# Source code for gemseo.mda.gauss_seidel

# -*- coding: utf-8 -*-
# Copyright 2021 IRT Saint Exupéry, https://www.irt-saintexupery.com
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# 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
#    OTHER AUTHORS   - MACROSCOPIC CHANGES
"""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())