NonlinearBlockJac

NonlinearBlockJac is a nonlinear solver that uses the block Jacobi method to solve the system. When to choose this solver over NonlinearBlockGS is an advanced topic, but it is valid for systems that satisfy the same conditions:

  1. System (or subsystem) contains a cycle, though subsystems may.

  2. System does not contain any implicit states, though subsystems may.

Note that you may not know if you satisfy the second condition, so choosing a solver can be “trial and error.” If NonlinearBlockJac doesn’t work, then you will need to use NewtonSolver

The main difference over NonlinearBlockGS is that data passing is delayed until after all subsystems have been executed.

Here, we choose NonlinearBlockJac to solve the Sellar problem, which has two components with a cyclic dependency, has no implicit states, and works very well with Jacobi.

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.linear_solver = om.LinearBlockGS()
model.nonlinear_solver = om.NonlinearBlockJac()

prob.setup()

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

prob.run_model()

print(prob['y1'])
print(prob['y2'])
NL: NLBJSolver 'NL: NLBJ' on system '' failed to converge in 10 iterations.
[25.58830249]
[12.05848818]

This solver runs all of the subsystems each iteration, but just passes the data along all connections simultaneously once per iteration. After each iteration, the iteration count and the residual norm are checked to see if termination has been satisfied.

You can control the termination criteria for the solver using the following options:

NonlinearBlockJac Options

Option DefaultAcceptable Values Acceptable Types Description
atol 1e-10N/A N/A absolute error tolerance
debug_print 0 [True, False] ['bool'] If true, the values of input and output variables at the start of iteration are printed and written to a file after a failure to converge.
err_on_non_converge 0 [True, False] ['bool'] When True, AnalysisError will be raised if we don't converge.
iprint 1 N/A ['int'] whether to print output
maxiter 10 N/A ['int'] maximum number of iterations
rtol 1e-10N/A N/A relative error tolerance
stall_limit 0 N/A N/A Number of iterations after which, if the residual norms are identical within the stall_tol, then terminate as if max iterations were reached. Default is 0, which disables this feature.
stall_tol 1e-12N/A N/A When stall checking is enabled, the threshold below which the residual norm is considered unchanged.

NonlinearBlockJac Constructor

The call signature for the NonlinearBlockJac constructor is:

NonlinearBlockJac.__init__(**kwargs)

Initialize all attributes.

NonlinearBlockJac Option Examples

maxiter

This lets you specify the maximum number of Jacobi iterations to apply. In this example, we cut it back from the default, ten, down to two, so that it terminates a few iterations earlier and doesn’t reach the specified absolute or relative tolerance.

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.linear_solver = om.LinearBlockGS()

nlbgs = model.nonlinear_solver = om.NonlinearBlockJac()
nlbgs.options['maxiter'] = 4

prob.setup()

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

prob.run_model()

print(prob['y1'])
print(prob['y2'])
NL: NLBJSolver 'NL: NLBJ' on system '' failed to converge in 4 iterations.
[25.57238139]
[12.05425424]

atol

Here, we set the absolute tolerance to a looser value that will trigger an earlier termination. After each iteration, the norm of the residuals is calculated by calling apply_nonlinear on all of the components. If this norm value is lower than the absolute tolerance atol, the iteration will terminate.

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.linear_solver = om.LinearBlockGS()

nlbgs = model.nonlinear_solver = om.NonlinearBlockJac()
nlbgs.options['atol'] = 1e-2

prob.setup()

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

prob.run_model()

print(prob['y1'])
print(prob['y2'])
NL: NLBJ Converged in 6 iterations
[25.58861716]
[12.05857185]

rtol

Here, we set the relative tolerance to a looser value that will trigger an earlier termination. After each iteration, the norm of the residuals is calculated by calling apply_nonlinear on all of the components. If the ratio of the currently calculated norm to the initial residual norm is lower than the relative tolerance rtol, the iteration will terminate.

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.linear_solver = om.LinearBlockGS()

nlbgs = model.nonlinear_solver = om.NonlinearBlockJac()
nlbgs.options['rtol'] = 1e-3

prob.setup()

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

prob.run_model()

print(prob.get_val('y1'))
print(prob.get_val('y2'))
NL: NLBJ Converged in 5 iterations
[25.58914915]
[12.05691422]