# LinearBlockGS¶

LinearBlockGS uses Block Gauss-Seidel to solve the linear system. LinearBlockGS iterates until the linear residual is below a tolerance, or the maximum number of iterations has been exceeded. As such, it is generally usable for any system topology, and can handle cycles and implicit states alike. It is not always the best solver to choose, however, and is known to diverge or plateau on some problems. In such a case, you may need to use a solver such as ScipyKrylov.

LinearBlockGS is a block solver, so you can specify different linear solvers in the subsystems and they will be utilized to solve the subsystem linear problem.

Note that systems without cycles or implicit states will converge in one iteration of Block Gauss-Seidel.

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('obj_cmp', om.ExecComp('obj = x**2 + z + 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.NonlinearBlockGS()

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
LN: LNBGS Converged in 7 iterations
LN: LNBGS Converged in 7 iterations

print(J['obj', 'z'])
print(J['obj', 'z'])

9.610010557013222
1.7844853356442276


## LinearBlockGS Options¶

OptionDefaultAcceptable ValuesAcceptable TypesDescription
aitken_initial_factor1.0N/AN/Ainitial value for Aitken relaxation factor
aitken_max_factor1.5N/AN/Aupper limit for Aitken relaxation factor
aitken_min_factor0.1N/AN/Alower limit for Aitken relaxation factor
assemble_jacFalse[True, False]['bool']Activates use of assembled jacobian by this solver.
atol1e-10N/AN/Aabsolute error tolerance
err_on_non_convergeFalse[True, False]['bool']When True, AnalysisError will be raised if we don't converge.
iprint1N/A['int']whether to print output
maxiter10N/A['int']maximum number of iterations
rtol1e-10N/AN/Arelative error tolerance
use_aitkenFalse[True, False]['bool']set to True to use Aitken relaxation

## LinearBlockGS Constructor¶

The call signature for the LinearBlockGS constructor is:

LinearBlockGS.__init__(**kwargs)[source]

Initialize all attributes.

## Aitken relaxation¶

This solver implements Aitken relaxation, as described in Algorithm 1 of this paper on aerostructual design optimization. The relaxation is turned off by default, but it may help convergence for more tightly coupled models.

## LinearBlockGS Option Examples¶

maxiter

This lets you specify the maximum number of Gauss-Seidel 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 either of the specified absolute or relative tolerances.

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('obj_cmp', om.ExecComp('obj = x**2 + z + 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.LinearBlockGS()
model.linear_solver.options['maxiter'] = 2

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
LN: LNBGSSolver 'LN: LNBGS' on system '' failed to converge in 2 iterations.
LN: LNBGSSolver 'LN: LNBGS' on system '' failed to converge in 2 iterations.

print(J['obj', 'z'])
print(J['obj', 'z'])

9.602301180035953
1.780225005469842


atol

Here, we set the absolute tolerance to a looser value that will trigger an earlier termination. After each iteration, the norm of the linear residuals is calculated by calling apply_linear. 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('obj_cmp', om.ExecComp('obj = x**2 + z + 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.LinearBlockGS()
model.linear_solver.options['atol'] = 1.0e-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
LN: LNBGS Converged in 3 iterations
LN: LNBGS Converged in 3 iterations

print(J['obj', 'z'])
print(J['obj', 'z'])

9.610162961754462
1.7845695570436657


rtol

Here, we set the relative tolerance to a looser value that will trigger an earlier termination. After each iteration, the norm of the linear residuals is calculated by calling apply_linear. 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('obj_cmp', om.ExecComp('obj = x**2 + z + 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.LinearBlockGS()
model.linear_solver.options['rtol'] = 1.0e-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
LN: LNBGS Converged in 3 iterations
LN: LNBGS Converged in 3 iterations

print(J['obj', 'z'])
print(J['obj', 'z'])

9.610162961754462
1.7845695570436657