41 cmake build with mumps (#47)

CMakeLists.txt

@@ -3,6 +3,7 @@ cmake_minimum_required(VERSION 3.16.3)
 project(GMGPolar VERSION 1.0.0)
 
 option(GMGPOLAR_BUILD_TESTS "Build GMGPolar unit tests." ON)
+option(GMGPOLAR_USE_MUMPS "Use MUMPS to compute matrix factorizations." OFF)
 
 set(CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/cmake" ${CMAKE_MODULE_PATH})
 set(CMAKE_POSITION_INDEPENDENT_CODE ON)
@@ -36,9 +37,48 @@ add_library(GMGPolar ${SOURCES_SRC})
 
 add_executable(gmgpolar_simulation ./src/main.cpp)
 
+configure_file(${CMAKE_SOURCE_DIR}/include/config_internal.h.in ${CMAKE_SOURCE_DIR}/include/config_internal.h)
+
 target_include_directories(gmgpolar_simulation PRIVATE ${CMAKE_SOURCE_DIR}/include ${CMAKE_SOURCE_DIR}/include/test_cases )
 target_include_directories(GMGPolar PRIVATE ${CMAKE_SOURCE_DIR}/include ${CMAKE_SOURCE_DIR}/include/test_cases )
 
+if(GMGPOLAR_USE_MUMPS)
+
+    set(INC_DIRS
+        /home/kueh_mj/.spack/rev.23.05/install/linux-rocky8-zen2/gcc-10.4.0/mumps-5.4.1-fftqkl/include
+        /sw/rev/23.05/linux-rocky8-zen2/gcc-10.4.0/metis-5.1.0-akhgsf/include
+    )
+
+    set(LIB_DIRS
+        /home/kueh_mj/.spack/rev.23.05/install/linux-rocky8-zen2/gcc-10.4.0/mumps-5.4.1-fftqkl/lib
+        /sw/rev/23.05/linux-rocky8-zen2/gcc-10.4.0/metis-5.1.0-akhgsf/lib
+    )
+
+    include_directories(
+        ${INC_DIRS}
+    )
+
+    target_link_directories(
+        GMGPolar
+        PUBLIC
+        ${LIB_DIRS}
+    )
+
+    set(LIBS
+        mpiseq
+        dmumps
+        mumps_common
+        metis
+    )
+
+    target_link_libraries(
+        GMGPolar
+        PUBLIC
+        ${LIBS}
+    ) 
+endif()
+
+
 find_package(OpenMP)
 
 #currently works on GNU compiler
@@ -58,8 +98,10 @@ else()
     message(FATAL_ERROR "Please use GNU compiler or change CMakeLists manually")
 endif()
 
+
 target_link_libraries(gmgpolar_simulation PRIVATE GMGPolar)
 
+
 include(thirdparty/CMakeLists.txt)
 add_subdirectory(tests)
 

README.md

@@ -21,8 +21,10 @@ Obtaining the source code
 The GmgPolar Solver does not require any external libraries.
 It is possible to link the code with the sparse direct solver ``MUMPS``.
 
-* ``MUMPS`` is optionnal. To use it, compile the code with option -DUSE_MUMPS. It is 
-  recommended to use the latest version (currently 5.4.1) but any version ulterior 
+* ``MUMPS`` is optional. However, it is absolutely recommended if large grids are considered.
+  Otherwise, the nonoptimal backup solver will be used for factorization of the matrices and will
+  slow down the setup phase significantly. To use it, compile the code with option -DGMGPOLAR_USE_MUMPS. 
+  It is recommended to use the latest version (currently 5.4.1) but any version ulterior 
   to 5.1.0 should be okay. MUMPS is available freely on demand on the MUMPS consortium 
   website "mumps-solver.org".
 

include/config.h

@@ -0,0 +1,32 @@
+/* 
+* Copyright (C) 2019-2023 The GMGPolar Development Team
+*
+* Authors: Martin J. Kühn
+*
+* Contact: 
+*    Carola Kruse <Carola.Kruse@CERFACS.fr>
+*    Martin J. Kuehn <Martin.Kuehn@DLR.de>
+*
+* Licensed under the Apache License, Version 2.0 (the "License");
+* you may not use this file except in compliance with the License.
+* You may obtain a copy of the License at
+*
+*     http://www.apache.org/licenses/LICENSE-2.0
+*
+* Unless required by applicable law or agreed to in writing, software
+* distributed under the License is distributed on an "AS IS" BASIS,
+* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+* See the License for the specific language governing permissions and
+* limitations under the License.
+*/
+
+/**
+ * Configuration of GMGPolar library.
+ */
+
+#ifndef GMGPOLAR_CONFIG_H
+#define GMGPOLAR_CONFIG_H
+
+#include "config_internal.h"
+
+#endif // GMGPOLAR_CONFIG_H

include/config_internal.h

@@ -0,0 +1,32 @@
+/* 
+* Copyright (C) 2019-2023 The GMGPolar Development Team
+*
+* Authors: Martin J. Kühn
+*
+* Contact: 
+*    Carola Kruse <Carola.Kruse@CERFACS.fr>
+*    Martin J. Kuehn <Martin.Kuehn@DLR.de>
+*
+* Licensed under the Apache License, Version 2.0 (the "License");
+* you may not use this file except in compliance with the License.
+* You may obtain a copy of the License at
+*
+*     http://www.apache.org/licenses/LICENSE-2.0
+*
+* Unless required by applicable law or agreed to in writing, software
+* distributed under the License is distributed on an "AS IS" BASIS,
+* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+* See the License for the specific language governing permissions and
+* limitations under the License.
+*/
+
+/**
+ * Configured by cmake.
+ */
+
+#ifndef GMGPOLAR_CONFIG_INTERNAL_H
+#define GMGPOLAR_CONFIG_INTERNAL_H
+
+#define GMGPOLAR_USE_MUMPS
+
+#endif // GMGPOLAR_CONFIG_INTERNAL_H

include/config_internal.h.in

@@ -0,0 +1,32 @@
+/* 
+* Copyright (C) 2019-2023 The GMGPolar Development Team
+*
+* Authors: Martin J. Kühn
+*
+* Contact: 
+*    Carola Kruse <Carola.Kruse@CERFACS.fr>
+*    Martin J. Kuehn <Martin.Kuehn@DLR.de>
+*
+* Licensed under the Apache License, Version 2.0 (the "License");
+* you may not use this file except in compliance with the License.
+* You may obtain a copy of the License at
+*
+*     http://www.apache.org/licenses/LICENSE-2.0
+*
+* Unless required by applicable law or agreed to in writing, software
+* distributed under the License is distributed on an "AS IS" BASIS,
+* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+* See the License for the specific language governing permissions and
+* limitations under the License.
+*/
+
+/**
+ * Configured by cmake.
+ */
+
+#ifndef GMGPOLAR_CONFIG_INTERNAL_H
+#define GMGPOLAR_CONFIG_INTERNAL_H
+
+#cmakedefine GMGPOLAR_USE_MUMPS
+
+#endif // GMGPOLAR_CONFIG_INTERNAL_H

include/gmgpolar.h

@@ -41,14 +41,15 @@
 #include <sys/stat.h>
 #include <iomanip>
 #include <ctime>
+#include "config.h"
 #include "constants.h"
 #include "level.h"
 #include "gyro.h"
 
 class gmgpolar
 {
 public:
-    /*******************************************************************************
+/*******************************************************************************
  * Attributes
  ******************************************************************************/
     /***************************************************************************

include/level.h

@@ -44,9 +44,11 @@
 #include <fstream>
 #include <iomanip>
 #include <sys/stat.h>
+
+#include "config.h"
 #include "gyro.h"
 
-#ifdef USE_MUMPS
+#ifdef GMGPOLAR_USE_MUMPS
     #include "mpi.h"
     #include "dmumps_c.h"
 
@@ -91,12 +93,11 @@ class level
     // - using in-house solver
     std::vector<int> row_Ac_LU, col_Ac_LU;
     std::vector<double> vals_Ac_LU;
-#ifdef USE_MUMPS
+#ifdef GMGPOLAR_USE_MUMPS
     // - using MUMPS
     DMUMPS_STRUC_C mumps_Ac;
     DMUMPS_STRUC_C mumps_across;
 #endif
-
     /*! Beta coefficient update */
     std::vector<double> betaVec;
 
@@ -154,7 +155,7 @@ class level
     std::vector<std::vector<int>> A_Zebra_r_LU_row[4];
     std::vector<std::vector<int>> A_Zebra_c_LU_row[4];
     std::vector<std::vector<double>> A_Zebra_v_LU_row[4];
-#ifdef USE_MUMPS
+#ifdef GMGPOLAR_USE_MUMPS
     // - using MUMPS
     DMUMPS_STRUC_C mumps_A_Zebra[4];
 
@@ -197,7 +198,6 @@ class level
     /***************************************************************************
      * gmgpolar
      **************************************************************************/
-    void build_bound();
     void define_coarse_nodes_onelevel(level* finer);
     void store_theta_n_co();
     // void define_colors();
@@ -287,7 +287,7 @@ class level
                                     std::vector<double>& A_vals, int m_solution);
     std::vector<double> solve_gaussian_elimination(std::vector<int> A_row_indices, std::vector<int> A_col_indices,
                                                    std::vector<double> A_vals, std::vector<double> f);
-#ifdef USE_MUMPS
+#ifdef GMGPOLAR_USE_MUMPS
     void init_mumps(DMUMPS_STRUC_C& mumps);
     void facto_mumps(DMUMPS_STRUC_C& mumps, std::vector<int> A_row_indices, std::vector<int> A_col_indices,
                      std::vector<double> A_vals, int m_solution);

src/RHS.cpp

@@ -46,8 +46,9 @@ double gyro::eval_def_rhs(double r, double theta, int verbose)
     else {
         rhs_val = functions->rho_pole(r, theta, kappa_eps, delta_e, Rmax);
     }
-    if (verbose)
+    if (verbose) {
         std::cout << "RHS(" << r << ", " << theta << "): " << rhs_val << "\n";
+    }
     return rhs_val;
 } /* ----- end of gyro::eval_def_rhs ----- */
 
@@ -71,8 +72,9 @@ std::vector<double> gyro::eval_def_rhs(double r, std::vector<double> theta, std:
     else {
         functions->rho_pole(r, theta, kappa_eps, delta_e, Rmax, rhs_val, sin_theta, cos_theta);
     }
-    if (verbose)
+    if (verbose) {
         for (int i = 0; i < ntheta; i++)
             std::cout << "RHS(" << r << ", " << theta[i] << "): " << rhs_val[i] << "\n";
+    }
     return rhs_val;
 } /* ----- end of gyro::eval_def_rhs ----- */

src/direct_solve.cpp

@@ -28,7 +28,7 @@
  */
 #include "level.h"
 
-#ifdef USE_MUMPS
+#ifdef GMGPOLAR_USE_MUMPS
 /*!
  *  \brief Initialization of a MUMPS structure
  *

src/level.cpp

@@ -43,7 +43,7 @@ level::level(int l_)
     reset_timers();
     delete_circles = 0;
 
-#ifdef USE_MUMPS
+#ifdef GMGPOLAR_USE_MUMPS
     init_mumps(mumps_Ac);
     if (gyro::icntl[Param::optimized] == 0) {
         for (int i = 0; i < 4; i++) {
@@ -64,7 +64,7 @@ level::level(int l_)
  */
 level::~level()
 {
-#ifdef USE_MUMPS
+#ifdef GMGPOLAR_USE_MUMPS
     finalize_mumps(mumps_Ac);
     if (gyro::icntl[Param::optimized] == 0) {
         for (int i = 0; i < 4; i++) {
@@ -124,36 +124,6 @@ void level::display_theta()
     gyro::disp(theta, "Coordinates theta");
 } /* ----- end of level::display_theta ----- */
 
-/*!
- *  \brief Defines which nodes are on the boundary
- *
- *  returns ndistance to Dirichlet boundary in binary (0: o boundary, >0: not on boundary)
- *  uses tolerance of 1e-10 to define boundary
- *
- */
-void level::build_bound()
-{
-    double tol_bound_check = gyro::dcntl[Param::tol_bound_check];
-    double r0_DB           = gyro::dcntl[Param::r0_DB];
-    double R               = gyro::dcntl[Param::R];
-
-    is_bound = std::vector<int>(m);
-
-    for (int j = 0; j < nr; j++) {
-        int is_DB = fabs((r[j] - r0_DB) * (r[j] - R)) < tol_bound_check;
-        for (int i = 0; i < ntheta_int; i++) {
-            is_bound[j * ntheta_int + i] = is_DB;
-        }
-    }
-
-    if (gyro::icntl[Param::verbose] > 5)
-        std::cout << "Printing for (r, theta) if the node is on the boundary.\n";
-        for (int j = 0; j < nr; j++)
-            for (int i = 0; i < ntheta_int; i++)
-                std::cout << "(" << r[j] << ", " << theta[j] << "): " << is_bound[j * ntheta_int + i]
-                          << "\n";
-} /* ----- end of gyro::build_bound ----- */
-
 /*!
  *  \brief Defines the number of entries in A
  *

src/main.cpp

@@ -461,5 +461,6 @@ int main(int argc, char* argv[])
             }
         }
     }
+
     return error;
 } /* ----- end of main ----- */

src/multigrid_iter.cpp

@@ -477,12 +477,12 @@ void gmgpolar::multigrid_cycle_extrapol(int l)
     if (l == levels - 2) {
 //solve exactly on the coarsest level for the error (A * error = res) (use whole A from coarsest level)
 // check for the second coarsest level (levels-2), because we have no smoothing on the coarsest level (and thus no Asc and no Asc_ortho)
-#ifdef USE_MUMPS
+#ifdef GMGPOLAR_USE_MUMPS
         if (gyro::icntl[Param::optimized] == 0) {
 #endif
             error_coarse = v_level[l + 1]->solve_gaussian_elimination(
                 v_level[l + 1]->row_Ac_LU, v_level[l + 1]->col_Ac_LU, v_level[l + 1]->vals_Ac_LU, v_level[l + 1]->fVec);
-#ifdef USE_MUMPS
+#ifdef GMGPOLAR_USE_MUMPS
         }
         else
             error_coarse = v_level[l + 1]->solve_mumps(v_level[l + 1]->mumps_Ac, v_level[l + 1]->fVec);