Skip to content
GitLab
Commit 6be7933a authored by Lena Noack's avatar Lena Noack
Browse files

Add plate tectonics benchmark

parent 907346b2
Hide whitespace changes
Inline Side-by-side
Octave/read_snap.m
View file @ 6be7933a
......@@ -36,7 +36,7 @@ function [data,input,program_dir,fname,fpath] = read_snap(finput=' ')
data.solidus_diff_w=1;, data.liquidus=1;, data.ds_dz=1;, data.dl_dz=1;, data.DeltaF=1;
data.serp=1;, data.total_mass_H2O=1;, data.oc_vol_tot=1;, data.fc_vol_tot=1;, data.cc_vol_tot=1;
data.cr_age=1;, data.output_tracer=1;, data.ref_g=1;, data.shear=1;, data.bulk_mod=1;
data.ref_p=1;, data.T_old=1;, data.entropy=1;, data.FP=1;, data.FG=0;, data.CO2=1;
data.ref_p=1;, data.T_old=1;, data.entropy=1;, data.FP=1;, data.FG=0;, data.MP=0;, data.CO2=1;
data.U238=0;, data.U235=0;, data.Th232=0;, data.K40=0;
data.m_vel_u=0.0;,data.m_vel_v=0.0;,data.m_vel_w=0.0;,data.mprod=0.0;,data.poros=0.0;
......@@ -109,6 +109,7 @@ function [data,input,program_dir,fname,fpath] = read_snap(finput=' ')
elseif(name == "en"), hnl=nl+2;,hny=ny+2*n3D;,hnr=nr+2;,data.entropy=read_matrix(hnl,hny,hnr,data_snap);
elseif(name == "fp"), hnl=nl+2;,hny=ny+2*n3D;,hnr=nr+2;,data.FP=read_matrix(hnl,hny,hnr,data_snap);
elseif(name == "fg"), hnl=nl+2;,hny=ny+2*n3D;,hnr=nr+2;,data.FG=read_matrix(hnl,hny,hnr,data_snap);
elseif(name == "ph"), hnl=nl+2;,hny=ny+2*n3D;,hnr=nr+2;,data.MP=read_matrix(hnl,hny,hnr,data_snap);
elseif(name == "co"), hnl=nl+2;,hny=ny+2*n3D;,hnr=nr+2;,data.CO2=read_matrix(hnl,hny,hnr,data_snap);
elseif(name == "ul"), hnl=nl+2;,hny=ny+2*n3D;,hnr=nr+2;,data.U238=read_matrix(hnl,hny,hnr,data_snap);
elseif(name == "us"), hnl=nl+2;,hny=ny+2*n3D;,hnr=nr+2;,data.U235=read_matrix(hnl,hny,hnr,data_snap);
......
bin/input_PTBench.txt 0 → 100644
View file @ 6be7933a
! Geometry 0
geometry = 0 ! 0-box, 1-sph. annulus, 2-cylinder
dim = 2
nl = 375 !3000 ! number of grid points in lateral direction (2 ghost points are added when reflective boundary is used)
nr = 84 !670 ! number of grid points in radial direction (+2 ghost points)
lmin = 0.0 ! only used when geometry=0
lmax = 3000000 ! only used when geometry=0
rmin = 670000 !
rmax = 1340000 ! rmax=rmin+1.0
phi_factor = 0 ! only used when geometry>0, 1.0: full annulus, 0.5: halph annulus, ...
! Adaptivity 0
ad_factor = 0 !2 ! 0-uniform box grid, 1-refine area (ad_min/ad_max values needed), 2-smallest/largest resolution (res_min/res_max) is set at (ad_l=0.5*lmax, ad_r=rmax)
res_ratio_l = 0.2
res_ratio_r = 0.2
! Bnd_cond 0
Ttop = 0 ! temperature at surface, until now: has to be 0 or dimensional value (e.g. 288 for Earth)
Tbottom = 1440 ! temperature at CMB, until now: has to be 1 or dimensional value (e.g. 3900 for Earth)
noslip_t = 0 ! 0: free-slip, 1: no-slip at surface boundary
noslip_b = 0 ! 0: free-slip, 1: no-slip at CMB boundary
inner_bound = 0 ! 1: isolated
outer_bound = 0 ! 1: isolated
open_b_bound = 0 ! 1: open boundary at bottom
periodic = 0 ! 0: reflective boundary, 1: periodic boundary, is currently tested, normally: 0 for box and 1 for sphere
bnd_l = 2
bnd_r = 1
! Restart_TS 0
use_snap = 0 ! insert number of TS from the snapshot, from which the simulation shall re-start
last_snap = 0 ! set to 1, then the last snapshot is automatically used (use_snap is ignored)
Cdt = 1 ! Courant factor when positive value (e.g. 0.1), Delta criterium factor when negative (e.g. -10.0)
nbiter = 0 ! max number of time steps, not used if tmax>0
iter_out = 0 ! output every ... time steps
tmax = 1.5e+7 ! 15Myr ! maximal time, stop the simulation when reached; if tmax=0 then either nbiter or conv is used to stop the simulation
t_output = 5.0e+5 ! 0.5Myr ! output time
! Initial 0
ampl = 0.0 !0.1 ! amplitude of initially applied sphercial harmonics
sph_l = -1 ! mode of spherical harmonics, chose +/-1.0 times the aspect ratio for box geometry and number of preferred plumes for polar geometry
linear = 0 ! initial temperature profile: 0: T=Tini everywhere; 1: linear, 2: conductive profile, 3: convective profile with TBLs
Tini = 1300 ! initial mantle temperature for profile 3
s_bot = 67000 ! lower thermal boundary layer (TBL) for profile 3
s_top = 67000 ! upper thermal boundary layer (TBL) for profile 3
TinclTref = 1
! Viscosity 0
n_fac = 0 ! Arrhenius AND FK: non-Newtonian factor (1.0 - Newtonian, >1.0 non-Newtonian)
p_fac = 0 ! Newtonian: grain-size factor
iniStrRate = 1.0e-15 ! 10^-15 1/s ! Initial Strain rate
StartNN = 0 ! Number of timesteps after which non-Newtonian law is applied
Tref = 0.0 ! reference temperature for Rayleigh number
zref = 0.0 ! reference depth for Rayleigh number
gam = 0.0 ! FK Viscosity: gamma_T
gamP = 0.0 ! FK Viscosity: gamma_p
E = 0 ! Arrhenius: actvation energy
E_wet = 0
V = 0 ! Arrhenius: activation volume
V_wet = 0
E_dis = 0 ! Arrhenius, non-Newtonian: activation energy
E_dis_wet = 0
V_dis = 0 ! Arrhenius, non-Newtonian: activation volume
V_dis_wet = 0
A_dif = 0 ! Arrhenius, mixed Newtonian: prefactor Newtonian: A' = A * eta_ref / d_ref^p
A_dif_wet = 0
A_dis = 0 ! Arrhenius, mixed Newtonian: prefactor non-Newtonian: A' = A * eta_ref^n * (kappa/D^2)^(n-1)
A_dis_wet = 0
water_r_dif = 0 ! is used in Arrhenius when wet formulation used
water_r_dis = 0
T0 = 273.15 ! Arrhenius: surface temperature (Tsurf/DeltaT)
AdVisc = 0 ! adapt viscosity at boundaries, such that viscosity at ri/ro is correct
ViscMin = 1.0e+20
ViscMax = 1.0e+25
! Plasticity 0
plasticity = 0 ! 0-Bingham, 1-nonsmooth, 2-angle-dependent
friction = 0 !1.0 ! non-dim factor: rho*g*D^3 / kappa / eta
mantle_fr_angl = 0 ! plasticity = 2: friction angle in degrees
mantle_fr_angl_m = 0 ! plasticity = 2: minimal friction angle for strain weakening
mantle_fr_coh = 0 ! plasticity = 2: surface cohesion, 20 MPa
YieldStr0 = 0 ! plasticity = 0/1: surface yield stress (YS_0); plasticity = 2: surface cohesion
YieldStrZ = 0 ! plasticity = 0/1: increase of yield stress with depth (YS_z*z); plasticity = 2: non-dim. factor for pressure-dependent part
strain_min = 0 ! plasticity = 2: for linear strain-weakening of friction angle
strain_max = 0 ! plasticity = 2: for linear strain-weakening of friction angle
UseRefDensity = 2
moveRef = 1
weak_x0 = 1487796 ! = 1500km - 0.5*weak_w/sin(weak_angle)
weak_x1 = 1530000 ! = 1530km
weak_angle = -35 ! angle of weak zone in degrees
weak_d = 82000
weak_w = 14000
weak_visc = 1.0e+20 ! ref visc = 1e20
weak_fr_angl = 0
weak_fr_angl_m = 0
weak_fr_coh = 0
weak_rho = 3250 ! ref density = 3200
weak_T = 0
weak_cp = 1250 ! ref heat capacity = 1250
weak_k = 2.5 ! ref conductivity (mantle) = 4
weak_alpha = 2.5e-5
weak_h = 0
weak_w_H2O = 0
weak_w_r = 0 ! set either 0 or 1 for dry or wet...
vel_r = 5 ! 5cm/yr vel imposed on right side
vel_z_in = 110000 ! 110km, depth until which inflow occurs
vel_z_out = 130000 ! 130km, depth from which on outflow occurs, linear velocity inbetween
OP_BOC = 7000 ! 7km
OP_BOC_rho = 3250
OP_BOC_T = 0
OP_BOC_eta = 1.0e+23
OP_BOC_fr_angl = 0
OP_BOC_fr_angl_m = 0
OP_BOC_fr_coh = 0
OP_BOC_cp = 1250
OP_BOC_k = 2.5
OP_BOC_alpha = 2.5e-5
OP_BOC_h = 0
OP_BOC_w_H2O = 0
OP_BOC_w_r = 0 ! set either 0 or 1 for dry or wet...
OP_SHB = 32000 ! 32km
OP_SHB_rho = 3250
OP_SHB_T = 0
OP_SHB_eta = 1.0e+23
OP_SHB_fr_angl = 0
OP_SHB_fr_angl_m = 0
OP_SHB_fr_coh = 0
OP_SHB_cp = 1250
OP_SHB_k = 2.5
OP_SHB_alpha = 2.5e-5
OP_SHB_h = 0
OP_SHB_w_H2O = 0
OP_SHB_w_r = 0 ! set either 0 or 1 for dry or wet...
OP_ThL = 43000 ! 43km
OP_ThL_rho = 3250
OP_ThL_T = 0
OP_ThL_eta = 1.0e+23
OP_ThL_fr_angl = 0
OP_ThL_fr_angl_m = 0
OP_ThL_fr_coh = 0
OP_ThL_age = 4.0e+7 ! 40 Myr[s]*kappa/D^2
OP_ThL_cp = 1250
OP_ThL_k = 2.5
OP_ThL_alpha = 2.5e-5
OP_ThL_h = 0
OP_ThL_w_H2O = 0
OP_ThL_w_r = 0 ! set either 0 or 1 for dry or wet...
SP_BOC = 8000 ! 8km
SP_BOC_rho = 3250
SP_BOC_T = 0
SP_BOC_eta = 1.0e+20
SP_BOC_fr_angl = 0
SP_BOC_fr_angl_m = 0
SP_BOC_fr_coh = 0
SP_BOC_cp = 1250
SP_BOC_k = 2.5
SP_BOC_alpha = 2.5e-5
SP_BOC_h = 0
SP_BOC_w_H2O = 0
SP_BOC_w_r = 0 ! set either 0 or 1 for dry or wet...
SP_SHB = 32000 ! 32km
SP_SHB_rho = 3250
SP_SHB_T = 0
SP_SHB_eta = 1.0e+23
SP_SHB_fr_angl = 0
SP_SHB_fr_angl_m = 0
SP_SHB_fr_coh = 0
SP_SHB_cp = 1250
SP_SHB_k = 2.5
SP_SHB_alpha = 2.5e-5
SP_SHB_h = 0
SP_SHB_w_H2O = 0
SP_SHB_w_r = 0 ! set either 0 or 1 for dry or wet...
SP_ThL = 70000 ! 70km
SP_ThL_rho = 3250
SP_ThL_T = 0
SP_ThL_eta = 1.0e+23
SP_ThL_fr_angl = 0
SP_ThL_fr_angl_m = 0
SP_ThL_fr_coh = 0
SP_ThL_age = 7.0e+7 ! 70 Myr[s]*kappa/D^2
SP_ThL_cp = 1250
SP_ThL_k = 2.5
SP_ThL_alpha = 2.5e-5
SP_ThL_h = 0
SP_ThL_w_H2O = 0
SP_ThL_w_r = 0 ! set either 0 or 1 for dry or wet...
T_init_plate = 7 ! set an initial temperature field
adiabat = 0.25 ! in Degree/km
isotherm = 800
DiT = 0 !1
fieldsAll = 1 !0
fieldP = 1
! Ra_H0 0
Ra = 0.0 ! Rayleigh number at Tref/zref
eta_ref = 1.0e+20 ! set Ra automatically
H0 = 0.0 ! initial heat source density, set to H0>1 for internal heating!
lambda = 0.0 ! decay of internal heat sources with time: H=H0*exp(-t*lambda)
B = 1 ! Buoyancy ratio for chemical diffusion, needed in momentum equation
Le = 0 ! Lewis number for chemical diffusion, if 0, then no double-diffusive solver is used
LeV = 0
comp_init = 0 !-2 ! +-1 = linear, +-2 = rho value; negative: add cylinder, positive: add random perturbation
comp_ampl = 0 ! Amplitude of composition perturbation
C_ref = 0 !1 ! reference composition (should correspond to reference density)
! Particles 0
use_part = 1 ! use particles instead of compositional field approach
part_n = 0 ! total number of particles, ignored if part_n_cell>0
part_n_cell = 10 ! number of particles per cell
part_depth = 1.0
restart_p = 1
part_rk = 4
! Numerics 0
debug = 0 ! set between 1 and 5 to obtain detailed information
max_outer = 50 ! number of max outer iterations (coupled energy-momentum solver), used if e_solver>0
max_inner = 3 ! number of max inner iterations (coupled pressure-velocity solver), used if m_solver>0
noBouss = 1 !0 ! if set to one then no Boussinesq is used in momentum solver
e_solver = -1 !3 ! -1: no e-solver, 0: direct first-order explicit, 1: second-order iterative explicit, 2: first-oder iterative implicit, 3: second-order iterative implicit
m_solver = 1 !
c_solver = 0
convT = 1e-8 ! criterion for convergence of outer iteration, used if e_solver>0 and number of outer iterations < max_outer
convV = 1e-4 ! criterion for convergence of inner iteration, used if m_solver>0 and number of inner iterations < max_inner
conv = 0 ! stops the simulation if convergence of time steps is reached, not used if tmax>0 or nbiter is reached
test = 0 ! output of pressure, when used IDL keyword /test has to be used, as well
beta = 1.0 ! for flux limiter; 1: minmod, 2: superbee
average = 1 ! averaging scheme to be used for viscosity: 0=old, 1=harmonic, 2=arithmetic, 3=geometric
p_average = 2 ! averaging scheme to be used for particles: 1=harmonic, 2=arithmetic, 3=geometric
p_average2 = 2 ! for density averaging; if 0 then set to p_average
p_penalty = 1.0e-7 ! this value shoud be 1.0e-7 and is divided by the maximal viscosity in the stokes equation
nmesures = 0 ! currently not used... should be included again for the cases when quasi-steady-state is reached to gain average values for the last X timesteps
! Physics 0
g = 9.81 ! still has to be tested with new version if dimensional quantities really can be used...
alpha = 2.5e-5
rho = 3200
kappa = 1.0e-6
grain_size = 1
mantle_w_H2O = 0 ! reference water value
mantle_w_r = 0 ! set either 0 or 1 for dry or wet...
k = 183.33 ! 183.33, ref value = 4
Tm_ref = 0
Cp = 1250
EOF = 1
src/parameters.f90
View file @ 6be7933a
......@@ -50,7 +50,7 @@ module parameters
real(dp)::lat_heat=0.0_dp
real(dp)::melt_depth=1.0_dp
real(dp)::melt_z_change=1.0_dp
real(dp)::dehydr_melt=0.0_dp
real(dp)::dehydr_melt=1.0_dp ! if depletion reaches this value, then material is completely dehydrated (new version: dehydrate using dH2O)
real(dp)::serp_instable_m=0.0_dp,serp_instable_T=0.0_dp
real(dp)::serp_instable_P1=0.0_dp,serp_instable_P2=0.0_dp
real(dp)::serp_instable_w=0.0_dp,serp_instable_w_T=0.0_dp,serp_instable_w_max=0.0_dp
......
src/plasticity.f90
View file @ 6be7933a
......@@ -1101,6 +1101,9 @@ module plasticity
y=mesh%rc(k)
x_w = 0.5_dp*mesh%lmax
if (geom.gt.0) then
x_w = 0.0_dp
endif
if (x.le.x_w) then ! left half of box
z_l = pW%OP_BOC+pW%OP_SHB+pW%OP_ThL
z_r = pW%SP_BOC+pW%SP_SHB + pW%SP_ThL
......@@ -1138,6 +1141,9 @@ module plasticity
y=mesh%rc(k)
x_w = 0.5_dp*mesh%lmax
if (geom.gt.0) then
x_w = 0.0_dp
endif
if (x.le.x_w) then ! left half of box
z_l = pW%OP_BOC+pW%OP_SHB+pW%OP_ThL
z_r = pW%SP_BOC+pW%SP_SHB + pW%SP_ThL
......@@ -1174,6 +1180,9 @@ module plasticity
y=mesh%rc(k)
x_w = 0.5_dp*mesh%lmax
if (geom.gt.0) then
x_w = 0.0_dp
endif
if (x.le.x_w) then ! left half of box
z_l = pW%OP_BOC+pW%OP_SHB+pW%OP_ThL
z_r = pW%SP_BOC+pW%SP_SHB + pW%SP_ThL
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment