Hi!
I would like to use the viscoplastic model with a phase transition of course at the 410 and 660. However, when I have changed from compositing to visco plastic my compositional fields (LM->UM at 660) are not switching (see figure). I attach my .prm, do I have wrong the syntax please?
Thanks a lot!
Elodie
–
#!/usr/bin/env python3
Global parameters
set Dimension = 2
set Use years in output instead of seconds = true
set Start time = 0
set End time = 5000.e6 # 4.5e9 #400.e6
set Output directory = output
set Nonlinear solver scheme = iterated Stokes #single Advection, iterated Stokes
set Nonlinear solver tolerance = 5e-4 # -4 -> -5
set Max nonlinear iterations = 150
set Max nonlinear iterations in pre-refinement = 20
set Pressure normalization = no # free surface
set Maximum time step = 1e6
#set Number of cheap Stokes solver steps = 200 # aspect 2.2
set Adiabatic surface temperature = 1600 # need 4 initial adiabatic temp.
set CFL number = 0.5
set Resume computation = auto
subsection Checkpointing
set Steps between checkpoint = 100
set Time between checkpoint = 0
end
subsection Termination criteria
set Checkpoint on termination = false
set End step = 1
set Termination criteria = end time #end step
end
subsection Solver parameters
subsection Stokes solver parameters
#set GMRES solver restart length = 50
set Linear solver tolerance = 1e-5
set Maximum number of expensive Stokes solver steps = 1000
set Number of cheap Stokes solver steps = 200 # not using too much, only get a guess
end
end
Model Geometry
subsection Geometry model
set Model name = spherical shell
subsection Spherical shell
set Inner radius = 3485e3
set Outer radius = 6371e3 #5961e3
set Opening angle = 360
end
end
Mesh refinement
subsection Mesh refinement
set Coarsening fraction = 0.2
set Refinement fraction = 0.8
set Initial adaptive refinement = 2
set Initial global refinement = 5
set Time steps between mesh refinement = 0
set Strategy = minimum refinement function #velocity, temperature, minimum refinement function #temperature #, composition, temperature, minimum refinement function
subsection Minimum refinement function
set Variable names = depth,t
set Coordinate system = depth
set Function expression = if(depth<=760e3,7,5)
#if(depth<560e3,0,if(depth>760e3,0,6))
if(depth<560e3,5,if(depth>760e3,5,6))
end
end
subsection Boundary composition model
#set Allow fixed composition on outflow boundaries = false # false for models without melt
set Fixed composition boundary indicators = bottom, top
end
subsection Boundary temperature model
#set Allow fixed temperature on outflow boundaries = true # default
set Fixed temperature boundary indicators = bottom, top
end
subsection Boundary velocity model
set Tangential velocity boundary indicators = left, right # 360
set Zero velocity boundary indicators = inner
end
subsection Mesh deformation
set Mesh deformation boundary indicators = top: free surface
subsection Free surface
set Free surface stabilization theta = 0.5
end
end
Number and names of compositional fields
subsection Compositional fields
set Number of fields = 4
set Names of fields = mantle_l, mantle_u, viscous_strain, plastic_strain #add strain field
end
Spatial domain of different compositional fields
Horizontal layers
subsection Initial composition model
set Model name = function
subsection Function
set Coordinate system = depth
set Variable names = depth,t
set Function expression =
if(depth>=660e3, 1, 0);
if(depth<660e3, 1, 0);
if(depth>=660e3, 0, 0);
if(depth<660e3, 0, 0)
end
end
Composition fixed on bottom, free on sides and top
subsection Boundary composition model
set Model name = initial composition
end
Set temperature Model
subsection Boundary temperature model
set Fixed temperature boundary indicators = top, bottom
set List of model names = spherical constant
subsection Spherical constant
set Inner temperature = 4000 #3300
set Outer temperature = 273 #273 for surface ############
end
end
Temperatur model
subsection Initial temperature model
#set Model name = adiabatic
set List of model names = adiabatic, function
set List of model operators = add
set Model name = adiabatic
subsection Adiabatic
set Age bottom boundary layer = 100.e6
set Age top boundary layer = 100.e6
set Radius = 0
set Amplitude = 0
set Position = center
set Subadiabaticity = 0
subsection Function
set Variable names = x,y
set Function constants = R=6371.e3
set Function expression =
if(sqrt(xx+yy)<=R-660.e3, 1, 0);
if(sqrt(xx+yy)>R-660.e3, 1, 0);
if(sqrt(xx+yy)<=R-660.e3, 0, 0);
if(sqrt(xx+yy)>R-660.e3, 0, 0)
end
end
subsection Function # temp. anomalie
set Coordinate system = cartesian
set Variable names = x,y,z
set Function constants = R=6371.e3, shift=3.e6, pi=3.14159, const=200
set Function expression =
if(sqrt(pow(x-shift+1e6,2)+pow(y-shift-1e6,2))<200.e3,const,
if(sqrt(pow(x-shift-2e6,2)+pow(y-shift+2e6,2))<250.e3,const+100,
if(sqrt(pow(x-shift-1.5e6,2)+pow(y+shift-2e6,2))<200.e3,const+200,
if(sqrt(pow(x-shift+1e6,2)+pow(y+shift+1e6,2))<400.e3,const,
if(sqrt(pow(x+shift-1e6,2)+pow(y+shift+1e6,2))<200.e3,const+100,
if(sqrt(pow(x+shift+2e6,2)+pow(y+shift-1.5e6,2))<300.e3,const,
if(sqrt(pow(x+shift+1e6,2)+pow(y-shift+1e6,2))<400.e3,const,
if(sqrt(pow(x+shift-1e6,2)+pow(y-shift-1e6,2))<200.e3,const+100,
if(sqrt(pow(x-shift-0.9e6,2)+pow(y-shift-1.1e6,2))<150.e3,-50.,
if(sqrt(pow(x,2)+pow(y,2))>R-200.e3 &&
sqrt(pow(x,2)+pow(y,2))<R-100.e3 &&
atan(y/x)> 43pi/180 &&
atan(y/x)< 46pi/180
,+0,
0)))))))) )) # mesh artefacts #if(sqrt(pow(x-shift-1.5e6,2)+pow(y-shift-1.5e6,2))<100.e3,+300.,
#set Function expression = if(r>0.e6 && phi>0 && phi<3000 && theta>0 && theta<3000,3100,2900)
#if(sqrt(pow(x,2)+pow(y,2))<1e5,3000,2900)
end
end
subsection Adiabatic conditions model # need to initialize for initial adiabatic temp.
set Model name = compute profile
#set Model name = initial profile
subsection Compute profile #Initial profile
set Variable names = x,y
set Function constants = R=6371.e3
set Function expression =
if(sqrt(xx+yy)<=R-660.e3, 1, 0);
if(sqrt(xx+yy)>R-660.e3, 1, 0);
if(sqrt(xx+yy)>=R-660.e3, 0, 0);
if(sqrt(xx+yy)<R-660.e3, 0, 0)
end
end
subsection Heating model
set List of model names = adiabatic heating, compositional heating #, shear heating
#set List of model names = compositional heating
subsection Compositional heating
set Use compositional field for heat production averaging = 1,1,1,0,0
set Compositional heating values = 0.02e-6
end
end
############################ MATERIAL MODEL ######################################################################
subsection Material model
set Model name = visco plastic #visco plastic #melt simpe
set Material averaging = harmonic average
############################ MATERIAL MODEL ######################################################################
#subsection Compositing
set Viscosity = visco plastic
set Density = visco plastic
set Thermal expansion coefficient = visco plastic
set Specific heat = visco plastic
set Thermal conductivity = visco plastic
set Compressibility = visco plastic
set Entropy derivative pressure = visco plastic
set Entropy derivative temperature = visco plastic
set Reaction terms = composition reaction acg
#end
subsection Visco Plastic
set Define transition by depth instead of pressure = true
set Phase transition Clapeyron slopes = background:1.5e6|-2e6, mantle_l:1.5e6|-2e6, mantle_u:1.5e6|-2e6 #based on Faccenda and Dal Zilio 2017
set Phase transition depths = background: 410e3|660e3, mantle_l: 410e3|660e3, mantle_u:410e3|660e3
set Phase transition widths = background: 13e3|5e3, mantle_l: 13e3|5e3, mantle_u:13e3|5e3 #13km and 5km based on Arredondo and Billen 2016, 403.5e3|415.5e3, 657.5e3|662.5e3
set Phase transition temperatures = background: 1810|1940, mantle_l: 1810|1940, mantle_u: 1810|1940 #based on Faccenda and Dal Zilio 2017
#set Maximum thermal expansivity =
# Reference temperature and viscosity
set Reference temperature = 293
set Reference viscosity = 1e21
# The minimum strain-rate helps limit large viscosities values that arise
# as the strain-rate approaches zero.
# The reference strain-rate is used on the first non-linear iteration
# of the first time step when the velocity has not been determined yet.
set Minimum strain rate = 1.e-20
set Reference strain rate = 1.e-15
# Limit the viscosity with minimum and maximum values
set Minimum viscosity = 1e19
set Maximum viscosity = 1e24
# Harmonic viscosity averaging
set Viscosity averaging scheme = harmonic
# Choose to have the viscosity (pre-yield) follow a dislocation
# diffusion or composite flow law. Here, dislocation is selected
# so no need to specify diffusion creep parameters below, which are
# only used if "diffusion" or "composite" option is selected.
set Viscous flow law = diffusion
from cizkov
order: background , mantle_low , mantle_up
set Thermal diffusivities = 1.e-6
set Thermal expansivities = 3.e-5
set Heat capacities = 1250 # specific heat
set Densities = 3416., 3416. , 3416.,3416,3416 #3750
set Thermal expansivities = 3.0e-5 , 3.0e-5, 2.7e-5
Dislocation creep
set Prefactors for dislocation creep = background:1e-30|1e-22|1e-22, mantle_l:1e-30|1e-22|1e-22, mantle_u:1e-30|1e-22|1e-22, viscous_strain:1e-30, plastic_strain:1e-22
set Stress exponents for dislocation creep = 3.5
set Activation energies for dislocation creep = background:4.8e5|4.8e5|4.8e5, mantle_l:4.8e5|4.8e5|4.8e5, mantle_u:4.8e5|4.8e5|4.8e5, viscous_strain:4.8e5, plastic_strain:4.8e5
set Activation volumes for dislocation creep = background:11.e-6|11.e-6|11.e-6, mantle_l:11.e-6|11.e-6|11.e-6, mantle_u:11.e-6|11.e-6|11.e-6, viscous_strain:11.e-6,plastic_strain:11.e-6
Diffusion creep
set Prefactors for diffusion creep = background:2.5e-17|0.125e-9|0.125e-9,mantle_l:2.5e-17|0.125e-9|0.125e-9, mantle_u:2.5e-17|0.125e-9|0.125e-9,viscous_strain:2.5e-17,plastic_strain:0.125e-9
set Stress exponents for diffusion creep = 1.0
set Activation energies for diffusion creep = background:1.6e5|3.35e5|3.35e5, mantle_l:1.6e5|3.35e5|3.35e5, mantle_u:1.6e5|3.35e5|3.35e5, viscous_strain:1.6e5, plastic_strain:3.35e5
set Activation volumes for diffusion creep = background:2.5e-6|4.8e-6|4.8e-6, mantle_l:2.5e-6|4.8e-6|4.8e-6, mantle_u:2.5e-6|4.8e-6|4.8e-6, viscous_strain:2.5e-6, plastic_strain:4.8e-6
set Grain size = 1 #1e-3
set Grain size exponents for diffusion creep = 0 #default
# Plasticity parameters
set Angles of internal friction = 0,0,0,0,0 # 30 # Angle Mohr
set Cohesions = 1.e6,1000.e6,40.e6,1000e6,40e6 #40.e6 # 20.e6 # Displ. Mohr
add strain weakening
set Strain weakening mechanism = plastic weakening with plastic strain and viscous weakening with viscous strain
set Start plasticity strain weakening intervals = 0.5
set Start prefactor strain weakening intervals = 0.5
set Cohesion strain weakening factors = 0.5
set Friction strain weakening factors = 0.5
set Prefactor strain weakening factors = 0.5
set End prefactor strain weakening intervals = 1.5
set End plasticity strain weakening intervals = 1.5
end
end
#######################
Gravity model
subsection Gravity model
set Model name = radial constant
subsection Radial constant
set Magnitude = 10
end
end
#######################
Postprocessing
subsection Postprocess
set List of postprocessors = velocity statistics, basic statistics, temperature statistics,
visualization, topography, composition statistics, depth average, dynamic topography, heat flux statistics
subsection Visualization
set List of output variables = density, viscosity, strain rate, stress, nonadiabatic temperature, depth, shear stress, vertical heat flux, heat flux map, named additional outputs
set Output format = vtu
set Time between graphical output = 1e7 # each time step
set Interpolate output = true
set Number of grouped files = 0
subsection Heat flux map
set Output point wise heat flux = false # default=false
end
end
subsection Depth average
set Number of zones = 500
set Output format = vtu #, gnuplot, vtk #gnuplot
end
end