Strip pattern in viscosity when using visco plastic

Hi,

I am running a 2D subduction model using Visco plastic and I am using version 2.4.0-pre. The results I got are behaving weirdly. There are a few problems:

  1. There is a weird high viscosity at the bottom center. I limited the viscosity between 1e18 to 1e22. If I increase the limit, I will encounter a solver converging problem.

  2. Once the model starts to run, the viscosity field has this strip pattern. And the mantle viscosity is really high. I check temperature, pressure and strain rate profiles, they all look right. I don’t understand what could cause this problem.

This figure shows the beginning of the model (the first output) and the second output (0.5Ma).

  1. When I tried to change the dislocation creep parameters (e.g. activation volume) or thermal parameters (e.g. thermal expansion), or change refinement, I will run into a solver converging problem and the error shows a high residual which means my matrix has a null-space.

Does anyone have any idea what could cause this problem? What should I try to resolve this problem?

Thank you so much!

Xiaowen

Hi Xiaowen,

A few comments on your model, as we encountered similar problems.

  1. on the high viscosity at the bottom: are you using dislocation creep only (i.e. without diffusion creep)? Because if so, then it might be a stagnation point, where viscosity is high because strainrate is low. Adding diffusion creep (i.e. make rheology composite) should resolve this, as diffusion creep is insensitive to strain rate.
  2. We had these similar ‘stripy’ patterns, and the underlying problem, in our case, was the composition field that shows these stripy patterns, and on which the viscosity is dependent. Changing to the SUPG advection scheme helped there, but to get completely rid of it, I suggest to use particles to advect composition, instead of a field method.
  3. There are lots of different solver settings that may be tried, but I am not the best person to comment on this, so leave this to others to comment on.

Hope this helps.
Best wishes,
Jeroen

Hi Xiaowen ,

Thank you for the detailed post and accompanying image summarizing the issue. I’ll build on the excellent questions and tips from Jeroen below (thanks Jeroen!).

  1. Same thought and question as Jeroen, are you only using dislocation creep? Also, what is the basal boundary condition (free-slip or no-slip)?

  2. I second switching from fields to particles. Let us know if you need input on how to do this. Question - does increasing or decreasing the resolution changing the characteristics (magnitude, timing) of the “stripy patterns” in the viscosity or underlying compositional field?

  3. Indeed, there are many solver options to choose from, but I think for this problem the best nonlinear solver scheme to start with is 'iterated Advection and defect correction Stokes" (+ block GMG for the linear solver). Let me know if it would be helpful to post an example prm snippet using these schemes.

Cheers,
John

Hi Jeroen and John,

Thank you so much for both of your replies.

  1. Yes, I am using a dislocation creep only. I thought to limit the minimum strain rate could prevent this happening, I guess it is still too low. I am now using composite rheology (adding diffusion creep). The mantle rheology looks correct now. However, I only got the first output. The model takes a quite a long time with a time step of less than 100yr (CFL=0.5). The bottom boundary is a free slip, and the side boundaries have prescribed velocity and fixed composition (inflow from the top part of the side boundaries and balanced outflow is through the bottom part of the side boundaries). Top is a free surface. I think I will need to tune my solver parameters.
  2. I understand that using particles would solve this problem. I am just thinking for future converting into a 3D model, using particles would be computationally expensive. I will try Jeroen’s suggestion about using SUPG advection scheme.
  3. I am using single advection and iterated defect correction stokes. John, if you could send me an example prm for 'iterated Advection and defect correction Stokes" (+ block GMG for the linear solver), that would be great. Thank you.

Thanks,
Xiaowen

Hi Xiaowen,

  1. A 100 yr time step is indeed quite small, but not out of the realm of conceivable values if you have very high local velocities. What is your smallest mesh size and highest velocity? Do you see anomalous velocities anywhere in the model?

  2. Good news about particles - there are now particle interpolation schemes (linear least squares and bilinear linear least squares) that when combined with a limiter, can achieve very accuracy with far fewer particles per cell than something like an arithmetic average interpolation scheme. I recommend looking through the test suite in the ASPECT master branch for examples. FYI, these have been introduced into ASPECT recently and we (developers) are still testing them in complex model setups, but the results are very encouraging.

  3. A solver scheme example is below. Note that in the Material model section, you will need to turn on averaging over the element (default is none). I suggest using harmonic averaging:
    set Material averaging = harmonic average only viscosity

Please let us know how things go with the model changes!

Cheers,
John

set Nonlinear solver scheme = iterated Advection and defect correction Stokes

Try increasing to a stricter value at some point

set Nonlinear solver tolerance = 1e-5

If it consistently takes more than 200 nonlinear iterations, adjustments are likely needed

set Max nonlinear iterations = 200

subsection Solver parameters
  subsection Stokes solver parameters
    set Stokes solver type = block GMG
    set Number of cheap Stokes solver steps = 1000
    set Use full A block as preconditioner  = true
    set Linear solver tolerance             = 1e-7 # test the effect of this value
  end
  subsection Newton solver parameters
    set Maximum linear Stokes solver tolerance = 1e-2 # test the effect of this value
    set Use Eisenstat Walker method for Picard iterations = true
  end
end

Hi John,

  1. At the beginning of the model run, the maximum velocity is 80 m/yr in the upper mantle beneath the subducting slab (flow goes around the hanging slab). The high velocity is at the first few step and decrease to 7m/yr which is still pretty high. My smallest mesh is 2.5km.
    After I use composite flow law, I keep running into solver converging problem. For my subduction model, there are velocities imposed on the side boundaries (7cm/yr) and my boundary conditions are fixed composition on the left and right boundary and the boundary temperature is fixed as the initial temperature. The bottom is free slip and the top is free surface. I am not sure if it is my boundary condition cause the problem.

I tried to use 'iterated Advection and defect correction Stokes" (+ block GMG). I attached my boundary conditions and the Material model section. Could you spot any obvious error?

subsection Boundary composition model
  set Fixed composition boundary indicators   = right, left
  set List of model names = initial composition
end

subsection Boundary temperature model
  set Fixed temperature boundary indicators   = bottom, top, left, right
  set List of model names =  initial temperature
end

subsection Mesh deformation
  set Mesh deformation boundary indicators        = top: free surface, top: diffusion
  subsection Free surface
    set Surface velocity projection = normal
  end
  subsection Diffusion
    set Hillslope transport coefficient = 1.e-7
  end
end

subsection Boundary velocity model
  set Tangential velocity boundary indicators = bottom
  set Prescribed velocity boundary indicators = right x:function, left x:function

  subsection Function
    set Function constants  = w=2000e3, cm=0.01, year=1, oc_v=5, cont_v=2, blance_out=0.5926 #=(5*80+2*120)/(560+520)   
    set Function expression = if (x<=w/2 && y>580e3, (oc_v)*cm/year, \
                               if (x<=w/2 && y<560e3, -(blance_out)*cm/year, \
                                 if (x<=w/2 && y<=580e3 && y>=560e3, \
                                 ((y-560e3)*(oc_v+blance_out)/20e3-blance_out)*cm/year, \
                                   if (x>w/2 && y>540e3, -cont_v*cm/year, \
                                     if (x>w/2 && y<520e3, blance_out*cm/year, \
                                         ((y-540e3)*(blance_out+cont_v)/(-20e3)-cont_v)*cm/year \
                                         )\
                                       )\
                                     )\
                                   )\
                                 );\
                              0
    set Variable names      = x,y
  end
end

Above are my boundary conditions

subsection Material model
  set Model name = visco plastic
  set Material averaging = harmonic average only viscosity

  subsection Visco Plastic
    
    set Minimum strain rate = 1.e-20
    set Reference strain rate = 1.e-16

    # Limit the viscosity with minimum and maximum values
    set Minimum viscosity = 1e18
    set Maximum viscosity = 1e25
    
    set Reference temperature = 273

    # compositional field: background, OC, OML, UC, LC, CML
    set Define thermal conductivities = true
    set Thermal conductivities        = 54.2, 2.25, 2.25, 2.25, 2.25, 2.25
    set Heat capacities       =        750.
    set Densities             =        3337.4,        2950,        3337.4,         2811.2,       2929,       3337.4
    set Thermal expansivities =        2e-5

    # Harmonic viscosity averaging
    set Viscosity averaging scheme = harmonic

    set Viscous flow law =  composite
    
    # Dislocation creep parameters
   # Materials                                   background,       OC,       OML,         UC,        LC,       CML
    set Prefactors for diffusion creep          =   1.5e-16,    1e-29,    1.5e-16,   1.5e-16,   1.5e-16,   1.5e-16
    set Activation energies for diffusion creep =    375.e3,     0.e3,     375.e3,    375.e3,    375.e3,    375.e3
    set Activation volumes for diffusion creep  =     3.e-6,       0.,      3.e-6,     3.e-6,     3.e-6,     3.e-6
    # Dislocation creep parameters
    # Note that the viscous pre-factors below are scaled to plane strain from unixial strain experiments.
    set Prefactors for dislocation creep          = 6.52e-16,    1e-50,  6.52e-16,  8.57e-28,  7.13e-18,  6.52e-16
    set Stress exponents for dislocation creep    =      3.5,      1.0,       3.5,       4.0,       3.0,       3.5
    set Activation energies for dislocation creep =   530.e3,       0.,    530.e3,    223.e3,    345.e3,    530.e3
    set Activation volumes for dislocation creep  =   18.e-6,       0.,    18.e-6,        0.,        0.,    18.e-6

    # Plasticity parameters
    set Angles of internal friction =    30,  0,    30,    30,    30,    30
    set Cohesions                   = 20.e6,  0, 20.e6, 20.e6, 20.e6, 20.e6
  end
end

Above is the material model section.
The temperature and compositional field look fine. I think it’s either the boundary conditions or the material model cause the problem.

  1. That is good to know about particles. I will try to implement it in the models.

Thanks,
Xiaowen

Hi Xiaowen,

Thank you for posting all of that information, this is really helpful. As these are very complicated models solving the various issues may take quite a bit more debugging and testing, but we will get there eventually.

Follow up questions, ideas and responses:

  1. Is this model setup based on a prior published experiment? If so, are there prior results to compare against?

  2. I can’t spot any obvious errors in your boundary conditions, but these things are unfortunately hard to debug without seeing the results. Do the imposed velocities on the boundaries from the model solution match what you intended to put in? For these types of checks, I often run a quick model with constant material properties (ie., fixed viscosity) just to make sure the boundary velocities are correct. The material properties look reasonable at first glance.

  3. A period of very fast velocities at the beginning of the model run is not entirely unexpected if you have lateral variations in composition/rheology (ex: ocean-continent boundary, slab nose pre-imposed in the upper mantle) in combination with a free a surface. Effectively, the system is not in isostatic equilibrium and is rapidly adjusting at the beginning of the model run. For your next step, I recommend setting the left and right boundaries to free-slip velocity BCs, and seeing if you still obtain these high velocities at the beginning of the model run.

  4. Side note on thermal conductivity - I see you are using a very high value in the asthenosphere, presumably to maintain an adiabatic gradient? As an alternative, you can use simplified adiabatic heating and shear heating with the extended BA approximation, and this should produce an adiabatic gradient.

Cheers,
John

Hi John,

Thank you for your suggestions and sorry for the late reply.
I have been trying different things and was able to get rid of the strip pattern. I find that using a maximum composition for viscosity averaging is able to get rid of the problem. Thank you for suggesting using adiabatic heating and shear heating. Using these instead of a fake high thermal conductivity in the mantle makes the model look better and it seems shorten the running time.

—I am now able to run the subduction model with prescribed velocity for over 10 Mry and the model look good. However, I run into a new problem. After running for a certain time, I got the error message “solve_advection did not converge”. It is because some points in the model have a nan viscosity, where there are surface topography high. Please see the figure below.

This problem won’t happen if I use a slower prescribed convergence rate. Do you have any suggestion that I can try to avoid this problem? I am using free surface with diffusion. I tried to play with the free surface parameters, but it doesn’t seem to help much.

–I have another question about one parameter in ASPECT. In the visco-plastic section, for the dislocation parameter, there is one called “constant viscosity prefactors (multiplicative factor)”. Is this a scaling factor that apply to the viscosity formula to linear scale viscosity?

Thank you.

Xiaowen

Hi Xiaowen,

Thank you for the update and great to hear the models are now working!

NaN viscosities - Interesting, I have not seen that before. Are the NaN viscosities being reported in Paraview, or within ASPECT (error message, debugger)? Numerical instabilities can certainly develop near the free surface when the mesh is highly deformed, but I think the topography and mesh deformation in the attached image is not necessarily at that point yet. Would you mind posting the full error message from the log.txt (or other) file?

What value of diffusivity are using for the free surface diffusion parameter. I think going up to a value of 1e-7 is within reasonable physical limits.

This problem won’t happen if I use a slower prescribed convergence rate. Do you have any suggestion that I can try to avoid this problem? I am using free surface with diffusion. I tried to play with the free surface parameters, but it doesn’t seem to help much.

For advection solver errors, I normally first try decreasing the time step (via CFL value or maximum time step).

You can also try increasing the GMRES restart length to 100 or 150 (default is 50), although this will require more memory:
subsection Advection solver parameters
# This is the number of iterations that define the GMRES solver restart
# length. Increasing this parameter makes the solver more robust and
# decreases the number of iterations. Be aware that increasing this number
# increases the memory usage of the advection solver, and makes individual
# iterations more expensive.
set GMRES solver restart length = 50
end

I have another question about one parameter in ASPECT. In the visco-plastic section, for the dislocation parameter, there is one called “constant viscosity prefactors (multiplicative factor)”. Is this a scaling factor that apply to the viscosity formula to linear scale viscosity?

Yes, this parameter will multiply the viscosity by the factor supplied in the parameter file. This happens right before the total viscous stresses are calculated. For example, if the constant viscosity prefactor is 0.1 and the viscosity determined by dislocation creep (or composite creep) is 1e24 Pa s, then the final “creep” viscosity will be 1e23 Pa s.

I hope this helps and please keep us posted!

Cheers,
John

Hi John,

ASPECT report the error as “the residual is -nan”. I checked on paraview, the last output there are points where viscosity is nan as I showed in the previous post. Here is the error message:

An error occurred in line <2778> of file </home/r/russ/xiaowen4/aspect/source/utilities.cc> in function
void aspect::Utilities::linear_solver_failed(const string&, const string&, const std::vectordealii::SolverControl&, const std::exception&, ompi_communicator_t* const&, const string&)
The violated condition was:
false
Additional information:
The iterative advection solver in Simulator::solve_advection did not
converge.
The initial residual was: -nan
The final residual is: -nan
The required residual for convergence is: 1.000000e-50
See Sub_1-2-1-7/solver_history.txt for the full convergence history.
The solver reported the following error:
--------------------------------------------------------
An error occurred in line <1075> of file
</scinet/niagara/software/2022a/opt/gcc-11.3.0-openmpi-4.1.4+ucx-1.11.2/dealii/9.4.0/deal.II-v9.4.0/include/deal.II/lac/solver_gmres.h>
in function
void dealii::SolverGMRES::solve(const MatrixType&,
VectorType&, const VectorType&, const PreconditionerType&) [with
MatrixType = dealii::TrilinosWrappers::SparseMatrix;
PreconditionerType = dealii::TrilinosWrappers::PreconditionILU;
VectorType = dealii::TrilinosWrappers::MPI::Vector]
The violated condition was:
iteration_state == SolverControl::success
Additional information:
Iterative method reported convergence failure in step 0. The residual
in the last step was -nan.
This error message can indicate that you have simply not allowed a
sufficiently large number of iterations for your iterative solver to
converge. This often happens when you increase the size of your
problem. In such cases, the last residual will likely still be very
small, and you can make the error go away by increasing the allowed
number of iterations when setting up the SolverControl object that
determines the maximal number of iterations you allow.
The other situation where this error may occur is when your matrix is
not invertible (e.g., your matrix has a null-space), or if you try to
apply the wrong solver to a matrix (e.g., using CG for a matrix that
is not symmetric or not positive definite). In these cases, the
residual in the last iteration is likely going to be large.

For the free surface, I did use normal projection and diffusion with coefficient of 1.e-7.

I tried as you suggested GMRES, I still got the same error message after ~12 Myr. I haven’t try change the CFL. My current one is

set CFL number = 0.3
set Maximum time step = 20e3

Thank you for you the reply.

Cheers,
Xiaowen