guess_nonlinear_cycle.py

import openmdao.api as om 
import numpy as np 

class SellarDis1(om.ExplicitComponent):
    """
    Component containing Discipline 1 -- no derivatives version.
    """

    def setup(self):

        # Global Design Variable
        self.add_input('z', val=np.zeros(2))

        # Local Design Variable
        self.add_input('x', val=0.)

        # Coupling parameter
        self.add_input('y2', val=1.0)

        # Coupling output
        self.add_output('y1', val=1.0)

        # Finite difference all partials.
        self.declare_partials('*', '*', method='fd')

    def compute(self, inputs, outputs):
        """
        Evaluates the equation
        y1 = z1**2 + z2 + x1 - 0.2*y2
        """
        z1 = inputs['z'][0]
        z2 = inputs['z'][1]
        x1 = inputs['x']
        y2 = inputs['y2']
        outputs['y1'] = z1**2 + z2 + x1 - 0.2*y2
        print('Dis1: '+str(outputs['y1']))

class SellarDis2(om.ExplicitComponent):
    """
    Component containing Discipline 2 -- no derivatives version.
    """

    def setup(self):
        # Global Design Variable
        self.add_input('z', val=np.zeros(2))

        # Coupling parameter
        self.add_input('y1', val=1.0)

        # Coupling output
        self.add_output('y2', val=1.0)

        # Finite difference all partials.
        self.declare_partials('*', '*', method='fd')

    def compute(self, inputs, outputs):
        """
        Evaluates the equation
        y2 = y1**(.5) + z1 + z2
        """

        z1 = inputs['z'][0]
        z2 = inputs['z'][1]
        y1 = inputs['y1']

        # Note: this may cause some issues. However, y1 is constrained to be
        # above 3.16, so lets just let it converge, and the optimizer will
        # throw it out
        if y1.real < 0.0:
            y1 *= -1

        outputs['y2'] = y1**.5 + z1 + z2
        print('Dis2: '+str(outputs['y2']))

class SellarCycle(om.Group):
    def setup(self):
        self.add_subsystem('d1', SellarDis1(), promotes_inputs=['x', 'z', 'y2'],
                            promotes_outputs=['y1'])
        self.add_subsystem('d2', SellarDis2(), promotes_inputs=['z', 'y1'],
                            promotes_outputs=['y2'])

        self.set_input_defaults('x', 1.0)
        self.set_input_defaults('z', np.array([5.0, 2.0]))

    def guess_nonlinear(self, inputs, outputs, residuals):
        print('in guess nonlinear - these nonconverged residuals shouldnt be zero')
        print(residuals['y1'])
        print(residuals['y2'])

class SellarMDA(om.Group):
    """
    Group containing the Sellar MDA.
    """

    def setup(self):
        cycle = self.add_subsystem('cycle', SellarCycle(), promotes=['*'])
        # Nonlinear Block Gauss Seidel is a gradient free solver
        cycle.nonlinear_solver = om.NonlinearBlockGS()

        self.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=['x', 'z', 'y1', 'y2', 'obj'])

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



if __name__ == "__main__":
    prob = om.Problem(SellarMDA())
    prob.setup(check=True)
    prob.run_model()

    # INFO: checking out_of_order
    # INFO: checking system
    # INFO: checking solvers
    # WARNING: Group 'cycle' contains cycles [['d1', 'd2']], but does not have an iterative linear solver.
    # INFO: checking dup_inputs
    # INFO: checking missing_recorders
    # WARNING: The Problem has no recorder of any kind attached
    # INFO: checking comp_has_no_outputs

    # =====
    # cycle
    # =====
    # in guess nonlinear - these nonconverged residuals shouldnt be zero
    # [0.]
    # [0.]
    # Dis1: [27.8]
    # Dis2: [12.27257053]
    # Dis1: [25.54548589]
    # Dis2: [12.05425424]
    # Dis1: [25.58914915]
    # Dis2: [12.05857185]
    # Dis1: [25.58828563]
    # Dis2: [12.0584865]
    # Dis1: [25.5883027]
    # Dis2: [12.05848818]
    # Dis1: [25.58830236]
    # Dis2: [12.05848815]
    # Dis1: [25.58830237]
    # Dis2: [12.05848815]
    # Dis1: [25.58830237]
    # Dis2: [12.05848815]
    # NL: NLBGS Converged in 7 iterations