[CoSim] Adding a strong block solver and block MVQN acceleration

.gitignore

@@ -86,3 +86,6 @@ remote_*.md
 
 # External downloaded libs (TPL downloaded libs, etc...)
 external_libraries/tetgen/*
+
+# Mac folder info files
+*.DS_Store

applications/CoSimulationApplication/python_scripts/convergence_accelerators/blockMvqn.py

@@ -0,0 +1,128 @@
+## @module blockMvqn
+# This module contains the class MVQNConvergenceAccelerator
+# Author: Tiba Azzeddine
+# Date: Nov. 06, 2023
+
+# Importing the Kratos Library
+import KratosMultiphysics as KM
+
+# Importing the base class
+from KratosMultiphysics.CoSimulationApplication.base_classes.co_simulation_convergence_accelerator import CoSimulationConvergenceAccelerator
+
+# CoSimulation imports
+import KratosMultiphysics.CoSimulationApplication.co_simulation_tools as cs_tools
+
+# Other imports
+import numpy as np
+from copy import deepcopy
+from collections import deque
+
+
+def Create(settings):
+    cs_tools.SettingsTypeCheck(settings)
+    return BLOCKMVQNConvergenceAccelerator(settings)
+
+## Class BLOCKMVQNConvergenceAccelerator.
+# This class contains the implementation of the Block MVQN method and helper functions.
+# Reference: A.E.J. Bogaers et al. "Quasi-Newton methods for implicit black-box FSI coupling", Computational methods in applied mechanics and engineering. 279(2014) 113-132.
+class BLOCKMVQNConvergenceAccelerator(CoSimulationConvergenceAccelerator):
+    ## The constructor.
+    # @param horizon Maximum number of vectors to be stored in each time step.
+    # @param alpha Relaxation factor for computing the update, when no vectors available.
+    def __init__( self, settings):
+        super().__init__(settings)
+
+        horizon = self.settings["horizon"].GetInt()
+        self.alpha = self.settings["alpha"].GetDouble()
+        self.epsilon = self.settings["epsilon"].GetDouble()
+
+        self.X_tilde = {}
+        self.X = {}
+        self.J = {}
+        self.J_hat = {}
+        self.coupl_data_names = {}
+
+        for solverData in settings["solver_sequence"].values():
+            self.X_tilde[solverData["data_name"].GetString()] = deque( maxlen = horizon )
+            self.X[solverData["data_name"].GetString()] = deque( maxlen = horizon )
+            self.J[solverData["data_name"].GetString()] = None # size will be determined when first time get the input vector
+            self.J_hat[solverData["data_name"].GetString()] = None
+            self.coupl_data_names[solverData["data_name"].GetString()] = solverData["coupled_data_name"].GetString()
+
+    ## UpdateSolution(r, x, y, data_name, yResidual)
+    # @param r residual r_k
+    # @param x solution x_k
+    # @param y (coupled solver) solution y_k
+    # @param data_name coupling variable
+    # @param yResidual (coupled solver) residual yResidual
+    # Computes the approximated update in each iteration.
+    def UpdateSolution( self, r, x, y, data_name, yResidual,):
+
+        coupled_data_name = self.coupl_data_names[data_name]
+        self.X_tilde[data_name].appendleft( deepcopy(r + x) )
+        self.X[coupled_data_name].appendleft( deepcopy(y) )
+
+        col = len(self.X[coupled_data_name]) - 1
+        row = len(r)
+        rowY = len(y)
+        k = col
+        if self.echo_level > 3:
+            cs_tools.cs_print_info(self._ClassName(), "Number of new modes: ", col )
+
+        ## For the first iteration
+        if k == 0:
+            if self.J[data_name] is None or self.J[coupled_data_name] is None:
+                ## Zero initial Jacobians
+                self.J_hat[data_name] = np.zeros( (row, rowY) )
+                self.J_hat[coupled_data_name] = np.zeros( (rowY, row) )
+                self.J[coupled_data_name] = np.zeros( (rowY, row) )
+                self.J[data_name] = np.zeros( (row, rowY) )
+                return self.alpha * r # Initial acceleration to be a constant relaxation
+            else:
+
+                blockJacobian = (np.eye(row) - self.J[data_name] @ self.J[coupled_data_name])
+                b = r - self.J[data_name] @ yResidual
+                return np.linalg.solve( blockJacobian, b )
+
+        ## Construct matrix W(differences of intermediate solutions x)
+        W = np.empty( shape = (col, rowY) ) # will be transposed later
+        for i in range(0, col):
+            W[i] = self.X[coupled_data_name][i] - self.X[coupled_data_name][i + 1]
+        W = W.T
+
+        ## Construct matrix W(differences of intermediate solutions y~)
+        V = np.empty( shape = (col, row) ) # will be transposed later
+        for i in range(0, col):
+            V[i] = self.X_tilde[data_name][i] - self.X_tilde[data_name][i + 1]
+        V = V.T
+
+        self.J_hat[data_name] = self.J[data_name] + (V - self.J[data_name] @ W) @ (np.linalg.pinv(W, rcond=self.epsilon))
+
+        blockJacobian = (np.eye(row) - self.J_hat[data_name] @ self.J_hat[coupled_data_name])
+        b = r - self.J_hat[data_name] @ yResidual
+
+        return np.linalg.solve( blockJacobian, b )
+
+    def FinalizeSolutionStep( self ):
+
+        ## Assign J=J_hat
+        for _, data_name in enumerate(self.J.keys()):
+            self.J[data_name] = self.J_hat[data_name].copy()
+
+            ## Clear the buffer
+            self.X_tilde[data_name].clear()
+            self.X[data_name].clear()
+
+        if self.echo_level > 3:
+            cs_tools.cs_print_info(self._ClassName(), "Jacobian matrix updated!")
+
+    @classmethod
+    def _GetDefaultParameters(cls):
+        this_defaults = KM.Parameters("""{
+            "horizon" : 15,
+            "alpha"   : 1.0,
+            "epsilon" : 1e-7,
+            "solver_sequence" : []
+        }""")
+        this_defaults.AddMissingParameters(super()._GetDefaultParameters())
+        return this_defaults

applications/CoSimulationApplication/python_scripts/convergence_accelerators/convergence_accelerator_wrapper.py

@@ -20,7 +20,20 @@ def __init__(self,
                  settings: KratosMultiphysics.Parameters,
                  interface_data_dict: "dict[str,CouplingInterfaceData]",
                  parent_coupled_solver_data_communicator: KratosMultiphysics.DataCommunicator):
-        self.interface_data = interface_data_dict[settings["data_name"].GetString()]
+        if settings.Has("solver_sequence"):
+            self.interface_data_dict = interface_data_dict
+            self.solver_names_sequence = []
+            self.coupl_data_names = {}
+            self.input_data = {}
+            for _, data_name in enumerate(self.interface_data_dict.keys()):
+                self.solver_names_sequence.append(self.interface_data_dict[data_name].solver_name)
+            self.isBlockAccelerator = True
+            for solverData in settings["solver_sequence"].values():
+                self.coupl_data_names[solverData["data_name"].GetString()] = solverData["coupled_data_name"].GetString()
+            self.yResidual_computation = CreateResidualComputation(settings, interface_data_dict)
+        else:
+            self.isBlockAccelerator = False
+            self.interface_data = interface_data_dict[settings["data_name"].GetString()]
         self.residual_computation = CreateResidualComputation(settings, interface_data_dict)
 
         # Remove extra entries from accelerator parameters
@@ -31,16 +44,36 @@ def __init__(self,
         self.conv_acc = CreateConvergenceAccelerator(settings)
         self.data_communicator = parent_coupled_solver_data_communicator
 
-        if self.interface_data.IsDefinedOnThisRank():
-            conv_acc_supports_dist_data = self.conv_acc.SupportsDistributedData()
-            self.executing_rank = conv_acc_supports_dist_data or (self.interface_data.GetModelPart().GetCommunicator().MyPID() == 0)
-            self.gather_scatter_required = self.interface_data.IsDistributed() and not conv_acc_supports_dist_data
-            if self.gather_scatter_required:
-                self.data_comm = self.interface_data.GetModelPart().GetCommunicator().GetDataCommunicator()
-                self.sizes_from_ranks = np.cumsum(self.data_comm.GatherInts([self.interface_data.Size()], 0))
+        if not self.isBlockAccelerator:
+            if self.interface_data.IsDefinedOnThisRank():
+                conv_acc_supports_dist_data = self.conv_acc.SupportsDistributedData()
+                self.executing_rank = conv_acc_supports_dist_data or (self.interface_data.GetModelPart().GetCommunicator().MyPID() == 0)
+                self.gather_scatter_required = self.interface_data.IsDistributed() and not conv_acc_supports_dist_data
+                if self.gather_scatter_required:
+                    self.data_comm = self.interface_data.GetModelPart().GetCommunicator().GetDataCommunicator()
+                    self.sizes_from_ranks = np.cumsum(self.data_comm.GatherInts([self.interface_data.Size()], 0))
+        else:
+            self.executing_rank = {}
+            self.gather_scatter_required = {}
+            self.data_comm = {}
+            self.sizes_from_ranks = {}
+            for _, data_name in enumerate(interface_data_dict.keys()):
+                if self.interface_data_dict[data_name].IsDefinedOnThisRank():
+                    conv_acc_supports_dist_data = self.conv_acc.SupportsDistributedData()
+                    self.executing_rank[data_name] = conv_acc_supports_dist_data or (self.interface_data_dict[data_name].GetModelPart().GetCommunicator().MyPID() == 0)
+                    self.gather_scatter_required[data_name] = self.interface_data_dict[data_name].IsDistributed() and not conv_acc_supports_dist_data
+                    if self.gather_scatter_required[data_name]:
+                        self.data_comm[data_name] = self.interface_data_dict[data_name].GetModelPart().GetCommunicator().GetDataCommunicator()
+                        self.sizes_from_ranks[data_name] = np.cumsum(self.data_comm[data_name].GatherInts([self.interface_data_dict[data_name].Size()], 0))
 
     def Initialize(self):
         self.conv_acc.Initialize()
+        if self.isBlockAccelerator:
+            self.prev_input_data = {}
+            self.output_data = {}
+            for _, data_name in enumerate(self.interface_data_dict.keys()):
+                self.prev_input_data[data_name] = self.interface_data_dict[data_name].GetData()
+                self.output_data[data_name] = self.interface_data_dict[data_name].GetData()
 
     def Finalize(self):
         self.conv_acc.Finalize()
@@ -52,17 +85,23 @@ def FinalizeSolutionStep(self):
         self.conv_acc.FinalizeSolutionStep()
 
     def InitializeNonLinearIteration(self):
-        if self.interface_data.IsDefinedOnThisRank():
-            # Saving the previous data for the computation of the residual
-            # and the computation of the solution update
-            self.input_data = self.interface_data.GetData()
+        # Saving the previous data for the computation of the residual
+        # and the computation of the solution update
+        if self.isBlockAccelerator:
+            for _, data_name in enumerate(self.interface_data_dict.keys()):
+                interface_data = self.interface_data_dict[data_name]
+                if interface_data.IsDefinedOnThisRank():
+                    self.input_data[data_name] = interface_data.GetData()
+        else:
+            if self.interface_data.IsDefinedOnThisRank():
+                self.input_data = self.interface_data.GetData()
 
         self.conv_acc.InitializeNonLinearIteration()
 
     def FinalizeNonLinearIteration(self):
         self.conv_acc.FinalizeNonLinearIteration()
 
-    def ComputeAndApplyUpdate(self):
+    def _ComputeAndApplyUpdate(self,):
         if not self.interface_data.IsDefinedOnThisRank(): return
 
         residual = self.residual_computation.ComputeResidual(self.input_data)
@@ -85,6 +124,60 @@ def ComputeAndApplyUpdate(self):
 
         self.interface_data.SetData(updated_data)
 
+    def _ComputeAndApplyBlockUpdate(self, solver_name = None):
+        # Defaulting to the first solver in the sequence
+        if solver_name is None:
+            solver_name = self.solver_names_sequence[0]
+
+        # Retrieving solver and data names
+        solver_id = self.solver_names_sequence.index(solver_name)
+        data_name, interface_data = list(self.interface_data_dict.items())[solver_id]
+        interface_data = self.interface_data_dict[data_name]
+        coupled_data_name = self.coupl_data_names[data_name]
+
+        # Data Comm
+        executing_rank = self.executing_rank[data_name]
+        gather_scatter_required = self.gather_scatter_required[data_name]
+
+        # Retrieving data
+        input_data = self.input_data[data_name]
+        self.prev_input_data[data_name] = interface_data.GetData()
+        prev_input_data = self.prev_input_data[coupled_data_name]
+
+        if not interface_data.IsDefinedOnThisRank(): return
+
+        residual = self.residual_computation.ComputeResidual(input_data, data_name)
+        yResidual = self.yResidual_computation.ComputeResidual(prev_input_data, coupled_data_name)
+
+        input_data_for_acc = input_data
+        input_other_solver = self.output_data[coupled_data_name]
+
+        if gather_scatter_required:
+            residual = np.array(np.concatenate(self.data_comm[data_name].GathervDoubles(residual, 0)))
+            yResidual = np.array(np.concatenate(self.data_comm[coupled_data_name].GathervDoubles(yResidual, 0)))
+            input_data_for_acc = np.array(np.concatenate(self.data_comm[data_name].GathervDoubles(input_data_for_acc, 0)))
+            input_other_solver = np.array(np.concatenate(self.data_comm[coupled_data_name].GathervDoubles(input_other_solver, 0)))
+
+        if executing_rank:
+            updated_data = input_data_for_acc + self.conv_acc.UpdateSolution(residual, input_data_for_acc, input_other_solver, data_name, yResidual)
+
+        if gather_scatter_required:
+            if executing_rank:
+                data_to_scatter = np.split(updated_data, self.sizes_from_ranks[data_name][:-1])
+            else:
+                data_to_scatter = []
+
+            updated_data = self.data_comm[data_name].ScattervDoubles(data_to_scatter, 0)
+
+        interface_data.SetData(updated_data)
+        self.output_data[data_name] = updated_data
+
+    def ComputeAndApplyUpdate(self, solver_name = None):
+        if self.isBlockAccelerator:
+            self._ComputeAndApplyBlockUpdate(solver_name)
+        else:
+            self._ComputeAndApplyUpdate()
+
     def PrintInfo(self):
         self.conv_acc.PrintInfo()
 
@@ -93,16 +186,24 @@ def Check(self):
 
 class ConvergenceAcceleratorResidual(metaclass=ABCMeta):
     @abstractmethod
-    def ComputeResidual(self, input_data): pass
+    def ComputeResidual(self, input_data, data_name=None): pass
 
 class DataDifferenceResidual(ConvergenceAcceleratorResidual):
     def __init__(self,
                  settings: KratosMultiphysics.Parameters,
                  interface_data_dict: "dict[str,CouplingInterfaceData]"):
-        self.interface_data = interface_data_dict[settings["data_name"].GetString()]
-
-    def ComputeResidual(self, input_data):
-        return self.interface_data.GetData() - input_data
+        if settings.Has("solver_sequence"):
+            self.interface_data_dict = interface_data_dict
+            self.isBlockResidualComputer = True
+        else:
+            self.isBlockResidualComputer = False
+            self.interface_data = interface_data_dict[settings["data_name"].GetString()]
+
+    def ComputeResidual(self, input_data, data_name=None):
+        if self.isBlockResidualComputer:
+            return self.interface_data_dict[data_name].GetData() - input_data
+        else:
+            return self.interface_data.GetData() - input_data
 
 class DifferentDataDifferenceResidual(ConvergenceAcceleratorResidual):
     def __init__(self,