first commit

src/reaktoro_pse/core/reaktoro_inputs.py

@@ -167,24 +167,30 @@ def add_specs(self, specs_object, assert_charge_neutrality, dissolveSpeciesInRkt
         pressure_not_set = True
         temperature_not_set = True
         for input_name, _ in self.state.inputs.items():
-            if input_name == "temperature":
+            if input_name == RktInputTypes.temperature:
                 specs_object.temperature()
                 temperature_not_set = False
-                self.rkt_inputs["temperature"] = self.state.inputs["temperature"]
-                self.rkt_inputs["temperature"].set_lower_bound(0)
-            elif input_name == "pressure":
+                self.rkt_inputs[RktInputTypes.temperature] = self.state.inputs[
+                    RktInputTypes.temperature
+                ]
+                self.rkt_inputs[RktInputTypes.temperature].set_lower_bound(0)
+            elif input_name == RktInputTypes.pressure:
                 specs_object.pressure()
                 pressure_not_set = False
-                self.rkt_inputs["pressure"] = self.state.inputs["pressure"]
-                self.rkt_inputs["pressure"].set_lower_bound(0)
-            elif input_name == "enthalpy":
+                self.rkt_inputs[RktInputTypes.pressure] = self.state.inputs[
+                    RktInputTypes.pressure
+                ]
+                self.rkt_inputs[RktInputTypes.pressure].set_lower_bound(0)
+            elif input_name == RktInputTypes.enthalpy:
                 specs_object.enthalpy()
-                self.rkt_inputs["enthalpy"] = self.state.inputs["enthalpy"]
-                self.rkt_inputs["enthalpy"].set_lower_bound(None)
-            elif input_name == "pH":
+                self.rkt_inputs[RktInputTypes.enthalpy] = self.state.inputs[
+                    RktInputTypes.enthalpy
+                ]
+                self.rkt_inputs[RktInputTypes.enthalpy].set_lower_bound(None)
+            elif input_name == RktInputTypes.pH:
                 specs_object.pH()
-                self.rkt_inputs["pH"] = self.state.inputs["pH"]
-                self.rkt_inputs["pH"].set_lower_bound(0)
+                self.rkt_inputs[RktInputTypes.pH] = self.state.inputs[RktInputTypes.pH]
+                self.rkt_inputs[RktInputTypes.pH].set_lower_bound(0)
             else:
                 pass
         if pressure_not_set:
@@ -196,11 +202,14 @@ def add_specs(self, specs_object, assert_charge_neutrality, dissolveSpeciesInRkt
         if assert_charge_neutrality:
             specs_object.charge()
             if self.neutrality_ion is not None:
-                self.ignore_elements_for_constraints.append(self.neutrality_ion)
+                if self.neutrality_ion == RktInputTypes.pH:
+                    specs_object.openTo("H+")
+                else:
+                    self.ignore_elements_for_constraints.append(self.neutrality_ion)
 
-                if self.neutrality_ion not in specs_object.namesInputs():
-                    # needs to be a species!
-                    specs_object.openTo(self.neutrality_ion)
+                    if self.neutrality_ion not in specs_object.namesInputs():
+                        # needs to be a species!
+                        specs_object.openTo(self.neutrality_ion)
 
         self._find_element_sums()
         # add/check if vars in rkt Inputs
@@ -337,7 +346,7 @@ def default_speciation(self):
             "NaOH": {"Na": 1, "O": 1, "H": 1},
             "H": {"H": 1},
             "OH": {"O": 1, "H": 1},
-            "H2O_evaporation": {"O": 1, "H": 2},
+            "H2O_evaporation": {"O": -1, "H": -2},
         }
 
     def get_modifier_mw(self, elemental_composition):

src/reaktoro_pse/examples/reaktoro_pse_to_phreeqc_comparison/standard_model.py

@@ -62,7 +62,7 @@ def build_modification_example(water_comp):
     m.water_recovery.fix()
     m.eq_water_fow = Constraint(
         expr=m.water_recovery
-        == -m.modified_properties_water_removal / m.feed_composition["H2O"]
+        == m.modified_properties_water_removal / m.feed_composition["H2O"]
     )
     return m
 
@@ -94,11 +94,11 @@ def add_standard_properties(m):
             "H2O_evaporation": m.modified_properties_water_removal,
             "NaOH": m.base_addition,
         },
-        dissolve_species_in_reaktoro=False,
+        dissolve_species_in_reaktoro=True,
         # we can use default converter as its defined for default database (Phreeqc and pitzer)
         # we are modifying state and must speciate inputs before adding acid to find final prop state.
         build_speciation_block=True,
-        reaktoro_solve_options={"open_species_on_property_block": ["H+", "OH-"]},
+        # reaktoro_solve_options={"open_species_on_property_block": ["OH-", "H2O"]},
         jacobian_options={
             "user_scaling": {
                 ("saturationIndex", "Calcite"): 1,
@@ -116,7 +116,7 @@ def scale_model(m):
         iscale.set_scaling_factor(
             m.feed_composition[key], 1 / m.feed_composition[key].value
         )
-    iscale.set_scaling_factor(m.water_recovery, 1)
+    iscale.set_scaling_factor(m.water_recovery, 1 / 1)
     iscale.set_scaling_factor(m.acid_addition, 1 / 0.001)
     iscale.set_scaling_factor(m.base_addition, 1 / 0.001)
 

src/reaktoro_pse/examples/reaktoro_pse_to_phreeqc_comparison/water_removal_comparison.py

@@ -31,7 +31,9 @@ def main(save_fig=False, show_fig=True):
     standardModel.initialize(m)
     reaktoro_output_dict = {}
     reaktoro_output_dict["water_removal_percent_sweep"] = {}
+    print(phreeqc_config["water_removal_percent"])
     for wr in phreeqc_config["water_removal_percent"]:
+        print(wr)
         m.water_recovery.fix(wr / 100)
         standardModel.solve(m)
         compUtils.get_reaktoro_solved_outputs(

src/reaktoro_pse/examples/simple_desalination.py

@@ -139,7 +139,7 @@ def eq_desal_composition(fs, key):
 
         # However, this can result in incorrect speciation for some databases, so please use with caution.
 
-        species_to_open = ["H+", "OH-"]
+        species_to_open = ["OH-"]
     else:
         species_to_open = None
     m.eq_desal_properties = ReaktoroBlock(

src/reaktoro_pse/reaktoro_block.py

@@ -150,20 +150,23 @@ class ReaktoroBlockData(ProcessBlockData):
         ConfigValue(
             default="Cl",
             domain=str,
-            description="Ion to use for maintaining charge neutrality",
-            doc="""This will unfix specified ion during equilibrium calculations while enforcing charge==0 constraint in reaktoro""",
+            description="Ion to use for maintaining charge neutrality (pH, ion, or element)",
+            doc="""This will unfix specified ion during equilibrium calculations while enforcing charge==0 constraint
+              in reaktoro, if exact speciation is provided, pH will be used by default""",
         ),
     )
     CONFIG.declare(
         "assert_charge_neutrality_on_all_blocks",
         ConfigValue(
-            default=False,
+            default=True,
             domain=bool,
             description="Defines if charge neutrality should be applied to both speciation and property block",
             doc="""
             When user provides chemical inputs, its assumed that user wants to modify an equilibrated state, 
             as such a speciation and property block will be built. Charge neutrality would only be applied on speciation block if enabled
-            if this option is set to True then charge neutrality will also be applied on property block """,
+            if this option is set to True then charge neutrality will also be applied on property block, 
+            by default this is set to true, and property block is charge neutralized against pH (H+) if present 
+            in system""",
         ),
     )
 
@@ -511,26 +514,41 @@ def build_rkt_inputs(
         block.rkt_inputs.register_free_elements(self.config.aqueous_phase.free_element)
         block.rkt_inputs.register_free_elements(self.config.liquid_phase.free_element)
         # register charge neutrality
-        if (
-            speciation_block_built == False
-            or self.config.assert_charge_neutrality_on_all_blocks
-        ):
-            # only ensure charge neutrality when doing first calculation
+        # only ensure charge neutrality when doing first calculation
+        if speciation_block_built == False:
+            # if exact speciation - then charge balance should be done on pH
+            if self.config.exact_speciation == True:
+                # only do so if we have 'H+' in species
+                ion_for_balancing = "pH"
+                if "H+" in block.rkt_state.database_species:
+                    assert_charge_neutrality = True
+
+            else:
+                assert_charge_neutrality = self.config.assert_charge_neutrality
+                ion_for_balancing = self.config.charge_neutrality_ion
+
             block.rkt_inputs.register_charge_neutrality(
-                assert_neutrality=self.config.assert_charge_neutrality,
-                ion=self.config.charge_neutrality_ion,
+                assert_neutrality=assert_charge_neutrality,
+                ion=ion_for_balancing,
             )
             block.rkt_inputs.configure_specs(
                 dissolve_species_in_rkt=self.config.dissolve_species_in_reaktoro,
                 exact_speciation=self.config.exact_speciation,
             )
         else:
             # if we have built a speciation block, the feed should be charge neutral and
-            # exact speciation is provided
+            # exact speciation is provided, then we balance on pH only,
+            # user can disable this self.assert_charge_neutrality_on_all_blocks
+            if self.config.assert_charge_neutrality_on_all_blocks:
+                # only do so if we have 'H+' in species
+                if "H+" in block.rkt_state.database_species:
+                    assert_charge_neutrality = True
+            else:
+                assert_charge_neutrality = False
             block.rkt_inputs.register_charge_neutrality(
-                assert_neutrality=False, ion=self.config.charge_neutrality_ion
+                assert_neutrality=assert_charge_neutrality,
+                ion="pH",
             )
-
             block.rkt_inputs.configure_specs(
                 dissolve_species_in_rkt=self.config.dissolve_species_in_reaktoro,
                 exact_speciation=True,