Create a new MaterialModel constructor. #1537

pmbremner · 2017-05-09T17:34:49Z

This patch adds a new MaterialModel array constructor that both initializes
and populates the MaterialModel arrays. The old constructor only sizes
the arrays, then the arrays must be populated within the calling function.
This constructor streamlines how these array values are obtained by making
a single standard call to the constructor. This patch is now called in
nullspace.cc to test.

ian-r-rose

There are also a number of places where the code indentation seems off. Running astyle should fix that.
EDIT: nevermind, my browser rendering was acting up. Indentation seems fine.

I am very excited for this.

ian-r-rose · 2017-05-09T17:37:41Z

include/aspect/material_model/interface.h

@@ -180,6 +182,16 @@ namespace aspect
      MaterialModelInputs(const unsigned int n_points,
                          const unsigned int n_comp);

+
+      //**************************************


This comment block should be of the same style as the other comment blocks (especially starting with /**). It is with this formatting that the documentation generation is able to grab your description.

ian-r-rose · 2017-05-09T17:41:40Z

source/material_model/interface.cc

+      velocity(fe.n_quadrature_points, numbers::signaling_nan<Tensor<1,dim> >()),
+      composition(fe.n_quadrature_points, std::vector<double>(introspection.n_compositional_fields, numbers::signaling_nan<double>())),
+      strain_rate(fe.n_quadrature_points, numbers::signaling_nan<SymmetricTensor<2,dim> >()),
+      cell (NULL)


This should be initialized. I think that fe.get_cell() would work.

I get an error with fe.get_cell() . Trying to troubleshoot

What is the error?

I think I know and know how to work around, but would need to see the error first.

/home/paulbremner/install/ASPECT/aspect/devel/source/material_model/interface.cc:270:7: error: class ‘aspect::MaterialModel::MaterialModelInputs’ does not have any field named ‘fe’
fe.get_cell()
^
/home/paulbremner/install/ASPECT/aspect/devel/source/material_model/interface.cc:270:9: error: expected ‘(’ before ‘.’ token
fe.get_cell()
^
/home/paulbremner/install/ASPECT/aspect/devel/source/material_model/interface.cc:270:9: error: expected ‘{’ before ‘.’ token
/home/paulbremner/install/ASPECT/aspect/devel/source/material_model/interface.cc: At global scope:
/home/paulbremner/install/ASPECT/aspect/devel/source/material_model/interface.cc:270:9: error: expected unqualified-id before ‘.’ token
/home/paulbremner/install/ASPECT/aspect/devel/source/material_model/interface.cc:760:1: error: expected ‘}’ at end of input
}
^

That can't be it -- you have a function argument above that's called fe. Did you rename it?

I don't believe so. The line right above this one uses fe

So do you have any changes in your local version right now compared to the version in the pull request?

Yes, but only changes to address the comments from Ian.

The only code changes (as in not a comment change) are this one being discussed, and deleting the line discussed in Ian's next comment.

ian-r-rose · 2017-05-09T17:43:07Z

source/material_model/interface.cc

+      fe[introspection.extractors.pressure].get_function_values (solution_vector, this->pressure);
+      fe[introspection.extractors.pressure].get_function_gradients (solution_vector, this->pressure_gradient);
+
+      std::vector<SymmetricTensor<2,dim> > strain_rate (fe.n_quadrature_points);


Is this used?

That isn't used... deleting line

ian-r-rose · 2017-05-09T17:45:07Z

source/material_model/interface.cc

+          for (unsigned int c=0; c<introspection.n_compositional_fields; ++c)
+            this->composition[i][c] = composition_values[c][i];
+        }
+      //**********************************


Remove comment.

bangerth · 2017-05-09T20:13:32Z

/run-tests

bangerth · 2017-05-09T20:48:02Z

include/aspect/material_model/interface.h

+
+      /**
+       * Constructor. Both initialize and populate the various arrays of this
+       * structure with the finite element and introspection objects and the


it's not actually the "finite element", but the "FEValues" object

Done as suggested.

bangerth · 2017-05-09T20:48:12Z

include/aspect/material_model/interface.h

+      /**
+       * Constructor. Both initialize and populate the various arrays of this
+       * structure with the finite element and introspection objects and the
+       * the solution_vector.


the the -> the

bangerth · 2017-05-09T20:50:51Z

source/material_model/interface.cc

+      fe_values[introspection.extractors.velocities].get_function_values (solution_vector, this->velocity);
+      fe_values[introspection.extractors.pressure].get_function_values (solution_vector, this->pressure);
+      fe_values[introspection.extractors.pressure].get_function_gradients (solution_vector, this->pressure_gradient);
+      fe_values[introspection.extractors.velocities].get_function_symmetric_gradients (solution_vector,this->strain_rate);


the last line is going to be an interesting one to consider. I don't think that all potential calling sites for this function do populate the strain_rate variable.

We should watch this in the places where we replace existing code by your function to see whether we need to modify this function further. The same is probably true for the gradient of the pressure, which is also not necessary everywhere.

Is there a way to make this function smart enough to only populate the necessary arrays?

Only by passing a whole bunch of bool arguments that indicate what is or isn't necessary to do. I would suggest leaving it as is for the moment (but undoing the change in nullspace.cc) and then carefully looking at each place where we want to call this function whether we are doing more than necessary or not -- and if yes, think about what to do next.

We could add additional bool flags to the constructor. It would get fairly busy at that point. Is it strictly reasonable to assume that a material model does not require something like strain rates or pressure gradients? They were added to MaterialModelInputs for a reason.

I feel like when we exclude the strain rates, we are implicitly making the assumption that we don't require the viscosity in the outputs, and strain rate only has an effect on the viscosity. This is probably not true in general. For instance, there are some dilational materials that change density upon shearing.

I'm not clear what you mean. If the strain rate is not filled, then the calling place explicitly doesn't want the viscosity. And because the strain rate vector is filled with signaling nans or has size zero, its contents can also not be used by accident. So I can't see a case where a place forgets to fill the strain rate but then accidentally uses the computed viscosity anyway.

I am thinking of a case where some of the other material model outputs (such as density) may depend on strain rates. In which case those material models would fail if strain_rate is not passed in.

Oh, I misread thinking that you were talking about the viscosity depending on strain rate, not density. I agree. Rock for example is such a material: if you shear it enough, you get sand that takes more volume than the original rock.

Nothing that can be done about this right now -- this simply isn't supported. I seem to recall that the reason we did it that way was that it is expensive to (i) evaluate the strain rate tensors, and (ii) to evaluate the viscosities, so we avoided it when possible.

Yeah, this is an issue for a different PR :)

For now I left the function as is, but I reverted nullspace.cc back to it's original form and I now test the function call in source/simulator/lateral_averaging.cc instead.

bangerth · 2017-05-09T20:51:05Z

source/material_model/interface.cc

+      fe_values[introspection.extractors.pressure].get_function_gradients (solution_vector, this->pressure_gradient);
+      fe_values[introspection.extractors.velocities].get_function_symmetric_gradients (solution_vector,this->strain_rate);
+
+      //Vectors for evaluating the finite element solution


add a space: // vectors...

also maybe say "...the compositional field parts of the finite element solution"

bangerth · 2017-05-09T20:51:55Z

source/material_model/interface.cc

+        }
+
+
+    }


no need for these empty lines

bangerth · 2017-05-09T20:57:10Z

source/simulator/nullspace.cc

-                }
-              in.cell = &cell;
-
-              material_model->evaluate(in, out);


this change doesn't look right. you used to call material_model->evaluate() before, but now you no longer do.

In fact, this may not be the very best place to test your changes. We really only need anything other than the velocities in in if use_constant_density == false. You now do substantially more work here by populating all of the fields (and adding back the now missing call to evaluate).

Maybe the best way to deal with this would be to leave this one place as it was before (can also easily be done with git) and pick a different place where we initialize a MaterialModelInputs object to test your new code.

nullspace.cc has been reverted back to it's original form. I have tested this new call in source/simulator/lateral_averaging.cc instead. The test passed. I have updated the commit message accordingly.

This patch adds a new MaterialModel array constructor that both initializes and populates the MaterialModel arrays. The old constructor only sizes the arrays, then the arrays must be populated within the calling function. This constructor streamlines how these array values are obtained by making a single standard call to the constructor. This patch is now called in source/simulator/lateral_averaging.cc to test.

bangerth · 2017-05-10T02:56:46Z

OK to merge if the tester is happy.

ian-r-rose reviewed May 9, 2017

View reviewed changes

pmbremner force-pushed the material_constructor branch from 004dbf5 to 9bb7b00 Compare May 9, 2017 20:43

bangerth reviewed May 9, 2017

View reviewed changes

pmbremner force-pushed the material_constructor branch from 9bb7b00 to cdcb02c Compare May 10, 2017 02:03

bangerth added the ready to merge label May 10, 2017

gassmoeller merged commit ca3e368 into geodynamics:master May 10, 2017

This was referenced May 10, 2017

Use new MaterialModelInputs constructor #1558

Closed

Add DataPostprocessor flavor of MaterialModelInputs constructor. #1722

Merged

Create a new MaterialModel constructor. #1537

Create a new MaterialModel constructor. #1537

Conversation

pmbremner commented May 9, 2017

ian-r-rose left a comment • edited

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

bangerth commented May 9, 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

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

bangerth commented May 10, 2017

ian-r-rose left a comment •

edited