Skip to content

Commit

Permalink
Merge pull request #129 from PALEOtoolkit/check_units
Browse files Browse the repository at this point in the history
Add option to check units for linked variables
  • Loading branch information
sjdaines authored Jun 20, 2024
2 parents 10a156f + 6d90b50 commit b0598ac
Show file tree
Hide file tree
Showing 12 changed files with 92 additions and 29 deletions.
55 changes: 52 additions & 3 deletions src/VariableDomain.jl
Original file line number Diff line number Diff line change
Expand Up @@ -214,22 +214,28 @@ end
[, eltypemap::Dict{String, DataType}],
[, default_host_dependent_field_data=nothing],
[, allow_base_link=true],
[, use_base_vars=String[]])
[. use_base_transfer_jacobian=true],
[, use_base_vars=String[]],
[, check_units_opt=:no])
Allocate or link memory for [`VariableDomain`](@ref)s `vars` in `modeldata` arrays `arrays_idx`
Allocate or link memory for [`VariableDomain`](@ref)s `vars` in `modeldata` array set `arrays_idx`
Element type of allocated Arrays is determined by `eltype(modeldata, arrays_idx)` (the usual case, allowing use of AD types),
which can be overridden by Variable `:datatype` attribute if present (allowing Variables to ignore AD types).
`:datatype` may be either a Julia `DataType` (eg Float64), or a string to be looked up in `eltypemap`.
If `allow_base_link==true`, and any of the following are true a link is made to the base array, instead of allocating:
If `allow_base_link==true`, and any of the following are true a link is made to the base array (`arrays_idx=1`),
instead of allocating a new array in array set `arrays_idx`:
- Variable element type matches `modeldata` base eltype (arrays_idx=1)
- `use_base_transfer_jacobian=true` and Variable `:transfer_jacobian` attribute is set
- Variable full name is in `use_base_vars`
Field data type is determined by Variable `:field_data` attribute, optionally this can take a
`default_host_dependent_field_data` default for Variables with `host_dependent(v)==true` (these are Variables with no Target or no Property linked,
intended to be external dependencies supplied by the solver).
If `check_units_opt != :no` then the `:units` field of linked variable is checked, resulting in either a warning (if `check_units_opt=:warn`)
or error (if `check_units_opt=:error`).
"""
function allocate_variables!(
vars, modeldata::AbstractModelData, arrays_idx::Int;
Expand All @@ -238,11 +244,14 @@ function allocate_variables!(
allow_base_link=true,
use_base_transfer_jacobian=true,
use_base_vars=String[],
check_units_opt=:no,
)

for v in vars
check_lengths(v)

check_units(v; check_units_opt)

data_dims = Tuple(
get_data_dimension(v.domain, dimname)
for dimname in get_attribute(v, :data_dims)
Expand Down Expand Up @@ -355,6 +364,46 @@ function check_lengths(var::VariableDomain)
return nothing
end

"""
check_units(var::VariableDomain; check_units_opt=:warn)
Check that units of all linked Variables match
"""
function check_units(var::VariableDomain; check_units_opt=:warn)

check_units_opt in (:no, :warn, :error) ||
error("check_units(): unsupported option check_units_opt=$check_units_opt (allowed values are :no, :warn, :error)")

check_units_opt == :no && return

var_units = get_attribute(var, :units)

num_errors = 0

for lv in get_all_links(var)
lv_units = get_attribute(lv, :units)

if !_compare_units(var_units, lv_units)
num_errors += 1
@warn "check_units: VariableDomain $(fullname(var)), :units=\"$var_units\" (from master $(typename(var.master.method.reaction)) $(fullname(var.master)))"*
" != $(fullname(lv)), :units=\"$lv_units\" (created by $(typename(lv.method.reaction)) $(fullname(lv.method.reaction)))"
end
end

if check_units_opt == :error && !iszero(num_errors)
error("check_units: VariableDomain $(fullname(var)), :units=$var_units units of linked variables do not match")
end

return num_errors
end

function _compare_units(units1, units2)
# very crude regularization of unit strings: m^3 and m3 are both accepted
units1 = replace(units1, "^"=>"")
units2 = replace(units2, "^"=>"")
return (units1 == units2) || (units1 == "unknown") || (units2 == "unknown")
end

####################################################################
# Manage linked VariableReactions
##################################################################
Expand Down
6 changes: 3 additions & 3 deletions src/reactioncatalog/GridForcings.jl
Original file line number Diff line number Diff line change
Expand Up @@ -94,18 +94,18 @@ function PB.register_methods!(rj::ReactionForceGrid)

vars = [
PB.VarDepScalar("global.tforce", "yr", "historical time at which to apply forcings, present = 0 yr"),
PB.VarProp("F", "", "interpolated forcing"),
PB.VarProp("F", "unknown", "interpolated forcing"),
]
if rj.pars.scale_offset_var[] != 0.0
@info " adding scalar offset from Variable 'scalar_offset'"
push!(vars, PB.VarDepScalar("scalar_offset", "", "scalar offset"))
push!(vars, PB.VarDepScalar("scalar_offset", "unknown", "scalar offset"))
end
PB.setfrozen!(rj.pars.scale_offset_var)

interp_vars = []
for (vname, vlog) in PB.IteratorUtils.zipstrict(rj.pars.interp_vars, rj.pars.interp_log; errmsg="length(interp_vars) != length(interp_log)")
@info " adding interpolation Variable $vname log $vlog"
push!(interp_vars, PB.VarDepScalar(vname, "", "interpolation variable log $vlog"))
push!(interp_vars, PB.VarDepScalar(vname, "unknown", "interpolation variable log $vlog"))
end
PB.setfrozen!(rj.pars.interp_vars, rj.pars.interp_log)

Expand Down
4 changes: 2 additions & 2 deletions src/reactioncatalog/Reservoirs.jl
Original file line number Diff line number Diff line change
Expand Up @@ -740,9 +740,9 @@ function PB.register_methods!(rj::ReactionConst)
for varnameisotope in rj.pars.constnames
varname, IsotopeType = PB.split_nameisotope(varnameisotope, rj.external_parameters)
if rj.scalar
constvar = PB.VarPropScalarStateIndep(varname, "", "constant value", attributes=(:field_data=>IsotopeType, ))
constvar = PB.VarPropScalarStateIndep(varname, "unknown", "constant value", attributes=(:field_data=>IsotopeType, ))
else
constvar = PB.VarPropStateIndep(varname, "", "constant value", attributes=(:field_data=>IsotopeType, ))
constvar = PB.VarPropStateIndep(varname, "unknown", "constant value", attributes=(:field_data=>IsotopeType, ))
end
push!(vars_const, constvar)
end
Expand Down
27 changes: 17 additions & 10 deletions src/reactioncatalog/VariableStats.jl
Original file line number Diff line number Diff line change
Expand Up @@ -77,7 +77,7 @@ function PB.register_methods!(rj::ReactionSum)

# mark all vars_to_add as optional to help diagnose config errors
# default :field_data=>PB.UndefinedData to allow Variable to link to any data type (this is checked later)
v = PB.VarDep(localname => "("*varname*")", "", "")
v = PB.VarDep(localname => "("*varname*")", "unknown", "")

# no length check if scalar sum
if !rj.pars.vectorsum[]
Expand All @@ -88,10 +88,10 @@ function PB.register_methods!(rj::ReactionSum)

if rj.pars.vectorsum[]
methodfn = do_vectorsum
var_sum = PB.VarProp("sum", "", "sum of specified variables")
var_sum = PB.VarProp("sum", "unknown", "sum of specified variables")
else
methodfn = do_scalarsum
var_sum = PB.VarPropScalar("sum", "", "sum of specified variables")
var_sum = PB.VarPropScalar("sum", "unknown", "sum of specified variables")
end

PB.add_method_do!(
Expand All @@ -116,6 +116,13 @@ function PB.register_dynamic_methods!(rj::ReactionSum)
var_sum = PB.get_variable(method_sum, "sum")
vars_to_add = PB.get_variables(method_sum, filterfn = v->v.localname != "sum")

# set units from sum variable
# (may have been set explicitly in yaml file)
sum_units = PB.get_attribute(var_sum, :units)
for v in vars_to_add
PB.set_attribute!(v, :units, sum_units)
end

# check variable components match and update var_sum.components
if rj.pars.component_to_add[] == 0
# add all components of vars_to_add Variables
Expand Down Expand Up @@ -231,11 +238,11 @@ end
function PB.register_methods!(rj::ReactionWeightedMean)

vars = [
PB.VarDep( "var", "measure-1", "variable to calculate weighted mean from",
PB.VarDep( "var", "unknown", "variable to calculate weighted mean from",
attributes=(:field_data=>rj.pars.field_data[],)),
PB.VarDep( "measure", "", "cell area or volume"),
PB.VarDepScalar( "measure_total", "", "total Domain area or volume"),
PB.VarPropScalar("var_mean", "", "weighted mean over Domain area or volume",
PB.VarDep( "measure", "unknown", "cell area or volume"),
PB.VarDepScalar( "measure_total", "unknown", "total Domain area or volume"),
PB.VarPropScalar("var_mean", "unknown", "weighted mean over Domain area or volume",
attributes=(:field_data=>rj.pars.field_data[], :initialize_to_zero=>true, :atomic=>true)),
]

Expand Down Expand Up @@ -289,10 +296,10 @@ end
function PB.register_methods!(rj::ReactionAreaVolumeValInRange)

vars = [
PB.VarDep( "rangevar", "mol m-3", "variable to check within range";
PB.VarDep( "rangevar", "unknown", "variable to check within range";
attributes=(:field_data=>PB.ScalarData, )), # total only, not isotope
PB.VarDep( "measure", "", "cell area or volume"),
PB.VarDepScalar( "measure_total", "", "total Domain area or volume"),
PB.VarDep( "measure", "unknown", "cell area or volume"),
PB.VarDepScalar( "measure_total", "unknown", "total Domain area or volume"),
PB.VarPropScalar("frac", "", "fraction of Domain area or volume in specified range",
attributes=(:initialize_to_zero=>true, :atomic=>true)),
]
Expand Down
4 changes: 2 additions & 2 deletions src/reactionmethods/RateStoich.jl
Original file line number Diff line number Diff line change
Expand Up @@ -133,7 +133,7 @@ function create_ratestoich_method(
etapar = ParDouble("eta", tmpeta, units="per mil", description="constant eta")
end

delta_var = VarDep(delta_varname, "", "generated by RateStoich.rate="*ratevarlocalname,
delta_var = VarDep(delta_varname, "per mil", "generated by RateStoich.rate="*ratevarlocalname,
attributes=(:space=>space,))
else
ratestoich.isotope_data = ScalarData
Expand All @@ -156,7 +156,7 @@ function create_ratestoich_method(
# create variable
push!(
vars_sms,
VarContrib(smsname, "", "generated by RateStoich rate="*ratevarlocalname,
VarContrib(smsname, "mol yr-1", "generated by RateStoich rate="*ratevarlocalname,
attributes=(:field_data=>disotope, :space=>space,))
)

Expand Down
1 change: 1 addition & 0 deletions test/configratestoich.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@ model1:
constnames: ["volume"]
variable_attributes:
volume:initial_value: 10.0
volume:units: m3

reservoir_H2S:
class: ReactionReservoirTotal
Expand Down
8 changes: 7 additions & 1 deletion test/configreservoirs.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,8 @@ model1:
vars_to_add: [1e-4*ocean.vectorsum, ConstS] # 3 + 1 = 4
variable_links:
sum: scalarsum
variable_attributes:
sum%units: mol

ocean:

Expand All @@ -60,13 +62,15 @@ model1:
constnames: ["volume"]
variable_attributes:
volume:initial_value: 10.0
volume:units: m3

const_volume_total:
class: ReactionScalarConst
parameters:
constnames: ["volume_total"]
variable_attributes:
volume_total:initial_value: 10000.0
volume_total:units: m^3

reservoir_const:
class: ReactionReservoirConst
Expand Down Expand Up @@ -99,9 +103,11 @@ model1:
vector_sum:
class: ReactionVectorSum
parameters:
vars_to_add: [volume, 2*T]
vars_to_add: [0.5*C, 2*T]
variable_links:
sum: vectorsum # 30.0
variable_attributes:
sum%units: mol

weightedmean:
class: ReactionWeightedMean
Expand Down
2 changes: 1 addition & 1 deletion test/runbasetests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,7 @@ include("ReactionPaleoMockModule.jl")
modeldata = PB.create_modeldata(model)

@test length(PB.get_unallocated_variables(ocean_domain, modeldata, 1)) == 5
PB.allocate_variables!(ocean_domain, modeldata, 1, hostdep=false)
PB.allocate_variables!(ocean_domain, modeldata, 1; hostdep=false, check_units_opt=:error)
@test length(PB.get_unallocated_variables(ocean_domain, modeldata, 1)) == 1
@test PB.check_ready(ocean_domain, modeldata, throw_on_error=false) == false

Expand Down
2 changes: 1 addition & 1 deletion test/runfluxtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ import PALEOboxes as PB

modeldata = PB.create_modeldata(model)

PB.allocate_variables!(model, modeldata, 1)
PB.allocate_variables!(model, modeldata, 1; check_units_opt=:error)

@test PB.check_ready(model, modeldata) == true

Expand Down
6 changes: 3 additions & 3 deletions test/rungridstests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ include("ReactionPaleoMockModule.jl")

modeldata = PB.create_modeldata(model)

PB.allocate_variables!(ocean_domain, modeldata, 1, hostdep=false)
PB.allocate_variables!(ocean_domain, modeldata, 1; hostdep=false, check_units_opt=:error)
@test length(PB.get_unallocated_variables(ocean_domain, modeldata, 1)) == 1

@test length(hostvars) == 1
Expand Down Expand Up @@ -73,7 +73,7 @@ end
@test PB.get_length(surface_domain) == prod(dsize)

modeldata = PB.create_modeldata(model)
PB.allocate_variables!(model, modeldata, 1)
PB.allocate_variables!(model, modeldata, 1; check_units_opt=:error)
@test PB.check_ready(model, modeldata; throw_on_error=false) == true

modelcreated_vars_dict = Dict([(var.name, var) for var in PB.get_variables(surface_domain, hostdep=false)])
Expand Down Expand Up @@ -104,7 +104,7 @@ end
@test PB.get_length(surface_domain) == prod(dsize)

modeldata = PB.create_modeldata(model)
PB.allocate_variables!(model, modeldata, 1)
PB.allocate_variables!(model, modeldata, 1; check_units_opt=:error)
@test PB.check_ready(model, modeldata, throw_on_error=false) == true

modelcreated_vars_dict = Dict([(var.name, var) for var in PB.get_variables(surface_domain, hostdep=false)])
Expand Down
2 changes: 1 addition & 1 deletion test/runratestoichtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@ include("ReactionRateStoichMock.jl")
domain.grid = PB.Grids.UnstructuredVectorGrid(ncells=10)
@test PB.get_length(domain) == 10
modeldata = PB.create_modeldata(model)
PB.allocate_variables!(model, modeldata, 1)
PB.allocate_variables!(model, modeldata, 1; check_units_opt=:error)
PB.check_ready(model, modeldata)

all_vars = PB.VariableAggregatorNamed(modeldata)
Expand Down
4 changes: 2 additions & 2 deletions test/runreservoirtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -36,13 +36,13 @@ end

# allocate internal variables
modeldata = PB.create_modeldata(model)
PB.allocate_variables!(model, modeldata, 1; hostdep=false)
PB.allocate_variables!(model, modeldata, 1; hostdep=false, check_units_opt=:error)

@test length( PB.get_unallocated_variables(global_domain, modeldata, 1)) == 4
@test PB.check_ready(global_domain, modeldata, throw_on_error=false) == false

# allocate arrays for host dependencies and set data pointers
PB.allocate_variables!(model, modeldata, 1; hostdep=true)
PB.allocate_variables!(model, modeldata, 1; hostdep=true, check_units_opt=:error)
@test PB.check_ready(global_domain, modeldata) == true

# get all variable data arrays
Expand Down

0 comments on commit b0598ac

Please sign in to comment.