Add optional check that units for linked variables match

src/VariableDomain.jl

@@ -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; 
@@ -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)
@@ -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
 ##################################################################

src/reactioncatalog/GridForcings.jl

@@ -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)
 

src/reactioncatalog/Reservoirs.jl

@@ -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

src/reactioncatalog/VariableStats.jl

@@ -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[]
@@ -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!(
@@ -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
@@ -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)),
     ]
 
@@ -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)),
     ]

src/reactionmethods/RateStoich.jl

@@ -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
@@ -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,))
         )
 

test/configratestoich.yaml

@@ -12,6 +12,7 @@ model1:
                 constnames: ["volume"]
             variable_attributes:                       
                 volume:initial_value:        10.0
+                volume:units:                m3
 
         reservoir_H2S:
             class: ReactionReservoirTotal

test/configreservoirs.yaml

@@ -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:
 
@@ -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
@@ -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

test/runbasetests.jl

@@ -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
 

test/runfluxtests.jl

@@ -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    
 

test/rungridstests.jl

@@ -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
@@ -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)])
@@ -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)])

test/runratestoichtests.jl

@@ -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)

test/runreservoirtests.jl

@@ -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