fix freezing in melt simple material model and improve documentation #2131

jdannberg · 2018-03-17T05:28:29Z

This was really confusing before (and I noticed this when @joeschools showed me one of his model setups). What used to happen was that there was a freezing rate parameter, but even if that was set to zero, melt could freeze, but only under certain circumstances (if it had not moved away from its source region). This update makes this behaviour more consistent (melt does not freeze if the freezing rate is set to zero; and when the freezing rate is larger than zero, it determines the rate of all freezing) and clearly documents how the amount of freezing is computed. I also improved some variable names.

bangerth · 2018-03-19T00:22:53Z

There are three failing tests. Can you check whether they indicate something significant, or just need updated output files?

bangerth · 2018-03-19T00:24:25Z

I think @tjhei should probably look at this.

fix for fractional melting

jdannberg · 2018-03-20T18:01:14Z

I have updated the tests, and I've also tested the code with a simple setup for melting and freezing at a mid-ocean ridge that I'd like to add as a cookbook later (it only works with our new melt solver, so we have to wait until that is merged).

I believe anyone could have a look at this, not necessarily just @tjhei, as it doesn't have to do anything with the melt solver, it is just a fix for the parametrization of melting/freezing (so it goes more in the direction of thermodynamics).

gassmoeller

We had a lengthy discussion about this PR, and came to the following conclusion: This material model is already very complicated and the fix of this PR is correct, but adds another layer of complexity. In order to fix this, we will for now keep the changes of this PR, but require the user from now on to use operator splitting when using this material model. This is the only reasonable option anyway (without operator splitting the models become unstable and the nonlinear solver often does not converge), and by removing the code paths for no-operator splitting we reduce the complexity. So in short: We implement this fix, and remove an option that is unstable anyway to end up with a comparable complexity of the code.

gassmoeller · 2018-03-20T20:50:34Z

source/material_model/melt_simple.cc

@@ -258,7 +258,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())
+          if (this->include_melt_transport() && in.strain_rate.size())


gassmoeller · 2018-03-20T20:50:54Z

source/material_model/melt_simple.cc

+                    porosity_change += freezing_amount;
+
+                    // Adapt time scale of freezing with respect to melting.
+                    if (porosity_change < 0 && this->get_parameters().use_operator_splitting)


Remove the non-operator split versions.

jdannberg · 2018-03-21T03:53:09Z

Okay, I changed the material model so that it can only be used with operator splitting now and I could remove some lines of code. I adapted the tests accordingly.

gassmoeller

I like this much better now. A few small comments and questions. Let me know when you have addressed them, likely this is ready to merge after the small changes.

gassmoeller · 2018-03-21T17:00:18Z

doc/modules/changes/20180316_jdannberg

@@ -0,0 +1,9 @@
+Overhaul of the melt simple model. 


melt simple material model

gassmoeller · 2018-03-21T17:04:03Z

source/material_model/melt_simple.cc

+
+                    // Adapt time scale of freezing with respect to melting.
+                    if (porosity_change < 0 )
+                      porosity_change *= freezing_rate * melting_time_scale;


Can you add a comment here why you multiply by both? "E.g. we divide by the melting_time_scale later on in any case, but this is not necessary for this case, thus ...". That would make it easier to understand.

gassmoeller · 2018-03-21T17:05:09Z

source/material_model/melt_simple.cc

                    reaction_rate_out->reaction_rates[i][c] = 0.0;
                }
            }

-          out.entropy_derivative_pressure[i]    = entropy_change (in.temperature[i], this->get_adiabatic_conditions().pressure(in.position[i]), maximum_melt_fractions[i], NonlinearDependence::pressure);


What happened to these lines? Are they not important for latent heat?

I moved them further up into the individual branches (with and without melt transport). This was convenient because they need the depletion field, which we do not necessarily have in the case without melt, and so we compute it only in the part of the if statement that is executed with melt transport (otherwise, 0 is used, just as before these changes).

gassmoeller · 2018-03-21T17:06:53Z

source/material_model/melt_simple.cc

-                             "Units: $1/yr$.");
+                             "Freezing rate of melt when in subsolidus regions. "
+                             "If this parameter is set to a number larger than 0.0, it specifies the "
+                             "fraction of melt that will freeze per year, as soon as the porosity "


per year or per second

gassmoeller · 2018-03-21T17:07:15Z

source/material_model/melt_simple.cc

          prm.declare_entry ("Melting time scale for operator splitting", "1e3",
                             Patterns::Double (0),
-                             "In case the operator splitting scheme is used, the porosity field can not "
+                             "Because the operator splitting scheme is used, the porosity field can not "


Nice rewrite 😄

gassmoeller · 2018-03-21T17:07:59Z

source/material_model/melt_simple.cc

@@ -742,8 +759,7 @@ namespace aspect
    void
    MeltSimple<dim>::create_additional_named_outputs (MaterialModel::MaterialModelOutputs<dim> &out) const
    {
-      if (this->get_parameters().use_operator_splitting
-          && out.template get_additional_output<ReactionRateOutputs<dim> >() == NULL)
+      if (this->get_parameters().use_operator_splitting && out.template get_additional_output<ReactionRateOutputs<dim> >() == NULL)


do you still need the first condition?

Yes. I first removed it, but it caused some kind of problem. I don't remember what exactly it was, but I think it might have been related to the point that we still allow this material to be used without operator splitting if melt transport is switched off (in this case, there are no reactions, so it doesn't matter), and I got some error message in one of the tests that my material model asks for reaction outputs, but they are not filled, or something like that.

gassmoeller · 2018-03-21T17:09:10Z

tests/shear_heating_with_melt.prm

 end

 subsection Postprocess

-  set List of postprocessors = composition statistics,velocity statistics, mass flux statistics, heating statistics, temperature statistics


the visualization was likely added by accident?

I can remove it again, but if anyone needs to look at this test again in the future because something has changed, they will need the visualization output again, because there's no other way to really see what the test does and if its behaviour is still correct. So unless you feel strongly about removing it, I think it would be kind of useful to just leave it switched on.

jdannberg · 2018-03-21T18:22:28Z

Alright, I think I addressed all of your comments.

gassmoeller · 2018-03-22T02:05:19Z

Thanks, easier to understand now.

jdannberg force-pushed the fix_freezing_in_melt_simple branch 2 times, most recently from 3e8fa02 to 04f5c1a Compare March 20, 2018 04:06

fix freezing in melt simple material model and improve documentation

5dd5692

fix for fractional melting

jdannberg force-pushed the fix_freezing_in_melt_simple branch from fc59a8a to 5dd5692 Compare March 20, 2018 17:19

gassmoeller reviewed Mar 20, 2018

View reviewed changes

jdannberg force-pushed the fix_freezing_in_melt_simple branch 2 times, most recently from 8bbaf8f to fc33b89 Compare March 21, 2018 03:11

gassmoeller reviewed Mar 21, 2018

View reviewed changes

remove operator splitting option

e6b48cf

jdannberg force-pushed the fix_freezing_in_melt_simple branch from fc33b89 to e6b48cf Compare March 21, 2018 18:23

gassmoeller merged commit 870a37f into geodynamics:master Mar 22, 2018

jdannberg deleted the fix_freezing_in_melt_simple branch March 22, 2018 18:20

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

fix freezing in melt simple material model and improve documentation #2131

fix freezing in melt simple material model and improve documentation #2131

jdannberg commented Mar 17, 2018

bangerth commented Mar 19, 2018

bangerth commented Mar 19, 2018

jdannberg commented Mar 20, 2018

gassmoeller left a comment

gassmoeller Mar 20, 2018

gassmoeller Mar 20, 2018

jdannberg commented Mar 21, 2018

gassmoeller left a comment

gassmoeller Mar 21, 2018

gassmoeller Mar 21, 2018

gassmoeller Mar 21, 2018

jdannberg Mar 21, 2018

gassmoeller Mar 21, 2018

gassmoeller Mar 21, 2018

gassmoeller Mar 21, 2018

jdannberg Mar 21, 2018

gassmoeller Mar 21, 2018

jdannberg Mar 21, 2018

jdannberg commented Mar 21, 2018

gassmoeller commented Mar 22, 2018

fix freezing in melt simple material model and improve documentation #2131

fix freezing in melt simple material model and improve documentation #2131

Conversation

jdannberg commented Mar 17, 2018

bangerth commented Mar 19, 2018

bangerth commented Mar 19, 2018

jdannberg commented Mar 20, 2018

gassmoeller left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jdannberg commented Mar 21, 2018

gassmoeller left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jdannberg commented Mar 21, 2018

gassmoeller commented Mar 22, 2018