add drucker prager rheology module #3257

naliboff · 2019-10-10T19:15:58Z

This PR adds a new model in the Rheology namespace for Drucker Prager plasticity, which will replace the similar functionality in MaterialUtilities.

This is currently work in progress and I wanted to let @anne-glerum take a look before proceeding further.

Remaining tasks and topics for discussion:

Use the new functions in the drucker_prager and viscoelastic_plastic material models.
Remove the related functionality in MaterialUtilities
Consistently use "yield stress" or "yield strength"
Rename various functions to be more intuitive/descriptive
Simplify existing documentation and make code more readable

I have followed the instructions for indenting my code.

For new features/models or changes of existing features:

I have tested my new feature locally to ensure it is correct.
I have created a testcase for the new feature/benchmark in the tests/ directory.
I have added a changelog entry in the doc/modules/changes directory that will inform other users of my change.

anne-glerum · 2019-10-18T13:52:51Z

Hi @naliboff , thanks for this! I like the new structure. Concerning your discussions points:

I'd go with yield_stress. This goes with most of the other function and variable names, like viscous_stress and compute_stress. Also, I associate stress sooner with the symbol sigma that we use in the comments than strength I think. But I'm also fine with strength :)
Instead of compute_stress I'd name the function compute_yield_stress.

Looking through it, I thought it might also be nice to have the DruckerPrager rheology store the updated cohesions and friction angles, instead of recomputing them to fill the AdditionalOutputs using the StrainDependent rheology, but maybe that's too much.

naliboff · 2019-10-21T19:12:14Z

@anne-glerum - Thanks for looking this over! I agree that yield_stress fits in better with the other function and variable names.

I like the idea of outputting the modified cohesion and friction angles through the DruckerPrager rheology. Eventually, we could change the visco plastic AdditionalOutputs variable from is_yielding to a specific deformation mechanism name (Plastic, Diffusion, Composite, etc).

naliboff · 2019-10-23T17:12:36Z

The test drucker_prager_derivatives is failing and it appears to be to related to the the parsing of parameters in the DruckerPrager rheology module (see errors below).

Looking at tests/drucker_prager_derivatives.cc, it's not clear to me what the issue is as the EquationOfState parameters are also parsed in the drucker_prager material model in a similar fashion.

Suggestions for what might be going on? I'll continue digging.

from GDB:

[ltt:12081] [ 1] /data1/Software/aspect/drucker_prager_rheology/aspect/debug/aspect(_ZNK6aspect15SimulatorAccessILi2EE22n_compositional_fieldsEv+0x4)[0x55798b907da2]
[ltt:12081] [ 2] /data1/Software/aspect/drucker_prager_rheology/aspect/debug/aspect(_ZN6aspect13MaterialModel8Rheology13DruckerPragerILi2EE16parse_parametersERN6dealii16ParameterHandlerE+0x32)[0x55798b4b5a02]
[ltt:12081] [ 3] /data1/Software/aspect/drucker_prager_rheology/aspect/debug/aspect(_ZN6aspect13MaterialModel13DruckerPragerILi2EE16parse_parametersERN6dealii16ParameterHandlerE+0xe3)[0x55798b41f0dd]

jdannberg · 2019-10-26T19:05:21Z

source/material_model/drucker_prager.cc

@@ -224,6 +224,7 @@ namespace aspect
                             "Units: $W/m/K$.");
          prm.enter_subsection ("Viscosity");
          {
+


@naliboff, I don't know if this causes the problem, but you want to remove the Angle of internal friction/Cohesion below (and same in parse parameters). These things are now read in the Rheology::DruckerPrager::declare_parameters(prm) function, right?

@jdannberg - Thanks for taking a look! I should have added a comment on Friday after figuring out the solution.

I actually intentionally left Angle of internal friction and Cohesion in drucker_prager as using the equivalent parameters in Rheology::DruckerPrager would have required adding a decent amount of new code and adjusting tests. This is due to the parameters in Rheology::DruckerPrager being lists, rather than single values.

@gassmoeller and I talked briefly on Friday about this and the issue is due to not having simulator access in the drucker_prager_derivatives test, which we need to get the number of compositional fields when parsing parameters. I checked with @MFraters and the correct path forward is to adjust the test (or make new tests) following the visco_plastic_derivatives_ tests, which have simulator access. This does indeed work and I'll submit the new tests in a separate PR.

Okay, but now that we have the new tests shouldnt this be removed? Or did I misunderstand our conversation? We do not want to have two parameters with nearly the same name, but different meanings in this material model.

naliboff · 2019-10-30T14:42:46Z

@MFraters - oddly, the 2D drucker_prager_derivatives test is failing, while the 3D one passes.

The changes occur at two points, one for pressure and for strain-rate:
https://jenkins.tjhei.info/job/aspect/job/PR-3257/5/artifact/changes-test-results.diff

The changes in strain-rate are very small, but the pressure less so. Is this something to look carefully into or is this kind of behavior consistent with what you saw previously? If the change is due to a bug, it would have to be in the compute_derivative function of the DruckerPrager rheology module, but I'm not seeing it.

gassmoeller

Not quite through yet, but a first set of comments. I am not quite happy with the current flow of information beween the DruckerPrager module and the material models. I suggested a different layout, but lets talk about it if that does not work for some reason. It would be nice to simplify all of this though, so I support the aim of this PR!

gassmoeller · 2019-10-30T18:48:21Z

doc/modules/changes/20191010_naliboff

+New: There is a new model in the MaterialModel::Rheology namespace for plastic
+material behavior based on the Drucker Prager yield criterion. The functionality
+is largely based on existing functions in MaterialUtilities and is currently used 
+in the visco plastic material. *Note, will updated once PR is approved and further


Dont forget to update this. What do you want to add?

gassmoeller · 2019-10-30T18:48:31Z

doc/modules/changes/20191010_naliboff

+in the visco plastic material. *Note, will updated once PR is approved and further
+changes are made. 
+<br>
+(John Naliboff, 2019-10-10)


This is not our usual date format

gassmoeller · 2019-10-30T18:49:52Z

include/aspect/material_model/drucker_prager.h

+        /*
+         * Objects for computing plastic stresses, viscosities, and additional outputs
+         */
+        Rheology::DruckerPrager<dim> drucker_prager_plasticity;


Shouldnt the other parameters below (angle_of_internal_friction, cohesion) now be removed, because they are inside drucker_prager_plasticity?

gassmoeller · 2019-10-30T18:50:14Z

include/aspect/material_model/dynamic_friction.h

+        /*
+          * Objects for computing plastic stresses, viscosities, and additional outputs
+          */
+        Rheology::DruckerPrager<dim> drucker_prager_plasticity;


Same here, remove other parameters that are no longer needed.

gassmoeller · 2019-10-30T18:51:22Z

include/aspect/material_model/rheology/drucker_prager.h

+           * Compute the plastic yield stress based on the Drucker Prager yield criterion.
+           */
+          double
+          compute_yield_stress (const double cohesion,


Shouldnt this function not get any inputs? The parameters are stored inside this class already.

Ah, I saw below why you are doing this. Still I think this is not a good solution. You pass back the parameters of this class, then modify them outside (using the strain weakening module) and then hand them over again to this function. This breaks the concept of encapsulation though, either these parameters are mostly used inside this class (then they shouldnt be modified outside of it), or they should not be members of this class. I would suggest the following interface instead:

double compute_yield_stress (const double pressure, const unsigned int composition_index, const double weakening_cohesion = 0.0, const double weakening_angle_of_friction = 0.0) const;

This way you do not need to pass the parameters out and in again, you just hand over all the information that is needed inside this function. You can modify the default value for the weakening parameters from 0.0 to 1.0 if that fits better with the formula.
Did I understand this correctly?

gassmoeller · 2019-10-30T18:53:19Z

include/aspect/material_model/rheology/drucker_prager.h

+
+
+          /**
+           * Return the value of the maximum yield stress for each composition used in


I do not see how this function can do for each composition. It gets no input, and always returns one double. It would either need to get an input (which composition?) or need to return a vector (for all compositions).

gassmoeller · 2019-10-30T18:56:53Z

source/material_model/drucker_prager.cc

@@ -224,6 +224,7 @@ namespace aspect
                             "Units: $W/m/K$.");
          prm.enter_subsection ("Viscosity");
          {
+


Okay, but now that we have the new tests shouldnt this be removed? Or did I misunderstand our conversation? We do not want to have two parameters with nearly the same name, but different meanings in this material model.

gassmoeller · 2019-10-30T18:58:24Z

source/material_model/rheology/drucker_prager.cc

+        // increment by one for background:
+        const unsigned int n_fields = this->n_compositional_fields() + 1;
+
+        angles_internal_friction = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Angles of internal friction"))),


All the functions in this class should use this parameter, not the ones handed over as parameters. Or did I misunderstand how this works?

gassmoeller · 2019-10-30T19:12:35Z

include/aspect/material_model/rheology/drucker_prager.h

+           * rheology model.
+           */
+          std::vector<double>
+          get_cohesions () const;


If the above comment works, we might not need this function anymore.

gassmoeller · 2019-10-30T19:12:42Z

include/aspect/material_model/rheology/drucker_prager.h

+           * composition used in the rheology model.
+           */
+          std::vector<double>
+          get_angles_internal_friction () const;


same for this one

naliboff · 2019-10-30T21:12:13Z

@gassmoeller - thanks for the review! A quick recap from our discussion:

I agree with the suggestion for improving the interface.
We can handle the duplicate parameters in the DruckerPrager material model through the declare_alias() function (link)
For now, we will hold off on moving the ViscoPlastic AdditionalOutputs to the DruckerPrager rheology module and revisit this when we add a second plasticity rheology module (e.g., dynamic friction, etc).

naliboff · 2019-11-04T20:42:21Z

@gassmoeller - It turns out there is an issue with the new proposed interface. In the dynamic_friction material model we calculate the steady-state friction coefficient, which is used to compute the Drucker Prager yield stress. Other material models in the future may want to calculate various permutations of the friction and cohesion angle as well (e.g., beyond simple strain weakening).

With this in mind, I don't think the new proposed interface is a good long-term solution and the best format is to keep with previous practice and have the material model pass the cohesion and internal friction value to the DruckerPrager rheology module.

In order to deal with the encapsulation issue, perhaps the most straightforward path is to move the Cohesions and Angles of internal Friction variables back to ViscoPlastic for now? As an alternative, we could also make a new rheology module for declaring these variables? For now, the first option seems like the most straightforward one.

naliboff · 2019-11-06T00:49:53Z

👍 - thanks for the change @gassmoeller, this is a good solution. I made the associated changed to drucker prager and dynamic friction, and merged those changes into your commit. I think I've now addressed all of your previous comments, nearly all of which were fixed by the last commit. So, likely ready for a final review.

@MFraters - any thoughts on why the drucker_prager_derivatives_2d test is failing, but the 3D case is not?

gassmoeller · 2019-11-06T17:05:22Z

The test was failing because the max_yield_stress for that test was changed. Previously the default was infity (no max yield stress), but the default for the input parameter is 1e12. I reverted that in another commit. I also removed the use of the DruckerPragerParameters in the drucker_prager material model and the dynamic_friction material model, because they both have a different structure for their input parameters. We can later think about how to unify them.

gassmoeller · 2019-11-06T17:05:51Z

When the tester agrees this is ready from my side.

naliboff · 2019-11-06T17:14:57Z

@gassmoeller - yikes, I'm surprised the change in maximum yield stress value affected only 1 out of 4 of the drucker prager tests. Thanks for pushing the other changes and agree that it is best to deal with unification of the interface and parameter declarations later on.

naliboff force-pushed the drucker_prager_rheology branch 2 times, most recently from 86fd259 to 3cec9aa Compare October 23, 2019 00:13

naliboff changed the title ~~[WIP] add drucker prager rheology module~~ add drucker prager rheology module Oct 23, 2019

jdannberg reviewed Oct 26, 2019

View reviewed changes

naliboff mentioned this pull request Oct 29, 2019

add tests #3275

Merged

4 tasks

naliboff and others added 5 commits October 29, 2019 20:11

add drucker prager rheology module

bb3fdf8

consistently use the phrase 'yield stress'

85f3478

remove duplicated parameters

ec5db6e

use drucker prager rheology module

91d462f

remove DruckerPrager material utilities

53d7487

naliboff force-pushed the drucker_prager_rheology branch from 3cec9aa to 53d7487 Compare October 30, 2019 03:41

gassmoeller reviewed Oct 30, 2019

View reviewed changes

naliboff force-pushed the drucker_prager_rheology branch from cc2367d to 4461b76 Compare November 5, 2019 23:50

Export drucker prager parameters

bd9feb4

naliboff force-pushed the drucker_prager_rheology branch from 4461b76 to bd9feb4 Compare November 5, 2019 23:52

update change log

d9cf56a

Simplify some code. Restore old behavior.

88fcf56

gassmoeller force-pushed the drucker_prager_rheology branch from 39ce754 to 88fcf56 Compare November 6, 2019 17:01

gassmoeller merged commit 89c11bf into geodynamics:master Nov 6, 2019

naliboff deleted the drucker_prager_rheology branch April 4, 2020 22:18

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

add drucker prager rheology module #3257

add drucker prager rheology module #3257

naliboff commented Oct 10, 2019

anne-glerum commented Oct 18, 2019

naliboff commented Oct 21, 2019

naliboff commented Oct 23, 2019

jdannberg Oct 26, 2019

naliboff Oct 26, 2019

gassmoeller Oct 30, 2019

naliboff commented Oct 30, 2019

gassmoeller left a comment

gassmoeller Oct 30, 2019

gassmoeller Oct 30, 2019

gassmoeller Oct 30, 2019

gassmoeller Oct 30, 2019

gassmoeller Oct 30, 2019

gassmoeller Oct 30, 2019

gassmoeller Oct 30, 2019

gassmoeller Oct 30, 2019

gassmoeller Oct 30, 2019

gassmoeller Oct 30, 2019

gassmoeller Oct 30, 2019

naliboff commented Oct 30, 2019

naliboff commented Nov 4, 2019

naliboff commented Nov 6, 2019

gassmoeller commented Nov 6, 2019

gassmoeller commented Nov 6, 2019

naliboff commented Nov 6, 2019



		/**
		* Return the value of the maximum yield stress for each composition used in

add drucker prager rheology module #3257

add drucker prager rheology module #3257

Conversation

naliboff commented Oct 10, 2019

For new features/models or changes of existing features:

anne-glerum commented Oct 18, 2019

naliboff commented Oct 21, 2019

naliboff commented Oct 23, 2019

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

naliboff commented Oct 30, 2019

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

naliboff commented Oct 30, 2019

naliboff commented Nov 4, 2019

naliboff commented Nov 6, 2019

gassmoeller commented Nov 6, 2019

gassmoeller commented Nov 6, 2019

naliboff commented Nov 6, 2019