Skip to content

Commit

Permalink
update
Browse files Browse the repository at this point in the history
  • Loading branch information
TorkelE committed Nov 14, 2023
1 parent a68e820 commit d607736
Show file tree
Hide file tree
Showing 4 changed files with 49 additions and 5 deletions.
10 changes: 10 additions & 0 deletions docs/src/catalyst_functionality/chemistry_related_functionality.md
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,16 @@ Finally, it is possible to check whether a species is a compound or not using th
iscompound(CO2)
```

!!! note
Currently, compounds with components that depend on variables that are not `t` are not supported. E.g.
```julia
@variables x y
@species O(x, y)
@compound O2 ~ 2O
```
will currently throw an error.


Compound components that are also compounds are allowed, e.g. we can create a carbonic acid compound (H₂CO₃) that consists of CO₂ and H₂O:
```@example chem1
@species H(t)
Expand Down
9 changes: 8 additions & 1 deletion src/chemistry_functionality.jl
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,7 @@ function make_compound(expr)

# Extracts the composite species name, and a Vector{ReactantStruct} of its components.
species_expr = expr.args[2] # E.g. :(CO2(t))
species_expr = insert_independent_variable(species_expr, :t) # Add independent variable (e.g. goes from `CO2` to `CO2(t)`).
species_expr = insert_independent_variable(species_expr, [:t]) # Add independent variable (e.g. goes from `CO2` to `CO2(t)`).
species_name = find_varname_in_declaration(expr.args[2]) # E.g. :CO2
composition = Catalyst.recursive_find_reactants!(expr.args[3], 1, Vector{ReactantStruct}(undef, 0))

Expand All @@ -61,13 +61,20 @@ function make_compound(expr)

# Creates the found expressions that will create the compound species.
# The `Expr(:escape, :(...))` is required so that teh expressions are evaluated in the scope the users use the macro in (to e.g. detect already exiting species).
non_t_iv_error_check_expr = Expr(:escape, :(@species $species_expr)) # E.g. `@species CO2(t)`
species_declaration_expr = Expr(:escape, :(@species $species_expr)) # E.g. `@species CO2(t)`
compound_designation_expr = Expr(:escape, :($species_name = ModelingToolkit.setmetadata($species_name, Catalyst.CompoundSpecies, true))) # E.g. `CO2 = ModelingToolkit.setmetadata(CO2, Catalyst.CompoundSpecies, true)`
components_designation_expr = Expr(:escape, :($species_name = ModelingToolkit.setmetadata($species_name, Catalyst.CompoundComponents, $components))) # E.g. `CO2 = ModelingToolkit.setmetadata(CO2, Catalyst.CompoundSpecies, [C, O])`
coefficients_designation_expr = Expr(:escape, :($species_name = ModelingToolkit.setmetadata($species_name, Catalyst.CompoundCoefficients, $coefficients))) # E.g. `CO2 = ModelingToolkit.setmetadata(CO2, Catalyst.CompoundSpecies, [1, 2])`

# Currently, non-t independent variables are not supported for compounds. If there are any like these, we throw an error:
non_t_iv_error_check_expr = Expr(:escape, :(issetequal(unique(reduce(vcat, arguments.(ModelingToolkit.unwrap.($components)))), [t]) || error("Currently, compounds depending on components that are not \"t\" are not supported.")))

println(non_t_iv_error_check_expr)

# Returns the rephrased expression.
return quote
$non_t_iv_error_check_expr
$species_declaration_expr
$compound_designation_expr
$components_designation_expr
Expand Down
12 changes: 8 additions & 4 deletions src/expression_utils.jl
Original file line number Diff line number Diff line change
Expand Up @@ -50,19 +50,23 @@ end
# X(t), [metadata=true]
# X(t) = 1.0, [metadata=true]
# (In this example the independent variable t was inserted).
# Extra dispatch is to ensure so that iv is a vector (so that we can handle both single or multiple iv insertion in one function).
# - This way we don't have to make a check for how to create the final expression (which is different for symbol or vector of symbols).
function insert_independent_variable(expr_in, iv)
# If iv is a symbol (e.g. :t), we convert to vector (e.g. [:t]). This way we can handle single and multiple ivs using the same code, without needing to check at the end.
(iv isa Symbol) && (iv = [iv])
# If expr is a symbol, just attach the iv. If not we have to create a new expr and mutate it.
# Because Symbols (a possible input) cannot be mutated, this function cannot mutate the input (would have been easier if Expr input was guaranteed).
(expr_in isa Symbol) && (return :($(expr_in)($iv)))
(expr_in isa Symbol) && (return Expr(:call, expr_in, iv...))
expr = deepcopy(expr_in)

if expr.head == :(=) # Case: :(X = 1.0)
expr.args[1] = :($(expr.args[1])($iv))
expr.args[1] = Expr(:call, expr.args[1], iv...)
elseif expr.head == :tuple
if expr.args[1] isa Symbol # Case: :(X, [metadata=true])
expr.args[1] = :($(expr.args[1])($iv))
expr.args[1] = Expr(:call, expr.args[1], iv...)
elseif (expr.args[1].head == :(=)) && (expr.args[1].args[1] isa Symbol) # Case: :(X = 1.0, [metadata=true])
expr.args[1].args[1] = :($(expr.args[1].args[1])($iv))
expr.args[1].args[1] = Expr(:call, expr.args[1].args[1], iv...)
end
end
(expr == expr_in) && error("Failed to add independent variable $(iv) to expression: $expr_in")
Expand Down
23 changes: 23 additions & 0 deletions test/miscellaneous_tests/compound_macro.jl
Original file line number Diff line number Diff line change
Expand Up @@ -67,6 +67,29 @@ let
@test iscompound(rn.NH3_5)
end

### Test Various Independent Variables ###
let
@variables t x y z
@species C(t) H(x) N(x) O(t) P(t,x) S(x,y)

@compound CO2 ~ C + 2O
@test issetequal(arguments(ModelingToolkit.unwrap(CO2)), [t])

# These 4 tests checks currently non-supported feature.
@test_broken false
@test_broken false
@test_broken false
@test_broken false
# @compound NH4, [output=true] ~ N + 4H
# @compound (H2O = 2.0) ~ 2H + N
# @compound PH4 ~ P + 4H
# @compound SO2 ~ S + 2O
# @test issetequal(arguments(ModelingToolkit.unwrap(NH4)), [x])
# @test issetequal(arguments(ModelingToolkit.unwrap(H2O)), [t, x])
# @test issetequal(arguments(ModelingToolkit.unwrap(PH4)), [t, x])
# @test issetequal(arguments(ModelingToolkit.unwrap(SO2)), [t, x, y])
end

### Other Minor Tests ###

# Test base functionality in two cases.
Expand Down

0 comments on commit d607736

Please sign in to comment.