Source code for gemseo.mda.jacobi

# 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
# 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
# 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 Jacobi algorithm for solving MDAs."""

from __future__ import annotations

from typing import TYPE_CHECKING
from typing import Any
from typing import ClassVar

from gemseo.core.parallel_execution.disc_parallel_execution import DiscParallelExecution
from gemseo.mda.base_mda import BaseProcessFlow
from gemseo.mda.base_mda import _BaseMDAProcessFlow
from gemseo.mda.base_mda_solver import BaseMDASolver
from gemseo.mda.jacobi_settings import MDAJacobi_Settings
from gemseo.utils.constants import READ_ONLY_EMPTY_DICT

if TYPE_CHECKING:
    from collections.abc import Sequence
    from typing import ClassVar

    from gemseo.core.coupling_structure import DependencyGraph
    from gemseo.core.discipline import Discipline
    from gemseo.typing import StrKeyMapping


class _ProcessFlow(_BaseMDAProcessFlow):
    """The process data and execution flow."""

    def _get_disciplines_couplings(
        self, graph: DependencyGraph
    ) -> list[tuple[str, str, list[str]]]:
        couplings_results = []
        for disc in self._node.disciplines:
            in_data = graph.graph.get_edge_data(self._node, disc)
            if in_data:
                couplings_results.append((
                    self._node,
                    disc,
                    sorted(in_data["io"]),
                ))
            out_data = graph.graph.get_edge_data(disc, self._node)
            if out_data:
                couplings_results.append((
                    disc,
                    self._node,
                    sorted(out_data["io"]),
                ))

        return couplings_results


[docs] class MDAJacobi(BaseMDASolver): r"""Perform an MDA using the Jacobi algorithm. This algorithm is a fixed point iteration method to solve systems of non-linear equations of the form, .. math:: \left\{ \begin{matrix} F_1(x_1, x_2, \dots, x_n) = 0 \\ F_2(x_1, x_2, \dots, x_n) = 0 \\ \vdots \\ F_n(x_1, x_2, \dots, x_n) = 0 \end{matrix} \right. Beginning with :math:`x_1^{(0)}, \dots, x_n^{(0)}`, the iterates are obtained as the solution of the following :math:`n` **independent** non-linear equations: .. math:: \left\{ \begin{matrix} r_1\left( x_1^{(i+1)} \right) = F_1(x_1^{(i+1)}, x_2^{(i)}, \dots, x_n^{(i)}) = 0 \\ r_2\left( x_2^{(i+1)} \right) = F_2(x_1^{(i)}, x_2^{(i+1)}, \dots, x_n^{(i)}) = 0 \\ \vdots \\ r_n\left( x_n^{(i+1)} \right) = F_n(x_1^{(i)}, x_2^{(i)}, \dots, x_n^{(i+1)}) = 0 \end{matrix} \right. """ Settings: ClassVar[type[MDAJacobi_Settings]] = MDAJacobi_Settings """The pydantic model for the settings.""" parallel_execution: DiscParallelExecution | None """Either an executor of disciplines in parallel or ``None`` in serial mode.""" settings: MDAJacobi_Settings """The settings of the MDA""" _process_flow_class: ClassVar[type[BaseProcessFlow]] = _ProcessFlow def __init__( # noqa: D107 self, disciplines: Sequence[Discipline], settings_model: MDAJacobi_Settings | None = None, **settings: Any, ) -> None: super().__init__(disciplines, settings_model=settings_model, **settings) self._compute_input_coupling_names() self._set_resolved_variables(self._input_couplings) if self.settings.n_processes == 1: self._execute_disciplines = self._execute_disciplines_sequentially self.parallel_execution = None else: self._execute_disciplines = self._execute_disciplines_in_parallel self.parallel_execution = DiscParallelExecution( disciplines, self.settings.n_processes, self.settings.use_threading, exceptions_to_re_raise=(ValueError,), )
[docs] def get_process_flow(self) -> BaseProcessFlow: # noqa: D102 process_flow = super().get_process_flow() process_flow.is_parallel = self.settings.n_processes > 1 return process_flow
def _compute_input_coupling_names(self) -> None: """Compute the coupling variables that are inputs of the MDA. This must be overloaded here because the Jacobi algorithm induces a delay between the couplings, the strong couplings may be fully resolved but the weak ones may need one more iteration. The base MDA class uses strong couplings only, which is not satisfying here if some disciplines are not strongly coupled. """ if len(self.coupling_structure.strongly_coupled_disciplines) == len( self.disciplines ): return super()._compute_input_coupling_names() self._input_couplings = sorted( set(self.coupling_structure.all_couplings).intersection( self.io.input_grammar.names ) ) self._numeric_input_couplings = sorted( set(self._input_couplings).difference(self._non_numeric_array_variables) ) return None def _execute_disciplines_in_parallel(self) -> None: """Execute the disciplines in parallel.""" self.parallel_execution.execute([self.io.data] * len(self.disciplines)) def _execute_disciplines_sequentially(self) -> None: """Execute the disciplines sequentially.""" for discipline in self.disciplines: discipline.execute(self.io.data) def _execute_disciplines_and_update_local_data( self, input_data: StrKeyMapping = READ_ONLY_EMPTY_DICT ) -> None: self._execute_disciplines() for discipline in self.disciplines: self.io.data.update(discipline.io.get_output_data()) def _execute(self) -> None: super()._execute() while True: local_data_before_execution = self.io.data.copy() self._execute_disciplines_and_update_local_data() self._compute_residuals(local_data_before_execution) if self._stop_criterion_is_reached: break updated_couplings = self._sequence_transformer.compute_transformed_iterate( self.get_current_resolved_variables_vector(), self.get_current_resolved_residual_vector(), ) self._update_local_data_from_array(updated_couplings)