# How to plot diff/disl creep ratio using ASPECT

**URL:** <https://community.geodynamics.org/t/how-to-plot-diff-disl-creep-ratio-using-aspect/3520>\
**Category:** ASPECT\
**Created:** [June 25, 2024, 8:01pm UTC](https://community.geodynamics.org/t/how-to-plot-diff-disl-creep-ratio-using-aspect/3520 "2024-06-25T20:01:07Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![maziqi96](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/maziqi96/32/1205_2.png) [@maziqi96](https://community.geodynamics.org/u/maziqi96)\
**Post date:** [June 25, 2024, 8:01pm UTC](https://community.geodynamics.org/t/how-to-plot-diff-disl-creep-ratio-using-aspect/3520/1 "2024-06-25T20:01:07Z")

</div>

Hi all,  
I have been working on plume-lithosphere interaction in 3d models. We are using the Visco plastic material model and we consider both diffusion and dislocation creep in our models. However, the output/postprocessor can only plot the total viscosity in Paraview. We want to plot a diff/disl creep figure in Paraview shown below in the upper panel.

 ![image](https://canada1.discourse-cdn.com/flex036/uploads/geodynamics/original/2X/2/24f9578bd0fbcce9c70537c52205086232b72b96.jpeg)  
[https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2019JB018560](https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2019JB018560)  
I have compared the Postprocess set up in this paper where the figure comes from with mine:  
**Setup in the 2019 JGR paper:**

```auto
subsection Postprocess
 set List of postprocessors = visualization, velocity statistics, basic statistics, depth average
.
.
.
  subsection Visualization
    set List of output variables = viscosity, density, strain rate, gamma # default

```

**Below is my setup:**

```auto
subsection Postprocess
  set List of postprocessors = visualization, temperature statistics, depth average, topography
.
.
.
  subsection Visualization
    set List of output variables = material properties, vertical heat flux, strain rate, stress, stress second invariant
    subsection Material properties
      set List of material properties = density, viscosity, thermal conductivity, thermal expansivity, specific heat, thermal diffusivity
    end

```

The difference is that the 2019 JGR paper output the strain rate. However, even if I output the strain rate, I am not sure how to plot the dislocation creep viscosity in Paraview. Do I plot the dislocation creep first with the following equation and then plot the diff/disl using (η-η\_disl)/η\_disl?  
η\_disl=1/2_A\_disl^(-1/n)_(ε ̇\_disl)^((1-n)/n)_exp((E\_disl_+pV\_disl\*)/(nRT))  
I would greatly appreciate it if someone could help me figure out a way to plot the diff/disl in Paraview. Thanks!

Best,  
Ziqi

---

<div class="post-metadata">

**Author:** ![jbnaliboff](https://avatars.discourse-cdn.com/v4/letter/j/58956e/32.png) [@jbnaliboff](https://community.geodynamics.org/u/jbnaliboff)\
**Post date:** [June 26, 2024, 4:12pm UTC](https://community.geodynamics.org/t/how-to-plot-diff-disl-creep-ratio-using-aspect/3520/2 "2024-06-26T16:12:04Z")

</div>

Hi @maziqi96,

Unfortunately, the visco plastic material model does not output the individual deformation mechanism (diffusion, dislocation, etc) viscosities. However, the material model could be modified to output them as additional outputs.

@bobmyhill @anne-glerum - would adding the creep viscosities as outputs to [IsostrainViscosities](https://github.com/geodynamics/aspect/blob/main/include/aspect/material_model/rheology/visco_plastic.h#L90) and then creating a new additional outputs option like [PlasticAdditionalOutputs](https://github.com/geodynamics/aspect/blob/main/include/aspect/material_model/rheology/visco_plastic.h#L53) be a reasonable solution?

---

<div class="post-metadata">

**Author:** ![bobmyhill](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/bobmyhill/32/38_2.png) [@bobmyhill](https://community.geodynamics.org/u/bobmyhill)\
**Post date:** [June 26, 2024, 5:16pm UTC](https://community.geodynamics.org/t/how-to-plot-diff-disl-creep-ratio-using-aspect/3520/3 "2024-06-26T17:16:02Z")

</div>

Hi @jbnaliboff,

Yes, I think that would be fine.

As an aside (and nothing to do as it’s in the code already), we didn’t name “IsostrainViscosities” very well. Yes, the viscosities are isostrain in the sense that they are calculated using the same (total) strain rate, but the combined viscosity implies that the elements experience isostress conditions (and this is what we want to approximate). True isostrain viscosities would be arithmetically added together (producing very different behaviour!).

I’ve thought a couple of times about how one could rationalise this choice, but I think we probably have to just settle for it being a good approximation.

---

<div class="post-metadata">

**Author:** ![maziqi96](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/maziqi96/32/1205_2.png) [@maziqi96](https://community.geodynamics.org/u/maziqi96)\
**Post date:** [July 3, 2024, 3:35pm UTC](https://community.geodynamics.org/t/how-to-plot-diff-disl-creep-ratio-using-aspect/3520/4 "2024-07-03T15:35:26Z")

</div>

Hi John and Bob,  
Thanks for your reply. However, I am sorry that I am a bit confused about how to achieve that. Does that mean I need to write a plugin.cc file with the dislocation/diffusion creep ratio as a new output under the structure IsostrainViscosities inside the class PlasticAdditionalOutputs? Can I modify the visco\_plastic.h file to what I want? Or do you recommend I write a new plugin.cc creating a new additional output from scratch?

Cheers,  
Ziqi

---

<div class="post-metadata">

**Author:** ![bobmyhill](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/bobmyhill/32/38_2.png) [@bobmyhill](https://community.geodynamics.org/u/bobmyhill)\
**Post date:** [July 5, 2024, 10:44am UTC](https://community.geodynamics.org/t/how-to-plot-diff-disl-creep-ratio-using-aspect/3520/5 "2024-07-05T10:44:06Z")

</div>

Hi Ziqi,

Sorry for not being clear. I think we / you will first have to decide what you want to output, because:

- in general, there are diffusion and dislocation viscosities for _each_ compositional field: ([aspect/source/material\_model/rheology/visco\_plastic.cc at e7ff10c15a8b8ef4f7397e1f72f71c6348ee3d95 · geodynamics/aspect · GitHub](https://github.com/geodynamics/aspect/blob/e7ff10c15a8b8ef4f7397e1f72f71c6348ee3d95/source/material_model/rheology/visco_plastic.cc#L138)). So do you want to output a vector of viscosity ratios?
- The concept of a viscosity ratio only makes any sense when the `viscous_flow_law` is set to `composite` ([aspect/source/material\_model/rheology/visco\_plastic.cc at e7ff10c15a8b8ef4f7397e1f72f71c6348ee3d95 · geodynamics/aspect · GitHub](https://github.com/geodynamics/aspect/blob/e7ff10c15a8b8ef4f7397e1f72f71c6348ee3d95/source/material_model/rheology/visco_plastic.cc#L215)). So what would you output if this option wasn’t chosen?
- Peierls creep ([aspect/source/material\_model/rheology/visco\_plastic.cc at e7ff10c15a8b8ef4f7397e1f72f71c6348ee3d95 · geodynamics/aspect · GitHub](https://github.com/geodynamics/aspect/blob/e7ff10c15a8b8ef4f7397e1f72f71c6348ee3d95/source/material_model/rheology/visco_plastic.cc#L232)), strain weakening ([aspect/source/material\_model/rheology/visco\_plastic.cc at e7ff10c15a8b8ef4f7397e1f72f71c6348ee3d95 · geodynamics/aspect · GitHub](https://github.com/geodynamics/aspect/blob/e7ff10c15a8b8ef4f7397e1f72f71c6348ee3d95/source/material_model/rheology/visco_plastic.cc#L246)) and elasticity ([aspect/source/material\_model/rheology/visco\_plastic.cc at e7ff10c15a8b8ef4f7397e1f72f71c6348ee3d95 · geodynamics/aspect · GitHub](https://github.com/geodynamics/aspect/blob/e7ff10c15a8b8ef4f7397e1f72f71c6348ee3d95/source/material_model/rheology/visco_plastic.cc#L257)) potentially affect the overall viscosity, but viscoplastic does not specify how these affect the diffusion/dislocation creep viscosity ratio. What is the sensible choice here?
- When a composition is yielding ([aspect/source/material\_model/rheology/visco\_plastic.cc at e7ff10c15a8b8ef4f7397e1f72f71c6348ee3d95 · geodynamics/aspect · GitHub](https://github.com/geodynamics/aspect/blob/e7ff10c15a8b8ef4f7397e1f72f71c6348ee3d95/source/material_model/rheology/visco_plastic.cc#L314)), neither the diffusion nor the dislocation viscosity contribute to the overall viscosity. What should be output then?

So yes, you could add the vector of dislocation/diffusion viscosity ratios to the structure IsostrainViscosities and add that vector to PlasticAdditionalOutputs. But I don’t know whether this is useful for you.

Best wishes,  
Bob

---

<div class="post-metadata">

**Author:** ![jbnaliboff](https://avatars.discourse-cdn.com/v4/letter/j/58956e/32.png) [@jbnaliboff](https://community.geodynamics.org/u/jbnaliboff)\
**Post date:** [July 12, 2024, 4:10pm UTC](https://community.geodynamics.org/t/how-to-plot-diff-disl-creep-ratio-using-aspect/3520/6 "2024-07-12T16:10:07Z")

</div>

Hi @maziqi96 and @bobmyhill,

@maziqi96 - To follow up on @bobmyhill last post, what we can do is add a new feature that outputs the viscosities of the different deformation mechanisms. It would be up to the user to then use these output to calculate ratios of interest, and decide what when such ratios are relevant (i.e., diff/disl ratio is not meaningful when yielding). If this is of interest, I am happy to add this feature, but it may take me a few weeks to get to it.

Cheers,  
John

---

<div class="post-metadata">

**Author:** ![maziqi96](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/maziqi96/32/1205_2.png) [@maziqi96](https://community.geodynamics.org/u/maziqi96)\
**Post date:** [July 12, 2024, 6:01pm UTC](https://community.geodynamics.org/t/how-to-plot-diff-disl-creep-ratio-using-aspect/3520/7 "2024-07-12T18:01:11Z")

</div>

Hi John and Bob,  
@jbnaliboff That would be great! Thank you very much! Currently, I am still trying to get the plugin to work. If I can’t, I plan to load the vtu file directly and post-process them to get the dislocation and diffusion creep separately.  
@bobmyhill Thanks for your followup. To answer your questions:

> do you want to output a vector of viscosity ratios?

yes

> The concept of a viscosity ratio only makes any sense when the viscous\_flow\_law is set to composite (aspect/source/material\_model/rheology/visco\_plastic.cc at e7ff10c15a8b8ef4f7397e1f72f71c6348ee3d95 · geodynamics/aspect · GitHub 1). So what would you output if this option wasn’t chosen?

The composite flow law is chosen in my models.

```auto
    # Select what type of viscosity law to use between diffusion, dislocation,
    # frank kamenetskii, and composite options. Soon there will be an option
    # to select a specific flow law for each assigned composition
    set Viscous flow law = composite

```

> Peierls creep (aspect/source/material\_model/rheology/visco\_plastic.cc at e7ff10c15a8b8ef4f7397e1f72f71c6348ee3d95 · geodynamics/aspect · GitHub), strain weakening (aspect/source/material\_model/rheology/visco\_plastic.cc at e7ff10c15a8b8ef4f7397e1f72f71c6348ee3d95 · geodynamics/aspect · GitHub) and elasticity (aspect/source/material\_model/rheology/visco\_plastic.cc at e7ff10c15a8b8ef4f7397e1f72f71c6348ee3d95 · geodynamics/aspect · GitHub) potentially affect the overall viscosity, but viscoplastic does not specify how these affect the diffusion/dislocation creep viscosity ratio. What is the sensible choice here?

Peierls creep, strain weakening and elasticity are all excluded in my models. So they won’t affect the overall viscosity in my models, right?

```auto
   # Whether to include Peierls creep in the rheological formulation.
    set Include Peierls creep = false

```

```auto
    # \item ``none'': No strain weakening is applied.
    # \item ``default'': The default option has the same behavior as ``none'',
    # but is there to make sure that the original parameters for specifying
    # the strain weakening mechanism (``Use plastic/viscous strain
    # weakening'') are still allowed, but to guarantee that one uses either
    # the old parameter names or the new ones, never both.
    set Strain weakening mechanism = default

```

```auto
  # Whether to include the additional elastic terms on the right-hand side of
  # the Stokes equation.
  set Enable elasticity = false

```

> When a composition is yielding (aspect/source/material\_model/rheology/visco\_plastic.cc at e7ff10c15a8b8ef4f7397e1f72f71c6348ee3d95 · geodynamics/aspect · GitHub), neither the diffusion nor the dislocation viscosity contribute to the overall viscosity. What should be output then?

There’s no yielding in my models.  
Hope this helps. Anyway, I appreciate your help!

Cheers,  
Ziqi

---

<div class="post-metadata">

**Author:** ![jbnaliboff](https://avatars.discourse-cdn.com/v4/letter/j/58956e/32.png) [@jbnaliboff](https://community.geodynamics.org/u/jbnaliboff)\
**Post date:** [July 17, 2024, 3:27pm UTC](https://community.geodynamics.org/t/how-to-plot-diff-disl-creep-ratio-using-aspect/3520/8 "2024-07-17T15:27:38Z")

</div>

> [@maziqi96](#):
>
> Peierls creep, strain weakening and elasticity are all excluded in my models. So they won’t affect the overall viscosity in my models, right?

That’s correct. It sounds like you are also working on a plugin to output the creep viscosities, or would you like me put in a PR to do this? The PR would not output the ratios (this is a bit to specific), but it will output the individual viscosities which can then be easily used for calculating the ratios in a postprocessing step.

Cheers,  
John

---

<div class="post-metadata">

**Author:** ![maziqi96](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/maziqi96/32/1205_2.png) [@maziqi96](https://community.geodynamics.org/u/maziqi96)\
**Post date:** [July 17, 2024, 3:55pm UTC](https://community.geodynamics.org/t/how-to-plot-diff-disl-creep-ratio-using-aspect/3520/9 "2024-07-17T15:55:12Z")

</div>

Hi John,  
Thanks for your reply. Sorry, I am a bit confused about what PR means. Does it mean plugin or Postprocessor? Outputing the individual viscosities sounds good to me! I can calculate the ratios in a postprocessing step in Paraview. Thank you very much!

Cheers,  
Ziqi

---

<div class="post-metadata">

**Author:** ![maziqi96](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/maziqi96/32/1205_2.png) [@maziqi96](https://community.geodynamics.org/u/maziqi96)\
**Post date:** [August 20, 2024, 3:37pm UTC](https://community.geodynamics.org/t/how-to-plot-diff-disl-creep-ratio-using-aspect/3520/10 "2024-08-20T15:37:05Z")

</div>

Hi John,  
Thanks for your help on this so far. I want to plot the diff/disl ratio before Ada Lovelace, which is near. The quickest approach I can think of right now is to calculate the diffusion and dislocation creep separately. Now I am trying to load the pvtu data into Python using PyVista. However, I am not sure how to rearrange the property (e.g. temperature, pressure, and strain rate) array according to x,y,z coordinates because I only want to plot the diff/disl ratio in the upper mantle (z= 340 km to 1000km) on the 2D cross-section (y=0 km).

```auto
>>> mesh.point_data 
pyvista DataSetAttributes
Association : POINT
Active Scalars : None
Active Vectors : None
Active Texture : None
Active Normals : None
Contains arrays :
    velocity float32 (20783574, 3)
    stress float32 (20783574, 9)
    p float32 (20783574,)
    T float32 (20783574,)
    density float32 (20783574,)
    viscosity float32 (20783574,)
    thermal_conductivity float32 (20783574,)
    thermal_expansivity float32 (20783574,)
    specific_heat float32 (20783574,)
    thermal_diffusivity float32 (20783574,)
    vertical_heat_flux float32 (20783574,)
    strain_rate float32 (20783574,)
    stress_second_invariant float32 (20783574,)
>>> points = mesh.points
>>> temperature = mesh['T']
>>> pressure = mesh['p']
>>> strain_rate = mesh['strain_rate']
>>> strain_rate
pyvista_ndarray([1.580637e-16, 1.580637e-16, 1.580637e-16, ...,
                 2.752588e-18, 2.752588e-18, 2.752588e-18], dtype=float32)

```

Another question is that I want to make sure I can use the following equations to calculate the diff/disl ratio correctly.  
 ![image](https://canada1.discourse-cdn.com/flex036/uploads/geodynamics/original/2X/c/c8e58c18d8ba7c2c2c4e9d4ddb1c0ef681e01e25.png)  
 ![image](https://canada1.discourse-cdn.com/flex036/uploads/geodynamics/original/2X/8/8050f516e8e39034de3d49af0b49c6162eb1dfed.png)  
 ![image](https://canada1.discourse-cdn.com/flex036/uploads/geodynamics/original/2X/d/d13d24651556c41d3f241888b39895ba0355937c.png)

I would greatly appreciate it if you could help me answer these questions.

Cheers,  
Ziqi

---

<div class="post-metadata">

**Author:** ![bobmyhill](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/bobmyhill/32/38_2.png) [@bobmyhill](https://community.geodynamics.org/u/bobmyhill)\
**Post date:** [August 20, 2024, 4:08pm UTC](https://community.geodynamics.org/t/how-to-plot-diff-disl-creep-ratio-using-aspect/3520/11 "2024-08-20T16:08:50Z")

</div>

Dear Ziqi,

> I am not sure how to rearrange the property (e.g. temperature, pressure, and strain rate) array according to x,y,z coordinates because I only want to plot the diff/disl ratio in the upper mantle (z= 340 km to 1000km) on the 2D cross-section (y=0 km).

`mesh.points` contains all the point positions. If you only want a subset of points, you can mask the ones you don’t want, e.g. `mask = np.abs(mesh.points[:,1]) < 1.e-5`, then `strain_rate[mask]`. Before doing this, you might want to convert from pyvista\_ndarray to numpy array ([pyvista.convert\_array — PyVista 0.44.1 documentation](https://docs.pyvista.org/api/utilities/_autosummary/pyvista.convert_array.html)) and then use numpy’s `reshape` method to get the arrays in the right shape for plotting using matplotlib. numpy’s online documentation should be sufficient to tell you how to do this.

> Another question is that I want to make sure I can use the following equations to calculate the diff/disl ratio correctly.

Yes, your equations are correct. You should bear in mind that ViscoPlastic does not decompose the strain rates, so \dot{\varepsilon} in your second expression should be the total strain rate to be consistent with what is done in the material model.

Best wishes,  
Bob