From 949e27bda7835c4af7f03448a1dd39dc6537e4d8 Mon Sep 17 00:00:00 2001 From: Albert de Montserrat Date: Tue, 26 Mar 2024 17:58:55 +0100 Subject: [PATCH] more non dim --- .../Layered_convection2D.jl | 27 +++--- .../Particles2D_nonDim/Layered_rheology.jl | 95 ++++++++++++++----- src/rheology/Viscosity.jl | 8 +- 3 files changed, 88 insertions(+), 42 deletions(-) diff --git a/miniapps/convection/Particles2D_nonDim/Layered_convection2D.jl b/miniapps/convection/Particles2D_nonDim/Layered_convection2D.jl index dbdb0222..36341292 100644 --- a/miniapps/convection/Particles2D_nonDim/Layered_convection2D.jl +++ b/miniapps/convection/Particles2D_nonDim/Layered_convection2D.jl @@ -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 @@ -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 -------------------------------------------------------------- @@ -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 -------------------------------- @@ -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) # ---------------------------------------------------- @@ -392,4 +395,4 @@ end # run main script -main2D(igg; figdir = figdir, ar = ar, nx = nx, ny = ny, do_vtk = do_vtk); \ No newline at end of file +# main2D(igg; figdir = figdir, ar = ar, nx = nx, ny = ny, do_vtk = do_vtk); \ No newline at end of file diff --git a/miniapps/convection/Particles2D_nonDim/Layered_rheology.jl b/miniapps/convection/Particles2D_nonDim/Layered_rheology.jl index cea9572e..527fe95f 100644 --- a/miniapps/convection/Particles2D_nonDim/Layered_rheology.jl +++ b/miniapps/convection/Particles2D_nonDim/Layered_rheology.jl @@ -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 @@ -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)), @@ -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)), @@ -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)), @@ -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) diff --git a/src/rheology/Viscosity.jl b/src/rheology/Viscosity.jl index d9089086..e07ca4c1 100644 --- a/src/rheology/Viscosity.jl +++ b/src/rheology/Viscosity.jl @@ -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...) @@ -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...) @@ -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...) @@ -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...)