How to solve an equation system as an array structure in OpenMDAO?

Question

I tried to solve the implicit equation shown below using OpenMDAO. The equations are shown below,

x[i]*z[i] + z[i] - 4 = 0,  
y[i] = x[i] + 2*z[i]

The solution is (for i=2) z = [2.666667, 2], y = [5.833333, 5] for x = [0.5, 1.0].

For this case, I have used the code which is shown below,

from __future__ import print_function
from openmdao.api import Component, Group, Problem, Newton, ScipyGMRES, IndepVarComp

class SimpleEquationSystem(Component):
    """Solve the Equation 
            x[i]*z[i] + z[i] - 4 = 0
            y[i] = x[i] + 2*z[i]

       Solution: z = [2.666667, 2], y = [5.833333, 5] for x = [0.5, 1.0]
    """

    def __init__(self):
        super(SimpleEquationSystem, self).__init__()


        self.add_param('x', [0.5, 1.0])
        self.add_state('y', [0.0,0.0])
        self.add_state('z', [0.0,0.0])
        self.iter=0

    def solve_nonlinear(self, params, unknowns, resids):
        """This component does no calculation on its own. It mainly holds the
        initial value of the state. An OpenMDAO solver outside of this
        component varies it to drive the residual to zero."""
        pass

    def apply_nonlinear(self, params, unknowns, resids):
        """ Report the residual """
        self.iter+=1
        for i in range(2):
            x=params['x'][i]
            y = unknowns['y'][i]
            z = unknowns['z'][i]


            resids['y'][i] = x*z + z - 4
            resids['z'][i] = x + 2*z - y

        print('y_%d' % self.iter,'=%s' %resids['y'], 'z_%d' % self.iter, '=%s' %resids['z'])
        print('x' ,'=%s' %x, 'y', '=%s' %y, 'z', '=%s' %z)

top = Problem()
root = top.root = Group()
root.add('comp', SimpleEquationSystem())
root.add('p1', IndepVarComp('x', [0.5, 1.0]))
root.connect('p1.x', 'comp.x')
# Tell these components to finite difference
root.comp.deriv_options['type'] = 'fd'
root.comp.deriv_options['form'] = 'central'
root.comp.deriv_options['step_size'] = 1.0e-4
root.nl_solver = Newton()
root.nl_solver.options['maxiter']=int(200)
root.ln_solver = ScipyGMRES()

top.setup()
top.print_all_convergence(level=1, depth=2)
top.run()
print('Solution x=%s, y=%s, z=%s' % (top['comp.x'], top['comp.y'], top['comp.z']))

This code errors out with the message "RuntimeWarning: invalid value encountered in double_scalars". One should be able to get this error in OpenMDAO while using the above code fragment.

When the residual equations are implemented as scalars instead of vectors as shown below, it works fine.

resids['z1']= params['x1'] + 2*unknowns['z1'] - unknowns['y1']
resids['z2']= params['x2'] + 2*unknowns['z2'] - unknowns['y2']
resids['y1']= params['x1']*unknowns['z1'] + unknowns['z1'] - 4
resids['y2']= params['x2']*unknowns['z2'] + unknowns['z2'] - 4

But I would like to have residuals as vectors so that it is easier to handle with a for loop. Can you please help me to resolve this problem?

Answer 1

When you declare a variable, it is precisely the variable type that you declare it. You gave it a list, which OpenMDAO interprets as a list, which isn't a data type that support differentiation, so Newton can't do anything with it. You need to make them numpy arrays by making the following changes:

import numpy as np

and

self.add_param('x', np.array([0.5, 1.0]))
self.add_state('y', np.array([0.0,0.0]))
self.add_state('z', np.array([0.0,0.0]))

also:

root.add('p1', IndepVarComp('x', np.array([0.5, 1.0])))