PETScKrylov

PETScKrylov is an iterative linear solver that wraps the linear solution methods found in PETSc via petsc4py. The default method is “fgmres”, or the Flexible Generalized Minimal RESidual method, though you may choose any of the other methods in PETSc. This linear solver is capable of handling any system topology effectively. It also solves all subsystems below it in the hierarchy, so assigning different solvers to subsystems will have no effect on the solution at this level.

This solver works under MPI, so it is a good alternative to ScipyKrylov. This solver is also re-entrant, so there are no problems if it is nested during preconditioning.

Here, we calculate the total derivatives across the Sellar system.

import numpy as np 
import openmdao.api as om
from openmdao.test_suite.components.sellar import SellarDis1withDerivatives, SellarDis2withDerivatives

prob = om.Problem()
model = prob.model

model.add_subsystem('d1', SellarDis1withDerivatives(), promotes=['x', 'z', 'y1', 'y2'])
model.add_subsystem('d2', SellarDis2withDerivatives(), promotes=['z', 'y1', 'y2'])

model.add_subsystem('obj_cmp', om.ExecComp('obj = x**2 + z[1] + y1 + exp(-y2)',
                                           z=np.array([0.0, 0.0]), x=0.0),
                    promotes=['obj', 'x', 'z', 'y1', 'y2'])

model.add_subsystem('con_cmp1', om.ExecComp('con1 = 3.16 - y1'), promotes=['con1', 'y1'])
model.add_subsystem('con_cmp2', om.ExecComp('con2 = y2 - 24.0'), promotes=['con2', 'y2'])

model.nonlinear_solver = om.NonlinearBlockGS()

model.linear_solver = om.PETScKrylov()

prob.setup()

prob.set_val('x', 1.)
prob.set_val('z', np.array([5.0, 2.0]))

prob.run_model()

wrt = ['z']
of = ['obj']

J = prob.compute_totals(of=of, wrt=wrt, return_format='flat_dict')
NL: NLBGS Converged in 8 iterations
print(J['obj', 'z'][0][0], 9.61001056, .00001)
print(J['obj', 'z'][0][1], 1.78448534, .00001)
9.610010556989947 9.61001056 1e-05
1.7844853356313641 1.78448534 1e-05

PETScKrylov Options

Option Default Acceptable Values Acceptable Types Description
assemble_jac False [True, False] ['bool'] Activates use of assembled jacobian by this solver.
atol 1e-10 N/A N/A absolute error tolerance
err_on_non_convergeFalse [True, False] ['bool'] When True, AnalysisError will be raised if we don't converge.
iprint 1 N/A ['int'] whether to print output
ksp_type fgmres ['richardson', 'chebyshev', 'cg', 'groppcg', 'pipecg', 'pipecgrr', 'cgne', 'nash', 'stcg', 'gltr', 'fcg', 'pipefcg', 'gmres', 'pipefgmres', 'fgmres', 'lgmres', 'dgmres', 'pgmres', 'tcqmr', 'bcgs', 'ibcgs', 'fbcgs', 'fbcgsr', 'bcgsl', 'cgs', 'tfqmr', 'cr', 'pipecr', 'lsqr', 'preonly', 'qcg', 'bicg', 'minres', 'symmlq', 'lcd', 'python', 'gcr', 'pipegcr', 'tsirm', 'cgls']N/A KSP algorithm to use. Default is 'fgmres'.
maxiter 100 N/A ['int'] maximum number of iterations
precon_side right ['left', 'right'] N/A Preconditioner side, default is right.
restart 1000 N/A ['int'] Number of iterations between restarts. Larger values increase iteration cost, but may be necessary for convergence
rtol 1e-10 N/A N/A relative error tolerance

PETScKrylov Constructor

The call signature for the PETScKrylov constructor is:

PETScKrylov.__init__(**kwargs)[source]

Declare the solver options.

PETScKrylov Option Examples

maxiter

maxiter lets you specify the maximum number of GMRES (or other algorithm) iterations to apply. The default maximum is 100, which is much higher than the other linear solvers because each multiplication by the system Jacobian is considered to be an iteration. You may have to decrease this value if you have a coupled system that is converging very slowly. (Of course, in such a case, it may be better to add a preconditioner.) Alternatively, you may have to raise maxiter if you have an extremely large number of components in your system (a 1000-component ring would need 1000 iterations just to make it around once.)

This example shows what happens if you set maxiter too low (the derivatives should be nonzero, but it stops too soon.)

from openmdao.test_suite.components.sellar import SellarDis1withDerivatives, SellarDis2withDerivatives

prob = om.Problem()
model = prob.model

model.add_subsystem('d1', SellarDis1withDerivatives(), promotes=['x', 'z', 'y1', 'y2'])
model.add_subsystem('d2', SellarDis2withDerivatives(), promotes=['z', 'y1', 'y2'])

model.add_subsystem('obj_cmp', om.ExecComp('obj = x**2 + z[1] + y1 + exp(-y2)',
                                           z=np.array([0.0, 0.0]), x=0.0),
                    promotes=['obj', 'x', 'z', 'y1', 'y2'])

model.add_subsystem('con_cmp1', om.ExecComp('con1 = 3.16 - y1'), promotes=['con1', 'y1'])
model.add_subsystem('con_cmp2', om.ExecComp('con2 = y2 - 24.0'), promotes=['con2', 'y2'])

model.nonlinear_solver = om.NonlinearBlockGS()

model.linear_solver = om.PETScKrylov()
model.linear_solver.options['maxiter'] = 3

prob.setup()

prob.set_val('x', 1.)
prob.set_val('z', np.array([5.0, 2.0]))

prob.run_model()

wrt = ['z']
of = ['obj']

J = prob.compute_totals(of=of, wrt=wrt, return_format='flat_dict')
NL: NLBGS Converged in 8 iterations
print(J['obj', 'z'][0][0])
print(J['obj', 'z'][0][1])
4.932180270630899
1.7340645454399009

atol

The absolute convergence tolerance, the absolute size of the (possibly preconditioned) residual norm.

You may need to adjust this setting if you have abnormally large or small values in your global Jacobian.

from openmdao.test_suite.components.sellar import SellarDis1withDerivatives, SellarDis2withDerivatives

prob = om.Problem()
model = prob.model

model.add_subsystem('d1', SellarDis1withDerivatives(), promotes=['x', 'z', 'y1', 'y2'])
model.add_subsystem('d2', SellarDis2withDerivatives(), promotes=['z', 'y1', 'y2'])

model.add_subsystem('obj_cmp', om.ExecComp('obj = x**2 + z[1] + y1 + exp(-y2)',
                                           z=np.array([0.0, 0.0]), x=0.0),
                    promotes=['obj', 'x', 'z', 'y1', 'y2'])

model.add_subsystem('con_cmp1', om.ExecComp('con1 = 3.16 - y1'), promotes=['con1', 'y1'])
model.add_subsystem('con_cmp2', om.ExecComp('con2 = y2 - 24.0'), promotes=['con2', 'y2'])

model.nonlinear_solver = om.NonlinearBlockGS()

model.linear_solver = om.PETScKrylov()
model.linear_solver.options['atol'] = 1.0e-20

prob.setup()

prob.set_val('x', 1.)
prob.set_val('z', np.array([5.0, 2.0]))

prob.run_model()

wrt = ['z']
of = ['obj']

J = prob.compute_totals(of=of, wrt=wrt, return_format='flat_dict')
NL: NLBGS Converged in 8 iterations
print(J['obj', 'z'][0][0])
print(J['obj', 'z'][0][1])
9.610010556989947
1.7844853356313641

rtol

The relative convergence tolerance, the relative decrease in the (possibly preconditioned) residual norm.

from openmdao.test_suite.components.sellar import SellarDis1withDerivatives, SellarDis2withDerivatives

prob = om.Problem()
model = prob.model

model.add_subsystem('d1', SellarDis1withDerivatives(), promotes=['x', 'z', 'y1', 'y2'])
model.add_subsystem('d2', SellarDis2withDerivatives(), promotes=['z', 'y1', 'y2'])


model.add_subsystem('obj_cmp', om.ExecComp('obj = x**2 + z[1] + y1 + exp(-y2)',
                                           z=np.array([0.0, 0.0]), x=0.0),
                    promotes=['obj', 'x', 'z', 'y1', 'y2'])

model.add_subsystem('con_cmp1', om.ExecComp('con1 = 3.16 - y1'), promotes=['con1', 'y1'])
model.add_subsystem('con_cmp2', om.ExecComp('con2 = y2 - 24.0'), promotes=['con2', 'y2'])

model.nonlinear_solver = om.NonlinearBlockGS()

model.linear_solver = om.PETScKrylov()
model.linear_solver.options['rtol'] = 1.0e-20

prob.setup()

prob.set_val('x', 1.)
prob.set_val('z', np.array([5.0, 2.0]))

prob.run_model()

wrt = ['z']
of = ['obj']

J = prob.compute_totals(of=of, wrt=wrt, return_format='flat_dict')
NL: NLBGS Converged in 8 iterations
print(J['obj', 'z'][0][0])
print(J['obj', 'z'][0][1])
9.610010556989947
1.7844853356313641

ksp_type

You can specify which PETSc algorithm to use in place of ‘fgmres’ by settng the “ksp_type” in the options dictionary. Here, we use ‘gmres’ instead.

from openmdao.test_suite.components.sellar import SellarDis1withDerivatives, SellarDis2withDerivatives

prob = om.Problem()
model = prob.model

model.add_subsystem('d1', SellarDis1withDerivatives(), promotes=['x', 'z', 'y1', 'y2'])
model.add_subsystem('d2', SellarDis2withDerivatives(), promotes=['z', 'y1', 'y2'])

model.add_subsystem('obj_cmp', om.ExecComp('obj = x**2 + z[1] + y1 + exp(-y2)',
                                           z=np.array([0.0, 0.0]), x=0.0),
                    promotes=['obj', 'x', 'z', 'y1', 'y2'])

model.add_subsystem('con_cmp1', om.ExecComp('con1 = 3.16 - y1'), promotes=['con1', 'y1'])
model.add_subsystem('con_cmp2', om.ExecComp('con2 = y2 - 24.0'), promotes=['con2', 'y2'])

model.nonlinear_solver = om.NonlinearBlockGS()

model.linear_solver = om.PETScKrylov()
model.linear_solver.options['ksp_type'] = 'gmres'

prob.setup()

prob.set_val('x', 1.)
prob.set_val('z', np.array([5.0, 2.0]))

prob.run_model()

wrt = ['z']
of = ['obj']

J = prob.compute_totals(of=of, wrt=wrt, return_format='flat_dict')
NL: NLBGS Converged in 8 iterations
print(J['obj', 'z'][0][0])
print(J['obj', 'z'][0][1])
9.610010556989947
1.7844853356313641

Specifying a Preconditioner

You can specify a preconditioner to improve the convergence of the iterative linear solution by setting the precon attribute. The motivation for using a preconditioner is the observation that iterative methods have better convergence properties if the linear system has a smaller condition number, so the goal of the preconditioner is to improve the condition number in part or all of the Jacobian.

Here, we add a Gauss-Seidel preconditioner to the simple Sellar solution with Newton. Note that the number of GMRES iterations is lower when using the preconditioner.

from openmdao.test_suite.components.sellar import SellarDis1withDerivatives, SellarDis2withDerivatives

prob = om.Problem()
model = prob.model

model.add_subsystem('d1', SellarDis1withDerivatives(), promotes=['x', 'z', 'y1', 'y2'])
model.add_subsystem('d2', SellarDis2withDerivatives(), promotes=['z', 'y1', 'y2'])

model.add_subsystem('obj_cmp', om.ExecComp('obj = x**2 + z[1] + y1 + exp(-y2)',
                                           z=np.array([0.0, 0.0]), x=0.0),
                    promotes=['obj', 'x', 'z', 'y1', 'y2'])

model.add_subsystem('con_cmp1', om.ExecComp('con1 = 3.16 - y1'), promotes=['con1', 'y1'])
model.add_subsystem('con_cmp2', om.ExecComp('con2 = y2 - 24.0'), promotes=['con2', 'y2'])

model.nonlinear_solver = om.NewtonSolver(solve_subsystems=False)
model.linear_solver = om.PETScKrylov()

model.linear_solver.precon = om.LinearBlockGS()
model.linear_solver.precon.options['maxiter'] = 2

prob.setup()

prob.set_val('x', 1.)
prob.set_val('z', np.array([5.0, 2.0]))

prob.run_model()
|  | precon:LN: LNBGSSolver 'LN: LNBGS' on system '' failed to converge in 2 iterations.
|  | precon:LN: LNBGSSolver 'LN: LNBGS' on system '' failed to converge in 2 iterations.
|  | precon:LN: LNBGSSolver 'LN: LNBGS' on system '' failed to converge in 2 iterations.
|  | precon:LN: LNBGSSolver 'LN: LNBGS' on system '' failed to converge in 2 iterations.
|  | precon:LN: LNBGSSolver 'LN: LNBGS' on system '' failed to converge in 2 iterations.
|  | precon:LN: LNBGSSolver 'LN: LNBGS' on system '' failed to converge in 2 iterations.
NL: Newton Converged in 3 iterations
print(prob.get_val('y1'))
print(prob.get_val('y2'))
[25.58830237]
[12.05848815]

While the default preconditioning “side” is right-preconditioning, you can also use left-preconditioning provided that you choose a “ksp_type” that supports it. Here we solve the same problem with left-preconditioning using the Richardson method and a DirectSolver.

from openmdao.test_suite.components.sellar import SellarDis1withDerivatives, SellarDis2withDerivatives

prob = om.Problem()
model = prob.model

model.add_subsystem('d1', SellarDis1withDerivatives(), promotes=['x', 'z', 'y1', 'y2'])
model.add_subsystem('d2', SellarDis2withDerivatives(), promotes=['z', 'y1', 'y2'])

model.add_subsystem('obj_cmp', om.ExecComp('obj = x**2 + z[1] + y1 + exp(-y2)',
                                           z=np.array([0.0, 0.0]), x=0.0),
                    promotes=['obj', 'x', 'z', 'y1', 'y2'])

model.add_subsystem('con_cmp1', om.ExecComp('con1 = 3.16 - y1'), promotes=['con1', 'y1'])
model.add_subsystem('con_cmp2', om.ExecComp('con2 = y2 - 24.0'), promotes=['con2', 'y2'])

model.nonlinear_solver = om.NewtonSolver(solve_subsystems=False)
model.linear_solver = om.PETScKrylov()

model.linear_solver.precon = om.DirectSolver()
model.linear_solver.options['precon_side'] = 'left'
model.linear_solver.options['ksp_type'] = 'richardson'

prob.setup()

prob.set_val('x', 1.)
prob.set_val('z', np.array([5.0, 2.0]))

prob.run_model()
NL: Newton Converged in 3 iterations
print(prob.get_val('y1'))
print(prob.get_val('y2'))
[25.58830237]
[12.05848815]