From cac378a7b07963203e0dd12527a4dfc4f6f93e56 Mon Sep 17 00:00:00 2001 From: Iddingsite Date: Tue, 2 Dec 2025 12:13:44 +0100 Subject: [PATCH 1/7] pressure_kernel --- src/stokes/PressureKernels.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/stokes/PressureKernels.jl b/src/stokes/PressureKernels.jl index fb1556f15..1d0ed22c6 100644 --- a/src/stokes/PressureKernels.jl +++ b/src/stokes/PressureKernels.jl @@ -178,7 +178,7 @@ function _compute_P!(P, P0, ∇V, Q, η, K, dt, r, θ_dτ) _Kdt = inv(K * dt) _dt = inv(dt) RP = muladd(-(P - P0), _Kdt, (-∇V + (Q * _dt))) - ψ = inv(inv(r / θ_dτ * η) + _Kdt) + ψ = inv(inv(η) + _Kdt) * r / θ_dτ P = ((muladd(P0, _Kdt, (-∇V + (Q * _dt)))) * ψ + P) / (1 + _Kdt * ψ) return RP, P @@ -188,7 +188,7 @@ function _compute_P!(P, P0, ∇V, Q, ΔTc, α, η, K, dt, r, θ_dτ) _Kdt = inv(K * dt) _dt = inv(dt) RP = muladd(-(P - P0), _Kdt, (-∇V + (α * (ΔTc * _dt)) + (Q * _dt))) - ψ = inv(inv(r / θ_dτ * η) + _Kdt) + ψ = inv(inv(η) + _Kdt) * r / θ_dτ P = ((muladd(P0, _Kdt, (-∇V + (α * (ΔTc * _dt)) + (Q * _dt)))) * ψ + P) / (1 + _Kdt * ψ) return RP, P From 394e0b40fdafce41cd2f6336b78f62bbe51181ea Mon Sep 17 00:00:00 2001 From: Iddingsite Date: Tue, 2 Dec 2025 15:28:08 +0100 Subject: [PATCH 2/7] correct expression --- src/stokes/PressureKernels.jl | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/src/stokes/PressureKernels.jl b/src/stokes/PressureKernels.jl index 1d0ed22c6..166f0e593 100644 --- a/src/stokes/PressureKernels.jl +++ b/src/stokes/PressureKernels.jl @@ -95,7 +95,8 @@ end ::Nothing, ) where {N, C <: JustRelax.CellArray} K = fn_ratio(get_bulk_modulus, rheology, @cell(phase_ratio[I...])) - @inbounds RP[I...], P[I...] = _compute_P!(P[I...], P0[I...], ∇V[I...], Q[I...], η[I...], K, dt, r, θ_dτ) + Gdt = fn_ratio(get_shear_modulus, rheology, @cell(phase_ratio[I...])) * dt + @inbounds RP[I...], P[I...] = _compute_P!(P[I...], P0[I...], ∇V[I...], Q[I...], η[I...], K, Gdt, dt, r, θ_dτ) return nothing end @@ -115,7 +116,8 @@ end melt_fraction, ) where {N, C <: JustRelax.CellArray} K = fn_ratio(get_bulk_modulus, rheology, @cell(phase_ratio[I...])) - @inbounds RP[I...], P[I...] = _compute_P!(P[I...], P0[I...], ∇V[I...], Q[I...], η[I...], K, dt, r, θ_dτ) + Gdt = fn_ratio(get_shear_modulus, rheology, @cell(phase_ratio[I...])) * dt + @inbounds RP[I...], P[I...] = _compute_P!(P[I...], P0[I...], ∇V[I...], Q[I...], η[I...], K, Gdt, dt, r, θ_dτ) return nothing end @@ -136,9 +138,10 @@ end ) where {N, C <: JustRelax.CellArray} @inbounds phase_ratio_I = phase_ratio[I...] @inbounds K = fn_ratio(get_bulk_modulus, rheology, phase_ratio_I) + @inbounds Gdt = fn_ratio(get_shear_modulus, rheology, phase_ratio_I) * dt @inbounds α = fn_ratio(get_thermal_expansion, rheology, phase_ratio_I) @inbounds RP[I...], P[I...] = _compute_P!( - P[I...], P0[I...], ∇V[I...], Q[I...], ΔTc[I...], α, η[I...], K, dt, r, θ_dτ + P[I...], P0[I...], ∇V[I...], Q[I...], ΔTc[I...], α, η[I...], K, Gdt, dt, r, θ_dτ ) return nothing end @@ -159,9 +162,10 @@ end melt_fraction, ) where {N, C <: JustRelax.CellArray} @inbounds K = fn_ratio(get_bulk_modulus, rheology, @cell(phase_ratio[I...])) + @inbounds Gdt = fn_ratio(get_shear_modulus, rheology, @cell(phase_ratio[I...])) * dt @inbounds α = fn_ratio(get_thermal_expansion, rheology, @cell(phase_ratio[I...]), (; ϕ = melt_fraction[I...])) @inbounds RP[I...], P[I...] = _compute_P!( - P[I...], P0[I...], ∇V[I...], Q[I...], ΔTc[I...], α, η[I...], K, dt, r, θ_dτ + P[I...], P0[I...], ∇V[I...], Q[I...], ΔTc[I...], α, η[I...], K, Gdt, dt, r, θ_dτ ) return nothing end @@ -174,21 +178,21 @@ function _compute_P!(P, ∇V, Q, η, dt, r, θ_dτ) return RP, P end -function _compute_P!(P, P0, ∇V, Q, η, K, dt, r, θ_dτ) +function _compute_P!(P, P0, ∇V, Q, η, K, Gdt, dt, r, θ_dτ) _Kdt = inv(K * dt) _dt = inv(dt) RP = muladd(-(P - P0), _Kdt, (-∇V + (Q * _dt))) - ψ = inv(inv(η) + _Kdt) * r / θ_dτ + ψ = inv(inv(η) + inv(Gdt)) * r / θ_dτ P = ((muladd(P0, _Kdt, (-∇V + (Q * _dt)))) * ψ + P) / (1 + _Kdt * ψ) return RP, P end -function _compute_P!(P, P0, ∇V, Q, ΔTc, α, η, K, dt, r, θ_dτ) +function _compute_P!(P, P0, ∇V, Q, ΔTc, α, η, K, Gdt, dt, r, θ_dτ) _Kdt = inv(K * dt) _dt = inv(dt) RP = muladd(-(P - P0), _Kdt, (-∇V + (α * (ΔTc * _dt)) + (Q * _dt))) - ψ = inv(inv(η) + _Kdt) * r / θ_dτ + ψ = inv(inv(η) + inv(Gdt)) * r / θ_dτ P = ((muladd(P0, _Kdt, (-∇V + (α * (ΔTc * _dt)) + (Q * _dt)))) * ψ + P) / (1 + _Kdt * ψ) return RP, P From 4fb2560bf5b4599cd344ad416d7b446a1d360ce4 Mon Sep 17 00:00:00 2001 From: Iddingsite Date: Tue, 2 Dec 2025 15:46:42 +0100 Subject: [PATCH 3/7] small fix --- src/stokes/PressureKernels.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/stokes/PressureKernels.jl b/src/stokes/PressureKernels.jl index 166f0e593..84c1e1cde 100644 --- a/src/stokes/PressureKernels.jl +++ b/src/stokes/PressureKernels.jl @@ -7,9 +7,9 @@ end ## Compressible -@parallel_indices (I...) function compute_P!(P, P0, RP, ∇V, Q, η, K, dt, r, θ_dτ) +@parallel_indices (I...) function compute_P!(P, P0, RP, ∇V, Q, η, K, Gdt, dt, r, θ_dτ) @inbounds RP[I...], P[I...] = _compute_P!( - P[I...], P0[I...], ∇V[I...], Q[I...], η[I...], K[I...], dt, r, θ_dτ + P[I...], P0[I...], ∇V[I...], Q[I...], η[I...], K[I...], Gdt[I...], dt, r, θ_dτ ) return nothing end @@ -20,7 +20,8 @@ end P, P0, RP, ∇V, Q, η, rheology::NTuple{N, MaterialParams}, phase, dt, r, θ_dτ, args ) where {N} K = get_bulk_modulus(rheology, phase[I...]) - @inbounds RP[I...], P[I...] = _compute_P!(P[I...], P0[I...], ∇V[I...], Q[I...], η[I...], K, dt, r, θ_dτ) + Gdt = get_shear_modulus(rheology, phase[I...]) * dt + @inbounds RP[I...], P[I...] = _compute_P!(P[I...], P0[I...], ∇V[I...], Q[I...], η[I...], K, Gdt, dt, r, θ_dτ) return nothing end From 239b6ad96c180bd707ced5bbdabba9fdfc78fb2a Mon Sep 17 00:00:00 2001 From: Iddingsite Date: Tue, 2 Dec 2025 16:09:49 +0100 Subject: [PATCH 4/7] more fix --- src/stokes/Stokes2D.jl | 2 +- src/stokes/Stokes3D.jl | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/src/stokes/Stokes2D.jl b/src/stokes/Stokes2D.jl index 3af252bd8..b1ef62ca8 100644 --- a/src/stokes/Stokes2D.jl +++ b/src/stokes/Stokes2D.jl @@ -203,7 +203,7 @@ function _solve!( @parallel (@idx ni) compute_∇V!(stokes.∇V, @velocity(stokes), _di) @parallel compute_P!( - stokes.P, stokes.P0, stokes.R.RP, stokes.∇V, stokes.Q, ητ, K, dt, r, θ_dτ + stokes.P, stokes.P0, stokes.R.RP, stokes.∇V, stokes.Q, ητ, K, G*dt, dt, r, θ_dτ ) @parallel (@idx ni .+ 1) compute_strain_rate!( diff --git a/src/stokes/Stokes3D.jl b/src/stokes/Stokes3D.jl index 661e0301b..95ceab0da 100644 --- a/src/stokes/Stokes3D.jl +++ b/src/stokes/Stokes3D.jl @@ -81,6 +81,7 @@ function _solve!( stokes.Q, η, K, + G*dt, dt, pt_stokes.r, pt_stokes.θ_dτ, From 418bbd8bcd865f55307515840ff39896118aa8bb Mon Sep 17 00:00:00 2001 From: Iddingsite Date: Tue, 2 Dec 2025 17:28:22 +0100 Subject: [PATCH 5/7] pass G and not Gdt --- src/stokes/PressureKernels.jl | 36 ++++++++++++++++++----------------- src/stokes/Stokes2D.jl | 2 +- src/stokes/Stokes3D.jl | 2 +- 3 files changed, 21 insertions(+), 19 deletions(-) diff --git a/src/stokes/PressureKernels.jl b/src/stokes/PressureKernels.jl index 84c1e1cde..9907486e4 100644 --- a/src/stokes/PressureKernels.jl +++ b/src/stokes/PressureKernels.jl @@ -7,21 +7,21 @@ end ## Compressible -@parallel_indices (I...) function compute_P!(P, P0, RP, ∇V, Q, η, K, Gdt, dt, r, θ_dτ) +@parallel_indices (I...) function compute_P!(P, P0, RP, ∇V, Q, η, K, G, dt, r, θ_dτ) @inbounds RP[I...], P[I...] = _compute_P!( - P[I...], P0[I...], ∇V[I...], Q[I...], η[I...], K[I...], Gdt[I...], dt, r, θ_dτ + P[I...], P0[I...], ∇V[I...], Q[I...], η[I...], K[I...], G[I...], dt, r, θ_dτ ) return nothing end -# With GeoParams + # With GeoParams @parallel_indices (I...) function compute_P!( P, P0, RP, ∇V, Q, η, rheology::NTuple{N, MaterialParams}, phase, dt, r, θ_dτ, args ) where {N} K = get_bulk_modulus(rheology, phase[I...]) - Gdt = get_shear_modulus(rheology, phase[I...]) * dt - @inbounds RP[I...], P[I...] = _compute_P!(P[I...], P0[I...], ∇V[I...], Q[I...], η[I...], K, Gdt, dt, r, θ_dτ) + G = get_shear_modulus(rheology, phase[I...]) + @inbounds RP[I...], P[I...] = _compute_P!(P[I...], P0[I...], ∇V[I...], Q[I...], η[I...], K, G, dt, r, θ_dτ) return nothing end @@ -96,8 +96,8 @@ end ::Nothing, ) where {N, C <: JustRelax.CellArray} K = fn_ratio(get_bulk_modulus, rheology, @cell(phase_ratio[I...])) - Gdt = fn_ratio(get_shear_modulus, rheology, @cell(phase_ratio[I...])) * dt - @inbounds RP[I...], P[I...] = _compute_P!(P[I...], P0[I...], ∇V[I...], Q[I...], η[I...], K, Gdt, dt, r, θ_dτ) + G = fn_ratio(get_shear_modulus, rheology, @cell(phase_ratio[I...])) + @inbounds RP[I...], P[I...] = _compute_P!(P[I...], P0[I...], ∇V[I...], Q[I...], η[I...], K, G, dt, r, θ_dτ) return nothing end @@ -117,8 +117,8 @@ end melt_fraction, ) where {N, C <: JustRelax.CellArray} K = fn_ratio(get_bulk_modulus, rheology, @cell(phase_ratio[I...])) - Gdt = fn_ratio(get_shear_modulus, rheology, @cell(phase_ratio[I...])) * dt - @inbounds RP[I...], P[I...] = _compute_P!(P[I...], P0[I...], ∇V[I...], Q[I...], η[I...], K, Gdt, dt, r, θ_dτ) + G = fn_ratio(get_shear_modulus, rheology, @cell(phase_ratio[I...])) + @inbounds RP[I...], P[I...] = _compute_P!(P[I...], P0[I...], ∇V[I...], Q[I...], η[I...], K, G, dt, r, θ_dτ) return nothing end @@ -139,10 +139,10 @@ end ) where {N, C <: JustRelax.CellArray} @inbounds phase_ratio_I = phase_ratio[I...] @inbounds K = fn_ratio(get_bulk_modulus, rheology, phase_ratio_I) - @inbounds Gdt = fn_ratio(get_shear_modulus, rheology, phase_ratio_I) * dt + @inbounds G = fn_ratio(get_shear_modulus, rheology, phase_ratio_I) @inbounds α = fn_ratio(get_thermal_expansion, rheology, phase_ratio_I) @inbounds RP[I...], P[I...] = _compute_P!( - P[I...], P0[I...], ∇V[I...], Q[I...], ΔTc[I...], α, η[I...], K, Gdt, dt, r, θ_dτ + P[I...], P0[I...], ∇V[I...], Q[I...], ΔTc[I...], α, η[I...], K, G, dt, r, θ_dτ ) return nothing end @@ -163,10 +163,10 @@ end melt_fraction, ) where {N, C <: JustRelax.CellArray} @inbounds K = fn_ratio(get_bulk_modulus, rheology, @cell(phase_ratio[I...])) - @inbounds Gdt = fn_ratio(get_shear_modulus, rheology, @cell(phase_ratio[I...])) * dt + @inbounds G = fn_ratio(get_shear_modulus, rheology, @cell(phase_ratio[I...])) @inbounds α = fn_ratio(get_thermal_expansion, rheology, @cell(phase_ratio[I...]), (; ϕ = melt_fraction[I...])) @inbounds RP[I...], P[I...] = _compute_P!( - P[I...], P0[I...], ∇V[I...], Q[I...], ΔTc[I...], α, η[I...], K, Gdt, dt, r, θ_dτ + P[I...], P0[I...], ∇V[I...], Q[I...], ΔTc[I...], α, η[I...], K, G, dt, r, θ_dτ ) return nothing end @@ -179,21 +179,23 @@ function _compute_P!(P, ∇V, Q, η, dt, r, θ_dτ) return RP, P end -function _compute_P!(P, P0, ∇V, Q, η, K, Gdt, dt, r, θ_dτ) +function _compute_P!(P, P0, ∇V, Q, η, K, G, dt, r, θ_dτ) _Kdt = inv(K * dt) + _Gdt = inv(G * dt) _dt = inv(dt) RP = muladd(-(P - P0), _Kdt, (-∇V + (Q * _dt))) - ψ = inv(inv(η) + inv(Gdt)) * r / θ_dτ + ψ = inv(inv(η) + _Gdt) * r / θ_dτ P = ((muladd(P0, _Kdt, (-∇V + (Q * _dt)))) * ψ + P) / (1 + _Kdt * ψ) return RP, P end -function _compute_P!(P, P0, ∇V, Q, ΔTc, α, η, K, Gdt, dt, r, θ_dτ) +function _compute_P!(P, P0, ∇V, Q, ΔTc, α, η, K, G, dt, r, θ_dτ) _Kdt = inv(K * dt) + _Gdt = inv(G * dt) _dt = inv(dt) RP = muladd(-(P - P0), _Kdt, (-∇V + (α * (ΔTc * _dt)) + (Q * _dt))) - ψ = inv(inv(η) + inv(Gdt)) * r / θ_dτ + ψ = inv(inv(η) + _Gdt) * r / θ_dτ P = ((muladd(P0, _Kdt, (-∇V + (α * (ΔTc * _dt)) + (Q * _dt)))) * ψ + P) / (1 + _Kdt * ψ) return RP, P diff --git a/src/stokes/Stokes2D.jl b/src/stokes/Stokes2D.jl index b1ef62ca8..467daef09 100644 --- a/src/stokes/Stokes2D.jl +++ b/src/stokes/Stokes2D.jl @@ -203,7 +203,7 @@ function _solve!( @parallel (@idx ni) compute_∇V!(stokes.∇V, @velocity(stokes), _di) @parallel compute_P!( - stokes.P, stokes.P0, stokes.R.RP, stokes.∇V, stokes.Q, ητ, K, G*dt, dt, r, θ_dτ + stokes.P, stokes.P0, stokes.R.RP, stokes.∇V, stokes.Q, ητ, K, G, dt, r, θ_dτ ) @parallel (@idx ni .+ 1) compute_strain_rate!( diff --git a/src/stokes/Stokes3D.jl b/src/stokes/Stokes3D.jl index 95ceab0da..a3fdee728 100644 --- a/src/stokes/Stokes3D.jl +++ b/src/stokes/Stokes3D.jl @@ -81,7 +81,7 @@ function _solve!( stokes.Q, η, K, - G*dt, + G, dt, pt_stokes.r, pt_stokes.θ_dτ, From f769c98b08641a8d8546a62d4c4bf33086266218 Mon Sep 17 00:00:00 2001 From: Iddingsite Date: Tue, 2 Dec 2025 17:29:33 +0100 Subject: [PATCH 6/7] runic --- src/stokes/PressureKernels.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stokes/PressureKernels.jl b/src/stokes/PressureKernels.jl index 9907486e4..23410927a 100644 --- a/src/stokes/PressureKernels.jl +++ b/src/stokes/PressureKernels.jl @@ -14,7 +14,7 @@ end return nothing end - # With GeoParams +# With GeoParams @parallel_indices (I...) function compute_P!( P, P0, RP, ∇V, Q, η, rheology::NTuple{N, MaterialParams}, phase, dt, r, θ_dτ, args From d65bfbcf422cccc2b5de8df4ef655ac1f9bb0f4f Mon Sep 17 00:00:00 2001 From: Iddingsite Date: Tue, 2 Dec 2025 17:44:20 +0100 Subject: [PATCH 7/7] last fix? --- src/rheology/GeoParams.jl | 6 +++--- src/stokes/Stokes2D.jl | 3 ++- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/src/rheology/GeoParams.jl b/src/rheology/GeoParams.jl index ca65c7401..e8cfad774 100644 --- a/src/rheology/GeoParams.jl +++ b/src/rheology/GeoParams.jl @@ -7,11 +7,11 @@ function get_bulk_modulus(args::Vararg{Any, N}) where {N} end function get_shear_modulus(args::Vararg{Any, N}) where {N} - Kb = GeoParams.get_G(args...) - if isnan(Kb) || iszero(Kb) + G = GeoParams.get_G(args...) + if isnan(G) || iszero(G) return Inf end - return Kb + return G end get_thermal_expansion(args::Vararg{Any, N}) where {N} = get_α(args...) diff --git a/src/stokes/Stokes2D.jl b/src/stokes/Stokes2D.jl index 467daef09..81213852e 100644 --- a/src/stokes/Stokes2D.jl +++ b/src/stokes/Stokes2D.jl @@ -328,6 +328,7 @@ function _solve!( # end Kb = get_Kb(rheology) + G = get_G(rheology) # errors err_it1 = 1.0 @@ -368,7 +369,7 @@ function _solve!( @parallel (@idx ni) compute_∇V!(stokes.∇V, @velocity(stokes), _di) @parallel compute_P!( - stokes.P, stokes.P0, stokes.R.RP, stokes.∇V, stokes.Q, η, Kb, dt, r, θ_dτ + stokes.P, stokes.P0, stokes.R.RP, stokes.∇V, stokes.Q, η, Kb, G, dt, r, θ_dτ ) update_ρg!(ρg[2], rheology, args)