Refactor ThermalDiffusionLithiumDensity to remove hardcoded parameters

Project.toml

@@ -8,7 +8,6 @@ ArraysOfArrays = "65a8f2f4-9b39-5baf-92e2-a9cc46fdf018"
 Clustering = "aaaa29a8-35af-508c-8bc3-b662a17a0fe5"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
-SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 Format = "1fa38f19-a742-5d3f-a2b9-30dd87b9d5f8"
 GPUArrays = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7"
@@ -31,6 +30,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"

examples/example_config_files/ThermalDiffusionLithiumDensity/Thermal_Diffusion_Config.yaml

@@ -0,0 +1,17 @@
+model: ThermalDiffusionLithiumDensity
+material: HPGe
+
+
+annealing_temperature_ranges:
+  - T_min: 473K
+    T_max: 873K
+    D0: 2.5e-7m^2/s
+    H: 11800cal
+  - T_min: 873K
+    T_max: 1273K
+    D0: 1.3e-7m^2/s
+    H: 10700cal    
+
+experimental_parameters:
+  a: 21.27
+  b: 2610

src/ImpurityDensities/ThermalDiffusionLithiumDensity.jl

@@ -31,24 +31,55 @@ struct ThermalDiffusionLithiumDensity{T <: SSDFloat} <: AbstractImpurityDensity{
     lithium_diffusivity_in_germanium::T
 end
 
-function calculate_lithium_diffusivity_in_germanium(lithium_annealing_temperature::T)::T where {T <: SSDFloat}
+struct LithiumDiffusionParameters{T <: SSDFloat}
+    T_min::T
+    T_max::T
+    D0::T
+    H::T
+end
+
+struct LithiumSaturationParameters{T <: SSDFloat}
+    a::T
+    b::T
+end
+
+struct ThermalDiffusionLithiumDensityParameters{T <: SSDFloat,N}
+    diffusion::NTuple{N, LithiumDiffusionParameters{T}}
+    saturation::LithiumSaturationParameters{T}
+end
+
+include("ThermalDiffusionLithiumDensityParameters.jl")
+
+function calculate_lithium_diffusivity_in_germanium(lithium_annealing_temperature::T, parameters::NTuple{N, LithiumDiffusionParameters{T}})::T where {T <: SSDFloat, N}
     # D0 [m^2*s^-1]
     # H [cal]
-    D0::T, H::T = ifelse(lithium_annealing_temperature <= 873, T.((2.5e-7, 11800)), T.((1.3e-7, 10700)))
+    D0 = parameters[end].D0
+    H  = parameters[end].H
+    if !(parameters[1].T_min ≤ lithium_annealing_temperature ≤ parameters[end].T_max)
+        throw(ArgumentError("Invalid lithium_annealing_temperature=$(lithium_annealing_temperature): expected $(parameters[1].T_min) ≤ lithium_annealing_temperature ≤ $(parameters[end].T_max)."))
+    end
+    for parameter in parameters
+        if lithium_annealing_temperature <= parameter.T_max
+            D0::T, H::T =  T.((parameter.D0, parameter.H))
+            break
+        end    
+    end
     D0 * exp(-H/(R_gas*lithium_annealing_temperature))
 end
-function calculate_lithium_saturated_density(lithium_annealing_temperature::T)::T where {T <: SSDFloat}
-    exp10(27.27 - 2610.0/lithium_annealing_temperature)
+function calculate_lithium_saturated_density(lithium_annealing_temperature::T, parameters::LithiumSaturationParameters)::T where {T <: SSDFloat}
+    exp10(parameters.a - parameters.b/lithium_annealing_temperature + 6)
 end
 
-function ThermalDiffusionLithiumDensity{T}(
+# N is the number of ranges included in the yaml file
+function ThermalDiffusionLithiumDensity{T}( 
     lithium_annealing_temperature::T,
     lithium_annealing_time::T,
     contact_with_lithium_doped::G,
     inactive_contact_id::Int;
+    model_parameters::ThermalDiffusionLithiumDensityParameters{T} = ThermalDiffusionLithiumParameters(T=T), 
     distance_to_contact::Function = pt::AbstractCoordinatePoint{T} -> ConstructiveSolidGeometry.distance_to_surface(pt, contact_with_lithium_doped),
-    lithium_density_on_contact::T = calculate_lithium_saturated_density(lithium_annealing_temperature),
-    lithium_diffusivity_in_germanium::T = calculate_lithium_diffusivity_in_germanium(lithium_annealing_temperature),
+    lithium_density_on_contact::T = calculate_lithium_saturated_density(lithium_annealing_temperature, model_parameters.saturation),
+    lithium_diffusivity_in_germanium::T = calculate_lithium_diffusivity_in_germanium(lithium_annealing_temperature, model_parameters.diffusion),
 ) where {T <: SSDFloat, G <: Union{<:AbstractGeometry, Nothing}}
     ThermalDiffusionLithiumDensity{T}(lithium_annealing_temperature, lithium_annealing_time, inactive_contact_id, distance_to_contact, lithium_density_on_contact, lithium_diffusivity_in_germanium)
 end
@@ -59,7 +90,8 @@ function ImpurityDensity(T::DataType, t::Val{:li_diffusion}, dict::AbstractDict,
     contact_with_lithium_doped = haskey(dict, "contact_with_lithium_doped") ? dict["contact_with_lithium_doped"] : nothing # you don't have to pass the geometry of doped contact only when the distance_to_contact is passed
     inactive_contact_id = get(dict, "doped_contact_id", -1)
     inactive_contact_id < 1 && error("Invalid doped_contact_id: missing or misspelled key")
-    ThermalDiffusionLithiumDensity{T}(lithium_annealing_temperature, lithium_annealing_time, contact_with_lithium_doped, inactive_contact_id)
+    model_parameters = haskey(dict,"model_parameters") ? ThermalDiffusionLithiumParameters(dict["model_parameters"]; T) : ThermalDiffusionLithiumParameters(; T)
+    ThermalDiffusionLithiumDensity{T}(lithium_annealing_temperature, lithium_annealing_time, contact_with_lithium_doped, inactive_contact_id, model_parameters)
 end
 
 function get_impurity_density(li_diffusion::ThermalDiffusionLithiumDensity{T}, pt::AbstractCoordinatePoint{T})::T where {T <: SSDFloat}

src/ImpurityDensities/ThermalDiffusionLithiumDensityParameters.jl

@@ -0,0 +1,93 @@
+function _ThermalDiffusionLithiumParameters(
+        config::AbstractDict;
+        T::Type = Float32,
+        temperature_ranges::AbstractVector{<:AbstractDict} = config["annealing_temperature_ranges"] ,
+        a::Union{Real, String} = config["experimental_parameters"]["a"],
+        b::Union{Real, String} = config["experimental_parameters"]["b"],
+        input_units::Union{Missing, NamedTuple} = missing
+    )
+
+    temperature_unit = !ismissing(input_units) ? input_units.temperature : internal_temperature_unit
+    diffusivity_unit = !ismissing(input_units) ?
+        input_units.length^2 / internal_time_unit : internal_length_unit^2/internal_time_unit
+
+    temp_ranges_parameters = Vector{LithiumDiffusionParameters{T}}()
+
+    for r in temperature_ranges
+        push!(temp_ranges_parameters,
+            LithiumDiffusionParameters{T}(
+                _parse_value(T, r["T_min"], temperature_unit),
+                _parse_value(T, r["T_max"], temperature_unit),
+                _parse_value(T, r["D0"], diffusivity_unit),
+                _parse_value(T, r["H"], u"cal") / ustrip(Unitful.NoUnits, one(T) * u"cal" / internal_energy_unit) # convert to cal internally
+            )
+        )
+    end
+
+    diffusion_params = Tuple(
+        LithiumDiffusionParameters{T}(
+            _parse_value(T, r["T_min"], temperature_unit),
+            _parse_value(T, r["T_max"], temperature_unit),
+            _parse_value(T, r["D0"], diffusivity_unit),
+            _parse_value(T, r["H"], u"cal") / ustrip(Unitful.NoUnits, one(T) * u"cal" / internal_energy_unit) # convert to cal internally
+        )
+        for r in temperature_ranges
+    )
+
+    for i in eachindex(diffusion_params)
+        if diffusion_params[i].T_max < diffusion_params[i].T_min
+            throw(ConfigFileError("Invalid annealing temperature range $i: T_max must be ≥ T_min."))
+        end
+        if i>1 && diffusion_params[i-1].T_max != diffusion_params[i].T_min
+            throw(ConfigFileError("Annealing temperature ranges must be contiguous and increasing: range $(i-1) ends at a different temperature than range $i starts."))
+        end
+    end
+
+
+    saturation_params = LithiumSaturationParameters{T}(
+        _parse_value(T, a, Unitful.NoUnits), _parse_value(T, b, Unitful.NoUnits)
+    )
+
+    return ThermalDiffusionLithiumDensityParameters(
+        diffusion_params,
+        saturation_params
+    )
+end
+
+
+function ThermalDiffusionLithiumParameters(config::AbstractDict, input_units::Union{Missing, NamedTuple} = missing; kwargs...)
+    if !haskey(config, "annealing_temperature_ranges") 
+        throw(ConfigFileError("ThermalDiffusionLithiumDensity config file needs entry 'annealing_temperature_ranges'.")) 
+    elseif isempty(config["annealing_temperature_ranges"])
+        throw(ConfigFileError("ThermalDiffusionLithiumDensity config file needs at least $minimum_num_of_temp_ranges entries in 'annealing_temperature_ranges'."))
+    end 
+
+    for (i, temperature_range) in enumerate(config["annealing_temperature_ranges"])
+        for axis in ("T_min", "T_max", "D0", "H")
+            if !haskey(temperature_range, axis)
+                throw(ConfigFileError("ThermalDiffusionLithiumDensity config file needs entry 'annealing_temperature_ranges[$i]/$(axis)'."))
+            end
+        end
+        for k in keys(temperature_range)
+            if !(k in ("T_min", "T_max", "D0", "H"))
+                @warn "Encountered unexpected key `$(k)` in 'annealing_temperature_ranges[$i]': Allowed keys are `T_min`, `T_max`, `D0` and `H`. Key `$(k)` will be ignored."
+            end
+        end
+    end
+
+
+    if !haskey(config, "experimental_parameters") 
+        throw(ConfigFileError("ThermalDiffusionLithiumDensity config file needs entry 'experimental_parameters'.")) 
+    elseif !haskey(config["experimental_parameters"], "a") 
+        throw(ConfigFileError("ThermalDiffusionLithiumDensity config file needs entry 'experimental_parameters/a'.")) 
+    elseif !haskey(config["experimental_parameters"], "b") 
+        throw(ConfigFileError("ThermalDiffusionLithiumDensity config file needs entry 'experimental_parameters/b'.")) 
+    end
+
+    _ThermalDiffusionLithiumParameters(config; input_units, kwargs...)
+end
+
+const default_ThermalDiffusionLithiumDensity_config_file = joinpath(get_path_to_example_config_files(), "ThermalDiffusionLithiumDensity/Thermal_Diffusion_Config.yaml")
+function ThermalDiffusionLithiumParameters(config_filename::AbstractString = default_ThermalDiffusionLithiumDensity_config_file, input_units::Union{Missing, NamedTuple} = missing; kwargs...)
+    ThermalDiffusionLithiumParameters(parse_config_file(config_filename), input_units; kwargs...)
+end