geodynamics / aspect

A parallel, extensible finite element code to simulate convection in both 2D and 3D models.
https://aspect.geodynamics.org/
Other
217 stars 232 forks source link

Add Darcy field convection timestep #5856

Closed danieldouglas92 closed 1 week ago

danieldouglas92 commented 3 weeks ago

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

danieldouglas92 commented 3 weeks ago

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.

danieldouglas92 commented 2 weeks ago

Got the melt outputs working thanks @gassmoeller!

jdannberg commented 2 weeks ago

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 commented 2 weeks ago

@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 commented 2 weeks ago

@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 commented 2 weeks ago

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 commented 2 weeks ago

@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 commented 2 weeks ago

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 commented 2 weeks ago

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

gassmoeller commented 2 weeks ago

/rebuild

gassmoeller commented 2 weeks ago

Is this ready for a review?

danieldouglas92 commented 2 weeks ago

@gassmoeller yes!