Skip to content

Commit

Permalink
more non dim
Browse files Browse the repository at this point in the history
  • Loading branch information
@albert-de-montserrat
albert-de-montserrat committed Mar 26, 2024
1 parent abf0b39 commit 949e27b
Show file tree
Hide file tree
Showing 3 changed files with 88 additions and 42 deletions.
27 changes: 15 additions & 12 deletions miniapps/convection/Particles2D_nonDim/Layered_convection2D.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,17 +52,18 @@ end
T[i + 1, j] = nondimensionalize(273e0K, CharDim)

elseif nondimensionalize(0e0km, CharDim) ≤ (depth) < nondimensionalize(35e3km, CharDim)
dTdZ = (923-273)/35e3
dTdZ = nondimensionalize((923-273)/35e3 * K/km, CharDim)

offset = nondimensionalize(273e0K, CharDim)
T[i + 1, j] = (depth) * dTdZ + offset

elseif nondimensionalize(110e3km, CharDim) > (depth) ≥ nondimensionalize(35e3km, CharDim)
dTdZ = (1492-923)/75e3
dTdZ = nondimensionalize((1492-923)/75e3 * K/km, CharDim)
offset = nondimensionalize(923K, CharDim)
T[i + 1, j] = (depth - nondimensionalize(35e3km, CharDim)) * dTdZ + offset

elseif (depth) ≥ nondimensionalize(110e3km, CharDim)
dTdZ = (1837 - 1492)/590e3
dTdZ = nondimensionalize((1837 - 1492)/590e3 * K/km, CharDim)
offset = nondimensionalize(1492e0K, CharDim)
T[i + 1, j] = (depth - nondimensionalize(110e3km, CharDim)) * dTdZ + offset

Expand All @@ -72,22 +73,24 @@ end
end

# Thermal rectangular perturbation
function rectangular_perturbation!(T, xc, yc, r, xvi, thick_air)
function rectangular_perturbation!(T, xc, yc, r, xvi, thick_air, CharDim)

@parallel_indices (i, j) function _rectangular_perturbation!(T, xc, yc, r, x, y)
@parallel_indices (i, j) function _rectangular_perturbation!(T, xc, yc, r, CharDim, x, y)
@inbounds if ((x[i]-xc)^2 ≤ r^2) && ((y[j] - yc - thick_air)^2 ≤ r^2)
depth = -y[j] - thick_air
dTdZ = (2047 - 2017) / 50e3
offset = 2017
T[i + 1, j] = (depth - 585e3) * dTdZ + offset
dTdZ = nondimensionalize((2047 - 2017)K / 50e3km, CharDim)
offset = nondimensionalize(2017e0K, CharDim)
T[i + 1, j] = (depth - nondimensionalize(585e3km, CharDim)) * dTdZ + offset
end
return nothing
end

ni = length.(xvi)
@parallel (@idx ni) _rectangular_perturbation!(T, xc, yc, r, xvi...)
@parallel (@idx ni) _rectangular_perturbation!(T, xc, yc, r, CharDim, xvi...)

return nothing
end

## END OF HELPER FUNCTION ------------------------------------------------------------

## BEGIN OF MAIN SCRIPT --------------------------------------------------------------
Expand Down Expand Up @@ -139,7 +142,7 @@ function main2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", do_vtk =false)
# STOKES ---------------------------------------------
# Allocate arrays needed for every Stokes problem
stokes = StokesArrays(ni, ViscoElastic)
pt_stokes = PTStokesCoeffs(li, di; ϵ=1e-4, CFL = 0.75 / √2.1)
pt_stokes = PTStokesCoeffs(li, di; ϵ=1e-7, CFL = 0.9 / √2.1)
# ----------------------------------------------------

# TEMPERATURE PROFILE --------------------------------
Expand All @@ -153,7 +156,7 @@ function main2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", do_vtk =false)
thermal_bcs!(thermal.T, thermal_bc)
Tbot = thermal.T[1, 1]
Ttop = thermal.T[1, end]
rectangular_perturbation!(thermal.T, xc_anomaly, yc_anomaly, r_anomaly, xvi, thick_air)
rectangular_perturbation!(thermal.T, xc_anomaly, yc_anomaly, r_anomaly, xvi, thick_air, CharDim)
@parallel (JustRelax.@idx size(thermal.Tc)...) temperature2center!(thermal.Tc, thermal.T)
# ----------------------------------------------------

Expand Down Expand Up @@ -392,4 +395,4 @@ end

# run main script

main2D(igg; figdir = figdir, ar = ar, nx = nx, ny = ny, do_vtk = do_vtk);
# main2D(igg; figdir = figdir, ar = ar, nx = nx, ny = ny, do_vtk = do_vtk);
95 changes: 69 additions & 26 deletions miniapps/convection/Particles2D_nonDim/Layered_rheology.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,26 +3,70 @@
function init_rheologies(CharDim; is_plastic = true)

# Dislocation and Diffusion creep
disl_upper_crust = DislocationCreep(A=5.07e-18, n=2.3, E=154e3, V=6e-6, r=0.0, R=8.3145)
disl_lower_crust = DislocationCreep(A=2.08e-23, n=3.2, E=238e3, V=6e-6, r=0.0, R=8.3145)
disl_lithospheric_mantle = DislocationCreep(A=2.51e-17, n=3.5, E=530e3, V=6e-6, r=0.0, R=8.3145)
disl_sublithospheric_mantle = DislocationCreep(A=2.51e-17, n=3.5, E=530e3, V=6e-6, r=0.0, R=8.3145)
diff_lithospheric_mantle = DislocationCreep(A=2.51e-17, n=1.0, E=530e3, V=6e-6, r=0.0, R=8.3145)
diff_sublithospheric_mantle = DislocationCreep(A=2.51e-17, n=1.0, E=530e3, V=6e-6, r=0.0, R=8.3145)
disl_upper_crust = DislocationCreep(
A=5.07e-18Pa^(-23//10) / s, # units are Pa^(-n) / s
n=2.3,
E=154e3J/mol,
V=6e-6m^3/mol,
r=0.0,
R=8.3145J/mol/K
)
disl_lower_crust = DislocationCreep(
A=2.08e-23Pa^(-32//10) / s, # units are Pa^(-n) / s
n=3.2,
E=238e3J/mol,
V=6e-6m^3/mol,
r=0.0,
R=8.3145J/mol/K
)
disl_lithospheric_mantle = DislocationCreep(
A=2.51e-17Pa^(-35//10) / s, # units are Pa^(-n) / s
n=3.5,
E=530e3J/mol,
V=6e-6m^3/mol,
r=0.0,
R=8.3145J/mol/K
)
disl_sublithospheric_mantle = DislocationCreep(
A=2.51e-17Pa^(-35//10) / s, # units are Pa^(-n) / s
n=3.5,
E=530e3J/mol,
V=6e-6m^3/mol,
r=0.0,
R=8.3145J/mol/K
)
diff_lithospheric_mantle = DiffusionCreep(
A=2.51e-17Pa^(-1 - 0) / s * m^(-0), # units are Pa^(-n - r) / s * m^(-p)
n=1.0,
E=530e3J/mol,
V=6e-6m^3/mol,
r=0.0,
p=0.0,
R=8.3145J/mol/K
)
diff_sublithospheric_mantle = DiffusionCreep(
A=2.51e-17Pa^(-1 - 0) / s * m^(-0), # units are Pa^(-n - r) / s * m^(-p)
n=1.0,
E=530e3J/mol,
V=6e-6m^3/mol,
r=0.0,
p=0.0,
R=8.3145J/mol/K
)

# Elasticity
el_upper_crust = SetConstantElasticity(; G=25e9, ν=0.5)
el_lower_crust = SetConstantElasticity(; G=25e9, ν=0.5)
el_lithospheric_mantle = SetConstantElasticity(; G=67e9, ν=0.5)
el_sublithospheric_mantle = SetConstantElasticity(; G=67e9, ν=0.5)
el_upper_crust = SetConstantElasticity(; G=25e9Pa, ν=0.5)
el_lower_crust = SetConstantElasticity(; G=25e9Pa, ν=0.5)
el_lithospheric_mantle = SetConstantElasticity(; G=67e9Pa, ν=0.5)
el_sublithospheric_mantle = SetConstantElasticity(; G=67e9Pa, ν=0.5)
β_upper_crust = inv(get_Kb(el_upper_crust))
β_lower_crust = inv(get_Kb(el_lower_crust))
β_lithospheric_mantle = inv(get_Kb(el_lithospheric_mantle))
β_sublithospheric_mantle = inv(get_Kb(el_sublithospheric_mantle))

# Physical properties using GeoParams ----------------
η_reg = 1e16
cohesion = 3e6
η_reg = 1e16 * Pa * s
cohesion = 3e6 * Pa
friction = asind(0.2)
pl_crust = if is_plastic
DruckerPrager_regularised(; C = cohesion, ϕ=friction, η_vp=η_reg, Ψ=0.0) # non-regularized plasticity
Expand Down Expand Up @@ -61,20 +105,20 @@ function init_rheologies(CharDim; is_plastic = true)
# Name = "UpperCrust",
SetMaterialParams(;
Phase = 1,
Density = PT_Density(; ρ0=2.75e3, β=β_upper_crust, T0=0.0, α = 3.5e-5),
HeatCapacity = ConstantHeatCapacity(; Cp=7.5e2),
Density = PT_Density(; ρ0=2.75e3kg / m^3, β=β_upper_crust, T0=273K, α = 3.5e-5/ K),
HeatCapacity = ConstantHeatCapacity(; Cp=7.5e2J / kg / K),
Conductivity = K_crust,
CompositeRheology = CompositeRheology((disl_upper_crust, el_upper_crust, pl_crust)),
Elasticity = el_upper_crust,
RadioactiveHeat = ConstantRadioactiveHeat(0.0),
Gravity = ConstantGravity(; g=9.81),
Gravity = ConstantGravity(; g=9.81m/s^2),
CharDim = CharDim,
),
# Name = "LowerCrust",
SetMaterialParams(;
Phase = 2,
Density = PT_Density(; ρ0=3e3, β=β_lower_crust, T0=0.0, α = 3.5e-5),
HeatCapacity = ConstantHeatCapacity(; Cp=7.5e2),
Density = PT_Density(; ρ0=3e3kg / m^3, β=β_upper_crust, T0=273K, α = 3.5e-5/ K),
HeatCapacity = ConstantHeatCapacity(; Cp=7.5e2J / kg / K),
Conductivity = K_crust,
RadioactiveHeat = ConstantRadioactiveHeat(0.0),
CompositeRheology = CompositeRheology((disl_lower_crust, el_lower_crust, pl_crust)),
Expand All @@ -84,8 +128,8 @@ function init_rheologies(CharDim; is_plastic = true)
# Name = "LithosphericMantle",
SetMaterialParams(;
Phase = 3,
Density = PT_Density(; ρ0=3.3e3, β=β_lithospheric_mantle, T0=0.0, α = 3e-5),
HeatCapacity = ConstantHeatCapacity(; Cp=1.25e3),
Density = PT_Density(; ρ0=3.3e3kg / m^3, β=β_upper_crust, T0=273K, α = 3.5e-5/ K),
HeatCapacity = ConstantHeatCapacity(; Cp=1.25e3J / kg / K),
Conductivity = K_mantle,
RadioactiveHeat = ConstantRadioactiveHeat(0.0),
CompositeRheology = CompositeRheology((disl_lithospheric_mantle, diff_lithospheric_mantle, el_lithospheric_mantle, pl)),
Expand All @@ -94,8 +138,8 @@ function init_rheologies(CharDim; is_plastic = true)
),
SetMaterialParams(;
Phase = 4,
Density = PT_Density(; ρ0=3.3e3-50, β=β_sublithospheric_mantle, T0=0.0, α = 3e-5),
HeatCapacity = ConstantHeatCapacity(; Cp=1.25e3),
Density = PT_Density(; ρ0=(3.3e3-50)kg / m^3, β=β_upper_crust, T0=273K, α = 3.5e-5/ K),
HeatCapacity = ConstantHeatCapacity(; Cp=1.25e3J / kg / K),
Conductivity = K_mantle,
RadioactiveHeat = ConstantRadioactiveHeat(0.0),
CompositeRheology = CompositeRheology((disl_sublithospheric_mantle, diff_sublithospheric_mantle, el_sublithospheric_mantle)),
Expand All @@ -105,17 +149,16 @@ function init_rheologies(CharDim; is_plastic = true)
# Name = "StickyAir",
SetMaterialParams(;
Phase = 5,
Density = ConstantDensity(; ρ=1e3), # water density
HeatCapacity = ConstantHeatCapacity(; Cp=3e3),
Density = ConstantDensity(; ρ=1e3kg / m^3), # water density
HeatCapacity = ConstantHeatCapacity(; Cp=3e3J / kg / K),
RadioactiveHeat = ConstantRadioactiveHeat(0.0),
Conductivity = ConstantConductivity(; k=1.0),
CompositeRheology = CompositeRheology((LinearViscous(; η=1e19),)),
Conductivity = ConstantConductivity(; k=1.0Watt / m / K),
CompositeRheology = CompositeRheology((LinearViscous(; η=1e19Pa*s),)),
CharDim = CharDim,
),
)
end


function init_phases!(phases, particles, Lx, d, r, thick_air, CharDim)
ni = size(phases)

Expand Down
8 changes: 4 additions & 4 deletions src/rheology/Viscosity.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@
ε = εxx[I...], εyy[I...]

# we need strain rate not to be zero, otherwise we get NaNs
εII_0 = allzero(ε...) * 1e-15
εII_0 = allzero(ε...) * eps()

# argument fields at local index
args_ij = local_viscosity_args(args, I...)
Expand Down Expand Up @@ -59,7 +59,7 @@ end
ε = εxx[I...], εyy[I...]

# we need strain rate not to be zero, otherwise we get NaNs
εII_0 = allzero(ε...) * 1e-18
εII_0 = allzero(ε...) * eps()

# argument fields at local index
args_ij = local_viscosity_args(args, I...)
Expand Down Expand Up @@ -95,7 +95,7 @@ end
εij_normal = εxx[I...], εyy[I...], εzz[I...]

# we need strain rate not to be zero, otherwise we get NaNs
εII_0 = allzero(εij_normal...) * 1e-18
εII_0 = allzero(εij_normal...) * eps()

# # argument fields at local index
args_ijk = local_viscosity_args(args, I...)
Expand Down Expand Up @@ -145,7 +145,7 @@ end
εij_normal = εxx[I...], εyy[I...], εzz[I...]

# we need strain rate not to be zero, otherwise we get NaNs
εII_0 = allzero(εij_normal...) * 1e-18
εII_0 = allzero(εij_normal...) * eps()

# # argument fields at local index
args_ijk = local_viscosity_args(args, I...)
Expand Down

0 comments on commit 949e27b

Please sign in to comment.