Fix operator splitting with latent heat #2106
Conversation
jdannberg
force-pushed
the
fix_operator_splitting_with_latent_heat
branch
from
0b94d78
to
10d391c
Compare
| @@ -131,6 +131,7 @@ namespace aspect | |||
| double compressibility; | |||
| double melt_compressibility; | |||
| bool include_melting_and_freezing; | |||
| double melting_rate_op_split; | |||
There was a problem hiding this comment.
I am not happy about this name. If at all, it should be melting_rate_for_operator_splitting, but that is not good either. What about melting_timescale? That would be the inverse of the current parameter. Also the fact that in the description for the parameter file you constantly have to refer to the inverse of this parameter hints to the fact that we maybe should store the inverse instead. What do you think?
There was a problem hiding this comment.
I could also just call it melting_rate, I guess.
I just though it would be more intuitive to give a rate, but maybe you're right and a time scale would make more sense. I guess it would be similar to a half life, so yes, I'll change it to melting_time_scale.
| @@ -133,6 +133,7 @@ namespace aspect | |||
| bool model_is_compressible; | |||
| bool fractional_melting; | |||
| double freezing_rate; | |||
| double melting_rate_op_split; | |||
| @@ -158,7 +158,7 @@ namespace aspect | |||
| out.densities[i] = (reference_rho_s + delta_rho) * temperature_dependence | |||
| * std::exp(compressibility * (in.pressure[i] - this->get_surface_pressure())); | |||
|
|
|||
| if (this->include_melt_transport() && include_melting_and_freezing) | |||
| if (this->include_melt_transport() && include_melting_and_freezing && in.strain_rate.size()) | |||
There was a problem hiding this comment.
Is this related to this PR or a separate improvement? Not saying you have to open another PR, just trying to understand if I am missing something, why this is now required.
There was a problem hiding this comment.
It is not really related, it is more of a bugfix. If you look below, you'll see that it was always required to compute melting rates, and I just moved it further up (so the code doesn't get executed if we don't have the strain rate).
In addition, before, we always computed the reaction rate outputs for the porosity, but the one for peridotite only when the strain rate was available. This didn't make sense (as both should always change in the same way), so the new syntax makes it more consistent.
source/material_model/melt_global.cc
Outdated
| @@ -198,9 +200,9 @@ namespace aspect | |||
| if (reaction_rate_out != NULL) | |||
| { | |||
| if (c == peridotite_idx && this->get_timestep_number() > 0) | |||
| reaction_rate_out->reaction_rates[i][c] = out.reaction_terms[i][c] / this->get_timestep() ; | |||
| reaction_rate_out->reaction_rates[i][c] = depletion_change * melting_rate_op_split; | |||
There was a problem hiding this comment.
If the melting_rate_op_split would be a timescale this would look more similar to the old thing just with a changed timescale.
source/material_model/melt_global.cc
Outdated
| "to be larger or equal to the reaction time step used in the operator splitting " | ||
| "scheme, otherwise reactions can not be computed. If the model does not use " | ||
| "operator splitting, this parameter is not used. " | ||
| "If this parameter is set to 0, no melting and freezing occurs in the model. " |
There was a problem hiding this comment.
Ah, is this the reason you want to use the inverse? Is there no better way to allow for no freezing / no melting?
There was a problem hiding this comment.
You could still set the time scale to 1e100 and have essentially no melting/freezing (and this is something that only comes up in test cases anyway), so I guess this is fine.
jdannberg
force-pushed
the
fix_operator_splitting_with_latent_heat
branch
from
10d391c
to
198f11f
Compare
|
Okay, I rebased and addressed your comment. Can you have another look? |
jdannberg
force-pushed
the
fix_operator_splitting_with_latent_heat
branch
from
198f11f
to
f3916bf
Compare
source/material_model/melt_global.cc
Outdated
| "Also note that the melting time scale has to be larger than or equal to the reaction " | ||
| "time step used in the operator splitting scheme, otherwise reactions can not be " | ||
| "computed. If the model does not use operator splitting, this parameter is not used. " | ||
| "Units: 1/yr or 1/s, depending on the ``Use years " |
source/material_model/melt_simple.cc
Outdated
| "Also note that the melting time scale has to be larger than or equal to the reaction " | ||
| "time step used in the operator splitting scheme, otherwise reactions can not be " | ||
| "computed. If the model does not use operator splitting, this parameter is not used. " | ||
| "Units: 1/yr or 1/s, depending on the ``Use years " |
fix tests address comments
jdannberg
force-pushed
the
fix_operator_splitting_with_latent_heat
branch
from
f3916bf
to
3012a8e
Compare
When looking at a model setup from @joeschools, I noticed that there is a problem when the melt material models are used together with the "latent heat melt" heating model in a model that uses operator splitting. What happens here is that a combination of the time step size and the number of operator splitting steps decides about how fast the porosity reaches the equilibrium melt fraction; in other words, how fast melting and freezing takes places; and there is no way to set the parameters in such a way that equilibrium is reached immediately.
This is not really great, as changing the time step size should not change the physics of the problem, so we thought about an alternative, and @gassmoeller suggested to introduce a separate parameter that controls the melting and freezing rate. This allows the user to control in the input file how fast melting and freezing occurs, independent of the number of operator splitting steps. It also means that the operator splitting time step has to be adapted to that melting rate. This allows it to make the melting rate so large that melt is basically always in equilibrium, and it also allows it to use melting rates inferred from experiments or observations. It also makes more sense, because the reaction rate outputs of the material models are now really "rates", and not just changes of something divided by the time step size.
This pull request introduces the new functionality, and adds a test for each material model to check that the latent heat of melting now works correctly.
This is based on #2104.