Skip to content

Commit

Permalink
Merge branch 'main' into adm-freesurface
Browse files Browse the repository at this point in the history
  • Loading branch information
@aelligp
aelligp committed Feb 14, 2024
2 parents cf65a1c + 872c178 commit 55795a8
Show file tree
Hide file tree
Showing 25 changed files with 104 additions and 84 deletions.
2 changes: 1 addition & 1 deletion miniapps/benchmarks/stokes2D/shear_heating/Shearheating2D.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,7 @@ function main2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", save_vtk =false)

# Physical properties using GeoParams ----------------
rheology = init_rheologies(; is_TP_Conductivity=false)
κ = (4 / (rheology[1].HeatCapacity[1].cp * rheology[1].Density[1].ρ))
κ = (4 / (rheology[1].HeatCapacity[1].Cp * rheology[1].Density[1].ρ))
dt = dt_diff = 0.5 * min(di...)^2 / κ / 2.01 # diffusive CFL timestep limiter
# ----------------------------------------------------

Expand Down
4 changes: 2 additions & 2 deletions miniapps/benchmarks/stokes2D/shear_heating/Shearheating_rheology.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@ function init_rheologies(; is_TP_Conductivity=true)
SetMaterialParams(;
Phase = 1,
Density = ConstantDensity(ρ = 2700),
HeatCapacity = ConstantHeatCapacity(; cp=1050.0),
HeatCapacity = ConstantHeatCapacity(; Cp=1050.0),
Conductivity = K_Matrix,
ShearHeat = ConstantShearheating(1.0NoUnits),
CompositeRheology = CompositeRheology((Matrix, )),
Expand All @@ -45,7 +45,7 @@ function init_rheologies(; is_TP_Conductivity=true)
SetMaterialParams(;
Phase = 2,
Density = ConstantDensity(ρ = 2700),
HeatCapacity = ConstantHeatCapacity(; cp=1050.0),
HeatCapacity = ConstantHeatCapacity(; Cp=1050.0),
Conductivity = K_Inclusion,
ShearHeat = ConstantShearheating(1.0NoUnits),
CompositeRheology = CompositeRheology((Inclusion, )),
Expand Down
2 changes: 1 addition & 1 deletion miniapps/benchmarks/stokes3D/shear_heating/Shearheating3D.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,7 @@ function main3D(igg; ar=8, ny=16, nx=ny*8, figdir="figs3D", save_vtk =false)

# Physical properties using GeoParams ----------------
rheology = init_rheologies(; is_TP_Conductivity=false)
κ = (4 / (rheology[1].HeatCapacity[1].cp * rheology[1].Density[1].ρ))
κ = (4 / (rheology[1].HeatCapacity[1].Cp * rheology[1].Density[1].ρ))
dt = dt_diff = 0.5 * min(di...)^3 / κ / 3.01 # diffusive CFL timestep limiter
# ----------------------------------------------------

Expand Down
4 changes: 2 additions & 2 deletions miniapps/benchmarks/stokes3D/shear_heating/Shearheating_rheology.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@ function init_rheologies(; is_TP_Conductivity=true)
SetMaterialParams(;
Phase = 1,
Density = ConstantDensity(ρ = 2700),
HeatCapacity = ConstantHeatCapacity(; cp=1050.0),
HeatCapacity = ConstantHeatCapacity(; Cp=1050.0),
Conductivity = K_Matrix,
ShearHeat = ConstantShearheating(1.0NoUnits),
CompositeRheology = CompositeRheology((Matrix, )),
Expand All @@ -45,7 +45,7 @@ function init_rheologies(; is_TP_Conductivity=true)
SetMaterialParams(;
Phase = 2,
Density = ConstantDensity(ρ = 2700),
HeatCapacity = ConstantHeatCapacity(; cp=1050.0),
HeatCapacity = ConstantHeatCapacity(; Cp=1050.0),
Conductivity = K_Inclusion,
ShearHeat = ConstantShearheating(1.0NoUnits),
CompositeRheology = CompositeRheology((Inclusion, )),
Expand Down
30 changes: 15 additions & 15 deletions miniapps/benchmarks/thermal_diffusion/diffusion/diffusion2D.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,16 +9,16 @@
return nothing
end

function elliptical_perturbation!(T, δT, xc, yc, r, xvi)
@parallel_indices (i, j) function _elliptical_perturbation!(T, δT, xc, yc, r, x, y)
@inbounds if (((x[i]-xc ))^2 + ((y[j] - yc))^2) ≤ r^2
T[i, j] += δT
end
return nothing
end
@parallel _elliptical_perturbation!(T, δT, xc, yc, r, xvi...)
function elliptical_perturbation!(T, δT, xc, yc, r, xvi)

@parallel_indices (i, j) function _elliptical_perturbation!(T, δT, xc, yc, r, x, y)
@inbounds if (((x[i]-xc ))^2 + ((y[j] - yc))^2) ≤ r^2
T[i, j] += δT
end
return nothing
end

@parallel _elliptical_perturbation!(T, δT, xc, yc, r, xvi...)
end

function diffusion_2D(; nx=32, ny=32, lx=100e3, ly=100e3, ρ0=3.3e3, Cp0=1.2e3, K0=3.0)
Expand All @@ -32,14 +32,14 @@ function diffusion_2D(; nx=32, ny=32, lx=100e3, ly=100e3, ρ0=3.3e3, Cp0=1.2e3,
li = (lx, ly) # domain length in x- and y-
di = @. li / ni # grid step in x- and -y
origin = 0, -ly
grid = Geometry(ni, li; origin = origin)
grid = Geometry(ni, li; origin = origin)
(; xci, xvi) = grid # nodes at the center and vertices of the cells

# Define the thermal parameters with GeoParams
rheology = SetMaterialParams(;
Phase = 1,
Density = PT_Density(; ρ0=3.1e3, β=0.0, T0=0.0, α = 1.5e-5),
HeatCapacity = ConstantHeatCapacity(; cp=Cp0),
HeatCapacity = ConstantHeatCapacity(; Cp=Cp0),
Conductivity = ConstantConductivity(; k=K0),
)
# fields needed to compute density on the fly
Expand All @@ -56,8 +56,8 @@ function diffusion_2D(; nx=32, ny=32, lx=100e3, ly=100e3, ρ0=3.3e3, Cp0=1.2e3,
ρCp = @. Cp * ρ

pt_thermal = PTThermalCoeffs(K, ρCp, dt, di, li)
thermal_bc = TemperatureBoundaryConditions(;
no_flux = (left = true, right = true, top = false, bot = false),
thermal_bc = TemperatureBoundaryConditions(;
no_flux = (left = true, right = true, top = false, bot = false),
)
@parallel (@idx size(thermal.T)) init_T!(thermal.T, xvi[2])

Expand Down
2 changes: 1 addition & 1 deletion miniapps/benchmarks/thermal_diffusion/diffusion/diffusion2D_MPI.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,7 @@ function diffusion_2D(;
rheology = SetMaterialParams(;
Phase = 1,
Density = PT_Density(; ρ0=3.1e3, β=0.0, T0=0.0, α = 1.5e-5),
HeatCapacity = ConstantHeatCapacity(; cp=Cp0),
HeatCapacity = ConstantHeatCapacity(; Cp=Cp0),
Conductivity = ConstantConductivity(; k=K0),
)
# fields needed to compute density on the fly
Expand Down
4 changes: 2 additions & 2 deletions miniapps/benchmarks/thermal_diffusion/diffusion/diffusion2D_multiphase.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -94,14 +94,14 @@ function diffusion_2D(; nx=32, ny=32, lx=100e3, ly=100e3, Cp0=1.2e3, K0=3.0)
SetMaterialParams(;
Phase = 1,
Density = PT_Density(; ρ0=3e3, β=0.0, T0=0.0, α = 1.5e-5),
HeatCapacity = ConstantHeatCapacity(; cp=Cp0),
HeatCapacity = ConstantHeatCapacity(; Cp=Cp0),
Conductivity = ConstantConductivity(; k=K0),
RadioactiveHeat = ConstantRadioactiveHeat(1e-6),
),
SetMaterialParams(;
Phase = 1,
Density = PT_Density(; ρ0=3.3e3, β=0.0, T0=0.0, α = 1.5e-5),
HeatCapacity = ConstantHeatCapacity(; cp=Cp0),
HeatCapacity = ConstantHeatCapacity(; Cp=Cp0),
Conductivity = ConstantConductivity(; k=K0),
RadioactiveHeat = ConstantRadioactiveHeat(1e-7),
),
Expand Down
10 changes: 5 additions & 5 deletions miniapps/benchmarks/thermal_diffusion/diffusion/diffusion3D.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,14 +48,14 @@ function diffusion_3D(;
di = @. li / ni # grid step in x- and -y
origin = 0, 0, -lz # nodes at the center and vertices of the cells
igg = IGG(init_global_grid(nx, ny, nz; init_MPI=init_MPI)...) # init MPI
grid = Geometry(ni, li; origin = origin)
grid = Geometry(ni, li; origin = origin)
(; xci, xvi) = grid # nodes at the center and vertices of the cells

# Define the thermal parameters with GeoParams
rheology = SetMaterialParams(;
Phase = 1,
Density = PT_Density(; ρ0=3.1e3, β=0.0, T0=0.0, α = 1.5e-5),
HeatCapacity = ConstantHeatCapacity(; cp=Cp0),
HeatCapacity = ConstantHeatCapacity(; Cp=Cp0),
Conductivity = ConstantConductivity(; k=K0),
)

Expand All @@ -75,8 +75,8 @@ function diffusion_3D(;

# Boundary conditions
pt_thermal = PTThermalCoeffs(K, ρCp, dt, di, li; CFL = 0.75 / √3.1)
thermal_bc = TemperatureBoundaryConditions(;
no_flux = (left = true , right = true , top = false, bot = false, front = true , back = true),
thermal_bc = TemperatureBoundaryConditions(;
no_flux = (left = true , right = true , top = false, bot = false, front = true , back = true),
periodicity = (left = false, right = false, top = false, bot = false, front = false, back = false),
)

Expand Down Expand Up @@ -104,7 +104,7 @@ function diffusion_3D(;
di,;
igg
)

t += dt
it += 1
end
Expand Down
2 changes: 1 addition & 1 deletion miniapps/benchmarks/thermal_diffusion/diffusion/diffusion3D_MPI.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,7 @@ function diffusion_3D(;
rheology = SetMaterialParams(;
Phase = 1,
Density = PT_Density(; ρ0=3.1e3, β=0.0, T0=0.0, α = 1.5e-5),
HeatCapacity = ConstantHeatCapacity(; cp=Cp0),
HeatCapacity = ConstantHeatCapacity(; Cp=Cp0),
Conductivity = ConstantConductivity(; k=K0),
)

Expand Down
6 changes: 3 additions & 3 deletions miniapps/convection/GlobalConvection2D_Upwind.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -114,7 +114,7 @@ function thermal_convection2D(; ar=8, ny=16, nx=ny*8, figdir="figs2D", thermal_p
Name = "Mantle",
Phase = 1,
Density = PT_Density(; ρ0=3.1e3, β=β, T0=0.0, α = 1.5e-5),
HeatCapacity = ConstantHeatCapacity(; cp=1.2e3),
HeatCapacity = ConstantHeatCapacity(; Cp=1.2e3),
Conductivity = ConstantConductivity(; k=3.0),
CompositeRheology = CompositeRheology((creep, el, )),
Elasticity = el,
Expand All @@ -124,14 +124,14 @@ function thermal_convection2D(; ar=8, ny=16, nx=ny*8, figdir="figs2D", thermal_p
Name = "Mantle",
Phase = 1,
Density = PT_Density(; ρ0=3.5e3, β=β, T0=0.0, α = 1.5e-5),
HeatCapacity = ConstantHeatCapacity(; cp=1.2e3),
HeatCapacity = ConstantHeatCapacity(; Cp=1.2e3),
Conductivity = ConstantConductivity(; k=3.0),
CompositeRheology = CompositeRheology((creep, el, pl)),
Elasticity = el,
Gravity = ConstantGravity(; g=9.81),
)
# heat diffusivity
κ = (rheology.Conductivity[1].k / (rheology.HeatCapacity[1].cp * rheology.Density[1].ρ0)).val
κ = (rheology.Conductivity[1].k / (rheology.HeatCapacity[1].Cp * rheology.Density[1].ρ0)).val
dt = dt_diff = 0.5 * min(di...)^2 / κ / 2.01 # diffusive CFL timestep limiter
# ----------------------------------------------------

Expand Down
6 changes: 3 additions & 3 deletions miniapps/convection/GlobalConvection2D_WENO5.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -115,7 +115,7 @@ function thermal_convection2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", therma
Name = "Mantle",
Phase = 1,
Density = PT_Density(; ρ0=3.1e3, β=β, T0=0.0, α = 1.5e-5),
HeatCapacity = ConstantHeatCapacity(; cp=1.2e3),
HeatCapacity = ConstantHeatCapacity(; Cp=1.2e3),
Conductivity = ConstantConductivity(; k=3.0),
CompositeRheology = CompositeRheology((creep, el, )),
Elasticity = el,
Expand All @@ -125,14 +125,14 @@ function thermal_convection2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", therma
Name = "Mantle",
Phase = 1,
Density = PT_Density(; ρ0=3.5e3, β=β, T0=0.0, α = 1.5e-5),
HeatCapacity = ConstantHeatCapacity(; cp=1.2e3),
HeatCapacity = ConstantHeatCapacity(; Cp=1.2e3),
Conductivity = ConstantConductivity(; k=3.0),
CompositeRheology = CompositeRheology((creep, el, pl)),
Elasticity = el,
Gravity = ConstantGravity(; g=9.81),
)
# heat diffusivity
κ = (rheology.Conductivity[1].k / (rheology.HeatCapacity[1].cp * rheology.Density[1].ρ0)).val
κ = (rheology.Conductivity[1].k / (rheology.HeatCapacity[1].Cp * rheology.Density[1].ρ0)).val
dt = dt_diff = 0.5 * min(di...)^2 / κ / 2.01 # diffusive CFL timestep limiter
# ----------------------------------------------------

Expand Down
6 changes: 3 additions & 3 deletions miniapps/convection/GlobalConvection2D_WENO5_MPI.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -115,7 +115,7 @@ function thermal_convection2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", therma
Name = "Mantle",
Phase = 1,
Density = PT_Density(; ρ0=3.1e3, β=β, T0=0.0, α = 1.5e-5),
HeatCapacity = ConstantHeatCapacity(; cp=1.2e3),
HeatCapacity = ConstantHeatCapacity(; Cp=1.2e3),
Conductivity = ConstantConductivity(; k=3.0),
CompositeRheology = CompositeRheology((creep, el, )),
Elasticity = el,
Expand All @@ -125,14 +125,14 @@ function thermal_convection2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", therma
Name = "Mantle",
Phase = 1,
Density = PT_Density(; ρ0=3.5e3, β=β, T0=0.0, α = 1.5e-5),
HeatCapacity = ConstantHeatCapacity(; cp=1.2e3),
HeatCapacity = ConstantHeatCapacity(; Cp=1.2e3),
Conductivity = ConstantConductivity(; k=3.0),
CompositeRheology = CompositeRheology((creep, el, pl)),
Elasticity = el,
Gravity = ConstantGravity(; g=9.81),
)
# heat diffusivity
κ = (rheology.Conductivity[1].k / (rheology.HeatCapacity[1].cp * rheology.Density[1].ρ0)).val
κ = (rheology.Conductivity[1].k / (rheology.HeatCapacity[1].Cp * rheology.Density[1].ρ0)).val
dt = dt_diff = 0.5 * min(di...)^2 / κ / 2.01 # diffusive CFL timestep limiter
# ----------------------------------------------------

Expand Down
6 changes: 3 additions & 3 deletions miniapps/convection/GlobalConvection3D_Upwind.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -117,7 +117,7 @@ function thermal_convection3D(; ar=8, nz=16, nx=ny*8, ny=nx, figdir="figs3D", th
Name = "Mantle",
Phase = 1,
Density = PT_Density(; ρ0=3.1e3, β=β, T0=0.0, α = 1.5e-5),
HeatCapacity = ConstantHeatCapacity(; cp=1.2e3),
HeatCapacity = ConstantHeatCapacity(; Cp=1.2e3),
Conductivity = ConstantConductivity(; k=3.0),
CompositeRheology = CompositeRheology((creep, el, )),
Elasticity = el,
Expand All @@ -127,14 +127,14 @@ function thermal_convection3D(; ar=8, nz=16, nx=ny*8, ny=nx, figdir="figs3D", th
Name = "Mantle",
Phase = 1,
Density = PT_Density(; ρ0=3.5e3, β=β, T0=0.0, α = 1.5e-5),
HeatCapacity = ConstantHeatCapacity(; cp=1.2e3),
HeatCapacity = ConstantHeatCapacity(; Cp=1.2e3),
Conductivity = ConstantConductivity(; k=3.0),
CompositeRheology = CompositeRheology((creep, el, pl)),
Elasticity = el,
Gravity = ConstantGravity(; g=9.81),
)
# heat diffusivity
κ = (rheology.Conductivity[1].k / (rheology.HeatCapacity[1].cp * rheology.Density[1].ρ0)).val
κ = (rheology.Conductivity[1].k / (rheology.HeatCapacity[1].Cp * rheology.Density[1].ρ0)).val
dt = dt_diff = min(di...)^2 / κ / 3.01 # diffusive CFL timestep limiter
# ----------------------------------------------------

Expand Down
2 changes: 1 addition & 1 deletion miniapps/convection/Particles2D/Layered_convection2D.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -109,7 +109,7 @@ function main2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", save_vtk =false)

# Physical properties using GeoParams ----------------
rheology = init_rheologies(; is_plastic = true)
κ = (4 / (rheology[4].HeatCapacity[1].cp * rheology[4].Density[1].ρ0))
κ = (4 / (rheology[4].HeatCapacity[1].Cp * rheology[4].Density[1].ρ0))
dt = dt_diff = 0.5 * min(di...)^2 / κ / 2.01 # diffusive CFL timestep limiter
# ----------------------------------------------------

Expand Down
18 changes: 9 additions & 9 deletions miniapps/convection/Particles2D/Layered_rheology.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,18 +24,18 @@ function init_rheologies(; is_plastic = true)
η_reg = 1e16
cohesion = 3e6
friction = asind(0.2)
pl_crust = if is_plastic
pl_crust = if is_plastic
DruckerPrager_regularised(; C = cohesion, ϕ=friction, η_vp=η_reg, Ψ=0.0) # non-regularized plasticity
else
DruckerPrager_regularised(; C = Inf, ϕ=friction, η_vp=η_reg, Ψ=0.0) # non-regularized plasticity
end
friction = asind(0.3)
pl = if is_plastic
pl = if is_plastic
DruckerPrager_regularised(; C = cohesion, ϕ=friction, η_vp=η_reg, Ψ=0.0) # non-regularized plasticity
else
DruckerPrager_regularised(; C = Inf, ϕ=friction, η_vp=η_reg, Ψ=0.0) # non-regularized plasticity
end
pl_wz = if is_plastic
pl_wz = if is_plastic
DruckerPrager_regularised(; C = 2e6, ϕ=2.0, η_vp=η_reg, Ψ=0.0) # non-regularized plasticity
else
DruckerPrager_regularised(; C = Inf, ϕ=friction, η_vp=η_reg, Ψ=0.0) # non-regularized plasticity
Expand All @@ -62,7 +62,7 @@ function init_rheologies(; is_plastic = true)
SetMaterialParams(;
Phase = 1,
Density = PT_Density(; ρ0=2.75e3, β=β_upper_crust, T0=0.0, α = 3.5e-5),
HeatCapacity = ConstantHeatCapacity(; cp=7.5e2),
HeatCapacity = ConstantHeatCapacity(; Cp=7.5e2),
Conductivity = K_crust,
CompositeRheology = CompositeRheology((disl_upper_crust, el_upper_crust, pl_crust)),
Elasticity = el_upper_crust,
Expand All @@ -73,7 +73,7 @@ function init_rheologies(; is_plastic = true)
SetMaterialParams(;
Phase = 2,
Density = PT_Density(; ρ0=3e3, β=β_lower_crust, T0=0.0, α = 3.5e-5),
HeatCapacity = ConstantHeatCapacity(; cp=7.5e2),
HeatCapacity = ConstantHeatCapacity(; Cp=7.5e2),
Conductivity = K_crust,
RadioactiveHeat = ConstantRadioactiveHeat(0.0),
CompositeRheology = CompositeRheology((disl_lower_crust, el_lower_crust, pl_crust)),
Expand All @@ -83,7 +83,7 @@ function init_rheologies(; is_plastic = true)
SetMaterialParams(;
Phase = 3,
Density = PT_Density(; ρ0=3.3e3, β=β_lithospheric_mantle, T0=0.0, α = 3e-5),
HeatCapacity = ConstantHeatCapacity(; cp=1.25e3),
HeatCapacity = ConstantHeatCapacity(; Cp=1.25e3),
Conductivity = K_mantle,
RadioactiveHeat = ConstantRadioactiveHeat(0.0),
CompositeRheology = CompositeRheology((disl_lithospheric_mantle, diff_lithospheric_mantle, el_lithospheric_mantle, pl)),
Expand All @@ -92,7 +92,7 @@ function init_rheologies(; 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),
HeatCapacity = ConstantHeatCapacity(; Cp=1.25e3),
Conductivity = K_mantle,
RadioactiveHeat = ConstantRadioactiveHeat(0.0),
CompositeRheology = CompositeRheology((disl_sublithospheric_mantle, diff_sublithospheric_mantle, el_sublithospheric_mantle)),
Expand All @@ -102,7 +102,7 @@ function init_rheologies(; is_plastic = true)
SetMaterialParams(;
Phase = 5,
Density = ConstantDensity(; ρ=1e3), # water density
HeatCapacity = ConstantHeatCapacity(; cp=3e3),
HeatCapacity = ConstantHeatCapacity(; Cp=3e3),
RadioactiveHeat = ConstantRadioactiveHeat(0.0),
Conductivity = ConstantConductivity(; k=1.0),
CompositeRheology = CompositeRheology((LinearViscous(; η=1e19),)),
Expand Down Expand Up @@ -132,7 +132,7 @@ function init_phases!(phases, particles, Lx, d, r, thick_air)
elseif depth > 90e3
@cell phases[ip, i, j] = 3.0

elseif depth < 0e0
elseif depth < 0e0
@cell phases[ip, i, j] = 5.0

end
Expand Down
Loading

0 comments on commit 55795a8

Please sign in to comment.