Add derivatives to visco plastic material model #2030

MFraters · 2017-12-13T17:39:49Z

Adding finite difference derivatives to the visco-plastic material model, so it can be used with the Newton solver. With this finite difference model, it should be relatively easy to add derivatives to other material models.

The finite difference could actually be done somewhere higher up or abstracted into a single function.

naliboff · 2017-12-18T20:47:05Z

@MFraters - Great! Given that there will likely be a few more material models added that use the Newton solver, perhaps it would be worth moving the finite difference section higher up as you suggested? Regardless, I'll take go through this today or tomorrow.

gassmoeller

In general I like the changes a lot, and think that they will be helpful for a number of other material models as well. However, I have a few technical issues with the changes. There are some pieces I do not quite understand and that need to be better documented. Moreover, I do not like to add checks that are only important, because of the way we test this material model. Is it possible to remove these get_simulator tests and still be able to test the model? Or would that be very complicated?

gassmoeller · 2018-01-30T05:55:34Z

source/material_model/visco_plastic.cc

@@ -182,7 +183,7 @@ namespace aspect
      // The first time this function is called (first iteration of first time step)
      // a specified "reference" strain rate is used as the returned value would
      // otherwise be zero.
-      const double edot_ii = ( (this->get_timestep_number() == 0 && strain_rate.norm() <= std::numeric_limits<double>::min())
+      const double edot_ii = ( (&(this->get_simulator()) != nullptr && this->get_timestep_number() == 0 && strain_rate.norm() <= std::numeric_limits<double>::min())


Why is this necessary? I was under the impression that we should never get here without a simulator object. Is this just something for the tests?

Yes, during the test of this material model, there is no simulator. So this is to make sure that the test actually works.

gassmoeller · 2018-01-30T05:58:29Z

source/material_model/visco_plastic.cc

+      MaterialModel::MaterialModelDerivatives<dim> *derivatives;
+      derivatives = out.template get_additional_output<MaterialModel::MaterialModelDerivatives<dim> >();
+
+      const double derivative_scaling_factor = &(this->get_simulator()) != nullptr ? (this->get_parameters().nonlinear_solver == Parameters<dim>::NonlinearSolver::Newton_Stokes ? this->get_newton_handler().get_newton_derivative_scaling_factor() : 0) : 1;


This is complicated to read, could you spell that out in actual if-else-loops. Again, why do we need to test for the simulator pointer?

If I remove the const I can spell it out in if-else statements. I could also add a comment explaining the line. Something like: If there is no simulator, it is a test environment and we want to enable the derivatives, so we set this variable to one. If we do not use the Newton solver, we do not need the derivatives, so we set this variable to zero. Otherwise, we use the value from the newton handler.

If having the const leads to code that is hard to read, then get rid of it. The comment sounds good to me!

gassmoeller · 2018-01-30T06:04:17Z

source/material_model/visco_plastic.cc

+              // compute derivatives if nessesary
+              if (derivative_scaling_factor != 0 && derivatives != NULL)
+                {
+                  const double finite_difference_accuracy = 1e-7;


This parts needs significantly more explanation. Why do you split up the components of the strain rate? I would prefer if we can move a lot of this into a separate function (either private in this material model or in aspect::Utilities). Is this currently only implemented for 2D (you only access elements 0 or 1)? What happens if strain_rate is close to zero?

Good point. I have moved the splitting of the components into a separate function, which allows for much less code :) It works now both in 2d and 3d. I also added more comments in the code.

gassmoeller · 2018-01-30T06:05:14Z

source/material_model/visco_plastic.cc

@@ -544,6 +666,14 @@ namespace aspect
    void
    ViscoPlastic<dim>::declare_parameters (ParameterHandler &prm)
    {
+      prm.enter_subsection("Compositional fields");


You do not need to declare a parameter again, this just needs to be done by the plugin owning this parameter.

So this, and the next comment are the same problem as above, the this->n_compositional_fields() is not available in when it is being tested.

gassmoeller · 2018-01-30T06:06:39Z

source/material_model/visco_plastic.cc

@@ -730,7 +860,13 @@ namespace aspect
    ViscoPlastic<dim>::parse_parameters (ParameterHandler &prm)
    {
      // increment by one for background:
-      const unsigned int n_fields = this->n_compositional_fields() + 1;
+      unsigned int n_fields = 0;


Why does the above not work? In theory you should have access to that variable already.

MFraters · 2018-02-06T23:39:07Z

It will be difficult to check the things I would like to check in a way where there is a simulator. We can of course discuss to see if there is a way to do it.

MFraters · 2018-02-10T00:04:31Z

@gassmoeller I guess I will have to make different tests for them after all, because I can't change the material model member variable Viscosity averaging scheme through creating a new parameter object and use simulator_access.get_material_model().parse_parameters(prm);, because get_material_model() returns a constant reference, so I get as an error:

error: passing ‘const aspect::MaterialModel::Interface<2>’ as ‘this’ argument discards qualifiers [-fpermissive]

Or do you know a way around this?

gassmoeller · 2018-02-12T16:28:59Z

Check by my office some time today, I am sure we can find a solution to this.

MFraters · 2018-02-12T23:11:30Z

@gassmoeller I think that I have addressed all your comments.

bangerth · 2018-02-14T00:00:00Z

@gassmoeller -- are you working this through with @MFraters or do you want me to also take a look?

gassmoeller

Much better, though I still have a number of comments. I also feel like the two loops over the compositions (separate for strain rate and pressure) in the evaluate function should be mergeable into a single one (with one inner loop for the strain rate components), but I can only see that clearly when my current comments are addressed. @bangerth: I have marked two questions for you, I can handle the rest.

gassmoeller · 2018-02-14T17:56:30Z

include/aspect/utilities.h

+     * Create matrix with unit at independent indices
+     */
+    template <int dim>
+    SymmetricTensor<2,dim> symmetric_independent_component_matrix (const unsigned int k);


I would name this function differently, though I am not sure if there is a convention for this. Maybe index_selector? It is something that has only zeroes as entries, except for a one in the kth index (or two 1/2 for the symmetric case). @bangerth: Do you have a good name for this?

The documentation could be the following: "A function that returns a SymmetricTensor, which entries are zero, except for the k'th component, which is set to one. If k is not on the main diagonal the resulting tensor is symmetrized."

For the text, it should be "...WHOSE entries are zero..."

Fundamentally, this function returns the elements of a basis for the set of symmetric tensors since every tensor can be written as a linear combination of the tensors returned by this function. So maybe something like nth_basis_for_symmetric_tensors()?

gassmoeller · 2018-02-14T17:59:17Z

source/material_model/visco_plastic.cc

+                  const double finite_difference_accuracy = 1e-7;
+
+                  // For each independent component, compute the derivative.
+                  for (unsigned int independent_component_k = 0; independent_component_k < SymmetricTensor<2,dim>::n_independent_components; independent_component_k++)


Lets rename the index to component, that makes it shorter, and the k is really not necessary if the name already explains what it is.

Also ++component instead of component++ for good style

gassmoeller · 2018-02-14T18:01:01Z

source/material_model/visco_plastic.cc

+                  for (unsigned int independent_component_k = 0; independent_component_k < SymmetricTensor<2,dim>::n_independent_components; independent_component_k++)
+                    {
+                      const TableIndices<2> indices_ij = SymmetricTensor<2,dim>::unrolled_to_component_indices (independent_component_k);
+                      SymmetricTensor<2,dim> strain_rate_component = strain_rate + std::max(std::fabs(strain_rate[indices_ij]), min_strain_rate) * finite_difference_accuracy


lets call this strain_rate_difference to stay consistent with the pressure_difference below.

also can this be const?

gassmoeller · 2018-02-14T18:05:47Z

source/material_model/visco_plastic.cc

+                      std::vector<double> eta_component = calculate_isostrain_viscosities(volume_fractions, pressure, temperature, composition, strain_rate_component,viscous_flow_law,yield_mechanism);
+
+                      // For each composition of the independent component, compute the derivative.
+                      for (unsigned int composition_i = 0; composition_i < eta_component.size(); composition_i++)


call this composition_index and do ++composition_index

gassmoeller · 2018-02-14T18:09:01Z

source/material_model/visco_plastic.cc

+                          if (derivative_component != 0)
+                            {
+                              // when the difference is non-zero, divide by the difference.
+                              derivative_component /= std::max(std::fabs(strain_rate_component[indices_ij]), min_strain_rate) * finite_difference_accuracy;


How is the term on the right hand side different from what you add in line 431? If it is the same: Give it a name and define it in line 430.

one uses the strain-rate and the other one uses the strain-rate plus the finite difference of the strain-rate.

gassmoeller · 2018-02-14T18:11:40Z

source/material_model/visco_plastic.cc

+                  // For each independent component, compute the derivative.
+                  for (unsigned int independent_component_k = 0; independent_component_k < SymmetricTensor<2,dim>::n_independent_components; independent_component_k++)
+                    {
+                      const TableIndices<2> indices_ij = SymmetricTensor<2,dim>::unrolled_to_component_indices (independent_component_k);


rename this to strain_rate_indices

gassmoeller · 2018-02-14T18:11:56Z

source/material_model/visco_plastic.cc

+                  /**
+                   * Now compute the derivative of the viscosity to the pressure
+                   */
+                  double pressure_difference = in.pressure[i] + (std::fabs(in.pressure[i]) * finite_difference_accuracy);


gassmoeller · 2018-02-14T18:12:32Z

source/material_model/visco_plastic.cc

+                   */
+                  double pressure_difference = in.pressure[i] + (std::fabs(in.pressure[i]) * finite_difference_accuracy);
+
+                  std::vector<double> pressure_difference_eta = calculate_isostrain_viscosities(volume_fractions, pressure_difference, temperature, composition, strain_rate,viscous_flow_law,yield_mechanism);


viscosity_difference and const?

gassmoeller · 2018-02-14T18:12:55Z

source/material_model/visco_plastic.cc

+
+                  for (unsigned int composition_i = 0; composition_i < pressure_difference_eta.size(); composition_i++)
+                    {
+                      double deriv_pressure = pressure_difference_eta[composition_i] - composition_viscosities[composition_i];


rename to viscosity_derivative

gassmoeller · 2018-02-14T18:23:29Z

source/utilities.cc

+    symmetric_independent_component_matrix (const unsigned int k)
+    {
+      // somehow I need this, because otherwise it complains that I am passing 3 arguments and it only takes 2...
+      const bool boolian = k < SymmetricTensor<2,dim>::n_independent_components;


@bangerth: I have tested this and found the same strange behavior. The following code does not work with deal.II 8.5, although it should:

const unsigned int k=0; Assert(k< SymmetricTensor<2,dim>::n_independent_components,ExcMessage(""))

The compiler error message is:

/home/rengas/Software/Aspect-Versionen/aspect/source/material_model/visco_plastic.cc:815:80: error: macro "Assert" passed 3 arguments, but takes just 2 Assert(k< SymmetricTensor<2,dim>::n_independent_components,ExcMessage("")); ^

As discussed on the way to lunch, you'll have to put a set of (parentheses) around the condition or the preprocessor is going to interpret the comma in 2,dim as one that separates macro arguments.

MFraters · 2018-02-16T02:13:50Z

It is ready for an new round of reviews :)

bangerth

I have only a couple small comments. Please address and squash your patches. Hopefully your testcase of today then becomes simpler to look at as well!

bangerth · 2018-02-19T21:55:07Z

source/material_model/visco_plastic.cc

+      // otherwise set it to zero.
+      double derivative_scaling_factor = 0.0;
+      if (this->get_parameters().nonlinear_solver == Parameters<dim>::NonlinearSolver::Newton_Stokes)
+        derivative_scaling_factor = this->get_newton_handler().get_newton_derivative_scaling_factor();


Hm, why do you actually need this in the material model? Shouldn't the scaling factor be applied only at the one or two places in the assembly where we compute the Newton matrix? If that were the case, you wouldn't have to deal with it in every single material model.

So the factor itself is actually not used to compute something here. It is only used to determine wheter the derivatives should be computed. Because when this scaling factor is zero, which may happen even when derivatives != NULL, there is no use in computing the derivatives.

bangerth · 2018-02-19T21:56:17Z

source/material_model/visco_plastic.cc

@@ -396,14 +408,83 @@ namespace aspect
              // TODO: This is only consistent with viscosity averaging if the arithmetic averaging
              // scheme is chosen. It would be useful to have a function to calculate isostress viscosities.
              const std::vector<double> composition_viscosities =
-                calculate_isostrain_viscosities(volume_fractions, pressure, temperature, composition, in.strain_rate[i],viscous_flow_law,yield_mechanism);
+                calculate_isostrain_viscosities(volume_fractions, pressure, temperature, composition, strain_rate,viscous_flow_law,yield_mechanism);


while you're here, break the line into shorter ones around 80 characters of length, and put a space after every comma separating the arguments

bangerth · 2018-02-19T21:57:50Z

source/material_model/visco_plastic.cc

+                  /**
+                   * Now compute the derivative of the viscosity to the pressure
+                   */
+                  const double pressure_difference = in.pressure[i] + (std::fabs(in.pressure[i]) * finite_difference_accuracy);


there are several long lines here. can you break them as well?

bangerth · 2018-02-19T21:59:57Z

source/utilities.cc

+    {
+      // somehow I need this, because otherwise it complains that I am passing 3 arguments and it only takes 2...
+      const bool boolian = k < SymmetricTensor<2,dim>::n_independent_components;
+      Assert(boolian, ExcMessage("The component is larger then the amount of independent components in the matrix.") );


You can write this as

Assert((k < SymmetricTensor<2,dim>::n_independent_components), ExcMessage("The component is larger then the amount of independent components in the matrix.") );

This has the advantage that the violated condition appears in the screen output.

(Unrelated: The correct spelling is "boolean" because it is named after Mr. Boole.)

bangerth · 2018-02-20T13:56:59Z

About the derivative scaling factor, as I mentioned yesterday my preference would be if the assembly loop simply didn't ask about information it doesn't need, rather than every material model separately trying to figure out whether the information it is being asked to deliver is actually used.

bangerth · 2018-02-20T17:11:50Z

We convinced ourselves that the code is correct as is. So I'll merge.

gassmoeller reviewed Jan 30, 2018

View reviewed changes

gassmoeller mentioned this pull request Feb 9, 2018

Simplify derivative test MFraters/aspect#16

Merged

MFraters force-pushed the add_derivatives_to_visco_plastic_material_model branch 2 times, most recently from a71f80c to e657541 Compare February 12, 2018 22:37

gassmoeller reviewed Feb 14, 2018

View reviewed changes

MFraters mentioned this pull request Feb 19, 2018

Fix some small Newton related problems #2100

Merged

bangerth reviewed Feb 19, 2018

View reviewed changes

MFraters force-pushed the add_derivatives_to_visco_plastic_material_model branch 2 times, most recently from f4866c7 to 4539f3c Compare February 20, 2018 00:01

MFraters and others added 7 commits February 20, 2018 01:35

Add derivatives to visco-plastic material model

9c9e5f1

Add tests for derivatives of visco-plastic material model derivatives.

63c7169

add symmetric_independent_component_matrix function to utilities.

5b41992

Simplify derivative test

1d38e6f

Let the test test everything in on go after the change.

48504b1

rename variables in visco_plastic.cc

79478b3

rename function in utilities.

56c55b6

MFraters force-pushed the add_derivatives_to_visco_plastic_material_model branch from 4539f3c to 56c55b6 Compare February 20, 2018 00:38

bangerth merged commit 1d912a2 into geodynamics:master Feb 20, 2018

Add derivatives to visco plastic material model #2030

Add derivatives to visco plastic material model #2030

Conversation

MFraters commented Dec 13, 2017

naliboff commented Dec 18, 2017

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

MFraters commented Feb 6, 2018

MFraters commented Feb 10, 2018

gassmoeller commented Feb 12, 2018

MFraters commented Feb 12, 2018

bangerth commented Feb 14, 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

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

MFraters commented Feb 16, 2018

bangerth 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

bangerth commented Feb 20, 2018

bangerth commented Feb 20, 2018