diff --git a/CMakeLists.txt b/CMakeLists.txt index 34c7e8ad4b7d4d8e62fce5227d0412f0a58225fc..756eaa23b503c432275595aed7adc759cec4e4e7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -143,7 +143,8 @@ SET(NEKTAR_LIBRARY_TYPE "SHARED") # Set up RPATH SET(CMAKE_SKIP_BUILD_RPATH FALSE) -SET(CMAKE_BUILD_RPATH "${TP_LIB_DIR}") +SET(CMAKE_BUILD_RPATH "${TP_LIB_DIR}" "${TPDIST}/lib") +SET(CMAKE_INSTALL_RPATH_USE_LINK_PATH FALSE) LIST(FIND CMAKE_PLATFORM_IMPLICIT_LINK_DIRECTORIES "${CMAKE_INSTALL_PREFIX}/${LIB_DIR}" isSystemDir) IF("${isSystemDir}" STREQUAL "-1") @@ -151,6 +152,7 @@ IF("${isSystemDir}" STREQUAL "-1") "${CMAKE_INSTALL_PREFIX}/${LIB_DIR}" "${TP_LIB_DIR}/") ELSE() SET(CMAKE_INSTALL_RPATH "${TP_LIB_DIR}") + message(STATUS $CMAKE_INSTALL_RPATH) ENDIF() # Enable the use of @rpath in macOS install names so that we can use multiple @@ -259,6 +261,7 @@ INCLUDE (ThirdPartyTetGen) INCLUDE (ThirdPartyCCM) INCLUDE (ThirdPartyBlasLapack) INCLUDE (ThirdPartyCwipi) +INCLUDE (ThirdPartySaena) INCLUDE (FindCFI) INCLUDE (FindLikwid) diff --git a/cmake/ThirdPartySaena.cmake b/cmake/ThirdPartySaena.cmake new file mode 100644 index 0000000000000000000000000000000000000000..2b3ffc4ea853df2a8f0ab1aff12051456a58ce1b --- /dev/null +++ b/cmake/ThirdPartySaena.cmake @@ -0,0 +1,59 @@ +######################################################################## +# +# ThirdParty configuration for Nektar++ +# +# PETSc +# +######################################################################## + +OPTION(NEKTAR_USE_SAENA "Enable Saena parallel matrix solver support." OFF) + +# Deactivate if forced ON but NEKTAR_USE_MKL or NEKTAR_USE_MPI are OFF +CMAKE_DEPENDENT_OPTION (NEKTAR_USE_SAENA + "Enable saena parallel matrix solver support." ${NEKTAR_USE_SAENA} + "NEKTAR_USE_MKL;NEKTAR_USE_MPI" OFF) + +IF (NEKTAR_USE_SAENA) + SET(BUILD_SAENA ON) + + CMAKE_DEPENDENT_OPTION(THIRDPARTY_BUILD_SAENA + "Build Saena if needed" ${BUILD_SAENA} + "NEKTAR_USE_SAENA" OFF) + + IF (THIRDPARTY_BUILD_SAENA) + INCLUDE(ExternalProject) + + EXTERNALPROJECT_ADD( + saena + PREFIX ${TPSRC} + STAMP_DIR ${TPBUILD}/stamp + GIT_REPOSITORY https://github.com/paralab/Saena_Public.git + GIT_TAG 8f49d89347d931ef21a26979246c955d1d1de4fb + DOWNLOAD_DIR ${TPSRC} + SOURCE_DIR ${TPBUILD}/saena + TMP_DIR ${TPBUILD}/saena-tmp + INSTALL_DIR ${TPDIST} + BINARY_DIR ${TPBUILD}/saena + CONFIGURE_COMMAND ${CMAKE_COMMAND} + -G ${CMAKE_GENERATOR} + -DCMAKE_BUILD_TYPE:STRING=${CMAKE_BUILD_TYPE} + -DCMAKE_C_COMPILER:FILEPATH=${CMAKE_C_COMPILER} + -DCMAKE_CXX_COMPILER:FILEPATH=${CMAKE_CXX_COMPILER} + -DCMAKE_INSTALL_PREFIX:PATH=${TPDIST} + ${TPBUILD}/saena + ) + + THIRDPARTY_LIBRARY(SAENA_LIBRARY SHARED saena + DESCRIPTION "Saena library") + SET(SAENA_INCLUDE_DIR ${TPDIST}/include CACHE FILEPATH + "Saenae includes" FORCE) + MESSAGE(STATUS "Build Saena: ${SAENA_LIBRARY}") + SET(SAENA_CONFIG_INCLUDE_DIR ${TPINC}) + ELSE() + MESSAGE(FATAL "Saena only available through third-party install") + ENDIF() + + ADD_DEFINITIONS(-DNEKTAR_USING_SAENA) + INCLUDE_DIRECTORIES(${SAENA_INCLUDE_DIR}) + MARK_AS_ADVANCED(SAENA_LIBRARY SAENA_INCLUDE_DIR) +ENDIF() diff --git a/library/MultiRegions/AssemblyMap/AssemblyMapCG.cpp b/library/MultiRegions/AssemblyMap/AssemblyMapCG.cpp index 016fa36f6923a52a8ce0b75ac5e236f92589b4a9..e3b971672c97b74c7d608c69f671851aee3562ed 100644 --- a/library/MultiRegions/AssemblyMap/AssemblyMapCG.cpp +++ b/library/MultiRegions/AssemblyMap/AssemblyMapCG.cpp @@ -1247,6 +1247,8 @@ int AssemblyMapCG::CreateGraph( case eIterativeStaticCond: case ePETScStaticCond: case ePETScFullMatrix: + case eSaenaStaticCond: + case eSaenaFullMatrix: case eXxtFullMatrix: case eXxtStaticCond: { @@ -1263,6 +1265,7 @@ int AssemblyMapCG::CreateGraph( case ePETScMultiLevelStaticCond: case eDirectMultiLevelStaticCond: case eIterativeMultiLevelStaticCond: + case eSaenaMultiLevelStaticCond: case eXxtMultiLevelStaticCond: { MultiLevelBisectionReordering(boostGraphObj, perm, iperm, diff --git a/library/MultiRegions/CMakeLists.txt b/library/MultiRegions/CMakeLists.txt index 67daccb0e6e6b8e1954cb4b4b522f2719b20beb7..9cd109a7810d052c735c2e1df0e8a23b5df6664a 100644 --- a/library/MultiRegions/CMakeLists.txt +++ b/library/MultiRegions/CMakeLists.txt @@ -113,6 +113,18 @@ IF(NEKTAR_USE_PETSC) ) ENDIF(NEKTAR_USE_PETSC) +IF(NEKTAR_USE_SAENA) + SET(MULTI_REGIONS_HEADERS ${MULTI_REGIONS_HEADERS} + GlobalLinSysSaena.h + GlobalLinSysSaenaFull.h + GlobalLinSysSaenaStaticCond.h + ) + SET(MULTI_REGIONS_SOURCES ${MULTI_REGIONS_SOURCES} + GlobalLinSysSaena.cpp + GlobalLinSysSaenaFull.cpp + GlobalLinSysSaenaStaticCond.cpp + ) +ENDIF(NEKTAR_USE_SAENA) ADD_DEFINITIONS(-DMULTI_REGIONS_EXPORTS) @@ -133,6 +145,11 @@ IF( NEKTAR_USE_PETSC ) ADD_DEPENDENCIES(MultiRegions petsc-3.11.4) ENDIF( NEKTAR_USE_PETSC ) +IF (NEKTAR_USE_SAENA) + TARGET_LINK_LIBRARIES(MultiRegions LINK_PRIVATE ${SAENA_LIBRARY}) + ADD_DEPENDENCIES(MultiRegions saena) +ENDIF() + IF (NEKTAR_BUILD_PYTHON) SUBDIRS(Python) ENDIF() diff --git a/library/MultiRegions/GlobalLinSys.cpp b/library/MultiRegions/GlobalLinSys.cpp index 7ac5137699609c31a7452947c69c60a5feb798a3..75257be73246d4f2ecb8031f29c9bddaaeed1ed4 100644 --- a/library/MultiRegions/GlobalLinSys.cpp +++ b/library/MultiRegions/GlobalLinSys.cpp @@ -47,7 +47,7 @@ namespace Nektar { namespace MultiRegions { -std::string GlobalLinSys::lookupIds[12] = { +std::string GlobalLinSys::lookupIds[15] = { LibUtilities::SessionReader::RegisterEnumValue( "GlobalSysSoln", "DirectFull", MultiRegions::eDirectFullMatrix), LibUtilities::SessionReader::RegisterEnumValue( @@ -76,7 +76,14 @@ std::string GlobalLinSys::lookupIds[12] = { "GlobalSysSoln", "PETScStaticCond", MultiRegions::ePETScStaticCond), LibUtilities::SessionReader::RegisterEnumValue( "GlobalSysSoln", "PETScMultiLevelStaticCond", - MultiRegions::ePETScMultiLevelStaticCond)}; + MultiRegions::ePETScMultiLevelStaticCond), + LibUtilities::SessionReader::RegisterEnumValue( + "GlobalSysSoln", "SaenaFull", MultiRegions::eSaenaFullMatrix), + LibUtilities::SessionReader::RegisterEnumValue( + "GlobalSysSoln", "SaenaStaticCond", MultiRegions::eSaenaStaticCond), + LibUtilities::SessionReader::RegisterEnumValue( + "GlobalSysSoln", "SaenaMultiLevelStaticCond", + MultiRegions::eSaenaMultiLevelStaticCond)}; #ifdef NEKTAR_USE_SCOTCH std::string GlobalLinSys::def = diff --git a/library/MultiRegions/GlobalLinSysSaena.cpp b/library/MultiRegions/GlobalLinSysSaena.cpp new file mode 100644 index 0000000000000000000000000000000000000000..47f61ec26ef78e3045a8bd3230224664eedd86ab --- /dev/null +++ b/library/MultiRegions/GlobalLinSysSaena.cpp @@ -0,0 +1,268 @@ +/////////////////////////////////////////////////////////////////////////////// +// +// File GlobalLinSys.cpp +// +// For more information, please see: http://www.nektar.info +// +// The MIT License +// +// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA), +// Department of Aeronautics, Imperial College London (UK), and Scientific +// Computing and Imaging Institute, University of Utah (USA). +// +// License for the specific language governing rights and limitations under +// Permission is hereby granted, free of charge, to any person obtaining a +// copy of this software and associated documentation files (the "Software"), +// to deal in the Software without restriction, including without limitation +// the rights to use, copy, modify, merge, publish, distribute, sublicense, +// and/or sell copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included +// in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS +// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +// DEALINGS IN THE SOFTWARE. +// +// Description: GlobalLinSys definition +// +/////////////////////////////////////////////////////////////////////////////// + +#include +#include +#include + +using namespace std; + +namespace Nektar +{ +namespace MultiRegions +{ +/** + * @class GlobalLinSysSaena + * + * Solves a linear system using Saena. + */ +GlobalLinSysSaena::GlobalLinSysSaena( + const GlobalLinSysKey &pKey, const std::weak_ptr &pExp, + const std::shared_ptr &pLocToGloMap, + const unsigned int pPolynomialOrder) + : GlobalLinSys(pKey, pExp, pLocToGloMap) +{ + if (pPolynomialOrder) + { + // setup and use supplied + SetPolyOrder(pPolynomialOrder); + } +} + +/** + * @brief Clean up Saena objects. + * + * Note that if SessionReader::Finalize is called before the end of the + * program, Saena may have been finalized already, at which point we + * cannot deallocate our objects. If that's the case we do nothing and + * let the kernel clear up after us. + */ +GlobalLinSysSaena::~GlobalLinSysSaena() +{ +} + +/** + * @brief Solve linear system using Saena. + * + * The general strategy being a Saena solve is to: + * + * - Copy values into the Saena vector #m_b + * - Solve the system #m_ksp and place result into #m_x. + * - Scatter results back into #m_locVec using #m_ctx scatter object. + * - Copy from #m_locVec to output array #pOutput. + */ +void GlobalLinSysSaena::v_SolveLinearSystem( + const int pNumRows, const Array &pInput, + Array &pOutput, const AssemblyMapSharedPtr &locToGloMap, + const int pNumDir) +{ + boost::ignore_unused(locToGloMap); + + // @TODO: shouldn't need to but we require a new RHS vector every + // time this is called. + saena::vector m_rhs; + m_rhs.set_comm(m_comm); + + const int nHomDofs = pNumRows - pNumDir; + + m_rhs.set(&m_reorderedMap[0], &pInput[pNumDir], nHomDofs); + m_rhs.assemble(); + m_amg.set_rhs(m_rhs); + + // Temporary solution storage? + NekDouble *sol = nullptr; + + // Solve with pCG method + m_amg.solve_pCG(sol, &m_opts); + + Vmath::Vcopy(nHomDofs, sol, 1, &pOutput[pNumDir], 1); + + if (sol != nullptr) + { + free(sol); + sol = nullptr; + } +} + +/** + * @brief Calculate a reordering of universal IDs for Saena. + * + * Saena requires a unique, contiguous index of all global and universal + * degrees of freedom which represents its position inside the + * matrix. Presently Gs does not guarantee this, so this routine + * constructs a new universal mapping. + * + * @param glo2uniMap Global to universal map + * @param glo2unique Global to unique map + * @param pLocToGloMap Assembly map for this system + */ +void GlobalLinSysSaena::CalculateReordering( + const Array &glo2uniMap, + const Array &glo2unique, + const AssemblyMapSharedPtr &pLocToGloMap) +{ + LibUtilities::CommSharedPtr vComm = + m_expList.lock()->GetSession()->GetComm(); + + const int nDirDofs = pLocToGloMap->GetNumGlobalDirBndCoeffs(); + const int nHomDofs = glo2uniMap.size() - nDirDofs; + const int nProc = vComm->GetSize(); + const int rank = vComm->GetRank(); + + int n, cnt; + + // Count number of unique degrees of freedom on each process. + m_nLocal = Vmath::Vsum(nHomDofs, glo2unique + nDirDofs, 1); + m_reorderedMap.resize(nHomDofs); + + // Reduce coefficient counts across all processors. + Array localCounts(nProc, 0), localOffset(nProc, 0); + localCounts[rank] = nHomDofs; + vComm->AllReduce(localCounts, LibUtilities::ReduceSum); + + for (n = 1; n < nProc; ++n) + { + localOffset[n] = localOffset[n - 1] + localCounts[n - 1]; + } + + int totHomDofs = Vmath::Vsum(nProc, localCounts, 1); + vector allUniIds(totHomDofs, 0); + + // Assemble list of universal IDs + for (n = 0; n < nHomDofs; ++n) + { + int gid = n + nDirDofs; + allUniIds[n + localOffset[rank]] = glo2uniMap[gid]; + } + + // Reduce this across processors so that each process has a list of + // all universal IDs. + vComm->AllReduce(allUniIds, LibUtilities::ReduceSum); + std::sort(allUniIds.begin(), allUniIds.end()); + map uniIdReorder; + + // Renumber starting from 0. + for (cnt = n = 0; n < allUniIds.size(); ++n) + { + if (uniIdReorder.count(allUniIds[n]) > 0) + { + continue; + } + + uniIdReorder[allUniIds[n]] = cnt++; + } + + // Populate reordering map. + for (n = 0; n < nHomDofs; ++n) + { + int gid = n + nDirDofs; + int uniId = glo2uniMap[gid]; + ASSERTL0(uniIdReorder.count(uniId) > 0, "Error in ordering"); + m_reorderedMap[n] = uniIdReorder[uniId]; + } + + m_bdydof = nDirDofs; +} + +/** + * @brief Construct Saena matrix and vector handles. + * + * @todo Preallocation should be done at this point, since presently + * matrix allocation takes a significant amount of time. + * + * @param nGlobal Number of global degrees of freedom in the system (on + * this processor) + * @param nDir Number of Dirichlet degrees of freedom (on this + * processor). + */ +void GlobalLinSysSaena::SetUpMatVec() +{ + LibUtilities::CommSharedPtr comm = + m_expList.lock()->GetSession()->GetComm(); + auto mpiComm = std::dynamic_pointer_cast(comm); + + m_comm = mpiComm->GetComm(); + m_matrix.set_comm(m_comm); + m_matrix.add_duplicates(true); + m_rhs.set_comm(m_comm); +} + +/** + * @brief Set up KSP solver object. + * + * This is reasonably generic setup -- most solver types can be changed + * using the ? file. + * + * @param tolerance Residual tolerance to converge to. + */ +void GlobalLinSysSaena::SetUpSolver(NekDouble tolerance) +{ + m_scale = false; + m_opts.set_relative_tolerance(tolerance); + // m_opts.set_dynamic_levels(false); + // m_opts.set_max_lev(5); + // m_opts.set_vcycle_num(400); + // m_opts.set_smoother("chebyshev"); // chebyshev, jacobi + // m_opts.set_preSmooth(3); + // m_opts.set_postSmooth(3); +} + +void GlobalLinSysSaena::SetUpMultigrid() +{ + int nummodes = m_expList.lock()->GetFieldDefinitions()[0]->m_numModes[0]; + int p_order = m_polyOrder == 0 ? nummodes - 1 : m_polyOrder; + int prodim = m_expList.lock()->GetCoordim(0); + + m_matrix.set_p_order(p_order); + m_matrix.set_prodim(prodim); + + // set p_coarsen levels computation. subtract by a constant. + vector order_dif; + for (int i = 0; i < p_order - 1; ++i) + { + order_dif.emplace_back(1); + } + + // set number of multigrid levels + int max_h_level = 1; // h-multigrid levels + m_amg.set_multigrid_max_level(static_cast(order_dif.size()) + + max_h_level); + + m_amg.set_scale(m_scale); + m_amg.set_matrix(&m_matrix, &m_opts, m_l2g, m_reorderedMap, m_bdydof, + order_dif); +} +} // namespace MultiRegions +} // namespace Nektar diff --git a/library/MultiRegions/GlobalLinSysSaena.h b/library/MultiRegions/GlobalLinSysSaena.h new file mode 100644 index 0000000000000000000000000000000000000000..15e696fdb1181c2f52f716c8a3953fd1c47d540c --- /dev/null +++ b/library/MultiRegions/GlobalLinSysSaena.h @@ -0,0 +1,110 @@ +/////////////////////////////////////////////////////////////////////////////// +// +// File GlobalLinSys.h +// +// For more information, please see: http://www.nektar.info +// +// The MIT License +// +// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA), +// Department of Aeronautics, Imperial College London (UK), and Scientific +// Computing and Imaging Institute, University of Utah (USA). +// +// License for the specific language governing rights and limitations under +// Permission is hereby granted, free of charge, to any person obtaining a +// copy of this software and associated documentation files (the "Software"), +// to deal in the Software without restriction, including without limitation +// the rights to use, copy, modify, merge, publish, distribute, sublicense, +// and/or sell copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included +// in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS +// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +// DEALINGS IN THE SOFTWARE. +// +// Description: GlobalLinSysSaena header +// +/////////////////////////////////////////////////////////////////////////////// +#ifndef NEKTAR_LIB_MULTIREGIONS_GLOBALLINSYSSAENA_H +#define NEKTAR_LIB_MULTIREGIONS_GLOBALLINSYSSAENA_H + +#include +#include + +#include +#include + +namespace Nektar +{ +namespace MultiRegions +{ +// Forward declarations +class ExpList; + +/// A Saena global linear system. +class GlobalLinSysSaena : virtual public GlobalLinSys +{ +public: + /// Constructor for full direct matrix solve. + MULTI_REGIONS_EXPORT GlobalLinSysSaena( + const GlobalLinSysKey &pKey, const std::weak_ptr &pExp, + const std::shared_ptr &pLocToGloMap, + const unsigned int pPolynomialOrder = 0); + + MULTI_REGIONS_EXPORT virtual ~GlobalLinSysSaena(); + + virtual void v_SolveLinearSystem(const int pNumRows, + const Array &pInput, + Array &pOutput, + const AssemblyMapSharedPtr &locToGloMap, + const int pNumDir); + + void SetPolyOrder(int p) + { + m_polyOrder = p; + } + +protected: + /// Saena matrix object. + saena::matrix m_matrix; + /// Saena vector to store rhs + saena::vector m_rhs; + /// Saena object for options + saena::options m_opts; + /// Saena object that represents solver system. + saena::amg m_amg; + /// Reordering that takes universal IDs to a unique row in the Saena + /// matrix. @see GlobalLinSysSaena::CalculateReordering + std::vector m_reorderedMap; + /// MPI communicator + MPI_Comm m_comm; + /// Number of unique degrees of freedom on this process. + int m_nLocal; + /// Number of boundary degrees of freedom + int m_bdydof; + /// Mesh information + std::vector> m_l2g; + /// flag to set the linear system to be scaled + bool m_scale; + int m_polyOrder = 0; + + PreconditionerSharedPtr m_precon; + + void SetUpMatVec(); + void SetUpSolver(NekDouble tolerance); + void SetUpMultigrid(); + void CalculateReordering(const Array &glo2uniMap, + const Array &glo2unique, + const AssemblyMapSharedPtr &pLocToGloMap); +}; +} // namespace MultiRegions +} // namespace Nektar + +#endif diff --git a/library/MultiRegions/GlobalLinSysSaenaFull.cpp b/library/MultiRegions/GlobalLinSysSaenaFull.cpp new file mode 100644 index 0000000000000000000000000000000000000000..1cb8944a3f061a922554a981349c066c7baaa026 --- /dev/null +++ b/library/MultiRegions/GlobalLinSysSaenaFull.cpp @@ -0,0 +1,263 @@ +/////////////////////////////////////////////////////////////////////////////// +// +// File GlobalLinSysSaenaFull.cpp +// +// For more information, please see: http://www.nektar.info +// +// The MIT License +// +// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA), +// Department of Aeronautics, Imperial College London (UK), and Scientific +// Computing and Imaging Institute, University of Utah (USA). +// +// License for the specific language governing rights and limitations under +// Permission is hereby granted, free of charge, to any person obtaining a +// copy of this software and associated documentation files (the "Software"), +// to deal in the Software without restriction, including without limitation +// the rights to use, copy, modify, merge, publish, distribute, sublicense, +// and/or sell copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included +// in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS +// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +// DEALINGS IN THE SOFTWARE. +// +// Description: GlobalLinSysSaenaFull definition +// +/////////////////////////////////////////////////////////////////////////////// + +#include +#include + +using namespace std; + +namespace Nektar +{ +namespace MultiRegions +{ +/** + * @class GlobalLinSysSaenaFull + */ + +/** + * Registers the class with the Factory. + */ +string GlobalLinSysSaenaFull::className = + GetGlobalLinSysFactory().RegisterCreatorFunction( + "SaenaFull", GlobalLinSysSaenaFull::create, "Saena Full Matrix."); + +/// Constructor for full direct matrix solve. +GlobalLinSysSaenaFull::GlobalLinSysSaenaFull( + const GlobalLinSysKey &pLinSysKey, const std::weak_ptr &pExp, + const std::shared_ptr &pLocToGloMap, + const unsigned int pPolynomialOrder) + : GlobalLinSys(pLinSysKey, pExp, pLocToGloMap), + GlobalLinSysSaena( + pLinSysKey, pExp, pLocToGloMap, + pPolynomialOrder) // @hari - Change this constructor instead +{ + // SET UP VECTORS AND MATRIX + SetUpMatVec(); + + int rank = 0, nprocs = 0; + MPI_Comm_size(m_comm, &nprocs); + MPI_Comm_rank(m_comm, &rank); + + auto tbegin = clock(); + + const int nDirDofs = pLocToGloMap->GetNumGlobalDirBndCoeffs(); + + int i, j, n, cnt, gid1, gid2, loc_lda; + NekDouble sign1, sign2, value; + DNekScalMatSharedPtr loc_mat; + + // CALCULATE REORDERING MAPPING + CalculateReordering(pLocToGloMap->GetGlobalToUniversalMap(), + pLocToGloMap->GetGlobalToUniversalMapUnique(), + pLocToGloMap); + + // STORE MESH INFO TO BE PASSED TO SAENA + // int total_elm = this->GetExp()->size(); + auto ExpTmp = m_expList.lock()->GetExp(); + int total_elm = ExpTmp->size(); + // std::cout << total_elm << "\n"; + + int counter = 0; + vector dof_elems; + for (i = 0; i < total_elm; ++i) + { + // std::cout << ExpTmp->at(i)->GetNcoeffs() << std::endl; + for (j = 0; j < ExpTmp->at(i)->GetNcoeffs(); ++j) + { + // printf("%i\t", + // pLocToGloMap->GetLocalToGlobalMap()[counter]); + dof_elems.emplace_back( + pLocToGloMap->GetLocalToGlobalMap()[counter] + 1); + ++counter; + } + m_l2g.emplace_back(dof_elems); + dof_elems.clear(); + } + + auto tend = clock(); + auto t = double(tend - tbegin) / CLOCKS_PER_SEC; + double t_ave = 0.0; + MPI_Reduce(&t, &t_ave, 1, MPI_DOUBLE, MPI_SUM, 0, m_comm); + if (!rank) + printf("Saena mesh info generation time: %f\n", t_ave / nprocs); + + // CONSTRUCT KSP OBJECT + SetUpSolver(pLocToGloMap->GetIterativeTolerance()); + + tbegin = clock(); + + m_matrix.erase_no_shrink_to_fit(); + + // POPULATE MATRIX + for (n = cnt = 0; n < m_expList.lock()->GetNumElmts(); ++n) + { + loc_mat = GetBlock(n); + loc_lda = loc_mat->GetRows(); + + for (i = 0; i < loc_lda; ++i) + { + gid1 = pLocToGloMap->GetLocalToGlobalMap(cnt + i) - nDirDofs; + sign1 = pLocToGloMap->GetLocalToGlobalSign(cnt + i); + if (gid1 >= 0) + { + int gid1ro = m_reorderedMap[gid1]; + for (j = 0; j < loc_lda; ++j) + { + gid2 = + pLocToGloMap->GetLocalToGlobalMap(cnt + j) - nDirDofs; + sign2 = pLocToGloMap->GetLocalToGlobalSign(cnt + j); + if (gid2 >= 0) + { + int gid2ro = m_reorderedMap[gid2]; + value = sign1 * sign2 * (*loc_mat)(i, j); + m_matrix.set(gid1ro, gid2ro, value); + } + } + } + } + cnt += loc_lda; + } + + // timing + tend = clock(); + t = double(tend - tbegin) / CLOCKS_PER_SEC; + MPI_Reduce(&t, &t_ave, 1, MPI_DOUBLE, MPI_SUM, 0, m_comm); + if (!rank) + printf("nektar assembly time: %f\n", t_ave / nprocs); + + tbegin = clock(); + + // ASSEMBLE MATRIX + // m_matrix.set_num_threads(1); + m_matrix.assemble(m_scale); + // m_matrix.assemble_writeToFile("matrix_folder"); + + // timing + tend = clock(); + t = double(tend - tbegin) / CLOCKS_PER_SEC; + MPI_Reduce(&t, &t_ave, 1, MPI_DOUBLE, MPI_SUM, 0, m_comm); + if (!rank) + printf("Saena matrix assembly time: %f\n", t_ave / nprocs); + + SetUpMultigrid(); +} + +GlobalLinSysSaenaFull::~GlobalLinSysSaenaFull() +{ +} + +/** + * Solve the linear system using a full global matrix system. + */ +void GlobalLinSysSaenaFull::v_Solve( + const Array &pLocInput, + Array &pLocOutput, + const AssemblyMapSharedPtr &pLocToGloMap, + const Array &pDirForcing) +{ + std::shared_ptr expList = m_expList.lock(); + bool dirForcCalculated = (bool)pDirForcing.size(); + int nDirDofs = pLocToGloMap->GetNumGlobalDirBndCoeffs(); + int nGlobDofs = pLocToGloMap->GetNumGlobalCoeffs(); + int nLocDofs = pLocToGloMap->GetNumLocalCoeffs(); + + // m_locToGloMap = pLocToGloMap; // required for DoMatrixMultiply + + Array tmp(nLocDofs); + Array tmp1(nLocDofs); + Array global(nGlobDofs, 0.0); + + int nDirTotal = nDirDofs; + expList->GetComm()->GetRowComm()->AllReduce(nDirTotal, + LibUtilities::ReduceSum); + + if (nDirTotal) + { + // calculate the dirichlet forcing + if (dirForcCalculated) + { + // assume pDirForcing is in local space + ASSERTL0( + pDirForcing.size() >= nLocDofs, + "DirForcing is not of sufficient size. Is it in local space?"); + Vmath::Vsub(nLocDofs, pLocInput, 1, pDirForcing, 1, tmp1, 1); + } + else + { + // Calculate the dirichlet forcing and substract it + // from the rhs + expList->GeneralMatrixOp(m_linSysKey, pLocOutput, tmp); + + // Apply robin boundary conditions to the solution. + for (auto &r : m_robinBCInfo) // add robin mass matrix + { + RobinBCInfoSharedPtr rBC; + Array tmploc; + + int n = r.first; + + int offset = expList->GetCoeff_Offset(n); + LocalRegions::ExpansionSharedPtr vExp = expList->GetExp(n); + + // add local matrix contribution + for (rBC = r.second; rBC; rBC = rBC->next) + { + vExp->AddRobinTraceContribution( + rBC->m_robinID, rBC->m_robinPrimitiveCoeffs, + pLocOutput + offset, tmploc = tmp + offset); + } + } + + Vmath::Vsub(nLocDofs, pLocInput, 1, tmp, 1, tmp1, 1); + } + + pLocToGloMap->Assemble(tmp1, tmp); + + SolveLinearSystem(nGlobDofs, tmp, global, pLocToGloMap, nDirDofs); + + pLocToGloMap->GlobalToLocal(global, tmp); + + // Add back initial and boundary condition + Vmath::Vadd(nLocDofs, tmp, 1, pLocOutput, 1, pLocOutput, 1); + } + else + { + pLocToGloMap->Assemble(pLocInput, tmp); + SolveLinearSystem(nGlobDofs, tmp, global, pLocToGloMap); + pLocToGloMap->GlobalToLocal(global, pLocOutput); + } +} +} // namespace MultiRegions +} // namespace Nektar diff --git a/library/MultiRegions/GlobalLinSysSaenaFull.h b/library/MultiRegions/GlobalLinSysSaenaFull.h new file mode 100644 index 0000000000000000000000000000000000000000..84364d85a7b0062915b63f81dedfeafa282f0c21 --- /dev/null +++ b/library/MultiRegions/GlobalLinSysSaenaFull.h @@ -0,0 +1,90 @@ +/////////////////////////////////////////////////////////////////////////////// +// +// File GlobalLinSysDirectXxt.h +// +// For more information, please see: http://www.nektar.info +// +// The MIT License +// +// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA), +// Department of Aeronautics, Imperial College London (UK), and Scientific +// Computing and Imaging Institute, University of Utah (USA). +// +// License for the specific language governing rights and limitations under +// Permission is hereby granted, free of charge, to any person obtaining a +// copy of this software and associated documentation files (the "Software"), +// to deal in the Software without restriction, including without limitation +// the rights to use, copy, modify, merge, publish, distribute, sublicense, +// and/or sell copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included +// in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS +// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +// DEALINGS IN THE SOFTWARE. +// +// Description: GlobalLinSysDirectXxt header +// +/////////////////////////////////////////////////////////////////////////////// + +#ifndef NEKTAR_LIB_MULTIREGIONS_GLOBALLINSYSSAENAFULL_H +#define NEKTAR_LIB_MULTIREGIONS_GLOBALLINSYSSAENAFULL_H + +#include +#include +#include + +namespace Nektar +{ +namespace MultiRegions +{ +// Forward declarations + +// class AssemblyMapDG; +class ExpList; + +/// A global linear system. +class GlobalLinSysSaenaFull : public GlobalLinSysSaena +{ +public: + /// Creates an instance of this class + static GlobalLinSysSharedPtr create( + const GlobalLinSysKey &pLinSysKey, + const std::weak_ptr &pExpList, + const std::shared_ptr &pLocToGloMap) + { + return MemoryManager::AllocateSharedPtr( + pLinSysKey, pExpList, + pLocToGloMap); // check for def args for poly order + } + + /// Name of class + MULTI_REGIONS_EXPORT static std::string className; + + /// Constructor for full direct matrix solve. + MULTI_REGIONS_EXPORT GlobalLinSysSaenaFull( + const GlobalLinSysKey &pLinSysKey, + const std::weak_ptr &pExpList, + const std::shared_ptr &pLocToGloMap, + const unsigned int pPolynomialOrder = 0); + + MULTI_REGIONS_EXPORT virtual ~GlobalLinSysSaenaFull(); + +private: + /// Solve the linear system for given input and output vectors + /// using a specified local to global map. + virtual void v_Solve( + const Array &in, Array &out, + const AssemblyMapSharedPtr &locToGloMap, + const Array &dirForcing = NullNekDouble1DArray); +}; +} // namespace MultiRegions +} // namespace Nektar + +#endif diff --git a/library/MultiRegions/GlobalLinSysSaenaStaticCond.cpp b/library/MultiRegions/GlobalLinSysSaenaStaticCond.cpp new file mode 100644 index 0000000000000000000000000000000000000000..c30a7810d9db5bf1d86bc2b5ea18e76bfa04b37d --- /dev/null +++ b/library/MultiRegions/GlobalLinSysSaenaStaticCond.cpp @@ -0,0 +1,223 @@ +/////////////////////////////////////////////////////////////////////////////// +// +// File GlobalLinSysSaenaStaticCond.cpp +// +// For more information, please see: http://www.nektar.info +// +// The MIT License +// +// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA), +// Department of Aeronautics, Imperial College London (UK), and Scientific +// Computing and Imaging Institute, University of Utah (USA). +// +// License for the specific language governing rights and limitations under +// Permission is hereby granted, free of charge, to any person obtaining a +// copy of this software and associated documentation files (the "Software"), +// to deal in the Software without restriction, including without limitation +// the rights to use, copy, modify, merge, publish, distribute, sublicense, +// and/or sell copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included +// in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS +// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +// DEALINGS IN THE SOFTWARE. +// +// Description: GlobalLinSysSaenaStaticCond definition +// +/////////////////////////////////////////////////////////////////////////////// + +#include + +//#include +//#include +//#include + +using namespace std; + +namespace Nektar +{ +namespace MultiRegions +{ +/** + * @class GlobalLinSysSaena + * + * Solves a linear system using single- or multi-level static + * condensation. + */ + +/** + * Registers the class with the Factory. + */ +string GlobalLinSysSaenaStaticCond::className = + GetGlobalLinSysFactory().RegisterCreatorFunction( + "SaenaStaticCond", GlobalLinSysSaenaStaticCond::create, + "Saena static condensation."); + +string GlobalLinSysSaenaStaticCond::className2 = + GetGlobalLinSysFactory().RegisterCreatorFunction( + "SaenaMultiLevelStaticCond", GlobalLinSysSaenaStaticCond::create, + "Saena multi-level static condensation."); + +/** + * For a matrix system of the form @f[ + * \left[ \begin{array}{cc} + * \boldsymbol{A} & \boldsymbol{B}\\ + * \boldsymbol{C} & \boldsymbol{D} + * \end{array} \right] + * \left[ \begin{array}{c} \boldsymbol{x_1}\\ \boldsymbol{x_2} + * \end{array}\right] + * = \left[ \begin{array}{c} \boldsymbol{y_1}\\ \boldsymbol{y_2} + * \end{array}\right], + * @f] + * where @f$\boldsymbol{D}@f$ and + * @f$(\boldsymbol{A-BD^{-1}C})@f$ are invertible, store and assemble + * a static condensation system, according to a given local to global + * mapping. #m_linSys is constructed by AssembleSchurComplement(). + * @param mKey Associated matrix key. + * @param pLocMatSys LocalMatrixSystem + * @param locToGloMap Local to global mapping. + */ +GlobalLinSysSaenaStaticCond::GlobalLinSysSaenaStaticCond( + const GlobalLinSysKey &pKey, const std::weak_ptr &pExpList, + const std::shared_ptr &pLocToGloMap) + : GlobalLinSys(pKey, pExpList, pLocToGloMap), + GlobalLinSysSaena(pKey, pExpList, pLocToGloMap), + GlobalLinSysStaticCond(pKey, pExpList, pLocToGloMap) +{ + std::cout << __func__ << std::endl; + + ASSERTL1((pKey.GetGlobalSysSolnType() == eSaenaStaticCond) || + (pKey.GetGlobalSysSolnType() == eSaenaMultiLevelStaticCond), + "This constructor is only valid when using static " + "condensation"); + ASSERTL1(pKey.GetGlobalSysSolnType() == + pLocToGloMap->GetGlobalSysSolnType(), + "The local to global map is not set up for the requested " + "solution type"); +} + +/** + * + */ +GlobalLinSysSaenaStaticCond::GlobalLinSysSaenaStaticCond( + const GlobalLinSysKey &pKey, const std::weak_ptr &pExpList, + const DNekScalBlkMatSharedPtr pSchurCompl, + const DNekScalBlkMatSharedPtr pBinvD, const DNekScalBlkMatSharedPtr pC, + const DNekScalBlkMatSharedPtr pInvD, + const std::shared_ptr &pLocToGloMap, + const PreconditionerSharedPtr pPrecon) + : GlobalLinSys(pKey, pExpList, pLocToGloMap), + GlobalLinSysSaena(pKey, pExpList, pLocToGloMap), + GlobalLinSysStaticCond(pKey, pExpList, pLocToGloMap) +{ + std::cout << __func__ << std::endl; + + m_schurCompl = pSchurCompl; + m_BinvD = pBinvD; + m_C = pC; + m_invD = pInvD; + m_precon = pPrecon; +} + +/** + * + */ +GlobalLinSysSaenaStaticCond::~GlobalLinSysSaenaStaticCond() +{ +} + +/** + * Assemble the schur complement matrix from the block matrices stored + * in #m_blkMatrices and the given local to global mapping information. + * @param locToGloMap Local to global mapping information. + */ +void GlobalLinSysSaenaStaticCond::v_AssembleSchurComplement( + AssemblyMapSharedPtr pLocToGloMap) +{ + std::cout << __func__ << std::endl; + + int i, j, n, cnt, gid1, gid2, loc_lda; + NekDouble sign1, sign2, value; + + const int nDirDofs = pLocToGloMap->GetNumGlobalDirBndCoeffs(); + + DNekScalBlkMatSharedPtr SchurCompl = m_schurCompl; + DNekScalBlkMatSharedPtr BinvD = m_BinvD; + DNekScalBlkMatSharedPtr C = m_C; + DNekScalBlkMatSharedPtr invD = m_invD; + DNekScalMatSharedPtr loc_mat; + + // CALCULATE REORDERING MAPPING + CalculateReordering(pLocToGloMap->GetGlobalToUniversalBndMap(), + pLocToGloMap->GetGlobalToUniversalBndMapUnique(), + pLocToGloMap); + + // SET UP VECTORS AND MATRIX + // SetUpMatVec(pLocToGloMap->GetNumGlobalBndCoeffs(), nDirDofs); + SetUpMatVec(); + + // CONSTRUCT KSP OBJECT + SetUpSolver(pLocToGloMap->GetIterativeTolerance()); + + // POPULATE MATRIX + for (n = cnt = 0; n < m_schurCompl->GetNumberOfBlockRows(); ++n) + { + loc_mat = m_schurCompl->GetBlock(n, n); + loc_lda = loc_mat->GetRows(); + + for (i = 0; i < loc_lda; ++i) + { + gid1 = pLocToGloMap->GetLocalToGlobalBndMap(cnt + i) - nDirDofs; + sign1 = pLocToGloMap->GetLocalToGlobalBndSign(cnt + i); + if (gid1 >= 0) + { + int gid1ro = m_reorderedMap[gid1]; + for (j = 0; j < loc_lda; ++j) + { + gid2 = pLocToGloMap->GetLocalToGlobalBndMap(cnt + j) - + nDirDofs; + sign2 = pLocToGloMap->GetLocalToGlobalBndSign(cnt + j); + if (gid2 >= 0) + { + int gid2ro = m_reorderedMap[gid2]; + value = sign1 * sign2 * (*loc_mat)(i, j); + m_matrix.set(gid1ro, gid2ro, value); + } + } + } + } + cnt += loc_lda; + } + + m_matrix.assemble(); +} + +GlobalLinSysStaticCondSharedPtr GlobalLinSysSaenaStaticCond::v_Recurse( + const GlobalLinSysKey &mkey, const std::weak_ptr &pExpList, + const DNekScalBlkMatSharedPtr pSchurCompl, + const DNekScalBlkMatSharedPtr pBinvD, const DNekScalBlkMatSharedPtr pC, + const DNekScalBlkMatSharedPtr pInvD, + const std::shared_ptr &l2gMap) +{ + // GlobalLinSysSaenaStaticCondSharedPtr sys = MemoryManager< + // GlobalLinSysSaenaStaticCond>::AllocateSharedPtr( + // mkey, pExpList, pSchurCompl, pBinvD, pC, pInvD, + // l2gMap, m_precon); + GlobalLinSysSaenaStaticCondSharedPtr sys = + MemoryManager::AllocateSharedPtr( + mkey, pExpList, pSchurCompl, pBinvD, pC, pInvD, l2gMap); + + std::cout << __func__ << std::endl; + + sys->Initialise(l2gMap); + return sys; +} +} // namespace MultiRegions +} // namespace Nektar diff --git a/library/MultiRegions/GlobalLinSysSaenaStaticCond.h b/library/MultiRegions/GlobalLinSysSaenaStaticCond.h new file mode 100644 index 0000000000000000000000000000000000000000..e550fc4f8bb4091fb3c71989ab138ad2b7de0ed4 --- /dev/null +++ b/library/MultiRegions/GlobalLinSysSaenaStaticCond.h @@ -0,0 +1,107 @@ +/////////////////////////////////////////////////////////////////////////////// +// +// File GlobalLinSysSaenaStaticCond.h +// +// For more information, please see: http://www.nektar.info +// +// The MIT License +// +// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA), +// Department of Aeronautics, Imperial College London (UK), and Scientific +// Computing and Imaging Institute, University of Utah (USA). +// +// License for the specific language governing rights and limitations under +// Permission is hereby granted, free of charge, to any person obtaining a +// copy of this software and associated documentation files (the "Software"), +// to deal in the Software without restriction, including without limitation +// the rights to use, copy, modify, merge, publish, distribute, sublicense, +// and/or sell copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included +// in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS +// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +// DEALINGS IN THE SOFTWARE. +// +// Description: GlobalLinSysStaticCond header +// +/////////////////////////////////////////////////////////////////////////////// + +#ifndef NEKTAR_LIB_MULTIREGIONS_GLOBALLINSYSSAENASTATICCOND_H +#define NEKTAR_LIB_MULTIREGIONS_GLOBALLINSYSSAENASTATICCOND_H + +#include +#include +#include + +namespace Nektar +{ +namespace MultiRegions +{ +// Forward declarations +class ExpList; +class GlobalLinSysSaenaStaticCond; + +typedef std::shared_ptr + GlobalLinSysSaenaStaticCondSharedPtr; + +/// A global linear system. +class GlobalLinSysSaenaStaticCond : virtual public GlobalLinSysSaena, + virtual public GlobalLinSysStaticCond +{ +public: + /// Creates an instance of this class + static GlobalLinSysSharedPtr create( + const GlobalLinSysKey &pLinSysKey, + const std::weak_ptr &pExpList, + const std::shared_ptr &pLocToGloMap) + { + GlobalLinSysSharedPtr p = + MemoryManager::AllocateSharedPtr( + pLinSysKey, pExpList, pLocToGloMap); + p->InitObject(); + return p; + } + + /// Name of class + MULTI_REGIONS_EXPORT static std::string className; + static std::string className2; + + /// Constructor for full direct matrix solve. + MULTI_REGIONS_EXPORT GlobalLinSysSaenaStaticCond( + const GlobalLinSysKey &mkey, const std::weak_ptr &pExpList, + const std::shared_ptr &locToGloMap); + + /// Constructor for full direct matrix solve. + MULTI_REGIONS_EXPORT GlobalLinSysSaenaStaticCond( + const GlobalLinSysKey &mkey, const std::weak_ptr &pExpList, + const DNekScalBlkMatSharedPtr pSchurCompl, + const DNekScalBlkMatSharedPtr pBinvD, const DNekScalBlkMatSharedPtr pC, + const DNekScalBlkMatSharedPtr pInvD, + const std::shared_ptr &locToGloMap, + const PreconditionerSharedPtr pPrecon = PreconditionerSharedPtr()); + + MULTI_REGIONS_EXPORT virtual ~GlobalLinSysSaenaStaticCond(); + +protected: + /// Assemble the Schur complement matrix. + virtual void v_AssembleSchurComplement( + std::shared_ptr locToGloMap); + + virtual GlobalLinSysStaticCondSharedPtr v_Recurse( + const GlobalLinSysKey &mkey, const std::weak_ptr &pExpList, + const DNekScalBlkMatSharedPtr pSchurCompl, + const DNekScalBlkMatSharedPtr pBinvD, const DNekScalBlkMatSharedPtr pC, + const DNekScalBlkMatSharedPtr pInvD, + const std::shared_ptr &locToGloMap); +}; +} // namespace MultiRegions +} // namespace Nektar + +#endif diff --git a/library/MultiRegions/MultiRegions.hpp b/library/MultiRegions/MultiRegions.hpp index 6680bc685ed938ea7e0451d53b667ad6eea41fc3..b69817afca9999aaa6a8dc09ecbed1f54b567d6a 100644 --- a/library/MultiRegions/MultiRegions.hpp +++ b/library/MultiRegions/MultiRegions.hpp @@ -176,41 +176,195 @@ struct PeriodicEntity { } - PeriodicEntity() - { - } + // Orientation of adjacent edge for use with boundary + // constraints + enum AdjacentTraceOrientation + { + eAdjacentEdgeIsForwards, + eAdjacentEdgeIsBackwards + }; - /// Geometry ID of entity. - int id; - /// Orientation of entity within higher dimensional entity. - StdRegions::Orientation orient; - /// Flag specifying if this entity is local to this partition. - bool isLocal; -}; + // Orientation of adjacent face for use with boundary + // constraints + enum AdjacentFaceOrientation + { + eAdjacentFaceDir1FwdDir1_Dir2FwdDir2, + eAdjacentFaceDir1FwdDir1_Dir2BwdDir2, + eAdjacentFaceDir1BwdDir1_Dir2FwdDir2, + eAdjacentFaceDir1BwdDir1_Dir2BwdDir2, + eAdjacentFaceDir1FwdDir2_Dir2FwdDir1, + eAdjacentFaceDir1FwdDir2_Dir2BwdDir1, + eAdjacentFaceDir1BwdDir2_Dir2FwdDir1, + eAdjacentFaceDir1BwdDir2_Dir2BwdDir1 + }; -typedef std::map> PeriodicMap; -static PeriodicMap NullPeriodicMap; + enum GlobalSysSolnType + { + eNoSolnType, ///< No Solution type specified + eDirectFullMatrix, + eDirectStaticCond, + eDirectMultiLevelStaticCond, + eIterativeFull, + eIterativeStaticCond, + eIterativeMultiLevelStaticCond, + eXxtFullMatrix, + eXxtStaticCond, + eXxtMultiLevelStaticCond, + ePETScFullMatrix, + ePETScStaticCond, + ePETScMultiLevelStaticCond, + eSaenaFullMatrix, + eSaenaStaticCond, + eSaenaMultiLevelStaticCond, + eSIZE_GlobalSysSolnType + }; -struct RotPeriodicInfo -{ - RotPeriodicInfo(const int dir, const NekDouble angle, const NekDouble tol) - : m_dir(dir), m_angle(angle), m_tol(tol) - { - } - RotPeriodicInfo() - { - } + const char* const GlobalSysSolnTypeMap[] = + { + "No Solution Type", + "DirectFull", + "DirectStaticCond", + "DirectMultiLevelStaticCond", + "IterativeFull", + "IterativeStaticCond", + "IterativeMultiLevelStaticCond", + "XxtFull", + "XxtStaticCond", + "XxtMultiLevelStaticCond", + "PETScFull", + "PETScStaticCond", + "PETScMultiLevelStaticCond", + "SaenaFull", + "SaenaStaticCond", + "SaenaMultiLevelStaticCond" + }; - /// Axis of rotation. 0 = 'x', 1 = 'y', 2 = 'z' - int m_dir; - /// Angle of rotation in radians - NekDouble m_angle; - /// Tolerance to rotation is considered identical - NekDouble m_tol; -}; + /// Type of Galerkin projection. + enum ProjectionType + { + eGalerkin, + eDiscontinuous, + eMixed_CG_Discontinuous + }; + + enum PreconditionerType + { + eNull, ///< No Solution type specified + eDiagonal, + eLinearWithDiagonal, + eLinear, + eLowEnergy, + eLinearWithLowEnergy, + eBlock, + eLinearWithBlock + }; + + const char* const PreconditionerTypeMap[] = + { + "Null", + "Diagonal", + "FullLinearSpaceWithDiagonal", + "FullLinearSpace", + "LowEnergyBlock", + "FullLinearSpaceWithLowEnergyBlock", + "Block", + "FullLinearSpaceWithBlock" + }; + + + enum LinSysIterSolver + { + eNoLinSysIterSolver,///< No Solution type specified + eConjugateGradient, + eGMRES + }; + + const char* const LinSysIterSolverMap[] = + { + "NoLinSysIterSolver", + "ConjugateGradient", + "GMRES" + }; + + + // let's keep this for linking to external + // sparse libraries + enum MatrixStorageType + { + eSmvBSR + }; + + const char* const MatrixStorageTypeMap[] = + { + "SmvBSR" + }; + + + typedef std::vector BndTypesVector; + typedef std::vector::iterator BndTypesVectorIter; + + + // structure to hold information about robin boundary conditions + + struct RobinBCInfo + { + RobinBCInfo(const int id, const Array &primCoeffs): + m_robinID(id), + m_robinPrimitiveCoeffs(primCoeffs) + { + } + + virtual ~RobinBCInfo() + {}; + + int m_robinID; /// id of which edge/face is robin condition + Array< OneD, const NekDouble > m_robinPrimitiveCoeffs; + std::shared_ptr next; + }; + + typedef std::shared_ptr RobinBCInfoSharedPtr; + + struct PeriodicEntity + { + PeriodicEntity( + const int id, + const StdRegions::Orientation orient, + const bool isLocal) : + id(id), orient(orient), isLocal(isLocal) {} + + PeriodicEntity() {} + + /// Geometry ID of entity. + int id; + /// Orientation of entity within higher dimensional entity. + StdRegions::Orientation orient; + /// Flag specifying if this entity is local to this partition. + bool isLocal; + }; + + typedef std::map > PeriodicMap; + static PeriodicMap NullPeriodicMap; + + struct RotPeriodicInfo + { + RotPeriodicInfo( + const int dir, + const NekDouble angle, + const NekDouble tol) : + m_dir(dir), m_angle(angle), m_tol(tol) {} + + RotPeriodicInfo() {} + + /// Axis of rotation. 0 = 'x', 1 = 'y', 2 = 'z' + int m_dir; + /// Angle of rotation in radians + NekDouble m_angle; + /// Tolerance to rotation is considered identical + NekDouble m_tol; + }; -} // namespace MultiRegions -} // namespace Nektar + }// end of namespace +}// end of namespace #endif diff --git a/library/MultiRegions/PreconditionerLinear.cpp b/library/MultiRegions/PreconditionerLinear.cpp index bbbc7ff3b9adb19929838ee9482ed244fa4ba42c..60302ab516d5343ede2b0ee1c882f3d825d48a56 100644 --- a/library/MultiRegions/PreconditionerLinear.cpp +++ b/library/MultiRegions/PreconditionerLinear.cpp @@ -43,6 +43,10 @@ #include #endif +#ifdef NEKTAR_USING_SAENA +#include +#endif + #include #include @@ -67,6 +71,8 @@ std::string PreconditionerLinear::solveType = std::string PreconditionerLinear::solveTypeIds[] = { LibUtilities::SessionReader::RegisterEnumValue( "LinearPreconSolver", "PETSc", MultiRegions::eLinearPreconPETSc), + LibUtilities::SessionReader::RegisterEnumValue( + "LinearPreconSolver", "Saena", MultiRegions::eLinearPreconSaena), LibUtilities::SessionReader::RegisterEnumValue( "LinearPreconSolver", "Xxt", MultiRegions::eLinearPreconXxt)}; @@ -112,6 +118,15 @@ void PreconditionerLinear::v_BuildPreconditioner() #ifndef NEKTAR_USING_PETSC NEKERROR(ErrorUtil::efatal, "Nektar++ has not been compiled with " "PETSc support."); +#endif + break; + } + case eLinearPreconSaena: + { + linSolveType = eSaenaFullMatrix; +#ifndef NEKTAR_USING_SAENA + NEKERROR(ErrorUtil::efatal, "Nektar++ has not been compiled with " + "Saena support."); #endif break; } @@ -154,6 +169,21 @@ void PreconditionerLinear::v_BuildPreconditioner() ASSERTL0(false, "Nektar++ has not been compiled with " "PETSc support."); #endif + break; + } + case eLinearPreconSaena: + { +#ifdef NEKTAR_USING_SAENA + auto vertLinsys = + MemoryManager::AllocateSharedPtr( + preconKey, expList, m_vertLocToGloMap, 1); + // vertLinsys->SetPolyOrder(1); + m_vertLinsys = vertLinsys; +#else + ASSERTL0(false, "Nektar++ has not been compiled with " + "Saena support."); +#endif + break; } } } diff --git a/library/MultiRegions/PreconditionerLinear.h b/library/MultiRegions/PreconditionerLinear.h index 890c0bf37ae2f140a5595ad2b230db232ace6778..224f2f32899e3fbfef16f5616ed82c94f97d6877 100644 --- a/library/MultiRegions/PreconditionerLinear.h +++ b/library/MultiRegions/PreconditionerLinear.h @@ -47,7 +47,8 @@ namespace MultiRegions enum LinearPreconSolver { eLinearPreconXxt, - eLinearPreconPETSc + eLinearPreconPETSc, + eLinearPreconSaena }; class PreconditionerLinear;