diff --git a/idaes/core/solvers/petsc.py b/idaes/core/solvers/petsc.py index 7129ccbee5..8cff031efb 100644 --- a/idaes/core/solvers/petsc.py +++ b/idaes/core/solvers/petsc.py @@ -42,7 +42,10 @@ import idaes import idaes.logger as idaeslog import idaes.config as icfg -from idaes.core.util.model_statistics import degrees_of_freedom, number_activated_constraints +from idaes.core.util.model_statistics import ( + degrees_of_freedom, + number_activated_constraints, +) # Importing a few things here so that they are cached @@ -252,7 +255,7 @@ def _copy_time(time_vars, t_from, t_to): def find_discretization_equations(m, time): """ This returns a list of continuity equations (if Lagrange-Legendre collocation is used) - and time discretization equations. Since we aren't solving the whole time period + and time discretization equations. Since we aren't solving the whole time period simultaneously, we'll want to deactivate these constraints. Args: @@ -290,6 +293,7 @@ def find_discretization_equations(m, time): return disc_eqns + def get_discretization_constraint_data(m, time): """ This returns a list of continuity equations (if Lagrange-Legendre collocation is used) @@ -567,10 +571,8 @@ def petsc_dae_by_time_element( elif representative_time not in between: raise RuntimeError("representative_time is not element of between.") - flattened_problem = _get_flattened_problem( - model=m, - time=time, - representative_time=representative_time + flattened_problem = _get_flattened_problem( + model=m, time=time, representative_time=representative_time ) solver_dae = pyo.SolverFactory("petsc_ts", options=ts_options) @@ -604,10 +606,13 @@ def petsc_dae_by_time_element( ) except RuntimeError as err: if "Zero constraints" in err.message: - init_subsystem=None + init_subsystem = None else: raise err - if init_subsystem is not None and number_activated_constraints(init_subsystem) > 0: + if ( + init_subsystem is not None + and number_activated_constraints(init_subsystem) > 0 + ): dof = degrees_of_freedom(init_subsystem) if dof > 0: raise RuntimeError( @@ -639,7 +644,9 @@ def petsc_dae_by_time_element( if t == between.first(): # t == between.first() was handled above continue - constraints = [con[t] for con in flattened_problem["time_constraints"] if t in con] + constraints = [ + con[t] for con in flattened_problem["time_constraints"] if t in con + ] variables = [var[t] for var in flattened_problem["time_variables"]] # Create a temporary block with references to original constraints # and variables so we can integrate this "subsystem" without @@ -1057,6 +1064,7 @@ def from_json(self, pth): self.vecs = json.load(fp) self.time = self.vecs["_time"] + def _get_flattened_problem(model, time, representative_time): """ Helper function for petsc_dae_by_time_element and get_initial_condition_problem. @@ -1067,11 +1075,11 @@ def _get_flattened_problem(model, time, representative_time): time (ContinuousSet): Time set representative_time (Element of time): A timepoint at which the DAE system is at its "normal" state after the constraints and variables associated with the initial time point - have passed. + have passed. Returns (dictionary): Dictionary containing lists of variables and constraints indexed by times and lists - of those unindexed by time. Also a list of discretization equations so they can be + of those unindexed by time. Also a list of discretization equations so they can be deactivated in the problem passed to PETSc. """ regular_vars, time_vars = flatten_dae_components( @@ -1092,15 +1100,15 @@ def _get_flattened_problem(model, time, representative_time): def get_initial_condition_problem( - model, - time, - initial_time, - representative_time=None, - initial_constraints=None, - initial_variables=None, - detect_initial=True, - flattened_problem=None, - ): + model, + time, + initial_time, + representative_time=None, + initial_constraints=None, + initial_variables=None, + detect_initial=True, + flattened_problem=None, +): """ Solve a DAE problem step by step using the PETSc DAE solver. This integrates from one time point to the next. @@ -1139,11 +1147,11 @@ def get_initial_condition_problem( if flattened_problem is None: if representative_time is None: - raise RuntimeError("The user must supply either the flattened problem or a representative time.") + raise RuntimeError( + "The user must supply either the flattened problem or a representative time." + ) flattened_problem = _get_flattened_problem( - model=model, - time=time, - representative_time=representative_time + model=model, time=time, representative_time=representative_time ) if detect_initial: @@ -1156,17 +1164,14 @@ def get_initial_condition_problem( const_init_set = ComponentSet(initial_constraints) initial_constraints = list(const_no_t_set | const_init_set) - - constraints = [ - con[initial_time] + con[initial_time] for con in flattened_problem["time_constraints"] if initial_time in con and con[initial_time] not in flattened_problem["discretization_equations"] ] + initial_constraints variables = [ - var[initial_time] - for var in flattened_problem["time_variables"] + var[initial_time] for var in flattened_problem["time_variables"] ] + initial_variables if len(constraints) <= 0: @@ -1181,4 +1186,4 @@ def get_initial_condition_problem( variables, ) _sub_problem_scaling_suffix(model, subsystem) - return subsystem \ No newline at end of file + return subsystem diff --git a/idaes/core/solvers/tests/test_petsc.py b/idaes/core/solvers/tests/test_petsc.py index 846ef067c4..13ac59531e 100644 --- a/idaes/core/solvers/tests/test_petsc.py +++ b/idaes/core/solvers/tests/test_petsc.py @@ -812,6 +812,7 @@ def test_petsc_skip_initial_solve(): assert pytest.approx(y5, rel=1e-3) == pyo.value(m.y[m.t.last(), 5]) assert pytest.approx(y6, rel=1e-3) == pyo.value(m.y[m.t.last(), 6]) + @pytest.mark.unit @pytest.mark.skipif(not petsc.petsc_available(), reason="PETSc solver not available") def test_petsc_collocation_lagrange_radau(): @@ -819,9 +820,7 @@ def test_petsc_collocation_lagrange_radau(): Ensure that PETSc works with Lagrange-Radau collocation """ m, y1, y2, y3, y4, y5, y6 = dae_with_non_time_indexed_constraint( - transformation_method="dae.collocation", - scheme="LAGRANGE-RADAU", - ncp=3 + transformation_method="dae.collocation", scheme="LAGRANGE-RADAU", ncp=3 ) res = petsc.petsc_dae_by_time_element( @@ -843,6 +842,7 @@ def test_petsc_collocation_lagrange_radau(): assert pytest.approx(y5, rel=1e-3) == pyo.value(m.y[m.t.last(), 5]) assert pytest.approx(y6, rel=1e-3) == pyo.value(m.y[m.t.last(), 6]) + @pytest.mark.unit @pytest.mark.skipif(not petsc.petsc_available(), reason="PETSc solver not available") def test_petsc_collocation_lagrange_legendre(): @@ -850,9 +850,7 @@ def test_petsc_collocation_lagrange_legendre(): Ensure that PETSc works with Lagrange-Legendre collocation """ m, y1, y2, y3, y4, y5, y6 = dae_with_non_time_indexed_constraint( - transformation_method="dae.collocation", - scheme="LAGRANGE-LEGENDRE", - ncp=3 + transformation_method="dae.collocation", scheme="LAGRANGE-LEGENDRE", ncp=3 ) res = petsc.petsc_dae_by_time_element( @@ -932,11 +930,11 @@ def test_petsc_traj_previous(): initial_solver="ipopt", initial_solver_options={ # 'bound_push' : 1e-22, - 'linear_solver': 'ma57', - 'OF_ma57_automatic_scaling': 'yes', - 'max_iter': 300, - 'tol': 1e-6, - 'halt_on_ampl_error': 'yes', + "linear_solver": "ma57", + "OF_ma57_automatic_scaling": "yes", + "max_iter": 300, + "tol": 1e-6, + "halt_on_ampl_error": "yes", # 'mu_strategy': 'monotone', }, # skip_initial=True, @@ -1353,6 +1351,7 @@ def diff_eq_rule(m, t): calculate_derivatives=True, ) + @pytest.mark.unit def test_get_initial_condition_problem_no_active_constraints(): m = pyo.ConcreteModel() @@ -1378,16 +1377,21 @@ def diff_eq_rule(m, t): model=m, time=m.time, initial_time=m.time.at(1), - representative_time=m.time.at(2) + representative_time=m.time.at(2), ) assert mstat.number_variables_in_activated_constraints(init_subsystem) == 0 - local_var_name_list = [var.name for var in init_subsystem.component_data_objects(ctype=pyo.Var)] - assert set(local_var_name_list) == {'x[0.0]', 'u[0.0]', 'dxdt[0.0]'} + local_var_name_list = [ + var.name for var in init_subsystem.component_data_objects(ctype=pyo.Var) + ] + assert set(local_var_name_list) == {"x[0.0]", "u[0.0]", "dxdt[0.0]"} assert mstat.number_activated_constraints(init_subsystem) == 0 - local_con_name_list = [con.name for con in init_subsystem.component_data_objects(ctype=pyo.Constraint)] + local_con_name_list = [ + con.name for con in init_subsystem.component_data_objects(ctype=pyo.Constraint) + ] assert local_con_name_list[0] == "diff_eq[0.0]" + @pytest.mark.unit def test_get_initial_condition_problem_no_active_constraints_not_t0(): m = pyo.ConcreteModel() @@ -1413,16 +1417,21 @@ def diff_eq_rule(m, t): model=m, time=m.time, initial_time=m.time.at(2), - representative_time=m.time.at(3) + representative_time=m.time.at(3), ) assert mstat.number_variables_in_activated_constraints(init_subsystem) == 0 - local_var_name_list = [var.name for var in init_subsystem.component_data_objects(ctype=pyo.Var)] - assert set(local_var_name_list) == {'x[1.0]', 'u[1.0]', 'dxdt[1.0]'} + local_var_name_list = [ + var.name for var in init_subsystem.component_data_objects(ctype=pyo.Var) + ] + assert set(local_var_name_list) == {"x[1.0]", "u[1.0]", "dxdt[1.0]"} assert mstat.number_activated_constraints(init_subsystem) == 0 - local_con_name_list = [con.name for con in init_subsystem.component_data_objects(ctype=pyo.Constraint)] + local_con_name_list = [ + con.name for con in init_subsystem.component_data_objects(ctype=pyo.Constraint) + ] assert local_con_name_list[0] == "diff_eq[1.0]" + @pytest.mark.unit def test_get_initial_condition_problem_at_time_no_active_constraints(): m = pyo.ConcreteModel() @@ -1445,34 +1454,103 @@ def diff_eq_rule(m, t): model=m, time=m.time, initial_time=m.time.at(1), - representative_time=m.time.at(2) + representative_time=m.time.at(2), ) assert mstat.number_variables_in_activated_constraints(init_subsystem) == 0 - local_var_name_list = [var.name for var in init_subsystem.component_data_objects(ctype=pyo.Var)] - assert set(local_var_name_list) == {'x[0.0]', 'u[0.0]', 'dxdt[0.0]'} + local_var_name_list = [ + var.name for var in init_subsystem.component_data_objects(ctype=pyo.Var) + ] + assert set(local_var_name_list) == {"x[0.0]", "u[0.0]", "dxdt[0.0]"} assert mstat.number_activated_constraints(init_subsystem) == 0 - local_con_name_list = [con.name for con in init_subsystem.component_data_objects(ctype=pyo.Constraint)] + local_con_name_list = [ + con.name for con in init_subsystem.component_data_objects(ctype=pyo.Constraint) + ] assert local_con_name_list[0] == "diff_eq[0.0]" + @pytest.mark.unit def test_get_initial_condition_problem_non_time_indexed_constraint(): m, y1, y2, y3, y4, y5, y6 = dae_with_non_time_indexed_constraint(nfe=10) init_subsystem = petsc.get_initial_condition_problem( model=m, time=m.t, initial_time=m.t.at(1), representative_time=m.t.at(2) ) - - unfixed_var_name_set = set(var.name for var in mstat.unfixed_variables_set(init_subsystem)) - assert unfixed_var_name_set == {'y[0,6]', 'ydot[0,1]', 'ydot[0,2]', 'ydot[0,3]', 'ydot[0,4]', 'ydot[0,5]', 'r[0,1]', 'r[0,2]', 'r[0,3]', 'r[0,4]', 'r[0,5]', 'Fin[0]', 'H'} - total_var_name_set = set(var.name for var in mstat.variables_set(init_subsystem)) - assert total_var_name_set == {'r[0,1]', 'r[0,2]', 'ydot[0,1]', 'ydot[0,4]', 'H', 'y[0,2]', 'ydot[0,3]', 'y[0,3]', 'y[0,5]', 'r[0,3]', 'r[0,5]', 'Fin[0]', 'y[0,1]', 'ydot[0,2]', 'ydot[0,6]', 'r[0,4]', 'y[0,4]', 'y[0,6]', 'ydot[0,5]'} + unfixed_var_name_set = set( + var.name for var in mstat.unfixed_variables_set(init_subsystem) + ) + assert unfixed_var_name_set == { + "y[0,6]", + "ydot[0,1]", + "ydot[0,2]", + "ydot[0,3]", + "ydot[0,4]", + "ydot[0,5]", + "r[0,1]", + "r[0,2]", + "r[0,3]", + "r[0,4]", + "r[0,5]", + "Fin[0]", + "H", + } - active_con_name_set = set(con.name for con in mstat.activated_constraints_set(init_subsystem)) - assert active_con_name_set == {'eq_y6[0]', 'eq_r1[0]', 'eq_r2[0]', 'eq_r3[0]', 'eq_r4[0]', 'eq_r5[0]', 'eq_Fin[0]', 'H_eqn'} + total_var_name_set = set(var.name for var in mstat.variables_set(init_subsystem)) + assert total_var_name_set == { + "r[0,1]", + "r[0,2]", + "ydot[0,1]", + "ydot[0,4]", + "H", + "y[0,2]", + "ydot[0,3]", + "y[0,3]", + "y[0,5]", + "r[0,3]", + "r[0,5]", + "Fin[0]", + "y[0,1]", + "ydot[0,2]", + "ydot[0,6]", + "r[0,4]", + "y[0,4]", + "y[0,6]", + "ydot[0,5]", + } + + active_con_name_set = set( + con.name for con in mstat.activated_constraints_set(init_subsystem) + ) + assert active_con_name_set == { + "eq_y6[0]", + "eq_r1[0]", + "eq_r2[0]", + "eq_r3[0]", + "eq_r4[0]", + "eq_r5[0]", + "eq_Fin[0]", + "H_eqn", + } + + total_con_name_set = set( + con.name for con in mstat.total_constraints_set(init_subsystem) + ) + assert total_con_name_set == { + "eq_r5[0]", + "eq_y6[0]", + "eq_r4[0]", + "eq_Fin[0]", + "H_eqn", + "eq_ydot5[0]", + "eq_ydot1[0]", + "eq_r1[0]", + "eq_ydot4[0]", + "eq_r3[0]", + "eq_r2[0]", + "eq_ydot2[0]", + "eq_ydot3[0]", + } - total_con_name_set = set(con.name for con in mstat.total_constraints_set(init_subsystem)) - assert total_con_name_set == {'eq_r5[0]', 'eq_y6[0]', 'eq_r4[0]', 'eq_Fin[0]', 'H_eqn', 'eq_ydot5[0]', 'eq_ydot1[0]', 'eq_r1[0]', 'eq_ydot4[0]', 'eq_r3[0]', 'eq_r2[0]', 'eq_ydot2[0]', 'eq_ydot3[0]'} @pytest.mark.unit def test_get_initial_condition_problem_non_time_indexed_constraint_not_t0(): @@ -1489,17 +1567,80 @@ def test_get_initial_condition_problem_non_time_indexed_constraint_not_t0(): m.eq_ydot4[t4].deactivate() m.eq_ydot5[t4].deactivate() - unfixed_var_name_set = set(var.name for var in mstat.unfixed_variables_set(init_subsystem)) - assert unfixed_var_name_set == {'y[54.0,6]', 'ydot[54.0,1]', 'ydot[54.0,2]', 'ydot[54.0,3]', 'ydot[54.0,4]', 'ydot[54.0,5]', 'r[54.0,1]', 'r[54.0,2]', 'r[54.0,3]', 'r[54.0,4]', 'r[54.0,5]', 'Fin[54.0]', 'H'} + unfixed_var_name_set = set( + var.name for var in mstat.unfixed_variables_set(init_subsystem) + ) + assert unfixed_var_name_set == { + "y[54.0,6]", + "ydot[54.0,1]", + "ydot[54.0,2]", + "ydot[54.0,3]", + "ydot[54.0,4]", + "ydot[54.0,5]", + "r[54.0,1]", + "r[54.0,2]", + "r[54.0,3]", + "r[54.0,4]", + "r[54.0,5]", + "Fin[54.0]", + "H", + } total_var_name_set = set(var.name for var in mstat.variables_set(init_subsystem)) - assert total_var_name_set == {'r[54.0,1]', 'r[54.0,2]', 'ydot[54.0,1]', 'ydot[54.0,4]', 'H', 'y[54.0,2]', 'ydot[54.0,3]', 'y[54.0,3]', 'y[54.0,5]', 'r[54.0,3]', 'r[54.0,5]', 'Fin[54.0]', 'y[54.0,1]', 'ydot[54.0,2]', 'ydot[54.0,6]', 'r[54.0,4]', 'y[54.0,4]', 'y[54.0,6]', 'ydot[54.0,5]'} - - active_con_name_set = set(con.name for con in mstat.activated_constraints_set(init_subsystem)) - assert active_con_name_set == {'eq_y6[54.0]', 'eq_r1[54.0]', 'eq_r2[54.0]', 'eq_r3[54.0]', 'eq_r4[54.0]', 'eq_r5[54.0]', 'eq_Fin[54.0]', 'H_eqn'} - - total_con_name_set = set(con.name for con in mstat.total_constraints_set(init_subsystem)) - assert total_con_name_set == {'eq_r5[54.0]', 'eq_y6[54.0]', 'eq_r4[54.0]', 'eq_Fin[54.0]', 'H_eqn', 'eq_ydot5[54.0]', 'eq_ydot1[54.0]', 'eq_r1[54.0]', 'eq_ydot4[54.0]', 'eq_r3[54.0]', 'eq_r2[54.0]', 'eq_ydot2[54.0]', 'eq_ydot3[54.0]'} + assert total_var_name_set == { + "r[54.0,1]", + "r[54.0,2]", + "ydot[54.0,1]", + "ydot[54.0,4]", + "H", + "y[54.0,2]", + "ydot[54.0,3]", + "y[54.0,3]", + "y[54.0,5]", + "r[54.0,3]", + "r[54.0,5]", + "Fin[54.0]", + "y[54.0,1]", + "ydot[54.0,2]", + "ydot[54.0,6]", + "r[54.0,4]", + "y[54.0,4]", + "y[54.0,6]", + "ydot[54.0,5]", + } + + active_con_name_set = set( + con.name for con in mstat.activated_constraints_set(init_subsystem) + ) + assert active_con_name_set == { + "eq_y6[54.0]", + "eq_r1[54.0]", + "eq_r2[54.0]", + "eq_r3[54.0]", + "eq_r4[54.0]", + "eq_r5[54.0]", + "eq_Fin[54.0]", + "H_eqn", + } + + total_con_name_set = set( + con.name for con in mstat.total_constraints_set(init_subsystem) + ) + assert total_con_name_set == { + "eq_r5[54.0]", + "eq_y6[54.0]", + "eq_r4[54.0]", + "eq_Fin[54.0]", + "H_eqn", + "eq_ydot5[54.0]", + "eq_ydot1[54.0]", + "eq_r1[54.0]", + "eq_ydot4[54.0]", + "eq_r3[54.0]", + "eq_r2[54.0]", + "eq_ydot2[54.0]", + "eq_ydot3[54.0]", + } if __name__ == "__main__": @@ -1518,7 +1659,10 @@ def test_get_initial_condition_problem_non_time_indexed_constraint_not_t0(): m.eq_ydot3[t4].deactivate() m.eq_ydot4[t4].deactivate() m.eq_ydot5[t4].deactivate() - init_subsystem = petsc.get_initial_condition_problem(model=m, time=m.t, initial_time=m.t.at(4), representative_time=m.t.at(5)) + init_subsystem = petsc.get_initial_condition_problem( + model=m, time=m.t, initial_time=m.t.at(4), representative_time=m.t.at(5) + ) from idaes.core.util.model_diagnostics import DiagnosticsToolbox - dt=DiagnosticsToolbox(init_subsystem) + + dt = DiagnosticsToolbox(init_subsystem) dt.report_structural_issues()