From 6255add479b95c5ada0dc7002345b486f3eed633 Mon Sep 17 00:00:00 2001 From: maorz1998 Date: Sat, 8 Apr 2023 18:52:02 +0800 Subject: [PATCH 1/3] use amgx to solve dfMatrix --- .../solvers/dfLowMachFoam/CMakeLists.txt | 10 + .../solvers/dfLowMachFoam/Make/options | 7 +- applications/solvers/dfLowMachFoam/UEqn.H | 28 +- .../solvers/dfLowMachFoam/createdfSolver.H | 3 +- src_gpu/AmgXSolver.H | 309 ++++++++++++++++++ src_gpu/AmgXSolver.cu | 296 +++++++++++++++++ src_gpu/CMakeLists.txt | 15 +- src_gpu/dfMatrix.H | 19 +- src_gpu/dfMatrix.cu | 209 +++++++++++- 9 files changed, 867 insertions(+), 29 deletions(-) create mode 100644 src_gpu/AmgXSolver.H create mode 100644 src_gpu/AmgXSolver.cu diff --git a/applications/solvers/dfLowMachFoam/CMakeLists.txt b/applications/solvers/dfLowMachFoam/CMakeLists.txt index 3c489e081..645289a64 100644 --- a/applications/solvers/dfLowMachFoam/CMakeLists.txt +++ b/applications/solvers/dfLowMachFoam/CMakeLists.txt @@ -1,6 +1,8 @@ cmake_minimum_required(VERSION 3.5) project(dfLowMachFoam LANGUAGES CXX) FIND_PACKAGE(MPI REQUIRED) +FIND_PACKAGE(OpenMP REQUIRED) +FIND_PACKAGE(CUDA REQUIRED) # Check valid thirdParty if(DEFINED ENV{WM_PROJECT_DIR}) @@ -26,6 +28,8 @@ SET(SRC_ORIG $ENV{SRC_ORIG}) # set compilation options SET(CMAKE_EXE_LINKER_FLAGS "-fuse-ld=bfd -Xlinker --add-needed -Xlinker --no-as-needed") +SET (CMAKE_C_FLAGS ${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}) +SET (CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}) SET(CMAKE_C_COMPILER g++) SET(PATH_LIB_OPENMPI "openmpi-system") # Foundation version @@ -83,6 +87,9 @@ include_directories( ${CANTERA_ROOT}/include ${MPI_INCLUDE_PATH} ${PROJECT_SOURCE_DIR} + ${CUDA_INCLUDE_DIRS} + /home/runze/AmgX/AMGX/include + /home/runze/deepflame-dev/src_gpu ) # add execution @@ -98,6 +105,9 @@ target_link_libraries(${PROJECT_NAME} ${DF_ROOT}/lib/libdfCombustionModels.so $ENV{FOAM_LIBBIN}/openmpi-system/libPstream.so ${MPI_LIBRARIES} + ${CUDA_LIBRARIES} + /home/runze/AmgX/AMGX/build/libamgxsh.so + /home/runze/deepflame-dev/src_gpu/build/libdfMatrix.so ) if(DEFINED ENV{PYTHON_INC_DIR}) diff --git a/applications/solvers/dfLowMachFoam/Make/options b/applications/solvers/dfLowMachFoam/Make/options index 49190c6b2..43ea7958c 100644 --- a/applications/solvers/dfLowMachFoam/Make/options +++ b/applications/solvers/dfLowMachFoam/Make/options @@ -1,6 +1,7 @@ -include $(GENERAL_RULES)/mplibType EXE_INC = -std=c++14 \ + -fopenmp \ -Wno-unused-variable \ -Wno-unused-but-set-variable \ -Wno-old-style-cast \ @@ -25,7 +26,8 @@ EXE_INC = -std=c++14 \ $(if $(LIBTORCH_ROOT),-I$(LIBTORCH_ROOT)/include/torch/csrc/api/include,) \ $(PYTHON_INC_DIR) \ -I/home/runze/deepflame-dev/src_gpu \ - -I/usr/local/cuda-11.6/include + -I/usr/local/cuda-11.6/include \ + -I/home/runze/AmgX/AMGX/include EXE_LIBS = \ -lcompressibleTransportModels \ @@ -48,5 +50,6 @@ EXE_LIBS = \ $(if $(PYTHON_LIB_DIR),-L$(PYTHON_LIB_DIR),) \ $(if $(PYTHON_LIB_DIR),-lpython3.8,) \ /home/runze/deepflame-dev/src_gpu/build/libdfMatrix.so \ - /usr/local/cuda-11.6/lib64/libcudart.so + /usr/local/cuda-11.6/lib64/libcudart.so \ + /home/runze/AmgX/AMGX/build/libamgxsh.so diff --git a/applications/solvers/dfLowMachFoam/UEqn.H b/applications/solvers/dfLowMachFoam/UEqn.H index 983c1aaa0..6d9cb4c0c 100644 --- a/applications/solvers/dfLowMachFoam/UEqn.H +++ b/applications/solvers/dfLowMachFoam/UEqn.H @@ -7,14 +7,20 @@ tmp tUEqn ); fvVectorMatrix& UEqn = tUEqn.ref(); -UEqn.relax(); +fvVectorMatrix turb_source +( + turbulence->divDevRhoReff(U) +); -if (pimple.momentumPredictor()) -{ - solve(UEqn == -fvc::grad(p)); - K = 0.5*magSqr(U); -} +// UEqn.relax(); + +// if (pimple.momentumPredictor()) +// { +// solve(UEqn == -fvc::grad(p)); + +// K = 0.5*magSqr(U); +// } int offset = 0; @@ -44,13 +50,13 @@ forAll(U.boundaryField(), patchi) offset += patchSize; } -fvVectorMatrix turb_source -( - turbulence->divDevRhoReff(U) -); UEqn_GPU.fvm_ddt(&rho.oldTime()[0], &rho[0], &mesh.V()[0], &U.oldTime()[0][0]); UEqn_GPU.fvm_div(&mesh.surfaceInterpolation::weights()[0], &phi[0], ueqn_internalCoeffs_init, ueqn_boundaryCoeffs_init); UEqn_GPU.fvc_grad(&mesh.Sf()[0][0], &p[0], boundary_face_vector_init, boundary_pressure_init); UEqn_GPU.add_fvMatrix(&turb_source.lower()[0], &turb_source.diag()[0], &turb_source.upper()[0], &turb_source.source()[0][0]); -UEqn_GPU.print(); \ No newline at end of file +// UEqn_GPU.print(); +UEqn_GPU.solve(); +UEqn_GPU.updatePsi(&U[0][0]); +K = 0.5*magSqr(U); + diff --git a/applications/solvers/dfLowMachFoam/createdfSolver.H b/applications/solvers/dfLowMachFoam/createdfSolver.H index f76b6e953..7b9bc7696 100644 --- a/applications/solvers/dfLowMachFoam/createdfSolver.H +++ b/applications/solvers/dfLowMachFoam/createdfSolver.H @@ -15,7 +15,8 @@ forAll(mesh.boundary(), patchi) } int num_boundary_cells; -dfMatrix UEqn_GPU(num_surfaces, num_cells, num_boundary_faces, num_boundary_cells, &neighbour[0], &owner[0], boundaryCellIndex); +dfMatrix UEqn_GPU(num_surfaces, num_cells, num_boundary_faces, num_boundary_cells, &neighbour[0], &owner[0], boundaryCellIndex, "dDDI", + "/home/runze/deepflame-dev/examples/dfLowMachFoam/threeD_reactingTGV/H2/cvodeIntegrator/system/amgxUOptions"); std::vector ueqn_internalCoeffs_init(3*num_boundary_faces), ueqn_boundaryCoeffs_init(3*num_boundary_faces); std::vector boundary_pressure_init(num_boundary_faces), boundary_face_vector_init(3*num_boundary_faces); \ No newline at end of file diff --git a/src_gpu/AmgXSolver.H b/src_gpu/AmgXSolver.H new file mode 100644 index 000000000..28baf006f --- /dev/null +++ b/src_gpu/AmgXSolver.H @@ -0,0 +1,309 @@ +/** + * \file AmgXSolver.hpp + * \brief Definition of class AmgXSolver. + * \author Pi-Yueh Chuang (pychuang@gwu.edu) + * \author Matt Martineau (mmartineau@nvidia.com) + * \date 2015-09-01 + * \copyright Copyright (c) 2015-2019 Pi-Yueh Chuang, Lorena A. Barba. + * \copyright Copyright (c) 2019-2021, NVIDIA CORPORATION. All rights reserved. + * This project is released under MIT License. + */ + + +#ifndef __AMGX_SOLVER_H__ +#define __AMGX_SOLVER_H__ + +// CUDA +#include + +// STL +# include +# include + +// AmgX +# include + +// PETSc +// # include + + +/** \brief A macro to check the returned CUDA error code. + * + * \param call [in] Function call to CUDA API. + */ +# define CHECK(call) \ +do \ +{ \ + const cudaError_t error_code = call; \ + if (error_code != cudaSuccess) \ + { \ + printf("CUDA Error:\n"); \ + printf(" File: %s\n", __FILE__); \ + printf(" Line: %d\n", __LINE__); \ + printf(" Error code: %d\n", error_code); \ + printf(" Error text: %s\n", \ + cudaGetErrorString(error_code)); \ + exit(1); \ + } \ +} while (0) + + + + + + +/** \brief A wrapper class for coupling PETSc and AmgX. + * + * This class is a wrapper of AmgX library for PETSc. PETSc users only need to + * pass a PETSc matrix and vectors into an AmgXSolver instance to solve their + * linear systems. The class is designed specifically for the situation where + * the number of MPI processes is more than the number of GPU devices. + * + * Eaxmple usage: + * \code + * int main(int argc, char **argv) + * { + * // initialize matrix A, RHS, etc using PETSc + * ... + * + * // create an instance of the solver wrapper + * AmgXSolver solver; + * // initialize the instance with communicator, executation mode, and config file + * solver.initialize(comm, mode, file); + * // set matrix A. Currently it only accept PETSc AIJ matrix + * solver.setA(A); + * // solve. x and rhs are PETSc vectors. unkns will be the final result in the end + * solver.solve(unks, rhs); + * // get number of iterations + * int iters; + * solver.getIters(iters); + * // get residual at the last iteration + * double res; + * solver.getResidual(iters, res); + * // finalization + * solver.finalize(); + * + * // other codes + * .... + * + * return 0; + * } + * \endcode + */ +class AmgXSolver +{ + public: + + /** \brief Default constructor. */ + AmgXSolver() = default; + + /** \brief Construct a AmgXSolver instance. + * + * \param comm [in] MPI communicator. + * \param modeStr [in] A string; target mode of AmgX (e.g., dDDI). + * \param cfgFile [in] A string; the path to AmgX configuration file. + */ + AmgXSolver + ( + const std::string &modeStr, + const std::string &cfgFile + ); + + /** \brief Destructor. */ + ~AmgXSolver(); + + /** \brief Initialize a AmgXSolver instance. + * + * \param comm [in] MPI communicator. + * \param modeStr [in] A string; target mode of AmgX (e.g., dDDI). + * \param cfgFile [in] A string; the path to AmgX configuration file. + * + */ + void initialize + ( + const std::string &modeStr, + const std::string &cfgFile + ); + + + /** \brief Finalize this instance. + * + * This function destroys AmgX data. When there are more than one + * AmgXSolver instances, the last one destroyed is also in charge of + * destroying the shared resource object and finalizing AmgX. + * + */ + void finalize(); + + /** \brief Set up the matrix used by AmgX. + * + * This function sets up the AmgX matrix from the provided CSR data + * structures and partition data. + * + * \param nGlobalRows [in] The number of global rows. + * \param nLocalRows [in] The number of local rows on this rank. + * \param nLocalNz [in] The total number of non zero entries locally. + * \param rowOffsets [in] The local CSR matrix row offsets. + * \param colIndicesGlobal [in] The global CSR matrix column indices. + * \param values [in] The local CSR matrix values. + * id of the owning rank for each row. + * + */ + void setOperator + ( + const int nRows, + const int nNz, + const int *rowIndex, + const int *colIndex, + const double *value + ); + + /** \brief Re-sets up an existing AmgX matrix. + * + * Replaces the matrix coefficients with the provided values and performs + * a resetup for the AmgX matrix. + * + * \param nLocalRows [in] The number of local rows on this rank. + * \param nLocalNz [in] The total number of non zero entries locally. + * \param values [in] The local CSR matrix values. + * + */ + void updateOperator + ( + const int nRows, + const int nNz, + const double *value + ); + + /** \brief Solve the linear system. + * + * \p p vector will be used as an initial guess and will be updated to the + * solution by the end of solving. + * + * For cases that use more MPI processes than the number of GPUs, this + * function will do data gathering before solving and data scattering + * after the solving. + * + * \param nLocalRows [in] The number of rows owned by this rank. + * \param pscalar [in, out] The unknown array. + * \param bscalar [in] The RHS array. + * \param matrix [in,out] The AmgX CSR matrix, A. + * + */ + void solve + ( + int nRows, + double* psi, + const double* rhs + ); + + /** \brief Solve the linear system. + * + * \p p vector will be used as an initial guess and will be updated to the + * solution by the end of solving. + * + * For cases that use more MPI processes than the number of GPUs, this + * function will do data gathering before solving and data scattering + * after the solving. + * + * \param nLocalRows [in] The number of rows owned by this rank. + * \param p [in, out] The unknown vector. + * \param b [in] The RHS vector. + * \param matrix [in,out] The AmgX CSR matrix, A. + * + */ + // void solve + // ( + // int nLocalRows, + // Vec& p, + // Vec& b, + // AmgXCSRMatrix& matrix + // ); + + + /** \brief Get the number of iterations of the last solving. + * + * \param iter [out] Number of iterations. + * + */ + void getIters + ( + int &iter + ); + + /** \brief Get the residual at a specific iteration during the last solving. + * + * \param iter [in] Target iteration. + * \param res [out] Returned residual. + * + */ + void getResidual + ( + const int &iter, + double &res + ); + + + private: + + /** \brief Current count of AmgXSolver instances. + * + * This static variable is used to count the number of instances. The + * fisrt instance is responsable for initializing AmgX library and the + * resource instance. + */ + static int count; + + /** \brief A flag indicating if this instance has been initialized. */ + bool isInitialised = false; + + /** \brief A parameter used by AmgX. */ + int ring; + + /** \brief AmgX solver mode. */ + AMGX_Mode mode; + + /** \brief AmgX config object. */ + AMGX_config_handle cfg = nullptr; + + /** \brief AmgX matrix object. */ + AMGX_matrix_handle AmgXA = nullptr; + + /** \brief AmgX vector object representing unknowns. */ + AMGX_vector_handle AmgXP = nullptr; + + /** \brief AmgX vector object representing RHS. */ + AMGX_vector_handle AmgXRHS = nullptr; + + /** \brief AmgX solver object. */ + AMGX_solver_handle solver = nullptr; + + /** \brief AmgX resource object. + * + * Due to the design of AmgX library, using more than one resource + * instance may cause some problems. So we make the resource instance + * as a static member to keep only one instance. + */ + static AMGX_resources_handle rsrc; + + /** \brief Set AmgX solver mode based on the user-provided string. + * + * Available modes are: dDDI, dDFI, dFFI, hDDI, hDFI, hFFI. + * + * \param modeStr [in] a std::string. + */ + void setMode(const std::string &modeStr); + + /** \brief Perform necessary initialization of AmgX. + * + * This function initializes AmgX for current instance. Based on + * \ref AmgXSolver::count "count", only the instance initialized first + * is in charge of initializing AmgX and the resource instance. + * + * \param cfgFile [in] Path to AmgX solver configuration file. + */ + void initAmgX(const std::string &cfgFile); +}; + +#endif + diff --git a/src_gpu/AmgXSolver.cu b/src_gpu/AmgXSolver.cu new file mode 100644 index 000000000..b0076e5c3 --- /dev/null +++ b/src_gpu/AmgXSolver.cu @@ -0,0 +1,296 @@ +/** + * \file AmgXSolver.cpp + * \brief Definition of member functions of the class AmgXSolver. + * \author Pi-Yueh Chuang (pychuang@gwu.edu) + * \author Matt Martineau (mmartineau@nvidia.com) + * \date 2015-09-01 + * \copyright Copyright (c) 2015-2019 Pi-Yueh Chuang, Lorena A. Barba. + * \copyright Copyright (c) 2019-2021, NVIDIA CORPORATION. All rights reserved. + * This project is released under MIT License. + */ + +// AmgXWrapper +#include "AmgXSolver.H" +#include +#include + +// initialize AmgXSolver::count to 0 +int AmgXSolver::count = 0; + +// initialize AmgXSolver::rsrc to nullptr; +AMGX_resources_handle AmgXSolver::rsrc = nullptr; + + +/* \implements AmgXSolver::AmgXSolver */ +AmgXSolver::AmgXSolver(const std::string &modeStr, const std::string &cfgFile) +{ + initialize(modeStr, cfgFile); +} + + +/* \implements AmgXSolver::~AmgXSolver */ +AmgXSolver::~AmgXSolver() +{ + if (isInitialised) finalize(); +} + + +/* \implements AmgXSolver::initialize */ +void AmgXSolver::initialize(const std::string &modeStr, const std::string &cfgFile) +{ + + // if this instance has already been initialized, skip + if (isInitialised) { + fprintf(stderr, + "This AmgXSolver instance has been initialized on this process.\n"); + exit(0); + } + + // increase the number of AmgXSolver instances + count += 1; + + // get the mode of AmgX solver + setMode(modeStr); + + // initialize AmgX + initAmgX(cfgFile); + + // a bool indicating if this instance is initialized + isInitialised = true; + + return; +} + +/* \implements AmgXSolver::setMode */ +void AmgXSolver::setMode(const std::string &modeStr) +{ + if (modeStr == "dDDI") + mode = AMGX_mode_dDDI; + else if (modeStr == "dDFI") + mode = AMGX_mode_dDFI; + else if (modeStr == "dFFI") + mode = AMGX_mode_dFFI; + else if (modeStr[0] == 'h') { + printf("CPU mode, %s, is not supported in this wrapper!", + modeStr.c_str()); + exit(0); + } + else { + printf("%s is not an available mode! Available modes are: " + "dDDI, dDFI, dFFI.\n", modeStr.c_str()); + exit(0); + } +} + + +/* \implements AmgXSolver::initAmgX */ + void AmgXSolver::initAmgX(const std::string &cfgFile) +{ + // only the first instance (AmgX solver) is in charge of initializing AmgX + if (count == 1) + { + // initialize AmgX + AMGX_SAFE_CALL(AMGX_initialize()); + + // intialize AmgX plugings + AMGX_SAFE_CALL(AMGX_initialize_plugins()); + + // let AmgX to handle errors returned + AMGX_SAFE_CALL(AMGX_install_signal_handler()); + } + + // create an AmgX configure object + AMGX_SAFE_CALL(AMGX_config_create_from_file(&cfg, cfgFile.c_str())); + + // let AmgX handle returned error codes internally + AMGX_SAFE_CALL(AMGX_config_add_parameters(&cfg, "exception_handling=1")); + + // create an AmgX resource object, only the first instance is in charge + if (count == 1) AMGX_resources_create_simple(&rsrc, cfg); + + // create AmgX vector object for unknowns and RHS + AMGX_vector_create(&AmgXP, rsrc, mode); + AMGX_vector_create(&AmgXRHS, rsrc, mode); + + // create AmgX matrix object for unknowns and RHS + AMGX_matrix_create(&AmgXA, rsrc, mode); + + // create an AmgX solver object + AMGX_solver_create(&solver, rsrc, mode, cfg); + + // obtain the default number of rings based on current configuration + AMGX_config_get_default_number_of_rings(cfg, &ring); +} + +/* \implements AmgXSolver::finalize */ +void AmgXSolver::finalize() +{ + // skip if this instance has not been initialised + if (!isInitialised) + { + fprintf(stderr, + "This AmgXWrapper has not been initialised. " + "Please initialise it before finalization.\n"); + exit(0); + } + + // destroy solver instance + AMGX_solver_destroy(solver); + + // destroy matrix instance + AMGX_matrix_destroy(AmgXA); + + // destroy RHS and unknown vectors + AMGX_vector_destroy(AmgXP); + AMGX_vector_destroy(AmgXRHS); + + // only the last instance need to destroy resource and finalizing AmgX + if (count == 1) + { + AMGX_resources_destroy(rsrc); + AMGX_SAFE_CALL(AMGX_config_destroy(cfg)); + + AMGX_SAFE_CALL(AMGX_finalize_plugins()); + AMGX_SAFE_CALL(AMGX_finalize()); + } + else + { + AMGX_config_destroy(cfg); + } + + // decrease the number of instances + count -= 1; + + // change status + isInitialised = false; +} + +/* \implements AmgXSolver::setOperator */ +void AmgXSolver::setOperator +( + const int nRows, + const int nNz, + const int *rowIndex, + const int *colIndex, + const double *value +) +{ + + // Check the matrix size is not larger than tolerated by AmgX + if(nRows > std::numeric_limits::max()) + { + fprintf(stderr, + "AmgX does not support a global number of rows greater than " + "what can be stored in 32 bits (nGlobalRows = %d).\n", + nRows); + exit(0); + } + + if (nNz > std::numeric_limits::max()) + { + fprintf(stderr, + "AmgX does not support non-zeros per (consolidated) rank greater than" + "what can be stored in 32 bits (nLocalNz = %d).\n", + nNz); + exit(0); + } + + // upload matrix A to AmgX + AMGX_matrix_upload_all( + AmgXA, nRows, nNz, 1, 1, rowIndex, colIndex, value, nullptr); + + // bind the matrix A to the solver + AMGX_solver_setup(solver, AmgXA); + + // connect (bind) vectors to the matrix + AMGX_vector_bind(AmgXP, AmgXA); + AMGX_vector_bind(AmgXRHS, AmgXA); +} + + +/* \implements AmgXSolver::updateOperator */ +void AmgXSolver::updateOperator +( + const int nRows, + const int nNz, + const double *value +) +{ + + // Replace the coefficients for the CSR matrix A within AmgX + AMGX_matrix_replace_coefficients(AmgXA, nRows, nNz, value, nullptr); + + // Re-setup the solver (a reduced overhead setup that accounts for consistent matrix structure) + AMGX_solver_resetup(solver, AmgXA); +} + +/* \implements AmgXSolver::solve */ +// void AmgXSolver::solve( +// int nLocalRows, Vec& p, Vec& b, AmgXCSRMatrix& matrix) +// { +// double* pscalar; +// double* bscalar; + +// // get pointers to the raw data of local vectors +// VecGetArray(p, &pscalar); +// VecGetArray(b, &bscalar); + +// solve(nLocalRows, pscalar, bscalar, matrix); + +// VecRestoreArray(p, &pscalar); +// VecRestoreArray(b, &bscalar); +// } + + +/* \implements AmgXSolver::solve */ +void AmgXSolver::solve( + int nRows, double* psi, const double* rhs) +{ + // Upload potentially consolidated vectors to AmgX + AMGX_vector_upload(AmgXP, nRows, 1, psi); + AMGX_vector_upload(AmgXRHS, nRows, 1, rhs); + + // Solve + AMGX_solver_solve(solver, AmgXRHS, AmgXP); + + // Get the status of the solver + AMGX_SOLVE_STATUS status; + AMGX_solver_get_status(solver, &status); + + // Check whether the solver successfully solved the problem + if (status != AMGX_SOLVE_SUCCESS) + { + fprintf(stderr, "AmgX solver failed to solve the system! " + "The error code is %d.\n", + status); + } + + // Download data from device + AMGX_vector_download(AmgXP, psi); + + // get norm and iteration number + double irnorm = 0., rnorm = 0.; + int nIters = 0; + getResidual(0, irnorm); + getIters(nIters); + getResidual(nIters, rnorm); + printf("Initial residual = %.10lf, Final residual = %.5e, No Iterations %d\n", irnorm, rnorm, nIters); + +} + + +/* \implements AmgXSolver::getIters */ +void AmgXSolver::getIters(int &iter) +{ + // only processes using AmgX will try to get # of iterations + AMGX_solver_get_iterations_number(solver, &iter); +} + + +/* \implements AmgXSolver::getResidual */ +void AmgXSolver::getResidual(const int &iter, double &res) +{ + // only processes using AmgX will try to get residual + AMGX_solver_get_iteration_residual(solver, iter, 0, &res); +} + diff --git a/src_gpu/CMakeLists.txt b/src_gpu/CMakeLists.txt index d8e6674fb..2c3207c22 100644 --- a/src_gpu/CMakeLists.txt +++ b/src_gpu/CMakeLists.txt @@ -1,15 +1,26 @@ cmake_minimum_required(VERSION 3.5) project(dfMatrix) FIND_PACKAGE(CUDA REQUIRED) +FIND_PACKAGE(MPI REQUIRED) enable_language(CUDA) +set(AMGX_DIR $ENV{AMGX_DIR}) include_directories( + ${MPI_INCLUDE_PATH} ${CUDA_INCLUDE_DIRS} + ${AMGX_DIR}/include ) SET(CMAKE_C_COMPILER g++) add_compile_options(-std=c++17) -add_library(${PROJECT_NAME} SHARED dfMatrix.cu) +add_library(${PROJECT_NAME} + SHARED + dfMatrix.cu + AmgXSolver.cu) -target_link_libraries(${PROJECT_NAME} ${CUDA_LIBRARIES}) +target_link_libraries(${PROJECT_NAME} + ${MPI_LIBRARIES} + ${CUDA_LIBRARIES} + ${AMGX_DIR}/build/libamgxsh.so +) diff --git a/src_gpu/dfMatrix.H b/src_gpu/dfMatrix.H index ac822693a..140488ae5 100644 --- a/src_gpu/dfMatrix.H +++ b/src_gpu/dfMatrix.H @@ -10,6 +10,8 @@ #include #include +#include "AmgXSolver.H" + static const char *_cudaGetErrorEnum(cudaError_t error) { return cudaGetErrorName(error); } @@ -126,23 +128,32 @@ private: size_t csr_value_vec_bytes; // extra matrix information - double *h_A_csr = nullptr, *h_b = nullptr; - double *d_A_csr = nullptr, *d_b = nullptr; double *d_turbSrc_A = nullptr, *d_turbSrc_b = nullptr; std::vector h_turbSrc_init_mtx_vec, h_turbSrc_init_1mtx; std::vector h_turbSrc_init_src_vec, h_turbSrc_src_vec; std::vector tmpPermutatedList; + double *h_A_csr = nullptr, *h_b = nullptr, *h_psi = nullptr; + double *d_A_csr = nullptr, *d_b = nullptr, *d_psi = nullptr; + /** * @brief cuda functions */ cudaStream_t stream; + /** + * @brief AmgX functions + */ + AmgXSolver* UxSolver = nullptr, *UySolver = nullptr, *UzSolver = nullptr; + int num_iteration; + + public: dfMatrix(); dfMatrix(int num_surfaces, int num_cells, int num_boundary_faces, int &num_boundary_cells_output, - const int *neighbour, const int *owner, std::vector boundary_cell_id_init); + const int *neighbour, const int *owner, std::vector boundary_cell_id_init, + const std::string &modeStr, const std::string &cfgFile); ~dfMatrix(); void print(); @@ -157,7 +168,9 @@ public: void add_fvMatrix(double* turbSrc_low, double* turbSrc_diag, double* turbSrc_upp, double* turbSrc_source); + void solve(const std::string &cfgFile); void solve(); + void updatePsi(double* Psi); }; #endif \ No newline at end of file diff --git a/src_gpu/dfMatrix.cu b/src_gpu/dfMatrix.cu index 9702c063a..e56d858f5 100644 --- a/src_gpu/dfMatrix.cu +++ b/src_gpu/dfMatrix.cu @@ -3,7 +3,7 @@ __global__ void fvm_ddt_kernel(int num_cells, int num_faces, const double rdelta_t, const int* csr_row_index, const int* csr_diag_index, const double* rho_old, const double* rho_new, const double* volume, const double* velocity_old, - const double* A_csr_input, const double* b_input, double* A_csr_output, double* b_output) { + const double* A_csr_input, const double* b_input, double* A_csr_output, double* b_output, double* psi) { int index = blockDim.x * blockIdx.x + threadIdx.x; if (index >= num_cells) return; @@ -22,6 +22,10 @@ __global__ void fvm_ddt_kernel(int num_cells, int num_faces, const double rdelta b_output[num_cells * 0 + index] = b_input[num_cells * 0 + index] + ddt_part_term * velocity_old[index * 3 + 0]; b_output[num_cells * 1 + index] = b_input[num_cells * 1 + index] + ddt_part_term * velocity_old[index * 3 + 1]; b_output[num_cells * 2 + index] = b_input[num_cells * 2 + index] + ddt_part_term * velocity_old[index * 3 + 2]; + + psi[num_cells * 0 + index] = velocity_old[index * 3 + 0]; + psi[num_cells * 1 + index] = velocity_old[index * 3 + 1]; + psi[num_cells * 2 + index] = velocity_old[index * 3 + 2]; } __global__ void fvm_div_internal(int num_cells, int num_faces, @@ -239,11 +243,12 @@ __global__ void add_fvMatrix_kernel(int num_cells, int num_faces, // constructor (construct mesh variable) dfMatrix::dfMatrix(){} dfMatrix::dfMatrix(int num_surfaces, int num_cells, int num_boundary_faces, int & num_boundary_cells_output, - const int *neighbour, const int *owner, std::vector boundary_cell_id_init) -: num_cells(num_cells), num_faces(num_surfaces*2), num_surfaces(num_surfaces), + const int *neighbour, const int *owner, std::vector boundary_cell_id_init, const std::string &modeStr, + const std::string &cfgFile) +: num_cells(num_cells), num_faces(num_surfaces*2), num_surfaces(num_surfaces), num_iteration(0), num_boundary_faces(num_boundary_faces) -{ - // resize vector +{ + // resize vectors h_weight_vec_init.resize(num_faces); h_weight_vec.resize(num_faces); h_phi_vec_init.resize(num_faces); @@ -269,6 +274,10 @@ dfMatrix::dfMatrix(int num_surfaces, int num_cells, int num_boundary_faces, int csr_value_bytes = (num_cells + num_faces) * sizeof(double); csr_value_vec_bytes = (num_cells + num_faces) * 3 * sizeof(double); + h_A_csr = new double[csr_value_vec_bytes]; + h_b = new double[cell_vec_bytes]; + h_psi = new double[cell_vec_bytes]; + /************************construct mesh variables****************************/ /** * 1. h_csr_row_index & h_csr_diag_index @@ -494,6 +503,7 @@ dfMatrix::dfMatrix(int num_surfaces, int num_cells, int num_boundary_faces, int checkCudaErrors(cudaMalloc((void**)&d_A_csr, csr_value_vec_bytes)); checkCudaErrors(cudaMalloc((void**)&d_b, cell_vec_bytes)); + checkCudaErrors(cudaMalloc((void**)&d_psi, cell_vec_bytes)); total_bytes += (boundary_face_bytes + boundary_face_vec_bytes * 3); checkCudaErrors(cudaMalloc((void**)&d_turbSrc_A, csr_value_bytes)); @@ -512,6 +522,11 @@ dfMatrix::dfMatrix(int num_surfaces, int num_cells, int num_boundary_faces, int checkCudaErrors(cudaMemcpyAsync(d_A_csr_diag_index, h_A_csr_diag_index, cell_index_bytes, cudaMemcpyHostToDevice, stream)); checkCudaErrors(cudaMemcpyAsync(d_boundary_cell_offset, h_boundary_cell_offset, (num_boundary_cells+1) * sizeof(int), cudaMemcpyHostToDevice, stream)); checkCudaErrors(cudaMemcpyAsync(d_boundary_cell_id, h_boundary_cell_id, boundary_face_index_bytes, cudaMemcpyHostToDevice, stream)); + + // construct AmgXSolver + UxSolver = new AmgXSolver(modeStr, cfgFile); + UySolver = new AmgXSolver(modeStr, cfgFile); + UzSolver = new AmgXSolver(modeStr, cfgFile); } dfMatrix::~dfMatrix() @@ -539,7 +554,7 @@ void dfMatrix::fvm_ddt(double *rho_old, double *rho_new, const double* volume, printf("CUDA kernel fvm_ddt launch with %d blocks of %d threads\n", blocks_per_grid, threads_per_block); fvm_ddt_kernel<<>>(num_cells, num_faces, rdelta_t, d_A_csr_row_index, d_A_csr_diag_index, - d_rho_old, d_rho_new, d_volume, d_velocity_old, d_A_csr, d_b, d_A_csr, d_b); + d_rho_old, d_rho_new, d_volume, d_velocity_old, d_A_csr, d_b, d_A_csr, d_b, d_psi); // Synchronize stream checkCudaErrors(cudaStreamSynchronize(stream)); } @@ -680,9 +695,6 @@ void dfMatrix::add_fvMatrix(double* turbSrc_low, double* turbSrc_diag, double* t void dfMatrix::print() { - h_A_csr = new double[csr_value_vec_bytes]; - h_b = new double[cell_vec_bytes]; - checkCudaErrors(cudaMemcpyAsync(h_A_csr, d_A_csr, csr_value_vec_bytes, cudaMemcpyDeviceToHost, stream)); checkCudaErrors(cudaMemcpyAsync(h_b, d_b, cell_vec_bytes, cudaMemcpyDeviceToHost, stream)); @@ -694,7 +706,184 @@ void dfMatrix::print() for (int i = 0; i < num_cells * 3; i++) fprintf(stderr, "h_b[%d]: %.15lf\n", i, h_b[i]); + + char *input_file = "of_output.txt"; + FILE *fp = fopen(input_file, "rb+"); + if (fp == NULL) { + fprintf(stderr, "Failed to open input file: %s!\n", input_file); + } + int readfile = 0; + double *of_b = new double[3*num_cells]; + double *of_A = new double[2*num_surfaces + num_cells]; + readfile = fread(of_b, num_cells * 3 * sizeof(double), 1, fp); + readfile = fread(of_A, (2*num_surfaces + num_cells) * sizeof(double), 1, fp); + + std::vector h_A_of_init_vec(num_cells+2*num_surfaces); + std::copy(of_A, of_A + num_cells+2*num_surfaces, h_A_of_init_vec.begin()); + + std::vector h_A_of_vec_1mtx(2*num_surfaces + num_cells, 0); + for (int i = 0; i < 2*num_surfaces + num_cells; i++) + { + h_A_of_vec_1mtx[i] = h_A_of_init_vec[tmpPermutatedList[i]]; + } + for (int i = 0; i < (2*num_surfaces + num_cells); i++) + fprintf(stderr, "h_A_of_vec_1mtx[%d]: %.15lf\n", i, h_A_of_vec_1mtx[i]); + + std::vector h_A_of_vec((2*num_surfaces + num_cells)*3); + for (int i =0; i < 3; i ++) + { + std::copy(h_A_of_vec_1mtx.begin(), h_A_of_vec_1mtx.end(), h_A_of_vec.begin()+i*(2*num_surfaces + num_cells)); + } + + // b + std::vector h_b_of_init_vec(3*num_cells); + std::copy(of_b, of_b + 3*num_cells, h_b_of_init_vec.begin()); + std::vector h_b_of_vec; + for (int i = 0; i < 3*num_cells; i+=3) + { + h_b_of_vec.push_back(h_b_of_init_vec[i]); + } + // fill RHS_y + for (int i = 1; i < 3*num_cells; i+=3) + { + h_b_of_vec.push_back(h_b_of_init_vec[i]); + } + // fill RHS_z + for (int i = 2; i < 3*num_cells; i+=3) + { + h_b_of_vec.push_back(h_b_of_init_vec[i]); + } + for (int i = 0; i < 3*num_cells; i++) + fprintf(stderr, "h_b_of_vec[%d]: %.15lf\n", i, h_b_of_vec[i]); } -void dfMatrix::solve(){} +// void dfMatrix::solve(const std::string &cfgFile) +// { +// // checkCudaErrors(cudaMemcpyAsync(h_A_csr, d_A_csr, csr_value_vec_bytes, cudaMemcpyDeviceToHost, stream)); +// // checkCudaErrors(cudaMemcpyAsync(h_b, d_b, cell_vec_bytes, cudaMemcpyDeviceToHost, stream)); +// checkCudaErrors(cudaMemcpyAsync(h_psi, d_psi, cell_vec_bytes, cudaMemcpyDeviceToHost, stream)); +// /** +// * @brief initialize AmgXSolver +// */ +// // - setMode + +// // - initAmgX(cfgFile) +// // initialize AmgX +// AMGX_SAFE_CALL(AMGX_initialize()); +// // intialize AmgX plugings +// AMGX_SAFE_CALL(AMGX_initialize_plugins()); +// // let AmgX to handle errors returned +// AMGX_SAFE_CALL(AMGX_install_signal_handler()); +// // create an AmgX configure object +// AMGX_config_handle cfg = nullptr; +// AMGX_SAFE_CALL(AMGX_config_create_from_file(&cfg, cfgFile.c_str())); +// // let AmgX handle returned error codes internally +// AMGX_SAFE_CALL(AMGX_config_add_parameters(&cfg, "exception_handling=1")); +// // initialize rsrc to nullptr; +// AMGX_resources_handle rsrc = nullptr; +// // create an AmgX resource object, only the first instance is in charge +// // if (solver_instance_count == 1) +// AMGX_resources_create_simple(&rsrc, cfg); // remember to set deviceID + +// // create AmgX vector object for unknowns and RHS +// AMGX_vector_handle AmgXP = nullptr; +// AMGX_vector_handle AmgXRHS = nullptr; +// AMGX_Mode mode = AMGX_mode_dDDI; +// AMGX_vector_create(&AmgXP, rsrc, mode); +// AMGX_vector_create(&AmgXRHS, rsrc, mode); + +// // create AmgX matrix object for unknowns and RHS +// AMGX_matrix_handle AmgXA = nullptr; +// AMGX_matrix_create(&AmgXA, rsrc, mode); + +// // bind the matrix A to the solver +// AMGX_solver_handle solver = nullptr; +// AMGX_solver_create(&solver, rsrc, mode, cfg); + +// int ring; +// // obtain the default number of rings based on current configuration +// AMGX_config_get_default_number_of_rings(cfg, &ring); + +// AMGX_matrix_upload_all(AmgXA, num_cells, (num_cells + num_faces), 1, 1, d_A_csr_row_index, d_A_csr_col_index, +// d_A_csr, nullptr); + +// AMGX_solver_setup(solver, AmgXA); + +// // bind +// AMGX_vector_bind(AmgXP, AmgXA); +// AMGX_vector_bind(AmgXP, AmgXA); + +// // create an AmgX solver object + +// /** +// * @brief amgx.initialiseMatrixComms(Amat); +// */ + +// /** +// * @brief amgx.setOperator +// */ +// AMGX_vector_upload(AmgXRHS, num_cells, 1, d_b); +// AMGX_vector_upload(AmgXP, num_cells, 1, d_psi); + +// // Solve +// AMGX_solver_solve(solver, AmgXRHS, AmgXP); +// AMGX_SOLVE_STATUS status; +// AMGX_solver_get_status(solver, &status); + +// if (status != AMGX_SOLVE_SUCCESS) +// { +// fprintf(stderr, "AmgX solver failed to solve the system! " +// "The error code is %d.\n", +// status); +// } +// // Download data from device +// double * h_psi_solved; +// h_psi_solved = new double[cell_vec_bytes]; +// AMGX_vector_download(AmgXP, h_psi_solved); +// for (int i = 0; i < num_cells; i++) +// { +// fprintf(stderr, "h_psi_solved[%d] = %lf\n", i, h_psi_solved[i]); +// } + +// } + +void dfMatrix::solve() +{ + // for (size_t i = 0; i < num_cells; i++) + // fprintf(stderr, "h_velocity_old[%d]: (%.15lf, %.15lf, %.15lf)\n", i, h_velocity_old[3*i], + // h_velocity_old[3*i + 1], h_velocity_old[3*i + 2]); + // constructor AmgXSolver at first interation + int nNz = num_cells + num_faces; // matrix entries + if (num_iteration == 0) // first interation + { + printf("Initializing AmgX Linear Solver\n"); + UxSolver->setOperator(num_cells, nNz, d_A_csr_row_index, d_A_csr_col_index, d_A_csr); + UySolver->setOperator(num_cells, nNz, d_A_csr_row_index, d_A_csr_col_index, d_A_csr + nNz); + UzSolver->setOperator(num_cells, nNz, d_A_csr_row_index, d_A_csr_col_index, d_A_csr + 2*nNz); + } + else + { + UxSolver->updateOperator(num_cells, nNz, d_A_csr); + UySolver->updateOperator(num_cells, nNz, d_A_csr + nNz); + UzSolver->updateOperator(num_cells, nNz, d_A_csr + 2*nNz); + } + UxSolver->solve(num_cells, d_psi, d_b); + UxSolver->solve(num_cells, d_psi + num_cells, d_b + num_cells); + UxSolver->solve(num_cells, d_psi + 2*num_cells, d_b + 2*num_cells); + num_iteration ++; + + checkCudaErrors(cudaMemcpy(h_psi, d_psi, cell_vec_bytes, cudaMemcpyDeviceToHost)); + // for (size_t i = 0; i < num_cells; i++) + // fprintf(stderr, "h_velocity_after[%d]: (%.15lf, %.15lf, %.15lf)\n", i, h_psi[i], + // h_psi[num_cells + i], h_psi[num_cells*2 + i]); +} +void dfMatrix::updatePsi(double* Psi) +{ + for (size_t i = 0; i < num_cells; i++) + { + Psi[i*3] = h_psi[i]; + Psi[i*3 + 1] = h_psi[num_cells + i]; + Psi[i*3 + 2] = h_psi[num_cells*2 + i]; + } +} From e9855399d69761fffee8d6262ea13e07c9a108c4 Mon Sep 17 00:00:00 2001 From: maorz1998 Date: Sun, 9 Apr 2023 18:15:36 +0800 Subject: [PATCH 2/3] fix bugs and perform preliminary optimization --- applications/solvers/dfLowMachFoam/UEqn.H | 74 +++-- .../solvers/dfLowMachFoam/createdfSolver.H | 12 +- .../solvers/dfLowMachFoam/dfLowMachFoam.C | 32 ++- bashrc.in | 4 + src_gpu/AmgXSolver.H | 1 + src_gpu/dfMatrix.H | 15 +- src_gpu/dfMatrix.cu | 262 +++++++----------- 7 files changed, 211 insertions(+), 189 deletions(-) diff --git a/applications/solvers/dfLowMachFoam/UEqn.H b/applications/solvers/dfLowMachFoam/UEqn.H index 6d9cb4c0c..6a97c5205 100644 --- a/applications/solvers/dfLowMachFoam/UEqn.H +++ b/applications/solvers/dfLowMachFoam/UEqn.H @@ -1,28 +1,45 @@ // Solve the Momentum equation - +start1 = std::clock(); +// Info << "U_of = " << U << endl; tmp tUEqn ( fvm::ddt(rho, U) + fvm::div(phi, U) - + turbulence->divDevRhoReff(U) + + turbulence->divDevRhoReff(U) == -fvc::grad(p) ); fvVectorMatrix& UEqn = tUEqn.ref(); +end1 = std::clock(); +time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC); +time_monitor_UEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC); -fvVectorMatrix turb_source -( - turbulence->divDevRhoReff(U) -); +// start1 = std::clock(); +// fvVectorMatrix turb_source +// ( +// turbulence->divDevRhoReff(U) +// ); +// end1 = std::clock(); +// time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC); +// time_monitor_UEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC); // UEqn.relax(); +start1 = std::clock(); +if (pimple.momentumPredictor()) +{ + solve(UEqn); + + K = 0.5*magSqr(U); +} +// Info << "U_of = " << U << endl; +end1 = std::clock(); +time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC); +time_monitor_UEqn_Solve += double(end1 - start1) / double(CLOCKS_PER_SEC); -// if (pimple.momentumPredictor()) -// { -// solve(UEqn == -fvc::grad(p)); -// K = 0.5*magSqr(U); -// } +start1 = std::clock(); +UEqn_GPU.fvm_ddt(&rho.oldTime()[0], &rho[0], &U.oldTime()[0][0]); +start2 = std::clock(); int offset = 0; forAll(U.boundaryField(), patchi) { @@ -45,18 +62,41 @@ forAll(U.boundaryField(), patchi) // boundary pressure std::copy(&patchP[0], &patchP[0]+patchSize, boundary_pressure_init.begin()+offset); - // boundary face vector - std::copy(&pSf[0][0], &pSf[0][0]+3*patchSize, boundary_face_vector_init.begin()+3*offset); offset += patchSize; } +end2 = std::clock(); +// time_monitor_CPU += double(end2 - start2) / double(CLOCKS_PER_SEC); + +UEqn_GPU.fvm_div(&phi[0], ueqn_internalCoeffs_init, ueqn_boundaryCoeffs_init); +UEqn_GPU.fvc_grad(&p[0], boundary_pressure_init); + +start2 = std::clock(); +fvVectorMatrix turb_source +( + turbulence->divDevRhoReff(U) +); +end2 = std::clock(); +// time_monitor_CPU += double(end2 - start2) / double(CLOCKS_PER_SEC); -UEqn_GPU.fvm_ddt(&rho.oldTime()[0], &rho[0], &mesh.V()[0], &U.oldTime()[0][0]); -UEqn_GPU.fvm_div(&mesh.surfaceInterpolation::weights()[0], &phi[0], ueqn_internalCoeffs_init, ueqn_boundaryCoeffs_init); -UEqn_GPU.fvc_grad(&mesh.Sf()[0][0], &p[0], boundary_face_vector_init, boundary_pressure_init); UEqn_GPU.add_fvMatrix(&turb_source.lower()[0], &turb_source.diag()[0], &turb_source.upper()[0], &turb_source.source()[0][0]); -// UEqn_GPU.print(); +end1 = std::clock(); +// time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC); +// time_monitor_UEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC); + +start1 = std::clock(); +// UEqn_GPU.checkValue(false); UEqn_GPU.solve(); +end1 = std::clock(); +// time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC); +// time_monitor_UEqn_Solve += double(end1 - start1) / double(CLOCKS_PER_SEC); + +start1 = std::clock(); UEqn_GPU.updatePsi(&U[0][0]); K = 0.5*magSqr(U); +end1 = std::clock(); +// time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC); +// time_monitor_UEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC); +// Info << "U_amgx = " << U << endl; +// time_monitor_UEqn_Solve += double(end1 - start1) / double(CLOCKS_PER_SEC); diff --git a/applications/solvers/dfLowMachFoam/createdfSolver.H b/applications/solvers/dfLowMachFoam/createdfSolver.H index 7b9bc7696..5cea30d05 100644 --- a/applications/solvers/dfLowMachFoam/createdfSolver.H +++ b/applications/solvers/dfLowMachFoam/createdfSolver.H @@ -4,19 +4,25 @@ int num_cells = mesh.nCells(); int num_surfaces = neighbour.size(); std::vector boundaryCellIndex; +std::vector boundary_face_vector_init; int num_boundary_faces = 0; int patchSize; +int offset = 0; forAll(mesh.boundary(), patchi) { labelUList sub_boundary = mesh.boundary()[patchi].faceCells(); patchSize = sub_boundary.size(); + const vectorField& pSf = mesh.Sf().boundaryField()[patchi]; + boundaryCellIndex.insert(boundaryCellIndex.end(), &sub_boundary[0], &sub_boundary[0]+patchSize); + boundary_face_vector_init.insert(boundary_face_vector_init.end()+offset, &pSf[0][0], &pSf[0][0]+3*patchSize); num_boundary_faces += patchSize; } int num_boundary_cells; -dfMatrix UEqn_GPU(num_surfaces, num_cells, num_boundary_faces, num_boundary_cells, &neighbour[0], &owner[0], boundaryCellIndex, "dDDI", - "/home/runze/deepflame-dev/examples/dfLowMachFoam/threeD_reactingTGV/H2/cvodeIntegrator/system/amgxUOptions"); +dfMatrix UEqn_GPU(num_surfaces, num_cells, num_boundary_faces, num_boundary_cells, &neighbour[0], &owner[0], &mesh.V()[0], &mesh.surfaceInterpolation::weights()[0], +&mesh.Sf()[0][0], boundary_face_vector_init, boundaryCellIndex, "dDDI", "/home/runze/deepflame-dev/examples/dfLowMachFoam/threeD_reactingTGV/H2/cvodeIntegrator/system/amgxUOptions"); std::vector ueqn_internalCoeffs_init(3*num_boundary_faces), ueqn_boundaryCoeffs_init(3*num_boundary_faces); -std::vector boundary_pressure_init(num_boundary_faces), boundary_face_vector_init(3*num_boundary_faces); \ No newline at end of file +std::vector boundary_pressure_init(num_boundary_faces); +// std::vector boundary_face_vector_init(3*num_boundary_faces); \ No newline at end of file diff --git a/applications/solvers/dfLowMachFoam/dfLowMachFoam.C b/applications/solvers/dfLowMachFoam/dfLowMachFoam.C index 28695533b..e1434b09f 100644 --- a/applications/solvers/dfLowMachFoam/dfLowMachFoam.C +++ b/applications/solvers/dfLowMachFoam/dfLowMachFoam.C @@ -83,13 +83,18 @@ int main(int argc, char *argv[]) #include "createRhoUfIfPresent.H" double time_monitor_flow=0; + double time_monitor_UEqn=0; + double time_monitor_UEqn_mtxAssembly=0; + double time_monitor_UEqn_Solve=0; + double time_monitor_UEqn_sum=0; double time_monitor_chem=0; double time_monitor_Y=0; double time_monitor_E=0; double time_monitor_corrThermo=0; double time_monitor_corrDiff=0; + double time_monitor_CPU=0; label timeIndex = 0; - clock_t start, end; + clock_t start, end, start1, end1, start2, end2; turbulence->validate(); @@ -99,7 +104,12 @@ int main(int argc, char *argv[]) #include "setInitialDeltaT.H" } + start1 = std::clock(); #include "createdfSolver.H" + end1 = std::clock(); + // time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC); + // time_monitor_UEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC); + // double time_UEqn_initial = time_monitor_UEqn; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -146,7 +156,7 @@ int main(int argc, char *argv[]) { #include "rhoEqn.H" } - + start = std::clock(); #include "UEqn.H" end = std::clock(); @@ -199,6 +209,7 @@ int main(int argc, char *argv[]) rho = thermo.rho(); runTime.write(); + time_monitor_UEqn_sum += time_monitor_UEqn; Info << "output time index " << runTime.timeIndex() << endl; @@ -209,11 +220,28 @@ int main(int argc, char *argv[]) Info<< "Energy Equations = " << time_monitor_E << " s" << endl; Info<< "thermo & Trans Properties = " << time_monitor_corrThermo << " s" << endl; Info<< "Diffusion Correction Time = " << time_monitor_corrDiff << " s" << endl; + Info<< "UEqn Time = " << time_monitor_UEqn << " s" << endl; + Info<< "UEqn Time assamble Mtx = " << time_monitor_UEqn_mtxAssembly << " s" << endl; + Info<< "UEqn Time solve = " << time_monitor_UEqn_Solve << " s" << endl; + // Info<< "UEqn sum Time = " << time_monitor_UEqn_sum << " s" << endl; + // Info<< "UEqn sum Time - overhead = " << time_monitor_UEqn_sum - time_UEqn_initial << " s" << endl; Info<< "sum Time = " << (time_monitor_chem + time_monitor_Y + time_monitor_flow + time_monitor_E + time_monitor_corrThermo + time_monitor_corrDiff) << " s" << endl; + Info<< "CPU Time (get turb souce) = " << time_monitor_CPU << " s" << endl; Info<< "============================================"< # include +# include // AmgX # include diff --git a/src_gpu/dfMatrix.H b/src_gpu/dfMatrix.H index 140488ae5..3b7f36cc1 100644 --- a/src_gpu/dfMatrix.H +++ b/src_gpu/dfMatrix.H @@ -9,6 +9,7 @@ #include #include #include +#include #include "AmgXSolver.H" @@ -147,24 +148,26 @@ private: AmgXSolver* UxSolver = nullptr, *UySolver = nullptr, *UzSolver = nullptr; int num_iteration; + double time_monitor_CPU, time_monitor_GPU; + public: dfMatrix(); dfMatrix(int num_surfaces, int num_cells, int num_boundary_faces, int &num_boundary_cells_output, - const int *neighbour, const int *owner, std::vector boundary_cell_id_init, + const int *neighbour, const int *owner, const double* volume, const double* weight, const double* face_vector, + std::vector boundary_face_vector_init, std::vector boundary_cell_id_init, const std::string &modeStr, const std::string &cfgFile); ~dfMatrix(); - void print(); + void checkValue(bool print); - void fvm_ddt(double *rho_old, double *rho_new, const double* volume, double* vector_old); + void fvm_ddt(double *rho_old, double *rho_new, double* vector_old); - void fvm_div(const double* weight, double* phi, std::vector ueqn_internalCoeffs_init, + void fvm_div(double* phi, std::vector ueqn_internalCoeffs_init, std::vector ueqn_boundaryCoeffs_init); - void fvc_grad(const double* face_vector, double* pressure, std::vector boundary_face_vector_init, - std::vector boundary_pressure_init); + void fvc_grad(double* pressure, std::vector boundary_pressure_init); void add_fvMatrix(double* turbSrc_low, double* turbSrc_diag, double* turbSrc_upp, double* turbSrc_source); diff --git a/src_gpu/dfMatrix.cu b/src_gpu/dfMatrix.cu index e56d858f5..e275ae361 100644 --- a/src_gpu/dfMatrix.cu +++ b/src_gpu/dfMatrix.cu @@ -1,4 +1,15 @@ #include "dfMatrix.H" + +void checkVectorEqual(int count, double* basevec, double* vec, double max_relative_error) { + for (size_t i = 0; i < count; ++i) + { + double abs_diff = fabs(basevec[i] - vec[i]); + double rel_diff = fabs(basevec[i] - vec[i]) / fabs(basevec[i]); + if (abs_diff > 1e-16 && rel_diff > max_relative_error) + fprintf(stderr, "mismatch index %d, cpu data: %.16lf, gpu data: %.16lf, relative error: %.16lf\n", i, basevec[i], vec[i], rel_diff); + } +} + // kernel functions __global__ void fvm_ddt_kernel(int num_cells, int num_faces, const double rdelta_t, const int* csr_row_index, const int* csr_diag_index, @@ -235,18 +246,19 @@ __global__ void add_fvMatrix_kernel(int num_cells, int num_faces, A_csr_output[csr_dim * 1 + i] = A_csr_input[csr_dim * 1 + i] + A_entry; A_csr_output[csr_dim * 2 + i] = A_csr_input[csr_dim * 2 + i] + A_entry; } - b_output[num_cells * 0 + index] = b_input[num_cells * 0 + index] + turbSrc_b[num_cells * 0 + index]; - b_output[num_cells * 1 + index] = b_input[num_cells * 1 + index] + turbSrc_b[num_cells * 1 + index]; - b_output[num_cells * 2 + index] = b_input[num_cells * 2 + index] + turbSrc_b[num_cells * 2 + index]; + b_output[num_cells * 0 + index] = b_input[num_cells * 0 + index] + turbSrc_b[index * 3 + 0]; + b_output[num_cells * 1 + index] = b_input[num_cells * 1 + index] + turbSrc_b[index * 3 + 1]; + b_output[num_cells * 2 + index] = b_input[num_cells * 2 + index] + turbSrc_b[index * 3 + 2]; } // constructor (construct mesh variable) dfMatrix::dfMatrix(){} dfMatrix::dfMatrix(int num_surfaces, int num_cells, int num_boundary_faces, int & num_boundary_cells_output, - const int *neighbour, const int *owner, std::vector boundary_cell_id_init, const std::string &modeStr, + const int *neighbour, const int *owner, const double* volume, const double* weight, const double* face_vector, + std::vector boundary_face_vector_init, std::vector boundary_cell_id_init, const std::string &modeStr, const std::string &cfgFile) : num_cells(num_cells), num_faces(num_surfaces*2), num_surfaces(num_surfaces), num_iteration(0), - num_boundary_faces(num_boundary_faces) + num_boundary_faces(num_boundary_faces), h_volume(volume), time_monitor_CPU(0.), time_monitor_GPU(0.) { // resize vectors h_weight_vec_init.resize(num_faces); @@ -473,6 +485,29 @@ dfMatrix::dfMatrix(int num_surfaces, int num_cells, int num_boundary_faces, int return pair.second; }); + // copy and permutate cell variables + std::copy(weight, weight + num_surfaces, h_weight_vec_init.begin()); + std::copy(weight, weight + num_surfaces, h_weight_vec_init.begin() + num_surfaces); + std::copy(face_vector, face_vector + 3*num_surfaces, h_face_vector_vec_init.begin()); + std::copy(face_vector, face_vector + 3*num_surfaces, h_face_vector_vec_init.begin() + 3*num_surfaces); + for (int i = 0; i < num_faces; i++) + { + h_weight_vec[i] = h_weight_vec_init[permedIndex[i]]; + h_face_vector_vec[i*3] = h_face_vector_vec_init[3*permedIndex[i]]; + h_face_vector_vec[i*3+1] = h_face_vector_vec_init[3*permedIndex[i]+1]; + h_face_vector_vec[i*3+2] = h_face_vector_vec_init[3*permedIndex[i]+2]; + } + h_weight = h_weight_vec.data(); + h_face_vector = h_face_vector_vec.data(); + + for (int i = 0; i < num_boundary_faces; i++) + { + boundary_face_vector[3*i] = boundary_face_vector_init[3*boundPermutationList[i]]; + boundary_face_vector[3*i+1] = boundary_face_vector_init[3*boundPermutationList[i]+1]; + boundary_face_vector[3*i+2] = boundary_face_vector_init[3*boundPermutationList[i]+2]; + } + h_boundary_face_vector = boundary_face_vector.data(); + /************************allocate memory on device****************************/ int total_bytes = 0; @@ -523,6 +558,11 @@ dfMatrix::dfMatrix(int num_surfaces, int num_cells, int num_boundary_faces, int checkCudaErrors(cudaMemcpyAsync(d_boundary_cell_offset, h_boundary_cell_offset, (num_boundary_cells+1) * sizeof(int), cudaMemcpyHostToDevice, stream)); checkCudaErrors(cudaMemcpyAsync(d_boundary_cell_id, h_boundary_cell_id, boundary_face_index_bytes, cudaMemcpyHostToDevice, stream)); + checkCudaErrors(cudaMemcpyAsync(d_volume, h_volume, cell_bytes, cudaMemcpyHostToDevice, stream)); + checkCudaErrors(cudaMemcpyAsync(d_weight, h_weight, face_bytes, cudaMemcpyHostToDevice, stream)); + checkCudaErrors(cudaMemcpyAsync(d_face_vector, h_face_vector, face_vec_bytes, cudaMemcpyHostToDevice, stream)); + checkCudaErrors(cudaMemcpyAsync(d_boundary_face_vector, h_boundary_face_vector, boundary_face_vec_bytes, cudaMemcpyHostToDevice, stream)); + // construct AmgXSolver UxSolver = new AmgXSolver(modeStr, cfgFile); UySolver = new AmgXSolver(modeStr, cfgFile); @@ -533,19 +573,17 @@ dfMatrix::~dfMatrix() { } -void dfMatrix::fvm_ddt(double *rho_old, double *rho_new, const double* volume, - double* vector_old) +void dfMatrix::fvm_ddt(double *rho_old, double *rho_new, double* vector_old) { + clock_t start = std::clock(); // copy cell variables directly h_rho_new = rho_new; h_rho_old = rho_old; - h_volume = volume; h_velocity_old = vector_old; // Copy the host input array in host memory to the device input array in device memory checkCudaErrors(cudaMemcpyAsync(d_rho_old, h_rho_old, cell_bytes, cudaMemcpyHostToDevice, stream)); checkCudaErrors(cudaMemcpyAsync(d_rho_new, h_rho_new, cell_bytes, cudaMemcpyHostToDevice, stream)); - checkCudaErrors(cudaMemcpyAsync(d_volume, h_volume, cell_bytes, cudaMemcpyHostToDevice, stream)); checkCudaErrors(cudaMemcpyAsync(d_velocity_old, h_velocity_old, cell_vec_bytes, cudaMemcpyHostToDevice, stream)); // launch cuda kernel @@ -555,25 +593,22 @@ void dfMatrix::fvm_ddt(double *rho_old, double *rho_new, const double* volume, fvm_ddt_kernel<<>>(num_cells, num_faces, rdelta_t, d_A_csr_row_index, d_A_csr_diag_index, d_rho_old, d_rho_new, d_volume, d_velocity_old, d_A_csr, d_b, d_A_csr, d_b, d_psi); - // Synchronize stream - checkCudaErrors(cudaStreamSynchronize(stream)); + clock_t end = std::clock(); + time_monitor_GPU += double(end - start) / double(CLOCKS_PER_SEC); } -void dfMatrix::fvm_div(const double* weight, double* phi, std::vector ueqn_internalCoeffs_init, +void dfMatrix::fvm_div(double* phi, std::vector ueqn_internalCoeffs_init, std::vector ueqn_boundaryCoeffs_init) { // copy and permutate face variables + clock_t start = std::clock(); std::copy(phi, phi + num_surfaces, h_phi_vec_init.begin()); std::copy(phi, phi + num_surfaces, h_phi_vec_init.begin() + num_surfaces); - std::copy(weight, weight + num_surfaces, h_weight_vec_init.begin()); - std::copy(weight, weight + num_surfaces, h_weight_vec_init.begin() + num_surfaces); for (int i = 0; i < num_faces; i++) { - h_weight_vec[i] = h_weight_vec_init[permedIndex[i]]; h_phi_vec[i] = h_phi_vec_init[permedIndex[i]]; } - h_weight = h_weight_vec.data(); h_phi = h_phi_vec.data(); // copy and permutate boundary variable @@ -588,9 +623,11 @@ void dfMatrix::fvm_div(const double* weight, double* phi, std::vector ue } h_internal_coeffs = ueqn_internalCoeffs.data(); h_boundary_coeffs = ueqn_boundaryCoeffs.data(); + clock_t end = std::clock(); + time_monitor_CPU += double(end - start) / double(CLOCKS_PER_SEC); + start = std::clock(); // Copy the host input array in host memory to the device input array in device memory - checkCudaErrors(cudaMemcpyAsync(d_weight, h_weight, face_bytes, cudaMemcpyHostToDevice, stream)); checkCudaErrors(cudaMemcpyAsync(d_phi, h_phi, face_bytes, cudaMemcpyHostToDevice, stream)); checkCudaErrors(cudaMemcpyAsync(d_internal_coeffs, h_internal_coeffs, boundary_face_vec_bytes, cudaMemcpyHostToDevice, stream)); checkCudaErrors(cudaMemcpyAsync(d_boundary_coeffs, h_boundary_coeffs, boundary_face_vec_bytes, cudaMemcpyHostToDevice, stream)); @@ -606,46 +643,29 @@ void dfMatrix::fvm_div(const double* weight, double* phi, std::vector ue d_A_csr_row_index, d_A_csr_diag_index, d_boundary_cell_offset, d_boundary_cell_id, d_internal_coeffs, d_boundary_coeffs, d_A_csr, d_b, d_A_csr, d_b); - - // Synchronize stream - checkCudaErrors(cudaStreamSynchronize(stream)); + end = std::clock(); + time_monitor_GPU += double(end - start) / double(CLOCKS_PER_SEC); } -void dfMatrix::fvc_grad(const double* face_vector, double* pressure, std::vector boundary_face_vector_init, - std::vector boundary_pressure_init) +void dfMatrix::fvc_grad(double* pressure, std::vector boundary_pressure_init) { // copy cell variables directly + clock_t start = std::clock(); h_pressure = pressure; - - // copy and permutate face variables - std::copy(face_vector, face_vector + 3*num_surfaces, h_face_vector_vec_init.begin()); - std::copy(face_vector, face_vector + 3*num_surfaces, h_face_vector_vec_init.begin() + 3*num_surfaces); - for (int i = 0; i < num_faces; i++) - { - h_face_vector_vec[i*3] = h_face_vector_vec_init[3*permedIndex[i]]; - h_face_vector_vec[i*3+1] = h_face_vector_vec_init[3*permedIndex[i]+1]; - h_face_vector_vec[i*3+2] = h_face_vector_vec_init[3*permedIndex[i]+2]; - } - h_face_vector = h_face_vector_vec.data(); - // copy and permutate boundary variable for (int i = 0; i < num_boundary_faces; i++) { - boundary_face_vector[3*i] = boundary_face_vector_init[3*boundPermutationList[i]]; - boundary_face_vector[3*i+1] = boundary_face_vector_init[3*boundPermutationList[i]+1]; - boundary_face_vector[3*i+2] = boundary_face_vector_init[3*boundPermutationList[i]+2]; boundary_pressure[i] = boundary_pressure_init[boundPermutationList[i]]; } - h_boundary_face_vector = boundary_face_vector.data(); h_boundary_pressure = boundary_pressure.data(); + clock_t end = std::clock(); + time_monitor_CPU += double(end - start) / double(CLOCKS_PER_SEC); // Copy the host input array in host memory to the device input array in device memory - checkCudaErrors(cudaMemcpyAsync(d_face_vector, h_face_vector, face_vec_bytes, cudaMemcpyHostToDevice, stream)); + start = std::clock(); checkCudaErrors(cudaMemcpyAsync(d_pressure, h_pressure, cell_bytes, cudaMemcpyHostToDevice, stream)); checkCudaErrors(cudaMemcpyAsync(d_boundary_pressure, h_boundary_pressure, boundary_face_bytes, cudaMemcpyHostToDevice, stream)); - checkCudaErrors(cudaMemcpyAsync(d_boundary_face_vector, h_boundary_face_vector, boundary_face_vec_bytes, cudaMemcpyHostToDevice, stream)); - - + // launch cuda kernel size_t threads_per_block = 1024; size_t blocks_per_grid = (num_cells + threads_per_block - 1) / threads_per_block; @@ -656,14 +676,14 @@ void dfMatrix::fvc_grad(const double* face_vector, double* pressure, std::vector fvc_grad_boundary_face<<>>(num_cells, num_boundary_cells, d_boundary_cell_offset, d_boundary_cell_id, d_boundary_face_vector, d_boundary_pressure, d_b, d_b); - - // Synchronize stream - checkCudaErrors(cudaStreamSynchronize(stream)); + end = std::clock(); + time_monitor_GPU += double(end - start) / double(CLOCKS_PER_SEC); } void dfMatrix::add_fvMatrix(double* turbSrc_low, double* turbSrc_diag, double* turbSrc_upp, double* turbSrc_source) { // copy and permutate matrix variables + clock_t start = std::clock(); std::copy(turbSrc_low, turbSrc_low + num_surfaces, h_turbSrc_init_mtx_vec.begin()); std::copy(turbSrc_diag, turbSrc_diag + num_cells, h_turbSrc_init_mtx_vec.begin() + num_surfaces); std::copy(turbSrc_upp, turbSrc_upp + num_surfaces, h_turbSrc_init_mtx_vec.begin() + num_surfaces + num_cells); @@ -671,41 +691,35 @@ void dfMatrix::add_fvMatrix(double* turbSrc_low, double* turbSrc_diag, double* t // permutate for (int i = 0; i < num_cells+2*num_surfaces; i++) h_turbSrc_init_1mtx[i] = h_turbSrc_init_mtx_vec[tmpPermutatedList[i]]; - // fill RHS_x - for (int i = 0; i < 3*num_cells; i+=3) - h_turbSrc_src_vec.push_back(h_turbSrc_init_src_vec[i]); - // fill RHS_y - for (int i = 1; i < 3*num_cells; i+=3) - h_turbSrc_src_vec.push_back(h_turbSrc_init_src_vec[i]); - // fill RHS_z - for (int i = 2; i < 3*num_cells; i+=3) - h_turbSrc_src_vec.push_back(h_turbSrc_init_src_vec[i]); + clock_t end = std::clock(); + time_monitor_CPU += double(end - start) / double(CLOCKS_PER_SEC); // Copy the host input array in host memory to the device input array in device memory + start = std::clock(); checkCudaErrors(cudaMemcpyAsync(d_turbSrc_A, h_turbSrc_init_1mtx.data(), csr_value_bytes, cudaMemcpyHostToDevice, stream)); - checkCudaErrors(cudaMemcpyAsync(d_turbSrc_b, h_turbSrc_src_vec.data(), cell_vec_bytes, cudaMemcpyHostToDevice, stream)); + checkCudaErrors(cudaMemcpyAsync(d_turbSrc_b, h_turbSrc_init_src_vec.data(), cell_vec_bytes, cudaMemcpyHostToDevice, stream)); size_t threads_per_block = 1024; size_t blocks_per_grid = (num_cells + threads_per_block - 1) / threads_per_block; add_fvMatrix_kernel<<>>(num_cells, num_faces, d_A_csr_row_index, d_turbSrc_A, d_turbSrc_b, d_A_csr, d_b, d_A_csr, d_b); - - // Synchronize stream - checkCudaErrors(cudaStreamSynchronize(stream)); + end = std::clock(); + time_monitor_GPU += double(end - start) / double(CLOCKS_PER_SEC); } -void dfMatrix::print() +void dfMatrix::checkValue(bool print) { checkCudaErrors(cudaMemcpyAsync(h_A_csr, d_A_csr, csr_value_vec_bytes, cudaMemcpyDeviceToHost, stream)); checkCudaErrors(cudaMemcpyAsync(h_b, d_b, cell_vec_bytes, cudaMemcpyDeviceToHost, stream)); // Synchronize stream checkCudaErrors(cudaStreamSynchronize(stream)); - - for (int i = 0; i < (2*num_surfaces + num_cells); i++) - fprintf(stderr, "h_A_csr[%d]: %.15lf\n", i, h_A_csr[i]); - - for (int i = 0; i < num_cells * 3; i++) - fprintf(stderr, "h_b[%d]: %.15lf\n", i, h_b[i]); + if (print) + { + for (int i = 0; i < (2*num_surfaces + num_cells); i++) + fprintf(stderr, "h_A_csr[%d]: %.15lf\n", i, h_A_csr[i]); + for (int i = 0; i < num_cells * 3; i++) + fprintf(stderr, "h_b[%d]: %.15lf\n", i, h_b[i]); + } char *input_file = "of_output.txt"; FILE *fp = fopen(input_file, "rb+"); @@ -726,8 +740,6 @@ void dfMatrix::print() { h_A_of_vec_1mtx[i] = h_A_of_init_vec[tmpPermutatedList[i]]; } - for (int i = 0; i < (2*num_surfaces + num_cells); i++) - fprintf(stderr, "h_A_of_vec_1mtx[%d]: %.15lf\n", i, h_A_of_vec_1mtx[i]); std::vector h_A_of_vec((2*num_surfaces + num_cells)*3); for (int i =0; i < 3; i ++) @@ -753,100 +765,21 @@ void dfMatrix::print() { h_b_of_vec.push_back(h_b_of_init_vec[i]); } - for (int i = 0; i < 3*num_cells; i++) - fprintf(stderr, "h_b_of_vec[%d]: %.15lf\n", i, h_b_of_vec[i]); -} -// void dfMatrix::solve(const std::string &cfgFile) -// { -// // checkCudaErrors(cudaMemcpyAsync(h_A_csr, d_A_csr, csr_value_vec_bytes, cudaMemcpyDeviceToHost, stream)); -// // checkCudaErrors(cudaMemcpyAsync(h_b, d_b, cell_vec_bytes, cudaMemcpyDeviceToHost, stream)); -// checkCudaErrors(cudaMemcpyAsync(h_psi, d_psi, cell_vec_bytes, cudaMemcpyDeviceToHost, stream)); -// /** -// * @brief initialize AmgXSolver -// */ -// // - setMode - -// // - initAmgX(cfgFile) -// // initialize AmgX -// AMGX_SAFE_CALL(AMGX_initialize()); -// // intialize AmgX plugings -// AMGX_SAFE_CALL(AMGX_initialize_plugins()); -// // let AmgX to handle errors returned -// AMGX_SAFE_CALL(AMGX_install_signal_handler()); -// // create an AmgX configure object -// AMGX_config_handle cfg = nullptr; -// AMGX_SAFE_CALL(AMGX_config_create_from_file(&cfg, cfgFile.c_str())); -// // let AmgX handle returned error codes internally -// AMGX_SAFE_CALL(AMGX_config_add_parameters(&cfg, "exception_handling=1")); -// // initialize rsrc to nullptr; -// AMGX_resources_handle rsrc = nullptr; -// // create an AmgX resource object, only the first instance is in charge -// // if (solver_instance_count == 1) -// AMGX_resources_create_simple(&rsrc, cfg); // remember to set deviceID - -// // create AmgX vector object for unknowns and RHS -// AMGX_vector_handle AmgXP = nullptr; -// AMGX_vector_handle AmgXRHS = nullptr; -// AMGX_Mode mode = AMGX_mode_dDDI; -// AMGX_vector_create(&AmgXP, rsrc, mode); -// AMGX_vector_create(&AmgXRHS, rsrc, mode); - -// // create AmgX matrix object for unknowns and RHS -// AMGX_matrix_handle AmgXA = nullptr; -// AMGX_matrix_create(&AmgXA, rsrc, mode); - -// // bind the matrix A to the solver -// AMGX_solver_handle solver = nullptr; -// AMGX_solver_create(&solver, rsrc, mode, cfg); - -// int ring; -// // obtain the default number of rings based on current configuration -// AMGX_config_get_default_number_of_rings(cfg, &ring); - -// AMGX_matrix_upload_all(AmgXA, num_cells, (num_cells + num_faces), 1, 1, d_A_csr_row_index, d_A_csr_col_index, -// d_A_csr, nullptr); - -// AMGX_solver_setup(solver, AmgXA); - + if (print) + { + for (int i = 0; i < (2*num_surfaces + num_cells); i++) + fprintf(stderr, "h_A_of_vec_1mtx[%d]: %.15lf\n", i, h_A_of_vec_1mtx[i]); + for (int i = 0; i < 3*num_cells; i++) + fprintf(stderr, "h_b_of_vec[%d]: %.15lf\n", i, h_b_of_vec[i]); + } -// // bind -// AMGX_vector_bind(AmgXP, AmgXA); -// AMGX_vector_bind(AmgXP, AmgXA); - -// // create an AmgX solver object - -// /** -// * @brief amgx.initialiseMatrixComms(Amat); -// */ - -// /** -// * @brief amgx.setOperator -// */ -// AMGX_vector_upload(AmgXRHS, num_cells, 1, d_b); -// AMGX_vector_upload(AmgXP, num_cells, 1, d_psi); - -// // Solve -// AMGX_solver_solve(solver, AmgXRHS, AmgXP); -// AMGX_SOLVE_STATUS status; -// AMGX_solver_get_status(solver, &status); - -// if (status != AMGX_SOLVE_SUCCESS) -// { -// fprintf(stderr, "AmgX solver failed to solve the system! " -// "The error code is %d.\n", -// status); -// } -// // Download data from device -// double * h_psi_solved; -// h_psi_solved = new double[cell_vec_bytes]; -// AMGX_vector_download(AmgXP, h_psi_solved); -// for (int i = 0; i < num_cells; i++) -// { -// fprintf(stderr, "h_psi_solved[%d] = %lf\n", i, h_psi_solved[i]); -// } - -// } + // check + fprintf(stderr, "check of h_A_csr\n"); + checkVectorEqual(2*num_surfaces + num_cells, h_A_of_vec_1mtx.data(), h_A_csr, 1e-5); + fprintf(stderr, "check of h_b\n"); + checkVectorEqual(3*num_cells, h_b_of_vec.data(), h_b, 1e-5); +} void dfMatrix::solve() { @@ -854,6 +787,13 @@ void dfMatrix::solve() // fprintf(stderr, "h_velocity_old[%d]: (%.15lf, %.15lf, %.15lf)\n", i, h_velocity_old[3*i], // h_velocity_old[3*i + 1], h_velocity_old[3*i + 2]); // constructor AmgXSolver at first interation + // Synchronize stream + checkCudaErrors(cudaStreamSynchronize(stream)); + printf("CPU Time (copy&permutate) = %.6lf s\n", time_monitor_CPU); + printf("GPU Time = %.6lf s\n", time_monitor_GPU); + time_monitor_CPU = 0; + time_monitor_GPU = 0; + int nNz = num_cells + num_faces; // matrix entries if (num_iteration == 0) // first interation { @@ -869,8 +809,8 @@ void dfMatrix::solve() UzSolver->updateOperator(num_cells, nNz, d_A_csr + 2*nNz); } UxSolver->solve(num_cells, d_psi, d_b); - UxSolver->solve(num_cells, d_psi + num_cells, d_b + num_cells); - UxSolver->solve(num_cells, d_psi + 2*num_cells, d_b + 2*num_cells); + UySolver->solve(num_cells, d_psi + num_cells, d_b + num_cells); + UzSolver->solve(num_cells, d_psi + 2*num_cells, d_b + 2*num_cells); num_iteration ++; checkCudaErrors(cudaMemcpy(h_psi, d_psi, cell_vec_bytes, cudaMemcpyDeviceToHost)); From 0ffd0f2628e72aa3bc58a379f2ce41078b554dfd Mon Sep 17 00:00:00 2001 From: maorz1998 Date: Sun, 9 Apr 2023 22:51:47 +0800 Subject: [PATCH 3/3] use GPU to accelerate permulation --- applications/solvers/dfLowMachFoam/UEqn.H | 51 ++++----- src_gpu/dfMatrix.H | 14 ++- src_gpu/dfMatrix.cu | 126 ++++++++++++++-------- 3 files changed, 118 insertions(+), 73 deletions(-) diff --git a/applications/solvers/dfLowMachFoam/UEqn.H b/applications/solvers/dfLowMachFoam/UEqn.H index 6a97c5205..c90c7bdc0 100644 --- a/applications/solvers/dfLowMachFoam/UEqn.H +++ b/applications/solvers/dfLowMachFoam/UEqn.H @@ -1,5 +1,5 @@ // Solve the Momentum equation -start1 = std::clock(); +// start1 = std::clock(); // Info << "U_of = " << U << endl; tmp tUEqn ( @@ -7,9 +7,9 @@ tmp tUEqn + turbulence->divDevRhoReff(U) == -fvc::grad(p) ); fvVectorMatrix& UEqn = tUEqn.ref(); -end1 = std::clock(); -time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC); -time_monitor_UEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC); +// end1 = std::clock(); +// time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC); +// time_monitor_UEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC); // start1 = std::clock(); // fvVectorMatrix turb_source @@ -22,17 +22,17 @@ time_monitor_UEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC); // UEqn.relax(); -start1 = std::clock(); -if (pimple.momentumPredictor()) -{ - solve(UEqn); +// start1 = std::clock(); +// if (pimple.momentumPredictor()) +// { +// solve(UEqn); - K = 0.5*magSqr(U); -} +// K = 0.5*magSqr(U); +// } // Info << "U_of = " << U << endl; -end1 = std::clock(); -time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC); -time_monitor_UEqn_Solve += double(end1 - start1) / double(CLOCKS_PER_SEC); +// end1 = std::clock(); +// time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC); +// time_monitor_UEqn_Solve += double(end1 - start1) / double(CLOCKS_PER_SEC); start1 = std::clock(); @@ -66,10 +66,10 @@ forAll(U.boundaryField(), patchi) offset += patchSize; } end2 = std::clock(); -// time_monitor_CPU += double(end2 - start2) / double(CLOCKS_PER_SEC); +time_monitor_CPU += double(end2 - start2) / double(CLOCKS_PER_SEC); -UEqn_GPU.fvm_div(&phi[0], ueqn_internalCoeffs_init, ueqn_boundaryCoeffs_init); -UEqn_GPU.fvc_grad(&p[0], boundary_pressure_init); +UEqn_GPU.fvm_div(&phi[0], ueqn_internalCoeffs_init, ueqn_boundaryCoeffs_init, boundary_pressure_init); +UEqn_GPU.fvc_grad(&p[0]); start2 = std::clock(); fvVectorMatrix turb_source @@ -77,26 +77,27 @@ fvVectorMatrix turb_source turbulence->divDevRhoReff(U) ); end2 = std::clock(); -// time_monitor_CPU += double(end2 - start2) / double(CLOCKS_PER_SEC); +time_monitor_CPU += double(end2 - start2) / double(CLOCKS_PER_SEC); UEqn_GPU.add_fvMatrix(&turb_source.lower()[0], &turb_source.diag()[0], &turb_source.upper()[0], &turb_source.source()[0][0]); end1 = std::clock(); -// time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC); -// time_monitor_UEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC); +time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC); +time_monitor_UEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC); -start1 = std::clock(); // UEqn_GPU.checkValue(false); + +start1 = std::clock(); UEqn_GPU.solve(); end1 = std::clock(); -// time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC); -// time_monitor_UEqn_Solve += double(end1 - start1) / double(CLOCKS_PER_SEC); +time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC); +time_monitor_UEqn_Solve += double(end1 - start1) / double(CLOCKS_PER_SEC); start1 = std::clock(); UEqn_GPU.updatePsi(&U[0][0]); K = 0.5*magSqr(U); end1 = std::clock(); -// time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC); -// time_monitor_UEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC); +time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC); +time_monitor_UEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC); +time_monitor_CPU += double(end1 - start1) / double(CLOCKS_PER_SEC); // Info << "U_amgx = " << U << endl; -// time_monitor_UEqn_Solve += double(end1 - start1) / double(CLOCKS_PER_SEC); diff --git a/src_gpu/dfMatrix.H b/src_gpu/dfMatrix.H index 3b7f36cc1..6c94fae9a 100644 --- a/src_gpu/dfMatrix.H +++ b/src_gpu/dfMatrix.H @@ -93,7 +93,9 @@ private: double *h_internal_coeffs = nullptr, *h_boundary_coeffs = nullptr, *h_boundary_pressure = nullptr, *h_boundary_face_vector = nullptr, *d_internal_coeffs = nullptr, *d_boundary_coeffs = nullptr, - *d_boundary_pressure = nullptr, *d_boundary_face_vector = nullptr; + *d_internal_coeffs_init = nullptr, *d_boundary_coeffs_init = nullptr, + *d_boundary_pressure = nullptr, *d_boundary_face_vector = nullptr, + *d_boundary_pressure_init = nullptr; std::vector boundPermutationList; std::vector ueqn_internalCoeffs, ueqn_boundaryCoeffs; std::vector boundary_face_vector; @@ -102,6 +104,8 @@ private: // - the device pointer to the permutated index list std::vector permedIndex; int *d_permedIndex=nullptr; + int *d_bouPermedIndex = nullptr; + // bytesize // - bytes of diagnal entries @@ -113,6 +117,7 @@ private: // - bytes of diagnal index size_t face_bytes; size_t face_vec_bytes; + size_t face_index_bytes; size_t boundary_cell_bytes; size_t boundary_cell_vec_bytes; @@ -129,10 +134,11 @@ private: size_t csr_value_vec_bytes; // extra matrix information - double *d_turbSrc_A = nullptr, *d_turbSrc_b = nullptr; + double *d_turbSrc_A = nullptr, *d_turbSrc_b = nullptr, *d_turbSrc_A_init = nullptr; std::vector h_turbSrc_init_mtx_vec, h_turbSrc_init_1mtx; std::vector h_turbSrc_init_src_vec, h_turbSrc_src_vec; std::vector tmpPermutatedList; + int * d_tmpPermutatedList = nullptr; double *h_A_csr = nullptr, *h_b = nullptr, *h_psi = nullptr; double *d_A_csr = nullptr, *d_b = nullptr, *d_psi = nullptr; @@ -165,9 +171,9 @@ public: void fvm_ddt(double *rho_old, double *rho_new, double* vector_old); void fvm_div(double* phi, std::vector ueqn_internalCoeffs_init, - std::vector ueqn_boundaryCoeffs_init); + std::vector ueqn_boundaryCoeffs_init, std::vector boundary_pressure_init); - void fvc_grad(double* pressure, std::vector boundary_pressure_init); + void fvc_grad(double* pressure); void add_fvMatrix(double* turbSrc_low, double* turbSrc_diag, double* turbSrc_upp, double* turbSrc_source); diff --git a/src_gpu/dfMatrix.cu b/src_gpu/dfMatrix.cu index e275ae361..91b3aebd5 100644 --- a/src_gpu/dfMatrix.cu +++ b/src_gpu/dfMatrix.cu @@ -251,6 +251,44 @@ __global__ void add_fvMatrix_kernel(int num_cells, int num_faces, b_output[num_cells * 2 + index] = b_input[num_cells * 2 + index] + turbSrc_b[index * 3 + 2]; } +__global__ void offdiagPermutation(const int num_faces, const int *permedIndex, + const double *d_phi_init, double *d_phi) +{ + int index = blockDim.x * blockIdx.x + threadIdx.x; + if (index >= num_faces) return; + + int p = permedIndex[index]; + d_phi[index] = d_phi_init[p]; +} + +__global__ void boundaryPermutation(const int num_boundary_faces, const int *bouPermedIndex, + const double *ueqn_internalCoeffs_init, const double *ueqn_boundaryCoeffs_init, + const double *boundary_pressure_init, double *ueqn_internalCoeffs, + double *ueqn_boundaryCoeffs, double *boundary_pressure) +{ + int index = blockDim.x * blockIdx.x + threadIdx.x; + if (index >= num_boundary_faces) return; + + int p = bouPermedIndex[index]; + ueqn_internalCoeffs[3*index] = ueqn_internalCoeffs_init[3*p]; + ueqn_internalCoeffs[3*index + 1] = ueqn_internalCoeffs_init[3*p + 1]; + ueqn_internalCoeffs[3*index + 2] = ueqn_internalCoeffs_init[3*p + 2]; + ueqn_boundaryCoeffs[3*index] = ueqn_boundaryCoeffs_init[3*p]; + ueqn_boundaryCoeffs[3*index + 1] = ueqn_boundaryCoeffs_init[3*p + 1]; + ueqn_boundaryCoeffs[3*index + 2] = ueqn_boundaryCoeffs_init[3*p + 2]; + boundary_pressure[index] = boundary_pressure_init[p]; +} + +__global__ void csrPermutation(const int num_Nz, const int *csrPermedIndex, + const double *d_turbSrc_A_init, double *d_turbSrc_A) +{ + int index = blockDim.x * blockIdx.x + threadIdx.x; + if (index >= num_Nz) return; + + int p = csrPermedIndex[index]; + d_turbSrc_A[index] = d_turbSrc_A_init[p]; +} + // constructor (construct mesh variable) dfMatrix::dfMatrix(){} dfMatrix::dfMatrix(int num_surfaces, int num_cells, int num_boundary_faces, int & num_boundary_cells_output, @@ -279,6 +317,7 @@ dfMatrix::dfMatrix(int num_surfaces, int num_cells, int num_boundary_faces, int face_bytes = num_faces * sizeof(double); face_vec_bytes = num_faces * 3 * sizeof(double); + face_index_bytes = num_faces * sizeof(int); // A_csr has one more element in each row: itself csr_row_index_bytes = (num_cells + 1) * sizeof(int); @@ -523,18 +562,22 @@ dfMatrix::dfMatrix(int num_surfaces, int num_cells, int num_boundary_faces, int checkCudaErrors(cudaMalloc((void**)&d_velocity_old, cell_vec_bytes)); checkCudaErrors(cudaMalloc((void**)&d_weight, face_bytes)); checkCudaErrors(cudaMalloc((void**)&d_phi, face_bytes)); + checkCudaErrors(cudaMalloc((void**)&d_phi_init, face_bytes)); checkCudaErrors(cudaMalloc((void**)&d_face_vector, face_vec_bytes)); - total_bytes += (cell_bytes * 4 + face_bytes * 2 + cell_vec_bytes + face_vec_bytes); + total_bytes += (cell_bytes * 4 + face_bytes * 3 + cell_vec_bytes + face_vec_bytes); checkCudaErrors(cudaMalloc((void**)&d_boundary_cell_offset, (num_boundary_cells+1) * sizeof(int))); checkCudaErrors(cudaMalloc((void**)&d_boundary_cell_id, boundary_face_index_bytes)); total_bytes += (boundary_face_index_bytes + (num_boundary_cells+1) * sizeof(int)); + checkCudaErrors(cudaMalloc((void**)&d_boundary_pressure_init, boundary_face_bytes)); checkCudaErrors(cudaMalloc((void**)&d_boundary_pressure, boundary_face_bytes)); checkCudaErrors(cudaMalloc((void**)&d_boundary_face_vector, boundary_face_vec_bytes)); + checkCudaErrors(cudaMalloc((void**)&d_internal_coeffs_init, boundary_face_vec_bytes)); + checkCudaErrors(cudaMalloc((void**)&d_boundary_coeffs_init, boundary_face_vec_bytes)); checkCudaErrors(cudaMalloc((void**)&d_internal_coeffs, boundary_face_vec_bytes)); checkCudaErrors(cudaMalloc((void**)&d_boundary_coeffs, boundary_face_vec_bytes)); - total_bytes += (boundary_face_bytes + boundary_face_vec_bytes * 3); + total_bytes += (boundary_face_bytes*2 + boundary_face_vec_bytes * 5); checkCudaErrors(cudaMalloc((void**)&d_A_csr, csr_value_vec_bytes)); checkCudaErrors(cudaMalloc((void**)&d_b, cell_vec_bytes)); @@ -542,8 +585,14 @@ dfMatrix::dfMatrix(int num_surfaces, int num_cells, int num_boundary_faces, int total_bytes += (boundary_face_bytes + boundary_face_vec_bytes * 3); checkCudaErrors(cudaMalloc((void**)&d_turbSrc_A, csr_value_bytes)); + checkCudaErrors(cudaMalloc((void**)&d_turbSrc_A_init, csr_value_bytes)); checkCudaErrors(cudaMalloc((void**)&d_turbSrc_b, cell_vec_bytes)); - total_bytes += (csr_value_bytes + cell_vec_bytes); + total_bytes += (2*csr_value_bytes + cell_vec_bytes); + + checkCudaErrors(cudaMalloc((void**)&d_permedIndex, face_index_bytes)); + checkCudaErrors(cudaMalloc((void**)&d_bouPermedIndex, boundary_face_index_bytes)); + checkCudaErrors(cudaMalloc((void**)&d_tmpPermutatedList, csr_col_index_bytes)); + total_bytes += (face_index_bytes + boundary_face_index_bytes + csr_col_index_bytes); fprintf(stderr, "Total bytes malloc on GPU: %.2fMB\n", total_bytes * 1.0 / 1024 / 1024); @@ -563,6 +612,10 @@ dfMatrix::dfMatrix(int num_surfaces, int num_cells, int num_boundary_faces, int checkCudaErrors(cudaMemcpyAsync(d_face_vector, h_face_vector, face_vec_bytes, cudaMemcpyHostToDevice, stream)); checkCudaErrors(cudaMemcpyAsync(d_boundary_face_vector, h_boundary_face_vector, boundary_face_vec_bytes, cudaMemcpyHostToDevice, stream)); + checkCudaErrors(cudaMemcpyAsync(d_permedIndex, permedIndex.data(), face_index_bytes, cudaMemcpyHostToDevice, stream)); + checkCudaErrors(cudaMemcpyAsync(d_bouPermedIndex, boundPermutationList.data(), boundary_face_index_bytes, cudaMemcpyHostToDevice, stream)); + checkCudaErrors(cudaMemcpyAsync(d_tmpPermutatedList, tmpPermutatedList.data(), csr_col_index_bytes, cudaMemcpyHostToDevice, stream)); + // construct AmgXSolver UxSolver = new AmgXSolver(modeStr, cfgFile); UySolver = new AmgXSolver(modeStr, cfgFile); @@ -598,43 +651,32 @@ void dfMatrix::fvm_ddt(double *rho_old, double *rho_new, double* vector_old) } void dfMatrix::fvm_div(double* phi, std::vector ueqn_internalCoeffs_init, - std::vector ueqn_boundaryCoeffs_init) + std::vector ueqn_boundaryCoeffs_init, std::vector boundary_pressure_init) { // copy and permutate face variables clock_t start = std::clock(); std::copy(phi, phi + num_surfaces, h_phi_vec_init.begin()); std::copy(phi, phi + num_surfaces, h_phi_vec_init.begin() + num_surfaces); - for (int i = 0; i < num_faces; i++) - { - h_phi_vec[i] = h_phi_vec_init[permedIndex[i]]; - } - h_phi = h_phi_vec.data(); - - // copy and permutate boundary variable - for (int i = 0; i < num_boundary_faces; i++) - { - ueqn_internalCoeffs[3*i] = ueqn_internalCoeffs_init[3*boundPermutationList[i]]; - ueqn_internalCoeffs[3*i+1] = ueqn_internalCoeffs_init[3*boundPermutationList[i]+1]; - ueqn_internalCoeffs[3*i+2] = ueqn_internalCoeffs_init[3*boundPermutationList[i]+2]; - ueqn_boundaryCoeffs[3*i] = ueqn_boundaryCoeffs_init[3*boundPermutationList[i]]; - ueqn_boundaryCoeffs[3*i+1] = ueqn_boundaryCoeffs_init[3*boundPermutationList[i]+1]; - ueqn_boundaryCoeffs[3*i+2] = ueqn_boundaryCoeffs_init[3*boundPermutationList[i]+2]; - } - h_internal_coeffs = ueqn_internalCoeffs.data(); - h_boundary_coeffs = ueqn_boundaryCoeffs.data(); clock_t end = std::clock(); time_monitor_CPU += double(end - start) / double(CLOCKS_PER_SEC); start = std::clock(); - // Copy the host input array in host memory to the device input array in device memory - checkCudaErrors(cudaMemcpyAsync(d_phi, h_phi, face_bytes, cudaMemcpyHostToDevice, stream)); - checkCudaErrors(cudaMemcpyAsync(d_internal_coeffs, h_internal_coeffs, boundary_face_vec_bytes, cudaMemcpyHostToDevice, stream)); - checkCudaErrors(cudaMemcpyAsync(d_boundary_coeffs, h_boundary_coeffs, boundary_face_vec_bytes, cudaMemcpyHostToDevice, stream)); + checkCudaErrors(cudaMemcpyAsync(d_phi_init, h_phi_vec_init.data(), face_bytes, cudaMemcpyHostToDevice, stream)); + size_t threads_per_block = 1024; + size_t blocks_per_grid = (num_faces + threads_per_block - 1) / threads_per_block; + offdiagPermutation<<>>(num_faces, d_permedIndex, d_phi_init, d_phi); + + // copy and permutate boundary variable + checkCudaErrors(cudaMemcpyAsync(d_internal_coeffs_init, ueqn_internalCoeffs_init.data(), boundary_face_vec_bytes, cudaMemcpyHostToDevice, stream)); + checkCudaErrors(cudaMemcpyAsync(d_boundary_coeffs_init, ueqn_boundaryCoeffs_init.data(), boundary_face_vec_bytes, cudaMemcpyHostToDevice, stream)); + checkCudaErrors(cudaMemcpyAsync(d_boundary_pressure_init, boundary_pressure_init.data(), boundary_face_bytes, cudaMemcpyHostToDevice, stream)); + blocks_per_grid = (num_boundary_faces + threads_per_block - 1) / threads_per_block; + boundaryPermutation<<>>(num_boundary_faces, d_bouPermedIndex, d_internal_coeffs_init, d_boundary_coeffs_init, + d_boundary_pressure_init, d_internal_coeffs, d_boundary_coeffs, d_boundary_pressure); // launch cuda kernel - size_t threads_per_block = 1024; - size_t blocks_per_grid = (num_cells + threads_per_block - 1) / threads_per_block; + blocks_per_grid = (num_cells + threads_per_block - 1) / threads_per_block; fvm_div_internal<<>>(num_cells, num_faces, d_A_csr_row_index, d_A_csr_diag_index, d_weight, d_phi, d_A_csr, d_b, d_A_csr, d_b); @@ -647,24 +689,17 @@ void dfMatrix::fvm_div(double* phi, std::vector ueqn_internalCoeffs_init time_monitor_GPU += double(end - start) / double(CLOCKS_PER_SEC); } -void dfMatrix::fvc_grad(double* pressure, std::vector boundary_pressure_init) +void dfMatrix::fvc_grad(double* pressure) { // copy cell variables directly clock_t start = std::clock(); h_pressure = pressure; - // copy and permutate boundary variable - for (int i = 0; i < num_boundary_faces; i++) - { - boundary_pressure[i] = boundary_pressure_init[boundPermutationList[i]]; - } - h_boundary_pressure = boundary_pressure.data(); clock_t end = std::clock(); time_monitor_CPU += double(end - start) / double(CLOCKS_PER_SEC); // Copy the host input array in host memory to the device input array in device memory start = std::clock(); checkCudaErrors(cudaMemcpyAsync(d_pressure, h_pressure, cell_bytes, cudaMemcpyHostToDevice, stream)); - checkCudaErrors(cudaMemcpyAsync(d_boundary_pressure, h_boundary_pressure, boundary_face_bytes, cudaMemcpyHostToDevice, stream)); // launch cuda kernel size_t threads_per_block = 1024; @@ -688,18 +723,21 @@ void dfMatrix::add_fvMatrix(double* turbSrc_low, double* turbSrc_diag, double* t std::copy(turbSrc_diag, turbSrc_diag + num_cells, h_turbSrc_init_mtx_vec.begin() + num_surfaces); std::copy(turbSrc_upp, turbSrc_upp + num_surfaces, h_turbSrc_init_mtx_vec.begin() + num_surfaces + num_cells); std::copy(turbSrc_source, turbSrc_source + 3*num_cells, h_turbSrc_init_src_vec.begin()); - // permutate - for (int i = 0; i < num_cells+2*num_surfaces; i++) - h_turbSrc_init_1mtx[i] = h_turbSrc_init_mtx_vec[tmpPermutatedList[i]]; - clock_t end = std::clock(); time_monitor_CPU += double(end - start) / double(CLOCKS_PER_SEC); - // Copy the host input array in host memory to the device input array in device memory + + // permutate start = std::clock(); - checkCudaErrors(cudaMemcpyAsync(d_turbSrc_A, h_turbSrc_init_1mtx.data(), csr_value_bytes, cudaMemcpyHostToDevice, stream)); - checkCudaErrors(cudaMemcpyAsync(d_turbSrc_b, h_turbSrc_init_src_vec.data(), cell_vec_bytes, cudaMemcpyHostToDevice, stream)); + checkCudaErrors(cudaMemcpyAsync(d_turbSrc_A_init, h_turbSrc_init_mtx_vec.data(), csr_value_bytes, cudaMemcpyHostToDevice, stream)); size_t threads_per_block = 1024; - size_t blocks_per_grid = (num_cells + threads_per_block - 1) / threads_per_block; + size_t blocks_per_grid = (num_cells + num_faces + threads_per_block - 1) / threads_per_block; + // for (int i = 0; i < num_cells+2*num_surfaces; i++) + // h_turbSrc_init_1mtx[i] = h_turbSrc_init_mtx_vec[tmpPermutatedList[i]]; + csrPermutation<<>>(num_cells + num_faces, d_tmpPermutatedList, d_turbSrc_A_init, d_turbSrc_A); + + // Copy the host input array in host memory to the device input array in device memory + checkCudaErrors(cudaMemcpyAsync(d_turbSrc_b, h_turbSrc_init_src_vec.data(), cell_vec_bytes, cudaMemcpyHostToDevice, stream)); + blocks_per_grid = (num_cells + threads_per_block - 1) / threads_per_block; add_fvMatrix_kernel<<>>(num_cells, num_faces, d_A_csr_row_index, d_turbSrc_A, d_turbSrc_b, d_A_csr, d_b, d_A_csr, d_b); end = std::clock();