How to plot diff/disl creep ratio using ASPECT

ASPECT
1

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
image1920Ă—1266 303 KB

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:

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:

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/2A_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

2

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 and then creating a new additional outputs option like PlasticAdditionalOutputs be a reasonable solution?

3

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.

4

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

5

Hi Ziqi,

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

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

6

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

1 Like
7

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.

    # 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?

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

    # \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

  # 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

8

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

9

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

10

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).

>>> 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
image
image

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

Cheers,
Ziqi

11

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) 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