Hello ASPECT community!
I’m a grad student working on modeling convection in an icy shell, I have a question about a 2D visco plastic model with variable viscosity that my advisor and I have been working on. Specifically about how the viscosity is calculated. Referring to Section 2.2.1 of this paper, equations 5 and 14:
Equation 5 is the same form described in the ASPECT documentation for our visco plastic material model.
We mainly reference values from this work (https://doi.org/10.1029/2000JB900336) which models the deformation of ice to set the value for constants. In the subsection Material model, model name visco platic, we set the following parameters:
Minimum strain rate = 1.e-20
Reference strain rate = 1.e-16
Minimum viscosity = 1e13
Maximum viscosity = 1e14
Viscosity averaging scheme = harmonic
Viscous flow law = diffusion
A = Prefactors for diffusion creep = 6.181e-14
E = Activation energies for diffusion creep = 49.0e3
V = Activation volumes for diffusion creep = 0
n = Stress exponents for diffusion creep = 1.8 (from Han and Showman 2011)
m = Grain size exponents for diffusion creep = 1.4
d = Grain size = 1e-2
We were expecting to achieve variable-viscosity within our flows, however the viscosity remains constant at 1e13 (the minimum viscosity limit we set in the parameters):
Going back to equation 14 referenced earlier, our total viscosity being limited to the minimum could be due to the maximum viscosity (μ_max) and effective viscoplastic viscosity (μ_vp / eff) in the second half of equation 14 being much less than the minimum viscosity (μ_min).
Are we missing something about how viscosity is calculated within the visco plastic material model? We are expecting a viscosity range between 1e13 and 1e14, as set by the min/max parameters we set. Does it have to do with the harmonic viscosity averaging scheme we chose? We also ran tests with “harmonic averaging only viscosity”, and our resulting viscosity is still the constant minimum value we set (1e13).
I can also share the full parameter file privately if it would be useful for troubleshooting. Thanks for your time reading this, and please let me know if there are any other details I can include!
Dear Sarah,
It is hard to tell what exactly causes the problem without plugging all the numbers into your equation, but I have some general ideas for how you might go about debugging the problem: In general, it is good to make the problem as simple as possible, so that you can see what is going on.
You could adapt your model for testing so that the temperature and pressure are constant (you can achieve that for the pressure by setting gravity to zero). This way, you can rule out everything to do with spatial averaging. The Viscosity averaging scheme = harmonic is for different compositions, so as long as you do not have compositional fields, that should not make a difference.
Then you could set you minimum viscosity to 1e-100, and you maximum viscosity to 1e100, that way, you would see what viscosity is actually computed. Since large viscosity contrasts cause problems for the linear solver, you could set the linear solver tolerance to 1 so that ASPECT does not actually solve the linear system, since you only want to look at the viscosity anyway.
If this still does not help, you can then make the problem even simpler by just starting to set all your parameters in your flow law to 1 (or whatever the appropriate value is so that they no longer influence the resulting viscosity). At some point, the problem should be so simple that you see why the viscosity is different from what you expect.
Best regards,
Juliane
Hello Juliane,
Thank you for your detailed response, and tips for debugging! I didn’t realize that the viscosity averaging scheme was for different compositions, that’s useful to know. We took your advice in trying to simplify our model: setting gravity to zero, min/max viscosity to 1e-100 and 1e100, linear Stokes solver tolerance to 1. After that, we set the additional parameters to 0 and 1 in the flow law to try to isolate the viscosity.
We were expecting these limited parameters to make our viscosity =1, but we found a resulting 0 viscosity. We did some additional tests to see if there’s an obvious problem I’m missing, but can’t seem to narrow down the cause of the viscosity.
Our current plan is to install another version of ASPECT that we can use as a debugger, modifying the source code to include some print statements so that we can see the values of viscosity that the code is calculating, and hopefully we’ll be able to figure things out from there.
I appreciate your help, and if you have further advice in terms of compiling a debug version of ASPECT, please let me know.
Thanks and best,
Sarah
Hi Sarah,
Just to make sure: You have compiled ASPECT in Debug model (instead of Release mode)? That will give you more information if values are different from what they should be.
And then you did something like
Viscous flow law = diffusion
A = Prefactors for diffusion creep = 1
E = Activation energies for diffusion creep = 0
V = Activation volumes for diffusion creep = 0
n = Stress exponents for diffusion creep = 1
m = Grain size exponents for diffusion creep = 1
d = Grain size = 1
which I guess should give you a viscosity of 0.5, but then the viscosity was 0 instead?
The only reason I can think of would be that the pressure, temperature or strain rate could have a value of infinity or something similarly unreasonable (which should give you an error in debug mode).
Can you post a minimal example that reproduces your problem?
Best regards,
Juliane
Dear Sarah,
Another couple of thoughts:
- what is the strain rate in your models? For reference, I’ve put what your flow law should produce below.
- why not try your diffusion creep flow law with a simpler material model first? For example, you could try the Diffusion dislocation model and set the Dislocation activation energy very high or your Dislocation prefactor very low.
I wouldn’t modify ASPECT until you’ve tried this.
Bob
Hello Bob and Juliane,
Thanks for your responses!
Juliane:
-
Yes, the version of ASPECT we are running is in the default Debug mode. It seems some users use GDB to debug, and while neither my advisor or I have experience with this method, it seems like a better option than modifying the ASPECT source files, as it has many features including the ability to print a stored variable.
-
Yes, we set the visco plastic material model parameters as you described. Here is a simplified version (with as many 0’s and 1’s for parameter values to isolate viscosity) of the parameter file we’re using:
## Global parameters
set Dimension = 2
set Start time = 0
set End time = 10000
set Use years in output instead of seconds = true
set Nonlinear solver scheme = single Advection, iterated defect correction Stokes
set Nonlinear solver tolerance = 1
set Max nonlinear iterations = 2
set CFL number = 0.5
set Maximum time step = 25
set Output directory = 2dviscoplastic
set Pressure normalization = no
#### Parameters describing the model
subsection Formulation
set Formulation = Boussinesq approximation
end
subsection Geometry model
set Model name = box
subsection Box
set X repetitions = 35
set X extent = 35000
set Y repetitions = 25
set Y extent = 10000
end
end
subsection Mesh refinement
set Initial adaptive refinement = 0
set Initial global refinement = 0
set Time steps between mesh refinement = 25
end
subsection Solver parameters
subsection Stokes solver parameters
set Number of cheap Stokes solver steps = 2000
set Stokes solver type = block GMG
end
subsection Newton solver parameters
set Maximum linear Stokes solver tolerance = 1
set Use Eisenstat Walker method for Picard iterations = true
end
end
subsection Boundary velocity model
set Tangential velocity boundary indicators = left, right, top, bottom
end
subsection Compositional fields
set Number of fields = 1
set Names of fields = plastic_strain
end
subsection Initial composition model
set Model name = function
subsection Function
set Variable names = x,y
set Function expression = 0;
end
end
# Composition: fixed on bottom (inflow boundary), free on sides and top
subsection Boundary composition model
set Fixed composition boundary indicators = bottom
set List of model names = initial composition
end
subsection Boundary temperature model
set Fixed temperature boundary indicators = bottom, top
set List of model names = box
subsection Box
set Bottom temperature = 273
set Top temperature = 110
end
end
subsection Initial temperature model
set Model name = function
subsection Function
set Variable names = x,y
set Function constants = pi=3.14159, T1= 273, T2=110, y1=0, y2=10e3, p=10, N=8, Ly=10e3, Lx=35e3
set Function expression = T1 + ((T2-T1)/(y2-y1))*(y-y1) + p*sin(2*pi*y*N / Ly) * sin(2*pi*x*N/Lx)
end end
subsection Heating model
set List of model names = compositional heating
subsection Compositional heating
set Use compositional field for heat production averaging = 1, 0
set Compositional heating values = 1.e-6, 0
end
end
subsection Material model
set Model name = visco plastic
subsection Visco Plastic
set Reference temperature = 273
set Minimum strain rate = 1
set Reference strain rate = 1
set Minimum viscosity = 1e-100
set Maximum viscosity = 1e100
set Define thermal conductivities = true
set Thermal conductivities = 1
set Heat capacities = 1
set Densities = 1
set Thermal expansivities = 1
set Viscosity averaging scheme = harmonic
set Viscous flow law = diffusion
set Prefactors for diffusion creep = 1
set Activation energies for diffusion creep = 0
set Activation volumes for diffusion creep = 0
set Stress exponents for diffusion creep = 1
set Grain size exponents for diffusion creep = 1
set Grain size = 1
set Angles of internal friction = 0
set Cohesions = 0
end
end
subsection Gravity model
set Model name = vertical
subsection Vertical
set Magnitude = 0
end
end
# Post processing
subsection Postprocess
set List of postprocessors = basic statistics, composition statistics, heat flux densities, heat flux statistics, mass flux statistics, matrix statistics, pressure statistics, temperature statistics, topography, velocity statistics, visualization
subsection Visualization
set List of output variables = density, heat flux map, named additional outputs, strain rate, viscosity
set Output format = vtu
set Time between graphical output = 0
set Interpolate output = true
set Output format = vtu
end
end
And Bob:
- the strain rate were were initially using were a minimum strain rate 1e-20 and a reference strain rate 1e-16, which should work for our temperature range of ~100-300K. In attempts to simplify the model, we changed both values to 1 in the example above.
- We didn’t consider trying a simpler material model first, but I like that idea! We’ll take your suggestion and try a Diffusion dislocation model with a very low value for Dislocation prefactor.
I’ll post an update with Bob’s suggestions, and we’ll start trying to use GDB to debug rather than having an alternate version of ASPECT installed and modifying the source code. If anything seems incorrect about the simplified code I posted, please let me know.
I greatly appreciate your time and wisdom!
Best,
Sarah
Dear Sarah,
A couple of issues:
- Your Cohesions are set to zero, so the model is always in plastic failure. This is why your viscosity is zero
- Your gravity is equal to zero and you have free slip boundary conditions. Nothing is driving flow.
May I suggest starting with visco_plastic_diffusion.prm in tests, changing the prm piecewise until you get to your model?
Best wishes,
Bob
@bobmyhill I was the one to suggest to set the gravity to zero, exactly because then nothing drives the flow and you can have a constant temperature, pressure, and strain rate, which is easier for debugging.
But I agree that the cohesion being zero is a problem.
Oops, sorry @jdannberg , I forgot you suggested that. Of course you’re right that it’s a useful thing to do if the viscosity is too low.
@SarahChinski , good luck 
Hi Sarah,
First thing I noticed: The file you sent crashes in debug mode; it’s missing the “set Material averaging = project to Q1 only viscosity” . The part that you can remove is the “set Viscosity averaging scheme = harmonic” not the material averaging.
Second thing: You set the nonlinear solver tolerance to 1, not the linear solver tolerance, so ASPECT still solves the equations. When I run it in debug mode, the model crashes.
Third: I set the Cohesion to 1e100 (along the lines of what Bob suggested) and that gave me the expected viscosity of 0.5. So that is the result we expect. From there, you can step by step go back to the actual viscosity you want changing the parameters one by one, and along the way you’ll hopefully figure out what went wrong.
Best,
Juliane
