From 227f7136277bac6ec053709d99650fa88406f5ce Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Tue, 1 Oct 2024 06:31:25 +0000 Subject: [PATCH 01/33] Bump crate-ci/typos from 1.24.3 to 1.25.0 Bumps [crate-ci/typos](https://github.com/crate-ci/typos) from 1.24.3 to 1.25.0. - [Release notes](https://github.com/crate-ci/typos/releases) - [Changelog](https://github.com/crate-ci/typos/blob/master/CHANGELOG.md) - [Commits](https://github.com/crate-ci/typos/compare/v1.24.3...v1.25.0) --- updated-dependencies: - dependency-name: crate-ci/typos dependency-type: direct:production update-type: version-update:semver-minor ... Signed-off-by: dependabot[bot] --- .github/workflows/SpellCheck.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml index fedf097e..50e06406 100644 --- a/.github/workflows/SpellCheck.yml +++ b/.github/workflows/SpellCheck.yml @@ -10,4 +10,4 @@ jobs: - name: Checkout Actions Repository uses: actions/checkout@v4 - name: Check spelling - uses: crate-ci/typos@v1.24.3 + uses: crate-ci/typos@v1.25.0 From a9c53d461a6a2f6d46464d9ad7771c760a663aca Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Wed, 9 Oct 2024 12:53:14 +0200 Subject: [PATCH 02/33] init commit permeability Co-authored by: Jacob Frasunkiewicz --- src/GeoParams.jl | 23 +++-- src/MaterialParameters.jl | 97 ++++++++++--------- src/Permeability/Permeability.jl | 156 +++++++++++++++++++++++++++++++ 3 files changed, 225 insertions(+), 51 deletions(-) create mode 100644 src/Permeability/Permeability.jl diff --git a/src/GeoParams.jl b/src/GeoParams.jl index dea42709..aadc4929 100644 --- a/src/GeoParams.jl +++ b/src/GeoParams.jl @@ -335,6 +335,17 @@ export compute_meltfraction, Vector_MeltingParam, SmoothMelting + +using .MaterialParameters.Permeability +export compute_permeability, + compute_permeability!, + param_info, + AbstractPermeability, + ConstantPermeability, + HazenPermeability, + PowerLawPermeability, + CarmanKozenyPermeability + include("Traits/rheology.jl") export RheologyTrait export islinear, LinearRheologyTrait, NonLinearRheologyTrait @@ -363,7 +374,7 @@ function creeplaw_list(m::Module) out = string.(names(m; all=true, imported=true)) filter!(x -> !startswith(x, "#"), out) return [getfield(m, Symbol(x)) for x in out if !isnothing(tryparse(Int, string(x[end]))) || endswith(x, "a") || endswith(x, "b")] -end +end diffusion_law_list() = creeplaw_list(Diffusion) dislocation_law_list() = creeplaw_list(Dislocation) @@ -371,10 +382,10 @@ grainboundarysliding_law_list() = creeplaw_list(GBS) nonlinearpeierls_law_list() = creeplaw_list(NonLinearPeierls) peierls_law_list() = creeplaw_list(Peierls) -export diffusion_law_list, - dislocation_law_list, - grainboundarysliding_law_list, - nonlinearpeierls_law_list, +export diffusion_law_list, + dislocation_law_list, + grainboundarysliding_law_list, + nonlinearpeierls_law_list, peierls_law_list # Define Table output functions @@ -454,4 +465,4 @@ export get_G, get_Kb const get_shearmodulus = get_G const get_bulkmodulus = get_Kb -end # module GeoParams \ No newline at end of file +end # module GeoParams diff --git a/src/MaterialParameters.jl b/src/MaterialParameters.jl index 6fada711..2af74f4b 100644 --- a/src/MaterialParameters.jl +++ b/src/MaterialParameters.jl @@ -11,7 +11,7 @@ using Static import Base.show, Base.convert using GeoParams: - AbstractMaterialParam, AbstractMaterialParamsStruct, AbstractPhaseDiagramsStruct, AbstractComposite, ptr2string + AbstractMaterialParam, AbstractMaterialParamsStruct, AbstractPhaseDiagramsStruct, AbstractComposite, ptr2string # Define an "empty" Material parameter structure struct No_MaterialParam{_T} <: AbstractMaterialParam end @@ -44,6 +44,7 @@ include("./PhaseDiagrams/PhaseDiagrams.jl") #include("./Elasticity/Elasticity.jl") include("./ConstitutiveRelationships.jl") include("./Density/Density.jl") +include("./Permeability/Permeability.jl") include("./GravitationalAcceleration/GravitationalAcceleration.jl") include("./Energy/HeatCapacity.jl") include("./Energy/Conductivity.jl") @@ -58,7 +59,7 @@ using .ConstitutiveRelationships: print_rheology_matrix """ MaterialParams - + Structure that holds all material parameters for a given phase """ @@ -74,6 +75,7 @@ Structure that holds all material parameters for a given phase Vradioact<:Tuple, Vlatent<:Tuple, Vshearheat<:Tuple, + Vpermeability<:Tuple, Vmelting<:Tuple, Vseismvel<:Tuple, } <: AbstractMaterialParamsStruct @@ -85,38 +87,40 @@ Structure that holds all material parameters for a given phase CreepLaws::Vcreep = () # Creep laws Elasticity::Velastic = () # Elastic parameters Plasticity::Vplastic = () # Plasticity - CompositeRheology::Vcomposite = () # Composite (combined) rheologies - Conductivity::Vcond = () # Parameters related to the energy equation - HeatCapacity::Vheatc = () # Heat capacity + CompositeRheology::Vcomposite = () # Composite (combined) rheologies + Conductivity::Vcond = () # Parameters related to the energy equation + HeatCapacity::Vheatc = () # Heat capacity RadioactiveHeat::Vradioact = () # Radioactive heating source terms in energy conservation equation LatentHeat::Vlatent = () # Latent heating source terms in energy conservation equation ShearHeat::Vshearheat = () # Shear heating source terms in energy conservation equation + Permeability::Vpermeability = () # Permeability Melting::Vmelting = () # Melting model SeismicVelocity::Vseismvel = () # Seismic velocity end """ SetMaterialParams(; Name::String="", Phase::Int64=1, - Density = nothing, + Density = nothing, Gravity = nothing, - CreepLaws = nothing, - Elasticity = nothing, - Plasticity = nothing, + CreepLaws = nothing, + Elasticity = nothing, + Plasticity = nothing, CompositeRheology = nothing, - Conductivity = nothing, - HeatCapacity = nothing, + Conductivity = nothing, + HeatCapacity = nothing, RadioactiveHeat = nothing, LatentHeat = nothing, ShearHeat = nothing, + Permeability = nothing, Melting = nothing, SeismicVelocity = nothing, CharDim::GeoUnits = nothing) -Sets material parameters for a given phase. +Sets material parameters for a given phase. -If `CharDim` is specified the input parameters are non-dimensionalized. +If `CharDim` is specified the input parameters are non-dimensionalized. Note that if `Density` is specified, we also set `Gravity` even if not explicitly listed - + # Examples Define two viscous creep laws & constant density: @@ -126,10 +130,10 @@ julia> Phase = SetMaterialParams(Name="Viscous Matrix", CreepLaws = (PowerlawViscous(), LinearViscous(η=1e21Pa*s))) Phase 1 : Viscous Matrix | [dimensional units] - | - |-- Density : Constant density: ρ=2900 kg m⁻³ - |-- Gravity : Gravitational acceleration: g=9.81 m s⁻² - |-- CreepLaws : Powerlaw viscosity: η0=1.0e18 Pa s, n=2.0, ε0=1.0e-15 s⁻¹ + | + |-- Density : Constant density: ρ=2900 kg m⁻³ + |-- Gravity : Gravitational acceleration: g=9.81 m s⁻² + |-- CreepLaws : Powerlaw viscosity: η0=1.0e18 Pa s, n=2.0, ε0=1.0e-15 s⁻¹ | Linear viscosity: η=1.0e21 Pa s ``` @@ -142,11 +146,11 @@ julia> Phase = SetMaterialParams(Name="Viscous Matrix", Phase=33, CharDim = CharUnits_GEO) Phase 33: Viscous Matrix | [non-dimensional units] - | - |-- Density : P/T-dependent density: ρ0=2.9e-16, α=0.038194500000000006, β=0.01, T0=0.21454659702313156, P0=0.0 - |-- Gravity : Gravitational acceleration: g=9.810000000000002e18 - |-- CreepLaws : Powerlaw viscosity: η0=0.1, n=3, ε0=0.001 - | Linear viscosity: η=10000.0 + | + |-- Density : P/T-dependent density: ρ0=2.9e-16, α=0.038194500000000006, β=0.01, T0=0.21454659702313156, P0=0.0 + |-- Gravity : Gravitational acceleration: g=9.810000000000002e18 + |-- CreepLaws : Powerlaw viscosity: η0=0.1, n=3, ε0=0.001 + | Linear viscosity: η=10000.0 ``` You can also create an array that holds several parameters: @@ -164,19 +168,19 @@ julia> MatParam 2-element Vector{MaterialParams}: Phase 1 : Upper Crust | [dimensional units] - | - |-- Density : Constant density: ρ=2900 kg m⁻³ - |-- Gravity : Gravitational acceleration: g=9.81 m s⁻² - |-- CreepLaws : Powerlaw viscosity: η0=1.0e18 Pa s, n=2.0, ε0=1.0e-15 s⁻¹ - | Linear viscosity: η=1.0e23 Pa s - + | + |-- Density : Constant density: ρ=2900 kg m⁻³ + |-- Gravity : Gravitational acceleration: g=9.81 m s⁻² + |-- CreepLaws : Powerlaw viscosity: η0=1.0e18 Pa s, n=2.0, ε0=1.0e-15 s⁻¹ + | Linear viscosity: η=1.0e23 Pa s + Phase 2 : Lower Crust | [dimensional units] - | - |-- Density : P/T-dependent density: ρ0=3000 kg m⁻³, α=3.0e-5 K⁻¹, β=1.0e-9 Pa⁻¹, T0=0 °C, P0=0 MPa - |-- Gravity : Gravitational acceleration: g=9.81 m s⁻² - |-- CreepLaws : Powerlaw viscosity: η0=1.0e18 Pa s, n=5, ε0=1.0e-15 s⁻¹ - | Linear viscosity: η=1.0e21 Pa s + | + |-- Density : P/T-dependent density: ρ0=3000 kg m⁻³, α=3.0e-5 K⁻¹, β=1.0e-9 Pa⁻¹, T0=0 °C, P0=0 MPa + |-- Gravity : Gravitational acceleration: g=9.81 m s⁻² + |-- CreepLaws : Powerlaw viscosity: η0=1.0e18 Pa s, n=5, ε0=1.0e-15 s⁻¹ + | Linear viscosity: η=1.0e21 Pa s ``` @@ -195,6 +199,7 @@ function SetMaterialParams(; RadioactiveHeat=nothing, LatentHeat=nothing, ShearHeat=nothing, + Permeability=nothing, Melting=nothing, SeismicVelocity=nothing, CharDim=nothing, @@ -213,6 +218,7 @@ function SetMaterialParams(; ConvField(RadioactiveHeat, :RadioactiveHeat; maxAllowedFields=1), ConvField(LatentHeat, :LatentHeat; maxAllowedFields=1), ConvField(ShearHeat, :ShearHeat; maxAllowedFields=1), + ConvField(Permeability, :Permeability; maxAllowedFields=1), ConvField(Melting, :Melting; maxAllowedFields=1), ConvField(SeismicVelocity, :SeismicVelocity; maxAllowedFields=1), CharDim, @@ -233,12 +239,13 @@ function SetMaterialParams( RadioactiveHeat, LatentHeat, ShearHeat, + Permeability, Melting, SeismicVelocity, CharDim, ) - # define struct for phase, while also specifying the maximum number of definitions for every field + # define struct for phase, while also specifying the maximum number of definitions for every field phase = MaterialParams( pointer(ptr2string(Name)), Phase, @@ -254,12 +261,13 @@ function SetMaterialParams( RadioactiveHeat, LatentHeat, ShearHeat, + Permeability, Melting, SeismicVelocity, ) # [optionally] non-dimensionalize the struct - phase_nd = nondimensionalize_phase(phase, CharDim) + phase_nd = nondimensionalize_phase(phase, CharDim) return phase_nd end @@ -274,11 +282,11 @@ end set_gravity(Gravity, Density) = Gravity # Helper function that converts a field to a Tuple, provided it is not nothing -# This also checks for the maximum allowed number of definitions -# (some rheological phases may allow for an arbitrary combination per phase; others like density EoS not) +# This also checks for the maximum allowed number of definitions +# (some rheological phases may allow for an arbitrary combination per phase; others like density EoS not) ConvField(::Nothing, fieldname::Symbol; maxAllowedFields=1e6) = () ConvField(field::AbstractMaterialParam, fieldname::Symbol; maxAllowedFields=1e6) = (field, ) -function ConvField(field::NTuple{N, Any}, fieldname::Symbol; maxAllowedFields=1e6) where N +function ConvField(field::NTuple{N, Any}, fieldname::Symbol; maxAllowedFields=1e6) where N if length(field) > maxAllowedFields error("Maximum $(maxAllowedFields) field allowed for: $fieldname") end @@ -286,11 +294,11 @@ function ConvField(field::NTuple{N, Any}, fieldname::Symbol; maxAllowedFields=1e end # Helper that prints info about each of the material parameters -# for this to look nice, you need to define a Base.show +# for this to look nice, you need to define a Base.show function Print_MaterialParam(io::IO, name::Symbol, Data) if length(Data) > 0 if Data isa Ptr - str = unsafe_string(Data) + str = unsafe_string(Data) print(io, " |-- $(rpad(name,18)): $str \n") elseif typeof(Data[1]) <: AbstractMaterialParam @@ -299,7 +307,7 @@ function Print_MaterialParam(io::IO, name::Symbol, Data) str = Data[i] if isa(str, AbstractComposite) # The CompositeRheology object is formatted a bit different - str = print_composite(Data[i],32) + str = print_composite(Data[i],32) end if i == 1 @@ -329,7 +337,7 @@ function Base.show(io::IO, phase::MaterialParams) end end -# Slightly nicer printout in case we have a tuple with material parameters +# Slightly nicer printout in case we have a tuple with material parameters function Base.show(io::IO, phase_tuple::NTuple{N,MaterialParams}) where {N} for i in 1:N Base.show(io, phase_tuple[i]) @@ -349,7 +357,7 @@ function print_composite(a, spaces=10) str = str.*"\n" for i=2:length(str) for j=1:spaces - str[i] = " "*str[i] + str[i] = " "*str[i] end end str = join(str) @@ -358,4 +366,3 @@ function print_composite(a, spaces=10) end end - diff --git a/src/Permeability/Permeability.jl b/src/Permeability/Permeability.jl new file mode 100644 index 00000000..86f03417 --- /dev/null +++ b/src/Permeability/Permeability.jl @@ -0,0 +1,156 @@ +module Permeability + +using Parameters, Unitful, LaTeXStrings, MuladdMacro +using ..Units +using ..PhaseDiagrams +using GeoParams: AbstractMaterialParam, AbstractMaterialParamsStruct, @extractors, add_extractor_functions +import ..Units: isdimensional +using ..MaterialParameters: No_MaterialParam, MaterialParamsInfo +import Base.show, GeoParams.param_info + +include("../Computations.jl") + +abstract type AbstractPermeability{T} <: AbstractMaterialParam end + +export compute_permeability, # calculation routines + compute_permeability!, # in place calculation + param_info, # info about the parameters + AbstractPermeability, + ConstantPermeability, # constant + HazenPermeability, # Hazen equation + PowerLawPermeability, # Power-law permeability + CarmanKozenyPermeability # Carman-Kozeny permeability + +# Define "empty" computational routines in case nothing is defined +function compute_permeability!( + k::_T, s::No_MaterialParam{_T}; ϕ::_T=zero(_T) +) where {_T} + return zero(_T) +end +function compute_permeability(s::No_MaterialParam{_T}; ϕ::_T=zero(_T)) where {_T} + return zero(_T) +end + +# Constant Permeability +@with_kw_noshow struct ConstantPermeability{_T,U} <: AbstractPermeability{_T} + k::GeoUnit{_T,U} = 1e-12m^2 # permeability +end +ConstantPermeability(args...) = ConstantPermeability(convert.(GeoUnit, args)...) +isdimensional(s::ConstantPermeability) = isdimensional(s.k) + +@inline compute_permeability(s::ConstantPermeability{_T}, args) where {_T} = s(; args...) +@inline compute_permeability(s::ConstantPermeability{_T}) where {_T} = s() + +function param_info(s::ConstantPermeability) + return MaterialParamsInfo(; Equation=L"k = cst") +end + +function compute_permeability!(k::AbstractArray, s::ConstantPermeability; kwargs...) + @unpack_val k_val = s + k[:] .= k_val + return nothing +end + +function compute_permeability!(k::AbstractArray, s::ConstantPermeability, args) + return compute_permeability!(k, s; args...) +end + +function show(io::IO, g::ConstantPermeability) + return print(io, "Constant permeability: k=$(UnitValue(g.k))") +end + +# Hazen Permeability +@with_kw_noshow struct HazenPermeability{_T,U1,U2} <: AbstractPermeability{_T} + C::GeoUnit{_T,U1} = 1.0 # Hazen constant + D10::GeoUnit{_T,U2} = 1e-4m # Effective grain size +end +HazenPermeability(args...) = HazenPermeability(convert.(GeoUnit, args)...) +isdimensional(s::HazenPermeability) = isdimensional(s.C) + +function param_info(s::HazenPermeability) + return MaterialParamsInfo(; Equation = L"k = C \cdot D_{10}^2") +end + +function (s::HazenPermeability{_T})(; kwargs...) where {_T} + if kwargs isa Quantity + @unpack_units C, D10 = s + else + @unpack_val C, D10 = s + end + + return C * D10^2 +end + +@inline (s::HazenPermeability)(args) = s(; args...) +@inline compute_permeability(s::HazenPermeability, args) = s(args) + +function show(io::IO, g::HazenPermeability) + return print(io, "Hazen permeability: k = C * D10^2; C=$(g.C); D10=$(g.D10)") +end + +# Power-law Permeability +@with_kw_noshow struct PowerLawPermeability{_T,U1,U2,U3} <: AbstractPermeability{_T} + c::GeoUnit{_T,U1} = 1.0 # Power-law constant + k0::GeoUnit{_T,U1} = 1e-12m^2 # reference permeability + n::GeoUnit{_T,U3} = 3.0 # exponent +end +PowerLawPermeability(args...) = PowerLawPermeability(convert.(GeoUnit, args)...) +isdimensional(s::PowerLawPermeability) = isdimensional(s.k0) + +function param_info(s::PowerLawPermeability) + return MaterialParamsInfo(; Equation = L"k = c * k_0 * (\phi^n)") +end + +function (s::PowerLawPermeability{_T})(; ϕ=0e0, kwargs...) where {_T} + if ϕ isa Quantity + @unpack_units c, k0, n = s + else + @unpack_val c, k0, n = s + end + + return c * k0 * ϕ^n +end + +@inline (s::PowerLawPermeability)(args) = s(; args...) +@inline compute_permeability(s::PowerLawPermeability, args) = s(args) + +function show(io::IO, g::PowerLawPermeability) + return print(io, "Power-law permeability: k = c* k0 * ϕ^n; c = $(g.c), k0=$(g.k0); n=$(g.n)") +end + +# Carman-Kozeny Permeability +@with_kw_noshow struct CarmanKozenyPermeability{_T,U1,U2,U3} <: AbstractPermeability{_T} + c::GeoUnit{_T,U1} = 1.0 # Carman-Kozeny constant + ϕ0::GeoUnit{_T,U2} = 0.01 # reference porosity + n::GeoUnit{_T,U3} = 3.0 # exponent +end +CarmanKozenyPermeability(args...) = CarmanKozenyPermeability(convert.(GeoUnit, args)...) +isdimensional(s::CarmanKozenyPermeability) = isdimensional(s.c) + +function param_info(s::CarmanKozenyPermeability) + return MaterialParamsInfo(; Equation = L"k = c \left(\frac{\phi}{\phi_0}\right)^n") +end + +function (s::CarmanKozenyPermeability{_T})(; ϕ=1e-2, kwargs...) where {_T} + if ϕ isa Quantity + @unpack_units c, ϕ0, n = s + else + @unpack_val c, ϕ0, n = s + end + + return c * (ϕ / ϕ0)^n +end + +@inline (s::CarmanKozenyPermeability)(args) = s(; args...) +@inline compute_permeability(s::CarmanKozenyPermeability, args) = s(args) + +function show(io::IO, g::CarmanKozenyPermeability) + return print(io, "Carman-Kozeny permeability: k = c * (ϕ / ϕ0)^n; c=$(g.c); ϕ0=$(g.ϕ0); n=$(g.n)") +end + +# extractor methods +for type in (ConstantPermeability, HazenPermeability, PowerLawPermeability, CarmanKozenyPermeability) + @extractors(type, :Permeability) +end + +end From 901742242a224a8e48f33d22d93acfb458ecd4dc Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Wed, 9 Oct 2024 12:53:14 +0200 Subject: [PATCH 03/33] init commit permeability Co-authored by: Jacob Frasunkiewicz --- src/GeoParams.jl | 23 +++-- src/MaterialParameters.jl | 97 ++++++++++--------- src/Permeability/Permeability.jl | 156 +++++++++++++++++++++++++++++++ 3 files changed, 225 insertions(+), 51 deletions(-) create mode 100644 src/Permeability/Permeability.jl diff --git a/src/GeoParams.jl b/src/GeoParams.jl index dea42709..aadc4929 100644 --- a/src/GeoParams.jl +++ b/src/GeoParams.jl @@ -335,6 +335,17 @@ export compute_meltfraction, Vector_MeltingParam, SmoothMelting + +using .MaterialParameters.Permeability +export compute_permeability, + compute_permeability!, + param_info, + AbstractPermeability, + ConstantPermeability, + HazenPermeability, + PowerLawPermeability, + CarmanKozenyPermeability + include("Traits/rheology.jl") export RheologyTrait export islinear, LinearRheologyTrait, NonLinearRheologyTrait @@ -363,7 +374,7 @@ function creeplaw_list(m::Module) out = string.(names(m; all=true, imported=true)) filter!(x -> !startswith(x, "#"), out) return [getfield(m, Symbol(x)) for x in out if !isnothing(tryparse(Int, string(x[end]))) || endswith(x, "a") || endswith(x, "b")] -end +end diffusion_law_list() = creeplaw_list(Diffusion) dislocation_law_list() = creeplaw_list(Dislocation) @@ -371,10 +382,10 @@ grainboundarysliding_law_list() = creeplaw_list(GBS) nonlinearpeierls_law_list() = creeplaw_list(NonLinearPeierls) peierls_law_list() = creeplaw_list(Peierls) -export diffusion_law_list, - dislocation_law_list, - grainboundarysliding_law_list, - nonlinearpeierls_law_list, +export diffusion_law_list, + dislocation_law_list, + grainboundarysliding_law_list, + nonlinearpeierls_law_list, peierls_law_list # Define Table output functions @@ -454,4 +465,4 @@ export get_G, get_Kb const get_shearmodulus = get_G const get_bulkmodulus = get_Kb -end # module GeoParams \ No newline at end of file +end # module GeoParams diff --git a/src/MaterialParameters.jl b/src/MaterialParameters.jl index 6fada711..2af74f4b 100644 --- a/src/MaterialParameters.jl +++ b/src/MaterialParameters.jl @@ -11,7 +11,7 @@ using Static import Base.show, Base.convert using GeoParams: - AbstractMaterialParam, AbstractMaterialParamsStruct, AbstractPhaseDiagramsStruct, AbstractComposite, ptr2string + AbstractMaterialParam, AbstractMaterialParamsStruct, AbstractPhaseDiagramsStruct, AbstractComposite, ptr2string # Define an "empty" Material parameter structure struct No_MaterialParam{_T} <: AbstractMaterialParam end @@ -44,6 +44,7 @@ include("./PhaseDiagrams/PhaseDiagrams.jl") #include("./Elasticity/Elasticity.jl") include("./ConstitutiveRelationships.jl") include("./Density/Density.jl") +include("./Permeability/Permeability.jl") include("./GravitationalAcceleration/GravitationalAcceleration.jl") include("./Energy/HeatCapacity.jl") include("./Energy/Conductivity.jl") @@ -58,7 +59,7 @@ using .ConstitutiveRelationships: print_rheology_matrix """ MaterialParams - + Structure that holds all material parameters for a given phase """ @@ -74,6 +75,7 @@ Structure that holds all material parameters for a given phase Vradioact<:Tuple, Vlatent<:Tuple, Vshearheat<:Tuple, + Vpermeability<:Tuple, Vmelting<:Tuple, Vseismvel<:Tuple, } <: AbstractMaterialParamsStruct @@ -85,38 +87,40 @@ Structure that holds all material parameters for a given phase CreepLaws::Vcreep = () # Creep laws Elasticity::Velastic = () # Elastic parameters Plasticity::Vplastic = () # Plasticity - CompositeRheology::Vcomposite = () # Composite (combined) rheologies - Conductivity::Vcond = () # Parameters related to the energy equation - HeatCapacity::Vheatc = () # Heat capacity + CompositeRheology::Vcomposite = () # Composite (combined) rheologies + Conductivity::Vcond = () # Parameters related to the energy equation + HeatCapacity::Vheatc = () # Heat capacity RadioactiveHeat::Vradioact = () # Radioactive heating source terms in energy conservation equation LatentHeat::Vlatent = () # Latent heating source terms in energy conservation equation ShearHeat::Vshearheat = () # Shear heating source terms in energy conservation equation + Permeability::Vpermeability = () # Permeability Melting::Vmelting = () # Melting model SeismicVelocity::Vseismvel = () # Seismic velocity end """ SetMaterialParams(; Name::String="", Phase::Int64=1, - Density = nothing, + Density = nothing, Gravity = nothing, - CreepLaws = nothing, - Elasticity = nothing, - Plasticity = nothing, + CreepLaws = nothing, + Elasticity = nothing, + Plasticity = nothing, CompositeRheology = nothing, - Conductivity = nothing, - HeatCapacity = nothing, + Conductivity = nothing, + HeatCapacity = nothing, RadioactiveHeat = nothing, LatentHeat = nothing, ShearHeat = nothing, + Permeability = nothing, Melting = nothing, SeismicVelocity = nothing, CharDim::GeoUnits = nothing) -Sets material parameters for a given phase. +Sets material parameters for a given phase. -If `CharDim` is specified the input parameters are non-dimensionalized. +If `CharDim` is specified the input parameters are non-dimensionalized. Note that if `Density` is specified, we also set `Gravity` even if not explicitly listed - + # Examples Define two viscous creep laws & constant density: @@ -126,10 +130,10 @@ julia> Phase = SetMaterialParams(Name="Viscous Matrix", CreepLaws = (PowerlawViscous(), LinearViscous(η=1e21Pa*s))) Phase 1 : Viscous Matrix | [dimensional units] - | - |-- Density : Constant density: ρ=2900 kg m⁻³ - |-- Gravity : Gravitational acceleration: g=9.81 m s⁻² - |-- CreepLaws : Powerlaw viscosity: η0=1.0e18 Pa s, n=2.0, ε0=1.0e-15 s⁻¹ + | + |-- Density : Constant density: ρ=2900 kg m⁻³ + |-- Gravity : Gravitational acceleration: g=9.81 m s⁻² + |-- CreepLaws : Powerlaw viscosity: η0=1.0e18 Pa s, n=2.0, ε0=1.0e-15 s⁻¹ | Linear viscosity: η=1.0e21 Pa s ``` @@ -142,11 +146,11 @@ julia> Phase = SetMaterialParams(Name="Viscous Matrix", Phase=33, CharDim = CharUnits_GEO) Phase 33: Viscous Matrix | [non-dimensional units] - | - |-- Density : P/T-dependent density: ρ0=2.9e-16, α=0.038194500000000006, β=0.01, T0=0.21454659702313156, P0=0.0 - |-- Gravity : Gravitational acceleration: g=9.810000000000002e18 - |-- CreepLaws : Powerlaw viscosity: η0=0.1, n=3, ε0=0.001 - | Linear viscosity: η=10000.0 + | + |-- Density : P/T-dependent density: ρ0=2.9e-16, α=0.038194500000000006, β=0.01, T0=0.21454659702313156, P0=0.0 + |-- Gravity : Gravitational acceleration: g=9.810000000000002e18 + |-- CreepLaws : Powerlaw viscosity: η0=0.1, n=3, ε0=0.001 + | Linear viscosity: η=10000.0 ``` You can also create an array that holds several parameters: @@ -164,19 +168,19 @@ julia> MatParam 2-element Vector{MaterialParams}: Phase 1 : Upper Crust | [dimensional units] - | - |-- Density : Constant density: ρ=2900 kg m⁻³ - |-- Gravity : Gravitational acceleration: g=9.81 m s⁻² - |-- CreepLaws : Powerlaw viscosity: η0=1.0e18 Pa s, n=2.0, ε0=1.0e-15 s⁻¹ - | Linear viscosity: η=1.0e23 Pa s - + | + |-- Density : Constant density: ρ=2900 kg m⁻³ + |-- Gravity : Gravitational acceleration: g=9.81 m s⁻² + |-- CreepLaws : Powerlaw viscosity: η0=1.0e18 Pa s, n=2.0, ε0=1.0e-15 s⁻¹ + | Linear viscosity: η=1.0e23 Pa s + Phase 2 : Lower Crust | [dimensional units] - | - |-- Density : P/T-dependent density: ρ0=3000 kg m⁻³, α=3.0e-5 K⁻¹, β=1.0e-9 Pa⁻¹, T0=0 °C, P0=0 MPa - |-- Gravity : Gravitational acceleration: g=9.81 m s⁻² - |-- CreepLaws : Powerlaw viscosity: η0=1.0e18 Pa s, n=5, ε0=1.0e-15 s⁻¹ - | Linear viscosity: η=1.0e21 Pa s + | + |-- Density : P/T-dependent density: ρ0=3000 kg m⁻³, α=3.0e-5 K⁻¹, β=1.0e-9 Pa⁻¹, T0=0 °C, P0=0 MPa + |-- Gravity : Gravitational acceleration: g=9.81 m s⁻² + |-- CreepLaws : Powerlaw viscosity: η0=1.0e18 Pa s, n=5, ε0=1.0e-15 s⁻¹ + | Linear viscosity: η=1.0e21 Pa s ``` @@ -195,6 +199,7 @@ function SetMaterialParams(; RadioactiveHeat=nothing, LatentHeat=nothing, ShearHeat=nothing, + Permeability=nothing, Melting=nothing, SeismicVelocity=nothing, CharDim=nothing, @@ -213,6 +218,7 @@ function SetMaterialParams(; ConvField(RadioactiveHeat, :RadioactiveHeat; maxAllowedFields=1), ConvField(LatentHeat, :LatentHeat; maxAllowedFields=1), ConvField(ShearHeat, :ShearHeat; maxAllowedFields=1), + ConvField(Permeability, :Permeability; maxAllowedFields=1), ConvField(Melting, :Melting; maxAllowedFields=1), ConvField(SeismicVelocity, :SeismicVelocity; maxAllowedFields=1), CharDim, @@ -233,12 +239,13 @@ function SetMaterialParams( RadioactiveHeat, LatentHeat, ShearHeat, + Permeability, Melting, SeismicVelocity, CharDim, ) - # define struct for phase, while also specifying the maximum number of definitions for every field + # define struct for phase, while also specifying the maximum number of definitions for every field phase = MaterialParams( pointer(ptr2string(Name)), Phase, @@ -254,12 +261,13 @@ function SetMaterialParams( RadioactiveHeat, LatentHeat, ShearHeat, + Permeability, Melting, SeismicVelocity, ) # [optionally] non-dimensionalize the struct - phase_nd = nondimensionalize_phase(phase, CharDim) + phase_nd = nondimensionalize_phase(phase, CharDim) return phase_nd end @@ -274,11 +282,11 @@ end set_gravity(Gravity, Density) = Gravity # Helper function that converts a field to a Tuple, provided it is not nothing -# This also checks for the maximum allowed number of definitions -# (some rheological phases may allow for an arbitrary combination per phase; others like density EoS not) +# This also checks for the maximum allowed number of definitions +# (some rheological phases may allow for an arbitrary combination per phase; others like density EoS not) ConvField(::Nothing, fieldname::Symbol; maxAllowedFields=1e6) = () ConvField(field::AbstractMaterialParam, fieldname::Symbol; maxAllowedFields=1e6) = (field, ) -function ConvField(field::NTuple{N, Any}, fieldname::Symbol; maxAllowedFields=1e6) where N +function ConvField(field::NTuple{N, Any}, fieldname::Symbol; maxAllowedFields=1e6) where N if length(field) > maxAllowedFields error("Maximum $(maxAllowedFields) field allowed for: $fieldname") end @@ -286,11 +294,11 @@ function ConvField(field::NTuple{N, Any}, fieldname::Symbol; maxAllowedFields=1e end # Helper that prints info about each of the material parameters -# for this to look nice, you need to define a Base.show +# for this to look nice, you need to define a Base.show function Print_MaterialParam(io::IO, name::Symbol, Data) if length(Data) > 0 if Data isa Ptr - str = unsafe_string(Data) + str = unsafe_string(Data) print(io, " |-- $(rpad(name,18)): $str \n") elseif typeof(Data[1]) <: AbstractMaterialParam @@ -299,7 +307,7 @@ function Print_MaterialParam(io::IO, name::Symbol, Data) str = Data[i] if isa(str, AbstractComposite) # The CompositeRheology object is formatted a bit different - str = print_composite(Data[i],32) + str = print_composite(Data[i],32) end if i == 1 @@ -329,7 +337,7 @@ function Base.show(io::IO, phase::MaterialParams) end end -# Slightly nicer printout in case we have a tuple with material parameters +# Slightly nicer printout in case we have a tuple with material parameters function Base.show(io::IO, phase_tuple::NTuple{N,MaterialParams}) where {N} for i in 1:N Base.show(io, phase_tuple[i]) @@ -349,7 +357,7 @@ function print_composite(a, spaces=10) str = str.*"\n" for i=2:length(str) for j=1:spaces - str[i] = " "*str[i] + str[i] = " "*str[i] end end str = join(str) @@ -358,4 +366,3 @@ function print_composite(a, spaces=10) end end - diff --git a/src/Permeability/Permeability.jl b/src/Permeability/Permeability.jl new file mode 100644 index 00000000..86f03417 --- /dev/null +++ b/src/Permeability/Permeability.jl @@ -0,0 +1,156 @@ +module Permeability + +using Parameters, Unitful, LaTeXStrings, MuladdMacro +using ..Units +using ..PhaseDiagrams +using GeoParams: AbstractMaterialParam, AbstractMaterialParamsStruct, @extractors, add_extractor_functions +import ..Units: isdimensional +using ..MaterialParameters: No_MaterialParam, MaterialParamsInfo +import Base.show, GeoParams.param_info + +include("../Computations.jl") + +abstract type AbstractPermeability{T} <: AbstractMaterialParam end + +export compute_permeability, # calculation routines + compute_permeability!, # in place calculation + param_info, # info about the parameters + AbstractPermeability, + ConstantPermeability, # constant + HazenPermeability, # Hazen equation + PowerLawPermeability, # Power-law permeability + CarmanKozenyPermeability # Carman-Kozeny permeability + +# Define "empty" computational routines in case nothing is defined +function compute_permeability!( + k::_T, s::No_MaterialParam{_T}; ϕ::_T=zero(_T) +) where {_T} + return zero(_T) +end +function compute_permeability(s::No_MaterialParam{_T}; ϕ::_T=zero(_T)) where {_T} + return zero(_T) +end + +# Constant Permeability +@with_kw_noshow struct ConstantPermeability{_T,U} <: AbstractPermeability{_T} + k::GeoUnit{_T,U} = 1e-12m^2 # permeability +end +ConstantPermeability(args...) = ConstantPermeability(convert.(GeoUnit, args)...) +isdimensional(s::ConstantPermeability) = isdimensional(s.k) + +@inline compute_permeability(s::ConstantPermeability{_T}, args) where {_T} = s(; args...) +@inline compute_permeability(s::ConstantPermeability{_T}) where {_T} = s() + +function param_info(s::ConstantPermeability) + return MaterialParamsInfo(; Equation=L"k = cst") +end + +function compute_permeability!(k::AbstractArray, s::ConstantPermeability; kwargs...) + @unpack_val k_val = s + k[:] .= k_val + return nothing +end + +function compute_permeability!(k::AbstractArray, s::ConstantPermeability, args) + return compute_permeability!(k, s; args...) +end + +function show(io::IO, g::ConstantPermeability) + return print(io, "Constant permeability: k=$(UnitValue(g.k))") +end + +# Hazen Permeability +@with_kw_noshow struct HazenPermeability{_T,U1,U2} <: AbstractPermeability{_T} + C::GeoUnit{_T,U1} = 1.0 # Hazen constant + D10::GeoUnit{_T,U2} = 1e-4m # Effective grain size +end +HazenPermeability(args...) = HazenPermeability(convert.(GeoUnit, args)...) +isdimensional(s::HazenPermeability) = isdimensional(s.C) + +function param_info(s::HazenPermeability) + return MaterialParamsInfo(; Equation = L"k = C \cdot D_{10}^2") +end + +function (s::HazenPermeability{_T})(; kwargs...) where {_T} + if kwargs isa Quantity + @unpack_units C, D10 = s + else + @unpack_val C, D10 = s + end + + return C * D10^2 +end + +@inline (s::HazenPermeability)(args) = s(; args...) +@inline compute_permeability(s::HazenPermeability, args) = s(args) + +function show(io::IO, g::HazenPermeability) + return print(io, "Hazen permeability: k = C * D10^2; C=$(g.C); D10=$(g.D10)") +end + +# Power-law Permeability +@with_kw_noshow struct PowerLawPermeability{_T,U1,U2,U3} <: AbstractPermeability{_T} + c::GeoUnit{_T,U1} = 1.0 # Power-law constant + k0::GeoUnit{_T,U1} = 1e-12m^2 # reference permeability + n::GeoUnit{_T,U3} = 3.0 # exponent +end +PowerLawPermeability(args...) = PowerLawPermeability(convert.(GeoUnit, args)...) +isdimensional(s::PowerLawPermeability) = isdimensional(s.k0) + +function param_info(s::PowerLawPermeability) + return MaterialParamsInfo(; Equation = L"k = c * k_0 * (\phi^n)") +end + +function (s::PowerLawPermeability{_T})(; ϕ=0e0, kwargs...) where {_T} + if ϕ isa Quantity + @unpack_units c, k0, n = s + else + @unpack_val c, k0, n = s + end + + return c * k0 * ϕ^n +end + +@inline (s::PowerLawPermeability)(args) = s(; args...) +@inline compute_permeability(s::PowerLawPermeability, args) = s(args) + +function show(io::IO, g::PowerLawPermeability) + return print(io, "Power-law permeability: k = c* k0 * ϕ^n; c = $(g.c), k0=$(g.k0); n=$(g.n)") +end + +# Carman-Kozeny Permeability +@with_kw_noshow struct CarmanKozenyPermeability{_T,U1,U2,U3} <: AbstractPermeability{_T} + c::GeoUnit{_T,U1} = 1.0 # Carman-Kozeny constant + ϕ0::GeoUnit{_T,U2} = 0.01 # reference porosity + n::GeoUnit{_T,U3} = 3.0 # exponent +end +CarmanKozenyPermeability(args...) = CarmanKozenyPermeability(convert.(GeoUnit, args)...) +isdimensional(s::CarmanKozenyPermeability) = isdimensional(s.c) + +function param_info(s::CarmanKozenyPermeability) + return MaterialParamsInfo(; Equation = L"k = c \left(\frac{\phi}{\phi_0}\right)^n") +end + +function (s::CarmanKozenyPermeability{_T})(; ϕ=1e-2, kwargs...) where {_T} + if ϕ isa Quantity + @unpack_units c, ϕ0, n = s + else + @unpack_val c, ϕ0, n = s + end + + return c * (ϕ / ϕ0)^n +end + +@inline (s::CarmanKozenyPermeability)(args) = s(; args...) +@inline compute_permeability(s::CarmanKozenyPermeability, args) = s(args) + +function show(io::IO, g::CarmanKozenyPermeability) + return print(io, "Carman-Kozeny permeability: k = c * (ϕ / ϕ0)^n; c=$(g.c); ϕ0=$(g.ϕ0); n=$(g.n)") +end + +# extractor methods +for type in (ConstantPermeability, HazenPermeability, PowerLawPermeability, CarmanKozenyPermeability) + @extractors(type, :Permeability) +end + +end From 7a58461c0ea409ff2fe82d7fd73c7ec0a60ebf07 Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Wed, 9 Oct 2024 17:39:14 +0200 Subject: [PATCH 04/33] add a few tests --- src/Permeability/Permeability.jl | 28 +++++------ src/aliases.jl | 1 + test/test_Permeability.jl | 79 ++++++++++++++++++++++++++++++++ 3 files changed, 95 insertions(+), 13 deletions(-) create mode 100644 test/test_Permeability.jl diff --git a/src/Permeability/Permeability.jl b/src/Permeability/Permeability.jl index 86f03417..1fbfa662 100644 --- a/src/Permeability/Permeability.jl +++ b/src/Permeability/Permeability.jl @@ -61,11 +61,11 @@ end # Hazen Permeability @with_kw_noshow struct HazenPermeability{_T,U1,U2} <: AbstractPermeability{_T} - C::GeoUnit{_T,U1} = 1.0 # Hazen constant - D10::GeoUnit{_T,U2} = 1e-4m # Effective grain size + C::GeoUnit{_T,U1} = 1.0 * NoUnits # Hazen constant + D10::GeoUnit{_T,U2} = 1e-4 * m # Effective grain size end HazenPermeability(args...) = HazenPermeability(convert.(GeoUnit, args)...) -isdimensional(s::HazenPermeability) = isdimensional(s.C) +isdimensional(s::HazenPermeability) = isdimensional(s.D10) function param_info(s::HazenPermeability) return MaterialParamsInfo(; Equation = L"k = C \cdot D_{10}^2") @@ -82,23 +82,25 @@ function (s::HazenPermeability{_T})(; kwargs...) where {_T} end @inline (s::HazenPermeability)(args) = s(; args...) -@inline compute_permeability(s::HazenPermeability, args) = s(args) +@inline compute_permeability(s::HazenPermeability, args) = s(; args...) + function show(io::IO, g::HazenPermeability) return print(io, "Hazen permeability: k = C * D10^2; C=$(g.C); D10=$(g.D10)") end # Power-law Permeability -@with_kw_noshow struct PowerLawPermeability{_T,U1,U2,U3} <: AbstractPermeability{_T} - c::GeoUnit{_T,U1} = 1.0 # Power-law constant - k0::GeoUnit{_T,U1} = 1e-12m^2 # reference permeability - n::GeoUnit{_T,U3} = 3.0 # exponent +@with_kw_noshow struct PowerLawPermeability{_T,U1,U2,U3, U4} <: AbstractPermeability{_T} + c::GeoUnit{_T,U1} = 1.0 * NoUnits # Power-law constant + k0::GeoUnit{_T,U2} = 1e-12 * m^2 # reference permeability + ϕ::GeoUnit{_T,U3} = 1e-2 * NoUnits # reference porosity + n::GeoUnit{_T,U4} = 3.0 * NoUnits # exponent end PowerLawPermeability(args...) = PowerLawPermeability(convert.(GeoUnit, args)...) isdimensional(s::PowerLawPermeability) = isdimensional(s.k0) function param_info(s::PowerLawPermeability) - return MaterialParamsInfo(; Equation = L"k = c * k_0 * (\phi^n)") + return MaterialParamsInfo(; Equation = L"k = c \cdot k_0 \cdot \phi^n") end function (s::PowerLawPermeability{_T})(; ϕ=0e0, kwargs...) where {_T} @@ -120,12 +122,12 @@ end # Carman-Kozeny Permeability @with_kw_noshow struct CarmanKozenyPermeability{_T,U1,U2,U3} <: AbstractPermeability{_T} - c::GeoUnit{_T,U1} = 1.0 # Carman-Kozeny constant - ϕ0::GeoUnit{_T,U2} = 0.01 # reference porosity - n::GeoUnit{_T,U3} = 3.0 # exponent + c::GeoUnit{_T,U1} = 1.0 * NoUnits # Carman-Kozeny constant + ϕ0::GeoUnit{_T,U2} = 0.01 * NoUnits # reference porosity + n::GeoUnit{_T,U3} = 3.0 * NoUnits # exponent end CarmanKozenyPermeability(args...) = CarmanKozenyPermeability(convert.(GeoUnit, args)...) -isdimensional(s::CarmanKozenyPermeability) = isdimensional(s.c) +# isdimensional(s::CarmanKozenyPermeability) = isdimensional(s.c) function param_info(s::CarmanKozenyPermeability) return MaterialParamsInfo(; Equation = L"k = c \left(\frac{\phi}{\phi_0}\right)^n") diff --git a/src/aliases.jl b/src/aliases.jl index 612ba50b..20282e7e 100644 --- a/src/aliases.jl +++ b/src/aliases.jl @@ -8,6 +8,7 @@ const VALID_KWARGS = ( :latent_heat, :radioactive_heat, :shearheating, + :permeability, :pwave_velocity, :swave_velocity, :meltfraction, diff --git a/test/test_Permeability.jl b/test/test_Permeability.jl new file mode 100644 index 00000000..34fefb96 --- /dev/null +++ b/test/test_Permeability.jl @@ -0,0 +1,79 @@ +using Test, GeoParams, Unitful, StaticArrays, LaTeXStrings + +@testset "Permeability.jl" begin + + # Set alias for permeability function + if !isdefined(Main, :GeoParamsAliases) + eval(:(@use GeoParamsAliases permeability = k)) + end + + # Make sure that structs are isbits + x = ConstantPermeability() + @test isbits(x) + @test param_info(x).Equation === L"k = cst" + @test isdimensional(x) === true + + x = HazenPermeability() + @test isbits(x) + @test param_info(x).Equation === L"k = C \cdot D_{10}^2" + @test isdimensional(x) === true + + x = PowerLawPermeability() + @test isbits(x) + @test param_info(x).Equation === L"k = c \cdot k_0 \cdot \phi^n" + @test isdimensional(x) === true + + x = CarmanKozenyPermeability() + @test isbits(x) + @test param_info(x).Equation === L"k = c \left(\frac{\phi}{\phi_0}\right)^n" + + # Test the permeability calculations with units + x1 = ConstantPermeability(; k=1e-12m^2) + @test x1.k.val == 1e-12 + @test GeoParams.get_k(x1) == 1e-12 + + x2 = HazenPermeability(; C=1.0, D10=1e-4m) + args = (;ϕ=0.4) + @test compute_permeability(x2, args) ≈ 1.0 * (0.1e-3)^2 + @test x2() ≈ 1.0 * (0.1e-3)^2 + + x3 = PowerLawPermeability(; c=1.0, k0=1e-12m^2, n=3.0) + args = (;ϕ=0.4) + @test compute_permeability(x3, args) ≈ 1.0 * 1.0e-12 * (0.4)^3 + @test x3(args) ≈ 1.0 * 1.0e-12 * (0.4)^3 + + x4 = CarmanKozenyPermeability(; c=1.0, ϕ0=0.3, n=3.0) + args = (;ϕ=0.4) + @test compute_permeability(x4, args) ≈ 1.0 * (0.4 / 0.3)^3 + @test x4(args) ≈ 1.0 * (0.4 / 0.3)^3 + + # Test the permeability calculations with non-dimensionalized units + CharUnits_GEO = GEO_units(; viscosity=1e19, length=1000km) + x1 = nondimensionalize(x1, CharUnits_GEO) + @test x1.k.val ≈ 1e-24 + + x2 = nondimensionalize(x2, CharUnits_GEO) + @test x2.C.val ≈ 1.0 + @test x2.D10.val ≈ 1e-10 + + x3 = nondimensionalize(x3, CharUnits_GEO) + @test x3.k0.val ≈ 1e-24 + + x4 = nondimensionalize(x4, CharUnits_GEO) + @test x4.c.val ≈ 1.0 + + # # Test allocations + # k = [0.0] + # ϕ = 0.4 + # args = (ϕ=ϕ) + + # # Test allocations using k alias + # k!(k, x3, args) + # num_alloc = @allocated k!(k, x3, args) + # @test num_alloc == 0 + + # k!(k, x4, args) + # num_alloc = @allocated k!(k, x4, args) + # @test num_alloc == 0 + +end From 6f6390a36f8820331a06f48baf60850b6b18e6c9 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Fri, 11 Oct 2024 06:57:18 +0200 Subject: [PATCH 05/33] update readme --- README.md | 60 +++++++++++++++++++++++++++++++++---------------------- 1 file changed, 36 insertions(+), 24 deletions(-) diff --git a/README.md b/README.md index 166c2f03..fb6d71dc 100644 --- a/README.md +++ b/README.md @@ -19,15 +19,20 @@ We also implement some typically used creep law parameters, together with tools NOTE: The package remains under development and the API is not yet fully fixed. Therefore feel free to look at it, but be aware that things may still change when you incorporate it into your codes. Comments/ideas/suggestions are highly apprecciated! ### Contents -* [1. Nondimensionalization](#1-nondimensionalization) -* [2. Material parameters](#2-material-parameters) -* [3. Plotting and output](#3-plotting-and-output) -* [4. Computational engine](#4-computational-engine) -* [5. Installation](#5-installation) -* [6. Documentation](#6-documentation) -* [7. Dependencies](#7-dependencies) -* [8. Contributing](#8-contributing) -* [9. Funding](#9-funding) +- [Contents](#contents) +- [1. Nondimensionalization](#1-nondimensionalization) +- [2. Material parameters](#2-material-parameters) + - [2.1 Constant density, constant linear viscosity](#21-constant-density-constant-linear-viscosity) + - [2.2 Nonlinear creep laws](#22-nonlinear-creep-laws) +- [3. Plotting and output](#3-plotting-and-output) + - [3.1 Plotting](#31-plotting) + - [3.2 Automatically create data tables](#32-automatically-create-data-tables) +- [4. Computational engine](#4-computational-engine) +- [5. Installation](#5-installation) +- [6. Documentation](#6-documentation) +- [7. Dependencies](#7-dependencies) +- [8. Contributing](#8-contributing) +- [9. Funding](#9-funding) ### 1. Nondimensionalization Typical geodynamic simulations involve dimensions on the order of 10's-1000's of kilometers, and viscosities on the order of ~1e20 Pas. If such values are directly employed in numerical solvers, they may result in roundoff errors. It is therefore common practice to nondimensionalize the input parameters by dividing them by typical values such that the result gives numbers that are closer to one. @@ -126,23 +131,30 @@ Phase 2 : Viscous Sinker You can add pre-defined non-linear creep laws as: ```julia julia> Phase = SetMaterialParams(Name="Viscous Matrix", Phase=2, - Density = ConstantDensity(), - CreepLaws = (SetDislocationCreep("Wet Olivine | Hirth & Kohlstedt (2003)"), - LinearViscous(η=1e23Pa*s)) ) + Density = ConstantDensity(), + CompositeRheology = CompositeRheology( + SetDislocationCreep(Dislocation.wet_olivine_Hirth_2003), + LinearViscous(η=1e23Pa*s)) + ) Phase 2 : Viscous Matrix | [dimensional units] - | - |-- Density : Constant density: ρ=2900.0 kg m⁻³·⁰ - |-- Gravity : Gravitational acceleration: g=9.81 m s⁻²·⁰ - |-- CreepLaws : DislocationCreep: Name = Wet Olivine | Hirth & Kohlstedt (2003), n=3.5, r=1.2, A=90.0, E=480.0, V=1.1e-5, Apparatus=1 - | Linear viscosity: η=1.0e23 -``` -Note that the dictionary `DislocationCreep_info` has all pre-defined creep laws, so for an overview type: -```julia -julia> DislocationCreep_info -Dict{String, DislocationCreep} with 2 entries: - "Dry Olivine | Hirth & Kohlstedt (2003)" => DislocationCreep: n=3.05, r=0, A=110000.0 MPa⁻³·⁰⁵ s⁻¹, E… - "Wet Olivine | Hirth & Kohlstedt (2003)" => DislocationCreep: n=3.5, r=1.2, A=90 MPa⁻³·⁵ s⁻¹, E=480 k… + | + |-- Name : Viscous Matrix + |-- Density : Constant density: ρ=2900.0 kg m⁻³·⁰ + |-- Gravity : Gravitational acceleration: g=9.81 m s⁻²·⁰ + |-- CompositeRheology : --⟦▪̲̅▫̲̅▫̲̅▫̲̅----⟦▪̲̅▫̲̅▫̲̅▫̲̅-- +``` +Note that the functions `dislocation_law_list()` and `diffusion_law_list()` list all pre-defined creep laws, so for an overview type: +```julia +julia> dislocation_law_list() +40-element Vector{Function}: + dry_anorthite_Rybacki_2000 (generic function with 1 method) + dry_olivine_Hirth_2003 (generic function with 1 method) + dry_olivine_Karato_2003 (generic function with 1 method) + dry_quartzite_Jaoul_1984 (generic function with 1 method) + ⋮ + wet_omphacite_Zhang_2006 (generic function with 1 method) + wet_quartzite_Hirth_2001 (generic function with 1 method) ``` ### 3. Plotting and output From 4f61db8f4f5b03de90a11a265e434aa8d1b9bb4e Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Fri, 11 Oct 2024 07:03:36 +0200 Subject: [PATCH 06/33] more fixes --- README.md | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index fb6d71dc..ef968ff1 100644 --- a/README.md +++ b/README.md @@ -176,9 +176,10 @@ When writing scientific papers that describes numerical modelling results, it is That is why we provide routines that fully automatize this process: First, we need to define a phase. ```julia -julia> MatParam = (SetMaterialParams(Name="Viscous Matrix", Phase=1, Density=ConstantDensity(),CreepLaws = SetDislocationCreep("Quartz Diorite | Hansen & Carter (1982)")), +julia> import GeoParams.Dislocation +julia> MatParam = (SetMaterialParams(Name="Viscous Matrix", Phase=1, Density=ConstantDensity(),CreepLaws = SetDislocationCreep(Dislocation.quartz_diorite_HansenCarter_1982)), SetMaterialParams(Name="Viscous Sinker", Phase=2, Density= PT_Density(),CreepLaws = LinearViscous(η=1e21Pa*s)), - SetMaterialParams(Name="Viscous Bottom", Phase=3, Density= PT_Density(),CreepLaws = SetDislocationCreep("Diabase | Caristan (1982)"))) + SetMaterialParams(Name="Viscous Bottom", Phase=3, Density= PT_Density(),CreepLaws = SetDislocationCreep(Dislocation.diabase_Caristan_1982))) ``` Next, you can create a LaTeX table for the defined phase ... ```julia From 72eb6664fc26a62c67a624bf2cb996d98cb98aff Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Fri, 11 Oct 2024 07:08:46 +0200 Subject: [PATCH 07/33] update CI --- .github/workflows/blank.yml | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/.github/workflows/blank.yml b/.github/workflows/blank.yml index 3514c836..75f5c71e 100644 --- a/.github/workflows/blank.yml +++ b/.github/workflows/blank.yml @@ -16,13 +16,13 @@ jobs: matrix: version: - '1.7' - - '1.8' - - '1.9' - '1.10' + - '1.11' - 'nightly' os: - ubuntu-latest - macOS-latest + - macOS-13 - windows-latest arch: - x64 @@ -48,6 +48,7 @@ jobs: - uses: codecov/codecov-action@v4 env: CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }} + DISPLAY: 0 build_docs: runs-on: ubuntu-20.04 steps: From ac8ba39d6ca18d4d67489175b0d2a1662d7ede13 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Fri, 11 Oct 2024 07:16:55 +0200 Subject: [PATCH 08/33] try deactivating DISPLAY in CI --- .github/workflows/blank.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/blank.yml b/.github/workflows/blank.yml index 75f5c71e..626ff629 100644 --- a/.github/workflows/blank.yml +++ b/.github/workflows/blank.yml @@ -43,7 +43,7 @@ jobs: ${{ runner.os }}-test- ${{ runner.os }}- - uses: julia-actions/julia-buildpkg@latest - - uses: julia-actions/julia-runtest@latest + - run: DISPLAY=:0 xvfb-run -s '-screen 0 1024x768x24' julia --project=@. -e 'using Pkg; Pkg.test(coverage=true)' - uses: julia-actions/julia-processcoverage@v1 - uses: codecov/codecov-action@v4 env: From 792c6dcb630558d3d5034c8c126663a0ea7143be Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Fri, 11 Oct 2024 07:20:34 +0200 Subject: [PATCH 09/33] update CI --- .github/workflows/blank.yml | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/.github/workflows/blank.yml b/.github/workflows/blank.yml index 626ff629..7d11cc28 100644 --- a/.github/workflows/blank.yml +++ b/.github/workflows/blank.yml @@ -7,6 +7,11 @@ on: tags: '*' pull_request: +# Cancel redundant CI tests automatically +concurrency: + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: true + jobs: run_tests: name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} @@ -43,7 +48,8 @@ jobs: ${{ runner.os }}-test- ${{ runner.os }}- - uses: julia-actions/julia-buildpkg@latest - - run: DISPLAY=:0 xvfb-run -s '-screen 0 1024x768x24' julia --project=@. -e 'using Pkg; Pkg.test(coverage=true)' + - name: Run tests + run: DISPLAY=:0 xvfb-run -s '-screen 0 1024x768x24' julia --project=@. -e 'using Pkg; Pkg.test(coverage=true)' - uses: julia-actions/julia-processcoverage@v1 - uses: codecov/codecov-action@v4 env: From 977affb276a945a30f730e67344252972e85f0a0 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Fri, 11 Oct 2024 07:25:31 +0200 Subject: [PATCH 10/33] add packages --- .github/workflows/blank.yml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.github/workflows/blank.yml b/.github/workflows/blank.yml index 7d11cc28..ecdb0bd7 100644 --- a/.github/workflows/blank.yml +++ b/.github/workflows/blank.yml @@ -48,6 +48,8 @@ jobs: ${{ runner.os }}-test- ${{ runner.os }}- - uses: julia-actions/julia-buildpkg@latest + - name: Add packages + run: sudo apt-get update && sudo apt-get install -y xorg-dev mesa-utils xvfb libgl1 freeglut3-dev libxrandr-dev libxinerama-dev libxcursor-dev libxi-dev libxext-dev - name: Run tests run: DISPLAY=:0 xvfb-run -s '-screen 0 1024x768x24' julia --project=@. -e 'using Pkg; Pkg.test(coverage=true)' - uses: julia-actions/julia-processcoverage@v1 From cd97cfc7a8e12ea336811287449ad97c1ccd40b4 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Fri, 11 Oct 2024 07:30:11 +0200 Subject: [PATCH 11/33] update CI --- .github/workflows/blank.yml | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/.github/workflows/blank.yml b/.github/workflows/blank.yml index ecdb0bd7..b5785673 100644 --- a/.github/workflows/blank.yml +++ b/.github/workflows/blank.yml @@ -48,15 +48,13 @@ jobs: ${{ runner.os }}-test- ${{ runner.os }}- - uses: julia-actions/julia-buildpkg@latest - - name: Add packages - run: sudo apt-get update && sudo apt-get install -y xorg-dev mesa-utils xvfb libgl1 freeglut3-dev libxrandr-dev libxinerama-dev libxcursor-dev libxi-dev libxext-dev - - name: Run tests - run: DISPLAY=:0 xvfb-run -s '-screen 0 1024x768x24' julia --project=@. -e 'using Pkg; Pkg.test(coverage=true)' + - uses: julia-actions/julia-runtest@latest + env: + DISPLAY: 0 - uses: julia-actions/julia-processcoverage@v1 - uses: codecov/codecov-action@v4 env: CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }} - DISPLAY: 0 build_docs: runs-on: ubuntu-20.04 steps: From aad2443aae3063f231eaf5cf38968aa056ced592 Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Sat, 12 Oct 2024 12:47:42 +0200 Subject: [PATCH 12/33] update --- src/GeoParams.jl | 4 +- src/Permeability/Permeability.jl | 116 +++++++++++++++++++++++++++++-- src/Viscosity/BulkViscosity.jl | 20 ++++++ test/test_Permeability.jl | 4 +- 4 files changed, 134 insertions(+), 10 deletions(-) create mode 100644 src/Viscosity/BulkViscosity.jl diff --git a/src/GeoParams.jl b/src/GeoParams.jl index aadc4929..116674b8 100644 --- a/src/GeoParams.jl +++ b/src/GeoParams.jl @@ -344,7 +344,9 @@ export compute_permeability, ConstantPermeability, HazenPermeability, PowerLawPermeability, - CarmanKozenyPermeability + CarmanKozenyPermeability, + BiotWillis, + SkemptonCoeff include("Traits/rheology.jl") export RheologyTrait diff --git a/src/Permeability/Permeability.jl b/src/Permeability/Permeability.jl index 1fbfa662..0d565266 100644 --- a/src/Permeability/Permeability.jl +++ b/src/Permeability/Permeability.jl @@ -4,6 +4,7 @@ using Parameters, Unitful, LaTeXStrings, MuladdMacro using ..Units using ..PhaseDiagrams using GeoParams: AbstractMaterialParam, AbstractMaterialParamsStruct, @extractors, add_extractor_functions +using GeoParams: fastpow, pow_check, @pow import ..Units: isdimensional using ..MaterialParameters: No_MaterialParam, MaterialParamsInfo import Base.show, GeoParams.param_info @@ -19,7 +20,9 @@ export compute_permeability, # calculation routines ConstantPermeability, # constant HazenPermeability, # Hazen equation PowerLawPermeability, # Power-law permeability - CarmanKozenyPermeability # Carman-Kozeny permeability + CarmanKozenyPermeability, # Carman-Kozeny permeability + BiotWillis, # Biot-Willis coefficient + SkemptonCoeff # Skempton coefficient # Define "empty" computational routines in case nothing is defined function compute_permeability!( @@ -32,6 +35,25 @@ function compute_permeability(s::No_MaterialParam{_T}; ϕ::_T=zero(_T)) where {_ end # Constant Permeability +""" + ConstantPermeability(k = 1e-12m^2) + +Defines a constant permeability value for a given material. + +# Arguments +- `k::Float64`: The permeability value in square meters (m^2). Default is `1e-12 m^2`. + +# Example +```julia +rheology = SetMaterialParams(; + Phase=1, + CreepLaws=(PowerlawViscous(), LinearViscous(; η=1e21Pa * s)), + Gravity=ConstantGravity(; g=9.81.0m / s^2), + Density= MeltDependent_Density(), + Permeability = ConstantPermeability(; k=1e-12m^2), + ) +``` +""" @with_kw_noshow struct ConstantPermeability{_T,U} <: AbstractPermeability{_T} k::GeoUnit{_T,U} = 1e-12m^2 # permeability end @@ -78,7 +100,7 @@ function (s::HazenPermeability{_T})(; kwargs...) where {_T} @unpack_val C, D10 = s end - return C * D10^2 + return @pow C * D10^2 end @inline (s::HazenPermeability)(args) = s(; args...) @@ -94,7 +116,7 @@ end c::GeoUnit{_T,U1} = 1.0 * NoUnits # Power-law constant k0::GeoUnit{_T,U2} = 1e-12 * m^2 # reference permeability ϕ::GeoUnit{_T,U3} = 1e-2 * NoUnits # reference porosity - n::GeoUnit{_T,U4} = 3.0 * NoUnits # exponent + n::GeoUnit{_T,U4} = 3 * NoUnits # exponent end PowerLawPermeability(args...) = PowerLawPermeability(convert.(GeoUnit, args)...) isdimensional(s::PowerLawPermeability) = isdimensional(s.k0) @@ -110,7 +132,7 @@ function (s::PowerLawPermeability{_T})(; ϕ=0e0, kwargs...) where {_T} @unpack_val c, k0, n = s end - return c * k0 * ϕ^n + return @pow c * k0 * ϕ^n end @inline (s::PowerLawPermeability)(args) = s(; args...) @@ -122,9 +144,9 @@ end # Carman-Kozeny Permeability @with_kw_noshow struct CarmanKozenyPermeability{_T,U1,U2,U3} <: AbstractPermeability{_T} - c::GeoUnit{_T,U1} = 1.0 * NoUnits # Carman-Kozeny constant + c::GeoUnit{_T,U1} = 1.0 * m^2 # Carman-Kozeny constant ϕ0::GeoUnit{_T,U2} = 0.01 * NoUnits # reference porosity - n::GeoUnit{_T,U3} = 3.0 * NoUnits # exponent + n::GeoUnit{_T,U3} = 3 * NoUnits # exponent end CarmanKozenyPermeability(args...) = CarmanKozenyPermeability(convert.(GeoUnit, args)...) # isdimensional(s::CarmanKozenyPermeability) = isdimensional(s.c) @@ -140,7 +162,7 @@ function (s::CarmanKozenyPermeability{_T})(; ϕ=1e-2, kwargs...) where {_T} @unpack_val c, ϕ0, n = s end - return c * (ϕ / ϕ0)^n + return @pow c * (ϕ / ϕ0)^n end @inline (s::CarmanKozenyPermeability)(args) = s(; args...) @@ -150,6 +172,86 @@ function show(io::IO, g::CarmanKozenyPermeability) return print(io, "Carman-Kozeny permeability: k = c * (ϕ / ϕ0)^n; c=$(g.c); ϕ0=$(g.ϕ0); n=$(g.n)") end +# ----------------------------------------------- +# This implements the methods described by Yarushina and Podladchikov, 2015 +""" +Biot-Willis coefficient +""" +@with_kw_noshow struct BiotWillis{_T,U} <: AbstractPermeability{_T} + Kd::GeoUnit{_T,U} = 1.0 * Pa # Drained bulk modulus + Ks::GeoUnit{_T,U} = 1.0 * Pa # Solid grain bulk modulus + α::GeoUnit{_T,U} = 1.0 * NoUnits # Biot-Willis coefficient to keep track of dimensionality +end + +BiotWillis(args...) = BiotWillis(convert.(GeoUnit, args)...) +isdimensional(s::BiotWillis) = isdimensional(s.Kd) + +function param_info(s::BiotWillis) + return MaterialParamsInfo(; Equation = L"\alpha = 1 - \frac{K_d}{K_s}") +end + +function (bw::BiotWillis{_T})(; kwargs...) where {_T} + if kwargs isa Quantity + @unpack_units Kd, Ks, α = bw + else + @unpack_val Kd, Ks, α = bw + end + + return 1 - (Kd * inv(Ks)) +end + +@inline (s::BiotWillis)(args) = s(; args...) +@inline compute_biot_willis(s::BiotWillis, args) = s(args) + +function show(io::IO, g::BiotWillis) + return print(io, "Biot-Willis coefficient: α = 1 - (Kd / Ks); Kd=$(g.Kd); Ks=$(g.Ks)") +end + +""" +Skempton coefficient +""" +@with_kw_noshow struct SkemptonCoeff{_T,U} <: AbstractPermeability{_T} + Kd::GeoUnit{_T,U} = 1.0 * Pa # Drained bulk modulus + Ks::GeoUnit{_T,U} = 1.0 * Pa # Solid grain bulk modulus + Kf::GeoUnit{_T,U} = 1.0 * Pa # Fluid bulk modulus + ϕ::GeoUnit{_T,U} = 1.0 * NoUnits # Porosity + B::GeoUnit{_T,U} = 1.0 * Pa # Skempton coefficient to keep track of dimensionality +end + +SkemptonCoeff(args...) = SkemptonCoeff(convert.(GeoUnit, args)...) +isdimensional(s::SkemptonCoeff) = isdimensional(s.Kd) + +function param_info(s::SkemptonCoeff) + return MaterialParamsInfo(; Equation = L"B = \frac{(1/K_d - 1/K_s)}{(1/K_d - 1/K_s) + \phi (1/K_f - 1/K_s)}") +end + +function (sc::SkemptonCoeff{_T})(; kwargs...) where {_T} + if kwargs isa Quantity + @unpack_units Kd, Ks, Kf, ϕ, B = sc + else + @unpack_val Kd, Ks, Kf, ϕ, B = sc + end + + return (inv(Kd) - inv(Ks)) * inv(inv(Kd) - inv(Ks) + ϕ * (inv(Kf) - inv(Ks))) +end + +@inline (s::SkemptonCoeff)(args) = s(; args...) +@inline compute_skempton_coeff(s::SkemptonCoeff, args) = s(args) + +function show(io::IO, g::SkemptonCoeff) + return print(io, "Skempton Coefficient: B = (1/Kd - 1/Ks) / ((1/Kd - 1/Ks) + ϕ * (1/Kf - 1/Ks)); Kd=$(g.Kd); Ks=$(g.Ks); Kf=$(g.Kf); ϕ=$(g.ϕ)") +end + +""" + compute_permeability!(k::AbstractArray{_T, N}, MatParam::NTuple{K,AbstractMaterialParamsStruct}, PhaseRatios::AbstractArray{_T, M}, P=nothing, T=nothing) + +In-place computation of permeability `k` for the whole domain and all phases, in case a vector with phase properties `MatParam` is provided, along with `P` and `T` arrays. +This assumes that the `PhaseRatio` of every point is specified as an Integer in the `PhaseRatios` array, which has one dimension more than the data arrays (and has a phase fraction between 0-1) +""" +@inline compute_permeability!(args::Vararg{Any, N}) where N = compute_param!(compute_permeability, args...) +@inline compute_permeability(args::Vararg{Any, N}) where N = compute_param(compute_permeability, args...) +@inline compute_permeability_ratio(args::Vararg{Any, N}) where N = compute_param_times_frac(compute_permeability, args...) + # extractor methods for type in (ConstantPermeability, HazenPermeability, PowerLawPermeability, CarmanKozenyPermeability) @extractors(type, :Permeability) diff --git a/src/Viscosity/BulkViscosity.jl b/src/Viscosity/BulkViscosity.jl new file mode 100644 index 00000000..7774c389 --- /dev/null +++ b/src/Viscosity/BulkViscosity.jl @@ -0,0 +1,20 @@ +""" +bulk_viscosity(η_shear, ϕ, C, R, λ, P_eff) + +# η_shear: shear viscosity +# ϕ: porosity +# C: shear viscosity ratio (value, dimensionless) +# R: compaction strength ratio (value, dimensionless) +# λ: effective pressure transition zone +# P_eff: effective Pressure +""" +function bulk_viscosity(η_shear, ϕ, C, R, λ, P_eff) + + η_bulk = η_shear / ϕ * C * (1.0 + 0.5 * (inv(R) - 1.0) * (1.0 + tanh(-P_eff * inv(λ)))) + + return η_bulk +end + + +### Georg Adjoint 2 pahse paper, +### Ludo's 2019 paper has it quite nicely also with a calc for R diff --git a/test/test_Permeability.jl b/test/test_Permeability.jl index 34fefb96..7704389e 100644 --- a/test/test_Permeability.jl +++ b/test/test_Permeability.jl @@ -37,12 +37,12 @@ using Test, GeoParams, Unitful, StaticArrays, LaTeXStrings @test compute_permeability(x2, args) ≈ 1.0 * (0.1e-3)^2 @test x2() ≈ 1.0 * (0.1e-3)^2 - x3 = PowerLawPermeability(; c=1.0, k0=1e-12m^2, n=3.0) + x3 = PowerLawPermeability(; c=1.0, k0=1e-12m^2, n=3) args = (;ϕ=0.4) @test compute_permeability(x3, args) ≈ 1.0 * 1.0e-12 * (0.4)^3 @test x3(args) ≈ 1.0 * 1.0e-12 * (0.4)^3 - x4 = CarmanKozenyPermeability(; c=1.0, ϕ0=0.3, n=3.0) + x4 = CarmanKozenyPermeability(; c=1.0, ϕ0=0.3, n=3) args = (;ϕ=0.4) @test compute_permeability(x4, args) ≈ 1.0 * (0.4 / 0.3)^3 @test x4(args) ≈ 1.0 * (0.4 / 0.3)^3 From adc98ff1e8df1d1af08cb962158a5adc6798db58 Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Sat, 12 Oct 2024 12:49:58 +0200 Subject: [PATCH 13/33] adapt melting to other similar functions --- src/MeltFraction/MeltingParameterization.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/MeltFraction/MeltingParameterization.jl b/src/MeltFraction/MeltingParameterization.jl index 1dfae569..6a688760 100644 --- a/src/MeltFraction/MeltingParameterization.jl +++ b/src/MeltFraction/MeltingParameterization.jl @@ -747,7 +747,7 @@ end # Computational routines needed for computations with the MaterialParams structure function compute_meltfraction(s::AbstractMaterialParamsStruct, args) if isempty(s.Melting) #in case there is a phase with no melting parametrization - return zero(typeof(args).types[1]) + return isempty(args) ? 0.0 : zero(typeof(args).types[1]) # return zero if not specified else return compute_meltfraction(s.Melting[1], args) end From 4a46439319535bcc5e33eb1a1a80db19510001c8 Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Sat, 12 Oct 2024 13:00:04 +0200 Subject: [PATCH 14/33] fix? --- src/MeltFraction/MeltingParameterization.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/MeltFraction/MeltingParameterization.jl b/src/MeltFraction/MeltingParameterization.jl index 6a688760..91210237 100644 --- a/src/MeltFraction/MeltingParameterization.jl +++ b/src/MeltFraction/MeltingParameterization.jl @@ -747,7 +747,7 @@ end # Computational routines needed for computations with the MaterialParams structure function compute_meltfraction(s::AbstractMaterialParamsStruct, args) if isempty(s.Melting) #in case there is a phase with no melting parametrization - return isempty(args) ? 0.0 : zero(typeof(args).types[1]) # return zero if not specified + return 0.0 # return zero if not specified else return compute_meltfraction(s.Melting[1], args) end From 6769629cc09ee3b36839ca64594bdcbc17f64a0e Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Sat, 12 Oct 2024 13:09:14 +0200 Subject: [PATCH 15/33] test more consisten option --- src/MeltFraction/MeltingParameterization.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/MeltFraction/MeltingParameterization.jl b/src/MeltFraction/MeltingParameterization.jl index 91210237..c5bcdd21 100644 --- a/src/MeltFraction/MeltingParameterization.jl +++ b/src/MeltFraction/MeltingParameterization.jl @@ -747,7 +747,7 @@ end # Computational routines needed for computations with the MaterialParams structure function compute_meltfraction(s::AbstractMaterialParamsStruct, args) if isempty(s.Melting) #in case there is a phase with no melting parametrization - return 0.0 # return zero if not specified + return zero(typeof(args)) # return zero if not specified else return compute_meltfraction(s.Melting[1], args) end From 48a8f0349ef546e61b1804966e2eeca2bae65359 Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Sat, 12 Oct 2024 13:13:55 +0200 Subject: [PATCH 16/33] revert back to hardcoded 0.0 as in SH --- src/MeltFraction/MeltingParameterization.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/MeltFraction/MeltingParameterization.jl b/src/MeltFraction/MeltingParameterization.jl index c5bcdd21..c0295f06 100644 --- a/src/MeltFraction/MeltingParameterization.jl +++ b/src/MeltFraction/MeltingParameterization.jl @@ -747,7 +747,7 @@ end # Computational routines needed for computations with the MaterialParams structure function compute_meltfraction(s::AbstractMaterialParamsStruct, args) if isempty(s.Melting) #in case there is a phase with no melting parametrization - return zero(typeof(args)) # return zero if not specified + return 0.0 # return zero if not specified else return compute_meltfraction(s.Melting[1], args) end @@ -755,7 +755,7 @@ end function compute_dϕdT(s::AbstractMaterialParamsStruct, args) if isempty(s.Melting) #in case there is a phase with no melting parametrization - return zero(typeof(args).types[1]) + return 0.0 # return zero if not specified else return compute_dϕdT(s.Melting[1], args) end From 33b9368b5c70ad9c958f733f34706b676e8eb88a Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Sat, 12 Oct 2024 13:29:47 +0200 Subject: [PATCH 17/33] number format --- src/Energy/Shearheating.jl | 2 +- src/MeltFraction/MeltingParameterization.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/Energy/Shearheating.jl b/src/Energy/Shearheating.jl index ed5df37c..876a8d81 100644 --- a/src/Energy/Shearheating.jl +++ b/src/Energy/Shearheating.jl @@ -172,7 +172,7 @@ end function compute_shearheating(s::AbstractMaterialParamsStruct, args::Vararg{Any,N}) where N if isempty(s.ShearHeat) - return 0.0 # return zero if not specified + return 0e0 # return zero if not specified else return compute_shearheating(s.ShearHeat[1], args...) end diff --git a/src/MeltFraction/MeltingParameterization.jl b/src/MeltFraction/MeltingParameterization.jl index c0295f06..ae7847e7 100644 --- a/src/MeltFraction/MeltingParameterization.jl +++ b/src/MeltFraction/MeltingParameterization.jl @@ -747,7 +747,7 @@ end # Computational routines needed for computations with the MaterialParams structure function compute_meltfraction(s::AbstractMaterialParamsStruct, args) if isempty(s.Melting) #in case there is a phase with no melting parametrization - return 0.0 # return zero if not specified + return 0e0 # return zero if not specified else return compute_meltfraction(s.Melting[1], args) end @@ -755,7 +755,7 @@ end function compute_dϕdT(s::AbstractMaterialParamsStruct, args) if isempty(s.Melting) #in case there is a phase with no melting parametrization - return 0.0 # return zero if not specified + return 0e0 # return zero if not specified else return compute_dϕdT(s.Melting[1], args) end From a174816b5a7502263dd081351d8c76311396c05a Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Sat, 12 Oct 2024 15:03:21 +0200 Subject: [PATCH 18/33] activate GLMakie as dependency --- Project.toml | 2 -- src/GeoParams.jl | 18 +++++++++--------- src/Plotting/Plotting.jl | 13 ++++++++++++- 3 files changed, 21 insertions(+), 12 deletions(-) diff --git a/Project.toml b/Project.toml index 5e3eb77a..1e9b5119 100644 --- a/Project.toml +++ b/Project.toml @@ -15,7 +15,6 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Loess = "4345ca2d-374a-55d4-8d30-97f9976e7612" MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a" -Requires = "ae029012-a4dd-5104-9daa-d747884805df" Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665" Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" @@ -46,7 +45,6 @@ LinearAlgebra = "1.0" Loess = "0.5,0.6" MuladdMacro = "0.2.4" Parameters = "0.12" -Requires = "0.5.0, 0.6, 0.7, 0.8, 1.0 - 1.3" Roots = "1.0, 2.0" Setfield = "0.5.1, 0.6, 0.7, 0.8, 1.0" SpecialFunctions = "2.0" diff --git a/src/GeoParams.jl b/src/GeoParams.jl index dea42709..3afa8178 100644 --- a/src/GeoParams.jl +++ b/src/GeoParams.jl @@ -15,7 +15,7 @@ module GeoParams using Parameters # helps setting default parameters in structures using Unitful # Units using BibTeX # references of creep laws -using Requires: @require # To only add plotting routines if GLMakie is loaded +#using Requires: @require # To only add plotting routines if GLMakie is loaded using StaticArrays using LinearAlgebra using ForwardDiff @@ -423,14 +423,14 @@ export PlotStrainrateStress, # We do not check `isdefined(Base, :get_extension)` as recommended since # Julia v1.9.0 does not load package extensions when their dependency is # loaded from the main environment. -function __init__() - @static if !(VERSION >= v"1.9.1") - @require GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" begin - print("Adding plotting routines of GeoParams through GLMakie \n") - @eval include("../ext/GeoParamsGLMakieExt.jl") - end - end -end +#function __init__() +# @static if !(VERSION >= v"1.9.1") +# @require GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" begin +# print("Adding plotting routines of GeoParams through GLMakie \n") +# @eval include("../ext/GeoParamsGLMakieExt.jl") +# end +# end +#end #Set functions aliases using @use include("aliases.jl") diff --git a/src/Plotting/Plotting.jl b/src/Plotting/Plotting.jl index 07c0d728..f55114b2 100644 --- a/src/Plotting/Plotting.jl +++ b/src/Plotting/Plotting.jl @@ -232,10 +232,21 @@ function customize_plot!(li, args) end """ - fig,ax,τII,εII = PlotStressStrainrate(x; args=(T=1000.0, P=0.0, d=1e-3, f=1.0), Stress=(1e0,1e8), plt=nothing) + fig,ax,τII,εII = PlotStressStrainrate(x; args=(T=1000.0, P=0.0, d=1e-3, f=1.0), Stress=(1e0,1e8)) Same as `PlotStrainrateStress` but with stress (in MPa) versus strainrate (in 1/s) instead. +Example +=== + +```julia +julia> import GeoParams.Dislocation, GeoParams.Diffusion; +julia> a1=SetDislocationCreep(Dislocation.wet_olivine_Hirth_2003); +julia> a2=SetDiffusionCreep(Diffusion.wet_olivine_Hirth_2003); +julia> x=CompositeRheology(a1,a2); +julia> fig,ax,τII,εII = PlotStressStrainrate(x; args=(T=1000.0, P=0.0, d=1e-3, f=1.0), Stress=(1e0,1e8)); +``` + """ function PlotStressStrainrate( x; From 9fecb803354d32ecf05f3fafaee1125f42e734d1 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Sat, 12 Oct 2024 15:19:00 +0200 Subject: [PATCH 19/33] next attempt --- .github/workflows/blank.yml | 2 -- Project.toml | 4 ++++ src/GeoParams.jl | 13 ------------- 3 files changed, 4 insertions(+), 15 deletions(-) diff --git a/.github/workflows/blank.yml b/.github/workflows/blank.yml index b5785673..51e36e74 100644 --- a/.github/workflows/blank.yml +++ b/.github/workflows/blank.yml @@ -49,8 +49,6 @@ jobs: ${{ runner.os }}- - uses: julia-actions/julia-buildpkg@latest - uses: julia-actions/julia-runtest@latest - env: - DISPLAY: 0 - uses: julia-actions/julia-processcoverage@v1 - uses: codecov/codecov-action@v4 env: diff --git a/Project.toml b/Project.toml index 1e9b5119..a1cac6b7 100644 --- a/Project.toml +++ b/Project.toml @@ -60,3 +60,7 @@ julia = "1.7" [extras] GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["Test"] diff --git a/src/GeoParams.jl b/src/GeoParams.jl index 3afa8178..61f14fbc 100644 --- a/src/GeoParams.jl +++ b/src/GeoParams.jl @@ -15,7 +15,6 @@ module GeoParams using Parameters # helps setting default parameters in structures using Unitful # Units using BibTeX # references of creep laws -#using Requires: @require # To only add plotting routines if GLMakie is loaded using StaticArrays using LinearAlgebra using ForwardDiff @@ -420,18 +419,6 @@ export PlotStrainrateStress, PlotPressureStressTime_0D, StrengthEnvelopePlot -# We do not check `isdefined(Base, :get_extension)` as recommended since -# Julia v1.9.0 does not load package extensions when their dependency is -# loaded from the main environment. -#function __init__() -# @static if !(VERSION >= v"1.9.1") -# @require GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" begin -# print("Adding plotting routines of GeoParams through GLMakie \n") -# @eval include("../ext/GeoParamsGLMakieExt.jl") -# end -# end -#end - #Set functions aliases using @use include("aliases.jl") export ntuple_idx From a5428f6a0ea5d6e3bd1c0324b9be108da9a5e8ce Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Sat, 12 Oct 2024 16:27:02 +0200 Subject: [PATCH 20/33] forgot to remove GLMakie from /test/GLMakie --- Project.toml | 2 +- ext/GeoParamsGLMakieExt.jl | 2 ++ test/Project.toml | 1 - test/runtests.jl | 4 +++- 4 files changed, 6 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index a1cac6b7..c7841a32 100644 --- a/Project.toml +++ b/Project.toml @@ -63,4 +63,4 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test"] +test = ["Test", "Statistics"] diff --git a/ext/GeoParamsGLMakieExt.jl b/ext/GeoParamsGLMakieExt.jl index e2597c3d..92633e66 100644 --- a/ext/GeoParamsGLMakieExt.jl +++ b/ext/GeoParamsGLMakieExt.jl @@ -1,6 +1,8 @@ # Package extension for adding GLMakie-based features to GeoParams.jl module GeoParamsGLMakieExt +using GeoParams + # We do not check `isdefined(Base, :get_extension)` as recommended since # Julia v1.9.0 does not load package extensions when their dependency is # loaded from the main environment. diff --git a/test/Project.toml b/test/Project.toml index c23eaddf..f8a39ef5 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -2,7 +2,6 @@ AbstractDifferentiation = "c29ec348-61ec-40c8-8164-b8c60e9d9f3d" DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" -GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" diff --git a/test/runtests.jl b/test/runtests.jl index 20a32f9b..ae410c4f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,10 +1,11 @@ using Test, Statistics + function runtests() files = readdir(@__DIR__) test_files = filter(startswith("test_"), files) - for f in test_files + for f in test_files[1] if !isdir(f) include(f) end @@ -12,3 +13,4 @@ function runtests() end runtests() + From 280416204cfd41d660df97551e81b43a785881a6 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Sat, 12 Oct 2024 16:31:49 +0200 Subject: [PATCH 21/33] remove debugging comment --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index ae410c4f..a1d1ac00 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,7 +5,7 @@ function runtests() files = readdir(@__DIR__) test_files = filter(startswith("test_"), files) - for f in test_files[1] + for f in test_files if !isdir(f) include(f) end From 90b5d5e85a750ef20f2d1a342b6066f05fb03fe7 Mon Sep 17 00:00:00 2001 From: Pascal Aellig <93140290+aelligp@users.noreply.github.com> Date: Sat, 12 Oct 2024 16:45:43 +0200 Subject: [PATCH 22/33] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index c7841a32..ec9df9ce 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "GeoParams" uuid = "e018b62d-d9de-4a26-8697-af89c310ae38" authors = ["Boris Kaus , Albert De Montserrat "] -version = "0.6.4" +version = "0.6.5" [deps] BibTeX = "7b0aa2c9-049f-5cec-894a-2b6b781bb25e" From d733d217eee1f633c39147ec2c465de9ffba6fd9 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Sat, 12 Oct 2024 17:14:28 +0200 Subject: [PATCH 23/33] fix plots --- src/Plotting/Plotting.jl | 102 ++++++++++++++++++++++++++++++--------- 1 file changed, 80 insertions(+), 22 deletions(-) diff --git a/src/Plotting/Plotting.jl b/src/Plotting/Plotting.jl index f55114b2..d621642b 100644 --- a/src/Plotting/Plotting.jl +++ b/src/Plotting/Plotting.jl @@ -45,12 +45,13 @@ Plots deviatoric stress versus deviatoric strain rate for a single or multiple c First, we retrieve the data for anorthite creeplaws ```julia-repl -julia> pp = SetDiffusionCreep("Dry Anorthite | Rybacki et al. (2006)"); -julia> pp1 = SetDislocationCreep("Dry Anorthite | Rybacki et al. (2006)"); +julia> import GeoParams.Diffusion, GeoParams.Dislocation +julia> pp = SetDiffusionCreep(Diffusion.dry_anorthite_Rybacki_2006); +julia> pp1 = SetDislocationCreep(Dislocation.dry_anorthite_Rybacki_2006); ``` -Next you can define each of the creeplaws inidvidually, plus a combined diffusion & dislocation creep law: +Next you can define each of the creeplaws individually, plus a combined diffusion & dislocation creep law: ```julia-repl -julia> v = (pp,pp1,(pp,pp1)); +julia> v = (pp,pp1,CompositeRheology(pp,pp1)); ``` Next, define temperature to be `900K` and grainsize to be `100 μm` and create a default plot of the 3 mechanisms: ```julia-repl @@ -531,7 +532,7 @@ function PlotStressViscosity( end """ - T,Cp,plt = PlotHeatCapacity(Cp::AbstractHeatCapacity; T=nothing, plt=nothing, lbl=nothing) + fig,ax,T,Cp_vec = PlotHeatCapacity(Cp::AbstractHeatCapacity; T=nothing, plt=nothing, lbl=nothing) Creates a plot of temperature `T` vs. heat capacity, as specified in Cp (which can be temperature-dependent). @@ -543,36 +544,93 @@ Creates a plot of temperature `T` vs. heat capacity, as specified in Cp (which c # Example ``` julia> Cp = T_HeatCapacity_Whittacker() -julia> T,Cp,plt = PlotHeatCapacity(Cp) +julia> fig,ax,T,Cp_vec = PlotHeatCapacity(Cp) ``` -you can now save the figure to disk with: -``` -julia> using Plots -julia> savefig(plt,"Tdependent_heatcapacity.png") -``` - """ -function PlotHeatCapacity(Cp::AbstractHeatCapacity; T=nothing, plt=nothing, lbl=nothing) +function PlotHeatCapacity(x; + args=(; ), + Stress=(1e0, 1e8), + linestyle=:solid, + linewidth=1, + color=nothing, + label=nothing, + title="", + fig=nothing, + filename=nothing, + res=(1200, 1200), + legendsize=15, + labelsize=35, + T=nothing) + if isnothing(T) T = collect(273.0:10:1250) * K end + n = 1 + if isa(x, Tuple) + n = length(x) + end + + if isnothing(fig) + fig = Figure(; fontsize=25, size=res) + end + ax = Axis( + fig[1, 1]; + xlabel="Temperature [K]", + ylabel="Heat Capacity [J kg⁻¹·⁰ K⁻¹·⁰]", + xlabelsize=labelsize, + ylabelsize=labelsize, + title=title, + ) + args = (; T=ustrip.(T)) Cp1 = zeros(size(T)) - compute_heatcapacity!(Cp1, Cp, args) - if length(Cp) == 1 - Cp1 = ones(size(T)) * Cp1 + #compute_heatcapacity!(Cp1, Cp, args) + #if length(Cp) == 1 + # Cp1 = ones(size(T)) * Cp1 + #end + + #if isnothing(plt) + # plt = plot(ustrip(T), ustrip(Cp); label=lbl) + #else + # plt = plot!(ustrip(T), ustrip(Cp); label=lbl) + #end + #lines(plt; xlabel="Temperature [$(unit(T[1]))]", ylabel="Cp [$(unit(Cp[1]))]") + #gui(plt) + + for i in 1:n + if isa(x, Tuple) + p = x[i] + else + p = x + end + + if isa(args, Tuple) + args_in = args[i] + else + args_in = args + end + + compute_heatcapacity!(Cp1, p, args) + # Retrieve plot arguments (label, color etc.) + plot_args = ObtainPlotArgs(i, p, args_in, linewidth, linestyle, color, label) + + # Create plot: + li = lines!(ustrip.(T), Cp1) # plot line + + # Customize plot: + customize_plot!(li, plot_args) end - if isnothing(plt) - plt = plot(ustrip(T), ustrip(Cp); label=lbl) + #axislegend(ax; labelsize=legendsize) + + if !isnothing(filename) + save(filename, fig) else - plt = plot!(ustrip(T), ustrip(Cp); label=lbl) + display(fig) end - plot!(plt; xlabel="Temperature [$(unit(T[1]))]", ylabel="Cp [$(unit(Cp[1]))]") - gui(plt) - return T, Cp1, plt + return fig, ax, T, Cp1 end """ From 7572cfc8029fd68c6749023f8ef83b64646dd402 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Sat, 12 Oct 2024 17:27:12 +0200 Subject: [PATCH 24/33] more fixes --- src/Plotting/Plotting.jl | 22 +++++++++++++++------- 1 file changed, 15 insertions(+), 7 deletions(-) diff --git a/src/Plotting/Plotting.jl b/src/Plotting/Plotting.jl index d621642b..e8d92bc4 100644 --- a/src/Plotting/Plotting.jl +++ b/src/Plotting/Plotting.jl @@ -633,6 +633,8 @@ function PlotHeatCapacity(x; return fig, ax, T, Cp1 end + +# TO BE FIXED """ T,Kk,plt = PlotConductivity(Cp::AbstractConductivity; T=nothing, plt=nothing, lbl=nothing) @@ -688,6 +690,7 @@ function PlotConductivity( return T, Cond, plt end +# TO BE FIXED """ T,phi,plt = PlotMeltFraction(p::AbstractMeltingParam; T=nothing, plt=nothing, lbl=nothing) @@ -747,6 +750,7 @@ function PlotMeltFraction( return T, phi, dϕdT end +# BROKEN """ plt, data, Tvec, Pvec = PlotPhaseDiagram(p::PhaseDiagram_LookupTable; fieldname::Symbol, Tvec=nothing, Pvec=nothing) @@ -757,7 +761,7 @@ The return arguments are the plotting object `plt` (so you can modify properties Example ======= ```julia -julia> PD_data = Read_LaMEM_Perple_X_Diagram("Peridotite.in") +julia> PD_Data = PerpleX_LaMEM_Diagram("./test/test_data/Peridotite.in") Perple_X/LaMEM Phase Diagram Lookup Table: File : Peridotite.in T : 293.0 - 1573.000039 @@ -831,16 +835,16 @@ end """ - plt = Plot_TAS_diagram() + plt = Plot_TAS_diagram(; displayLabel=true,size=(1500,1500), fontsize=18) Creates a TAS diagram plot """ -function Plot_TAS_diagram(; displayLabel=true) +function Plot_TAS_diagram(; displayLabel=true,sz=(1500,1500), fontsz=18) # get TAS diagram data from TASclassification routine ClassTASdata = TASclassificationData() @unpack litho, n_ver, ver = ClassTASdata - f = Figure(size = (1100, 1100), fontsize = 18) + f = Figure(size = sz, fontsize = fontsz) p1 = GridLayout(f[1, 1]) ax1 = Axis( p1[1, 1], @@ -853,6 +857,7 @@ function Plot_TAS_diagram(; displayLabel=true) ) n_poly = size(litho, 2) shift = 1 + #= for poly in 1:n_poly shift_poly = shift:(shift + n_ver[poly] - 1) x = sum(ver[i, 1] for i in shift_poly) / n_ver[poly] @@ -888,6 +893,7 @@ function Plot_TAS_diagram(; displayLabel=true) rowsize!(p2, 1, 700) colsize!(p2, 1, 225) end + =# display(f) return f @@ -1001,9 +1007,10 @@ Creates a deformation mechanism map (T/εII vs. stress/viscosity or T/τII vs. s # Example ```julia -julia> v1 = SetDiffusionCreep("Dry Anorthite | Rybacki et al. (2006)") -julia> v2 = SetDislocationCreep("Dry Anorthite | Rybacki et al. (2006)") -julia> v=CompositeRheology(v1,v2) +julia> import GeoParams.Diffusion, GeoParams.Dislocation +julia> v1 = SetDiffusionCreep(Diffusion.dry_anorthite_Rybacki_2006); +julia> v2 = SetDislocationCreep(Dislocation.dry_anorthite_Rybacki_2006); +julia> v = CompositeRheology(v1,v2) julia> PlotDeformationMap(v, levels=100, colormap=:roma) ``` Next, let's plot viscosity and flip x & y axis: @@ -1055,6 +1062,7 @@ function PlotDeformationMap( n_components = length(v) end + d_vec = range(d[1], d[2], n+1) if strainrate # compute ε as a function of τ and T From 271dc8110f7cfe1a7abc40e429434f30689275af Mon Sep 17 00:00:00 2001 From: Pascal Aellig <93140290+aelligp@users.noreply.github.com> Date: Thu, 17 Oct 2024 11:16:17 +0200 Subject: [PATCH 25/33] Update blank.yml --- .github/workflows/blank.yml | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/.github/workflows/blank.yml b/.github/workflows/blank.yml index 51e36e74..27466de9 100644 --- a/.github/workflows/blank.yml +++ b/.github/workflows/blank.yml @@ -26,11 +26,23 @@ jobs: - 'nightly' os: - ubuntu-latest - - macOS-latest - macOS-13 - windows-latest arch: - x64 + include: + - os: macOS-latest + arch: aarch64 + version: '1.7' + - os: macOS-latest + arch: aarch64 + version: '1.10' + - os: macOS-latest + arch: aarch64 + version: '1.11' + - os: macOS-latest + arch: aarch64 + version: 'nightly' steps: - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@v2 From b25c0902a07c392486bc66e2ae3c99a2cf2c5b0b Mon Sep 17 00:00:00 2001 From: Pascal Aellig <93140290+aelligp@users.noreply.github.com> Date: Thu, 17 Oct 2024 11:38:00 +0200 Subject: [PATCH 26/33] Move to 1.9 as lower bound for both arch --- .github/workflows/blank.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/blank.yml b/.github/workflows/blank.yml index 27466de9..2bfb7a70 100644 --- a/.github/workflows/blank.yml +++ b/.github/workflows/blank.yml @@ -20,7 +20,7 @@ jobs: fail-fast: false matrix: version: - - '1.7' + - '1.9' - '1.10' - '1.11' - 'nightly' @@ -33,7 +33,7 @@ jobs: include: - os: macOS-latest arch: aarch64 - version: '1.7' + version: '1.9' - os: macOS-latest arch: aarch64 version: '1.10' From e8d5668a7f72cb0e4c1669dcc44747a9aadf62a1 Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Thu, 17 Oct 2024 17:28:52 +0200 Subject: [PATCH 27/33] update to working tests --- src/GeoParams.jl | 1 + src/Permeability/Permeability.jl | 84 +++++--------------------------- test/test_Permeability.jl | 52 +++++++++++++++----- 3 files changed, 52 insertions(+), 85 deletions(-) diff --git a/src/GeoParams.jl b/src/GeoParams.jl index 116674b8..eb9fda4e 100644 --- a/src/GeoParams.jl +++ b/src/GeoParams.jl @@ -339,6 +339,7 @@ export compute_meltfraction, using .MaterialParameters.Permeability export compute_permeability, compute_permeability!, + compute_permeability_ratio, param_info, AbstractPermeability, ConstantPermeability, diff --git a/src/Permeability/Permeability.jl b/src/Permeability/Permeability.jl index 0d565266..878ebe30 100644 --- a/src/Permeability/Permeability.jl +++ b/src/Permeability/Permeability.jl @@ -15,14 +15,13 @@ abstract type AbstractPermeability{T} <: AbstractMaterialParam end export compute_permeability, # calculation routines compute_permeability!, # in place calculation + compute_permeability_ratio, # calculation with phase ratios param_info, # info about the parameters AbstractPermeability, ConstantPermeability, # constant HazenPermeability, # Hazen equation PowerLawPermeability, # Power-law permeability - CarmanKozenyPermeability, # Carman-Kozeny permeability - BiotWillis, # Biot-Willis coefficient - SkemptonCoeff # Skempton coefficient + CarmanKozenyPermeability # Carman-Kozeny permeability # Define "empty" computational routines in case nothing is defined function compute_permeability!( @@ -60,6 +59,8 @@ end ConstantPermeability(args...) = ConstantPermeability(convert.(GeoUnit, args)...) isdimensional(s::ConstantPermeability) = isdimensional(s.k) +@inline (s::ConstantPermeability)(; args...) = s.k.val +@inline (s::ConstantPermeability)(args) = s(; args...) @inline compute_permeability(s::ConstantPermeability{_T}, args) where {_T} = s(; args...) @inline compute_permeability(s::ConstantPermeability{_T}) where {_T} = s() @@ -104,7 +105,7 @@ function (s::HazenPermeability{_T})(; kwargs...) where {_T} end @inline (s::HazenPermeability)(args) = s(; args...) -@inline compute_permeability(s::HazenPermeability, args) = s(; args...) +@inline compute_permeability(s::HazenPermeability, args) = s(args) function show(io::IO, g::HazenPermeability) @@ -125,7 +126,7 @@ function param_info(s::PowerLawPermeability) return MaterialParamsInfo(; Equation = L"k = c \cdot k_0 \cdot \phi^n") end -function (s::PowerLawPermeability{_T})(; ϕ=0e0, kwargs...) where {_T} +function (s::PowerLawPermeability{_T})(; ϕ=1e-2, kwargs...) where {_T} if ϕ isa Quantity @unpack_units c, k0, n = s else @@ -172,75 +173,12 @@ function show(io::IO, g::CarmanKozenyPermeability) return print(io, "Carman-Kozeny permeability: k = c * (ϕ / ϕ0)^n; c=$(g.c); ϕ0=$(g.ϕ0); n=$(g.n)") end -# ----------------------------------------------- -# This implements the methods described by Yarushina and Podladchikov, 2015 -""" -Biot-Willis coefficient -""" -@with_kw_noshow struct BiotWillis{_T,U} <: AbstractPermeability{_T} - Kd::GeoUnit{_T,U} = 1.0 * Pa # Drained bulk modulus - Ks::GeoUnit{_T,U} = 1.0 * Pa # Solid grain bulk modulus - α::GeoUnit{_T,U} = 1.0 * NoUnits # Biot-Willis coefficient to keep track of dimensionality -end - -BiotWillis(args...) = BiotWillis(convert.(GeoUnit, args)...) -isdimensional(s::BiotWillis) = isdimensional(s.Kd) - -function param_info(s::BiotWillis) - return MaterialParamsInfo(; Equation = L"\alpha = 1 - \frac{K_d}{K_s}") -end - -function (bw::BiotWillis{_T})(; kwargs...) where {_T} - if kwargs isa Quantity - @unpack_units Kd, Ks, α = bw - else - @unpack_val Kd, Ks, α = bw - end - - return 1 - (Kd * inv(Ks)) -end - -@inline (s::BiotWillis)(args) = s(; args...) -@inline compute_biot_willis(s::BiotWillis, args) = s(args) - -function show(io::IO, g::BiotWillis) - return print(io, "Biot-Willis coefficient: α = 1 - (Kd / Ks); Kd=$(g.Kd); Ks=$(g.Ks)") -end - -""" -Skempton coefficient -""" -@with_kw_noshow struct SkemptonCoeff{_T,U} <: AbstractPermeability{_T} - Kd::GeoUnit{_T,U} = 1.0 * Pa # Drained bulk modulus - Ks::GeoUnit{_T,U} = 1.0 * Pa # Solid grain bulk modulus - Kf::GeoUnit{_T,U} = 1.0 * Pa # Fluid bulk modulus - ϕ::GeoUnit{_T,U} = 1.0 * NoUnits # Porosity - B::GeoUnit{_T,U} = 1.0 * Pa # Skempton coefficient to keep track of dimensionality -end - -SkemptonCoeff(args...) = SkemptonCoeff(convert.(GeoUnit, args)...) -isdimensional(s::SkemptonCoeff) = isdimensional(s.Kd) - -function param_info(s::SkemptonCoeff) - return MaterialParamsInfo(; Equation = L"B = \frac{(1/K_d - 1/K_s)}{(1/K_d - 1/K_s) + \phi (1/K_f - 1/K_s)}") -end - -function (sc::SkemptonCoeff{_T})(; kwargs...) where {_T} - if kwargs isa Quantity - @unpack_units Kd, Ks, Kf, ϕ, B = sc - else - @unpack_val Kd, Ks, Kf, ϕ, B = sc - end - - return (inv(Kd) - inv(Ks)) * inv(inv(Kd) - inv(Ks) + ϕ * (inv(Kf) - inv(Ks))) -end - -@inline (s::SkemptonCoeff)(args) = s(; args...) -@inline compute_skempton_coeff(s::SkemptonCoeff, args) = s(args) - -function show(io::IO, g::SkemptonCoeff) - return print(io, "Skempton Coefficient: B = (1/Kd - 1/Ks) / ((1/Kd - 1/Ks) + ϕ * (1/Kf - 1/Ks)); Kd=$(g.Kd); Ks=$(g.Ks); Kf=$(g.Kf); ϕ=$(g.ϕ)") +#------------------------------------------------------------------------------------------------------------------# +# Computational routines needed for computations with the MaterialParams structure +function compute_permeability(s::AbstractMaterialParamsStruct, args) + return compute_permeability(s.Permeability[1], args) end +#------------------------------------------------------------------------------------------------------------- """ compute_permeability!(k::AbstractArray{_T, N}, MatParam::NTuple{K,AbstractMaterialParamsStruct}, PhaseRatios::AbstractArray{_T, M}, P=nothing, T=nothing) diff --git a/test/test_Permeability.jl b/test/test_Permeability.jl index 7704389e..bca04325 100644 --- a/test/test_Permeability.jl +++ b/test/test_Permeability.jl @@ -31,6 +31,7 @@ using Test, GeoParams, Unitful, StaticArrays, LaTeXStrings x1 = ConstantPermeability(; k=1e-12m^2) @test x1.k.val == 1e-12 @test GeoParams.get_k(x1) == 1e-12 + @test compute_permeability(x1) ≈ 1e-12 x2 = HazenPermeability(; C=1.0, D10=1e-4m) args = (;ϕ=0.4) @@ -62,18 +63,45 @@ using Test, GeoParams, Unitful, StaticArrays, LaTeXStrings x4 = nondimensionalize(x4, CharUnits_GEO) @test x4.c.val ≈ 1.0 - # # Test allocations - # k = [0.0] - # ϕ = 0.4 - # args = (ϕ=ϕ) - # # Test allocations using k alias - # k!(k, x3, args) - # num_alloc = @allocated k!(k, x3, args) - # @test num_alloc == 0 - - # k!(k, x4, args) - # num_alloc = @allocated k!(k, x4, args) - # @test num_alloc == 0 + # Define material parameters with different permeability parameterizations + Mat_tup = ( + SetMaterialParams(; + Name="Mantle", Phase=1, Permeability=ConstantPermeability(; k=1e-12m^2) + ), + SetMaterialParams(; Name="Crust", Phase=2, Permeability=HazenPermeability(; C=1.0, D10=1e-4m)), + SetMaterialParams(; + Name="UpperCrust", + Phase=3, + Permeability=PowerLawPermeability(; c=1.0, k0=1e-12m^2, n=3, ϕ=0.1), + Density=PT_Density(), + ), + SetMaterialParams(; Name="LowerCrust", Phase=4, Permeability=CarmanKozenyPermeability(; c=1.0, ϕ0=0.3, n=3), Density=PT_Density()), + ) + + n = 100 + Phases = ones(Int64, n, n, n) + Phases[:, :, 20:end] .= 2 + Phases[:, :, 50:end] .= 3 + Phases[:, :, 70:end] .= 4 + + ϕ = fill(1e-2, size(Phases)) + T = ones(size(Phases)) * (800 + 273.15) + P = ones(size(Phases)) * 10 + args = (P=P, T=T) + + @test compute_permeability(Mat_tup, Phases[1], args) == 1e-12 + + compute_permeability!(ϕ, Mat_tup, Phases, args) + + @test sum(ϕ) / n^3 ≈ 1.1484481671481687e-5 # Adjust this value based on expected results + + # Test PhaseRatio and StaticArrays PhaseRatios as input + args = (P=0.0, T=1000.0 + 273.15) + PhaseRatio = (0.25, 0.25, 0.25, 0.25) + @test compute_permeability_ratio(PhaseRatio, Mat_tup, args) == 9.26175950925951e-6 # Adjust this value based on expected results + + SvPhaseRatio = SA[0.25,0.25,0.25,0.25] + @test compute_permeability_ratio(SvPhaseRatio, Mat_tup, args) == 9.26175950925951e-6 # Adjust this value based on expected results end From eb2b85aa8f21b100a3b0be04b0a722689de5ae0a Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Thu, 17 Oct 2024 17:42:13 +0200 Subject: [PATCH 28/33] delete obsolete code --- src/GeoParams.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/GeoParams.jl b/src/GeoParams.jl index eb9fda4e..98f33e7d 100644 --- a/src/GeoParams.jl +++ b/src/GeoParams.jl @@ -345,9 +345,7 @@ export compute_permeability, ConstantPermeability, HazenPermeability, PowerLawPermeability, - CarmanKozenyPermeability, - BiotWillis, - SkemptonCoeff + CarmanKozenyPermeability include("Traits/rheology.jl") export RheologyTrait From 58bb9e7514335690897d0bf72098747721be6e85 Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Thu, 17 Oct 2024 17:58:44 +0200 Subject: [PATCH 29/33] rm Unitful --- test/test_Permeability.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_Permeability.jl b/test/test_Permeability.jl index bca04325..316cd6ee 100644 --- a/test/test_Permeability.jl +++ b/test/test_Permeability.jl @@ -1,4 +1,4 @@ -using Test, GeoParams, Unitful, StaticArrays, LaTeXStrings +using Test, GeoParams, StaticArrays, LaTeXStrings @testset "Permeability.jl" begin From 43ad526c5bf555195d922b12821a504f4f7a012d Mon Sep 17 00:00:00 2001 From: Pascal Aellig <93140290+aelligp@users.noreply.github.com> Date: Fri, 18 Oct 2024 16:18:02 +0200 Subject: [PATCH 30/33] Apply suggestions pt 1 Co-authored-by: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> --- src/Permeability/Permeability.jl | 20 ++++++++++---------- test/test_Permeability.jl | 4 ++-- 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/Permeability/Permeability.jl b/src/Permeability/Permeability.jl index 878ebe30..ef746746 100644 --- a/src/Permeability/Permeability.jl +++ b/src/Permeability/Permeability.jl @@ -84,8 +84,8 @@ end # Hazen Permeability @with_kw_noshow struct HazenPermeability{_T,U1,U2} <: AbstractPermeability{_T} - C::GeoUnit{_T,U1} = 1.0 * NoUnits # Hazen constant - D10::GeoUnit{_T,U2} = 1e-4 * m # Effective grain size + C::GeoUnit{_T,U1} = 1.0 * NoUnits # Hazen constant + D10::GeoUnit{_T,U2} = 1e-4 * m # Effective grain size end HazenPermeability(args...) = HazenPermeability(convert.(GeoUnit, args)...) isdimensional(s::HazenPermeability) = isdimensional(s.D10) @@ -101,7 +101,7 @@ function (s::HazenPermeability{_T})(; kwargs...) where {_T} @unpack_val C, D10 = s end - return @pow C * D10^2 + return C * D10^2 end @inline (s::HazenPermeability)(args) = s(; args...) @@ -114,10 +114,10 @@ end # Power-law Permeability @with_kw_noshow struct PowerLawPermeability{_T,U1,U2,U3, U4} <: AbstractPermeability{_T} - c::GeoUnit{_T,U1} = 1.0 * NoUnits # Power-law constant - k0::GeoUnit{_T,U2} = 1e-12 * m^2 # reference permeability - ϕ::GeoUnit{_T,U3} = 1e-2 * NoUnits # reference porosity - n::GeoUnit{_T,U4} = 3 * NoUnits # exponent + c::GeoUnit{_T,U1} = 1.0 * NoUnits # Power-law constant + k0::GeoUnit{_T,U2} = 1e-12 * m^2 # reference permeability + ϕ::GeoUnit{_T,U3} = 1e-2 * NoUnits # reference porosity + n::GeoUnit{_T,U4} = 3 * NoUnits # exponent end PowerLawPermeability(args...) = PowerLawPermeability(convert.(GeoUnit, args)...) isdimensional(s::PowerLawPermeability) = isdimensional(s.k0) @@ -146,8 +146,8 @@ end # Carman-Kozeny Permeability @with_kw_noshow struct CarmanKozenyPermeability{_T,U1,U2,U3} <: AbstractPermeability{_T} c::GeoUnit{_T,U1} = 1.0 * m^2 # Carman-Kozeny constant - ϕ0::GeoUnit{_T,U2} = 0.01 * NoUnits # reference porosity - n::GeoUnit{_T,U3} = 3 * NoUnits # exponent + ϕ0::GeoUnit{_T,U2} = 0.01 * NoUnits # reference porosity + n::GeoUnit{_T,U3} = 3 * NoUnits # exponent end CarmanKozenyPermeability(args...) = CarmanKozenyPermeability(convert.(GeoUnit, args)...) # isdimensional(s::CarmanKozenyPermeability) = isdimensional(s.c) @@ -166,7 +166,7 @@ function (s::CarmanKozenyPermeability{_T})(; ϕ=1e-2, kwargs...) where {_T} return @pow c * (ϕ / ϕ0)^n end -@inline (s::CarmanKozenyPermeability)(args) = s(; args...) +@inline (s::CarmanKozenyPermeability)(args) = s(; args...) @inline compute_permeability(s::CarmanKozenyPermeability, args) = s(args) function show(io::IO, g::CarmanKozenyPermeability) diff --git a/test/test_Permeability.jl b/test/test_Permeability.jl index 316cd6ee..ca27c20a 100644 --- a/test/test_Permeability.jl +++ b/test/test_Permeability.jl @@ -86,8 +86,8 @@ using Test, GeoParams, StaticArrays, LaTeXStrings Phases[:, :, 70:end] .= 4 ϕ = fill(1e-2, size(Phases)) - T = ones(size(Phases)) * (800 + 273.15) - P = ones(size(Phases)) * 10 + T = fill((800 + 273.15), size(Phases)) + P = fill(10, size(Phases)) args = (P=P, T=T) @test compute_permeability(Mat_tup, Phases[1], args) == 1e-12 From 6d679e20cef9950d17796cd87a56bfb98840e72a Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Fri, 18 Oct 2024 16:20:49 +0200 Subject: [PATCH 31/33] address suggestions pt 2 --- src/Permeability/Permeability.jl | 4 ++-- test/test_Permeability.jl | 5 ----- 2 files changed, 2 insertions(+), 7 deletions(-) diff --git a/src/Permeability/Permeability.jl b/src/Permeability/Permeability.jl index ef746746..53c45db2 100644 --- a/src/Permeability/Permeability.jl +++ b/src/Permeability/Permeability.jl @@ -4,7 +4,7 @@ using Parameters, Unitful, LaTeXStrings, MuladdMacro using ..Units using ..PhaseDiagrams using GeoParams: AbstractMaterialParam, AbstractMaterialParamsStruct, @extractors, add_extractor_functions -using GeoParams: fastpow, pow_check, @pow +using GeoParams: fastpow, pow_check, @pow, @fastpow import ..Units: isdimensional using ..MaterialParameters: No_MaterialParam, MaterialParamsInfo import Base.show, GeoParams.param_info @@ -133,7 +133,7 @@ function (s::PowerLawPermeability{_T})(; ϕ=1e-2, kwargs...) where {_T} @unpack_val c, k0, n = s end - return @pow c * k0 * ϕ^n + return @fastpow c * k0 * ϕ^n end @inline (s::PowerLawPermeability)(args) = s(; args...) diff --git a/test/test_Permeability.jl b/test/test_Permeability.jl index ca27c20a..068eb75a 100644 --- a/test/test_Permeability.jl +++ b/test/test_Permeability.jl @@ -2,11 +2,6 @@ using Test, GeoParams, StaticArrays, LaTeXStrings @testset "Permeability.jl" begin - # Set alias for permeability function - if !isdefined(Main, :GeoParamsAliases) - eval(:(@use GeoParamsAliases permeability = k)) - end - # Make sure that structs are isbits x = ConstantPermeability() @test isbits(x) From 035fe071e7e5dd129e664b5812caa70064db5662 Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Fri, 18 Oct 2024 16:33:09 +0200 Subject: [PATCH 32/33] correct syntax --- src/Permeability/Permeability.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Permeability/Permeability.jl b/src/Permeability/Permeability.jl index 53c45db2..36926553 100644 --- a/src/Permeability/Permeability.jl +++ b/src/Permeability/Permeability.jl @@ -4,7 +4,7 @@ using Parameters, Unitful, LaTeXStrings, MuladdMacro using ..Units using ..PhaseDiagrams using GeoParams: AbstractMaterialParam, AbstractMaterialParamsStruct, @extractors, add_extractor_functions -using GeoParams: fastpow, pow_check, @pow, @fastpow +using GeoParams: fastpow, pow_check, @pow import ..Units: isdimensional using ..MaterialParameters: No_MaterialParam, MaterialParamsInfo import Base.show, GeoParams.param_info @@ -133,7 +133,7 @@ function (s::PowerLawPermeability{_T})(; ϕ=1e-2, kwargs...) where {_T} @unpack_val c, k0, n = s end - return @fastpow c * k0 * ϕ^n + return c * k0 * fastpow(ϕ,n) end @inline (s::PowerLawPermeability)(args) = s(; args...) From ebfe9cfa6e71c125751314a684c8a4f67a86eb81 Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Fri, 18 Oct 2024 16:52:01 +0200 Subject: [PATCH 33/33] delete bulkviscosity file that slipped in --- src/Viscosity/BulkViscosity.jl | 20 -------------------- 1 file changed, 20 deletions(-) delete mode 100644 src/Viscosity/BulkViscosity.jl diff --git a/src/Viscosity/BulkViscosity.jl b/src/Viscosity/BulkViscosity.jl deleted file mode 100644 index 7774c389..00000000 --- a/src/Viscosity/BulkViscosity.jl +++ /dev/null @@ -1,20 +0,0 @@ -""" -bulk_viscosity(η_shear, ϕ, C, R, λ, P_eff) - -# η_shear: shear viscosity -# ϕ: porosity -# C: shear viscosity ratio (value, dimensionless) -# R: compaction strength ratio (value, dimensionless) -# λ: effective pressure transition zone -# P_eff: effective Pressure -""" -function bulk_viscosity(η_shear, ϕ, C, R, λ, P_eff) - - η_bulk = η_shear / ϕ * C * (1.0 + 0.5 * (inv(R) - 1.0) * (1.0 + tanh(-P_eff * inv(λ)))) - - return η_bulk -end - - -### Georg Adjoint 2 pahse paper, -### Ludo's 2019 paper has it quite nicely also with a calc for R