New melt solver rebase #2138

jdannberg · 2018-03-23T01:00:26Z

Supersedes #1600. For now, I mainly want to see which tests fail.

not part of pull request (yet) [temp prm file for testing] add missing flag melt 3.0, constraints and everything add test files add comment add kd1 files iteration numbers and indent enable force vector in examples switch to CG for S, hack around other bugs add simple testcase add rising_melt_blob melt_solver/ files fix and cleanup compute_initial_pressure_field() add a postprocessor for the volumetric strain rate output p_c_bar, constant k_d with cutoff, DG p_c update:is_melt_cell - output darcy coefficient and is_melt_cell - split is_melt_cell decision - growing of layer (disabled) add tests, many not working clean up and fix melt advection add average density as fluid boundary pressure

jdannberg · 2018-03-27T16:00:08Z

@tjhei, do you know what's going on with the tester? It seems like it doesn't run tests any more...

jdannberg · 2018-03-27T23:20:40Z

Okay, I have checked that all the tests are still working and the output makes sense (it seems the tester is currently not working, but @gassmoeller has set up a tester locally on his machine so that I could update the test output), so it would be great if someone could have a look now! Although the pull request changes many files, most of them are just tests (and because we now use a new type of solver for the melt transport problem, all melt tests have changed).

The main advantages are that this pull request should make the solver for melt models much faster (in some cases it reduces the iteration count by an order of magnitude), in particular in cases with low or zero porosity in a part of the domain. So that would be a big step forward.

There are a number of things that are not clear to me and that we should discuss:

Boundary terms: @tjhei wondered if we have to add special boundary terms between the regions with melt and without melt, but I've tried to do this and the results just didn't make sense any more (and without these terms, all tests I've done looked good to me)
We have implemented the cutoff between melt cells/cells without melt as a fixed number with respect to the reference permeability, and it basically replaces the parameter "Melt transport threshold" we used before. Do we want to keep it like this (and just remove the Melt transport threshold from all input files, as it is not used anymore), or do we instead want to reuse the melt transport threshold somehow and try to make the cutoff an input parameter?
How do we want to compute the nonlinear residual? We introduced a new scaled pressure (p_c_bar), which is now stored in the block of the solution vector where the compaction pressure was before. So we will use p_c_bar to compute the nonlinear residual. Is that what we want? Or is there a good way to use the real compaction pressure instead?

gassmoeller

First set of comments

gassmoeller · 2018-03-29T23:41:21Z

include/aspect/compat.h

@@ -59,8 +59,9 @@ namespace aspect
      {
        dealii::SolverControl::State return_value = dealii::SolverControl::check(step, check_value);

-        if (step == 0)
-          history_data.resize(history_data.size()+1);
+        if (history_data_enabled)


Please try removing this and see if tests fail

gassmoeller · 2018-03-29T23:42:04Z

include/aspect/melt.h

+        double p_c_scale (const MaterialModel::MaterialModelInputs<dim> &inputs,
+                          const MaterialModel::MaterialModelOutputs<dim> &outputs,
+                          const MeltHandler<dim> &handler,
+                          bool consider_is_melt_cell) const;


gassmoeller · 2018-03-29T23:43:49Z

include/aspect/melt.h

+         * The last input argument consider_is_melt_cell determines if
+         * this computation takes into account if a cell is a "melt cell"
+         * (cells where we solve the melt transport equations, as
+         * indicated by the entries stored in the is_melt_cell vector of


Please rewrite in shorter sentences what happens if consider_is_melt_cell is true or false.

gassmoeller · 2018-03-29T23:49:41Z

include/aspect/melt.h

+
+      /**
+        * Given the Darcy coefficient as computed by the material model, limit the
+        * coefficient to a minimum value based on the reference Darcy coefficient


... to a minimum value (1e-3 times the reference Darcy coefficient) in melt cells. If @p is_melt_cell is false it returns zero.

gassmoeller · 2018-03-29T23:51:55Z

include/aspect/melt.h

+      /**
+       * Whether to use a discontinuous element for the compaction pressure or not.
+       */
+      bool use_discontinuous_p_c;


make private

gassmoeller · 2018-03-30T00:16:56Z

source/simulator/melt.cc

-                              :
-                              0.0);
+
+          const double K_D = this->get_melt_handler().limited_darcy_coefficient(melt_outputs->permeabilities[q] / melt_outputs->fluid_viscosities[q], p_c_scale > 0);


gassmoeller · 2018-03-30T00:17:32Z

source/simulator/melt.cc

@@ -483,8 +535,9 @@ namespace aspect
      const FEValuesExtractors::Scalar ex_p_f = introspection.variable("fluid pressure").extractor_scalar();
      const unsigned int p_f_component_index = introspection.variable("fluid pressure").first_component_index;
      const unsigned int p_c_component_index = introspection.variable("compaction pressure").first_component_index;
+      const double p_c_scale = dynamic_cast<const MaterialModel::MeltInterface<dim>*>(&this->get_material_model())->p_c_scale(scratch.material_model_inputs, scratch.material_model_outputs, this->get_melt_handler(), true);


gassmoeller · 2018-03-30T00:17:57Z

source/simulator/melt.cc

-                              melt_outputs->permeabilities[q] / melt_outputs->fluid_viscosities[q]
-                              :
-                              0.0);
+          const double K_D = this->get_melt_handler().limited_darcy_coefficient(melt_outputs->permeabilities[q] / melt_outputs->fluid_viscosities[q], p_c_scale > 0);


gassmoeller · 2018-03-30T00:20:50Z

source/simulator/melt.cc

            this->get_material_model().evaluate(in, out);

+            const double p_c_scale = dynamic_cast<const MaterialModel::MeltInterface<dim>*>(&this->get_material_model())->p_c_scale(in,out,this->get_melt_handler(), true);


gassmoeller · 2018-03-30T00:24:01Z

source/simulator/melt.cc

@@ -1180,7 +1265,12 @@ namespace aspect
      FEValues<dim> fe_values (this->get_mapping(),
                               this->get_fe(),
                               quadrature,
-                               update_quadrature_points | update_values);
+                               update_quadrature_points | update_gradients | update_values);


check if you can remove the gradients and change compute_material_model_input_values argument to false

I can set use_strain_rate to false, but we still need the gradients (as the material model inputs always compute the pressure gradients).

Oh, and I changed the compute_material_model_input_values to material_model_inputs.reinit, as this is our new way to write this, correct?

Ok then do that.

Yes I think it is.

gassmoeller · 2018-04-03T23:35:56Z

We discussed the mentioned questions with @jdannberg and @naliboff and came to the following conclusions:

Since the benchmarks work without the flux terms, and do not with them (maybe wrong terms?) we go ahead without the terms for now. If we need them we can always add them later.
We will remove the outdated parameter and introduce a more appropriate one, something like 'Relative permeability threshold'.
Since the scaling factor remains the same during nonlinear iterations (although it is different per cell) it might make no difference which one to use. Lets stick with the simple assumption that it is sufficient to use the nonlinear residual of the variable we actually have in the equation (pcbar) unless we have evidence we would need to do something else.

gassmoeller

Given the large number of comments I would like to take another look after you addressed them. Please keep your updates in a separate commit.

gassmoeller · 2018-04-03T23:41:10Z

source/simulator/melt.cc

@@ -1211,8 +1310,8 @@ namespace aspect
                const double phi = std::max(0.0, porosity_values[j]);

                double p = p_f_values[j];
-                if (phi < 1.0-this->get_melt_handler().melt_transport_threshold)
-                  p = (p_c_values[j] - (phi-1.0) * p_f_values[j]) / (1.0-phi);
+                if (p_c_scale > 0)


this should also prevent division by zero in the next line

gassmoeller · 2018-04-03T23:44:43Z

source/simulator/melt.cc

+  namespace
+  {
+    template <int dim>
+    struct PcAssembleData


Try to find a more expressive name, e.g. PcConstraintAssembleData and add a comment, e.g. This is a scratch object for the assembly of the Pc constraints outside of the melt region.

gassmoeller · 2018-04-03T23:46:22Z

source/simulator/melt.cc

+        material_model_inputs(scratch.material_model_inputs),
+        material_model_outputs(scratch.material_model_outputs)
+      {
+


remove empty line

gassmoeller · 2018-04-03T23:46:48Z

source/simulator/melt.cc

+    template <int dim>
+    struct PcAssembleData
+    {
+      PcAssembleData (const FiniteElement<dim> &finite_element,


// standard constructor

gassmoeller · 2018-04-03T23:47:01Z

source/simulator/melt.cc

+      {}
+
+
+      PcAssembleData (const PcAssembleData &scratch)


// Copy constructor

gassmoeller · 2018-04-04T00:30:41Z

source/simulator/parameters.cc

@@ -1287,8 +1287,7 @@ namespace aspect
                      && !use_discontinuous_composition_discretization,
                      ExcMessage ("Using discontinuous elements for temperature "
                                  "or composition in models with melt transport is currently not implemented."));
-          // We can not have a DG p_f. While it would be possible to use a
-          // discontinuous p_c, this is not tested, so we disable it for now.
+          // We can not have a DG p_f.
          AssertThrow(!use_locally_conservative_discretization,
                      ExcMessage ("Discontinuous elements for the pressure "
                                  "in models with melt transport are not supported"));


Add dot at end of sentence

gassmoeller · 2018-04-04T00:31:20Z

source/simulator/solver.cc

@@ -280,10 +280,10 @@ namespace aspect
      // first solve with the bottom left block, which we have built
      // as a mass matrix with the inverse of the viscosity
      {
-        SolverControl solver_control(1000, src.block(1).l2_norm() * S_block_tolerance);
+        aspect::SolverControl solver_control(1000, src.block(1).l2_norm() * S_block_tolerance);


do you need this? should be fixed on master already

gassmoeller · 2018-04-04T00:32:52Z

source/simulator/solver.cc


 #ifdef ASPECT_USE_PETSC
-        SolverCG<LinearAlgebra::Vector> solver(solver_control);
+        SolverGMRES<LinearAlgebra::Vector> solver(solver_control);


revert this to CG as before

gassmoeller · 2018-04-04T00:41:20Z

source/simulator/solver.cc

@@ -835,6 +836,7 @@ namespace aspect
                  {
                    // output solver history
                    std::ofstream f((parameters.output_directory+"solver_history.txt").c_str());
+                    f << std::setprecision(16);


Leave that out

gassmoeller · 2018-04-04T00:48:08Z

tests/rising_melt_blob.prm

+# simple model, but with a different melting rate.
+
+# The melting rate has the form of a Gaussian function 
+# and every timestep a fixed amount of material is melting. 


Update description

jdannberg · 2018-04-05T00:43:01Z

source/simulator/melt.cc

+
+
+        void
+        local_assemble (const typename DoFHandler<dim>::active_cell_iterator &cell,


Now that I am looking over this, should I also give this a different name (something like local_save_nonzero_pc_dofs)?

sure, sounds good

gassmoeller

The only functional change I found was in the preconditioner. Not sure that is the reason, but give it a try, since it was correct before and is wrong now.

gassmoeller · 2018-04-05T20:15:54Z

source/simulator/melt.cc

-                       pressure_scaling)
-                      * scratch.phi_p_c[i] * scratch.phi_p[j]
-                    )
-                    * JxW;


The above lines might be the cause for the test changes. Previously you always added these line (if(true), although that condition is unnecessary), now you only add them if the condition above is true (component_index == component_index). Revert to the old state, but remove the if(true). Its only the preconditioner though, not sure this is the only reason.

Ah, I see, I'll try it. Thanks!

gassmoeller · 2018-04-05T20:20:18Z

include/aspect/melt.h

-      MeltHandler(ParameterHandler &prm);
-
+    template <int dim>
+    struct Parameters


This struct should only contain the input or other parameters for MeltHandler, not internal variables. In particular boundary_fluid_pressure, is_melt_cell_vector, and current_constraints should be private member variables of the class and not in this struct. If you need them outside, make get_ functions for them.

gassmoeller · 2018-04-05T20:28:18Z

source/simulator/melt.cc

-                         "should probably be set to 'false'.");
+      prm.enter_subsection ("Melt settings");
+      {
+        prm.declare_entry ("Darcy coefficient variation threshold", "1e-3",


I do not like the variation a lot. What about Relative Darcy coefficient threshold? Variation sounds like you are measuring the spatial derivatives of darcy coefficients.

I feel like it's not really clear what a relative threshold should be, so that name feels a bit weird to me. I thought one other way to express this would be to say, how much is K_D allowed to vary (or deviate) from the reference value? I could call it "Global Darcy coefficient variation threshold" to make clear it doesn't have anything to do with derivatives.

Hmm, I am not really happy with that either. Lets discuss it later today.

I looked for parameters with a similar meaning in the code, and the closest thing I could find was something like the "Maximum lateral viscosity variation" in the Steinberger material model. But sure, let's discuss this later.

gassmoeller · 2018-04-05T20:28:50Z

source/simulator/melt.cc

+    {
+      prm.enter_subsection ("Melt settings");
+      {
+        K_D_variation_threshold = prm.get_double("Darcy coefficient variation threshold");


relative_K_D_threshold?

gassmoeller · 2018-04-05T20:29:35Z

source/simulator/parameters.cc

-                                  "or composition in models with melt transport is currently not implemented."));
-          // We can not have a DG p_f.
-          AssertThrow(!use_locally_conservative_discretization,
-                      ExcMessage ("Discontinuous elements for the pressure "


I like this part 😄, less melt code in the general code.

gassmoeller · 2018-04-05T20:30:29Z

tests/compaction_length_refinement.prm

@@ -110,10 +110,6 @@ subsection Model settings
  set Include melt transport                  = true
 end

-subsection Melt settings
-  set Melt transport threshold                = 0


Could this be what is causing the changes? Or was the melt transport threshold not used before as well?

jdannberg · 2018-04-06T05:40:43Z

I think I have addressed all of your comments. Can you have another look?

gassmoeller · 2018-04-07T00:39:43Z

Thanks for working through the multiple reviews. I think if we missed anything, we need to address it as a fix later. Melt 2.0 done, 🎉

bangerth · 2018-04-07T04:28:23Z

yay -- good work y'all!

rrgrove6 and others added 11 commits March 22, 2018 15:11

fix rebase

129a893

hack to update system matrices

d92512f

add refinement criterion for the compaction pressure

2259b3d

fix error computation for tests

4d19d41

fix solitary wave benchmark

0c63db2

add option for discontinuous p_c

b1f1ee3

only update melt constraints instead of recomputing all constraints

dc1e84b

cleanup

7bdd47b

more cleanup

032a525

use geometric mean for melt velocity

5b0340f

jdannberg force-pushed the new-melt-solver-rebase branch from 228471e to 2b2ac56 Compare March 23, 2018 15:53

fix rebase

f3b9aac

jdannberg force-pushed the new-melt-solver-rebase branch 5 times, most recently from 6775a48 to 027a018 Compare March 24, 2018 05:37

some more cleanup

6b0d37d

jdannberg force-pushed the new-melt-solver-rebase branch 3 times, most recently from 33837fc to a1865a3 Compare March 27, 2018 05:09

jdannberg force-pushed the new-melt-solver-rebase branch from a1865a3 to 863839d Compare March 27, 2018 23:17

jdannberg force-pushed the new-melt-solver-rebase branch 2 times, most recently from 82c1b6b to 2b3aa5e Compare March 29, 2018 03:28

jdannberg changed the title ~~[WIP] New melt solver rebase~~ New melt solver rebase Mar 29, 2018

fix tests

069c31b

jdannberg force-pushed the new-melt-solver-rebase branch from 2b3aa5e to 069c31b Compare March 29, 2018 22:47

jdannberg added this to the 2.0 milestone Mar 29, 2018

gassmoeller reviewed Mar 30, 2018

View reviewed changes

jdannberg mentioned this pull request Apr 3, 2018

add cookbook for a mid-ocean ridge model #2149

Merged

gassmoeller requested changes Apr 4, 2018

View reviewed changes

jdannberg commented Apr 5, 2018

View reviewed changes

jdannberg force-pushed the new-melt-solver-rebase branch from e984d99 to 399fb84 Compare April 5, 2018 02:27

gassmoeller reviewed Apr 5, 2018

View reviewed changes

jdannberg force-pushed the new-melt-solver-rebase branch 2 times, most recently from cd89ec5 to 291e56f Compare April 5, 2018 22:09

address comments

3aa6ca3

jdannberg force-pushed the new-melt-solver-rebase branch from 291e56f to 3aa6ca3 Compare April 6, 2018 02:24

Merge branch 'master' into new-melt-solver-rebase

3f9944e

some more fixes

45849cf

jdannberg force-pushed the new-melt-solver-rebase branch from 4d97497 to 45849cf Compare April 6, 2018 22:48

gassmoeller approved these changes Apr 6, 2018

View reviewed changes

jdannberg mentioned this pull request Apr 6, 2018

New melt solver #1600

Closed

gassmoeller merged commit 014b29d into geodynamics:master Apr 7, 2018

gassmoeller mentioned this pull request Apr 20, 2018

Fix parameters.tex and manual #2169

Merged

		this->get_material_model().evaluate(in, out);

		const double p_c_scale = dynamic_cast<const MaterialModel::MeltInterface<dim>*>(&this->get_material_model())->p_c_scale(in,out,this->get_melt_handler(), true);



		void
		local_assemble (const typename DoFHandler<dim>::active_cell_iterator &cell,

New melt solver rebase #2138

New melt solver rebase #2138

Conversation

jdannberg commented Mar 23, 2018

jdannberg commented Mar 27, 2018

jdannberg commented Mar 27, 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

Choose a reason for hiding this comment

Choose a reason for hiding this comment

gassmoeller commented Apr 3, 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

Choose a reason for hiding this comment

Choose a reason for hiding this comment

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 Apr 6, 2018

gassmoeller commented Apr 7, 2018

bangerth commented Apr 7, 2018