From ad33dde2acd7a64fc2c173acf35f366560bf98dc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Martin=20J=2E=20K=C3=BChn?= <62713180+mknaranja@users.noreply.github.com> Date: Tue, 3 Oct 2023 15:35:29 +0200 Subject: [PATCH] 41 cmake build with mumps (#47) --- CMakeLists.txt | 42 ++++++++++++++++++++++++++++++++++++ README.md | 6 ++++-- include/config.h | 32 +++++++++++++++++++++++++++ include/config_internal.h | 32 +++++++++++++++++++++++++++ include/config_internal.h.in | 32 +++++++++++++++++++++++++++ include/gmgpolar.h | 3 ++- include/level.h | 12 +++++------ src/RHS.cpp | 6 ++++-- src/direct_solve.cpp | 2 +- src/level.cpp | 34 ++--------------------------- src/main.cpp | 1 + src/multigrid_iter.cpp | 4 ++-- src/polar_multigrid.cpp | 21 +++++++++--------- src/smoother.cpp | 8 +++---- src/smoother0.cpp | 4 ++-- 15 files changed, 176 insertions(+), 63 deletions(-) create mode 100644 include/config.h create mode 100644 include/config_internal.h create mode 100644 include/config_internal.h.in diff --git a/CMakeLists.txt b/CMakeLists.txt index dbaa177c..13143f65 100644 --- a/CMakeLists.txt +++ b/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) diff --git a/README.md b/README.md index 5f9091eb..8d057932 100644 --- a/README.md +++ b/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". diff --git a/include/config.h b/include/config.h new file mode 100644 index 00000000..ea43a7d3 --- /dev/null +++ b/include/config.h @@ -0,0 +1,32 @@ +/* +* Copyright (C) 2019-2023 The GMGPolar Development Team +* +* Authors: Martin J. Kühn +* +* Contact: +* Carola Kruse +* Martin J. Kuehn +* +* 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 diff --git a/include/config_internal.h b/include/config_internal.h new file mode 100644 index 00000000..38e6a59d --- /dev/null +++ b/include/config_internal.h @@ -0,0 +1,32 @@ +/* +* Copyright (C) 2019-2023 The GMGPolar Development Team +* +* Authors: Martin J. Kühn +* +* Contact: +* Carola Kruse +* Martin J. Kuehn +* +* 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 diff --git a/include/config_internal.h.in b/include/config_internal.h.in new file mode 100644 index 00000000..de594a47 --- /dev/null +++ b/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 +* Martin J. Kuehn +* +* 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 diff --git a/include/gmgpolar.h b/include/gmgpolar.h index f164e4e1..9de56530 100644 --- a/include/gmgpolar.h +++ b/include/gmgpolar.h @@ -41,6 +41,7 @@ #include #include #include +#include "config.h" #include "constants.h" #include "level.h" #include "gyro.h" @@ -48,7 +49,7 @@ class gmgpolar { public: - /******************************************************************************* +/******************************************************************************* * Attributes ******************************************************************************/ /*************************************************************************** diff --git a/include/level.h b/include/level.h index 31e90775..f8167002 100644 --- a/include/level.h +++ b/include/level.h @@ -44,9 +44,11 @@ #include #include #include + +#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 row_Ac_LU, col_Ac_LU; std::vector 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 betaVec; @@ -154,7 +155,7 @@ class level std::vector> A_Zebra_r_LU_row[4]; std::vector> A_Zebra_c_LU_row[4]; std::vector> 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& A_vals, int m_solution); std::vector solve_gaussian_elimination(std::vector A_row_indices, std::vector A_col_indices, std::vector A_vals, std::vector f); -#ifdef USE_MUMPS +#ifdef GMGPOLAR_USE_MUMPS void init_mumps(DMUMPS_STRUC_C& mumps); void facto_mumps(DMUMPS_STRUC_C& mumps, std::vector A_row_indices, std::vector A_col_indices, std::vector A_vals, int m_solution); diff --git a/src/RHS.cpp b/src/RHS.cpp index 285cf6ca..109d82fd 100644 --- a/src/RHS.cpp +++ b/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 gyro::eval_def_rhs(double r, std::vector 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 ----- */ diff --git a/src/direct_solve.cpp b/src/direct_solve.cpp index 160a38bb..4b021225 100644 --- a/src/direct_solve.cpp +++ b/src/direct_solve.cpp @@ -28,7 +28,7 @@ */ #include "level.h" -#ifdef USE_MUMPS +#ifdef GMGPOLAR_USE_MUMPS /*! * \brief Initialization of a MUMPS structure * diff --git a/src/level.cpp b/src/level.cpp index aa86bef4..42f08a2e 100644 --- a/src/level.cpp +++ b/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(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 * diff --git a/src/main.cpp b/src/main.cpp index 20a3d0fc..c64100b5 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -461,5 +461,6 @@ int main(int argc, char* argv[]) } } } + return error; } /* ----- end of main ----- */ diff --git a/src/multigrid_iter.cpp b/src/multigrid_iter.cpp index fd889afd..b33b39a5 100644 --- a/src/multigrid_iter.cpp +++ b/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); diff --git a/src/polar_multigrid.cpp b/src/polar_multigrid.cpp index d2070004..5cd6ed05 100644 --- a/src/polar_multigrid.cpp +++ b/src/polar_multigrid.cpp @@ -143,8 +143,7 @@ void gmgpolar::polar_multigrid() cycle_str = "V"; else if (gyro::icntl[Param::cycle] == 2) cycle_str = "W"; - if (gyro::icntl[Param::verbose] > 2) - { + if (gyro::icntl[Param::verbose] > 2) { std::cout << "\nProb: " << gyro::icntl[Param::prob] << ", alpha_coeff: " << gyro::icntl[Param::alpha_coeff] << ", beta_coeff: " << gyro::icntl[Param::beta_coeff] << " ***** Problem size " << m << " (" << v_level[0]->nr << ", " << v_level[0]->ntheta << "), " << levels @@ -265,9 +264,6 @@ void gmgpolar::prepare_op_levels() if (l < levels - 1) v_level[l]->mc = v_level[l + 1]->m; - if (gyro::icntl[Param::verbose] > 3) - std::cout << "Create boundary array\n"; - v_level[l]->build_bound(); if (l == 0 && !gyro::f_sol_in.empty()) { v_level[0]->read_sol(); @@ -300,8 +296,11 @@ void gmgpolar::prepare_op_levels() if (gyro::icntl[Param::verbose] > 3) std::cout << "Factorizing coarse operator...\n"; TIC; -#ifdef USE_MUMPS +#ifdef GMGPOLAR_USE_MUMPS + std::cout << "\n Using GMGPolar with MUMPS\n"; if (gyro::icntl[Param::optimized] == 0) { +#else + std::cout << "\n Attention: Using GMGPolar without MUMPS (Coarse solve is very slow)\n"; #endif v_level[l]->row_Ac_LU = std::vector(v_level[l]->row_indices); v_level[l]->col_Ac_LU = std::vector(v_level[l]->col_indices); @@ -314,7 +313,7 @@ void gmgpolar::prepare_op_levels() v_level[l]->vals.shrink_to_fit(); v_level[l]->facto_gaussian_elimination(v_level[l]->row_Ac_LU, v_level[l]->col_Ac_LU, v_level[l]->vals_Ac_LU, v_level[l]->m); -#ifdef USE_MUMPS +#ifdef GMGPOLAR_USE_MUMPS } else v_level[l]->facto_mumps(v_level[l]->mumps_Ac, v_level[l]->row_indices, v_level[l]->col_indices, @@ -437,7 +436,7 @@ void gmgpolar::prepare_op_levels() TIC; if (gyro::icntl[Param::optimized] == 0) { -#ifdef USE_MUMPS +#ifdef GMGPOLAR_USE_MUMPS if (gyro::icntl[Param::optimized] == 0) { #endif v_level[l]->A_Zebra_r_LU.assign(4, std::vector()); @@ -465,7 +464,7 @@ void gmgpolar::prepare_op_levels() v_level[l]->A_Zebra_r.shrink_to_fit(); v_level[l]->A_Zebra_c.shrink_to_fit(); v_level[l]->A_Zebra_v.shrink_to_fit(); -#ifdef USE_MUMPS +#ifdef GMGPOLAR_USE_MUMPS } else for (int smoother = 0; smoother < 4; smoother++) @@ -544,7 +543,7 @@ void gmgpolar::prepare_op_levels() std::vector(v_level[l]->A_Zebra_v_row[smoother][ij]); if (gyro::icntl[Param::verbose] > 3) std::cout << "Across factorization..."; -#ifdef USE_MUMPS +#ifdef GMGPOLAR_USE_MUMPS if (gyro::icntl[Param::optimized] == 0) { #endif if (gyro::icntl[Param::verbose] > 3) @@ -554,7 +553,7 @@ void gmgpolar::prepare_op_levels() v_level[l]->A_Zebra_c_LU_row[smoother][ij], v_level[l]->A_Zebra_v_LU_row[smoother][ij], v_level[l]->m_sc[smoother]); -#ifdef USE_MUMPS +#ifdef GMGPOLAR_USE_MUMPS } else { if (gyro::icntl[Param::verbose] > 2) diff --git a/src/smoother.cpp b/src/smoother.cpp index 126c9940..0c183373 100644 --- a/src/smoother.cpp +++ b/src/smoother.cpp @@ -247,12 +247,12 @@ void level::multigrid_smoothing(int smoother, int v, std::vector& f_Asc_ } // Across (direct solver) else { -#ifdef USE_MUMPS +#ifdef GMGPOLAR_USE_MUMPS if (gyro::icntl[Param::optimized] == 0) { #endif u_sc = solve_gaussian_elimination(A_Zebra_r_LU_row[smoother][k], A_Zebra_c_LU_row[smoother][k], A_Zebra_v_LU_row[smoother][k], f_sc); -#ifdef USE_MUMPS +#ifdef GMGPOLAR_USE_MUMPS } else u_sc = solve_mumps(mumps_across, f_sc); @@ -322,12 +322,12 @@ void level::multigrid_smoothing(int smoother, int v, std::vector& f_Asc_ } // Across (direct solver) else { -#ifdef USE_MUMPS +#ifdef GMGPOLAR_USE_MUMPS if (gyro::icntl[Param::optimized] == 0) { #endif u_sc = solve_gaussian_elimination(A_Zebra_r_LU_row[smoother][k], A_Zebra_c_LU_row[smoother][k], A_Zebra_v_LU_row[smoother][k], f_sc); -#ifdef USE_MUMPS +#ifdef GMGPOLAR_USE_MUMPS } else u_sc = solve_mumps(mumps_across, f_sc); diff --git a/src/smoother0.cpp b/src/smoother0.cpp index e6c21258..7d86c6fb 100644 --- a/src/smoother0.cpp +++ b/src/smoother0.cpp @@ -70,12 +70,12 @@ void level::multigrid_smoothing0(int smoother) t_Asc_ortho += TOC; TIC; -#ifdef USE_MUMPS +#ifdef GMGPOLAR_USE_MUMPS if (gyro::icntl[Param::optimized] == 0) { #endif u_sc = solve_gaussian_elimination(A_Zebra_r_LU[smoother], A_Zebra_c_LU[smoother], A_Zebra_v_LU[smoother], f_total); -#ifdef USE_MUMPS +#ifdef GMGPOLAR_USE_MUMPS } else u_sc = solve_mumps(mumps_A_Zebra[smoother], f_total);