From e8c2f36b0d9ef1bb5053c28e36c269ec666f6aa6 Mon Sep 17 00:00:00 2001 From: Pascal Aellig <93140290+aelligp@users.noreply.github.com> Date: Wed, 14 Feb 2024 15:28:42 +0100 Subject: [PATCH 1/2] rename cp (#102) --- .../stokes2D/shear_heating/Shearheating2D.jl | 2 +- .../shear_heating/Shearheating_rheology.jl | 4 +-- .../stokes3D/shear_heating/Shearheating3D.jl | 2 +- .../shear_heating/Shearheating_rheology.jl | 4 +-- .../diffusion/diffusion2D.jl | 30 +++++++++---------- .../diffusion/diffusion2D_MPI.jl | 2 +- .../diffusion/diffusion2D_multiphase.jl | 4 +-- .../diffusion/diffusion3D.jl | 10 +++---- .../diffusion/diffusion3D_MPI.jl | 2 +- .../convection/GlobalConvection2D_Upwind.jl | 6 ++-- .../convection/GlobalConvection2D_WENO5.jl | 6 ++-- .../GlobalConvection2D_WENO5_MPI.jl | 6 ++-- .../convection/GlobalConvection3D_Upwind.jl | 6 ++-- .../Particles2D/Layered_convection2D.jl | 2 +- .../Particles2D/Layered_rheology.jl | 18 +++++------ .../Particles3D/Layered_convection3D.jl | 2 +- .../Particles3D/Layered_rheology.jl | 22 +++++++------- miniapps/convection/WENO5/Layered_rheology.jl | 20 ++++++------- .../convection/WENO5/WENO_convection2D.jl | 2 +- 19 files changed, 75 insertions(+), 75 deletions(-) diff --git a/miniapps/benchmarks/stokes2D/shear_heating/Shearheating2D.jl b/miniapps/benchmarks/stokes2D/shear_heating/Shearheating2D.jl index 530407bd..7e81bd14 100644 --- a/miniapps/benchmarks/stokes2D/shear_heating/Shearheating2D.jl +++ b/miniapps/benchmarks/stokes2D/shear_heating/Shearheating2D.jl @@ -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 # ---------------------------------------------------- diff --git a/miniapps/benchmarks/stokes2D/shear_heating/Shearheating_rheology.jl b/miniapps/benchmarks/stokes2D/shear_heating/Shearheating_rheology.jl index 6eb20053..74b5694f 100644 --- a/miniapps/benchmarks/stokes2D/shear_heating/Shearheating_rheology.jl +++ b/miniapps/benchmarks/stokes2D/shear_heating/Shearheating_rheology.jl @@ -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, )), @@ -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, )), diff --git a/miniapps/benchmarks/stokes3D/shear_heating/Shearheating3D.jl b/miniapps/benchmarks/stokes3D/shear_heating/Shearheating3D.jl index 177ee532..16b5e0e3 100644 --- a/miniapps/benchmarks/stokes3D/shear_heating/Shearheating3D.jl +++ b/miniapps/benchmarks/stokes3D/shear_heating/Shearheating3D.jl @@ -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 # ---------------------------------------------------- diff --git a/miniapps/benchmarks/stokes3D/shear_heating/Shearheating_rheology.jl b/miniapps/benchmarks/stokes3D/shear_heating/Shearheating_rheology.jl index 535b9a52..5041a86a 100644 --- a/miniapps/benchmarks/stokes3D/shear_heating/Shearheating_rheology.jl +++ b/miniapps/benchmarks/stokes3D/shear_heating/Shearheating_rheology.jl @@ -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, )), @@ -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, )), diff --git a/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion2D.jl b/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion2D.jl index ec01bc77..020feec2 100644 --- a/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion2D.jl +++ b/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion2D.jl @@ -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) @@ -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 @@ -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]) diff --git a/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion2D_MPI.jl b/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion2D_MPI.jl index 096bb71a..add48909 100644 --- a/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion2D_MPI.jl +++ b/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion2D_MPI.jl @@ -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 diff --git a/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion2D_multiphase.jl b/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion2D_multiphase.jl index c77d7ad6..1f366f8f 100644 --- a/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion2D_multiphase.jl +++ b/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion2D_multiphase.jl @@ -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), ), diff --git a/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion3D.jl b/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion3D.jl index a3fb114f..e5e5cbbb 100644 --- a/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion3D.jl +++ b/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion3D.jl @@ -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), ) @@ -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), ) @@ -104,7 +104,7 @@ function diffusion_3D(; di,; igg ) - + t += dt it += 1 end diff --git a/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion3D_MPI.jl b/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion3D_MPI.jl index f2bdfbfe..a802c5d1 100644 --- a/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion3D_MPI.jl +++ b/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion3D_MPI.jl @@ -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), ) diff --git a/miniapps/convection/GlobalConvection2D_Upwind.jl b/miniapps/convection/GlobalConvection2D_Upwind.jl index b501b168..45d8725e 100644 --- a/miniapps/convection/GlobalConvection2D_Upwind.jl +++ b/miniapps/convection/GlobalConvection2D_Upwind.jl @@ -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, @@ -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 # ---------------------------------------------------- diff --git a/miniapps/convection/GlobalConvection2D_WENO5.jl b/miniapps/convection/GlobalConvection2D_WENO5.jl index 4f239f97..f78513ec 100644 --- a/miniapps/convection/GlobalConvection2D_WENO5.jl +++ b/miniapps/convection/GlobalConvection2D_WENO5.jl @@ -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, @@ -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 # ---------------------------------------------------- diff --git a/miniapps/convection/GlobalConvection2D_WENO5_MPI.jl b/miniapps/convection/GlobalConvection2D_WENO5_MPI.jl index 756e434e..56f0873e 100644 --- a/miniapps/convection/GlobalConvection2D_WENO5_MPI.jl +++ b/miniapps/convection/GlobalConvection2D_WENO5_MPI.jl @@ -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, @@ -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 # ---------------------------------------------------- diff --git a/miniapps/convection/GlobalConvection3D_Upwind.jl b/miniapps/convection/GlobalConvection3D_Upwind.jl index bee653e5..e1a20083 100644 --- a/miniapps/convection/GlobalConvection3D_Upwind.jl +++ b/miniapps/convection/GlobalConvection3D_Upwind.jl @@ -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, @@ -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 # ---------------------------------------------------- diff --git a/miniapps/convection/Particles2D/Layered_convection2D.jl b/miniapps/convection/Particles2D/Layered_convection2D.jl index 41b0f81b..28d383e3 100644 --- a/miniapps/convection/Particles2D/Layered_convection2D.jl +++ b/miniapps/convection/Particles2D/Layered_convection2D.jl @@ -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 # ---------------------------------------------------- diff --git a/miniapps/convection/Particles2D/Layered_rheology.jl b/miniapps/convection/Particles2D/Layered_rheology.jl index f4596b2a..c5a6f6e6 100644 --- a/miniapps/convection/Particles2D/Layered_rheology.jl +++ b/miniapps/convection/Particles2D/Layered_rheology.jl @@ -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 @@ -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, @@ -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)), @@ -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)), @@ -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)), @@ -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),)), @@ -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 diff --git a/miniapps/convection/Particles3D/Layered_convection3D.jl b/miniapps/convection/Particles3D/Layered_convection3D.jl index dc81c68b..90e9c28e 100644 --- a/miniapps/convection/Particles3D/Layered_convection3D.jl +++ b/miniapps/convection/Particles3D/Layered_convection3D.jl @@ -94,7 +94,7 @@ function main3D(igg; ar=1, nx=16, ny=16, nz=16, figdir="figs3D", do_vtk =false) # Physical properties using GeoParams ---------------- rheology = init_rheologies(; is_plastic = true) - κ = (10 / (rheology[1].HeatCapacity[1].cp * rheology[1].Density[1].ρ0)) + κ = (10 / (rheology[1].HeatCapacity[1].Cp * rheology[1].Density[1].ρ0)) dt = dt_diff = 0.5 * min(di...)^3 / κ / 3.01 # diffusive CFL timestep limiter # ---------------------------------------------------- diff --git a/miniapps/convection/Particles3D/Layered_rheology.jl b/miniapps/convection/Particles3D/Layered_rheology.jl index 96f49cf2..bff86b86 100644 --- a/miniapps/convection/Particles3D/Layered_rheology.jl +++ b/miniapps/convection/Particles3D/Layered_rheology.jl @@ -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 @@ -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, @@ -72,7 +72,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, CompositeRheology = CompositeRheology((disl_lower_crust, el_lower_crust, pl_crust)), Elasticity = el_lower_crust, @@ -81,7 +81,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, CompositeRheology = CompositeRheology((disl_lithospheric_mantle, diff_lithospheric_mantle, el_lithospheric_mantle, pl)), Elasticity = el_lithospheric_mantle, @@ -90,7 +90,7 @@ function init_rheologies(; is_plastic = true) # SetMaterialParams(; # Phase = 4, # Density = PT_Density(; ρ0=3.3e3, β=β_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)), @@ -100,7 +100,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, CompositeRheology = CompositeRheology((disl_sublithospheric_mantle, diff_sublithospheric_mantle, el_sublithospheric_mantle)), Elasticity = el_sublithospheric_mantle, @@ -109,7 +109,7 @@ function init_rheologies(; is_plastic = true) SetMaterialParams(; Phase = 5, Density = ConstantDensity(; ρ=1e3), - HeatCapacity = ConstantHeatCapacity(; cp=1.25e3), + HeatCapacity = ConstantHeatCapacity(; Cp=1.25e3), Conductivity = ConstantConductivity(; k=15.0), CompositeRheology = CompositeRheology((LinearViscous(; η=1e21),)), ), @@ -120,7 +120,7 @@ function init_phases!(phases, particles, Lx, Ly; d=650e3, r=50e3) ni = size(phases) @parallel_indices (I...) function init_phases!(phases, px, py, pz, index, r, Lx, Ly) - + @inbounds for ip in JustRelax.JustRelax.cellaxes(phases) # quick escape JustRelax.@cell(index[ip, I...]) == 0 && continue @@ -141,7 +141,7 @@ function init_phases!(phases, particles, Lx, Ly; d=650e3, r=50e3) elseif depth > 90e3 @cell phases[ip, I...] = 3.0 - elseif 0e0 > depth + elseif 0e0 > depth @cell phases[ip, I...] = 5.0 end diff --git a/miniapps/convection/WENO5/Layered_rheology.jl b/miniapps/convection/WENO5/Layered_rheology.jl index 2c935162..bf85511a 100644 --- a/miniapps/convection/WENO5/Layered_rheology.jl +++ b/miniapps/convection/WENO5/Layered_rheology.jl @@ -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 @@ -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, @@ -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)), @@ -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)), @@ -93,7 +93,7 @@ function init_rheologies(; is_plastic = true) # SetMaterialParams(; # Phase = 4, # Density = PT_Density(; ρ0=3.3e3, β=β_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)), @@ -103,7 +103,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)), @@ -113,7 +113,7 @@ function init_rheologies(; is_plastic = true) SetMaterialParams(; Phase = 5, Density = ConstantDensity(; ρ=1e3), - HeatCapacity = ConstantHeatCapacity(; cp=1.25e3), + HeatCapacity = ConstantHeatCapacity(; Cp=1.25e3), RadioactiveHeat = ConstantRadioactiveHeat(0.0), Conductivity = ConstantConductivity(; k=15.0), CompositeRheology = CompositeRheology((LinearViscous(; η=1e21),)), @@ -143,7 +143,7 @@ function init_phases!(phases, particles, Lx; d=650e3, r=50e3) elseif depth > 90e3 @cell phases[ip, i, j] = 3.0 - elseif depth < 0e0 + elseif depth < 0e0 @cell phases[ip, i, j] = 5.0 end diff --git a/miniapps/convection/WENO5/WENO_convection2D.jl b/miniapps/convection/WENO5/WENO_convection2D.jl index 7100d581..36b91803 100644 --- a/miniapps/convection/WENO5/WENO_convection2D.jl +++ b/miniapps/convection/WENO5/WENO_convection2D.jl @@ -97,7 +97,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) - κ = (10 / (rheology[1].HeatCapacity[1].cp * rheology[1].Density[1].ρ0)) + κ = (10 / (rheology[1].HeatCapacity[1].Cp * rheology[1].Density[1].ρ0)) dt = dt_diff = 0.5 * min(di...)^2 / κ / 2.01 # diffusive CFL timestep limiter # ---------------------------------------------------- From 872c1787b766fbb4273cf64987918a303555654b Mon Sep 17 00:00:00 2001 From: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> Date: Wed, 14 Feb 2024 15:36:54 +0100 Subject: [PATCH 2/2] Patch to adapt to changes done in GP#main (#104) * patch to adapt to changes done in GP#main * format * format --- src/Utils.jl | 7 +++++-- src/rheology/GeoParams.jl | 15 +++++++++++++++ src/stokes/PressureKernels.jl | 4 ++-- src/stokes/Stokes2D.jl | 1 + src/stokes/Stokes3D.jl | 3 ++- src/stokes/StressKernels.jl | 8 ++++---- 6 files changed, 29 insertions(+), 9 deletions(-) create mode 100644 src/rheology/GeoParams.jl diff --git a/src/Utils.jl b/src/Utils.jl index 15d71c48..109a52d9 100644 --- a/src/Utils.jl +++ b/src/Utils.jl @@ -461,8 +461,11 @@ take(fldr::String) = !isdir(fldr) && mkpath(fldr) Do a continuation step `exp((1-ν)*log(x_old) + ν*log(x_new))` with damping parameter `ν` """ -# @inline continuation_log(x_new, x_old, ν) = exp((1 - ν) * log(x_old) + ν * log(x_new)) -@inline continuation_log(x_new, x_old, ν) = muladd((1 - ν), x_old, ν * x_new) # (1 - ν) * x_old + ν * x_new +@inline function continuation_log(x_new::T, x_old::T, ν) where {T} + x_cont = exp((1 - ν) * log(x_old) + ν * log(x_new)) + return isnan(x_cont) ? 0.0 : x_cont +end +# @inline continuation_log(x_new, x_old, ν) = muladd((1 - ν), x_old, ν * x_new) # (1 - ν) * x_old + ν * x_new # Others diff --git a/src/rheology/GeoParams.jl b/src/rheology/GeoParams.jl new file mode 100644 index 00000000..90d68eda --- /dev/null +++ b/src/rheology/GeoParams.jl @@ -0,0 +1,15 @@ +function get_bulk_modulus(args...) + Kb = GeoParams.get_Kb(args...) + if isnan(Kb) || iszero(Kb) + return Inf + end + return Kb +end + +function get_shear_modulus(args...) + Kb = GeoParams.get_G(args...) + if isnan(Kb) || iszero(Kb) + return Inf + end + return Kb +end diff --git a/src/stokes/PressureKernels.jl b/src/stokes/PressureKernels.jl index 6621b1d9..3adcd350 100644 --- a/src/stokes/PressureKernels.jl +++ b/src/stokes/PressureKernels.jl @@ -19,7 +19,7 @@ end @parallel_indices (I...) function compute_P!( P, P0, RP, ∇V, η, rheology::NTuple{N,MaterialParams}, phase, dt, r, θ_dτ ) where {N} - K = get_Kb(rheology, phase[I...]) + K = get_bulk_modulus(rheology, phase[I...]) RP[I...], P[I...] = _compute_P!(P[I...], P0[I...], ∇V[I...], η[I...], K, dt, r, θ_dτ) return nothing end @@ -27,7 +27,7 @@ end @parallel_indices (I...) function compute_P!( P, P0, RP, ∇V, η, rheology::NTuple{N,MaterialParams}, phase_ratio::C, dt, r, θ_dτ ) where {N,C<:JustRelax.CellArray} - K = fn_ratio(get_Kb, rheology, phase_ratio[I...]) + K = fn_ratio(get_bulk_modulus, rheology, phase_ratio[I...]) RP[I...], P[I...] = _compute_P!(P[I...], P0[I...], ∇V[I...], η[I...], K, dt, r, θ_dτ) return nothing end diff --git a/src/stokes/Stokes2D.jl b/src/stokes/Stokes2D.jl index 22d9c9df..434e2a0a 100644 --- a/src/stokes/Stokes2D.jl +++ b/src/stokes/Stokes2D.jl @@ -53,6 +53,7 @@ import JustRelax: mean_mpi, norm_mpi, maximum_mpi, minimum_mpi, backend @eval @init_parallel_stencil($backend, Float64, 2) +include("../rheology/GeoParams.jl") include("StressRotation.jl") include("PressureKernels.jl") include("VelocityKernels.jl") diff --git a/src/stokes/Stokes3D.jl b/src/stokes/Stokes3D.jl index 98f34f89..c76f1a5a 100644 --- a/src/stokes/Stokes3D.jl +++ b/src/stokes/Stokes3D.jl @@ -20,6 +20,7 @@ import JustRelax: mean_mpi, norm_mpi, minimum_mpi, maximum_mpi, backend @eval @init_parallel_stencil($backend, Float64, 3) +include("../rheology/GeoParams.jl") include("StressRotation.jl") include("PressureKernels.jl") include("VelocityKernels.jl") @@ -237,7 +238,7 @@ function JustRelax.solve!( norm_∇V = Float64[] Kb = get_Kb(rheology) - G = get_G(rheology) + G = get_shear_modulus(rheology) @copy stokes.P0 stokes.P λ = @zeros(ni...) diff --git a/src/stokes/StressKernels.jl b/src/stokes/StressKernels.jl index 55c9ff51..e5d22e23 100644 --- a/src/stokes/StressKernels.jl +++ b/src/stokes/StressKernels.jl @@ -276,14 +276,14 @@ end # numerics ηij = η[I...] - _Gdt = inv(get_G(rheology[1]) * dt) + _Gdt = inv(get_shear_modulus(rheology[1]) * dt) dτ_r = compute_dτ_r(θ_dτ, ηij, _Gdt) # get plastic parameters (if any...) is_pl, C, sinϕ, cosϕ, sinψ, η_reg = plastic_params_phase(rheology, EII[I...], 1) # plastic volumetric change K * dt * sinϕ * sinψ - K = get_bulkmodulus(rheology[1]) + K = get_bulk_modulus(rheology[1]) volume = isinf(K) ? 0.0 : K * dt * sinϕ * sinψ plastic_parameters = (; is_pl, C, sinϕ, cosϕ, η_reg, volume) @@ -318,14 +318,14 @@ end # numerics ηij = @inbounds η[I...] phase = @inbounds phase_center[I...] - _Gdt = inv(fn_ratio(get_G, rheology, phase) * dt) + _Gdt = inv(fn_ratio(get_shear_modulus, rheology, phase) * dt) dτ_r = compute_dτ_r(θ_dτ, ηij, _Gdt) # get plastic parameters (if any...) is_pl, C, sinϕ, cosϕ, sinψ, η_reg = plastic_params_phase(rheology, EII[I...], phase) # plastic volumetric change K * dt * sinϕ * sinψ - K = fn_ratio(get_bulkmodulus, rheology, phase) + K = fn_ratio(get_bulk_modulus, rheology, phase) volume = isinf(K) ? 0.0 : K * dt * sinϕ * sinψ plastic_parameters = (; is_pl, C, sinϕ, cosϕ, η_reg, volume)