fix P bug

PTsolvers · Dec 2, 2024 · f5efe36 · f5efe36
1 parent f084316
commit f5efe36
Show file tree

Hide file tree

Showing 2 changed files with 28 additions and 28 deletions.
diff --git a/src/stokes/PressureKernels.jl b/src/stokes/PressureKernels.jl
@@ -64,7 +64,7 @@ function compute_P!(
 ) where {N}
     ni = size(P)
     @parallel (@idx ni) compute_P_kernel!(
-        P, P0, RP, ∇V, η, rheology, phase_ratio, dt, r, θ_dτ, ΔTc, melt_fraction
+        P, P0, RP, ∇V, η, rheology, phase_ratio.center, dt, r, θ_dτ, ΔTc, melt_fraction
     )
     return nothing
 end

diff --git a/src/variational_stokes/PressureKernels.jl b/src/variational_stokes/PressureKernels.jl
@@ -45,30 +45,30 @@ end
     return nothing
 end
 
-@parallel_indices (I...) function compute_P_kernel!(
-    P,
-    P0,
-    RP,
-    ∇V,
-    η,
-    rheology::NTuple{N,MaterialParams},
-    phase_ratio::C,
-    ϕ::JustRelax.RockRatio,
-    dt,
-    r,
-    θ_dτ,
-    ΔTc,
-    ::Nothing,
-) where {N,C<:JustRelax.CellArray}
-    if isvalid_c(ϕ, I...)
-        phase_ratio_I = phase_ratio[I...]
-        K = fn_ratio(get_bulk_modulus, rheology, phase_ratio_I)
-        α = fn_ratio(get_thermal_expansion, rheology, phase_ratio_I)
-        RP[I...], P[I...] = _compute_P!(
-            P[I...], P0[I...], ∇V[I...], ΔTc[I...], α, η[I...], K, dt, r, θ_dτ
-        )
-    else
-        RP[I...] = P[I...] = zero(eltype(P))
-    end
-    return nothing
-end
+# @parallel_indices (I...) function compute_P_kernel!(
+#     P,
+#     P0,
+#     RP,
+#     ∇V,
+#     η,
+#     rheology::NTuple{N,MaterialParams},
+#     phase_ratio::C,
+#     ϕ::JustRelax.RockRatio,
+#     dt,
+#     r,
+#     θ_dτ,
+#     ΔTc,
+#     ::Nothing,
+# ) where {N,C<:JustRelax.CellArray}
+#     if isvalid_c(ϕ, I...)
+#         phase_ratio_I = phase_ratio[I...]
+#         K = fn_ratio(get_bulk_modulus, rheology, phase_ratio_I)
+#         α = fn_ratio(get_thermal_expansion, rheology, phase_ratio_I)
+#         RP[I...], P[I...] = _compute_P!(
+#             P[I...], P0[I...], ∇V[I...], ΔTc[I...], α, η[I...], K, dt, r, θ_dτ
+#         )
+#     else
+#         RP[I...] = P[I...] = zero(eltype(P))
+#     end
+#     return nothing
+# end