Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add perturbation kwarg to cohesion #202

Merged
merged 1 commit into from
Jun 14, 2024
Merged

Add perturbation kwarg to cohesion #202

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
19 changes: 12 additions & 7 deletions src/Plasticity/DruckerPrager.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -55,33 +55,37 @@ end

# Calculation routines
function (s::DruckerPrager{_T, U, U1, S, S})(;
P=zero(_T), τII=zero(_T), Pf=zero(_T), EII=zero(_T), kwargs...
P=0.0, τII=0.0, Pf=0.0, EII=0.0, perturbation_C = 1.0, kwargs...
) where {_T,U,U1,S<:AbstractSoftening}
@unpack_val sinϕ, cosϕ, ϕ, C = s
ϕ = s.softening_ϕ(EII, ϕ)
C = s.softening_C(EII, C)
C *= perturbation_C

cosϕ, sinϕ = iszero(EII) ? (cosϕ, sinϕ) : (cosd(ϕ), sind(ϕ))

F = τII - cosϕ * C - sinϕ * (P - Pf) # with fluid pressure (set to zero by default)
return F
end


function (s::DruckerPrager{_T, U, U1, NoSoftening, S})(;
P=zero(_T), τII=zero(_T), Pf=zero(_T), EII=zero(_T), kwargs...
P=0.0, τII=0.0, Pf=0.0, EII=0.0, perturbation_C = 1.0, kwargs...
) where {_T,U,U1,S}
@unpack_val sinϕ, cosϕ, ϕ, C = s
C = s.softening_C(EII, C)
C *= perturbation_C

F = τII - cosϕ * C - sinϕ * (P - Pf) # with fluid pressure (set to zero by default)
return F
end

function (s::DruckerPrager{_T, U, U1, S, NoSoftening})(;
P=zero(_T), τII=zero(_T), Pf=zero(_T), EII=zero(_T), kwargs...
P=0.0, τII=0.0, Pf=0.0, EII=0.0, perturbation_C = 1.0, kwargs...
) where {_T,U,U1,S}
@unpack_val sinϕ, cosϕ, ϕ, C = s
ϕ = s.softening_ϕ(EII, ϕ)
C *= perturbation_C

cosϕ, sinϕ = iszero(EII) ? (cosϕ, sinϕ) : (cosd(ϕ), sind(ϕ))

Expand All @@ -90,9 +94,10 @@ function (s::DruckerPrager{_T, U, U1, S, NoSoftening})(;
end

function (s::DruckerPrager{_T, U, U1, NoSoftening, NoSoftening})(;
P=zero(_T), τII=zero(_T), Pf=zero(_T), kwargs...
P=0.0, τII=0.0, Pf=0.0, perturbation_C = 1.0, kwargs...
) where {_T,U,U1}
@unpack_val sinϕ, cosϕ, ϕ, C = s
C *= perturbation_C

F = τII - cosϕ * C - sinϕ * (P - Pf) # with fluid pressure (set to zero by default)
return F
Expand All @@ -104,9 +109,9 @@ end
Computes the plastic yield function `F` for a given second invariant of the deviatoric stress tensor `τII`, `P` pressure, and `Pf` fluid pressure.
"""
function compute_yieldfunction(
s::DruckerPrager{_T}; P=zero(_T), τII=zero(_T), Pf=zero(_T), EII=zero(_T)
s::DruckerPrager{_T}; P=0.0, τII=0.0, Pf=0.0, EII=0.0, perturbation_C = 1.0
) where {_T}
return s(; P=P, τII=τII, Pf=Pf, EII=EII)
return s(; P=P, τII=τII, Pf=Pf, EII=EII, perturbation_C =perturbation_C)
end

"""
Expand Down Expand Up @@ -135,7 +140,7 @@ end
# Plastic Potential

# Derivatives w.r.t pressure
∂Q∂P(p::DruckerPrager, P=zero(_T); τII=zero(_T), kwargs...) = -NumValue(p.sinΨ)
∂Q∂P(p::DruckerPrager, P=0.0; τII=0.0, kwargs...) = -NumValue(p.sinΨ)

# Derivatives of yield function
∂F∂τII(p::DruckerPrager, τII::_T; P=zero(_T), kwargs...) where _T = _T(1)
Expand Down
13 changes: 6 additions & 7 deletions src/Plasticity/DruckerPrager_regularised.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -56,7 +56,7 @@

# Calculation routines
function (s::DruckerPrager_regularised{_T})(;
P=zero(_T), τII=zero(_T), Pf=zero(_T), λ= zero(_T), EII=zero(_T), kwargs...
P=0.0, τII=0.0, Pf=0.0, λ= 0.0, EII=0.0, perturbation_C = 1.0, kwargs...
) where {_T}
@unpack_val sinϕ, cosϕ, ϕ, C, η_vp = s
ϕ = s.softening_ϕ(EII, ϕ)
Expand All @@ -68,7 +68,7 @@
end

function (s::DruckerPrager_regularised{_T,U,U1,NoSoftening, S})(;
P=zero(_T), τII=zero(_T), Pf=zero(_T), λ= zero(_T), EII=zero(_T), kwargs...
P=0.0, τII=0.0, Pf=0.0, λ= 0.0, EII=0.0, perturbation_C = 1.0, kwargs...
) where {_T,U,U1,S}
@unpack_val sinϕ, cosϕ, ϕ, C, η_vp = s
C = s.softening_C(EII, C)
Expand All @@ -78,7 +78,7 @@
end

function (s::DruckerPrager_regularised{_T,U,U1,S, NoSoftening})(;
P=zero(_T), τII=zero(_T), Pf=zero(_T), λ= zero(_T), EII=zero(_T), kwargs...
P=0.0, τII=0.0, Pf=0.0, λ= 0.0, EII=0.0, perturbation_C = 1.0, kwargs...
) where {_T,U,U1,S}
@unpack_val sinϕ, cosϕ, ϕ, C, η_vp = s
ϕ = s.softening_ϕ(EII, ϕ)
Expand All @@ -89,7 +89,7 @@
end

function (s::DruckerPrager_regularised{_T,U,U1,NoSoftening,NoSoftening})(;
P=zero(_T), τII=zero(_T), Pf=zero(_T), λ= zero(_T), kwargs...
P=0.0, τII=0.0, Pf=0.0, λ= 0.0, perturbation_C = 1.0, kwargs...
) where {_T,U,U1}
@unpack_val sinϕ, cosϕ, ϕ, C, η_vp = s
ε̇II_pl = λ*∂Q∂τII(s, τII) # plastic strainrate
Expand All @@ -105,9 +105,9 @@
Computes the plastic yield function `F` for a given second invariant of the deviatoric stress tensor `τII`, `P` pressure, and `Pf` fluid pressure.
"""
function compute_yieldfunction(
s::DruckerPrager_regularised{_T}; P=zero(_T), τII=zero(_T), Pf=zero(_T), λ=zero(_T), EII=zero(_T)
s::DruckerPrager_regularised{_T}; P=0.0, τII=0.0, Pf=0.0, λ=0.0, EII=0.0, perturbation_C = 1.0
) where {_T}
return s(; P=P, τII=τII, Pf=Pf, λ=λ, EII=EII)
return s(; P=P, τII=τII, Pf=Pf, λ=λ, EII=EII, perturbation_C = perturbation_C)

Check warning on line 110 in src/Plasticity/DruckerPrager_regularised.jl

View check run for this annotation

Codecov / codecov/patch

src/Plasticity/DruckerPrager_regularised.jl#L110 

Added line #L110 was not covered by tests
end

"""
Expand Down Expand Up @@ -144,7 +144,6 @@
∂F∂P(p::DruckerPrager_regularised, P::_T; kwargs...) where _T = -NumValue(p.sinϕ)
∂F∂λ(p::DruckerPrager_regularised, τII::_T; P=zero(_T), kwargs...) where _T = -2 * NumValue(p.η_vp)*∂Q∂τII(p, τII, P=P)


# Derivatives w.r.t stress tensor

# Hard-coded partial derivatives of the plastic potential Q
Expand Down
102 changes: 50 additions & 52 deletions test/test_Plasticity.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,32 +7,38 @@ using GeoParams
CharUnits_GEO = GEO_units(; viscosity=1e19, length=10km)

# DruckerPrager ---------
p = DruckerPrager()
p = DruckerPrager()
info = param_info(p)
@test isbits(p)
@test NumValue(p.ϕ) == 30
@test NumValue(p.ϕ) == 30
@test isvolumetric(p) == false

p_nd = p
p_nd = nondimensionalize(p_nd, CharUnits_GEO)
@test p_nd.C.val ≈ 1

# Compute with dimensional units
τII = 20e6
P = 1e6
args = (P=P, τII=τII)
args1 = (τII=τII, P=P)
@test compute_yieldfunction(p, args) ≈ 1.0839745962155614e7 # no Pfluid
@test compute_yieldfunction(p, args) ≈ compute_yieldfunction(p, args1) # different order
τII = 20e6
P = 1e6
args = (P = P, τII = τII)
args1 = (τII = τII, P = P)
@test compute_yieldfunction(p, args) ≈ 1.0839745962155614e7 # no Pfluid
@test compute_yieldfunction(p, args) ≈ compute_yieldfunction(p, args1) # different order

args_f = (P=P, τII=τII, Pf=0.5e6)

# test cohesion perturbation
@test compute_yieldfunction(p, (τII = τII, P = P, perturbation_C = 1.0)) ≈ 1.0839745962155614e7 # no perturbation
@test compute_yieldfunction(p, (τII = τII, P = P, perturbation_C = 0.0)) == 1.95e7 # cohesion = 0.0
@test compute_yieldfunction(p, (τII = τII, P = P, perturbation_C = 0.5)) ≈ 1.5169872981077807e7 # midway

args_f = (P = P, τII = τII, Pf = 0.5e6)

@test compute_yieldfunction(p, args_f) ≈ 1.1089745962155614e7 # with Pfluid

# Test with arrays
P_array = ones(10) * 1e6
τII_array = ones(10) * 20e6
F_array = similar(P_array)
P_array = fill(1e6, 10)
τII_array = fill(20e6, 10)
F_array = similar(P_array)
compute_yieldfunction!(F_array, p, (; P=P_array, τII=τII_array))
@test F_array[1] ≈ 1.0839745962155614e7

Expand All @@ -51,16 +57,16 @@ using GeoParams
)

# test computing material properties
n = 100
Phases = ones(Int64, n, n, n)
n = 100
Phases = ones(Int64, n, n, n)
Phases[:, :, 20:end] .= 2
Phases[:, :, 60:end] .= 2

τII = ones(size(Phases)) * 10e6
P = ones(size(Phases)) * 1e6
Pf = ones(size(Phases)) * 0.5e6
F = zero(P)
args = (P=P, τII=τII)
τII = fill( 10e6, size(Phases)...)
P = fill( 1e6, size(Phases)...)
Pf = fill(0.5e6, size(Phases)...)
F = zero(P)
args = (P = P, τII = τII)
compute_yieldfunction!(F, MatParam, Phases, args) # computation routine w/out P (not used in most heat capacity formulations)
@test maximum(F[1, 1, :]) ≈ 839745.962155614

Expand All @@ -73,8 +79,8 @@ using GeoParams
# test if we provide phase ratios
PhaseRatio = zeros(n, n, n, 3)
for i in CartesianIndices(Phases)
iz = Phases[i]
I = CartesianIndex(i, iz)
iz = Phases[i]
I = CartesianIndex(i, iz)
PhaseRatio[I] = 1.0
end
compute_yieldfunction!(F, MatParam, PhaseRatio, args)
Expand All @@ -84,10 +90,10 @@ using GeoParams

# Test plastic potential derivatives
## 2D
τij = (1.0, 2.0, 3.0)
fxx(τij) = 0.5 * τij[1] / second_invariant(τij)
fyy(τij) = 0.5 * τij[2] / second_invariant(τij)
fxy(τij) = τij[3] / second_invariant(τij)
τij = (1.0, 2.0, 3.0)
fxx(τij) = 0.5 * τij[1] / second_invariant(τij)
fyy(τij) = 0.5 * τij[2] / second_invariant(τij)
fxy(τij) = τij[3] / second_invariant(τij)
solution2D = [fxx(τij), fyy(τij), fxy(τij)]

# # using StaticArrays
Expand All @@ -98,27 +104,24 @@ using GeoParams

# using tuples
τij_tuple = (1.0, 2.0, 3.0)
out2 = ∂Q∂τ(p, τij_tuple)
out2 = ∂Q∂τ(p, τij_tuple)
@test out2 == Tuple(solution2D)
@test compute_plasticpotentialDerivative(p, τij_tuple) == ∂Q∂τ(p, τij_tuple)

# using AD
Q = second_invariant # where second_invariant is a function
# ad1 = ∂Q∂τ(Q, τij_static)
# @test out1 == solution2D
# @test compute_plasticpotentialDerivative(p, τij_static) == ∂Q∂τ(p, τij_static)
Q = second_invariant # where second_invariant is a function
ad2 = ∂Q∂τ(Q, τij_tuple)
@test out2 == Tuple(solution2D)
@test compute_plasticpotentialDerivative(p, τij_tuple) == ∂Q∂τ(p, τij_tuple)

## 3D
τij = (1.0, 2.0, 3.0, 4.0, 5.0, 6.0)
gxx(τij) = 0.5 * τij[1] / second_invariant(τij)
gyy(τij) = 0.5 * τij[2] / second_invariant(τij)
gzz(τij) = 0.5 * τij[3] / second_invariant(τij)
gyz(τij) = τij[4] / second_invariant(τij)
gxz(τij) = τij[5] / second_invariant(τij)
gxy(τij) = τij[6] / second_invariant(τij)
τij = (1.0, 2.0, 3.0, 4.0, 5.0, 6.0)
gxx(τij) = 0.5 * τij[1] / second_invariant(τij)
gyy(τij) = 0.5 * τij[2] / second_invariant(τij)
gzz(τij) = 0.5 * τij[3] / second_invariant(τij)
gyz(τij) = τij[4] / second_invariant(τij)
gxz(τij) = τij[5] / second_invariant(τij)
gxy(τij) = τij[6] / second_invariant(τij)
solution3D = [gxx(τij), gyy(τij), gzz(τij), gyz(τij), gxz(τij), gxy(τij)]

# # using StaticArrays
Expand All @@ -129,42 +132,37 @@ using GeoParams

# using tuples
τij_tuple = (1.0, 2.0, 3.0, 4.0, 5.0, 6.0)
out4 = ∂Q∂τ(p, τij_tuple)
out4 = ∂Q∂τ(p, τij_tuple)
@test out4 == Tuple(solution3D)
@test compute_plasticpotentialDerivative(p, τij_tuple) == ∂Q∂τ(p, τij_tuple)

# using AD
Q = second_invariant # where second_invariant is a function
# ad3 = ∂Q∂τ(Q, τij_static)
# @test out3 == solution3D
# @test compute_plasticpotentialDerivative(p, τij_static) == ∂Q∂τ(p, τij_static)
Q = second_invariant # where second_invariant is a function
ad4 = ∂Q∂τ(Q, τij_tuple)
@test out4 == Tuple(solution3D)
@test compute_plasticpotentialDerivative(p, τij_tuple) == ∂Q∂τ(p, τij_tuple)

# -----------------------

# composite rheology with plasticity
η, G = 10, 1
t_M = η / G
εII = 1.0
args = (;)
pl2 = DruckerPrager(; C=η, ϕ=0) # plasticity
c_pl = CompositeRheology(
LinearViscous(; η=η * Pa * s), ConstantElasticity(; G=G * Pa), pl2
) # linear VEP
η, G = 10, 1
t_M = η / G
εII = 1.0
args = (;)
pl2 = DruckerPrager(; C=η, ϕ=0) # plasticity
c_pl = CompositeRheology(LinearViscous(; η=η * Pa * s), ConstantElasticity(; G=G * Pa), pl2) # linear VEP
c_pl2 = CompositeRheology(ConstantElasticity(; G=G * Pa), pl2) # linear VEP

# case where old stress is below yield & new stress is above
args = (τII_old=9.8001101017963, P=0.0, τII=9.8001101017963)
args = (τII_old=9.8001101017963, P=0.0, τII=9.8001101017963)
F_old = compute_yieldfunction(c_pl.elements[3], args)

#
τ1, = local_iterations_εII(c_pl, εII, args; verbose=false, max_iter=10)
τ2, = compute_τII(c_pl, εII, args; verbose=false)
@test τ1 == τ2

args = merge(args, (τII=τ1,))
args = merge(args, (τII=τ1,))
F_check = compute_yieldfunction(c_pl.elements[3], args)
@test abs(F_check) < 1e-12
end
Loading