Adds a parallel manager, allowing parallel execution of reaktoro blocks.

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):
@@ -198,8 +223,12 @@ def output_constraints(fs, prop, prop_index):
 
     def initialize(self, presolve_during_initialization=False):
         self.initialize_input_variables_and_constraints()
-        self.solver.state.equilibrate_state()
-        self.solver.solve_reaktoro_block(presolve=presolve_during_initialization)
+
+        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")
 
@@ -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
@@ -111,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

@@ -35,10 +35,30 @@ def __init__(self):
         self.exact_speciation = None
 
     def copy_chem_inputs(self, chem_inputs):
-        self.rkt_chemical_inputs = copy.deepcopy(chem_inputs)
-        for key, obj in self.rkt_chemical_inputs.items():
-            # self.inputs[key] = copy.deepcopy(obj)
-            self.rkt_chemical_inputs[key].delete_pyomo_var()
+        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:
@@ -146,7 +166,6 @@ def build_input_specs(self):
             self.assert_charge_neutrality,
             self.dissolve_species_in_rkt,
         )
-
         # get input name order!
         for idx, spec in enumerate(self.equilibrium_specs.namesInputs()):
             if spec == "T":

src/reaktoro_pse/core/reaktoro_outputs.py

@@ -89,6 +89,9 @@ def set_option_function(self, property_type, get_function):
     def set_poyomo_build_option(self, func):
         self.pyomo_build_options = func
 
+    def set_pyo_value(self, value):
+        self.pyomo_var.set_value(value)
+
     def set_get_function(self, func):
         self.get_function = func
 

src/reaktoro_pse/core/reaktoro_solver.py

@@ -87,6 +87,7 @@ def __init__(
         self._sequential_fails = 0
         self._max_fails = 30
         self._input_params = {}
+        self.jacobian_scaling_values = None
 
     def export_config(self):
         export_object = ReaktoroSolverExport()
@@ -234,7 +235,9 @@ def solve_reaktoro_block(
                 _log.warning(f"{key}: {value}")
             self._sequential_fails += 1
             if self._sequential_fails > self._max_fails:
-                assert False
+                raise RuntimeError(
+                    "Number of failed solves exceed maximum allowed number"
+                )
             raise cyipopt.CyIpoptEvaluationError
         else:
             self._sequential_fails = 0
@@ -259,3 +262,6 @@ def try_solve(self, presolve=False):
         )
         self.output_specs.update_supported_props()
         return result
+
+    def get_jacobian_scaling(self):
+        return self.jacobian_scaling_values

src/reaktoro_pse/core/reaktoro_state.py

@@ -124,7 +124,7 @@ def set_activity_model(
             activity_model, state_of_matter
         )
 
-    def get_registered_phases(self, activate_database):
+    def get_registered_phases(self, active_database):
         """
         creates listof phases with applied activity models for generation of reaktoro state
         args:
@@ -139,7 +139,7 @@ def get_registered_phases(self, activate_database):
                 or phase_object.non_speciate_phase_list is not None
             ):
                 rkt_phase_object = self.create_rkt_phase(
-                    activate_database,
+                    active_database,
                     phase_object.phase_function,
                     phase_object.phase_list,
                     phase_object.non_speciate_phase_list,
@@ -149,27 +149,29 @@ def get_registered_phases(self, activate_database):
                     for rpo in rkt_phase_object:
                         self.apply_activity_model(
                             rpo,
+                            active_database,
                             phase_object.activity_model,
                             phase_object.state_of_matter,
                         )
                         activate_phase_list.append(rpo)
                 else:
                     self.apply_activity_model(
                         rkt_phase_object,
+                        active_database,
                         phase_object.activity_model,
                         phase_object.state_of_matter,
                     )
                     activate_phase_list.append(rkt_phase_object)
         return activate_phase_list
 
     def apply_activity_model(
-        self, rkt_phase_object, activity_model, state_of_matter=None
+        self, rkt_phase_object, active_database, activity_model, state_of_matter=None
     ):
         """sets activity mode"""
         if activity_model is None:
             raise TypeError(f"Activity model for {rkt_phase_object} is not set.")
         rkt_activity_model_object = self._process_activity(
-            activity_model, activity_model, state_of_matter
+            active_database, activity_model, state_of_matter
         )
         rkt_phase_object.set(rkt_activity_model_object)
 
@@ -235,10 +237,9 @@ def get_func(activity_model, state_of_matter):
                 if state_of_matter is None:
                     try:
                         return getattr(rkt, activity_model)()
-                    except RuntimeError:
+                    except TypeError:
                         # might require database as input (e.g ActivityModelPhreeqc)
                         return getattr(rkt, activity_model)(active_database)
-
                 else:
                     return getattr(rkt, activity_model)(state_of_matter)
             elif isinstance(activity_model, tuple):
@@ -276,10 +277,26 @@ def __init__(self):
         self.database_file = None
 
     def copy_rkt_inputs(self, inputs):
-        self.inputs = copy.deepcopy(inputs)
-        for key, obj in self.inputs.items():
-            # self.inputs[key] = copy.deepcopy(obj)
-            self.inputs[key].delete_pyomo_var()
+        self.inputs = RktInputs.RktInputs()
+        for key, obj in inputs.items():
+            obj.update_values(True)
+            self.inputs[key] = None
+            self.inputs[key].time_unit = obj.time_unit
+            self.inputs[key].main_unit = obj.main_unit
+            self.inputs[key].conversion_unit = obj.conversion_unit
+            self.inputs[key].conversion_value = obj.conversion_value
+            self.inputs[key].required_unit = obj.required_unit
+            self.inputs[key].lower_bound = obj.lower_bound
+            self.inputs[key].input_type = obj.input_type
+            self.inputs[key].value = obj.value
+            self.inputs[key].converted_value = obj.converted_value
+            self.inputs.registered_phases = inputs.registered_phases
+        self.inputs.all_species = inputs.all_species
+        self.inputs.species_list = inputs.species_list
+        self.inputs.convert_to_rkt_species = inputs.convert_to_rkt_species
+        self.inputs.composition_is_elements = inputs.composition_is_elements
+        self.inputs.conversion_method = inputs.conversion_method
+        self.inputs.rkt_input_list = inputs.rkt_input_list
 
 
 # base class for configuring reaktoro states and solver
@@ -792,3 +809,4 @@ def load_from_export_object(self, export_object):
         self.phase_manager = export_object.phase_manager
         self.database_file = export_object.database_file
         self.database_type = export_object.database_type
+        self.exclude_species_list = export_object.exclude_species_list

src/reaktoro_pse/core/tests/test_reaktoro_builder.py

@@ -143,7 +143,6 @@ def test_build_with_rkt_dissolution(build_with_dissolve_in_rkt):
     builder.initialize()
     m.display()
     m.rkt_block.reaktoro_model.display()
-    m.rkt_block.output_constraints.pprint()
     # will have as many DOFs as outputs due to pyomo not
     # knowing tha graybox exists.
     assert len(m.rkt_block.reaktoro_model.inputs) == len(
@@ -175,6 +174,7 @@ def test_build_with_pyomo_dissolution(build_with_dissolve_in_pyomo):
     builder.initialize()
     # will have as many DOFs as outputs due to pyomo not
     # knowing tha graybox exists.
+    print(rkt_solver.output_specs.rkt_outputs)
     assert degrees_of_freedom(m) == len(rkt_solver.output_specs.rkt_outputs)
     cy_solver = get_solver(solver="cyipopt-watertap")
     cy_solver.options["max_iter"] = 20