Newton merge 3, part 2: Newton drucker prager material #1684

MFraters · 2017-05-15T18:14:36Z

split up of #1653, part 2. This pull request changes the drucker prager material model to the new material model style and adds derivatives needed for the Newton method.

bangerth · 2017-05-15T18:17:23Z

Can anyone with better knowledge of material models look at this? @gassmoeller , @anne-glerum ?

bangerth · 2017-05-15T18:17:29Z

/run-tests

tjhei · 2017-05-15T20:43:14Z

source/material_model/drucker_prager.cc

+          const double pressure=std::max(in.pressure[i],0.0);
+
+          // calculate effective viscosity
+          if (in.strain_rate.size())


tjhei · 2017-05-15T20:43:31Z

source/material_model/drucker_prager.cc

+
+              const SymmetricTensor<2,dim> strain_rate_deviator = deviator(in.strain_rate[i]);
+
+              const double edot_ii_strict = (&(this->get_simulator()) == NULL


please add the long comment that you deleted in this place.

tjhei · 2017-05-15T20:44:07Z

source/material_model/drucker_prager.cc

+                  out.viscosities[i] = effective_viscosity;
+                }
+
+              Assert(dealii::numbers::is_finite(out.viscosities[i]),ExcMessage ("Error: Averaged viscosity is not finite."));


space after comma, maybe break line

tjhei · 2017-05-15T20:44:20Z

source/material_model/drucker_prager.cc

+
+                }
+            }
+          out.densities[i] = reference_rho * (1 - thermal_expansivity * (temperature - reference_T));


tjhei · 2017-05-15T20:44:27Z

source/material_model/drucker_prager.cc

+          out.thermal_conductivities[i] = thermal_conductivities;
+          // Compressibility at the given positions.
+          // The compressibility is given as
+          // $\frac 1\rho \frac{\partial\rho}{\partial p}$.


comment wrong

what should it be?

tjhei · 2017-05-15T20:45:34Z

tests/drucker_prager_derivatives.cc

+
+  if (Error)
+    {
+      std::cout << "Some parts of the test where not succesfull." << std::endl;


tjhei · 2017-05-15T20:45:44Z

tests/drucker_prager_derivatives.prm

@@ -0,0 +1 @@
+set Additional shared libraries = tests/libdrucker_prager_derivatives.so


add a comment please

tjhei · 2017-05-15T20:45:55Z

tests/drucker_prager_derivatives/screen-output

+oneone 3: Finite difference = 0. Analytical derivative = -0
+oneone 4: Finite difference = -5.03316e+25. Analytical derivative = -0
+   Error: The derivative of the viscosity to the strain rate is too different from the analitical value.
+Some parts of the test where not succesfull.


so the test fails?

ah, right. I had written a bit about that in the big pull request, but forgot to copy it to here to discuss.

I had written something on this in a email to wolfgang a long time ago:

In short the problem is that for all common cases I found that the finite difference and the analytical solution are exactly the same. But with a strain rate in which the diagonal entries are the same the analytical solution gives -0 and the finite difference gives results on the order of the other strain rates (~-1e+20).

Here is the relevant code (There is also a minimum to edot_ii in the real code.):

const double edot_ii = 2 * std::sqrt(0.5*deviator(in.strain_rate[i])*deviator(in.strain_rate[i])); const double stress_exponent_inv = (1./stress_exponent[c]); composition_viscosities[c] = std::pow(prefactor[c],-stress_exponent_inv) * std::pow(edot_ii,stress_exponent_inv-1); composition_viscosities_derivatives[c] = 2 * (stress_exponent_inv-1) * composition_viscosities[c] * (1/(edot_ii*edot_ii)) * deviator(in.strain_rate[i]);

I think this happens for the following reason: In the analytical equation, the deviator make the diagonal then zero, and the analytical solution then also ends up with zero entries on the diagonal. The finite difference method first gives the original strain rate tensor to the material model, and the deviator makes the diagonal zero. But then the finite difference method perturbs only one of the diagonal entries, which make the diagonal entries not equal any more and as a consequence the diagonal of the deviator will not contain zero entries any more, and I think this lead to a viscosity change on the order of what normally happens.

Thoughts?

gassmoeller

Looks ok except for some style issues. Thanks for unifying the functions. How rigorously did you test that the new viscosity computation is identical to the old one? What about corner cases (like epsilon=0)?

gassmoeller · 2017-05-15T20:44:20Z

source/material_model/drucker_prager.cc

+                  out.viscosities[i] = effective_viscosity;
+                }
+
+              Assert(dealii::numbers::is_finite(out.viscosities[i]),ExcMessage ("Error: Averaged viscosity is not finite."));


Please break this line

gassmoeller · 2017-05-15T20:45:48Z

source/material_model/drucker_prager.cc

+                    for (int y = 0; y < dim; y++)
+                      Assert(dealii::numbers::is_finite(derivatives->viscosity_derivative_wrt_strain_rate[i][x][y]),
+                             ExcMessage ("Error: Averaged dviscosities_dstrain_rate is not finite."));
+


Why do you separate this if-block from the one above? If you could unify them it would be easier to read the code, because all derivative stuff would be in one place.

gassmoeller · 2017-05-15T20:46:09Z

source/material_model/drucker_prager.cc

          reference_specific_heat    = prm.get_double ("Reference specific heat");
-          thermal_alpha              = prm.get_double ("Thermal expansion coefficient");
+          thermal_expansivity              = prm.get_double ("Thermal expansion coefficient");


align equal signs

gassmoeller · 2017-05-15T20:46:43Z

tests/drucker_prager_derivatives.cc

+  in_base.pressure[4] = 2e12;
+
+  /**
+   * We cant take to slow values, because then the difference in the


what are slow values? Do you mean small strain rates?

gassmoeller · 2017-05-15T20:48:17Z

tests/drucker_prager_derivatives.cc

+// run this function by initializing a global variable by it
+int ii = f(-1000); // Testing min function
+int iz = f(-2); // Testing min function
+int ij = f(-1.5); // Testing generalized p norm mean with negative p


I really dont think these comments apply to these function calls 😄. Please fix all of them.

gassmoeller · 2017-05-15T20:49:19Z

include/aspect/material_model/drucker_prager.h

+         * computes the volume fractions of the compositional fields, and returns it so that the sum
+         * of the volume fractions is one.
+         */
+        std::vector<double> compute_volume_fractions( const std::vector<double> &compositional_fields) const;


Is this function used in the .cc file? I do not see it.

you are right, this comes from code when there where different compositions present. I will remove this.

gassmoeller · 2017-05-15T20:50:08Z

include/aspect/material_model/drucker_prager.h

@@ -76,43 +76,17 @@ namespace aspect
     * @ingroup MaterialModels
     */
    template <int dim>
-    class DruckerPrager : public MaterialModel::InterfaceCompatibility<dim>, public ::aspect::SimulatorAccess<dim>
+    class DruckerPrager : public MaterialModel::Interface<dim>, public ::aspect::SimulatorAccess<dim>


Please read through the class documentation above and make sure it is still accurate after this change. If not please update, and do not forget to update the string at the bottom of the .cc file.

I don't see anything what should change based on this. Did you have anything in mind? I could add that the derivatives are added to the material model.

gassmoeller · 2017-05-15T20:50:51Z

include/aspect/material_model/drucker_prager.h

-                                             const double pressure,
-                                             const std::vector<double> &compositional_fields,
-                                             const Point<dim> &position) const;
+        virtual void evaluate(const MaterialModel::MaterialModelInputs<dim> &in,


Thanks for the cleanup 👍

gassmoeller · 2017-05-15T20:52:02Z

source/material_model/drucker_prager.cc

+                }
+              else
+                {
+                  out.viscosities[i] = effective_viscosity;


I suppose since the tests pass the viscosity calculation is unchanged? I did not go through all of the formulas, but did you intentionally change something, or should everything be as before? Except for the addition of derivatives of course.

The viscosity should be exactly the same as before.

bangerth mentioned this pull request May 15, 2017

Newton merge 3, part 1: Add newton simple nonlinear material #1683

Merged

MFraters mentioned this pull request May 15, 2017

Newton merge 3: To the spd factor visualisation and (a bit too far) beyond #1653

Closed

MFraters force-pushed the newton_drucker_prager_material branch from aa24b5e to 44e1246 Compare May 15, 2017 18:31

tjhei requested changes May 15, 2017

View reviewed changes

gassmoeller reviewed May 15, 2017

View reviewed changes

MFraters added 2 commits May 16, 2017 05:23

rewrite drucker prager material model and add derivates.

7eb2b4b

Add test for the derivatives of the drucker prager material model.

f28e722

MFraters force-pushed the newton_drucker_prager_material branch from f00f38e to 7049917 Compare May 16, 2017 03:28

add changelog for drucker prager derivative changes.

4159efa

MFraters force-pushed the newton_drucker_prager_material branch from 7049917 to 4159efa Compare May 16, 2017 04:19

Merge branch 'master' into newton_drucker_prager_material

f2fe0ee

gassmoeller merged commit 33d570a into geodynamics:master May 16, 2017

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Newton merge 3, part 2: Newton drucker prager material #1684

Newton merge 3, part 2: Newton drucker prager material #1684

MFraters commented May 15, 2017

bangerth commented May 15, 2017

bangerth commented May 15, 2017

tjhei May 15, 2017

tjhei May 15, 2017

tjhei May 15, 2017

tjhei May 15, 2017

tjhei May 15, 2017

MFraters May 16, 2017

tjhei May 15, 2017

tjhei May 15, 2017

tjhei May 15, 2017

MFraters May 15, 2017

gassmoeller left a comment

gassmoeller May 15, 2017

gassmoeller May 15, 2017

gassmoeller May 15, 2017

gassmoeller May 15, 2017

gassmoeller May 15, 2017

gassmoeller May 15, 2017

MFraters May 15, 2017

gassmoeller May 15, 2017

MFraters May 16, 2017

gassmoeller May 15, 2017

gassmoeller May 15, 2017

MFraters May 15, 2017


		const SymmetricTensor<2,dim> strain_rate_deviator = deviator(in.strain_rate[i]);

		const double edot_ii_strict = (&(this->get_simulator()) == NULL

		@@ -0,0 +1 @@
		set Additional shared libraries = tests/libdrucker_prager_derivatives.so

Newton merge 3, part 2: Newton drucker prager material #1684

Newton merge 3, part 2: Newton drucker prager material #1684

Conversation

MFraters commented May 15, 2017

bangerth commented May 15, 2017

bangerth commented May 15, 2017

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

Choose a reason for hiding this comment

Choose a reason for hiding this comment