Implementation of plastic dilation #6373

YiminJin · 2025-06-12T02:21:03Z

This PR is an implementation of plastic dilation in the framework of ASPECT. Currently, plastic dilation is added to the Drucker Prager rheology, but is only available in the ViscoPlastic material model. It solves the oscillation problem shown in #5953.

I have tested the Picard iteration method, the defect correction method (DCM) and Newton method with the strip-footing experiment, and they all gave the correct result:

The convergence rates of different nonlinear solvers are as expected: Newton > DCM > Picard (either with AMG or with GMG).

VEP model with dilation angle ($\psi$) identical to the friction angle ($\phi$) has two distinctive features:

The shear bands develop along the Coulomb angle;
The width of shear bands are independent on the grid size.

The features can be observed in the numerical results of Kaus' (2010) benchmark (without viscoplastic damper):

The reason of feature 1 has been widely discussed and verified, while the reason of feature 2 is not so clear (the theory of associated/non-associated plastic flow can explain it, but I do not fully understand the theory).

I am sorry for modifying over 1000 lines of codes in a single PR. Please let me know if you find something confusing in my codes.

tjhei

I don't understand the reasoning for splitting dilation into LHS and RHS and the change to the -2/3 term for incompressible simulations, but the rest of the code looks clean and correct!
@naliboff can you take a look and then we can maybe discuss together?

source/simulator/assemblers/newton_stokes.cc

tjhei · 2025-06-12T16:06:30Z

source/simulator/assemblers/newton_stokes.cc

-
-              // Only assemble this term if we are running incompressible, otherwise this term
-              // is already included on the LHS of the equation.
-              if (enable_prescribed_dilation && !material_model_is_compressible)


why is this no longer needed?

The boolean material_model_is_compressible determines whether the assembler StokesCompressibleStrainRateViscosityTerm is executed. I have forced the assembler to be executed when plastic dilation is enabled (see line 76 in source/simulator/assembly.cc and line 55 in source/simulator/newton.cc), so the volumetric part of strain rate is subtracted from the LHS, and this term is no longer needed.

tjhei · 2025-06-12T16:09:55Z

source/simulator/solver/matrix_free_operators.cc

+            if (cell_data->is_compressible ||
+                cell_data->enable_prescribed_dilation)
              for (unsigned int d=0; d<dim; ++d)
                velocity_terms[d][d] -= viscosity_x_2 / 3. * div_u;


I did not see this change in the matrix-based version. Is this there as well?

Yes, the same change is made in the matrix-based assembler (line 76 in source/simulator/assembly.cc and line 55 in source/simulator/newton.cc).

tjhei · 2025-06-12T16:11:42Z

source/simulator/solver/stokes_matrix_free_local_smoothing.cc

-            this->get_newton_handler().parameters.newton_derivative_scaling_factor;
-
-          active_cell_data.enable_newton_derivatives = true;
+          active_cell_data.enable_newton_derivatives = (Parameters<dim>::is_defect_correction(this->get_parameters().nonlinear_solver)


With your change we now create the table for newtion_derivatives if dilation=true but we are not using a Newton solver. Is this necessary?

I think the Newton derivative tables are not created if the condition on the RHS of line 433 is not fulfilled. I have added if (active_cell_data.enable_newton_derivatives) before reinitializing those tables (line 461).

tjhei · 2025-06-12T16:14:33Z

Finally, it would be good to have some simple tests that exercise the new code. In fact the tests are currently broken:

/__w/aspect/aspect/tests/prescribed_dilation.cc:211:40: error: 'class aspect::MaterialModel::PrescribedPlasticDilation<2>' has no member named 'dilation'

naliboff · 2025-06-12T21:06:58Z

@YiminJin - Thank you for continuing to work on this and open the PR!

A few initial comments and thoughts from your initial post (more to come in the review):

I agree the total width of the shear band looks very similar across grid refinement levels, but the strain rate within the shear band is certainly changing as seen in the Brick Benchmark case. As such, I think the quasi-length scale introduced by the plastic dilation is different than the one introduced by the plastic damper term (i.e., the strain rate should not change once the damper length scale is resolved).
Following on point 1, does the strain-rate within the strip-footing experiment change at all as a function of resolution?

naliboff

@YiminJin - Nice work on this PR, it is in very good shape in terms of the structure, readability, and formatting. I confess I am having a very hard time evaluating the correctness of specific details, without seeing a higher-level summary of the equations and approach. Maybe we can discuss this in person? In any case, thank you again for pushing this forward and this looks like great progress!

include/aspect/material_model/interface.h

YiminJin · 2025-06-12T22:36:42Z

As such, I think the quasi-length scale introduced by the plastic dilation is different than the one introduced by the plastic damper term

I think I misunderstood the theory about associated/non-associated plastic flow. I reviewed the theory and it did not say that associated plastic flow is essentially mesh-independent. On the contrary, since we use a local method for plastic yielding (not global methods such as the strain gradient method), the thickness of the shear bands should be mesh-dependent. I will try to figure out why the shear bands behave like this in my results.

YiminJin · 2025-06-13T02:16:02Z

Following on point 1, does the strain-rate within the strip-footing experiment change at all as a function of resolution?

The strain rate field of strip-footing experiment does not change with refinement level. Here is the comparison:

bobmyhill · 2025-06-13T14:55:32Z

@YiminJin Nice work! I won't add another review, but I wanted to ask something. When you say "they all gave the correct result", do you mean the correct result, or just the same result? In #5953 you suggested that this was in fact the incorrect result (or at least different from your elASPECT code).

Whatever your answer, I don't think this should impact the merging of this PR. It would just be good to know what will still need fixing after this PR is merged (and why, if that's become clearer since last year).

YiminJin · 2025-06-14T05:17:06Z

I have found a simple way to change the 2D models in ASPECT to plane-strain models. After the modification, the results now look like this:

However, the Newton solver has some bugs currently, and the GMG version has not been implemented, so the new modifications are not pushed. Besides, my modifications conflict with some other PRs that use the orginal PrescribedPlasticDilation, so I think I had better put the "standard" plastic dilation terms to somewhere else.

bobmyhill · 2025-06-14T14:19:24Z

@YiminJin thanks for your answer! It's brilliant to see that you have isolated the cause of the difference between elASPECT and ASPECT results.

Could you clarify the difference you have identified between "2D models in ASPECT" and "plane-strain models"? It is my understanding that 2D models in ASPECT are explicitly using the plane-strain approximation: https://aspect-documentation.readthedocs.io/en/latest/user/methods/basic-equations/2d-models.html

Would it be equivalent to say that one of the "plastic dilation" formulations implicitly assumes isotropic dilation? If so, it is interesting that the two approximations produce such different results, and we should discuss as a larger team what we want in our final implementation.

YiminJin · 2025-06-14T16:00:33Z

Could you clarify the difference you have identified between "2D models in ASPECT" and "plane-strain models"? It is my understanding that 2D models in ASPECT are explicitly using the plane-strain approximation

There are two functions in ASPECT that are inconsistent with plane strain:

Function deviator(), which is frequently used in material models. It is defined as $\text{dev}(A) := A - \text{trace}(A) / \text{dim}$, while in plane-strain assumption it should always be $\text{dev}(A) = A - \text{trace}(A) / 3$.
Function second_invariant(), which is defined as $-A_{ij}A_{ij}/2$, but in plane strain case the contribution of $A_{33}$ should also be included.

After replacing the two functions by the plane strain version, the shear bands become sharp.

Would it be equivalent to say that one of the "plastic dilation" formulations implicitly assumes isotropic dilation?

I do not understand this question. Could you specify the "plastic dilation formulation" and the definition of "isotropic dilation" ? (Is there anisotropic dilation?)

we should discuss as a larger team what we want in our final implementation.

Sure! I hope we could make progress together.

bobmyhill · 2025-06-14T17:45:30Z

There are two functions in ASPECT that are inconsistent with plane strain

Can you raise a GitHub Issue describing the problem? I agree with your assessment, but it would be good to ask the other developers (a) whether they see something we don't, and (b) how they think we should address the problem given that any change would be a breaking change.

Could you specify the "plastic dilation formulation" and the definition of "isotropic dilation" ? (Is there anisotropic dilation?)

I think we can ignore my original question for the purposes of this PR (because I agree with your plane strain comment). But to answer your question, in plane strain any dilation must be anisotropic, because $\varepsilon_z=0$. Locally, shear bands are pretty much planar, and so any dilation should be uniaxial to first order. In terms of implementation, I don't think this matters, but before you pointed out the deviator() and second_invariant() functions I was reaching for alternative explanations.

gassmoeller

I will have to continue this review later, but just a few comments for now.

include/aspect/material_model/interface.h

gassmoeller · 2025-06-14T22:22:22Z

include/aspect/material_model/interface.h

+     * $\bar\alpha\gamma$ should be added to the right-hand side of the
+     * mass conservation equation, where $\bar\alpha$ is the negative
+     * derivative of plastic potential with respect to the pressure
+     * ($\sin\psi$ in 2D case), and $\gamma$ is the plastic multiplier.
+     * The plastic multiplier is given by
+     * $\gamma = (\tau_{II} - \alpha p - k) / \eta^{ve}$,
+     * where $\tau_{II}$ is the second invariant of the deviatoric stress,
+     * $\alpha$ is the negative derivative of yield function with respect to
+     * the pressure ($\sin\phi$ in 2D case), $k$ is cohesion, and $\eta^{ve}$,
+     * is the pre-yielding viscosity. When the Picard method or Defect Correction
+     * Method is applied, the term $\bar\alpha\gamma$ should be split into two
+     * terms:
+     * $\bar\alpha\gamma = \bar\alpha\alpha p / \eta^{ve} +
+     * \bar\alpha(\tau_{II} - k) / \eta^{ve}$,
+     * the former of which should be moved to the left-hand side in order to
+     * guarantee the stability of the nonlinear solver. Therefore, this output
+     * provides two terms: dilation_lhs_term corresponds to
+     * $\bar\alpha\alpha / \eta^{ve}$ (p is replaced by the shape function),
+     * and dilation_rhs_term corresponds to $(\tau_{II} - k) / \eta^{ve}$.


All of this description is pretty specific to the case of plastic dilation, but not the general case of someone prescribing delation because of other processes. Maybe it would be better to limit this documentation to the discussion of where the individual terms will go into the equation, and move the specific parts to drucker_prager.h:145

gassmoeller · 2025-06-14T22:23:42Z

include/aspect/material_model/rheology/drucker_prager.h

+          compute_dilation_terms_for_stokes_system(const double angle_friction,
+                                                   const double angle_dilation,
+                                                   const double cohesion,
+                                                   const double non_yielding_viscosity,
+                                                   const double effective_strain_rate) const;


could you use the following interface instead, this is more in line with that I did in #6384.

Suggested change

compute_dilation_terms_for_stokes_system(const double angle_friction,

const double angle_dilation,

const double cohesion,

const double non_yielding_viscosity,

const double effective_strain_rate) const;

compute_dilation_terms_for_stokes_system(const double angle_friction,

const double angle_dilation,

const double cohesion,

const double non_yielding_viscosity,

const double effective_strain_rate) const;

gassmoeller · 2025-06-14T22:25:23Z

include/aspect/material_model/rheology/visco_plastic.h

+      /**
+       * The LHS term corresponding to plastic dilation in the
+       * Stokes system. For details, see the comments of
+       * MaterialModel::PrescribedDilation::dilation_lhs_term.
+       */
+      std::vector<double> dilation_lhs_terms;
+
+      /**
+       * The RHS term corresponding to plastic dilation in the
+       * Stokes system. For details, see the comments of
+       * MaterialModel::PrescribedDilation::dilation_rhs_term.
+       */
+      std::vector<double> dilation_rhs_terms;


this code part looks different since #6384. Please rebase, the new version will be simpler, because it includes the DruckerPragerParameters structure already, which will automatically include the dilation angle.

gassmoeller · 2025-06-14T22:29:17Z

source/material_model/rheology/drucker_prager.cc

            // Average among phases
            drucker_prager_parameters.angle_internal_friction = MaterialModel::MaterialUtilities::phase_average_value(phase_function_values, n_phase_transitions_per_composition,
                                                                angles_internal_friction, composition);
+            drucker_prager_parameters.angle_dilation = angles_dilation[composition];


wouldnt you also need to use a call to MaterialModel::MaterialUtilities::phase_average_value to allow different dilation angles for different phases?

I forgot that. It is corrected in the newest commit.

source/material_model/rheology/drucker_prager.cc

gassmoeller · 2025-06-14T22:56:22Z

source/material_model/rheology/visco_plastic.cc

+            // Compute the dilation terms if necessary.
+            if (this->get_parameters().enable_prescribed_dilation == true)
+              {
+                if (output_parameters.composition_yielding[j] == true)
+                  {
+                    const double current_dilation = drucker_prager_parameters.angle_dilation * weakening_factors[1];
+                    const std::pair<double,double> dilation_terms = drucker_prager_plasticity.compute_dilation_terms_for_stokes_system (
+                                                                      current_friction,
+                                                                      current_dilation,
+                                                                      current_cohesion,
+                                                                      non_yielding_viscosity,
+                                                                      effective_edot_ii);
+                    output_parameters.dilation_lhs_terms[j] = dilation_terms.first;
+                    output_parameters.dilation_rhs_terms[j] = dilation_terms.second;
+                  }
+                else
+                  {
+                    output_parameters.dilation_lhs_terms[j] = 0;
+                    output_parameters.dilation_rhs_terms[j] = 0;
+                  }
+              }


we may need to think where to put this code compared to the code above. @naliboff any suggestions?

YiminJin · 2025-06-14T23:17:27Z

Can you raise a GitHub Issue describing the problem?

Sure. I have raised an issue (#6434 ) about this problem.

gassmoeller

Thanks for addressing my earlier comments. I think this code is already in pretty good shape and almost ready to merge. I left just a few comments concerned with optimizing assembly speed and a few requests for more documentation. Let me know when you have addressed these.

include/aspect/simulator/solver/matrix_free_operators.h

source/material_model/rheology/visco_plastic.cc

gassmoeller · 2025-06-16T20:32:54Z

source/material_model/rheology/visco_plastic.cc

+                                                                   volume_fractions,
+                                                                   dilation_values,
+                                                                   composition_dilation_derivatives_wrt_pressure,
+                                                                   1);


gassmoeller · 2025-06-16T20:37:38Z

source/simulator/assemblers/newton_stokes.cc

+                                                 ( one_over_eta +
+                                                   // add the dilation LHS term if plastic dilation
+                                                   // is enabled
+                                                   ( prescribed_dilation != nullptr ?
+                                                     prescribed_dilation->dilation_lhs_term[q] :
+                                                     0.0 )
+                                                 )


please do this computation up in line 146 where we compute one_over_eta. E.g.:
const double one_over_eta_plus_dilation = ...

The current loop is one of the most often executed lines in ASPECT, so it is better to move computations and if-branches outside of it.

gassmoeller · 2025-06-16T20:38:08Z

source/simulator/assemblers/newton_stokes.cc

                           // consider the derivatives deta/deps
                           // here, but we leave this as a TODO
-                           one_over_eta
+                           (one_over_eta - dilation_newton_factor)


yes, here you already precomputated dilation_newton_factor so this can stay as it is.

gassmoeller · 2025-06-16T20:40:29Z

source/simulator/assemblers/newton_stokes.cc

+                                       * (prescribed_dilation->dilation_rhs_term[q] -
+                                          prescribed_dilation->dilation_lhs_term[q] *


why is this the difference between rhs and lhs term? Since we are assembling the rhs here, I would expect to only see the rhs side term. Please add a comment that explains this above line 464.

gassmoeller · 2025-06-16T20:41:01Z

source/simulator/assemblers/newton_stokes.cc

+                                                - (prescribed_dilation == nullptr ? 0.0 :
+                                                   pressure_scaling * pressure_scaling *
+                                                   prescribed_dilation->dilation_lhs_term[q] *


precompute this term outside the loop over i and j.

gassmoeller · 2025-06-16T20:46:40Z

source/simulator/assemblers/stokes.cc

+                                                - (prescribed_dilation == nullptr ? 0.0 :
+                                                   pressure_scaling * pressure_scaling *
+                                                   prescribed_dilation->dilation_lhs_term[q] *


precompute this term outside the loop over i and j

gassmoeller · 2025-06-16T20:47:46Z

source/simulator/assemblers/stokes.cc

      Assert(!this->get_parameters().enable_prescribed_dilation
             ||
-             outputs.template get_additional_output_object<MaterialModel::PrescribedPlasticDilation<dim>>()->dilation.size()
+             outputs.template get_additional_output_object<MaterialModel::PrescribedPlasticDilation<dim>>()->dilation_lhs_term.size()


do you need to also add an assert for the dilation_rhs_term here? just in case someone accidentally sets them to different sizes.

YiminJin · 2025-06-16T22:05:35Z

@gassmoeller Thank you for the close examination on my codes! The problems you mentioned are addressed now.

gassmoeller

Thank you! This looks good to go to me, but I will leave it open for a day to give @naliboff, @bangerth, and @tjhei the chance to comment if they want to.

Can you also add a changelog entry that describes the changes that users will have to do if they wrote their own dilation terms (i.e. that the variable has been renamed and there are now two terms, lhs and rhs).

bangerth

I don't have enough of an understanding of the material models any more to be useful here. I'll let the experts be in charge for this PR :-)

gassmoeller · 2025-06-17T15:38:37Z

Since this is blocking #6419 let's merge it.

tjhei approved these changes Jun 12, 2025

View reviewed changes

tjhei added the ready to test label Jun 12, 2025

naliboff reviewed Jun 12, 2025

View reviewed changes

include/aspect/material_model/interface.h Outdated Show resolved Hide resolved

gassmoeller mentioned this pull request Jun 14, 2025

[WIP] [Request of initial comments] Add direction dilation. #6419

Open

gassmoeller reviewed Jun 14, 2025

View reviewed changes

YiminJin mentioned this pull request Jun 14, 2025

Make the 2D formulations consistent with plane strain assumption #6434

Open

YiminJin force-pushed the plastic-dilation branch 2 times, most recently from 8fda76d to 35d7af0 Compare June 16, 2025 17:58

gassmoeller reviewed Jun 16, 2025

View reviewed changes

YiminJin force-pushed the plastic-dilation branch from 697f16a to a832cf4 Compare June 16, 2025 22:02

gassmoeller approved these changes Jun 16, 2025

View reviewed changes

bangerth reviewed Jun 17, 2025

View reviewed changes

YiminJin force-pushed the plastic-dilation branch from e0b7451 to b25c571 Compare June 17, 2025 00:28

Implementation of plastic dilation

f283ae8

YiminJin force-pushed the plastic-dilation branch from b25c571 to f283ae8 Compare June 17, 2025 01:13

gassmoeller merged commit 7bc9ac6 into geodynamics:main Jun 17, 2025
8 checks passed

This was referenced Jun 17, 2025

Make deviator() and second_invariant() consistent with plane strain assumption in 2D #6471

Open

Newton method for plastic dilation #5953

Closed

YiminJin mentioned this pull request Jun 21, 2025

Modify the plastic dilation terms to reduce noises in effective viscosity #6546

Merged

		* (prescribed_dilation->dilation_rhs_term[q] -
		prescribed_dilation->dilation_lhs_term[q] *

Implementation of plastic dilation #6373

Implementation of plastic dilation #6373

Uh oh!

Conversation

YiminJin commented Jun 12, 2025

Uh oh!

tjhei left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

tjhei commented Jun 12, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

naliboff commented Jun 12, 2025

Uh oh!

naliboff left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

YiminJin commented Jun 12, 2025

Uh oh!

YiminJin commented Jun 13, 2025

Uh oh!

bobmyhill commented Jun 13, 2025

Uh oh!

YiminJin commented Jun 14, 2025

Uh oh!

bobmyhill commented Jun 14, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

YiminJin commented Jun 14, 2025

Uh oh!

bobmyhill commented Jun 14, 2025

Uh oh!

gassmoeller left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

YiminJin commented Jun 14, 2025

Uh oh!

gassmoeller left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

tjhei commented Jun 12, 2025 •

edited

Loading

bobmyhill commented Jun 14, 2025 •

edited

Loading