diff --git a/Project.toml b/Project.toml index bcf0ff196..9301a80e5 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "GeoParams" uuid = "e018b62d-d9de-4a26-8697-af89c310ae38" authors = ["Boris Kaus "] -version = "0.2.3" +version = "0.2.4" [deps] BibTeX = "7b0aa2c9-049f-5cec-894a-2b6b781bb25e" diff --git a/docs/src/man/melting.md b/docs/src/man/melting.md index fef6da182..8515010ff 100644 --- a/docs/src/man/melting.md +++ b/docs/src/man/melting.md @@ -16,7 +16,14 @@ GeoParams.MeltingParam.compute_meltfraction! GeoParams.MeltingParam.compute_meltfraction ``` -Also note that phase diagrams can be imported using `PerpleX_LaMEM_Diagram`. +You can also obtain the derivative of melt fraction versus temperature with (useful to compute latent heat effects): +```@docs +GeoParams.MeltingParam.compute_dϕdT! +GeoParams.MeltingParam.compute_dϕdT +``` + +Also note that phase diagrams can be imported using `PerpleX_LaMEM_Diagram`, which may also have melt content information. +The computational routines work with that as well. # Plotting routines You can use the routine `PlotMeltFraction` to create a plot, provided that the `Plots` package has been loaded diff --git a/src/Energy/Conductivity.jl b/src/Energy/Conductivity.jl index 0549cbcf9..8fb0c1ab4 100644 --- a/src/Energy/Conductivity.jl +++ b/src/Energy/Conductivity.jl @@ -87,16 +87,16 @@ Their parameterization is originally given for the thermal diffusivity, together Cp = a + b T - c/T^2 ``` ```math - \\kappa = d/T - e, if T<=846K + \\kappa = d/T - e \\textrm{ if } T<=846K ``` ```math - \\kappa = f - g*T, if T>846K + \\kappa = f - g*T \\textrm{ if } T>846K ``` ```math - \\rho = 2700 kg/m3 + \\rho = 2700 kg/m^3 ``` ```math - k = \\kappa Cp \\rho + k = \\kappa \\rho Cp ``` where ``Cp`` is the heat capacity [``J/mol/K``], and ``a,b,c`` are parameters that dependent on the temperature `T`: @@ -197,7 +197,7 @@ end Sets a temperature-dependent conductivity that is parameterization after *Whittington, et al. 2009* -The original parameterization involves quite a few parameters; this is a polynomial fit that is roughly valid from 0-1000C +The original parameterization involves quite a few parameters; this is a polynomial fit that is roughly valid from 0-1000Celcius ```math k [W/m/K] = -2 10^{-9} (T-Ts)^3 + 6 10^{-6} (T-Ts)^2 - 0.0062 (T-Ts) + 4 ``` @@ -206,8 +206,6 @@ The original parameterization involves quite a few parameters; this is a polynom ``` where `T[K]` is the temperature in Kelvin (or the nondimensional equivalent of it). """ - - @with_kw_noshow struct T_Conductivity_Whittington_parameterised{T,U1,U2,U3,U4,U5} <: AbstractConductivity{T} # Note: the resulting curve of k was visually compared with Fig. 2 of the paper a::GeoUnit{T,U1} = -2e-09Watt/m/K^4 # @@ -269,7 +267,6 @@ end #------------------------------------------------------------------------- - # Temperature (& Pressure) dependent conductivity ------------------------------- """ TP_Conductivity() diff --git a/src/GeoParams.jl b/src/GeoParams.jl index 67ce5ebc8..bb23521cc 100644 --- a/src/GeoParams.jl +++ b/src/GeoParams.jl @@ -111,16 +111,17 @@ export compute_zirconsaturation, compute_zirconsaturation!, # calculation # Seismic velocities using .MaterialParameters.SeismicVelocity -export compute_pwave_velocity, compute_swave_velocity, - compute_pwave_velocity!, compute_swave_velocity!, +export compute_pwave_velocity, compute_swave_velocity, + compute_pwave_velocity!, compute_swave_velocity!, ConstantSeismicVelocity # Add melting parameterizations include("./MeltFraction/MeltingParameterization.jl") using .MeltingParam -export compute_meltfraction, compute_meltfraction!, # calculation routines - MeltingParam_Caricchi, MeltingParam_4thOrder, - MeltingParam_5thOrder, MeltingParam_Quadratic +export compute_meltfraction, compute_meltfraction!, # calculation routines + compute_dϕdT, compute_dϕdT!, + MeltingParam_Caricchi, MeltingParam_4thOrder, + MeltingParam_5thOrder, MeltingParam_Quadratic # Add plotting routines - only activated if the "Plots.jl" package is loaded diff --git a/src/MeltFraction/MeltingParameterization.jl b/src/MeltFraction/MeltingParameterization.jl index 3268f8d03..ff934b1f4 100644 --- a/src/MeltFraction/MeltingParameterization.jl +++ b/src/MeltFraction/MeltingParameterization.jl @@ -13,6 +13,8 @@ abstract type AbstractMeltingParam{T} <: AbstractMaterialParam end export compute_meltfraction, compute_meltfraction!, # calculation routines + compute_dϕdT, # derivative of melt fraction versus + compute_dϕdT!, param_info, MeltingParam_Caricchi, MeltingParam_4thOrder, @@ -28,10 +30,10 @@ include("../Computations.jl") Implements the T-dependent melting parameterisation used by Caricchi et al ```math - \\theta = (800.0 .- (T + 273.15))./23.0 + \\theta = {(800.0 - (T + 273.15)) \\over 23.0} ``` ```math - \\phi_{solid} = 1.0 - {1.0 \\over (1.0 + e^\\theta)}; + \\phi_{melt} = {1.0 \\over (1.0 + e^\\theta)} ``` Note that T is in Kelvin. @@ -48,7 +50,7 @@ function param_info(s::MeltingParam_Caricchi) # info about the struct return MaterialParamsInfo(Equation = L"\phi = {1 \over {1 + \exp( {800-T[^oC] \over 23})}}") end -# Calculation routine +# Calculation routines function compute_meltfraction(p::MeltingParam_Caricchi{_T}, P::Quantity, T::Quantity) where _T @unpack_units a,b,c = p @@ -75,6 +77,31 @@ function compute_meltfraction!(ϕ::AbstractArray{_T}, p::MeltingParam_Caricchi{_ return nothing end +function compute_dϕdT(p::MeltingParam_Caricchi{_T}, P::Quantity, T::Quantity) where _T + @unpack_units a,b,c = p + + # analytically computed derivative (using Symbolics.jl) + dϕdT = exp((a + c - T) / b) / (b*((1.0 + exp((a + c - T) / b))^2)) + + return dϕdT +end + +function compute_dϕdT(p::MeltingParam_Caricchi{_T}, P::_T, T::_T ) where _T + @unpack_val a,b,c = p + + dϕdT = exp((a + c - T) / b) / (b*((1.0 + exp((a + c - T) / b))^2)) + + return dϕdT +end + +function compute_dϕdT!(dϕdT::AbstractArray{_T}, p::MeltingParam_Caricchi{_T}, P::AbstractArray{_T}, T::AbstractArray{_T}) where _T + @unpack_val a,b,c = p + + @. dϕdT = exp((a + c - T) / b) / (b*((1.0 + exp((a + c - T) / b))^2)) + + return nothing +end + # Print info function show(io::IO, g::MeltingParam_Caricchi) print(io, "Caricchi et al. melting parameterization") @@ -88,13 +115,13 @@ end Uses a 5th order polynomial to describe the melt fraction `phi` between solidus temperature `T_s` and liquidus temperature `T_l` ```math - \\phi = a T^5 + b T^4 + c T^3 + d T^2 + e T + f \\textrm{ if } T_s ≤ T ≤ T_l + \\phi = a T^5 + b T^4 + c T^3 + d T^2 + e T + f \\textrm{ for } T_s ≤ T ≤ T_l ``` ```math - \\phi = 1 \\textrm{ if } T>T_l + \\phi = 1 \\textrm{ if } T>T_l ``` ```math - \\phi = 0 \\textrm{ if } TT_l + dϕdT = 0. + end + + return dϕdT +end + +function compute_dϕdT(p::MeltingParam_5thOrder{_T}, P::_T, T::_T ) where _T + @unpack_val a,b,c,d,e,T_s,T_l = p + + dϕdT = 5*a*T^4 + 4*b*T^3 + 3*c*T^2 + 2*d*T + e + if TT_l + dϕdT = 0. + end + + return dϕdT +end + +function compute_dϕdT!(dϕdT::AbstractArray{_T}, p::MeltingParam_5thOrder{_T}, P::AbstractArray{_T}, T::AbstractArray{_T}) where _T + @unpack_val a,b,c,d,e,f,T_s,T_l = p + + @. dϕdT = 5*a*T^4 + 4*b*T^3 + 3*c*T^2 + 2*d*T + e + + dϕdT[T.T_l] .= 0. + + return nothing +end + # Print info function show(io::IO, g::MeltingParam_5thOrder) print(io, "5th order polynomial melting curve: phi = $(NumValue(g.a)) T^5 + $(NumValue(g.b))T^4 + $(NumValue(g.c))T^3 + $(NumValue(g.d))T^2 + $(NumValue(g.e))T + $(NumValue(g.f)) $(Value(g.T_s)) ≤ T ≤ $(Value(g.T_l))") @@ -171,13 +230,13 @@ end Uses a 4th order polynomial to describe the melt fraction `phi` between solidus temperature `T_s` and liquidus temperature `T_l` ```math - \\phi = b T^4 + c T^3 + d T^2 + e T + f \\textrm{ if } T_s ≤ T ≤ T_l + \\phi = b T^4 + c T^3 + d T^2 + e T + f \\textrm{ for } T_s ≤ T ≤ T_l ``` ```math - \\phi = 1 \\textrm{ if } T>T_l + \\phi = 1 \\textrm{ if } T>T_l ``` ```math - \\phi = 0 \\textrm{ if } TT_l + dϕdT = 0. + end + + return dϕdT +end + + +function compute_dϕdT(p::MeltingParam_4thOrder{_T}, P::_T, T::_T ) where _T + @unpack_val b,c,d,e,T_s,T_l = p + + dϕdT = 4*b*T^3 + 3*c*T^2 + 2*d*T + e + if TT_l + dϕdT = 0. + end + return dϕdT + +end + +function compute_dϕdT!(dϕdT::AbstractArray{_T}, p::MeltingParam_4thOrder{_T}, P::AbstractArray{_T}, T::AbstractArray{_T}) where _T + @unpack_val b,c,d,e,T_s,T_l = p + + @. dϕdT = 4*b*T^3 + 3*c*T^2 + 2*d*T + e + + dϕdT[T.T_l] .= 0. + + return nothing +end + + # Print info function show(io::IO, g::MeltingParam_4thOrder) print(io, "4th order polynomial melting curve: phi = $(NumValue(g.b))T^4 + $(NumValue(g.c))T^3 + $(NumValue(g.d))T^2 + $(NumValue(g.e))T + $(NumValue(g.f)) $(Value(g.T_s)) ≤ T ≤ $(Value(g.T_l))") @@ -273,7 +366,7 @@ function param_info(s::MeltingParam_Quadratic) # info about the struct return MaterialParamsInfo(Equation = L"\phi = 1.0 - ((T_l - T)/(T_l - T_s))^2") end -# Calculation routine +# Calculation routines function compute_meltfraction(p::MeltingParam_Quadratic{_T}, P::Quantity, T::Quantity) where _T @unpack_units T_s,T_l = p @@ -309,6 +402,36 @@ function compute_meltfraction!(ϕ::AbstractArray{_T}, p::MeltingParam_Quadratic{ return nothing end +function compute_dϕdT(p::MeltingParam_Quadratic{_T}, P::Quantity, T::Quantity) where _T + @unpack_units T_s,T_l = p + + dϕdT = (2T_l - 2T) / ((T_l - T_s)^2) + if T>T_l || TT_l || TT_l] .= 0. + + return nothing +end + # Print info function show(io::IO, g::MeltingParam_Quadratic) print(io, "Quadratic melting curve: ϕ = 1.0 - ((Tₗ-T)/(Tₗ-Tₛ))² with Tₛ=$(Value(g.T_s)), Tₗ=$(Value(g.T_l)) ") @@ -317,18 +440,18 @@ end """ - ComputeMeltingParam(P,T, p::AbstractPhaseDiagramsStruct) + compute_meltfraction(P,T, p::AbstractPhaseDiagramsStruct) -Computes melt fraction in case we use a phase diagram lookup table. The table should have the collum `:meltFrac` specified. +Computes melt fraction in case we use a phase diagram lookup table. The table should have the column `:meltFrac` specified. """ function compute_meltfraction(p::PhaseDiagram_LookupTable, P::_T,T::_T) where _T return p.meltFrac.(T,P) end """ - ComputeMeltingParam!(ϕ::AbstractArray{<:AbstractFloat}, P::AbstractArray{<:AbstractFloat},T:AbstractArray{<:AbstractFloat}, p::PhaseDiagram_LookupTable) + compute_meltfraction!(ϕ::AbstractArray{<:AbstractFloat}, P::AbstractArray{<:AbstractFloat},T:AbstractArray{<:AbstractFloat}, p::PhaseDiagram_LookupTable) -In-place computation of melt fraction in case we use a phase diagram lookup table. The table should have the collum `:meltFrac` specified. +In-place computation of melt fraction in case we use a phase diagram lookup table. The table should have the column `:meltFrac` specified. """ function compute_meltfraction!(ϕ::AbstractArray{_T}, p::PhaseDiagram_LookupTable, P::AbstractArray{_T}, T::AbstractArray{_T}) where _T ϕ[:] = p.meltFrac.(T,P) @@ -336,68 +459,82 @@ function compute_meltfraction!(ϕ::AbstractArray{_T}, p::PhaseDiagram_LookupTabl return nothing end -# Computational routines needed for computations with the MaterialParams structure -function compute_meltfraction(s::AbstractMaterialParamsStruct, P::_T=zero(_T),T::_T=zero(_T)) where {_T} - if isempty(s.Melting) #in case there is a phase with no melting parametrization - return zero(_T) - end - return compute_meltfraction(s.Melting[1], P,T) -end - """ - ComputeMeltingParam!(ϕ::AbstractArray{<:AbstractFloat}, Phases::AbstractArray{<:Integer}, P::AbstractArray{<:AbstractFloat},T::AbstractArray{<:AbstractFloat}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct}) + compute_dϕdT(P,T, p::AbstractPhaseDiagramsStruct) -In-place computation of density `rho` for the whole domain and all phases, in case a vector with phase properties `MatParam` is provided, along with `P` and `T` arrays. +Computes derivative of melt fraction vs T in case we use a phase diagram lookup table. The table should have the column `:meltFrac` specified. +The derivative is computed by finite differencing. """ -compute_meltfraction(args...) = compute_param(compute_meltfraction, args...) -compute_meltfraction!(args...) = compute_param!(compute_meltfraction, args...) - -#= -function compute_meltfraction!(ϕ::AbstractArray{<:AbstractFloat, N}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct, 1}, Phases::AbstractArray{<:Integer, N}, P::AbstractArray{<:AbstractFloat, N},T::AbstractArray{<:AbstractFloat, N}) where N +function compute_dϕdT(p::PhaseDiagram_LookupTable, P::_T,T::_T) where _T + + dT = (maximum(T) - minimum(T))/2.0*1e-6 + 1e-6 # 1e-6 of the average T value + ϕ1 = p.meltFrac.(T .+ dT ,P) + ϕ0 = p.meltFrac.(T ,P) + dϕdT = (ϕ1-ϕ0)/dT - for i = 1:length(MatParam) + return dϕdT +end - if length(MatParam[i].Melting)>0 +""" + compute_dϕdT!(dϕdT::AbstractArray{<:AbstractFloat}, P::AbstractArray{<:AbstractFloat},T:AbstractArray{<:AbstractFloat}, p::PhaseDiagram_LookupTable) - # Create views into arrays (so we don't have to allocate) - ind = Phases .== MatParam[i].Phase; - ϕ_local = view(ϕ, ind ) - P_local = view(P , ind ) - T_local = view(T , ind ) +In-place computation of melt fraction in case we use a phase diagram lookup table. The table should have the column `:meltFrac` specified. +The derivative is computed by finite differencing. +""" +function compute_meltfraction!(dϕdT::AbstractArray{_T}, p::PhaseDiagram_LookupTable, P::AbstractArray{_T}, T::AbstractArray{_T}) where _T + + dT = (maximum(T) - minimum(T))/2.0*1e-6 + 1e-6 # 1e-6 of the average T value + ϕ1 = p.meltFrac.(T .+ dT ,P) + ϕ0 = p.meltFrac.(T ,P) + dϕdT = (ϕ1-ϕ0)/dT - compute_meltfraction!(ϕ_local, MatParam[i].Melting[1], P_local, T_local) - end + return nothing +end +# Computational routines needed for computations with the MaterialParams structure +function compute_meltfraction(s::AbstractMaterialParamsStruct, P::_T=zero(_T),T::_T=zero(_T)) where {_T} + if isempty(s.Melting) #in case there is a phase with no melting parametrization + return zero(_T) end + return compute_meltfraction(s.Melting[1], P,T) +end - return nothing +function compute_dϕdT(s::AbstractMaterialParamsStruct, P::_T=zero(_T),T::_T=zero(_T)) where {_T} + if isempty(s.Melting) #in case there is a phase with no melting parametrization + return zero(_T) + end + return compute_dϕdT(s.Melting[1], P,T) end +""" + ϕ = compute_meltfraction(Phases::AbstractArray{<:Integer}, P::AbstractArray{<:AbstractFloat},T::AbstractArray{<:AbstractFloat}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct}) +Computation of melt fraction ϕ for the whole domain and all phases, in case an array with phase properties `MatParam` is provided, along with `P` and `T` arrays. """ - ComputeMeltingParam!(ϕ::AbstractArray{<:AbstractFloat}, Phases::AbstractArray{<:AbstractFloat}, P::AbstractArray{<:AbstractFloat},T::AbstractArray{<:AbstractFloat}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct}) +compute_meltfraction(args...) = compute_param(compute_meltfraction, args...) -In-place computation of density `rho` for the whole domain and all phases, in case a vector with phase properties `MatParam` is provided, along with `P` and `T` arrays. """ -function compute_meltfraction!(ϕ::AbstractArray{<:AbstractFloat, N}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct, 1}, PhaseRatios::AbstractArray{<:AbstractFloat, M}, P::AbstractArray{<:AbstractFloat, N},T::AbstractArray{<:AbstractFloat, N}) where {N,M} + compute_meltfraction(ϕ::AbstractArray{<:AbstractFloat}, Phases::AbstractArray{<:Integer}, P::AbstractArray{<:AbstractFloat},T::AbstractArray{<:AbstractFloat}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct}) - ϕ .= 0.0 - for i = 1:length(MatParam) - - ϕ_local = zeros(size(ϕ)) - Fraction = selectdim(PhaseRatios,M,i); - if (maximum(Fraction)>0.0) & (length(MatParam[i].Melting)>0) +In-place computation of melt fraction ϕ for the whole domain and all phases, in case an array with phase properties `MatParam` is provided, along with `P` and `T` arrays. +""" +compute_meltfraction!(args...) = compute_param!(compute_meltfraction, args...) - compute_meltfraction!(ϕ_local, MatParam[i].Melting[1] , P, T) +""" + ϕ = compute_dϕdT(Phases::AbstractArray{<:Integer}, P::AbstractArray{<:AbstractFloat},T::AbstractArray{<:AbstractFloat}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct}) - ϕ .= ϕ .+ ϕ_local.*Fraction - end +Computates the derivative of melt fraction ϕ versus temperature `T` for the whole domain and all phases, in case an array with phase properties `MatParam` is provided, along with `P` and `T` arrays. +This is employed in computing latent heat terms in an implicit manner, for example +""" +compute_dϕdT(args...) = compute_param(compute_dϕdT, args...) - end - return nothing -end -=# +""" + compute_dϕdT!(ϕ::AbstractArray{<:AbstractFloat}, Phases::AbstractArray{<:Integer}, P::AbstractArray{<:AbstractFloat},T::AbstractArray{<:AbstractFloat}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct}) +Computates the derivative of melt fraction `ϕ` versus temperature `T`, ``\\partial \\phi} \\over {\\partial T}`` for the whole domain and all phases, in case an array with phase properties `MatParam` is provided, along with `P` and `T` arrays. +This is employed in computing latent heat terms in an implicit manner, for example +""" +compute_dϕdT!(args...) = compute_param!(compute_dϕdT, args...) end \ No newline at end of file diff --git a/test/test_MeltingParam.jl b/test/test_MeltingParam.jl index 74104955c..b087a798d 100644 --- a/test/test_MeltingParam.jl +++ b/test/test_MeltingParam.jl @@ -16,7 +16,7 @@ T_nd = Float64.(T/CharUnits_GEO.Temperature) # Caricchi parameterization [in ND numbers, which is anyways the typical use case] p = MeltingParam_Caricchi() phi_dim = zeros(size(T)) -compute_meltfraction!(phi_dim, p, zeros(size(T)), ustrip.(T)) +compute_meltfraction!(phi_dim, p, zeros(size(T)), ustrip.(T)) phi_dim1 = zeros(size(phi_dim)) compute_meltfraction!(phi_dim1, p, zeros(size(T)), ustrip.(T)) # in-place routine @@ -35,6 +35,12 @@ Phi_anal = 1.0 .- Phi_solid @test sum(phi_dim1 - Phi_anal) < 1e-12 @test sum(phi_nd - Phi_anal) < 1e-12 +# test derivative vs T +dϕdT_dim = zeros(size(T)) +compute_dϕdT!(dϕdT_dim, p, zeros(size(T)), ustrip.(T)) +@test sum(dϕdT_dim) ≈ 0.008102237679214096 + + #------------------------------ # 5th order polynomial p = MeltingParam_5thOrder(); @@ -65,6 +71,12 @@ phi = zeros(size(Tdata)) compute_meltfraction!(phi, p, zeros(size(Tdata)), ustrip.(Tdata)) @test norm(data[:,2]-phi) ≈ 0.07151515017819135 + +# test derivative vs T +dϕdT_dim = zeros(size(T)) +compute_dϕdT!(dϕdT_dim, p, zeros(size(T)), ustrip.(T)) +@test sum(dϕdT_dim) ≈ 0.006484458453421382 + #------------------------------ @@ -99,6 +111,12 @@ phi = zeros(size(Tdata)) compute_meltfraction!(phi, p, zeros(size(Tdata)), ustrip.(Tdata)) @test norm(data[:,2]-phi) ≈ 0.0678052542705406 + +# test derivative vs T +dϕdT_dim = zeros(size(T)) +compute_dϕdT!(dϕdT_dim, p, zeros(size(T)), ustrip.(T)) +@test sum(dϕdT_dim) ≈ 0.00830985782591842 + #------------------------------ @@ -107,6 +125,10 @@ compute_meltfraction!(phi, p, zeros(size(Tdata)), ustrip.(Tdata)) p = MeltingParam_Quadratic(); compute_meltfraction!(phi_dim, p, zeros(size(T)), ustrip.(T)) @test sum(phi_dim) ≈ 5.0894901144641 + +dϕdT_dim = zeros(size(T)) +compute_dϕdT!(dϕdT_dim, p, zeros(size(T)), ustrip.(T)) +@test sum(dϕdT_dim) ≈ 0.009365244536940681 #------------------------------ @@ -137,13 +159,16 @@ Phases[:,:,80:end] .= 3 Phases[:,:,90:end] .= 4 ϕ = zeros(size(Phases)) +dϕdT = zeros(size(Phases)) T = ones(size(Phases))*1500 P = ones(size(Phases))*10 compute_meltfraction!(ϕ, Mat_tup, Phases, P,T) #allocations coming from computing meltfraction using PhaseDiagram_LookupTable - @test sum(ϕ)/n^3 ≈ 0.7463001302812086 +compute_dϕdT!(dϕdT, Mat_tup, Phases, P,T) +@test sum(dϕdT)/n^3 ≈ 0.000176112129245805 + # test computing material properties when we have PhaseRatios, instead of Phase numbers PhaseRatio = zeros(n,n,n,4); @@ -156,6 +181,8 @@ end compute_meltfraction!(ϕ, Mat_tup, PhaseRatio, P,T) @test sum(ϕ)/n^3 ≈ 0.7463001302812086 +compute_dϕdT!(dϕdT, Mat_tup, PhaseRatio, P,T) +@test sum(dϕdT)/n^3 ≈ 0.000176112129245805 end diff --git a/test/test_ZirconSaturation.jl b/test/test_ZirconSaturation.jl index 655005a6e..e2258a9ef 100644 --- a/test/test_ZirconSaturation.jl +++ b/test/test_ZirconSaturation.jl @@ -68,7 +68,6 @@ compute_zirconsaturation!(Fzrs, Mat_tup, Phases, P, T) #allocations coming from @test sum(Fzrs)/n^2 ≈ 0.8028012642752728 - # test computing material properties when we have PhaseRatios, instead of Phase numbers PhaseRatio = zeros(n,n,3); for i in CartesianIndices(Phases)