Multiple variables: only one field being updated at a time in goal based indicate_errors

[ddundo]

Hey @jwallwork23 and @stephankramer. As GoalOrientedMeshSeq.indicate_errors() is currently written, we are updating fields used in computing error indicators by first looping through fields and then looping through solutions at exported timesteps (L194-L201). This means that if we have more than one field in our mesh sequence, the updated fields won't represent the same timestep since only one field is being updated at the current timestep in the loop. This is an issue, for example, in the Gray-Scott example that we have.

In this multiple variable case, is it enough to just update all fields at the current timestep, something like this (just adding another loop)

for f in self.fields:
    # Loop over each timestep
    for j in range(len(sols[f]["forward"][i])):
        # Update fields
        for ff in self.fields:
          transfer(sols[ff][FWD][i][j], u[ff])
          transfer(sols[ff][FWD_OLD][i][j], u_[ff])
          transfer(sols[ff][ADJ][i][j], u_star[ff])
          transfer(sols[ff][ADJ_NEXT][i][j], u_star_next[ff])

or do we also need to modify the other calculations further down the line?

[jwallwork23]

Hi @ddundo we just looked into this in the meeting.

Essentially, the problem is the code was designed for single equations with single prognostic fields, but was then incorrectly extended to the case of multiple equations with multiple prognostic fields. The first f loop isn't really looping over fields, but over equations. The second ff loop in your suggestion then loops over the prognostic fields. Provided that indicator_fn is something general like the default get_dwr_indicator then possibly your suggestion would work, although it does duplicate effort. Also, if the user customises it then maybe not because they would be forced to use the same custom indicator for each equation.

Perhaps an alternative approach is just to reorder the loops? i.e., put the j loop first and a single f loop inside to update all the fields. Then to address the custom indicator_fn issue, I guess we would need to allow the user to pass multiple such functions. However, since no one has asked for that as yet, maybe we just open a new issue and leave it for now.

Does this sound like a good plan to you? If so, would you be willing to give it a test and open a PR (and a new issue)?

[jwallwork23]

Oh and I think you would want to loop over the fields once to update them all and then loop again (effectively over equations) to compute the error indicators. Both of those loops on the same level within the j loop.

[stephankramer]

Yes, I think you want to swap the two loops, i.e. j loop outside, f-loop inside - but then split the j loop to first loop through all f to do the transfers, and then in a second f loop do call the indicator_fn. So:

            # Loop over each timestep
            for j in range(len(sols[f]["forward"][i])):
                for f in self.fields:
                    # Update fields
                    transfer(sols[f][FWD][i][j], u[f])
                    transfer(sols[f][FWD_OLD][i][j], u_[f])
                    transfer(sols[f][ADJ][i][j], u_star[f])
                    transfer(sols[f][ADJ_NEXT][i][j], u_star_next[f])

                    # Combine adjoint solutions as appropriate
                    u_star[f].assign(0.5 * (u_star[f] + u_star_next[f]))
                    u_star_e[f].assign(
                        0.5 * (sols_e[f][ADJ][i][j] + sols_e[f][ADJ_NEXT][i][j])
                    )
                    u_star_e[f] -= u_star[f]

                for f in self.fields:
                    # Evaluate error indicator
                    indi_e = indicator_fn(forms[f], u_star_e[f])

                    # Project back to the base space
                    indi = project(indi_e, P0_spaces[i])
                    indi.interpolate(abs(indi))
                    indicators[f][i][j].interpolate(ufl.max_value(indi, 1.0e-16))

The second loop needs to be separate because the forms[f] does not just contain u[f] and u_[f] but potentially any of the other fields, so they should have been transferred at that stage.

The other thing is you'd have to be quite careful about how you use "the other" field in your equation. So say in the split Gray-Scott example: In the second solve we are solving for b but using the value of a that we have just solved for, so it is using the end-of-timestep value of a - which is consistent with the above. However in the first solve, which solves for a but uses the value of b, the value of b it uses is actually the beginning-of-timestep value, so in the form it should be using b instead of b. It doesn't matter in the forward solver, because b=b at the beginning of the timestep when that first solve is done - but in the error indicator calculation it would substitute the end-of-timestep value for b for the indicator of that first solve which wasn't even available to it.

[jwallwork23]

Hm yeah good points in the final paragraph. I think tackling those kinds of things properly will mean rethinking how we work with different time integration schemes in general.

[jwallwork23]

@ddundo I can't remember, are you working on this? If so, I'll wait and we can work on getting the fix merged. If not, I'll transfer the issue over to Goalie.

[ddundo]

Hey @jwallwork23 I've done this in a quick and dirty way just to get it working for my case (dirty referring to the way I save the solutions at the timestep prior to the exported timestep, re Stephan's comment). I was thinking of coming back to it in a month or so properly

[jwallwork23]

Okay, thanks for the update. If you are planning to revisit then I will transfer it over to Goalie, along with the other Issues I'm porting.