Add interface linear solvers

CMakeLists.txt

@@ -83,7 +83,6 @@ include(SimVascularFunctionCheckCompilerFlags)
 # These variables must later be add to 'ExternalProject_Add(svFSI' as -D options. 
 #
 set(SV_USE_TRILINOS OFF CACHE BOOL "Build with the Trilinos linear algebra package")
-#set(SV_USE_PETSC OFF CACHE BOOL "Build with the PETSc linear algebra package")
 set(SV_PETSC_DIR "" CACHE STRING "Path to a local install of the PETSc linear algebra package")
 set(ENABLE_COVERAGE OFF CACHE BOOL "Enable code coverage")
 set(ENABLE_ARRAY_INDEX_CHECKING OFF CACHE BOOL "Enable Array index checking")
@@ -144,7 +143,6 @@ ExternalProject_Add(svFSI
     -DSV_EXTERNALS_USE_TOPLEVEL_DIR:BOOL=ON
     -DSV_EXTERNALS_TOPLEVEL_DIR:PATH=${SV_EXTERNALS_TOPLEVEL_DIR}
     -DSV_USE_TRILINOS:BOOL=${SV_USE_TRILINOS}
-    #-DSV_USE_PETSC:BOOL=${SV_USE_PETSC}
     -DSV_PETSC_DIR:STRING=${SV_PETSC_DIR}
     -DENABLE_COVERAGE:BOOL=${ENABLE_COVERAGE}
     -DENABLE_ARRAY_INDEX_CHECKING:BOOL=${ENABLE_ARRAY_INDEX_CHECKING}

Code/Source/svFSI/CMakeLists.txt

@@ -1,3 +1,27 @@
+# Copyright (c) 2014-2015 The Regents of the University of California.
+# All Rights Reserved.
+#
+# Permission is hereby granted, free of charge, to any person obtaining
+# a copy of this software and associated documentation files (the
+# "Software"), to deal in the Software without restriction, including
+# without limitation the rights to use, copy, modify, merge, publish,
+# distribute, sublicense, and/or sell copies of the Software, and to
+# permit persons to whom the Software is furnished to do so, subject
+# to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included
+# in all copies or substantial portions of the Software.
+#
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+# TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
+# OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 
 set(CMAKE_CXX_STANDARD 17)
 set(CMAKE_CXX_STANDARD_REQUIRED True)
@@ -84,11 +108,12 @@ if(NOT "${SV_PETSC_DIR}" STREQUAL "")
   message("  PETSC_LIBRARY_DIRS = ${PETSC_LIBRARY_DIRS}")
   message("  PETSC_INCLUDE_DIRS = ${PETSC_INCLUDE_DIRS}")
 
-  # Set PETSc include directory.
+  # Set PETSc include and library directories.
   include_directories(${PETSC_INCLUDE_DIRS})
+  #link_directories(${PETSC_LIBRARY_DIRS})
 
-  # Set C++ directive to use PETSc in svFSIplus code.
-  ADD_DEFINITIONS(-DWITH_PETSC)
+  # Set C++ directive to use PETSc.
+  ADD_DEFINITIONS(-DUSE_PETSC)
   set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -pedantic -std=c99")
 
   set(USE_PETSC 1)
@@ -123,6 +148,10 @@ set(lib ${SV_LIB_SVFSI_NAME})
 set(CSRCS 
   Array3.h Array3.cpp 
   Array.h Array.cpp
+  LinearAlgebra.h LinearAlgebra.cpp
+  FsilsLinearAlgebra.h FsilsLinearAlgebra.cpp
+  PetscLinearAlgebra.h PetscLinearAlgebra.cpp
+  TrilinosLinearAlgebra.h TrilinosLinearAlgebra.cpp
   Tensor4.h Tensor4.cpp
   Vector.h Vector.cpp 
 
@@ -192,12 +221,12 @@ set(CSRCS
 )
 
   # Set PETSc interace code.
-  if(USE_PETSC)
-     #set(CSRCS ${CSRCS} petsc_linear_solver.c)
-  endif()
+  #if(USE_PETSC)
+  #   set(CSRCS ${CSRCS} petsc_linear_solver.h petsc_linear_solver.c)
+  #endif()
 
 if(USE_TRILINOS)
-  set(CSRCS ${CSRCS} trilinos_linear_solver.cpp)
+  #set(CSRCS ${CSRCS} trilinos_linear_solver.cpp)
 
   # trilinos directories and libraries
   include_directories(${Trilinos_TPL_INCLUDE_DIRS})

Code/Source/svFSI/ComMod.h

@@ -56,6 +56,8 @@
 #include <string>
 #include <vector>
 
+class LinearAlgebra;
+
 /// @brief Fourier coefficients that are used to specify unsteady BCs
 //
 class fcType
@@ -613,7 +615,6 @@ class outputType
     std::string name;
 };
 
-
 /// @brief Linear system of equations solver type
 //
 class lsType
@@ -623,9 +624,6 @@ class lsType
     /// @brief LS solver                     (IN)
     consts::SolverType LS_type = consts::SolverType::lSolver_NA;
 
-    /// @brief Preconditioner                (IN)
-    consts::PreconditionerType PREC_Type = consts::PreconditionerType::PREC_NONE;
-
     /// @brief Successful solving            (OUT)
     bool suc = false;
 
@@ -667,6 +665,9 @@ class lsType
 
     /// @brief Calling duration              (OUT)
     double callD = 0.0;
+
+    //@brief Configuration file for linear solvers (Trilinos, PETSc)
+    std::string config;
 };
 
 
@@ -1067,6 +1068,18 @@ class eqType
     /// @brief type of linear solver
     lsType ls;
 
+    /// @brief The type of interface to a numerical linear algebra library.
+    consts::LinearAlgebraType linear_algebra_type;
+
+    /// @brief The type of assembly interface to a numerical linear algebra library.
+    consts::LinearAlgebraType linear_algebra_assembly_type;
+
+    /// @brief The type of preconditioner used by the interface to a numerical linear algebra library.
+    consts::PreconditionerType linear_algebra_preconditioner = consts::PreconditionerType::PREC_FSILS;
+
+    /// @brief Interface to a numerical linear algebra library.
+    LinearAlgebra* linear_algebra = nullptr;
+
     /// @brief FSILS type of linear solver
     fsi_linear_solver::FSILS_lsType FSILS;
 
@@ -1285,23 +1298,6 @@ class ibType
     ibCommType cm;
 };
 
-/// @brief Data type for Trilinos Linear Solver related arrays
-//
-class tlsType
-{
-  public:
-
-    /// @brief Local to global mapping
-    Vector<int> ltg;
-
-    /// @brief Factor for Dirichlet BCs
-    Array<double> W;
-
-    /// @brief Residual
-    Array<double> R;
-};
-
-
 /// @brief The ComMod class duplicates the data structures in the Fortran COMMOD module
 /// defined in MOD.f. 
 ///
@@ -1582,9 +1578,6 @@ class ComMod {
     /// @brief IB: Immersed boundary data structure
     ibType ib;
 
-    /// @brief Trilinos Linear Solver data type
-    tlsType  tls;
-
     bool debug_active = false;
 
     Timer timer;

Code/Source/svFSI/FsilsLinearAlgebra.cpp

@@ -0,0 +1,149 @@
+/* Copyright (c) Stanford University, The Regents of the University of California, and others.
+ *
+ * All Rights Reserved.
+ *
+ * See Copyright-SimVascular.txt for additional details.
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+ * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+ * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+ * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
+ * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+ * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+ * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+ * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+ * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#include "FsilsLinearAlgebra.h"
+#include "fsils_api.hpp"
+#include "lhsa.h"
+#include <iostream>
+
+/////////////////////////////////////////////////////////////////
+//             F s i l s L i n e a r A l g e b r a             //
+/////////////////////////////////////////////////////////////////
+// The following methods implement the FSILS LinearAlgebra interface.
+
+std::set<consts::LinearAlgebraType> FsilsLinearAlgebra::valid_assemblers = {
+  consts::LinearAlgebraType::none,
+  consts::LinearAlgebraType::fsils,
+};
+
+FsilsLinearAlgebra::FsilsLinearAlgebra()
+{
+  interface_type = consts::LinearAlgebraType::fsils; 
+  assembly_type = consts::LinearAlgebraType::fsils; 
+  preconditioner_type = consts::PreconditionerType::PREC_FSILS;
+}
+
+/// @brief Allocate data arrays.
+void FsilsLinearAlgebra::alloc(ComMod& com_mod, eqType& lEq)
+{
+  #define n_debug_alloc
+  #ifdef debug_alloc
+  std::cout << "[FsilsLinearAlgebra::alloc] ---------- alloc ---------- " << std::endl;
+  #endif
+  int dof = com_mod.dof;
+
+  com_mod.Val.resize(dof*dof, com_mod.lhs.nnz);
+}
+
+/// @brief Assemble local element arrays.
+void FsilsLinearAlgebra::assemble(ComMod& com_mod, const int num_elem_nodes, const Vector<int>& eqN,
+        const Array3<double>& lK, const Array<double>& lR)
+{
+  #define n_debug_assemble
+  #ifdef debug_assemble
+  std::cout << "[FsilsLinearAlgebra::assemble] ---------- assemble ---------- " << std::endl;
+  std::cout << "[FsilsLinearAlgebra::assemble] num_elem_nodes: " << num_elem_nodes << std::endl;
+  std::cout << "[FsilsLinearAlgebra::assemble] eqN.size(): " << eqN.size() << std::endl;
+  std::cout << "[FsilsLinearAlgebra::assemble] lK.size(): " << lK.size() << std::endl;
+  std::cout << "[FsilsLinearAlgebra::assemble] lR.size(): " << lR.size() << std::endl;
+  #endif
+
+  lhsa_ns::do_assem(com_mod, num_elem_nodes, eqN, lK, lR);
+}
+
+/// @brief Check the validity of the preconditioner and assembly types options. 
+void FsilsLinearAlgebra::check_options(const consts::PreconditionerType prec_cond_type, 
+  const consts::LinearAlgebraType assembly_type)
+{
+  using namespace consts;
+  auto prec_cond_type_name = consts::preconditioner_type_to_name.at(prec_cond_type);
+  auto assembly_type_name = LinearAlgebra::type_to_name.at(assembly_type);
+  std::string error_msg;
+
+  if (valid_assemblers.count(assembly_type) == 0) {
+    error_msg = "fsils linear algebra can't use '" + assembly_type_name + "' for assembly.";
+  }
+
+  if (fsils_preconditioners.count(prec_cond_type) == 0) { 
+    error_msg = "fsils linear algebra can't use '" + prec_cond_type_name + "' for a preconditioner.";
+  }
+
+  if (error_msg != "") { 
+    throw std::runtime_error("[svFSIplus] ERROR: " + error_msg);
+  }
+}
+
+/// @brief Initialize framework.
+void FsilsLinearAlgebra::initialize(ComMod& com_mod, eqType& lEq)
+{
+  // Nothing is needed to initialize FSILS.
+}
+
+/// @brief Set the linear algebra package for assmbly.
+void FsilsLinearAlgebra::set_assembly(consts::LinearAlgebraType atype)
+{
+  if (atype == consts::LinearAlgebraType::none) {
+    return;
+  }
+
+  if (valid_assemblers.count(atype) == 0) {
+    auto str_type = LinearAlgebra::type_to_name.at(atype);
+    throw std::runtime_error("[FsilsLinearAlgebra] ERROR: Can't set fsils linear algebra to use '" +
+      str_type + "' for assembly.");
+  }
+
+  assembly_type = atype;
+}
+
+/// @brief Set the preconditioner.
+void FsilsLinearAlgebra::set_preconditioner(consts::PreconditionerType prec_type)
+{
+  if (consts::fsils_preconditioners.count(prec_type) == 0) {
+    auto prec_cond_type_name = consts::preconditioner_type_to_name.at(prec_type);
+    throw std::runtime_error("[FsilsLinearAlgebra] ERROR: fsils linear algebra can't use '" + 
+        prec_cond_type_name + "' for a preconditioner.");
+    return;
+  }
+
+  preconditioner_type = prec_type;
+}
+
+/// @brief Solve a system of linear equations.
+void FsilsLinearAlgebra::solve(ComMod& com_mod, eqType& lEq, const Vector<int>& incL, const Vector<double>& res)
+{
+  auto& lhs = com_mod.lhs;
+  int dof = com_mod.dof;
+  auto& R = com_mod.R;      
+  auto& Val = com_mod.Val;
+  auto preconditioner = lEq.linear_algebra_preconditioner;
+
+  fsi_linear_solver::fsils_solve(lhs, lEq.FSILS, dof, R, Val, preconditioner, incL, res);
+}
+

Code/Source/svFSI/FsilsLinearAlgebra.h

@@ -0,0 +1,60 @@
+/* Copyright (c) Stanford University, The Regents of the University of California, and others.
+ *
+ * All Rights Reserved.
+ *
+ * See Copyright-SimVascular.txt for additional details.
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+ * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+ * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+ * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
+ * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+ * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+ * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+ * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+ * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#ifndef FSILS_LINEAR_ALGEBRA_H 
+#define FSILS_LINEAR_ALGEBRA_H 
+
+#include "LinearAlgebra.h"
+
+/// @brief The FsilsLinearAlgebra class implements the LinearAlgebra 
+/// interface for the FSILS numerical linear algebra included in svFSIplus.
+///
+class FsilsLinearAlgebra : public virtual LinearAlgebra {
+
+  public:
+    FsilsLinearAlgebra();
+    ~FsilsLinearAlgebra() { };
+
+    virtual void alloc(ComMod& com_mod, eqType& lEq);
+    virtual void assemble(ComMod& com_mod, const int num_elem_nodes, const Vector<int>& eqN,
+        const Array3<double>& lK, const Array<double>& lR);
+    virtual void check_options(const consts::PreconditionerType prec_cond_type, const consts::LinearAlgebraType assembly_type);
+    virtual void initialize(ComMod& com_mod, eqType& lEq);
+    virtual void solve(ComMod& com_mod, eqType& lEq, const Vector<int>& incL, const Vector<double>& res);
+    virtual void set_assembly(consts::LinearAlgebraType atype);
+    virtual void set_preconditioner(consts::PreconditionerType prec_type);
+
+  private:
+    /// @brief A list of linear algebra interfaces that can be used for assembly.
+    static std::set<consts::LinearAlgebraType> valid_assemblers; 
+};
+
+#endif
+