How to plot diff/disl creep ratio using ASPECT

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.


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

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?

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.

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

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

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

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