Vector density, melt fraction and heat capacity

src/Density/Density.jl

@@ -26,6 +26,7 @@
     PT_Density,             # P & T dependent density
     Compressible_Density,   # Compressible density
     T_Density,              # T dependent density
+    Vector_Density,         # Vector with density
     MeltDependent_Density   # Melt dependent density
 
 # Define "empty" computational routines in case nothing is defined
@@ -263,6 +264,41 @@
 end
 #-------------------------------------------------------------------------
 
+# MAGEMin DB density -------------------------------
+"""
+    Vector_Density(_T)
+
+Stores a vector with density data that can be retrieved by providing an `index`
+"""
+struct Vector_Density{_T, V <: AbstractVector} <: AbstractDensity{_T}
+    rho::V       # Density
+end
+Vector_Density(; rho=Vector{Float64}()) = Vector_Density{eltype(rho), typeof(rho)}(rho)
+
+# This assumes that density always has a single parameter. If that is not the case, we will have to extend this (to be done)
+function param_info(s::Vector_Density) # info about the struct
+    return MaterialParamsInfo(; Equation=L"\rho from a precomputed vector")
+end
+
+# Calculation routine
+"""
+    compute_density(s::Vector_Density; index::Int64, kwargs...)
+
+Pointwise calculation of density from a vector where `index` is the index of the point
+"""
+@inline function (s::Vector_Density)(; index::Int64, kwargs...)
+    return s.rho[index]
+end
+
+@inline (s::Vector_Density)(args)                = s(; args...)
+@inline compute_density(s::Vector_Density, args) = s(args)
+
+# Print info
+function show(io::IO, g::Vector_Density)
+    return print(io, "Density from precomputed vector with $(length(g.rho)) entries.")
+end
+#-------------------------------------------------------------------------
+
 #-------------------------------------------------------------------------
 # Phase diagrams
 function param_info(s::PhaseDiagram_LookupTable) # info about the struct
@@ -363,7 +399,7 @@
 @inline compute_density_ratio(args::Vararg{Any, N}) where N = compute_param_times_frac(compute_density, args...)
 
 # extractor methods
-for type in (ConstantDensity, PT_Density, Compressible_Density, T_Density, MeltDependent_Density)
+for type in (ConstantDensity, PT_Density, Compressible_Density, T_Density, MeltDependent_Density, Vector_Density)
     @extractors(type, :Density)
 end
 

src/Energy/HeatCapacity.jl

@@ -19,6 +19,7 @@
     ConstantHeatCapacity,                # constant
     T_HeatCapacity_Whittington,          # T-dependent heat capacity
     Latent_HeatCapacity,                 # Implement latent heat by modifying heat capacity
+    Vector_HeatCapacity,                 # Vector 
     param_info
 
 include("../Computations.jl")
@@ -163,6 +164,37 @@
 end
 #-------------------------------------------------------------------------
 
+#-------------------------------------------------------------------------
+"""
+    Vector_HeatCapacity(_T)
+
+Stores a vector with heat capacity data that can be retrieved by providing an `index`
+"""
+struct Vector_HeatCapacity{_T, V <: AbstractVector} <: AbstractHeatCapacity{_T}
+    Cp::V       # Heat capacity
+end
+Vector_HeatCapacity(; Cp=Vector{Float64}()) = Vector_HeatCapacity{eltype(Cp), typeof(Cp)}(Cp)
+
+function param_info(s::Vector_HeatCapacity) # info about the struct
+    return MaterialParamsInfo(; Equation=L"Cp from a precomputed vector")
+end
+
+"""
+    compute_heatcapacity(a::Vector_HeatCapacity; index::Int64, kwargs...)
+
+Pointwise calculation of heat capacity from a vector where `index` is the index of the point
+"""
+@inline function compute_heatcapacity(s::Vector_HeatCapacity; index::Int64=1, kwargs...)
+    return s.Cp[index]
+end
+
+# Print info
+function show(io::IO, g::Vector_HeatCapacity)
+    return print(io, "Heat capacity from precomputed vector with $(length(g.Cp)) entries")
+end
+
+#-------------------------------------------------------------------------
+
 
 #-------------------------------------------------------------------------
 """
@@ -173,7 +205,7 @@
 # function compute_heatcapacity!(Cp_array::AbstractArray{_T, N},s::T_HeatCapacity_Whittington{_T}; T::AbstractArray{_T, N}, kwargs...) where {_T,N} end
 
 # add methods programmatically
-for myType in (:ConstantHeatCapacity, :T_HeatCapacity_Whittington, :Latent_HeatCapacity)
+for myType in (:ConstantHeatCapacity, :T_HeatCapacity_Whittington, :Latent_HeatCapacity, :Vector_HeatCapacity)
     @eval begin
         #(s::$(myType))(args) = s(; args...)
         compute_heatcapacity(s::$(myType), args) = compute_heatcapacity(s; args...)
@@ -191,6 +223,7 @@
 
 
 
+
 #-------------------------------------------------------------------------
 # Heat capacity from phase diagram
 
@@ -267,7 +300,7 @@
 compute_heatcapacity!(args::Vararg{Any, N}) where N = compute_param!(compute_heatcapacity, args...)
 
 # extractor methods
-for type in (ConstantHeatCapacity, T_HeatCapacity_Whittington, Latent_HeatCapacity)
+for type in (ConstantHeatCapacity, T_HeatCapacity_Whittington, Latent_HeatCapacity, Vector_HeatCapacity)
     @extractors(type, :HeatCapacity)
 end
 

src/GeoParams.jl

@@ -127,6 +127,7 @@ export compute_density,                                # computational routines
     PT_Density,
     Compressible_Density,
     T_Density,
+    Vector_Density,
     PhaseDiagram_LookupTable,
     Read_LaMEM_Perple_X_Diagram,
     MeltDependent_Density,
@@ -275,7 +276,7 @@ export compute_gravity,                                # computational routines
 # Energy parameters: Heat Capacity, Thermal conductivity, latent heat, radioactive heat
 using .MaterialParameters.HeatCapacity
 export compute_heatcapacity,
-    compute_heatcapacity!, ConstantHeatCapacity, T_HeatCapacity_Whittington, Latent_HeatCapacity
+    compute_heatcapacity!, ConstantHeatCapacity, T_HeatCapacity_Whittington, Latent_HeatCapacity, Vector_HeatCapacity
 
 using .MaterialParameters.Conductivity
 export compute_conductivity,
@@ -336,6 +337,7 @@ export compute_meltfraction,
     MeltingParam_5thOrder,
     MeltingParam_Quadratic,
     MeltingParam_Assimilation,
+    Vector_MeltingParam,
     SmoothMelting
 
 include("CreepLaw/Data/DislocationCreep.jl")

src/MeltFraction/MeltingParameterization.jl

@@ -25,6 +25,7 @@
     MeltingParam_5thOrder,
     MeltingParam_Quadratic,
     MeltingParam_Assimilation,
+    Vector_MeltingParam,
     SmoothMelting
 
 include("../Utils.jl")
@@ -518,6 +519,37 @@
 
 #-------------------------------------------------------------------------
 
+# MAGEMin Database -------------------------------------------------------
+"""
+    Vector_MeltingParam(_T)
+
+Stores a vector with melt fraction that can be retrieved by providing an `index`
+"""
+struct Vector_MeltingParam{_T, V <: AbstractVector} <: AbstractMeltingParam{_T}
+    ϕ::V       # melt fraction
+end
+Vector_MeltingParam(; ϕ=Vector{Float64}()) = Vector_MeltingParam{eltype(ϕ), typeof(ϕ)}(ϕ)
+
+function param_info(s::Vector_MeltingParam) # info about the struct
+    return MaterialParamsInfo(; Equation=L"\phi from a precomputed vector")
+end
+
+# Calculation routines
+function (g::Vector_MeltingParam)(; index, kwargs...)
+    return g.ϕ[index]
+end
+
+function compute_dϕdT(g::Vector_MeltingParam; kwargs...)
+    return 0.0
+end
+
+# Print info
+function show(io::IO, g::Vector_MeltingParam)
+    return print(io, "Melt fraction from precomputed vector with $(length(g.ϕ)) entries")
+end
+#-------------------------------------------------------------------------
+
+
 # Smooth melting function ------------------------------------------------
 
 """
@@ -710,7 +742,8 @@
     :MeltingParam_4thOrder,
     :MeltingParam_Quadratic,
     :MeltingParam_Assimilation,
-    :SmoothMelting,
+    :Vector_MeltingParam,
+    :SmoothMelting, 
 )
     @eval begin
         (p::$(myType))(args) = p(; args...)

test/test_Density.jl

@@ -156,7 +156,7 @@ using Test, GeoParams, StaticArrays, LaTeXStrings
     @test ustrip(dimensionalize(Vs_ND, km / s, CharDim)) ≈ PD_data.Vs(1500.0, 1e8)
 
     # Test computation of density for the whole computational domain, using arrays
-    MatParam = Vector{AbstractMaterialParamsStruct}(undef, 4)
+    MatParam = Vector{AbstractMaterialParamsStruct}(undef, 5)
 
     MatParam[1] = SetMaterialParams(;
         Name="Mantle",
@@ -184,7 +184,12 @@ using Test, GeoParams, StaticArrays, LaTeXStrings
         CreepLaws=(PowerlawViscous(), LinearViscous(; η=1e23Pa * s)),
         Density=Compressible_Density(),
     )
-
+    MatParam[5] = SetMaterialParams(;
+        Name="LowerCrust",
+        Phase=4,
+        CreepLaws=(PowerlawViscous(), LinearViscous(; η=1e23Pa * s)),
+        Density=Vector_Density(rho=fill(2900.0,100)),
+    )
     Mat_tup = Tuple(MatParam)  # create a tuple to avoid allocations
 
     MatParam1 = Vector{AbstractMaterialParamsStruct}(undef, 4)
@@ -219,27 +224,28 @@ using Test, GeoParams, StaticArrays, LaTeXStrings
     @views Phases[:,  20:end] .= 1
     @views Phases[:, 200:end] .= 2
     @views Phases[:, 300:end] .= 3
+    @views Phases[:, 350:end] .= 4
 
     #Phases .= 2;
     rho = zeros(size(Phases))
     T = ones(size(Phases))
     P = fill(10.0, size(Phases))
 
-    args = (P=P, T=T)
+    args = (P=P, T=T, index=fill(10,size(T)))
 
     compute_density!(rho, MatParam, Phases, args)
 
     # Test computing density when Mat_tup1 is provided as a tuple
     compute_density!(rho, Mat_tup1, Phases, args)
     num_alloc = @allocated compute_density!(rho, Mat_tup1, Phases, args)   #      287.416 μs (0 allocations: 0 bytes)
-    @test sum(rho) / 400^2 ≈ 2945.000013499999
+    @test sum(rho) / 400^2 ≈ 2575.250013499998
     # @test num_alloc ≤ 32
 
     #Same test using function alias
     rho = zeros(size(Phases))
     ρ!(rho, Mat_tup1, Phases, args)
     num_alloc = @allocated compute_density!(rho, Mat_tup1, Phases, args)
-    @test sum(rho) / 400^2 ≈ 2945.000013499999
+    @test sum(rho) / 400^2 ≈ 2575.250013499998
     # @test num_alloc ≤ 32
 
     # Test for single phase
@@ -250,7 +256,7 @@ using Test, GeoParams, StaticArrays, LaTeXStrings
     @test sum(rho) / 400^2 ≈ 2895.5241895725003
 
     # test computing material properties when we have PhaseRatios, instead of Phase numbers
-    PhaseRatio = zeros(size(Phases)..., length(Mat_tup1))
+    PhaseRatio = zeros(size(Phases)..., length(Mat_tup))
     for i in CartesianIndices(Phases)
         iz = Phases[i]
         I = CartesianIndex(i, iz + 1)
@@ -260,19 +266,19 @@ using Test, GeoParams, StaticArrays, LaTeXStrings
     compute_density!(rho, Mat_tup1, PhaseRatio, args)
 
     num_alloc = @allocated compute_density!(rho, Mat_tup1, PhaseRatio, args) #   136.776 μs (0 allocations: 0 bytes)
-    @test sum(rho) / 400^2 ≈ 2945.000013499999
+    @test sum(rho) / 400^2 ≈ 2575.250013499998
     @test num_alloc == 0           # for some reason this does indicate allocations but @btime does not
 
     # Test calling the routine with only pressure as input.
     # This is ok for Mat_tup1, as it only has constant & P-dependent densities.
     # Note, however, that if you have P & T dependent densities and do this it will use 0 as default value for T
     compute_density!(rho, Mat_tup1, PhaseRatio, (; P=P))
-    @test sum(rho) / 400^2 ≈ 2945.000013499999
+    @test sum(rho) / 400^2 ≈ 2575.250013499998
 
     # In case we only want to compute with T, do this:
     #  NOTE that in this example the results are actually wrong (as some functions require P as well)
-    compute_density!(rho, Mat_tup, PhaseRatio, (P=zeros(size(T)), T=T))
-    @test sum(rho) / 400^2 ≈ 2895.5241749999996
+    compute_density!(rho, Mat_tup, PhaseRatio, (P=zeros(size(T)), T=T, index=fill(10,size(T)) ))
+    @test sum(rho) / 400^2 ≈ 2895.524175
 
     #Test computation of density given a single phase and P,T as scalars
     Phase, P, T = 0, 1.0, 1.0

test/test_Energy.jl

@@ -45,6 +45,15 @@ using StaticArrays
     # Dimensionalize again and double-check the results
     @test sum(abs.(ustrip.(Cp_nd * CharUnits_GEO.heatcapacity) - Cp)) < 1e-11
 
+    # heat capacity
+    Cp3 = Vector_HeatCapacity(Cp=fill(1500.0, 100))
+    index = fill(10,size(T))
+    args  = (; index=index)
+    Cp = zero(T)
+    @test compute_heatcapacity(Cp3, index=10) == 1500.0
+
+   # compute_heatcapacity!(Cp, Cp3, args)
+
     # Test with arrays
     T_array = T * ones(10)'
     Cp_array = similar(T_array)
@@ -61,7 +70,7 @@ using StaticArrays
     @test sum(Cp_array[i, 1] for i in axes(Cp_array,1)) ≈ 11667.035717418683
 
     # Check that it works if we give a phase array
-    MatParam = Array{MaterialParams,1}(undef, 2)
+    MatParam = Array{MaterialParams,1}(undef, 3)
     MatParam[1] = SetMaterialParams(;
         Name="Mantle", Phase=1, HeatCapacity=ConstantHeatCapacity()
     )
@@ -70,6 +79,10 @@ using StaticArrays
         Name="Crust", Phase=2, HeatCapacity=T_HeatCapacity_Whittington()
     )
 
+    MatParam[3] = SetMaterialParams(;
+        Name="Crust1", Phase=3, HeatCapacity=Vector_HeatCapacity(Cp=fill(1500.0, 100))
+    )
+
     Mat_tup = Tuple(MatParam)
 
     Mat_tup1 = (
@@ -83,22 +96,23 @@ using StaticArrays
     n = 100
     Phases = ones(Int64, n, n, n);
     @views Phases[:, :, 20:end] .= 2;
-
+    @views Phases[:, :, 50:end] .= 3;
+
     Cp = zeros(size(Phases));
     T = fill(1500e0, size(Phases));
     P = zeros(size(Phases));
 
-    args = (; T=T)
+    args = (; T=T, index=fill(10, size(T)))
     compute_heatcapacity!(Cp, Mat_tup, Phases, args)    # computation routine w/out P (not used in most heat capacity formulations)
-    @test sum(Cp[1, 1, k] for k in axes(Cp,3)) ≈ 121399.0486067196
+    @test sum(Cp[1, 1, k] for k in axes(Cp,3)) ≈ 134023.72170619245
 
     # check with array of constant properties (and no required input args)
     args1 = (;)
     compute_heatcapacity!(Cp, Mat_tup1, Phases, args1)    # computation routine w/out P (not used in most heat capacity formulations)
-    @test sum(Cp[1, 1, k] for k in axes(Cp,3)) ≈ 109050.0
+    @test sum(Cp[1, 1, k] for k in axes(Cp,3)) ≈ 52950.0
 
     num_alloc = @allocated compute_heatcapacity!(Cp, Mat_tup, Phases, args)
-    @test sum(Cp[1, 1, k] for k in axes(Cp,3)) ≈ 121399.0486067196
+    @test sum(Cp[1, 1, k] for k in axes(Cp,3)) ≈ 134023.72170619245
 
     @test num_alloc == 0
 
@@ -111,7 +125,7 @@ using StaticArrays
     end
     compute_heatcapacity!(Cp, Mat_tup, PhaseRatio, args)
     num_alloc = @allocated compute_heatcapacity!(Cp, Mat_tup, PhaseRatio, args)
-    @test sum(Cp[1, 1, k] for k in axes(Cp,3)) ≈ 121399.0486067196
+    @test sum(Cp[1, 1, k] for k in axes(Cp,3)) ≈ 134023.72170619245
     @test num_alloc == 0
 
 
@@ -130,7 +144,6 @@ using StaticArrays
     @test isdimensional(x_ND)==false
     @test isdimensional(x_ND1)==false
 
-
     dϕdT = 0.1
     dϕdT_ND = nondimensionalize(dϕdT / K, CharUnits_GEO)
     args = (; dϕdT=dϕdT, T=300.0+273)