Initialise fields outside of `get_solver`

[ddundo]

Instead of initialising fields inside get_solver, I think it would be better to do that somewhere inside MeshSeq. That way we do not require the user to properly name the functions ("field" and "field_old") and assigning initial conditions to field_old. I think it would be good to hide this inner working of goalie from the user.

So, for example, the burgers.py demo is currently:

def get_solver(mesh_seq):
    def solver(index, ic):
        function_space = mesh_seq.function_spaces["u"][index]
        u = Function(function_space, name="u")

        # Initialise 'lagged' solution
        u_ = Function(function_space, name="u_old")
        u_.assign(ic["u"])
        ...

but would be simplified to

def get_solver(mesh_seq):
    def solver(index):                  # no longer passing ic
        u, u_ = mesh_seq.fields["u"]    # now mesh_seq.fields returns a dictionary
        ...

What do you think @jwallwork23?

[jwallwork23]

Sounds like a great idea, thanks @ddundo!

[stephankramer]

Ok, so let's start with this one. I think you have to spell out a bit more what you are proposing. From the code above I gather you will have goalie prepare the contents of mesh_seq.fields, before goalie calls the solver, such that it contains the u and u_ that you require the user code to use as the end-of-timestep and beginning-of-timestep solution in each of the timesteps in the solver, I'm just trying to spell out the sort of assumptions you are making that you need to explain, in order to be able to discuss the changes you are proposing. So please correct me if this description is wrong. You need to be very explicit here and describe things in terms of user code, what goalie does and expects, and what the interface is like.

So that use of mesh_seq.fields seems a bit of a funny way to pass on some arguments to a function, why not just have these u and u_ passed in as argument?

My main concern is that you are now telling the user code which functions it has to use, rather than having the user create these and tell goalie which ones they are using. This doesn't work well with coupling goalie with existing models that already do this allocation for you - which is typically required because it (the existing model code) is also in charge of creating the solver objects (NonLinearVariationalSolver etc.) that have these wired in.

[ddundo]

Thanks @stephankramer, I did not consider the case where we might not want this - and I am still struggling to imagine the scenario which you are talking about without seeing an actual example :)

With the changes in #182, the dictionary MeshSeq.fields holds solutions at the current timestep (so we should probably rename it to something more informative, see https://github.com/mesh-adaptation/goalie/pull/182#discussion_r1580789499). So at the beginning of get_solver this is the dictionary of initial conditions at that subinterval. So yes, you are right - we can also pass these initial conditions as an argument to get_solver and get_form. I agree that it makes more sense that way. I only did it this way because I thought it's more in line with how we access other things, such as mesh_seq.function_spaces etc.

[ddundo]

And to answer more specifically...

From the code above I gather you will have goalie prepare the contents of mesh_seq.fields, before goalie calls the solver, such that it contains the u and u_ that you require the user code to use as the end-of-timestep and beginning-of-timestep solution in each of the timesteps in the solver (...) My main concern is that you are now telling the user code which functions it has to use, rather than having the user create these and tell goalie which ones they are using.

The new function I introduced, MeshSeq._reinitialise_fields prepares functions u and u_ based on the information the user provided in get_function_spaces - i.e. the name of the field and associated function space. So I wouldn't say that I am "requiring" anything or restricting the user. The user can still use whatever they want, they just need to specify it in get_function_spaces.

And this is how it is on the main branch too. We can't create a function in get_solver that has a different name and function space than that provided in get_function_spaces.

[ddundo]

So that use of mesh_seq.fields seems a bit of a funny way to pass on some arguments to a function, why not just have these u and u_ passed in as argument?

Actually the more I think about it, I'm not sure if I agree... if we do that, where do we stop? Should we also pass mesh_seq.function_spaces, mesh_seq.form, mesh_seq.time_partition, etc. as arguments? The final stop would be passing the subinterval mesh as an argument.

So if we take a look at burgers.py on the main branch:

def get_function_spaces(mesh):
    return {"u": VectorFunctionSpace(mesh, "CG", 2)}

def get_solver(mesh_seq):
    def solver(index, ic):
        function_space = mesh_seq.function_spaces["u"][index]

Now let's pretend that we are passing the subinterval mesh as an argument to the solver. Then we could do the following:

def get_solver(mesh_seq):
    def solver(index, ic, mesh):
        function_space = VectorFunctionSpace(mesh, "CG", 2)

which is identical to what is happening in the first snippet, but avoids unnecessarily redefining what function space we wanted. So mesh_seq.function_spaces also "prepares" function spaces in the same way that my proposed mesh_seq.fields would prepare the solution functions. It just follows whatever the user specified in get_function_spaces.