[NOX] Inherit matrix-free operators from SparseOperator

src/fsi/src/monolithic/nonlinear_solver/4C_fsi_nox_linearsystem_gcr.cpp

@@ -14,7 +14,6 @@
 #include "4C_utils_shared_ptr_from_ref.hpp"
 
 #include <Epetra_LinearProblem.h>
-#include <Epetra_Operator.h>
 #include <Teuchos_ParameterList.hpp>
 
 #include <vector>
@@ -26,8 +25,8 @@ NOX::FSI::LinearSystemGCR::LinearSystemGCR(Teuchos::ParameterList& printParams,
     Teuchos::ParameterList& linearSolverParams,
     const std::shared_ptr<NOX::Nln::Interface::RequiredBase> iReq,
     const std::shared_ptr<NOX::Nln::Interface::JacobianBase> iJac,
-    const std::shared_ptr<Epetra_Operator>& jacobian, const NOX::Nln::Vector& cloneVector,
-    const Teuchos::RCP<::NOX::Epetra::Scaling> s)
+    const std::shared_ptr<Core::LinAlg::SparseOperator>& jacobian,
+    const NOX::Nln::Vector& cloneVector, const Teuchos::RCP<::NOX::Epetra::Scaling> s)
     : utils(printParams),
       jacInterfacePtr(iJac),
       jacType(EpetraOperator),
@@ -63,19 +62,15 @@ void NOX::FSI::LinearSystemGCR::reset(Teuchos::ParameterList& linearSolverParams
 bool NOX::FSI::LinearSystemGCR::apply_jacobian(
     const NOX::Nln::Vector& input, NOX::Nln::Vector& result) const
 {
-  jacPtr->SetUseTranspose(false);
-  int status = jacPtr->Apply(input.get_linalg_vector(), result.get_linalg_vector());
+  int status = jacPtr->multiply(false, input.get_linalg_vector(), result.get_linalg_vector());
   return status == 0;
 }
 
 
 bool NOX::FSI::LinearSystemGCR::apply_jacobian_transpose(
     const NOX::Nln::Vector& input, NOX::Nln::Vector& result) const
 {
-  jacPtr->SetUseTranspose(true);
-  int status = jacPtr->Apply(input.get_linalg_vector(), result.get_linalg_vector());
-  jacPtr->SetUseTranspose(false);
-
+  int status = jacPtr->multiply(true, input.get_linalg_vector(), result.get_linalg_vector());
   return status == 0;
 }
 

src/fsi/src/monolithic/nonlinear_solver/4C_fsi_nox_linearsystem_gcr.hpp

@@ -10,6 +10,7 @@
 
 #include "4C_config.hpp"
 
+#include "4C_linalg_sparseoperator.hpp"
 #include "4C_solver_nonlin_nox_interface_jacobian_base.hpp"
 #include "4C_solver_nonlin_nox_interface_required_base.hpp"
 #include "4C_solver_nonlin_nox_linearsystem_base.hpp"
@@ -57,7 +58,8 @@ namespace NOX
           Teuchos::ParameterList& linearSolverParams,
           const std::shared_ptr<NOX::Nln::Interface::RequiredBase> iReq,
           const std::shared_ptr<NOX::Nln::Interface::JacobianBase> iJac,
-          const std::shared_ptr<Epetra_Operator>& J, const NOX::Nln::Vector& cloneVector,
+          const std::shared_ptr<Core::LinAlg::SparseOperator>& J,
+          const NOX::Nln::Vector& cloneVector,
           const Teuchos::RCP<::NOX::Epetra::Scaling> scalingObject = Teuchos::null);
 
       //! Reset the linear solver parameters.
@@ -152,7 +154,7 @@ namespace NOX
       OperatorType jacType;
 
       //! Pointer to the Jacobian operator.
-      mutable std::shared_ptr<Epetra_Operator> jacPtr;
+      std::shared_ptr<Core::LinAlg::SparseOperator> jacPtr;
 
       //! Scaling object supplied by the user
       Teuchos::RCP<::NOX::Epetra::Scaling> scaling;

src/fsi/src/partitioned/4C_fsi_partitioned.cpp

@@ -503,7 +503,7 @@ Teuchos::RCP<NOX::Nln::LinearSystemBase> FSI::Partitioned::create_linear_system(
 
   std::shared_ptr<NOX::Nln::Interface::JacobianBase> iJac;
 
-  std::shared_ptr<Epetra_Operator> J;
+  std::shared_ptr<Core::LinAlg::SparseOperator> J;
 
   Teuchos::RCP<NOX::Nln::LinearSystemBase> linSys;
 
@@ -549,7 +549,7 @@ Teuchos::RCP<NOX::Nln::LinearSystemBase> FSI::Partitioned::create_linear_system(
     auto MF = std::make_shared<NOX::Nln::MatrixFree>(
         printParams, epetra_rcp_interface, noxSoln, lambda, kelleyPerturbation);
     iJac = MF;
-    J = Core::Utils::shared_ptr_from_ref(MF->get_matrix_free());
+    J = Core::Utils::shared_ptr_from_ref(MF->get_operator());
   }
 
   // No Jacobian at all. Do a fix point iteration.

src/fsi/src/partitioned/nonlinear_solver/4C_fsi_nox_jacobian.cpp

@@ -52,23 +52,93 @@ NOX::FSI::FSIMatrixFree::FSIMatrixFree(Teuchos::ParameterList& printParams,
   }
 }
 
+Epetra_Operator& NOX::FSI::FSIMatrixFree::epetra_operator() { return *this; }
 
+void NOX::FSI::FSIMatrixFree::zero() { FOUR_C_THROW("Not implemented"); }
 
-int NOX::FSI::FSIMatrixFree::SetUseTranspose(bool UseTranspose)
+void NOX::FSI::FSIMatrixFree::reset() { FOUR_C_THROW("Not implemented"); }
+
+void NOX::FSI::FSIMatrixFree::assemble(int eid, const std::vector<int>& lmstride,
+    const Core::LinAlg::SerialDenseMatrix& Aele, const std::vector<int>& lmrow,
+    const std::vector<int>& lmrowowner, const std::vector<int>& lmcol)
+{
+  FOUR_C_THROW("Not implemented");
+}
+
+void NOX::FSI::FSIMatrixFree::assemble(double val, int rgid, int cgid)
+{
+  FOUR_C_THROW("Not implemented");
+}
+
+bool NOX::FSI::FSIMatrixFree::filled() const { FOUR_C_THROW("Not implemented"); }
+
+void NOX::FSI::FSIMatrixFree::complete(Core::LinAlg::OptionsMatrixComplete options_matrix_complete)
+{
+  FOUR_C_THROW("Not implemented");
+}
+
+void NOX::FSI::FSIMatrixFree::complete(const Core::LinAlg::Map& domainmap,
+    const Core::LinAlg::Map& rangemap, Core::LinAlg::OptionsMatrixComplete options_matrix_complete)
+{
+  FOUR_C_THROW("Not implemented");
+}
+
+void NOX::FSI::FSIMatrixFree::un_complete() { FOUR_C_THROW("Not implemented"); }
+
+void NOX::FSI::FSIMatrixFree::apply_dirichlet(
+    const Core::LinAlg::Vector<double>& dbctoggle, bool diagonalblock)
+{
+  FOUR_C_THROW("Not implemented");
+}
+
+void NOX::FSI::FSIMatrixFree::apply_dirichlet(const Core::LinAlg::Map& dbcmap, bool diagonalblock)
+{
+  FOUR_C_THROW("Not implemented");
+}
+
+bool NOX::FSI::FSIMatrixFree::is_dbc_applied(const Core::LinAlg::Map& dbcmap, bool diagonalblock,
+    const Core::LinAlg::SparseMatrix* trafo) const
+{
+  FOUR_C_THROW("Not implemented");
+}
+
+const Core::LinAlg::Map& NOX::FSI::FSIMatrixFree::domain_map() const
+{
+  FOUR_C_THROW("Not implemented");
+}
+
+void NOX::FSI::FSIMatrixFree::add(const Core::LinAlg::SparseOperator& A, const bool transposeA,
+    const double scalarA, const double scalarB)
+{
+  FOUR_C_THROW("Not implemented");
+}
+
+void NOX::FSI::FSIMatrixFree::add_other(Core::LinAlg::SparseMatrix& A, const bool transposeA,
+    const double scalarA, const double scalarB) const
 {
-  if (UseTranspose == true)
+  FOUR_C_THROW("Not implemented");
+}
+
+void NOX::FSI::FSIMatrixFree::add_other(Core::LinAlg::BlockSparseMatrixBase& A,
+    const bool transposeA, const double scalarA, const double scalarB) const
+{
+  FOUR_C_THROW("Not implemented");
+}
+
+int NOX::FSI::FSIMatrixFree::scale(double ScalarConstant) { FOUR_C_THROW("Not implemented"); }
+
+int NOX::FSI::FSIMatrixFree::multiply(bool TransA, const Core::LinAlg::MultiVector<double>& X,
+    Core::LinAlg::MultiVector<double>& Y) const
+{
+  if (TransA == true)
   {
     utils.out()
-        << "ERROR: FSIMatrixFree::SetUseTranspose() - Transpose is unavailable in Matrix-Free mode!"
+        << "ERROR: FSIMatrixFree::multiply() - Transpose is unavailable in Matrix-Free mode!"
         << std::endl;
     throw "NOX Error";
+    return -1;
   }
-  return (-1);
-}
 
-
-int NOX::FSI::FSIMatrixFree::Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
-{
   // Calculate the matrix-vector product:
   //
   // y = R' x = S'(F(d)) F'(d) x - x
@@ -82,15 +152,12 @@ int NOX::FSI::FSIMatrixFree::Apply(const Epetra_MultiVector& X, Epetra_MultiVect
   // Convert X and Y from an Epetra_MultiVector to a Core::LinAlg::Vectors
   // and NOX::Nln::Vectors.  This is done so we use a consistent
   // vector space for norms and inner products.
-  Core::LinAlg::View wrappedX(X);
-  Core::LinAlg::View wrappedY(Y);
 
   // There is a const_cast introduced - should be removed
-  NOX::Nln::Vector nevX(Core::Utils::shared_ptr_from_ref(
-                            const_cast<Core::LinAlg::Vector<double>&>(wrappedX.underlying()(0))),
-      NOX::Nln::Vector::MemoryType::View);
-  NOX::Nln::Vector nevY(Core::Utils::shared_ptr_from_ref(wrappedY.underlying()(0)),
+  NOX::Nln::Vector nevX(
+      Core::Utils::shared_ptr_from_ref(const_cast<Core::LinAlg::Vector<double>&>(X(0))),
       NOX::Nln::Vector::MemoryType::View);
+  NOX::Nln::Vector nevY(Core::Utils::shared_ptr_from_ref(Y(0)), NOX::Nln::Vector::MemoryType::View);
 
   // The trial vector x is not guaranteed to be a suitable interface
   // displacement. It might be much too large to fit the ALE
@@ -132,43 +199,74 @@ int NOX::FSI::FSIMatrixFree::Apply(const Epetra_MultiVector& X, Epetra_MultiVect
   return 0;
 }
 
+int NOX::FSI::FSIMatrixFree::SetUseTranspose(bool UseTranspose)
+{
+  FOUR_C_THROW("Not implemented");
+  return -1;
+}
+
+
+int NOX::FSI::FSIMatrixFree::Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
+{
+  FOUR_C_THROW("Not implemented");
+  return -1;
+}
+
 
 int NOX::FSI::FSIMatrixFree::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
 {
-  utils.out() << "ERROR: FSIMatrixFree::ApplyInverse - Not available for Matrix Free!" << std::endl;
-  throw "NOX Error";
-  return (-1);
+  FOUR_C_THROW("Not implemented");
+  return -1;
 }
 
 
 double NOX::FSI::FSIMatrixFree::NormInf() const
 {
-  utils.out() << "ERROR: FSIMatrixFree::NormInf() - Not Available for Matrix-Free mode!"
-              << std::endl;
-  throw "NOX Error";
+  FOUR_C_THROW("Not implemented");
   return 1.0;
 }
 
 
-const char* NOX::FSI::FSIMatrixFree::Label() const { return label.c_str(); }
+const char* NOX::FSI::FSIMatrixFree::Label() const
+{
+  FOUR_C_THROW("Not implemented");
+  return label.c_str();
+}
 
 
-bool NOX::FSI::FSIMatrixFree::UseTranspose() const { return false; }
+bool NOX::FSI::FSIMatrixFree::UseTranspose() const
+{
+  FOUR_C_THROW("Not implemented");
+  return false;
+}
 
 
-bool NOX::FSI::FSIMatrixFree::HasNormInf() const { return false; }
+bool NOX::FSI::FSIMatrixFree::HasNormInf() const
+{
+  FOUR_C_THROW("Not implemented");
+  return false;
+}
 
 
 const Epetra_Comm& NOX::FSI::FSIMatrixFree::Comm() const
 {
+  FOUR_C_THROW("Not implemented");
   return Core::Communication::as_epetra_comm(currentX.get_linalg_vector().get_map().get_comm());
 }
 
 
-const Epetra_Map& NOX::FSI::FSIMatrixFree::OperatorDomainMap() const { return *epetraMap; }
+const Epetra_Map& NOX::FSI::FSIMatrixFree::OperatorDomainMap() const
+{
+  FOUR_C_THROW("Not implemented");
+  return *epetraMap;
+}
 
 
-const Epetra_Map& NOX::FSI::FSIMatrixFree::OperatorRangeMap() const { return *epetraMap; }
+const Epetra_Map& NOX::FSI::FSIMatrixFree::OperatorRangeMap() const
+{
+  FOUR_C_THROW("Not implemented");
+  return *epetraMap;
+}
 
 
 bool NOX::FSI::FSIMatrixFree::computeJacobian(const Epetra_Vector& x, Epetra_Operator& Jac)

src/fsi/src/partitioned/nonlinear_solver/4C_fsi_nox_jacobian.hpp

@@ -10,6 +10,7 @@
 
 #include "4C_config.hpp"
 
+#include "4C_linalg_sparseoperator.hpp"
 #include "4C_solver_nonlin_nox_interface_jacobian_base.hpp"
 #include "4C_solver_nonlin_nox_interface_required_base.hpp"
 #include "4C_solver_nonlin_nox_vector.hpp"
@@ -35,7 +36,8 @@ namespace NOX
   namespace FSI
   {
     /// Matrix Free Newton Krylov based on an approximation of the residuum derivatives
-    class FSIMatrixFree : public Epetra_Operator, public virtual NOX::Nln::Interface::JacobianBase
+    class FSIMatrixFree : public Core::LinAlg::SparseOperator,
+                          public virtual NOX::Nln::Interface::JacobianBase
     {
      public:
       /*! \brief Constructor
@@ -45,7 +47,54 @@ namespace NOX
       FSIMatrixFree(Teuchos::ParameterList& printParams,
           const std::shared_ptr<NOX::Nln::Interface::RequiredBase> i, const NOX::Nln::Vector& x);
 
+      // Methods of Core::LinAlg::SparseOperator interface
+      Epetra_Operator& epetra_operator() override;
 
+      void zero() override;
+
+      void reset() override;
+
+      void assemble(int eid, const std::vector<int>& lmstride,
+          const Core::LinAlg::SerialDenseMatrix& Aele, const std::vector<int>& lmrow,
+          const std::vector<int>& lmrowowner, const std::vector<int>& lmcol) override;
+
+      void assemble(double val, int rgid, int cgid) override;
+
+      bool filled() const override;
+
+      void complete(Core::LinAlg::OptionsMatrixComplete options_matrix_complete = {}) override;
+
+      void complete(const Core::LinAlg::Map& domainmap, const Core::LinAlg::Map& rangemap,
+          Core::LinAlg::OptionsMatrixComplete options_matrix_complete = {}) override;
+
+      void un_complete() override;
+
+      void apply_dirichlet(
+          const Core::LinAlg::Vector<double>& dbctoggle, bool diagonalblock = true) override;
+
+      void apply_dirichlet(const Core::LinAlg::Map& dbcmap, bool diagonalblock = true) override;
+
+      bool is_dbc_applied(const Core::LinAlg::Map& dbcmap, bool diagonalblock = true,
+          const Core::LinAlg::SparseMatrix* trafo = nullptr) const override;
+
+      const Core::LinAlg::Map& domain_map() const override;
+
+      void add(const Core::LinAlg::SparseOperator& A, const bool transposeA, const double scalarA,
+          const double scalarB) override;
+
+      void add_other(Core::LinAlg::SparseMatrix& A, const bool transposeA, const double scalarA,
+          const double scalarB) const override;
+
+      void add_other(Core::LinAlg::BlockSparseMatrixBase& A, const bool transposeA,
+          const double scalarA, const double scalarB) const override;
+
+      int scale(double ScalarConstant) override;
+
+      int multiply(bool TransA, const Core::LinAlg::MultiVector<double>& X,
+          Core::LinAlg::MultiVector<double>& Y) const override;
+
+
+      // Methods of Epetra_Operator interface
       //! If set true, transpose of this operator will be applied.
       /*! This flag allows the transpose of the given operator to be used implicitly.  Setting this
         flag affects only the Apply() and ApplyInverse() methods.  If the implementation of this