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] 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)