Merge d40f0e4 into bb60dd4

openmdao/solvers/nonlinear/nonlinear_block_gs.py

@@ -26,11 +26,59 @@ def _setup_solvers(self, system, depth):
         if len(system._subsystems_allprocs) != len(system._subsystems_myproc):
             raise RuntimeError('Nonlinear Gauss-Seidel cannot be used on a parallel group.')
 
+    def _declare_options(self):
+        """
+        Declare options before kwargs are processed in the init method.
+        """
+        self.options.declare('use_aitken', type_=bool, default=False,
+                             desc='set to True to use Aitken relaxation')
+        self.options.declare('aitken_min_factor', default=0.1,
+                             desc='lower limit for Aitken relaxation factor')
+        self.options.declare('aitken_max_factor', default=1.5,
+                             desc='upper limit for Aitken relaxation factor')
+
+    def _iter_initialize(self):
+        """
+        Perform any necessary pre-processing operations.
+
+        Returns
+        -------
+        float
+            initial error.
+        float
+            error at the first iteration.
+        """
+        if self.options['use_aitken']:
+            self._aitken_work1 = self._system._outputs._clone()
+            self._aitken_work2 = self._system._outputs._clone()
+            self._aitken_work3 = self._system._outputs._clone()
+            self._aitken_work4 = self._system._outputs._clone()
+            self._theta_n_1 = 1.
+
+        return super(NonlinearBlockGS, self)._iter_initialize()
+
     def _iter_execute(self):
         """
         Perform the operations in the iteration loop.
         """
         system = self._system
+        use_aitken = self.options['use_aitken']
+
+        if use_aitken:
+            outputs = self._system._outputs
+            aitken_min_factor = self.options['aitken_min_factor']
+            aitken_max_factor = self.options['aitken_max_factor']
+
+            # some variables that are used for Aitken's relaxation
+            delta_outputs_n_1 = self._aitken_work1
+            delta_outputs_n = self._aitken_work2
+            outputs_n = self._aitken_work3
+            temp = self._aitken_work4
+            theta_n_1 = self._theta_n_1
+            # store a copy of the outputs, used to compute the change in outputs later
+            delta_outputs_n.set_vec(outputs)
+            # store a copy of the outputs
+            outputs_n.set_vec(outputs)
 
         self._solver_info.append_subsolver()
         for isub, subsys in enumerate(system._subsystems_myproc):
@@ -40,6 +88,37 @@ def _iter_execute(self):
 
         self._solver_info.pop()
 
+        if use_aitken:
+            # compute the change in the outputs after the NLBGS iteration
+            delta_outputs_n -= outputs
+            delta_outputs_n *= -1
+
+            if self._iter_count >= 2:
+                # Compute relaxation factor. This method is used by Kenway et al. in
+                # "Scalable Parallel Approach for High-Fidelity Steady-State Aero-
+                # elastic Analysis and Adjoint Derivative Computations" (ln 22 of Algo 1)
+
+                temp.set_vec(delta_outputs_n)
+                temp -= delta_outputs_n_1
+                temp_norm = temp.get_norm()
+                if temp_norm == 0.:
+                    temp_norm = 1e-12  # prevent division by 0 in the next line
+                theta_n = theta_n_1 * (1 - temp.dot(delta_outputs_n) / temp_norm ** 2)
+                # limit relaxation factor to the specified range
+                theta_n = max(aitken_min_factor, min(aitken_max_factor, theta_n))
+                # save relaxation factor for the next iteration
+                theta_n_1 = theta_n
+            else:
+                theta_n = 1.
+
+            outputs.set_vec(outputs_n)
+
+            # compute relaxed outputs
+            outputs.add_scal_vec(theta_n, delta_outputs_n)
+
+            # save update to use in next iteration
+            delta_outputs_n_1.set_vec(delta_outputs_n)
+
     def _mpi_print_header(self):
         """
         Print header text before solving.

openmdao/solvers/nonlinear/tests/test_nonlinear_block_gs.py

@@ -307,6 +307,19 @@ def test_run_apply(self):
         self.assertTrue(s1p1._run_apply != s2p1._run_apply)
         self.assertTrue(s1p2._run_apply != s2p2._run_apply)
 
+    def test_NLBGS_Aitken(self):
+
+        prob = Problem()
+        model = prob.model = SellarDerivatives()
+        model.nonlinear_solver = NonlinearBlockGS()
+
+        prob.setup()
+        model.nonlinear_solver.options['use_aitken'] = True
+        prob.run_model()
+
+        assert_rel_error(self, prob['y1'], 25.58830273, .00001)
+        assert_rel_error(self, prob['y2'], 12.05848819, .00001)
+        self.assertTrue(model.nonlinear_solver._iter_count == 5)
 
 if __name__ == "__main__":
     unittest.main()

openmdao/solvers/tests/test_solver_features.py

@@ -63,6 +63,5 @@ def test_specify_subgroup_solvers(self):
         assert_rel_error(self, prob['g2.y1'], 0.64, .00001)
         assert_rel_error(self, prob['g2.y2'], 0.80, .00001)
 
-
 if __name__ == "__main__":
     unittest.main()

openmdao/vectors/default_vector.py

@@ -426,6 +426,26 @@ def set_const(self, val):
         for data in itervalues(self._data):
             data[:] = val
 
+    def dot(self, vec):
+        """
+        Compute the dot product of the real parts of the current vec and the incoming vec.
+
+        Parameters
+        ----------
+        vec : <Vector>
+            The incoming vector being dotted with self.
+
+        Returns
+        -------
+        float
+            The computed dot product value.
+        """
+        global_sum = 0
+        for set_name, data in iteritems(self._data):
+            global_sum += np.dot(self._data[set_name], vec._data[set_name])
+
+        return global_sum
+
     def get_norm(self):
         """
         Return the norm of this vector.

openmdao/vectors/tests/test_vector.py

@@ -1,6 +1,7 @@
 import unittest
 
 from openmdao.api import Problem, IndepVarComp
+from openmdao.parallel_api import PETScVector
 
 class TestVector(unittest.TestCase):
 
@@ -19,7 +20,35 @@ def test_iter(self):
 
         self.assertListEqual(outputs, expected, msg='Iter is not returning the expected names')
 
+    def test_dot(self):
+        p = Problem()
+        comp = IndepVarComp()
+        comp.add_output('v1', val=1.0)
+        comp.add_output('v2', val=2.0)
+        p.model.add_subsystem('des_vars', comp, promotes=['*'])
+        p.setup()
+        p.final_setup()
+
+        new_vec = p.model._outputs._clone()
+        new_vec.set_const(3.)
+
+        self.assertEqual(new_vec.dot(p.model._outputs), 9.)
+
+    def test_dot_petsc(self):
+        p = Problem()
+        comp = IndepVarComp()
+        comp.add_output('v1', val=1.0)
+        comp.add_output('v2', val=2.0)
+        p.model.add_subsystem('des_vars', comp, promotes=['*'])
+        p.setup(vector_class=PETScVector)
+        p.final_setup()
+
+        new_vec = p.model._outputs._clone()
+        new_vec.set_const(3.)
+
+        self.assertEqual(new_vec.dot(p.model._outputs), 9.)
+
+
 if __name__ == '__main__':
 
     unittest.main()
-

openmdao/vectors/vector.py

@@ -491,6 +491,19 @@ def set_const(self, val):
         """
         pass
 
+    def dot(self, vec):
+        """
+        Compute the dot product of the current vec and the incoming vec.
+
+        Must be implemented by the subclass.
+
+        Parameters
+        ----------
+        vec : <Vector>
+            The incoming vector being dotted with self.
+        """
+        pass
+
     def get_norm(self):
         """
         Return the norm of this vector.