Merge branch 'main' into ph_charge_balance_fix

README.md

@@ -84,6 +84,7 @@ This option will force Ipopt to use least squares method to calculate dual infea
 
     cy_solver = get_solver(solver="cyipopt-watertap")
     cy_solver.options['recalc_y']='yes'
+    cy_solver.options["recalc_y_feas_tol"] = 1e-2
 
 
 B. Use exact derivatives instead of numeric

src/reaktoro_pse/core/pyomo_property_writer/property_functions.py

@@ -10,6 +10,7 @@
 # "https://github.com/watertap-org/reaktoro-pse/"
 #################################################################################
 from pyomo.environ import log10, log, exp
+from idaes.core.util.math import smooth_max
 
 
 def build_scaling_tendency_constraint(rkt_output_object):
@@ -24,6 +25,12 @@ def build_scaling_tendency_constraint(rkt_output_object):
     )
 
 
+# log10(user_output_var) == smooth_max(
+#     build_properties[("saturationIndex", rkt_output_object.property_index)].pyomo_var,
+#     -5,
+# )
+
+
 def build_ph_constraint(rkt_output_object):
     user_output_var = rkt_output_object.pyomo_var
     build_properties = rkt_output_object.pyomo_build_options.properties

src/reaktoro_pse/core/reaktoro_block_builder.py

@@ -54,14 +54,33 @@ def __init__(self, block, reaktoro_solver, build_on_init=True):
         if isinstance(self.solver, ReaktoroSolver) == False:
             raise TypeError("Reaktoro block builder requires a ReaktoroSolver class")
         self.configure_jacobian_scaling()
+        self.reaktoro_initialize_function = None  # used to provide external solve call
+        self.display_reaktoro_state_function = (
+            None  # used to specifying external function to display rkt state
+        )
+        self.build_output_vars()
         if build_on_init:  # option to support legacy implementation
             self.build_reaktoro_block()
 
-    def build_reaktoro_block(self):
+    def build_reaktoro_block(
+        self,
+        gray_box_model=None,
+        reaktoro_initialize_function=None,
+        display_reaktoro_state_function=None,
+    ):
         """build reaktoro model"""
-        external_model = ReaktoroGrayBox()
-        external_model.configure(self.solver)
-        self.block.reaktoro_model = ExternalGreyBoxBlock(external_model=external_model)
+        if gray_box_model is None:
+            external_model = ReaktoroGrayBox()
+            external_model.configure(self.solver)
+            self.block.reaktoro_model = ExternalGreyBoxBlock(
+                external_model=external_model
+            )
+        else:
+            self.block.reaktoro_model = gray_box_model
+        if reaktoro_initialize_function is not None:
+            self.reaktoro_initialize_function = reaktoro_initialize_function
+        if display_reaktoro_state_function is not None:
+            self.display_reaktoro_state_function = display_reaktoro_state_function
         self.build_input_constraints()
         self.build_output_constraints()
 
@@ -157,31 +176,37 @@ def input_constraints(fs, key):
                         ].get_pyomo_with_required_units()
                     )
 
+    def build_output_vars(self):
+        new_output_vars = {}
+
+        for key, obj in self.solver.output_specs.user_outputs.items():
+            # NOTE: We do not set rkt_outputs to reaktoro_model outputs as they
+            # same as user inputs - we want RKt model to update "user provided vars"
+            # rather then pyomo vars in reaktoro model (e.g. reaktor_block.outputs)
+            if obj.get_pyomo_var() is None:
+                new_output_vars[key] = obj
+        if new_output_vars != {}:
+            self.block.outputs = Var(new_output_vars.keys(), initialize=1)
+            for key, obj in new_output_vars.items():
+                obj.set_pyomo_var(self.block.outputs[key])
+        self.new_output_vars = new_output_vars
+
     def build_output_constraints(self):
         """first update rktOuptutObjects for pyomoBuildProperties with reaktoro pyomo variables as
         they will be used in construction of constraints
         The is will also check if user provided an output pyomo var and if not will
         add them to new_output_var dict, which will be used to create new output variables on the block
         """
-        new_output_vars = {}
         for key, obj in self.solver.output_specs.user_outputs.items():
             if PropTypes.pyomo_built_prop == obj.property_type:
                 for (
                     pyoPropKey,
                     pyoPropObj,
                 ) in obj.pyomo_build_options.properties.items():
-                    pyoPropObj.set_pyomo_var(
-                        self.block.reaktoro_model.outputs[pyoPropKey]
-                    )
-            # NOTE: We do not set rkt_outputs to reaktoro_model outputs as they
-            # same as user inputs - we want RKt model to update "user provided vars"
-            # rather then pyomo vars in reaktoro model (e.g. reaktor_block.outputs)
-            if obj.get_pyomo_var() is None:
-                new_output_vars[key] = obj
-        if new_output_vars != {}:
-            self.block.outputs = Var(new_output_vars.keys(), initialize=1)
-            for key, obj in new_output_vars.items():
-                obj.set_pyomo_var(self.block.outputs[key])
+                    if pyoPropObj.get_pyomo_var() is None:
+                        pyoPropObj.set_pyomo_var(
+                            self.block.reaktoro_model.outputs[pyoPropKey]
+                        )
 
         @self.block.Constraint(self.solver.output_specs.user_outputs)
         def output_constraints(fs, prop, prop_index):
@@ -196,10 +221,14 @@ def output_constraints(fs, prop, prop_index):
                     == self.block.reaktoro_model.outputs[(prop, prop_index)]
                 )
 
-    def initialize(self, presolveDuringInitialization=False):
+    def initialize(self, presolve_during_initialization=False):
         self.initialize_input_variables_and_constraints()
-        self.solver.state.equilibrate_state()
-        self.solver.solve_reaktoro_block(presolve=presolveDuringInitialization)
+
+        if self.reaktoro_initialize_function is None:
+            self.solver.state.equilibrate_state()
+            self.solver.solve_reaktoro_block(presolve=presolve_during_initialization)
+        else:
+            self.reaktoro_initialize_function(presolve=presolve_during_initialization)
         self.initialize_output_variables_and_constraints()
         _log.info(f"Initialized rkt block")
 
@@ -219,9 +248,9 @@ def calc_scale(value):
 
             sf = calc_scale(abs(pyo_var.value))
             if sf > 1e16:
-                _log.warning(f"Var {pyo_var} scale >1e16")
+                _log.warning(f"Var {pyo_var} scale {sf}>1e16")
             if sf < 1e-16:
-                _log.warning(f"Var {pyo_var} scale <1e-16")
+                _log.warning(f"Var {pyo_var} scale {sf}<1e-16")
             return sf
 
     def initialize_output_variables_and_constraints(self):
@@ -267,6 +296,7 @@ def initialize_output_variables_and_constraints(self):
 
         # update jacobian scaling
         self.get_jacobian_scaling()
+        self.set_user_jacobian_scaling()
 
     def get_jacobian_scaling(self):
         if self.jacobian_scaling_type == JacScalingTypes.no_scaling:
@@ -284,12 +314,11 @@ def get_jacobian_scaling(self):
             self.solver.jacobian_scaling_values = (
                 np.sum(np.abs(self.solver.jacobian_matrix) ** 2, axis=1) ** 0.5
             )
-        self.set_user_jacobian_scaling()
 
     def set_user_jacobian_scaling(self, user_scaling=None):
+
         if user_scaling is None:
             user_scaling = self.user_scaling
-
         for i, (key, obj) in enumerate(self.solver.output_specs.rkt_outputs.items()):
             if user_scaling.get(key) != None:
                 scale = user_scaling[key]
@@ -329,3 +358,9 @@ def initialize_input_variables_and_constraints(self):
                 sf = self.get_sf(self.block.reaktoro_model.inputs[key])
             iscale.set_scaling_factor(self.block.reaktoro_model.inputs[key], sf)
             iscale.constraint_scaling_transform(self.block.input_constraints[key], sf)
+
+    def display_state(self):
+        if self.display_reaktoro_state_function is None:
+            print(self.solver.state)
+        else:
+            self.display_reaktoro_state_function()

src/reaktoro_pse/core/reaktoro_gray_box.py

@@ -34,13 +34,32 @@ class HessTypes:
 class ReaktoroGrayBox(ExternalGreyBoxModel):
     ########################################################################################
     # custom Grey Box functions
-    def configure(self, reaktoro_solver):
+    def configure(
+        self,
+        reaktoro_solver=None,
+        inputs=None,
+        input_dict=None,
+        outputs=None,
+        hessian_type=None,
+    ):
         # assign a Reaktoro state object to instance
         self.reaktoro_solver = reaktoro_solver
-        self.hess_type = reaktoro_solver.hessian_type
-        self.inputs = reaktoro_solver.input_specs.rkt_inputs.rkt_input_list
-        self.input_dict = reaktoro_solver.input_specs.rkt_inputs
-        self.outputs = list(reaktoro_solver.output_specs.rkt_outputs.keys())
+        if hessian_type is None:
+            self.hess_type = reaktoro_solver.hessian_type
+        else:
+            self.hess_type = hessian_type
+        if inputs is None:
+            self.inputs = reaktoro_solver.input_specs.rkt_inputs.rkt_input_list
+        else:
+            self.inputs = inputs
+        if input_dict is None:
+            self.input_dict = reaktoro_solver.input_specs.rkt_inputs
+        else:
+            self.input_dict = input_dict
+        if outputs is None:
+            self.outputs = list(reaktoro_solver.output_specs.rkt_outputs.keys())
+        else:
+            self.outputs = outputs
         self._outputs_dual_multipliers = np.ones(len(self.outputs))
         self._hess = np.zeros((len(self.inputs), len(self.inputs)))
         self.header_saved = False
@@ -96,7 +115,6 @@ def get_last_output(self, new_params):
                 self.reaktoro_solver.solve_reaktoro_block(params=new_params)
             )
             self.step += 1
-
         self.old_params = copy.deepcopy(new_params)
 
     def evaluate_jacobian_outputs(self):
@@ -112,7 +130,7 @@ def set_output_constraint_multipliers(self, _outputs_dual_multipliers):
         np.copyto(self._outputs_dual_multipliers, _outputs_dual_multipliers)
 
     def get_output_constraint_scaling_factors(self):
-        return self.reaktoro_solver.jacobian_scaling_values
+        return self.reaktoro_solver.get_jacobian_scaling()
 
     def hessian_gauss_newton_version(self, sparse_jac, threshold=7):
 

src/reaktoro_pse/core/reaktoro_inputs.py

@@ -14,6 +14,7 @@
 from reaktoro_pse.core.reaktoro_state import ReaktoroState
 
 import idaes.logger as idaeslog
+import copy
 
 _log = idaeslog.getLogger(__name__)
 
@@ -22,6 +23,44 @@
 """ class to setup input constraints, and specs for reaktoro solver class"""
 
 
+class ReaktoroInputExport:
+    def __init__(self):
+        self.ignore_elements_for_constraints = []
+        self.fixed_solvent_specie = {}
+        self.fixed_solvent_speciation = {}
+        self.rkt_chemical_inputs = None
+        self.assert_charge_neutrality = None
+        self.neutrality_ion = None
+        self.dissolve_species_in_rkt = None
+        self.exact_speciation = None
+
+    def copy_chem_inputs(self, chem_inputs):
+        self.rkt_chemical_inputs = RktInputs()
+        for key, obj in chem_inputs.items():
+            obj.update_values(True)
+            self.rkt_chemical_inputs[key] = None
+            self.rkt_chemical_inputs[key].time_unit = obj.time_unit
+            self.rkt_chemical_inputs[key].main_unit = obj.main_unit
+            self.rkt_chemical_inputs[key].conversion_unit = obj.conversion_unit
+            self.rkt_chemical_inputs[key].conversion_value = obj.conversion_value
+            self.rkt_chemical_inputs[key].required_unit = obj.required_unit
+            self.rkt_chemical_inputs[key].lower_bound = obj.lower_bound
+            self.rkt_chemical_inputs[key].input_type = obj.input_type
+            self.rkt_chemical_inputs[key].value = obj.value
+            self.rkt_chemical_inputs[key].converted_value = obj.converted_value
+        self.rkt_chemical_inputs.registered_phases = chem_inputs.registered_phases
+        self.rkt_chemical_inputs.all_species = chem_inputs.all_species
+        self.rkt_chemical_inputs.species_list = chem_inputs.species_list
+        self.rkt_chemical_inputs.convert_to_rkt_species = (
+            chem_inputs.convert_to_rkt_species
+        )
+        self.rkt_chemical_inputs.composition_is_elements = (
+            chem_inputs.composition_is_elements
+        )
+        self.rkt_chemical_inputs.conversion_method = chem_inputs.conversion_method
+        self.rkt_chemical_inputs.rkt_input_list = chem_inputs.rkt_input_list
+
+
 class ReaktoroInputSpec:
     def __init__(self, reaktor_state):
         # initialize parameters needed to build reaktor solver
@@ -91,7 +130,6 @@ def register_fixed_solvent_specie(self, phase, specie):
         system speciation, so if we want to specify pH, we need to allow system to find eq. H/O and fix
         H2O"""
         self.fixed_solvent_specie[phase] = specie
-
         self.fixed_solvent_speciation[phase] = {}
 
     def register_free_elements(self, elements):
@@ -118,14 +156,16 @@ def configure_specs(
         """
         self.dissolve_species_in_rkt = dissolve_species_in_rkt
         self.exact_speciation = exact_speciation
+
+    def build_input_specs(self):
+        """function to build all the input specs"""
         self.breakdown_species_to_elements()
         self.equilibrium_specs = rkt.EquilibriumSpecs(self.state.state.system())
         self.add_specs(
             self.equilibrium_specs,
             self.assert_charge_neutrality,
-            dissolve_species_in_rkt,
+            self.dissolve_species_in_rkt,
         )
-
         # get input name order!
         for idx, spec in enumerate(self.equilibrium_specs.namesInputs()):
             if spec == "T":
@@ -436,3 +476,29 @@ def write_empty_constraints(self, spec_object):
         for specie in self.empty_constraints:
             spec_object.openTo(specie)
             self.write_empty_con(spec_object, specie)
+
+    def export_config(self):
+        export_object = ReaktoroInputExport()
+        export_object.copy_chem_inputs(self.rkt_chemical_inputs)
+        export_object.ignore_elements_for_constraints = (
+            self.ignore_elements_for_constraints
+        )
+        export_object.fixed_solvent_specie = self.fixed_solvent_specie
+        export_object.fixed_solvent_speciation = self.fixed_solvent_speciation
+        export_object.assert_charge_neutrality = self.assert_charge_neutrality
+        export_object.neutrality_ion = self.neutrality_ion
+        export_object.dissolve_species_in_rkt = self.dissolve_species_in_rkt
+        export_object.exact_speciation = self.exact_speciation
+        return export_object
+
+    def load_from_export_object(self, export_object):
+        self.ignore_elements_for_constraints = (
+            export_object.ignore_elements_for_constraints
+        )
+        self.fixed_solvent_specie = export_object.fixed_solvent_specie
+        self.fixed_solvent_speciation = export_object.fixed_solvent_speciation
+        self.rkt_chemical_inputs = export_object.rkt_chemical_inputs
+        self.assert_charge_neutrality = export_object.assert_charge_neutrality
+        self.neutrality_ion = export_object.neutrality_ion
+        self.dissolve_species_in_rkt = export_object.dissolve_species_in_rkt
+        self.exact_speciation = export_object.exact_speciation

src/reaktoro_pse/core/reaktoro_jacobian.py

@@ -175,6 +175,13 @@ def display_available(self):
         return available_dict
 
 
+class ReaktoroJacobianExport:
+    def __init__(self):
+        self.der_step_size = None
+        self.jacobian_type = None
+        self.numerical_order = None
+
+
 class ReaktoroJacobianSpec:
     def __init__(self, reaktor_state, reaktor_outputs):
         self.state = reaktor_state
@@ -188,6 +195,20 @@ def __init__(self, reaktor_state, reaktor_outputs):
         self.configure_numerical_jacobian()
         self.check_existing_jacobian_props()
 
+    def export_config(self):
+        export_object = ReaktoroJacobianExport()
+        export_object.der_step_size = self.der_step_size
+        export_object.jacobian_type = self.jacobian_type
+        export_object.numerical_order = self.numerical_order
+        return export_object
+
+    def load_from_export_object(self, export_object):
+        self.configure_numerical_jacobian(
+            jacobian_type=export_object.jacobian_type,
+            order=export_object.numerical_order,
+            step_size=export_object.der_step_size,
+        )
+
     def configure_numerical_jacobian(
         self, jacobian_type="average", order=4, step_size=1e-8
     ):
@@ -200,6 +221,7 @@ def configure_numerical_jacobian(
         """
         self.der_step_size = step_size
         self.jacobian_type = JacType.average
+        self.numerical_order = order
         if jacobian_type == JacType.average:
             self.jacobian_type = JacType.average
             assert order % 2 == 0