Time-dependent GO adaptation demo and time-dependent constants

[ddundo]

Closes #28.

In this PR I add a bubble shear goal-oriented adaptation demo. It involves a time-varying velocity field which is prescribed and not solved for. So to ease into the GO demo, I also add a simple solve_forward bubble shear demo to introduce the idea of "time-dependent" constants.

Related to #73 discussion: I added a check in GoalOrientedMeshSequence.indicate_errors to see whether the user flagged in their get_form that there are time-dependent constants. In that case, I recompile the form for each exported timestep, which should update these fields. By "flag" I mean that the form function contains the kwarg err_ind_time, i.e.: def form(index, fields, err_ind_time=None). I am not fully happy with this so I also didn't fully tidy up the go_mesh_seq.py code, as I hope to get feedback from you on how to improve this - but here is my reasoning behind it:

1. This does not change the interface when we do not have time-dependent constants,
2. We do not need to recompile the form for every timestep when we are not inside indicate_errors(). Instead we can (and should) update the time-dependent constants in get_solver.
3. The name err_ind_time is unique enough to be sure that the user won't randomly choose it.
4. This kind of "flag" is necessary to handle the 2nd part of Stephan's comment here https://github.com/pyroteus/goalie/issues/55#issuecomment-1702280655

Overall, I think it's a good solution but I am happy to keep #73 open and explore other ideas.

I will make 2 comments in the code myself where things are left undone.

And I tried to keep the demos cheap (to prevent the testing suite taking too long) so the results could be better, but they still look good I'd say!

[ddundo]

I cancelled the workflow because it looked like the GO demo was taking a while (takes about 2-3 minutes for me locally). We should cut it down even further (e.g. use only one fp iteration, decrease end_time...).

Edit: During testing, the mesh resolution, timestep and number of fp iterations are now modified to cut down the time. This helped but it's still quite long: 172s for the GO demo and 105s for the solve_forward one. Locally it takes 48s and 12s, respectively. I am surprised that it takes so long since it's solving it on a super coarse mesh (5x5). So feedback is welcome on how to get this down further :)

Edit2: I further reduced it (6 subintervals and increased timestep during testing). They take 18s and 6s for me locally. I cancelled the workflow now to not waste resources - will check the time here after your review when I push other changes.

[ddundo]

Thanks a lot for reviewing this @jwallwork23! I addressed the minor comments in the new commit.

My main comment is the one suggesting that we actually adopt time as a third (compulsory) argument to form. Interested to hear what you think about this.

I might be missing an obvious solution... but if we make this change, I don't see how we would separate out the problems where we do and do not have these time-dependent constants. So in error_indicators we would always pass time to get_form, which would mean recompiling the form unnecessarily in all other demos we have. That's what I meant when I said that the kwarg err_ind_time flags that we have these fields, so we need to call get_form in order to update them since they're not stored in the solution dictionaries. Do you have a suggestion how to deal with this please?

Edit: And to reply to this related comment here:

I wonder if we should just add time as an argument to form in general? For problems where the form doesn't vary with time, it will just be ignored. IMO it would be better to only pass solution fields (and their lagged counterparts) through the second argument and always determine other variables based on the current time. I guess this would imply some refactoring, including the other demos.

If I understood you right, this would mean calling get_form inside the timestepping loop which would be extremely inefficient, right? I.e. we would go from

def get_solver(mesh_seq):
    def solver(index, ic):
        ...
        form_fields = {"c": (c, c_), "u": (u, u_)}
        F = mesh_seq.form(index, form_fields)["c"]
        nlvp = NonlinearVariationalProblem(F, c, bcs=mesh_seq.bcs(index))
        nlvs = NonlinearVariationalSolver(nlvp, ad_block_tag="c")

        while t < t_end + 0.5 * dt:
            u.interpolate(velocity_expression(x, y, t))
            nlvs.solve()
            ...

to something like

def get_solver(mesh_seq):
    def solver(index, ic):
        ...
        sol_fields = {"c": (c, c_)}

        while t < t_end + 0.5 * dt:
            F = mesh_seq.form(index, sol_fields, time=t)["c"]
            solve(F==0, ...)
            ...

I'm trying to avoid this for the same reason as in the first part of this comment.

[jwallwork23]

Good point. I guess we need get_form to produce an update_form method, which could get stashed on the MeshSeq? Then - if it exists - it is called at each timestep?

[acse-ej321]

@ddundo @jwallwork23 - I don't think you would need to go to that extreme in changing the def solver function. For me, the current setup works,

where I pass t as a parameter to get_form:

 def get_form(self):
        def form(index, solutions,t,field="u"):

and define it in def solver outside the update_forcing or t stepping loop as a firedrake.Constant such that in the loop the fd.Constant is updated with assign:

 def get_solver(self):
        def solver(index, ic):

            # Time integrate from t_start to t_end
            t_start, t_end = self.subintervals[index]
            dt = self.time_partition.timesteps[index]
            # define and initialise t as a firedrake.Constant
            t = fd.Constant(t_start)
            # set the dummy t as an integer for the timestep loop
            t_= t_start

           # assemble the form outside the loop
           F = self.form(index, {"u": (u, u_)},t, field="u")["u"]

          while t_ < t_end - 1.0e-05:
                # solve
                fd.solve(F == 0, u,bcs=bcs, ad_block_tag="u")
                # reassign prior
                u_.assign(u)
                # update dummy time stepping int for loop
                t_ += dt
                # update firedrake.Constant time in form with 'assign'
                # updates the form without full reassemble
                t.assign(t_)

Firedrake will automatically assemble the form with the updated Constant as far as I am aware, as long as you assign it - so there would be no need to nest the assembly of the form in thet-stepping loop for that intention. That is not to say anupdate_form` might be more elegant?

[ddundo]

@acse-ej321 I may be misunderstanding (please correct me!), but if you take a look at L125-L180 in https://github.com/firedrakeproject/firedrake/blob/master/firedrake/solving.py, when you do fd.solve(F == 0, u,bcs=bcs, ad_block_tag="u"), firedrake assembles the NonlinearVariationalProblem. So when we do it inside the while loop, it does this for every timestep, which is very inefficient. So it's better to define the NonlinearVariationalProblem outside of the while loop and just call it inside the loop. I get about 10x speed increase when I do that.

[ddundo]

I don't see how we would separate out the problems where we do and do not have these time-dependent constants. So in error_indicators we would always pass time to get_form, which would mean recompiling the form unnecessarily in all other demos we have.

I went to test how expensive this is. I called get_form 10 000 times and it took only 21 seconds! So I think we can sacrifice this bit of inefficiency in indicate_errors and always recompile the form? The extra cost of this is dwarfed by everything else in indicate_errors :)

Good point. I guess we need get_form to produce an update_form method, which could get stashed on the MeshSeq? Then - if it exists - it is called at each timestep?

Was this referring to calling update_form inside the solver or indicate_errors? In get_solver it doesn't really matter for goalie, right? Users can do whatever they want there anyway. And I think we should keep it like this - I don't see a point in imposing some restrictions here.

So in summary: now I don't see a downside to making time a compulsory argument in form. What do you think @acse-ej321 and @jwallwork23? :)

[jwallwork23]

Yes, in general we should create the NVP and repeatedly its solve method rather than calling the solve driver function - I had just used that in the preliminary demos because its more intuitive for beginners.

[jwallwork23]

Okay how about this: we introduce a time attribute for MeshSeq, which is a list containing a Function from R-space for each subinterval, initialised to the start time of that subinterval. Inside the solver, we get this time and use it in the loop. Since it's a Function rather than a Constant, it will still support time += dt. The same time[subinterval] should be used inside the form, which would mean that updating the time in the solver loop would automatically update the form, too. Similar changes presumably required for the error indication code.

Some comments:

• Probably best to have the time attribute private and provide a public @property that grabs it.
• The time values would need to be reset to the start times of the corresponding subintervals before solving.
• Plausibly, this approach would avoid needing to introduce update_form.

[ddundo]

Thanks @jwallwork23, I'll do as you suggested :)

[ddundo]

Actually sorry, I'm not sure I see how that would allow us to address the last paragraph of https://github.com/pyroteus/goalie/issues/55#issuecomment-1702280655.

That is, in the gray_scott_split.py demo, in solver we have:

while t < t_end - 0.5 * dt:
    nlvs_a.solve()
    nlvs_b.solve()

which means that when we are in indicate_errors, we have to use the values of b and b_ from the previous timestep when compiling F_a:

# no changes needed for F_b

if inside_indicate_errors:
    b = b_             # b from the previous timestep
    b_ = b_old_old     # b_ from the previous timestep
F_a = ...

So in your suggestion @jwallwork23, what would this if inside_indicate_errors condition actually be?

E.g. if we had def form(index, solutions, time=None) then this condition could be if time is not None, since we do not pass time from solver but pass it from indicate_errors.

---

Edit: sorry, I jumped a step ahead... this is what I do locally in my glacier example, where I then completely removed

transfer(self.solutions[f][FWD][i][j], u[f])
transfer(self.solutions[f][FWD_OLD][i][j], u_[f])

but rather just pass the correct fields in forms = enriched_mesh_seq.form(i, mapping, time=time). I think this is a good solution because calling enriched_mesh_seq.form is surprisingly cheap.

[jwallwork23]

Sorry @ddundo I'm confused.

• I am no longer suggesting passing time as a kwarg. Instead I'm suggesting it to be an attribute of the MeshSeq, with one time Function for each subinterval.
• I don't like this conditional if inside_indicate_errors. Can we not avoid this?
• Let's not get ahead of ourselves - we can think about how to address other PRs later.

Edit: (I'm struggling to follow what you're saying without a code example...)

[ddundo]

Yeah, sorry about that... I realised I jumped ahead after I wrote it. I will have more time over the weekend so I will code up your suggestion and post an update :) I'm just trying to avoid refactoring in the future. Thanks again both for the input!

[ddundo]

I was trying out @jwallwork23's suggestion now and I have 2 ugly solutions to offer.

So to summarise, the suggestion was:

• introduce time = MeshSeq.time
• only pass c, c_ to get_form from the solver. Define u, u_ inside get_form

So this would look something like this:

def get_form(mesh_seq):
    def form(index, solutions):
        Q = mesh_seq.function_spaces["c"][index]
        V = VectorFunctionSpace(mesh_seq[index], "CG", 1)
        R = FunctionSpace(mesh_seq[index], "R", 0)
        dt = Function(R).assign(mesh_seq.time_partition.timesteps[index])
        theta = Function(R).assign(0.5)  # Crank-Nicolson implicitness parameter

        c, c_ = solutions["c"]

        # Velocity
        x, y = SpatialCoordinate(mesh_seq[index])
        u = Function(V).interpolate(velocity_expression(x, y, mesh_seq.time[index]))
        u_ = Function(V).interpolate(velocity_expression(x, y, mesh_seq.time[index] - dt))

        # SUPG stabilisation parameter
        D = Function(R).assign(0.1)  # diffusivity coefficient
        h = CellSize(mesh_seq[index])  # mesh cell size
        U = sqrt(dot(u, u))  # velocity magnitude
        tau = 0.5 * h / U
        tau = min_value(tau, U * h / (6 * D))

        phi = TestFunction(Q)
        phi += tau * dot(u, grad(phi))

        a = c * phi * dx + dt * theta * dot(u, grad(c)) * phi * dx
        L = c_ * phi * dx - dt * (1 - theta) * dot(u_, grad(c_)) * phi * dx
        F = a - L

        return {"c": F}
    return form

def get_solver(mesh_seq):
    def solver(index, ic):
        tp = mesh_seq.time_partition
        dt = tp.timesteps[index]
        num_timesteps = tp.num_timesteps_per_subinterval[index]

        t = mesh_seq.time[index]

        # Initialise the concentration fields
        Q = mesh_seq.function_spaces["c"][index]
        c = Function(Q, name="c")
        c_ = Function(Q, name="c_old").assign(ic["c"])

        F = mesh_seq.form(index, {"c": (c, c_)})["c"]
        nlvp = NonlinearVariationalProblem(F, c, bcs=mesh_seq.bcs(index))
        nlvs = NonlinearVariationalSolver(nlvp, ad_block_tag="c")

        # Time integrate
        for _ in range(num_timesteps):
            # TODO: Update u and u_ without reassembling nlvp

            # solve the advection equation
            nlvs.solve()

            # update the 'lagged' concentration
            c_.assign(c)
            t += dt

        return {"c": c}
    return solver

So we need to somehow get u, u_ from get_form without calling get_form repeatedly inside the timestepping loop. Here are the two solutions I can think of:

1. Get u, u_ from F.coefficients():

   
   def get_form(mesh_seq):
   def form(index, solutions):
       # same as above, just need to name u, u_
       u = Function(V, name="u").interpolate(...)
       u_ = Function(V, name="u_old").interpolate(...)
   
       return {"c": F}
   return form

def get_solver(mesh_seq): def solver(index, ic):

same as above up to this point

    F = mesh_seq.form(index, {"c": (c, c_)})["c"]
    for coeff in F.coefficients():
        coeff_name = coeff.name()
        if coeff_name == "u":
            u = coeff
        elif coeff_name == "u_old":
            u_ = coeff
    x, y = SpatialCoordinate(mesh_seq[index])

    nlvp = NonlinearVariationalProblem(F, c, bcs=mesh_seq.bcs(index))
    nlvs = NonlinearVariationalSolver(nlvp, ad_block_tag="c")

    # Time integrate
    for _ in range(num_timesteps):
        u.interpolate(velocity_expression(x, y, t))

        # solve the advection equation
        nlvs.solve()

        # update the 'lagged' concentration and velocity field
        c_.assign(c)
        u_.assign(u)
        t += dt

    return {"c": c}
return solver

~~2. Have `get_form` return `(u, u_)`:~~ No longer an option since #170.

Am I missing a more elegant solution? :)

[ddundo]

@ddundo Would you mind addressing my comment and rebasing on top of main?

I've done it but I didn't update the bubble shear demos. I'd prefer to wait for #164 to be merged before making further changes, since that properly takes care of #55, which overlaps with changes in this PR. So I think it would make sense to prioritise that PR before this one for that reason :)

[ddundo]

I converted this to ready for review to draw more attention to it :) as I mentioned in the meeting (and as discussed in the long thread above), it remains to decide how to communicate these "time-dependent constants" within Goalie.

The two new demos fail - I haven't updated them yet to yield etc.

[ddundo]

Hi @jwallwork23 @stephankramer @acse-ej321, here is the summary of the issue :) I will also start by saying that the PR is longish because there are 2 demos, but it is enough to just look at one of them (either one) for the purposes of the discussion here

In this specific bubble shear case example, we have a background velocity field that we do not solve for - we just update it at each timestep. How goalie is currently set up, when we call indicate_errors we only update the prognostic fields that we solved for (and which we saved at each export timestep). We also need to update this velocity field at each export timestep in order to compute the error indicators correctly.

So the problem we need to solve is how to update these fields in indicate_errors. The answer is not obvious since these fields are updated manually by the user in user-defined get_form or get_solver functions.

---

Here is the idea that I had how we could do this automatically, without requiring extra effort from the user.

We could identify these changing fields and extract them from the coefficients method of the variational form. I made a small gist to help demonstrate: https://gist.github.com/ddundo/92fc7d9fd24471a37c5a903ddd035554 - I labelled the important part towards the end. Here it's all inside the "user code", but it could all be done really nicely in the goalie code. I have a clear idea how to do this nicely :)

So we could extract and save a copy of these fields at each exported timestep, before then using them in indicate_errors. This is potentially memory-intensive, so I also suggest that we do the following. At the moment, goalie is set up so that we solve the forward and adjoint problem over all subintervals, and only then do we call indicate_errors. I think that we should instead call indicate_errors after each individual subinterval. That way:

• we only need to store these changing fields for the export timesteps in the subinterval, rather than over the whole time interval
• similarly for forward_old, adjoint and adjoint_next fields which we are now storing over the whole time interval: unless we want to return them (could add a kwarg for it), we could drop them after each call to indicate_errors and free up memory after each subinterval

[jwallwork23]

Apologies for the delay getting back to you @ddundo. This sounds like a great approach to me. It'd be preferable if you could break down the efforts into several small PRs, if possible. e.g., start by splitting up the indicate_errors call in one PR. That sound okay?

[ddundo]

Thanks @jwallwork23! I'll do it in steps as you suggested :)