Skip to content

Commit

Permalink
Merge pull request #225 from JuliaGeodynamics/pa-permea
Browse files Browse the repository at this point in the history 
Add Permeability
  • Loading branch information
@aelligp
aelligp authored Oct 18, 2024
2 parents 746ea96 + ebfe9cf commit 2b4f474
Show file tree
Hide file tree
Showing 5 changed files with 371 additions and 51 deletions.
24 changes: 18 additions & 6 deletions src/GeoParams.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -334,6 +334,18 @@ export compute_meltfraction,
Vector_MeltingParam,
SmoothMelting


using .MaterialParameters.Permeability
export compute_permeability,
compute_permeability!,
compute_permeability_ratio,
param_info,
AbstractPermeability,
ConstantPermeability,
HazenPermeability,
PowerLawPermeability,
CarmanKozenyPermeability

include("Traits/rheology.jl")
export RheologyTrait
export islinear, LinearRheologyTrait, NonLinearRheologyTrait
Expand Down Expand Up @@ -362,18 +374,18 @@ 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)
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
Expand Down Expand Up @@ -441,4 +453,4 @@ export get_G, get_Kb
const get_shearmodulus = get_G
const get_bulkmodulus = get_Kb

end # module GeoParams
end # module GeoParams
97 changes: 52 additions & 45 deletions src/MaterialParameters.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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")
Expand All @@ -58,7 +59,7 @@ using .ConstitutiveRelationships: print_rheology_matrix

"""
MaterialParams
Structure that holds all material parameters for a given phase
"""
Expand All @@ -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
Expand All @@ -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:
Expand All @@ -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
```
Expand All @@ -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:
Expand All @@ -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
```
Expand All @@ -195,6 +199,7 @@ function SetMaterialParams(;
RadioactiveHeat=nothing,
LatentHeat=nothing,
ShearHeat=nothing,
Permeability=nothing,
Melting=nothing,
SeismicVelocity=nothing,
CharDim=nothing,
Expand All @@ -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,
Expand All @@ -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,
Expand All @@ -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

Expand All @@ -274,23 +282,23 @@ 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
return field
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
Expand All @@ -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
Expand Down Expand Up @@ -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])
Expand All @@ -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)
Expand All @@ -358,4 +366,3 @@ function print_composite(a, spaces=10)
end

end

Loading

0 comments on commit 2b4f474

Please sign in to comment.