ScipyKrylov

ScipyKrylov is an iterative linear solver that wraps the methods found in scipy.sparse.linalg. The default method is “gmres”, or the Generalized Minimal RESidual method. Support for other scipy.sparse.linalg solvers will be added over time. This linear solver is capable of handling any system topology very 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 is a serial solver, so it should never be used under MPI; use PETScKrylov instead.

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

import numpy as np

from openmdao.api import Problem, Group, IndepVarComp, ScipyKrylov, \
     NonlinearBlockGS, ExecComp
from openmdao.test_suite.components.sellar import SellarDis1withDerivatives, \
     SellarDis2withDerivatives

prob = Problem()
model = prob.model

model.add_subsystem('px', IndepVarComp('x', 1.0), promotes=['x'])
model.add_subsystem('pz', IndepVarComp('z', np.array([5.0, 2.0])), promotes=['z'])

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

model.add_subsystem('obj_cmp', 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', ExecComp('con1 = 3.16 - y1'), promotes=['con1', 'y1'])
model.add_subsystem('con_cmp2', ExecComp('con2 = y2 - 24.0'), promotes=['con2', 'y2'])

model.nonlinear_solver = NonlinearBlockGS()

model.linear_solver = ScipyKrylov()

prob.setup()
prob.run_model()
NL: NLBGS Converged in 7 iterations
wrt = ['z']
of = ['obj']

J = prob.compute_totals(of=of, wrt=wrt, return_format='flat_dict')
print(J['obj', 'z'][0][0])
9.61001055699
print(J['obj', 'z'][0][1])
1.78448533563

ScipyKrylov Options

Option Default Acceptable Values Acceptable Types Description
assemble_jac False N/A [‘bool’] Activates use of assembled jacobian by this solver.
atol 1e-12 N/A N/A absolute error tolerance
err_on_maxiter False N/A [‘bool’] When True, AnalysisError will be raised if we don’t converge.
iprint 1 N/A [‘int’] whether to print output
maxiter 1000 N/A [‘int’] maximum number of iterations
restart 20 N/A [‘int’] Number of iterations between restarts. Larger values increase iteration cost, but may be necessary for convergence. This option applies only to gmres.
rtol 1e-10 N/A N/A relative error tolerance
solver gmres [‘gmres’] N/A function handle for actual solver

Note: Options can be passed as keyword arguments at initialization.

ScipyKrylov Option Examples

maxiter

maxiter lets you specify the maximum number of GMRES iterations to apply. The default maximum is 1000, 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 it 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.)

import numpy as np

from openmdao.api import Problem, Group, IndepVarComp, ScipyKrylov, NonlinearBlockGS, ExecComp
from openmdao.test_suite.components.sellar import SellarDis1withDerivatives, SellarDis2withDerivatives

prob = Problem()
model = prob.model

model.add_subsystem('px', IndepVarComp('x', 1.0), promotes=['x'])
model.add_subsystem('pz', IndepVarComp('z', np.array([5.0, 2.0])), promotes=['z'])

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

model.add_subsystem('obj_cmp', 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', ExecComp('con1 = 3.16 - y1'), promotes=['con1', 'y1'])
model.add_subsystem('con_cmp2', ExecComp('con2 = y2 - 24.0'), promotes=['con2', 'y2'])

model.nonlinear_solver = NonlinearBlockGS()

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

prob.setup()
prob.run_model()
NL: NLBGS Converged in 7 iterations
wrt = ['z']
of = ['obj']

J = prob.compute_totals(of=of, wrt=wrt, return_format='flat_dict')
print(J['obj', 'z'][0][0])
0.0
print(J['obj', 'z'][0][1])
0.0

atol

Here, we set the absolute tolerance to a much tighter value (default is 1.0e-12) to show what happens. In practice, the tolerance serves a dual role in GMRES. In addition to being a termination criteria, the tolerance also defines what GMRES considers to be tiny. Tiny numbers are replaced by zero when the argument vector is normalized at the start of each new matrix-vector product. The end result here is that we iterate longer to get a marginally better answer.

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

import numpy as np

from openmdao.api import Problem, Group, IndepVarComp, ScipyKrylov, NonlinearBlockGS, ExecComp
from openmdao.test_suite.components.sellar import SellarDis1withDerivatives, SellarDis2withDerivatives

prob = Problem()
model = prob.model

model.add_subsystem('px', IndepVarComp('x', 1.0), promotes=['x'])
model.add_subsystem('pz', IndepVarComp('z', np.array([5.0, 2.0])), promotes=['z'])

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

model.add_subsystem('obj_cmp', 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', ExecComp('con1 = 3.16 - y1'), promotes=['con1', 'y1'])
model.add_subsystem('con_cmp2', ExecComp('con2 = y2 - 24.0'), promotes=['con2', 'y2'])

model.nonlinear_solver = NonlinearBlockGS()

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

prob.setup()
prob.run_model()
NL: NLBGS Converged in 7 iterations
wrt = ['z']
of = ['obj']

J = prob.compute_totals(of=of, wrt=wrt, return_format='flat_dict')
print(J['obj', 'z'][0][0])
9.61001055699
print(J['obj', 'z'][0][1])
1.78448533563

rtol

The ‘rtol’ setting is not supported by Scipy GMRES.

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.

import numpy as np

from openmdao.api import Problem, Group, IndepVarComp, ScipyKrylov, NewtonSolver, \
     LinearBlockGS, ExecComp
from openmdao.test_suite.components.sellar import SellarDis1withDerivatives, \
     SellarDis2withDerivatives

prob = Problem()
model = prob.model

model.add_subsystem('px', IndepVarComp('x', 1.0), promotes=['x'])
model.add_subsystem('pz', IndepVarComp('z', np.array([5.0, 2.0])), promotes=['z'])

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

model.add_subsystem('obj_cmp', 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', ExecComp('con1 = 3.16 - y1'), promotes=['con1', 'y1'])
model.add_subsystem('con_cmp2', ExecComp('con2 = y2 - 24.0'), promotes=['con2', 'y2'])

model.nonlinear_solver = NewtonSolver()
model.linear_solver = ScipyKrylov()

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

prob.setup()
prob.run_model()
|  | precon:LN: LNBGS Failed to Converge in 2 iterations
|  | precon:LN: LNBGS Converged in 0 iterations
|  | precon:LN: LNBGS Failed to Converge in 2 iterations
|  | precon:LN: LNBGS Converged in 0 iterations
|  | precon:LN: LNBGS Failed to Converge in 2 iterations
|  | precon:LN: LNBGS Converged in 0 iterations
|  | precon:LN: LNBGS Failed to Converge in 2 iterations
|  | precon:LN: LNBGS Converged in 0 iterations
|  | precon:LN: LNBGS Failed to Converge in 2 iterations
|  | precon:LN: LNBGS Converged in 0 iterations
|  | precon:LN: LNBGS Failed to Converge in 2 iterations
|  | precon:LN: LNBGS Converged in 0 iterations
|  | precon:LN: LNBGS Failed to Converge in 2 iterations
|  | precon:LN: LNBGS Converged in 0 iterations
|  | precon:LN: LNBGS Failed to Converge in 2 iterations
|  | precon:LN: LNBGS Converged in 0 iterations
NL: Newton Converged in 8 iterations
print(prob['y1'])
[ 25.58830237]
print(prob['y2'])
[ 12.05848815]

A note on nesting ScipyKrylov under a preconditoner: The underlying GMRES module is not re-entrant, so it cannot be called as a new instance while it is running. If you need to use gmres under gmres in a preconditioner stack, you should use PETScKrylov at one (ore more) of the levels.