Migrate monolithic FSI to new solid time integration

src/adapter/4C_adapter_str_structure_new.cpp

@@ -820,21 +820,39 @@ void Adapter::StructureBaseAlgorithmNew::create_wrapper(
     case Core::ProblemType::thermo_fsi:
     case Core::ProblemType::fsi_xfem:
     {
+      const Teuchos::ParameterList& fsidyn = problem->fsi_dynamic_params();
+      auto coupling = Teuchos::getIntegralValue<FsiCoupling>(fsidyn, "COUPALGO");
+
+      if (Core::Communication::my_mpi_rank((actdis_->get_comm())) == 0)
+        Core::IO::cout << "Using StructureNOXCorrectionWrapper()..." << Core::IO::endl;
+
       // Are there any constraint conditions active?
       const std::set<Inpar::Solid::ModelType>& modeltypes =
           ti_strategy->get_data_sdyn().get_model_types();
       if (modeltypes.find(Inpar::Solid::model_lag_pen_constraint) != modeltypes.end())
       {
-        if (Core::Communication::my_mpi_rank((actdis_->get_comm())) == 0)
-          Core::IO::cout << "Using StructureNOXCorrectionWrapper()..." << Core::IO::endl;
-
         str_wrapper_ = std::make_shared<StructureConstrMerged>(
             std::make_shared<StructureNOXCorrectionWrapper>(ti_strategy));
       }
       else
       {
-        // case of partitioned fsi
-        str_wrapper_ = std::make_shared<FSIStructureWrapper>(ti_strategy);
+        // case of monolithic fsi
+        if (coupling == fsi_iter_mortar_monolithicfluidsplit or
+            coupling == fsi_iter_mortar_monolithicstructuresplit or
+            coupling == fsi_iter_monolithicfluidsplit or
+            coupling == fsi_iter_monolithicstructuresplit or
+            coupling == fsi_iter_sliding_monolithicfluidsplit or
+            coupling == fsi_iter_sliding_monolithicstructuresplit or
+            coupling == fsi_iter_mortar_monolithicfluidsplit_saddlepoint)
+        {
+          str_wrapper_ = std::make_shared<FSIStructureWrapper>(
+              std::make_shared<StructureNOXCorrectionWrapper>(ti_strategy));
+        }
+        else
+        {
+          // case of partitioned fsi
+          str_wrapper_ = std::make_shared<FSIStructureWrapper>(ti_strategy);
+        }
       }
       break;
     }

src/core/linalg/src/sparse/4C_linalg_blocksparsematrix.cpp

@@ -55,6 +55,76 @@ std::shared_ptr<Core::LinAlg::SparseMatrix> Core::LinAlg::BlockSparseMatrixBase:
   return sparse;
 }
 
+/*----------------------------------------------------------------------*
+ *----------------------------------------------------------------------*/
+std::shared_ptr<Core::LinAlg::BlockSparseMatrixBase>
+Core::LinAlg::copy_sparse_to_block_sparse_matrix(const Core::LinAlg::SparseMatrix& sparse_matrix,
+    const Core::LinAlg::MultiMapExtractor& domainmaps,
+    const Core::LinAlg::MultiMapExtractor& rangemaps)
+{
+  // Step 1: Precompute the number of max number of entries per row in the blocks
+
+  const int number_of_domain_maps = domainmaps.num_maps();
+
+  const auto& epetra_matrix = sparse_matrix.epetra_matrix();
+  const Map& sparse_matrix_col_map = sparse_matrix.col_map();
+  const int nummyrows = epetra_matrix.NumMyRows();
+  std::vector<int> max_num_entries_per_row_per_block(number_of_domain_maps, 0);
+  std::vector<int> num_entries_per_row_per_block(number_of_domain_maps, 0);
+  for (int iRow = 0; iRow < nummyrows; ++iRow)
+  {
+    int num_entries_per_row;
+    double* values;
+    int* indices;
+    epetra_matrix.ExtractMyRowView(iRow, num_entries_per_row, values, indices);
+    num_entries_per_row_per_block.assign(number_of_domain_maps, 0);
+    for (int iCol = 0; iCol < num_entries_per_row; ++iCol)
+    {
+      int col_gid = sparse_matrix_col_map.gid(indices[iCol]);
+      for (int m = 0; m < number_of_domain_maps; ++m)
+      {
+        if (domainmaps.map(m)->my_gid(col_gid))
+        {
+          ++num_entries_per_row_per_block[m];
+          break;
+        }
+      }
+    }
+    // check whether this works element wise ?!?!
+    max_num_entries_per_row_per_block =
+        std::max(num_entries_per_row_per_block, max_num_entries_per_row_per_block);
+  }
+
+  int max_num_entries_per_row = *std::ranges::max_element(
+      max_num_entries_per_row_per_block.begin(), max_num_entries_per_row_per_block.end());
+
+  // allocate block matrix with educated guess for number of non-zeros per row
+  auto block_matrix =
+      std::make_shared<Core::LinAlg::BlockSparseMatrix<Core::LinAlg::DefaultBlockMatrixStrategy>>(
+          domainmaps, rangemaps, max_num_entries_per_row, false, true);
+
+  // Step 2: Copy sparse matrix to block structure
+
+  for (int iRow = 0; iRow < nummyrows; ++iRow)
+  {
+    int num_entries_per_row;
+    double* values;
+    int* indices;
+
+    epetra_matrix.ExtractMyRowView(iRow, num_entries_per_row, values, indices);
+    const int row_gid = epetra_matrix.RowMap().GID(iRow);
+
+    for (int iCol = 0; iCol < num_entries_per_row; ++iCol)
+    {
+      int col_gid = sparse_matrix_col_map.gid(indices[iCol]);
+      block_matrix->assemble(values[iCol], row_gid, col_gid);
+    }
+  }
+
+  if (sparse_matrix.filled()) block_matrix->complete();
+
+  return block_matrix;
+}
 
 /*----------------------------------------------------------------------*
  *----------------------------------------------------------------------*/
@@ -348,6 +418,35 @@ const Epetra_Map& Core::LinAlg::BlockSparseMatrixBase::OperatorRangeMap() const
   return full_range_map().get_epetra_map();
 }
 
+void Core::LinAlg::BlockSparseMatrixBase::print(std::ostream& os) const
+{
+  for (int row = 0; row < rows(); row++)
+  {
+    for (int col = 0; col < cols(); col++)
+    {
+      for (int proc = 0; proc < Core::Communication::num_mpi_ranks(get_comm()); ++proc)
+      {
+        Core::Communication::barrier(get_comm());
+        if (proc == Core::Communication::my_mpi_rank(get_comm()))
+        {
+          if (matrix(row, col).num_my_nonzeros() == 0)
+          {
+            os << "\nBlockSparseMatrix row, col: " << row << ", " << col << " on rank "
+               << Core::Communication::my_mpi_rank(get_comm()) << " is empty" << std::endl;
+          }
+          else
+          {
+            os << "\nBlockSparseMatrix row, col: " << row << ", " << col << " on rank "
+               << Core::Communication::my_mpi_rank(get_comm()) << std::endl;
+
+            matrix(row, col).print(os);
+          }
+        }
+        Core::Communication::barrier(get_comm());
+      }
+    }
+  }
+}
 
 /*----------------------------------------------------------------------*
  *----------------------------------------------------------------------*/

src/core/linalg/src/sparse/4C_linalg_blocksparsematrix.hpp

@@ -18,6 +18,9 @@ FOUR_C_NAMESPACE_OPEN
 
 namespace Core::LinAlg
 {
+  // forward declaration
+  class DefaultBlockMatrixStrategy;
+
   /// Internal base class of BlockSparseMatrix that contains the non-template stuff
   /*!
 
@@ -210,6 +213,9 @@ namespace Core::LinAlg
     /// access to range map extractor in derived classes
     const MultiMapExtractor& range_extractor() const { return rangemaps_; }
 
+    //! Print to user-provided output stream
+    void print(std::ostream& os) const;
+
    private:
     /// the full domain map together with all partial domain maps
     MultiMapExtractor domainmaps_;
@@ -230,6 +236,11 @@ namespace Core::LinAlg
     bool usetranspose_;
   };
 
+  // split LinAlg::SparseMatrix into LinAlg::BlockSparseMatrix based on the given maps
+  std::shared_ptr<BlockSparseMatrixBase> copy_sparse_to_block_sparse_matrix(
+      const Core::LinAlg::SparseMatrix& sparse_matrix,
+      const Core::LinAlg::MultiMapExtractor& domainmaps,
+      const Core::LinAlg::MultiMapExtractor& rangemaps);
 
 
   /// Block matrix consisting of SparseMatrix blocks

src/fsi/src/monolithic/4C_fsi_monolithic.cpp

@@ -10,8 +10,10 @@
 #include "4C_adapter_ale.hpp"
 #include "4C_adapter_ale_fsi.hpp"
 #include "4C_adapter_fld_fluid_fsi.hpp"
+#include "4C_adapter_str_factory.hpp"
 #include "4C_adapter_str_fsi_timint_adaptive.hpp"
 #include "4C_adapter_str_fsiwrapper.hpp"
+#include "4C_adapter_str_structure_new.hpp"
 #include "4C_constraint_manager.hpp"
 #include "4C_coupling_adapter.hpp"
 #include "4C_fem_discretization.hpp"
@@ -55,6 +57,7 @@ FSI::MonolithicBase::MonolithicBase(MPI_Comm comm, const Teuchos::ParameterList&
       isadastructure_(false),
       isadafluid_(false),
       isadasolver_(false),
+      use_old_structure_(true),
       verbosity_(Teuchos::getIntegralValue<Inpar::FSI::Verbosity>(
           Global::Problem::instance()->fsi_dynamic_params(), "VERBOSITY"))
 {
@@ -98,18 +101,38 @@ void FSI::MonolithicBase::create_structure_time_integrator(
   // access structural dynamic params
   const Teuchos::ParameterList& sdyn = Global::Problem::instance()->structural_dynamic_params();
 
-  // ask base algorithm for the structural time integrator
-  std::shared_ptr<Adapter::StructureBaseAlgorithm> structure =
-      std::make_shared<Adapter::StructureBaseAlgorithm>(
-          timeparams, const_cast<Teuchos::ParameterList&>(sdyn), structdis);
-  structure_ =
-      std::dynamic_pointer_cast<Adapter::FSIStructureWrapper>(structure->structure_field());
-  structure_->setup();
+  if (Teuchos::getIntegralValue<Inpar::Solid::IntegrationStrategy>(sdyn, "INT_STRATEGY") ==
+      Inpar::Solid::int_standard)
+  {
+    use_old_structure_ = false;
 
-  if (structure_ == nullptr)
-    FOUR_C_THROW(
-        "Cast from Adapter::Structure to Adapter::FSIStructureWrapper "
-        "failed.");
+    std::shared_ptr<Adapter::StructureBaseAlgorithmNew> structure =
+        Adapter::build_structure_algorithm(sdyn);
+    structure->init(Global::Problem::instance()->fsi_dynamic_params(),
+        const_cast<Teuchos::ParameterList&>(sdyn), structdis);
+    structure->setup();
+
+    structure_ =
+        std::dynamic_pointer_cast<Adapter::FSIStructureWrapper>(structure->structure_field());
+
+    if (structure_ == nullptr)
+      FOUR_C_THROW("cast from ADAPTER::Structure to ADAPTER::FSIStructureWrapper failed");
+  }
+  else
+  {
+    // ask base algorithm for the structural time integrator
+    std::shared_ptr<Adapter::StructureBaseAlgorithm> structure =
+        std::make_shared<Adapter::StructureBaseAlgorithm>(
+            timeparams, const_cast<Teuchos::ParameterList&>(sdyn), structdis);
+    structure_ =
+        std::dynamic_pointer_cast<Adapter::FSIStructureWrapper>(structure->structure_field());
+    structure_->setup();
+
+    if (structure_ == nullptr)
+      FOUR_C_THROW(
+          "Cast from Adapter::Structure to Adapter::FSIStructureWrapper "
+          "failed.");
+  }
 
   return;
 }
@@ -199,7 +222,7 @@ void FSI::MonolithicBase::output()
   fluid_field()->output();
   ale_field()->output();
 
-  if (structure_field()->get_constraint_manager()->have_monitor())
+  if (use_old_structure_)
   {
     structure_field()->get_constraint_manager()->compute_monitor_values(
         *structure_field()->dispnp());
@@ -625,7 +648,7 @@ void FSI::Monolithic::time_step(const std::shared_ptr<NOX::Nln::Interface::Requi
   // initial guess. (The Dirichlet conditions are already build in!)
   std::shared_ptr<Core::LinAlg::Vector<double>> initial_guess_v =
       std::make_shared<Core::LinAlg::Vector<double>>(*dof_row_map(), true);
-  initial_guess(initial_guess_v);
+  if (use_old_structure_) initial_guess(initial_guess_v);
 
   NOX::Nln::Vector noxSoln(initial_guess_v, NOX::Nln::Vector::MemoryType::View);
 

src/fsi/src/monolithic/4C_fsi_monolithic.hpp

@@ -249,6 +249,8 @@ namespace FSI
     bool isadafluid_;      ///< Time step size adaptivity based on fluid?
     bool isadasolver_;     ///< Time step size adaptivity based on solver convergence?
 
+    bool use_old_structure_;  /// indicator if old or new time integration is used TEMPORARY
+
     //@}
 
     //! @name output related stuff

src/fsi/src/monolithic/model_evaluator/4C_fsi_monolithicfluidsplit.cpp

@@ -315,10 +315,13 @@ void FSI::MonolithicFluidSplit::setup_rhs_residual(Core::LinAlg::Vector<double>&
   const double fluidscale = fluid_field()->residual_scaling();
 
   // get single field residuals
-  const Core::LinAlg::Vector<double> sv(*structure_field()->rhs());
+  Core::LinAlg::Vector<double> sv(*structure_field()->rhs());
   const Core::LinAlg::Vector<double> fv(*fluid_field()->rhs());
   const Core::LinAlg::Vector<double> av(*ale_field()->rhs());
 
+  // NOX treats rhs different in new solid time integration, hence sign inversion is necessary here
+  if (!use_old_structure_) sv.scale(-1.0);
+
   // extract only inner DOFs from fluid (=slave) and ALE field
   std::shared_ptr<Core::LinAlg::Vector<double>> fov =
       fluid_field()->interface()->extract_other_vector(fv);
@@ -1273,12 +1276,15 @@ void FSI::MonolithicFluidSplit::output()
 
   ale_field()->output();
 
-  if (structure_field()->get_constraint_manager()->have_monitor())
+  if (use_old_structure_)
   {
-    structure_field()->get_constraint_manager()->compute_monitor_values(
-        *structure_field()->dispnp());
-    if (Core::Communication::my_mpi_rank(get_comm()) == 0)
-      structure_field()->get_constraint_manager()->print_monitor_values();
+    if (structure_field()->get_constraint_manager()->have_monitor())
+    {
+      structure_field()->get_constraint_manager()->compute_monitor_values(
+          *structure_field()->dispnp());
+      if (Core::Communication::my_mpi_rank(get_comm()) == 0)
+        structure_field()->get_constraint_manager()->print_monitor_values();
+    }
   }
 }