Move remaining scatra to teko

src/core/linear_solver/src/method/4C_linear_solver_method.hpp

@@ -62,8 +62,7 @@ namespace Core::LinearSolver
   {
     ilu,              ///< incomplete LU factorization with fill in levels (Ifpack package)
     multigrid_muelu,  ///< multigrid preconditioner (MueLu package, recommended!)
-    multigrid_nxn,  ///< multigrid preconditioner for a nxn block matrix (indirectly MueLu package)
-    block_teko      ///< block preconditioning (Teko package, recommended!)
+    block_teko        ///< block preconditioning (Teko package, recommended!)
   };
 
   /// linear solver type base class
@@ -90,11 +89,10 @@ namespace Core::LinearSolver
     virtual int solve() = 0;
 
     /// return number of iterations performed by solver
-    virtual int get_num_iters() const
+    [[nodiscard]] virtual int get_num_iters() const
     {
       FOUR_C_THROW("Not implemented in base class!");
-      return -1;
-    };
+    }
   };
 }  // namespace Core::LinearSolver
 

src/core/linear_solver/src/method/4C_linear_solver_method_input.cpp

@@ -47,7 +47,6 @@ namespace Core::LinearSolver
             {
                 {"ILU", Core::LinearSolver::PreconditionerType::ilu},
                 {"MueLu", Core::LinearSolver::PreconditionerType::multigrid_muelu},
-                {"AMGnxn", Core::LinearSolver::PreconditionerType::multigrid_nxn},
                 {"Teko", Core::LinearSolver::PreconditionerType::block_teko},
             },
             {.description =
@@ -129,15 +128,6 @@ namespace Core::LinearSolver
         // user-given name of solver block (just for beauty)
         parameter<std::string>("NAME",
             {.description = "User specified name for solver block", .default_value = "No_name"}),
-
-        // Parameters for AMGnxn Preconditioner
-        parameter<std::string>("AMGNXN_TYPE",
-            {.description = "Name of the pre-built preconditioner to be used. If set "
-                            "to\"XML\" the preconditioner is defined using a xml file",
-                .default_value = "AMG(BGS)"}),
-
-        Core::IO::InputSpecBuilders::parameter<std::optional<std::filesystem::path>>(
-            "AMGNXN_XML_FILE", {.description = "xml file defining the AMGnxn preconditioner"}),
     });
   }
 

src/core/linear_solver/src/method/4C_linear_solver_method_linalg.cpp

@@ -51,7 +51,7 @@ int Core::LinAlg::Solver::get_num_iters() const { return solver_->get_num_iters(
 /*----------------------------------------------------------------------*
  *----------------------------------------------------------------------*/
 void Core::LinAlg::Solver::adapt_tolerance(
-    const double desirednlnres, const double currentnlnres, const double better)
+    const double desirednlnres, const double currentnlnres, const double better) const
 {
   if (!params().isSublist("Belos Parameters")) FOUR_C_THROW("Adaptive tolerance only for Belos.");
 
@@ -112,7 +112,7 @@ void Core::LinAlg::Solver::adapt_tolerance(
 
 /*----------------------------------------------------------------------*
  *----------------------------------------------------------------------*/
-void Core::LinAlg::Solver::set_tolerance(const double tolerance)
+void Core::LinAlg::Solver::set_tolerance(const double tolerance) const
 {
   if (!params().isSublist("Belos Parameters"))
     FOUR_C_THROW("Set tolerance of linear solver only for Belos solver.");
@@ -123,15 +123,15 @@ void Core::LinAlg::Solver::set_tolerance(const double tolerance)
   if (!have_saved_value)
   {
     solver_params.set<double>(
-        "Convergence Tolerance Saved", solver_params.get<double>("Convergence Tolerance", 1.0e-8));
+        "Convergence Tolerance Saved", solver_params.get<double>("Convergence Tolerance"));
   }
 
   solver_params.set<double>("Convergence Tolerance", tolerance);
 }
 
 /*----------------------------------------------------------------------*
  *----------------------------------------------------------------------*/
-void Core::LinAlg::Solver::reset_tolerance()
+void Core::LinAlg::Solver::reset_tolerance() const
 {
   if (!params().isSublist("Belos Parameters")) return;
 
@@ -191,7 +191,9 @@ void Core::LinAlg::Solver::setup(std::shared_ptr<Epetra_Operator> matrix,
           solvertype);
     }
     else
+    {
       FOUR_C_THROW("Unknown type of solver");
+    }
   }
 
   solver_->setup(matrix, x, b, refactor, params.reset, params.projector);
@@ -290,6 +292,12 @@ Teuchos::ParameterList translate_four_c_to_belos(const Teuchos::ParameterList& i
   if (xmlfile)
   {
     beloslist.set("SOLVER_XML_FILE", xmlfile->string());
+
+    // required for adaptive linear solver tolerance
+    beloslist.set("Convergence Tolerance", inparams.get<double>("AZTOL"));
+    beloslist.set("Implicit Residual Scaling",
+        Belos::convertScaleTypeToString(
+            Teuchos::getIntegralValue<Belos::ScaleType>(inparams, "AZCONV")));
   }
   else
   {
@@ -341,7 +349,6 @@ Teuchos::ParameterList translate_four_c_to_belos(const Teuchos::ParameterList& i
         FOUR_C_THROW("Flag '{}'! \nUnknown solver for Belos.",
             Teuchos::getIntegralValue<Core::LinearSolver::IterativeSolverType>(
                 inparams, "AZSOLVE"));
-        break;
       }
     }
   }
@@ -358,15 +365,11 @@ Teuchos::ParameterList translate_four_c_to_belos(const Teuchos::ParameterList& i
     case Core::LinearSolver::PreconditionerType::multigrid_muelu:
       beloslist.set("Preconditioner Type", "ML");
       break;
-    case Core::LinearSolver::PreconditionerType::multigrid_nxn:
-      beloslist.set("Preconditioner Type", "AMGnxn");
-      break;
     case Core::LinearSolver::PreconditionerType::block_teko:
       beloslist.set("Preconditioner Type", "Teko");
       break;
     default:
       FOUR_C_THROW("Unknown preconditioner for Belos");
-      break;
   }
 
   // set parameters for Ifpack if used
@@ -387,14 +390,6 @@ Teuchos::ParameterList translate_four_c_to_belos(const Teuchos::ParameterList& i
     Teuchos::ParameterList& tekolist = outparams.sublist("Teko Parameters");
     tekolist = translate_four_c_to_teko(inparams, &beloslist);
   }
-  if (azprectype == Core::LinearSolver::PreconditionerType::multigrid_nxn)
-  {
-    Teuchos::ParameterList& amgnxnlist = outparams.sublist("AMGnxn Parameters");
-    auto amgnxn_xml = inparams.get<std::optional<std::filesystem::path>>("AMGNXN_XML_FILE");
-    amgnxnlist.set("AMGNXN_XML_FILE", amgnxn_xml);
-    std::string amgnxn_type = inparams.get<std::string>("AMGNXN_TYPE");
-    amgnxnlist.set<std::string>("AMGNXN_TYPE", amgnxn_type);
-  }
 
   return outparams;
 }
@@ -415,10 +410,8 @@ Teuchos::ParameterList Core::LinAlg::Solver::translate_solver_parameters(
   switch (Teuchos::getIntegralValue<Core::LinearSolver::SolverType>(inparams, "SOLVER"))
   {
     case Core::LinearSolver::SolverType::undefined:
-      std::cout << "undefined solver! Set " << inparams.name() << "  in your input file!"
-                << std::endl;
+      std::cout << "undefined solver! Set " << inparams.name() << "  in your input file!" << '\n';
       FOUR_C_THROW("fix your input file");
-      break;
     case Core::LinearSolver::SolverType::umfpack:
       outparams.set("solver", "umfpack");
       break;
@@ -430,7 +423,6 @@ Teuchos::ParameterList Core::LinAlg::Solver::translate_solver_parameters(
       break;
     default:
       FOUR_C_THROW("Unsupported type of solver");
-      break;
   }
 
   return outparams;

src/core/linear_solver/src/method/4C_linear_solver_method_linalg.hpp

@@ -148,9 +148,9 @@ namespace Core::LinAlg
     void reset();
 
     //! get tolerance from Belos solver
-    double get_tolerance() const
+    [[nodiscard]] double get_tolerance() const
     {
-      return params().sublist("Belos Parameters").get<double>("Convergence Tolerance", 1.0e-8);
+      return params().sublist("Belos Parameters").get<double>("Convergence Tolerance");
     }
 
     /*!
@@ -164,7 +164,7 @@ namespace Core::LinAlg
 
     \sa adapt_tolerance
     */
-    void reset_tolerance();
+    void reset_tolerance() const;
 
     //@}
     //! @name Input of parameters
@@ -270,10 +270,10 @@ namespace Core::LinAlg
 
     */
     void adapt_tolerance(
-        const double desirednlnres, const double currentnlnres, const double better);
+        const double desirednlnres, const double currentnlnres, const double better) const;
 
     //! set tolerance to Belos solver
-    void set_tolerance(double tolerance);
+    void set_tolerance(double tolerance) const;
 
     //! a communicator
     MPI_Comm comm_;

src/core/linear_solver/src/method/4C_linear_solver_method_parameters.cpp

@@ -45,18 +45,26 @@ void Core::LinearSolver::Parameters::compute_solver_parameters(
     else
     {
       // if a map is given, grab the block information of the first element in that map
+      bool exit_loops(false);
       for (int i = 0; i < dis.num_my_row_nodes(); ++i)
       {
-        Core::Nodes::Node* actnode = dis.l_row_node(i);
-        std::vector<int> dofs = dis.dof(0, actnode);
-
-        const int localIndex = nullspace_node_map->lid(dofs[0]);
-
-        if (localIndex == -1) continue;
-
-        Core::Elements::Element* dwele = dis.l_row_element(localIndex);
-        actnode->elements()[0]->element_type().nodal_block_information(dwele, numdf, dimns, nv, np);
-        break;
+        const auto* current_node = dis.l_row_node(i);
+        const auto node_gid = current_node->id();
+        for (int j = 0; j < dis.num_my_row_elements(); ++j)
+        {
+          auto* element = dis.l_row_element(j);
+
+          if (const auto* const element_node_gids = element->node_ids();
+              std::find(element_node_gids, element_node_gids + element->num_node(), node_gid) !=
+              element_node_gids + element->num_node())
+          {
+            current_node->elements()[0]->element_type().nodal_block_information(
+                element, numdf, dimns, nv, np);
+            exit_loops = true;
+            break;
+          }
+        }
+        if (exit_loops) break;
       }
     }
 
@@ -67,6 +75,18 @@ void Core::LinearSolver::Parameters::compute_solver_parameters(
     numdf = gdata[0];
     dimns = gdata[1];
 
+    // for dof split within one field the number of global elements in the dof map and the node map
+    // have to be equal
+    if (nullspace_node_map != nullptr and nullspace_dof_map != nullptr)
+    {
+      // TODO Discuss with Max if this should work
+      if (nullspace_dof_map->num_global_elements() == nullspace_node_map->num_global_elements())
+      {
+        numdf = 1;
+        dimns = 1;
+      }
+    }
+
     // store nullspace information in solver list
     solverlist.set("PDE equations", numdf);
     solverlist.set("null space: dimension", dimns);

src/core/linear_solver/src/preconditioner/4C_linear_solver_preconditioner_teko.cpp

@@ -123,6 +123,16 @@ void Core::LinearSolver::TekoPreconditioner::setup(Epetra_Operator* matrix,
                       .get<std::shared_ptr<Core::LinAlg::MultiVector<double>>>("Coordinates")
                       ->get_epetra_multi_vector()));
 
+          // TODO Do temporary printing here
+          /*Teuchos::RCP<Teuchos::FancyOStream> out =
+              Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
+          std::cout << "Teko: " << inverse << '\n';
+          std::cout << "Nullspace: size: " << nullspace->getGlobalLength() << '\n';
+          nullspace->describe(*out, Teuchos::VERB_EXTREME);
+          std::cout << "Coordinates: size: " << coordinates->getGlobalLength() << '\n';
+          coordinates->describe(*out, Teuchos::VERB_EXTREME);
+          std::cout << "Number of equations: " << number_of_equations << '\n' << "\n";*/
+
           tekoParams.sublist("Inverse Factory Library")
               .sublist(inverse)
               .set("number of equations", number_of_equations);