Add Darcy field convection timestep

[danieldouglas92]

Add the functionality to limit the time step based on the Darcy velocity. Currently, when the 'darcy field' advection method is used, a field is advected with the Darcy velocity but the time step is not limited by this velocity. This leads to convergence issues, since the Darcy velocity is always equal to or greater than the solid velocity. This adds the functionality to convection_time_step.cc to limit the time step based on the Darcy velocity.

The implementation I have in this PR is currently working in that it gives the proper limited time step, however the fluid parameters are currently hard-coded in because when I try to access MeltOutputs through a variable fluid_out, fluid_out is a nullptr. @jdannberg am I missing something obvious here? Once I can access the Melt variables, I think this PR is basically good to go.

This addresses one of the issues in #4748

[gassmoeller]

I left some comments that should fix the issue. Nice test. Can you also add a changelog entry?

[gassmoeller]

move this to the top of the file, where you also compute consider_darcy_timestep

[danieldouglas92]

Unless I'm misunderstanding your comments, I don't know if it's possible to not compute the porosity_idx here because you need to divide by the porosity to calculate the Darcy velocity

[gassmoeller]

But you can compute the porosirty_idx above, and use it here. See my comment above.

[danieldouglas92]

Thanks for the comments @gassmoeller, I still wasn't able to get fluid_out to not be a nullptr, and I'm not sure if I followed exactly all of your comments about what to calculate when determining whether consider_darcy_timestep should be true.

Just to outline my thought process, I think that the following line:
consider_darcy_timestep = this->get_parameters().compositional_field_methods[this->introspection().compositional_index_for_name("porosity")] == Parameters<dim>::AdvectionFieldMethod::fem_darcy_field

is important because the field type could be "porosity" but not be advected with the darcy field method. I think I was able to follow all your comments, but the only confusing part to me was bringing the line:

const unsigned int porosity_idx = this->introspection().compositional_index_for_name("porosity");

up to where we calculate consider_darcy_timestep.

[gassmoeller]

Looks much better. A few more comments.
I think your fluid_out is empty, because your material model does not actually create a MeltOutputs object in the create_additional_named_outputs function. You need to modify ReactiveFluidTransport to do that (currently it only fills the values in evaluate if the MeltOutputs exist, but they never exist, so it never fills them).

[gassmoeller]

But you can compute the porosirty_idx above, and use it here. See my comment above.

[danieldouglas92]

Got the melt outputs working thanks @gassmoeller!

[jdannberg]

A point on the MeltOutputs: They are not additional named outputs, they are a separate type of outputs (AdditionalMaterialOutputs). In all other places, it is the responsibility of the caller to create them, not the material model itself. None of the melt models create their own MeltOutputs in the create_additional_named_outputs function.

To be consistent with the rest of ASPECT, the melt outputs should be created in the place where you need them, with MeltHandler<dim>::create_material_model_outputs(out), before you call evaluate (see my comment on #5902).

[danieldouglas92]

@jdannberg I had tried this actually but for whatever reason at the time it didn't work and I clearly didn't try that hard to get it working, but now it seems to work!

[danieldouglas92]

@gassmoeller Unfortunately I'm running into a git issue, the test viscoplastic_melt_fraction is failing here because the .prm file asked for named additional outputs, but somehow the PR (#5799) that I got merged at the hackathon which actually allows this to be done is not here. My local main is up to date and does have all the code required for viscoplastic_melt_fraction to not fail, but when I try to rebase onto main it says that this branch (limit_timestep_by_darcy_velocity) is already up to date with main, which clearly isn't true. I'm not really sure how this is happening?

[gassmoeller]

A bit hard to debug remotely, but if I check out your branch from this PR, I clearly see the commit from #5799 in the commit history. Are you sure it is functionality from #5799 that is missing for this PR to pass the tests?
Here is the test error message: https://github.com/geodynamics/aspect/actions/runs/9532699918/job/26275098979?pr=5856#step:6:2060
The error message does not mention named additional outputs.

[danieldouglas92]

@gassmoeller Yes there are two tests that fail, the first is uncoupled_two_phase_tian_parameterization_kinematic_slab fails because of the changes that I'm adding here and I need to fix this test also, but the second test viscoplastic_melt_fraction fails because of the named additional outputs. But you're right, I do also see the merge in my git log but the actual code that was merged in that PR is not in this branch.

[gassmoeller]

But you're right, I do also see the merge in my git log but the actual code that was merged in that PR is not in this branch.

The only way I can see that happening is if either you (in this PR) or someone else (in another PR) undid the changes. Can you check here if you maybe undid your changes accidentally?

[danieldouglas92]

@gassmoeller It never once occurred to me that I might've accidentally deleted the code and committed it! Thanks for catching that

[gassmoeller]

/rebuild

[gassmoeller]

Is this ready for a review?

[danieldouglas92]

@gassmoeller yes!

[gassmoeller]

Looks good, just two small comments left.
@jdannberg could you take a look at the equation for the darcy flow? I am not qualified for that.

[gassmoeller]

This line would need to happen before the call to evaluate, right?

[danieldouglas92]

I just tried putting the line:

this->get_material_model().fill_additional_material_model_inputs(in, this->get_solution(), fe_values, this->introspection());

before evalute, after evalute, and not including at all, and it seems to give the same solution for all 3 cases so I'll remove it and see if it causes any of the other tests to fail.

[danieldouglas92]

the line:

this->get_material_model().create_additional_named_outputs(out);

also doesn't seem to be required, so I'll remove this line too.

[gassmoeller]

Sorry for the dealy, looks good to go.