Source code for gemseo.mda.gauss_seidel

# -*- coding: utf-8 -*-
# Copyright 2021 IRT Saint Exupéry,
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
# License version 3 as published by the Free Software Foundation.
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# Lesser General Public License for more details.
# You should have received a copy of the GNU Lesser General Public License
# along with this program; if not, write to the Free Software Foundation,
# Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.

# Contributors:
#    INITIAL AUTHORS - API and implementation and/or documentation
#        :author: Francois Gallard
"""A Gauss Seidel algorithm for solving MDAs."""
from __future__ import division, unicode_literals

from typing import Any, Mapping, Optional, Sequence

from gemseo.core.chain import MDOChain
from gemseo.core.coupling_structure import MDOCouplingStructure
from gemseo.core.discipline import MDODiscipline
from gemseo.mda.mda import MDA

[docs]class MDAGaussSeidel(MDA): """An MDA analysis based on the Gauss-Seidel algorithm. This algorithm is an iterative technique to solve the linear system: .. math:: Ax = b by decomposing the matrix :math:`A` into the sum of a lower triangular matrix :math:`L_*` and a strictly upper triangular matrix :math:`U`. The new iterate is given by: .. math:: x_{k+1} = L_*^{-1}(b-Ux_k) """ _ATTR_TO_SERIALIZE = MDA._ATTR_TO_SERIALIZE + ( "strong_couplings", "over_relax_factor", "normed_residual", ) def __init__( self, disciplines, # type: Sequence[MDODiscipline] name=None, # type: Optional[str] max_mda_iter=10, # type: int grammar_type=MDODiscipline.JSON_GRAMMAR_TYPE, # type: str tolerance=1e-6, # type: float linear_solver_tolerance=1e-12, # type: float warm_start=False, # type: bool use_lu_fact=False, # type: bool over_relax_factor=1.0, # type: float coupling_structure=None, # type: Optional[MDOCouplingStructure] log_convergence=False, # type: bool linear_solver="DEFAULT", # type: str linear_solver_options=None, # type: Mapping[str,Any] ): # type: (...) -> None """ Args: over_relax_factor: The relaxation coefficient, used to make the method more robust, if ``0<over_relax_factor<1`` or faster if ``1<over_relax_factor<=2``. If ``over_relax_factor =1.``, it is deactivated. """ self.chain = MDOChain(disciplines, grammar_type=grammar_type) super(MDAGaussSeidel, self).__init__( disciplines, max_mda_iter=max_mda_iter, name=name, grammar_type=grammar_type, tolerance=tolerance, linear_solver_tolerance=linear_solver_tolerance, warm_start=warm_start, use_lu_fact=use_lu_fact, coupling_structure=coupling_structure, log_convergence=log_convergence, linear_solver=linear_solver, linear_solver_options=linear_solver_options, ) assert over_relax_factor > 0.0 assert over_relax_factor <= 2.0 self.over_relax_factor = over_relax_factor self._set_default_inputs() self._compute_input_couplings() def _initialize_grammars(self): self.input_grammar.update_from(self.chain.input_grammar) self.output_grammar.update_from(self.chain.output_grammar) def _run(self): # Run the disciplines in a sequential way # until the difference between outputs is under tolerance. if self.warm_start: self._couplings_warm_start() current_couplings = 0.0 relax = self.over_relax_factor use_relax = relax != 1.0 # store initial residual current_iter = 0 while not self._termination(current_iter) or current_iter == 0: for discipline in self.disciplines: discipline.execute(self.local_data) outs = discipline.get_output_data() if use_relax: # First time this output is computed, update directly local data self.local_data.update( {k: v for k, v in outs.items() if k not in self.local_data} ) # The couplings already exist in the local data, # so the over relaxation can be applied self.local_data.update( { k: relax * v + (1.0 - relax) * self.local_data[k] for k, v in outs.items() if k in self.local_data } ) else: self.local_data.update(outs) new_couplings = self._current_strong_couplings() self._compute_residual( current_couplings, new_couplings, current_iter, first=current_iter == 0, log_normed_residual=self.log_convergence, ) # store current residuals current_iter += 1 current_couplings = new_couplings for discipline in self.disciplines: # Update all outputs without relax self.local_data.update(discipline.get_output_data())