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.

SellarDis1withDerivatives class definition

class SellarDis1withDerivatives(SellarDis1):
    """
    Component containing Discipline 1 -- derivatives version.
    """

    def setup_partials(self):
        # Analytic Derivs
        self.declare_partials(of='*', wrt='*')

    def compute_partials(self, inputs, partials):
        """
        Jacobian for Sellar discipline 1.
        """
        partials['y1', 'y2'] = -0.2
        partials['y1', 'z'] = np.array([[2.0 * inputs['z'][0], 1.0]])
        partials['y1', 'x'] = 1.0

SellarDis2withDerivatives class definition

class SellarDis2withDerivatives(SellarDis2):
    """
    Component containing Discipline 2 -- derivatives version.
    """

    def setup_partials(self):
        # Analytic Derivs
        self.declare_partials(of='*', wrt='*')

    def compute_partials(self, inputs, J):
        """
        Jacobian for Sellar discipline 2.
        """
        y1 = inputs['y1']
        if y1.real < 0.0:
            y1 *= -1
        if y1.real < 1e-8:
            y1 = 1e-8

        J['y2', 'y1'] = .5*y1**-.5
        J['y2', 'z'] = np.array([[1.0, 1.0]])

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.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

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

9.610010557013222
1.7844853356442276

LinearBlockGS Options

Option	Default	Acceptable Values	Acceptable Types	Description
aitken_initial_factor	1.0	N/A	N/A	initial value for Aitken relaxation factor
aitken_max_factor	1.5	N/A	N/A	upper limit for Aitken relaxation factor
aitken_min_factor	0.1	N/A	N/A	lower limit for Aitken relaxation factor
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_converge	False	[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-10	N/A	N/A	relative error tolerance
use_aitken	False	[True, False]	['bool']	set to True to use Aitken relaxation

LinearBlockGS Constructor

The call signature for the LinearBlockGS constructor is:

LinearBlockGS.__init__(**kwargs)

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('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.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.

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

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('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.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

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

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('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.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

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

9.610162961754462
1.7845695570436657