Skip to content

Commit

Permalink
Patch to adapt to changes done in GP#main (#104)
Browse files Browse the repository at this point in the history
* patch to adapt to changes done in GP#main

* format

* format
  • Loading branch information
albert-de-montserrat authored Feb 14, 2024
1 parent e8c2f36 commit 872c178
Show file tree
Hide file tree
Showing 6 changed files with 29 additions and 9 deletions.
7 changes: 5 additions & 2 deletions src/Utils.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
15 changes: 15 additions & 0 deletions src/rheology/GeoParams.jl
Original file line number Diff line number Diff line change
@@ -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
4 changes: 2 additions & 2 deletions src/stokes/PressureKernels.jl
Original file line number Diff line number Diff line change
Expand Up @@ -19,15 +19,15 @@ 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

@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
Expand Down
1 change: 1 addition & 0 deletions src/stokes/Stokes2D.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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")
Expand Down
3 changes: 2 additions & 1 deletion src/stokes/Stokes3D.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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")
Expand Down Expand Up @@ -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...)

Expand Down
8 changes: 4 additions & 4 deletions src/stokes/StressKernels.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down Expand Up @@ -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)

Expand Down

0 comments on commit 872c178

Please sign in to comment.