ImplicitComponent

Implicit variables are those that are computed as an implicit function of other variables. For instance, \(y\) would be an implicit variable, given that it is computed by solving \(\cos(x \cdot y) - z \cdot y = 0\). In OpenMDAO, implicit variables are defined as the outputs of components that inherit from ImplicitComponent.

In the above implict expression, \(y\) is the implicit variable while \(x\) and \(z\) would be considered inputs.

ImplicitComponent Methods

The implementation of each method will be illustrated using a simple implicit component that computes the output \(x\) implicitly via a quadratic equation, \(ax^2 + bx + c =0\), where \(a\), \(b\), and \(c\) are inputs to the component.

class QuadraticComp(ImplicitComponent):
    """
    A Simple Implicit Component representing a Quadratic Equation.

    R(a, b, c, x) = ax^2 + bx + c

    Solution via Quadratic Formula:
    x = (-b + sqrt(b^2 - 4ac)) / 2a
    """
  • setup() :

    Declare input and output variables via add_input and add_output. Information like variable names, sizes, units, and bounds are declared. Also, declare partial derivatives that this component provides. Here we use the wild card to say that this component provides derivatives of all implicit residuals with respect to all inputs and outputs.

    def setup(self):
        self.add_input('a', val=1.)
        self.add_input('b', val=1.)
        self.add_input('c', val=1.)
        self.add_output('x', val=0.)
    
        self.declare_partials(of='*', wrt='*')
    
  • apply_nonlinear(inputs, outputs, residuals) :

    Compute the residuals, given the inputs and outputs.

    def apply_nonlinear(self, inputs, outputs, residuals):
        a = inputs['a']
        b = inputs['b']
        c = inputs['c']
        x = outputs['x']
        residuals['x'] = a * x ** 2 + b * x + c
    
  • solve_nonlinear(inputs, outputs) :

    [Optional] Compute the outputs, given the inputs.

    def solve_nonlinear(self, inputs, outputs):
        a = inputs['a']
        b = inputs['b']
        c = inputs['c']
        outputs['x'] = (-b + (b ** 2 - 4 * a * c) ** 0.5) / (2 * a)
    
  • linearize(inputs, outputs, partials) :

    [Optional] The component’s partial derivatives of interest are those of the residuals with respect to the inputs and the outputs. If the user computes the partial derivatives explicitly, they are provided here. If the user wants to implement partial derivatives in a matrix-free way, this method provides a place to perform any necessary assembly or pre-processing for the matrix-vector products. Regardless of how the partial derivatives are computed, this method provides a place to perform any relevant factorizations for directly solving or preconditioning the linear system.

    def linearize(self, inputs, outputs, partials):
        a = inputs['a']
        b = inputs['b']
        c = inputs['c']
        x = outputs['x']
    
        partials['x', 'a'] = x ** 2
        partials['x', 'b'] = x
        partials['x', 'c'] = 1.0
        partials['x', 'x'] = 2 * a * x + b
    
        self.inv_jac = 1.0 / (2 * a * x + b)
    
  • apply_linear(inputs, outputs, d_inputs, d_outputs, d_residuals, mode) :

    [Optional] If the user wants to implement partial derivatives in a matrix-free way, this method performs the matrix-vector product. If mode is ‘fwd’, this method computes \(d\_{residuals} = J \cdot [ d\_{inputs} \quad d\_{outputs} ]^T\). If mode is ‘rev’, this method computes \([ d\_{inputs} \quad d\_{outputs} ]^T = J^T \cdot d\_{residuals}\).

    def apply_linear(self, inputs, outputs,
                     d_inputs, d_outputs, d_residuals, mode):
        a = inputs['a']
        b = inputs['b']
        c = inputs['c']
        x = outputs['x']
        if mode == 'fwd':
            if 'x' in d_residuals:
                if 'x' in d_outputs:
                    d_residuals['x'] += (2 * a * x + b) * d_outputs['x']
                if 'a' in d_inputs:
                    d_residuals['x'] += x ** 2 * d_inputs['a']
                if 'b' in d_inputs:
                    d_residuals['x'] += x * d_inputs['b']
                if 'c' in d_inputs:
                    d_residuals['x'] += d_inputs['c']
        elif mode == 'rev':
            if 'x' in d_residuals:
                if 'x' in d_outputs:
                    d_outputs['x'] += (2 * a * x + b) * d_residuals['x']
                if 'a' in d_inputs:
                    d_inputs['a'] += x ** 2 * d_residuals['x']
                if 'b' in d_inputs:
                    d_inputs['b'] += x * d_residuals['x']
                if 'c' in d_inputs:
                    d_inputs['c'] += d_residuals['x']
    
  • apply_multi_linear(inputs, outputs, d_inputs, d_outputs, d_residuals, mode) :

    [Optional] If the user wants to implement partial derivatives in a matrix-free way, this method performs the matrix-matrix product. If mode is ‘fwd’, this method computes \(d\_{residuals} = J \cdot [ d\_{inputs} \quad d\_{outputs} ]^T\) where d_outputs and d_residuals are both matrices. If mode is ‘rev’, this method similarly computes \([ d\_{inputs} \quad d\_{outputs} ]^T = J^T \cdot d\_{residuals}\).

    This method is only used when “vectorize_derivs” is set to True on a design variable or response.

    def apply_multi_linear(self, inputs, outputs, d_inputs, d_outputs, d_residuals, mode):
        a = inputs['a']
        b = inputs['b']
        c = inputs['c']
        x = outputs['x']
        if mode == 'fwd':
            if 'x' in d_residuals:
                if 'x' in d_outputs:
                    d_residuals['x'] += (2 * a * x + b) * d_outputs['x']
                if 'a' in d_inputs:
                    d_residuals['x'] += x ** 2 * d_inputs['a']
                if 'b' in d_inputs:
                    d_residuals['x'] += x * d_inputs['b']
                if 'c' in d_inputs:
                    d_residuals['x'] += d_inputs['c']
        elif mode == 'rev':
            if 'x' in d_residuals:
                if 'x' in d_outputs:
                    d_outputs['x'] += (2 * a * x + b) * d_residuals['x']
                if 'a' in d_inputs:
                    d_inputs['a'] += x ** 2 * d_residuals['x']
                if 'b' in d_inputs:
                    d_inputs['b'] += x * d_residuals['x']
                if 'c' in d_inputs:
                    d_inputs['c'] += d_residuals['x']
    
  • solve_linear(d_outputs, d_residuals, mode) :

    [Optional] Solves a linear system where the matrix is \(d\_{residuals} / d\_{outputs}\) or its transpose. If mode is ‘fwd’, the right-hand side vector is \(d\_{residuals}\) and the solution vector is \(d\_{outputs}\). If mode is ‘rev’, the right-hand side vector is \(d\_{outputs}\) and the solution vector is \(d\_{residuals}\).

    def solve_linear(self, d_outputs, d_residuals, mode):
        if mode == 'fwd':
            d_outputs['x'] = self.inv_jac * d_residuals['x']
        elif mode == 'rev':
            d_residuals['x'] = self.inv_jac * d_outputs['x']
    
  • guess_nonlinear(self, inputs, outputs, residuals) :

    [Optional] This method allows the user to calculate and specify an initial guess for implicit states. It is called prior to every call to solve_nonlinear, so it is useful for when you would like to “reset” the initial conditions on an inner-nested solve whenever an outer loop solver or driver changes other values. Since it is a hook for custom code, you could also use it to monitor variables in the input, output, or residual vectors and change the initial guess when some condition is met. Here is a simple example where we use NewtonSolver to find one of the roots of a second-order quadratic equation. Which root you get depends on the initial guess.

    from openmdao.api import Problem, Group, ImplicitComponent, IndepVarComp, NewtonSolver, ScipyKrylov
    
    class ImpWithInitial(ImplicitComponent):
    
        def setup(self):
            self.add_input('a', val=1.)
            self.add_input('b', val=1.)
            self.add_input('c', val=1.)
            self.add_output('x', val=0.)
    
            self.declare_partials(of='*', wrt='*')
    
        def apply_nonlinear(self, inputs, outputs, residuals):
            a = inputs['a']
            b = inputs['b']
            c = inputs['c']
            x = outputs['x']
            residuals['x'] = a * x ** 2 + b * x + c
    
        def solve_nonlinear(self, inputs, outputs):
            a = inputs['a']
            b = inputs['b']
            c = inputs['c']
            outputs['x'] = (-b + (b ** 2 - 4 * a * c) ** 0.5) / 2 / a
    
        def linearize(self, inputs, outputs, partials):
            a = inputs['a']
            b = inputs['b']
            c = inputs['c']
            x = outputs['x']
    
            partials['x', 'a'] = x ** 2
            partials['x', 'b'] = x
            partials['x', 'c'] = 1.0
            partials['x', 'x'] = 2 * a * x + b
    
            self.inv_jac = 1.0 / (2 * a * x + b)
    
        def solve_nonlinear(self, inputs, outputs):
            """ Do nothing. """
            pass
    
        def guess_nonlinear(self, inputs, outputs, resids):
            # Solution at 1 and 3. Default value takes us to -1 solution. Here
            # we set it to a value that will tke us to the 3 solution.
            outputs['x'] = 5.0
    
    prob = Problem()
    model = prob.model = Group()
    
    model.add_subsystem('pa', IndepVarComp('a', 1.0))
    model.add_subsystem('pb', IndepVarComp('b', 1.0))
    model.add_subsystem('pc', IndepVarComp('c', 1.0))
    model.add_subsystem('comp2', ImpWithInitial())
    model.connect('pa.a', 'comp2.a')
    model.connect('pb.b', 'comp2.b')
    model.connect('pc.c', 'comp2.c')
    
    model.nonlinear_solver = NewtonSolver()
    model.nonlinear_solver.options['solve_subsystems'] = True
    model.nonlinear_solver.options['max_sub_solves'] = 1
    model.linear_solver = ScipyKrylov()
    
    prob.setup(check=False)
    
    prob['pa.a'] = 1.
    prob['pb.b'] = -4.
    prob['pc.c'] = 3.
    
    prob.run_model()
    
    print(prob['comp2.x'])