# 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 gemseo.core.chain import MDOChain
from gemseo.core.discipline import MDODiscipline
from gemseo.mda.mda import MDA

[docs]class MDAGaussSeidel(MDA):
"""Perform a MDA analysis using a Gauss-Seidel algorithm, 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)
"""

def __init__(
self,
disciplines,
name=None,
max_mda_iter=10,
grammar_type=MDODiscipline.JSON_GRAMMAR_TYPE,
tolerance=1e-6,
linear_solver_tolerance=1e-12,
warm_start=False,
use_lu_fact=False,
over_relax_factor=1.0,
):
"""Constructor.

:param disciplines: the disciplines list
:type disciplines: list(MDODiscipline)
:param max_mda_iter: maximum number of iterations
:type max_mda_iter: int
:param name: the name of the chain
:type name: str
:param grammar_type: the type of grammar to use for IO declaration
either JSON_GRAMMAR_TYPE or SIMPLE_GRAMMAR_TYPE
:type grammar_type: str
:param tolerance: tolerance of the iterative direct coupling solver,
norm of the current residuals divided by initial residuals norm
shall be lower than the tolerance to stop iterating
:type tolerance: float
:param linear_solver_tolerance: Tolerance of the linear solver
:type linear_solver_tolerance: float
:param warm_start: if True, the second iteration and ongoing
start from the previous coupling solution
:type warm_start: bool
:param use_lu_fact: if True, when using adjoint/forward
differenciation, store a LU factorization of the matrix
to solve faster multiple RHS problem
:type use_lu_fact: bool
:param over_relax_factor: 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
:type over_relax_factor: float
"""
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,
)
assert over_relax_factor > 0.0
assert over_relax_factor <= 2.0
self.over_relax_factor = over_relax_factor
self._initialize_grammars()
self._set_default_inputs()
self._compute_input_couplings()

def _initialize_grammars(self):
"""Defines all inputs and outputs of the chain."""
# self.chain.initialize_grammars()
self.input_grammar.update_from(self.chain.input_grammar)
self.output_grammar.update_from(self.chain.output_grammar)

def _run(self):
"""Runs the disciplines in a sequential way until the difference between outputs
is under tolerance.

:returns: the local data
"""
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
)
# 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())