PETSc Time-stepping Solver – Chemical Akzo Nobel Example

PETSc Time-stepping Solver – Chemical Akzo Nobel Example#

Author: John Eslick
Maintainer: John Eslick
Updated: 2023-06-01

This example provides an overview of the PETSc time-stepping solver utilities in IDAES, which can be used to solve systems of differential algebraic equations (DAEs). PETSc is a solver suite developed primarily by Argonne National Lab (https://petsc.org/release/). IDAES provides a wrapper for PETSc (IDAES/idaes-ext) that uses the AMPL solver interface (https://ampl.com/resources/learn-more/hooking-your-solver-to-ampl/) and utility functions that allow Pyomo and Pyomo.DAE (https://pyomo.readthedocs.io/en/stable/modeling_extensions/dae.html) problems to be solved using PETSc.

This demonstration problem describes a set of chemical reactions in a reactor. A full description of the problem is available at https://archimede.dm.uniba.it/~testset/report/chemakzo.pdf. This is part of a test set which can be found at https://archimede.uniba.it/~testset/.

Prerequisites#

The PETSc solver is an extra download for IDAES, which can be downloaded using the command idaes get-extensions --extra petsc, if it is not installed already. See the IDAES solver documentation for more information (https://idaes-pse.readthedocs.io/en/stable/reference_guides/core/solvers.html).

Imports#

Import the modules that will be used. Numpy and matplotlib are used to make some plots, idaes.core.solvers.petsc contains the PETSc utilities, and idaes.core.solvers.features contains the example model used here.

import numpy as np
import matplotlib.pyplot as plt

import pyomo.environ as pyo
import idaes.core.solvers.petsc as petsc  # petsc utilities module
from idaes.core.solvers.features import dae  # DAE example/test problem

Set Up the Model#

The model in this example is used for basic solver testing, so it is provided as part of an IDAES solver testing module. The model implementation is standard Pyomo.DAE, and nothing special needs to be done in the model to use the PETSc solver. The IDAES utilities for the PETSc solver will take the discretized Pyomo model and integrate between discrete time points to fill in the solution. To integrate over the entire time domain (as we will do here), you can discretize time using one time element, in which case, the problem will just contain the initial and final points. The intermediate solutions can be read from the trajectory data saved by the solver. The trajectory data can be used for analysis, or interpolation can be used to initialize a Pyomo problem before solving the fully time discretized problem. Integrating over the entire time domain is fastest with a coarsely discretized model (ideally just a single finite element in time) because the model is smaller and there are fewer calls to the integrator. This can be a good way to start testing a new dynamic IDAES model.

# To see the example problem code, uncomment the line below and execute this cell.
# ??dae

# Get the model and known solution for y variables at t=180 minutes.
m, y1, y2, y3, y4, y5, y6 = dae(nfe=10)

The variables y1 to y6 represent concentrations of chemical species. The values returned by the function above are the correct solution at t = 180. These values can be used to verify the solver results. The Pyomo model is m. We are mainly interested in the y variables. The variables y1 to y5 are differential variables and y6 is an algebraic variable. Initial conditions are required for y1 through y5, and the initial values of the other variables can be calculated from there. The variables y1 through y5 at t = 0 are: m.y[0, 1] to m.y[0, 5] and y6 is m.y6[0]. The variables at the final state are m.y[180, 1] to m.y[180, 5] and m.y6[180]. The variable y6 is indexed differently because we want to treat it differently than the differential variables.

# See the initial conditions:
print("at t = 0:")
print(f"    y1 = {pyo.value(m.y[0, 1])}")
print(f"    y2 = {pyo.value(m.y[0, 2])}")
print(f"    y3 = {pyo.value(m.y[0, 3])}")
print(f"    y4 = {pyo.value(m.y[0, 4])}")
print(f"    y5 = {pyo.value(m.y[0, 5])}")

at t = 0:
    y1 = 0.444
    y2 = 0.00123
    y3 = 0.0
    y4 = 0.007
    y5 = 0.0

Solve#

The petsc_dae_by_time_element() function is used to solve Pyomo.DAE discretized Pyomo problem with the PETSc time-stepping solver by integrating between discrete time points. In this case there is only one time element.

# To see the docs, uncomment the line below and execute this cell.
# ?petsc.petsc_dae_by_time_element

# The command below will solve the problem.  In this case, we want to read the saved
# trajectory for each time element in the Pyomo.DAE problem (in this case there is
# only 1) so we will need to provide solver options to save the trajectory to the PETSc
# solver, a file name stub for variable information files, and a file stub for saving
# the trajectory information.  The options shown below will delete the trajectory
# information written by PETSc and resave it as json.  This allows us to cleanly read
# the trajectory data for multiple time elements.

result = petsc.petsc_dae_by_time_element(
    m,
    time=m.t,
    between=[m.t.first(), m.t.last()],
    ts_options={
        "--ts_type": "cn",  # Crank–Nicolson
        "--ts_adapt_type": "basic",
        "--ts_dt": 0.01,
        "--ts_save_trajectory": 1,
    },
)
tj = result.trajectory
res = result.results

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[7], line 9
# The command below will solve the problem.  In this case, we want to read the saved
# trajectory for each time element in the Pyomo.DAE problem (in this case there is
# only 1) so we will need to provide solver options to save the trajectory to the PETSc
   (...)
# information written by PETSc and resave it as json.  This allows us to cleanly read
# the trajectory data for multiple time elements.
----> 9 result = petsc.petsc_dae_by_time_element(
   m,
   time=m.t,
   between=[m.t.first(), m.t.last()],
   ts_options={
       "--ts_type": "cn",  # Crank–Nicolson
       "--ts_adapt_type": "basic",
       "--ts_dt": 0.01,
       "--ts_save_trajectory": 1,
   },
)
tj = result.trajectory
res = result.results

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/solvers/petsc.py:592, in petsc_dae_by_time_element(m, time, timevar, initial_constraints, initial_variables, detect_initial, skip_initial, initial_solver, initial_solver_options, ts_options, keepfiles, symbolic_solver_labels, between, interpolate, calculate_derivatives, previous_trajectory, representative_time, snes_options)
           _sub_problem_scaling_suffix(m, t_block)
           with idaeslog.solver_log(solve_log, idaeslog.INFO) as slc:
--> 592                 res = initial_solver_obj.solve(
                   t_block,
                   tee=slc.tee,
                   symbolic_solver_labels=symbolic_solver_labels,
               )
           res_list.append(res)
tprev = t0

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/opt/base/solvers.py:534, in OptSolver.solve(self, *args, **kwds)
def solve(self, *args, **kwds):
   """Solve the problem"""
--> 534     self.available(exception_flag=True)
   #
   # If the inputs are models, then validate that they have been
   # constructed! Collect suffix names to try and import from solution.
   #
   from pyomo.core.base.block import BlockData

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/solvers/plugins/solvers/ASL.py:118, in ASL.available(self, exception_flag)
def available(self, exception_flag=True):
--> 118     if not super().available(exception_flag):
       return False
   return self.version() is not None

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/opt/solver/shellcmd.py:134, in SystemCallSolver.available(self, exception_flag)
   cm = nullcontext() if exception_flag else LoggingIntercept()
   with cm:
--> 134         ans = self.executable()
except NotImplementedError:
   ans = None

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/opt/solver/shellcmd.py:205, in SystemCallSolver.executable(self)
def executable(self):
   """
   Returns the executable used by this solver.
   """
   return (
       self._user_executable
       if (self._user_executable is not None)
--> 205         else self._default_executable()
   )

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/solvers/petsc.py:124, in Petsc._default_executable(self)
executable = Executable("petsc")
if not executable:
--> 124     raise RuntimeError("No PETSc executable found.")
return executable.path()

RuntimeError: No PETSc executable found.

# Verify results
assert abs(y1 - pyo.value(m.y[180, 1])) / y1 < 1e-3
assert abs(y2 - pyo.value(m.y[180, 2])) / y2 < 1e-3
assert abs(y3 - pyo.value(m.y[180, 3])) / y3 < 1e-3
assert abs(y4 - pyo.value(m.y[180, 4])) / y4 < 1e-3
assert abs(y5 - pyo.value(m.y[180, 5])) / y5 < 1e-3
assert abs(y6 - pyo.value(m.y6[180])) / y6 < 1e-3

Plot Results Stored in Model#

In the problem above we used the PETSc to integrate between the first and last time point. By default, the PETSc DAE utility will interpolate from the PETSc trajectory to fill in the skipped points.

a = plt.plot(m.t, [pyo.value(m.y[t, 1]) for t in m.t], label="y1")
a = plt.plot(m.t, [pyo.value(m.y[t, 2]) for t in m.t], label="y2")
a = plt.plot(m.t, [pyo.value(m.y[t, 3]) for t in m.t], label="y3")
a = plt.plot(m.t, [pyo.value(m.y[t, 4]) for t in m.t], label="y4")
a = plt.plot(m.t, [pyo.value(m.y[t, 5]) for t in m.t], label="y5")
a = plt.plot(m.t, [pyo.value(m.y6[t]) for t in m.t], label="y6")
a = plt.legend()
a = plt.ylabel("Concentration (mol/l)")
a = plt.xlabel("time (min)")

../../_images/cf41289cc14910af14042926c9b78e36a4e7734d56934a9dc0ea3823007bd400.png

Plot Trajectory#

# First plot all y's on one plot

a = plt.plot(tj.time, tj.get_vec(m.y[180, 1]), label="y1")
a = plt.plot(tj.time, tj.get_vec(m.y[180, 2]), label="y2")
a = plt.plot(tj.time, tj.get_vec(m.y[180, 3]), label="y3")
a = plt.plot(tj.time, tj.get_vec(m.y[180, 4]), label="y4")
a = plt.plot(tj.time, tj.get_vec(m.y[180, 5]), label="y5")
a = plt.plot(tj.time, tj.get_vec(m.y6[180]), label="y6")
a = plt.legend()
a = plt.ylabel("Concentration (mol/l)")
a = plt.xlabel("time (min)")

../../_images/52392e1df1e6de9a12e4b4dd6de2691e84572a5ab98e4dbdfeb8169f04f9bc56.png

# 2 and 4 are pretty low concentration, so plot those so we can see better
a = plt.plot(tj.time, tj.get_vec(m.y[180, 2]), label="y2")
a = plt.plot(tj.time, tj.get_vec(m.y[180, 4]), label="y4")
a = plt.legend()
a = plt.ylabel("Concentration (mol/l)")
a = plt.xlabel("time (min)")

../../_images/4a7d4798f4ecf77c265540f437819ef1038494e5dcaedbc3149a826429bff80f.png

# 2 seems to have some fast dynamics so plot a shorter time
a = plt.plot(tj.vecs["_time"], tj.vecs[str(m.y[180, 2])], label="y2")
a = plt.legend()
a = plt.ylabel("Concentration (mol/l)")
a = plt.xlabel("time (min)")
a = plt.xlim(0, 2)

../../_images/7e7b3002472f393e537eeec90d04eac0ae406fe73f16dbbfbd6f4930ae87395f.png

Interpolate Trajectory#

For a number of reasons, such as initializing Pyomo problems or showing results at even time intervals it can be useful to show values at specific time points. The PetscTrajectory class has a method to use linear interpolation to produce a new dictionary of trajectory data at specified time points.

# This creates a new trajectory data set with data every minute.
tji = tj.interpolate(np.linspace(0, 180, 181))

# The plot of this new data should look the same as the original, although some of the
# fast dynamics of component 2 will be obscured.

a = plt.plot(tji.time, tji.get_vec(m.y[180, 1]), label="y1")
a = plt.plot(tji.time, tji.get_vec(m.y[180, 2]), label="y2")
a = plt.plot(tji.time, tji.get_vec(m.y[180, 3]), label="y3")
a = plt.plot(tji.time, tji.get_vec(m.y[180, 4]), label="y4")
a = plt.plot(tji.time, tji.get_vec(m.y[180, 5]), label="y5")
a = plt.plot(tji.time, tji.get_vec(m.y6[180]), label="y6")
a = plt.legend()
a = plt.ylabel("Concentration (mol/l)")
a = plt.xlabel("time (min)")

../../_images/8f589b5bd60f7570c79f8fcf1458f6bbe287cfbd2414280b60ad9618158456cd.png

a = plt.plot(tji.time, tji.get_vec(m.y[180, 2]), label="y2 interpolate dt=1")
a = plt.plot(tj.time, tj.get_vec(m.y[180, 2]), label="y2 original")
a = plt.legend()
a = plt.ylabel("Concentration (mol/l)")
a = plt.xlabel("time (min)")
a = plt.xlim(0, 2)

../../_images/84d64c5858f31a304c8e46d047898261efc1acdddca082e3f9bd06aeebe79fbf.png