init

Author: TorkelE

docs/src/model_creation/chemistry_related_functionality.md

@@ -136,4 +136,23 @@ balanced_reaction = balance_reaction(unbalanced_reaction)[1]
 Reactions declared as a part of a `ReactionSystem` (e.g. using the DSL) can be retrieved for balancing using the [`reactions`](@ref) function. Please note that balancing these will not mutate the `ReactionSystem`, but a new reaction system will need to be created using the balanced reactions.
 
 !!! note
-    Reaction balancing is currently not supported for reactions involving compounds of compounds.
+    Reaction balancing is currently not supported for reactions involving compounds of compounds.
+
+### Balancing full systems
+It is possible to balance all the reactions of a reaction system simultaneously using the `balance_system` function. Here, the output is a new system, where all reactions are balanced. E.g. We can use it to balance this system of Methane formation/combustion:
+```@example chem2
+rs = @reaction_network begin
+    @species C(t) O(t) H(t)
+    @compounds begin
+        H2(t) = 2H
+        CH4(t) = C + 4H
+        O2(t) = 2O
+        CO2(t) = C + 2O
+        H2O(t) = 2H + O
+    end
+    1.0, C + H2 --> CH4
+    2.0, CH4 + O2 --> CO2 + H2O
+end
+rs_balanced = balance_system(rs)
+```
+Except for the modified reaction stoichiometries, the new system is identical to the previous one.

src/chemistry_functionality.jl

@@ -346,3 +346,16 @@ function create_matrix(reaction::Catalyst.Reaction)
 
     return A
 end
+
+"""
+    balance_system(rs::ReactionSystem)
+
+From a system, creates a new system where each reaction is a balanced version of the corresponding
+reaction of the original system. For more information, consider the `balance_reaction` function
+(which is internally applied to each system reaction).
+"""
+function balance_system(rs::ReactionSystem)
+    @set! rs.eqs = [(eq isa Reaction) ? balance_reaction(eq) : eq for eq in get_eqs(rs)]
+    @set! rs.rxs = [balance_reaction(rx) for rx in get_rxs(rs)]
+    return rs
+end

test/miscellaneous_tests/reaction_balancing.jl

@@ -425,5 +425,58 @@ let
     @test isequal([rn.CO2,  rn.H2O], brxs.substrates)
     @test isequal([rn.C6H12O6,  rn.O2], brxs.products)
     @test isequal([6, 6], brxs.substoich)
-    @test isequal([1,  6], brxs.prodstoich)    
+    @test isequal([1, 6], brxs.prodstoich)    
+end
+
+
+### Reaction System Balancing ###
+
+# Simple example with multiple reactions.
+let
+    # Creates a reaction system and tis balanced version.
+    rs = @reaction_network begin
+        @species C(t) O(t) H(t)
+        @compounds begin
+            H2(t) = 2H
+            CH4(t) = C + 4H
+            O2(t) = 2O
+            CO2(t) = C + 2O
+            H2O(t) = 2H + O
+        end
+        1.0, C + H2 --> CH4
+        2.0, CH4 + O2 --> CO2 + H2O
+    end
+    rs_balanced = balance_system(rs)
+
+    # Checks that the system is correctly balanced (while not taking order of reactions or of
+    # substrates/products for granted).
+    rx1 = filter(rx -> isequal(rx.rate, 1.0), reactions(rs_balanced))[1]
+    rx2 = filter(rx -> isequal(rx.rate, 2.0), reactions(rs_balanced))[1]
+    @test issetequal(rx1.substoich, [1,2])
+    @test issetequal(rx1.prodstoich, [1])
+    @test issetequal(rx2.substoich, [1,2])
+    @test issetequal(rx2.prodstoich, [1,2])
+end
+
+# Checks for system with non-reaction equations.
+let
+    # Creates a reaction system and tis balanced version.
+    rs = @reaction_network begin
+        @species O(t) H(t)
+        @compounds begin
+            H2(t) = 2H
+            O2(t) = 2O
+            H2O(t) = 2H + O
+        end
+        @equation begin
+            D(V) ~ 1 - V
+        end
+        1.0, H2 + O2 --> H2O
+    end
+    rs_balanced = balance_system(rs)
+
+    # Checks that the system is correctly balanced (while not taking order of reactions or of
+    # substrates/products for granted).
+    @test issetequal(reactions(rs_balanced)[1].substoich, [2,1])
+    @test issetequal(reactions(rs_balanced)[1].prodstoich, [2])
 end