From a07166255e9f0a3d68bfff87f442a00492480cb4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Benjamin=20Menk=C3=BCc?= Date: Sat, 8 Jan 2022 18:18:59 +0100 Subject: [PATCH 1/7] add wrapper to interpolate for python api, because pybind11 cannot handle std::vector --- src/expression/expression.cpp | 9 +++++++++ src/expression/expression.h | 1 + 2 files changed, 10 insertions(+) diff --git a/src/expression/expression.cpp b/src/expression/expression.cpp index 3285c663..a2844300 100755 --- a/src/expression/expression.cpp +++ b/src/expression/expression.cpp @@ -556,6 +556,15 @@ std::vector expression::max(int physreg, expression* meshdeform, int ref return maxval; } +void expression::interpolate(int physreg, std::vector& xyzcoord, std::vector& interpolated, std::vector& isfound) +{ + std::vector isfoundBool; + interpolate(physreg, xyzcoord, interpolated, isfoundBool); + isfound = std::vector(isfoundBool.size()); + for (size_t i = 0; i < isfound.size(); ++i) + isfound[i] = isfoundBool[i]; +} + void expression::interpolate(int physreg, std::vector& xyzcoord, std::vector& interpolated, std::vector& isfound) { universe::getrawmesh()->getphysicalregions()->errorundefined({physreg}); diff --git a/src/expression/expression.h b/src/expression/expression.h index 822f1f6f..6fc8c430 100755 --- a/src/expression/expression.h +++ b/src/expression/expression.h @@ -122,6 +122,7 @@ class expression // Only the highest dimension elements in physical region 'physreg' are considered. // In case the ith coordinate is not in the physical region or there was any other // issue then 'isfound[i]' is false. + void interpolate(int physreg, std::vector& xyzcoord, std::vector& interpolated, std::vector& isfound); void interpolate(int physreg, std::vector& xyzcoord, std::vector& interpolated, std::vector& isfound); void interpolate(int physreg, expression meshdeform, std::vector& xyzcoord, std::vector& interpolated, std::vector& isfound); // These two functions are added for convenience and work only to interpolate at a single (x,y,z) coordinate. From 77cc66553d4cc5bf9aed8818a2f4978d4c1b626c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Benjamin=20Menk=C3=BCc?= Date: Mon, 10 Jan 2022 23:24:23 +0100 Subject: [PATCH 2/7] undo last commit, its not needed anymore --- src/expression/expression.cpp | 9 --------- src/expression/expression.h | 1 - 2 files changed, 10 deletions(-) diff --git a/src/expression/expression.cpp b/src/expression/expression.cpp index a2844300..3285c663 100755 --- a/src/expression/expression.cpp +++ b/src/expression/expression.cpp @@ -556,15 +556,6 @@ std::vector expression::max(int physreg, expression* meshdeform, int ref return maxval; } -void expression::interpolate(int physreg, std::vector& xyzcoord, std::vector& interpolated, std::vector& isfound) -{ - std::vector isfoundBool; - interpolate(physreg, xyzcoord, interpolated, isfoundBool); - isfound = std::vector(isfoundBool.size()); - for (size_t i = 0; i < isfound.size(); ++i) - isfound[i] = isfoundBool[i]; -} - void expression::interpolate(int physreg, std::vector& xyzcoord, std::vector& interpolated, std::vector& isfound) { universe::getrawmesh()->getphysicalregions()->errorundefined({physreg}); diff --git a/src/expression/expression.h b/src/expression/expression.h index 6fc8c430..822f1f6f 100755 --- a/src/expression/expression.h +++ b/src/expression/expression.h @@ -122,7 +122,6 @@ class expression // Only the highest dimension elements in physical region 'physreg' are considered. // In case the ith coordinate is not in the physical region or there was any other // issue then 'isfound[i]' is false. - void interpolate(int physreg, std::vector& xyzcoord, std::vector& interpolated, std::vector& isfound); void interpolate(int physreg, std::vector& xyzcoord, std::vector& interpolated, std::vector& isfound); void interpolate(int physreg, expression meshdeform, std::vector& xyzcoord, std::vector& interpolated, std::vector& isfound); // These two functions are added for convenience and work only to interpolate at a single (x,y,z) coordinate. From 6582b455ba93b17f8cacbcf7a1eaa57200d3d77b Mon Sep 17 00:00:00 2001 From: Benjamin Menkuec Date: Wed, 19 Jan 2022 14:34:12 +0100 Subject: [PATCH 3/7] update required cmake version because of FILE_LINK --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index aeeed1d0..e7d0d7da 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,6 +1,6 @@ # Thanks to S. Matsievskiy for bringing cmake to the project. -cmake_minimum_required(VERSION 3.0 FATAL_ERROR) +cmake_minimum_required(VERSION 3.14 FATAL_ERROR) project(Sparselizard LANGUAGES CXX) From f1949cf46e779bf5d8ad21b334fb8cec388a91e0 Mon Sep 17 00:00:00 2001 From: benjamin menkuec Date: Sun, 13 Mar 2022 21:33:38 +0100 Subject: [PATCH 4/7] added getavals() and getdvals() --- src/formulation/mat.cpp | 2 ++ src/formulation/mat.h | 2 ++ src/formulation/rawmat.cpp | 10 ++++++++++ src/formulation/rawmat.h | 2 ++ 4 files changed, 16 insertions(+) diff --git a/src/formulation/mat.cpp b/src/formulation/mat.cpp index e52ea8d2..4e4a0d7a 100755 --- a/src/formulation/mat.cpp +++ b/src/formulation/mat.cpp @@ -102,6 +102,8 @@ vec mat::eliminate(vec b) indexmat mat::getainds(void) { errorifpointerisnull(); errorifinvalidated(); return rawmatptr->getainds(); } indexmat mat::getdinds(void) { errorifpointerisnull(); errorifinvalidated(); return rawmatptr->getdinds(); } +densemat mat::getavals(void) { errorifpointerisnull(); errorifinvalidated(); return rawmatptr->getavals(); } +densemat mat::getdvals(void) { errorifpointerisnull(); errorifinvalidated(); return rawmatptr->getdvals(); } Mat mat::getapetsc(void) { errorifpointerisnull(); errorifinvalidated(); return rawmatptr->getapetsc(); } Mat mat::getdpetsc(void) { errorifpointerisnull(); errorifinvalidated(); return rawmatptr->getdpetsc(); } diff --git a/src/formulation/mat.h b/src/formulation/mat.h index 16211d23..6d961334 100755 --- a/src/formulation/mat.h +++ b/src/formulation/mat.h @@ -68,6 +68,8 @@ class mat indexmat getainds(void); indexmat getdinds(void); + densemat getavals(void); + densemat getdvals(void); Mat getapetsc(void); Mat getdpetsc(void); diff --git a/src/formulation/rawmat.cpp b/src/formulation/rawmat.cpp index 40688572..6842159d 100755 --- a/src/formulation/rawmat.cpp +++ b/src/formulation/rawmat.cpp @@ -346,6 +346,16 @@ indexmat rawmat::getdinds(void) return Dinds; } +densemat rawmat::getavals(void) +{ + return Avals; +} + +densemat rawmat::getdvals(void) +{ + return Dvals; +} + Mat rawmat::getapetsc(void) { return Amat; diff --git a/src/formulation/rawmat.h b/src/formulation/rawmat.h index a71efde8..e912efdd 100755 --- a/src/formulation/rawmat.h +++ b/src/formulation/rawmat.h @@ -84,6 +84,8 @@ class rawmat indexmat getainds(void); indexmat getdinds(void); + densemat getavals(void); + densemat getdvals(void); Mat getapetsc(void); Mat getdpetsc(void); From 2adc84ac7936269c1021c4bcbd25451751e6d05c Mon Sep 17 00:00:00 2001 From: benjamin menkuec Date: Mon, 14 Mar 2022 19:51:25 +0100 Subject: [PATCH 5/7] added getarows getdrows getdcols getacols functions --- src/formulation/mat.cpp | 4 ++++ src/formulation/mat.h | 4 ++++ src/formulation/rawmat.cpp | 20 ++++++++++++++++++++ src/formulation/rawmat.h | 4 ++++ 4 files changed, 32 insertions(+) diff --git a/src/formulation/mat.cpp b/src/formulation/mat.cpp index 4e4a0d7a..fe293a3a 100755 --- a/src/formulation/mat.cpp +++ b/src/formulation/mat.cpp @@ -102,6 +102,10 @@ vec mat::eliminate(vec b) indexmat mat::getainds(void) { errorifpointerisnull(); errorifinvalidated(); return rawmatptr->getainds(); } indexmat mat::getdinds(void) { errorifpointerisnull(); errorifinvalidated(); return rawmatptr->getdinds(); } +indexmat mat::getarows(void) { errorifpointerisnull(); errorifinvalidated(); return rawmatptr->getarows(); } +indexmat mat::getdrows(void) { errorifpointerisnull(); errorifinvalidated(); return rawmatptr->getdrows(); } +indexmat mat::getacols(void) { errorifpointerisnull(); errorifinvalidated(); return rawmatptr->getacols(); } +indexmat mat::getdcols(void) { errorifpointerisnull(); errorifinvalidated(); return rawmatptr->getdcols(); } densemat mat::getavals(void) { errorifpointerisnull(); errorifinvalidated(); return rawmatptr->getavals(); } densemat mat::getdvals(void) { errorifpointerisnull(); errorifinvalidated(); return rawmatptr->getdvals(); } diff --git a/src/formulation/mat.h b/src/formulation/mat.h index 6d961334..362178d5 100755 --- a/src/formulation/mat.h +++ b/src/formulation/mat.h @@ -68,6 +68,10 @@ class mat indexmat getainds(void); indexmat getdinds(void); + indexmat getarows(void); + indexmat getdrows(void); + indexmat getacols(void); + indexmat getdcols(void); densemat getavals(void); densemat getdvals(void); diff --git a/src/formulation/rawmat.cpp b/src/formulation/rawmat.cpp index 6842159d..70d66a4d 100755 --- a/src/formulation/rawmat.cpp +++ b/src/formulation/rawmat.cpp @@ -346,6 +346,26 @@ indexmat rawmat::getdinds(void) return Dinds; } +indexmat rawmat::getarows(void) +{ + return Arows; +} + +indexmat rawmat::getdrows(void) +{ + return Drows; +} + +indexmat rawmat::getacols(void) +{ + return Acols; +} + +indexmat rawmat::getdcols(void) +{ + return Dcols; +} + densemat rawmat::getavals(void) { return Avals; diff --git a/src/formulation/rawmat.h b/src/formulation/rawmat.h index e912efdd..b7c81929 100755 --- a/src/formulation/rawmat.h +++ b/src/formulation/rawmat.h @@ -84,6 +84,10 @@ class rawmat indexmat getainds(void); indexmat getdinds(void); + indexmat getarows(void); + indexmat getdrows(void); + indexmat getacols(void); + indexmat getdcols(void); densemat getavals(void); densemat getdvals(void); From 2466f72e5e563696f2f25e1a7d39d0111a5cc6b9 Mon Sep 17 00:00:00 2001 From: benjamin menkuec Date: Thu, 17 Mar 2022 18:41:38 +0100 Subject: [PATCH 6/7] added another solve function for debugging purpose --- install_external_libs/install_petsc.sh | 1 - src/expression/operation/sl.cpp | 44 ++++++++++++++++++++++++++ src/expression/operation/sl.h | 1 + 3 files changed, 45 insertions(+), 1 deletion(-) diff --git a/install_external_libs/install_petsc.sh b/install_external_libs/install_petsc.sh index c78fa3b3..1655407e 100755 --- a/install_external_libs/install_petsc.sh +++ b/install_external_libs/install_petsc.sh @@ -37,7 +37,6 @@ fi # Metis is recommended but not mandatory. It can provide a major speedup for MUMPS during resolution. ./configure --with-openmp --with-mpi=0 --with-shared-libraries=1 --with-mumps-serial=1 --download-mumps --download-openblas --download-metis --download-slepc --with-debugging=0 --with-scalar-type=real --with-x=0 COPTFLAGS='-O3' CXXOPTFLAGS='-O3' FOPTFLAGS='-O3'; - ########## COMPILE PETSC : # In case petsc appends a 2 make a link: diff --git a/src/expression/operation/sl.cpp b/src/expression/operation/sl.cpp index 27915672..4b90f839 100755 --- a/src/expression/operation/sl.cpp +++ b/src/expression/operation/sl.cpp @@ -1525,6 +1525,50 @@ expression sl::array3x3(expression term11, expression term12, expression term13, return expression(3,3, {term11, term12, term13, term21, term22, term23, term31, term32, term33}); } +std::vector sl::solve(std::vector idxptr, std::vector indices, std::vector data, std::vector b) +{ + Mat Apetsc = PETSC_NULL; + MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, b.size(), b.size(), idxptr.data(), indices.data(), data.data(), &Apetsc); + MatAssemblyBegin(Apetsc, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(Apetsc, MAT_FINAL_ASSEMBLY); + + Vec bpetsc = PETSC_NULL; + VecCreate(PETSC_COMM_SELF, &bpetsc); + VecSetSizes(bpetsc, PETSC_DECIDE, b.size()); + VecSetFromOptions(bpetsc); + std::vector addr(b.size(), 0); + for (int i = 0; i < b.size(); i++) + addr[i] = i; + VecSetValues(bpetsc, b.size(), addr.data(), b.data(), INSERT_VALUES); + VecAssemblyBegin(bpetsc); + VecAssemblyEnd(bpetsc); + + Vec solpetsc = PETSC_NULL; + VecCreate(PETSC_COMM_SELF, &solpetsc); + VecSetSizes(solpetsc, PETSC_DECIDE, idxptr.size()-1); + VecSetFromOptions(solpetsc); + VecSet(solpetsc, 0.0); + VecAssemblyBegin(solpetsc); + VecAssemblyEnd(solpetsc); + + KSP myksp = PETSC_NULL; + KSP* ksp = &myksp; + KSPCreate(PETSC_COMM_SELF, ksp); + KSPSetOperators(*ksp, Apetsc, Apetsc); + KSPSetFromOptions(*ksp); + + PC pc; + KSPGetPC(*ksp,&pc); + PCSetType(pc,PCLU); + PCFactorSetMatSolverType(pc,MATSOLVERMUMPS); + + KSPSolve(*ksp, bpetsc, solpetsc); + + std::vector sol(b.size(), 0); + VecGetValues(solpetsc, b.size(), addr.data(), sol.data()); + return sol; +} + vec sl::solve(mat A, vec b, std::string soltype, bool diagscaling) { if (soltype != "lu" && soltype != "cholesky") diff --git a/src/expression/operation/sl.h b/src/expression/operation/sl.h index 8c1ec5df..5ec587c0 100755 --- a/src/expression/operation/sl.h +++ b/src/expression/operation/sl.h @@ -251,6 +251,7 @@ namespace sl expression array3x3(expression term11, expression term12, expression term13, expression term21, expression term22, expression term23, expression term31, expression term32, expression term33); + std::vector solve(std::vector idxptr, std::vector indices, std::vector data, std::vector b); // Direct resolution (with or without diagonal scaling): vec solve(mat A, vec b, std::string soltype = "lu", bool diagscaling = false); // Multi-rhs direct resolution: From 6a36e21f6150d419a1e7194c4d54dfa3495c0179 Mon Sep 17 00:00:00 2001 From: Benjamin Menkuec Date: Sun, 8 Jan 2023 14:57:26 +0100 Subject: [PATCH 7/7] Revert "added another solve function for debugging purpose" This reverts commit c8ef791584b1c7d46a9294c4a3b4ab444fda473b. --- install_external_libs/install_petsc.sh | 1 + src/expression/operation/sl.cpp | 44 -------------------------- src/expression/operation/sl.h | 1 - 3 files changed, 1 insertion(+), 45 deletions(-) diff --git a/install_external_libs/install_petsc.sh b/install_external_libs/install_petsc.sh index 1655407e..c78fa3b3 100755 --- a/install_external_libs/install_petsc.sh +++ b/install_external_libs/install_petsc.sh @@ -37,6 +37,7 @@ fi # Metis is recommended but not mandatory. It can provide a major speedup for MUMPS during resolution. ./configure --with-openmp --with-mpi=0 --with-shared-libraries=1 --with-mumps-serial=1 --download-mumps --download-openblas --download-metis --download-slepc --with-debugging=0 --with-scalar-type=real --with-x=0 COPTFLAGS='-O3' CXXOPTFLAGS='-O3' FOPTFLAGS='-O3'; + ########## COMPILE PETSC : # In case petsc appends a 2 make a link: diff --git a/src/expression/operation/sl.cpp b/src/expression/operation/sl.cpp index 4b90f839..27915672 100755 --- a/src/expression/operation/sl.cpp +++ b/src/expression/operation/sl.cpp @@ -1525,50 +1525,6 @@ expression sl::array3x3(expression term11, expression term12, expression term13, return expression(3,3, {term11, term12, term13, term21, term22, term23, term31, term32, term33}); } -std::vector sl::solve(std::vector idxptr, std::vector indices, std::vector data, std::vector b) -{ - Mat Apetsc = PETSC_NULL; - MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, b.size(), b.size(), idxptr.data(), indices.data(), data.data(), &Apetsc); - MatAssemblyBegin(Apetsc, MAT_FINAL_ASSEMBLY); - MatAssemblyEnd(Apetsc, MAT_FINAL_ASSEMBLY); - - Vec bpetsc = PETSC_NULL; - VecCreate(PETSC_COMM_SELF, &bpetsc); - VecSetSizes(bpetsc, PETSC_DECIDE, b.size()); - VecSetFromOptions(bpetsc); - std::vector addr(b.size(), 0); - for (int i = 0; i < b.size(); i++) - addr[i] = i; - VecSetValues(bpetsc, b.size(), addr.data(), b.data(), INSERT_VALUES); - VecAssemblyBegin(bpetsc); - VecAssemblyEnd(bpetsc); - - Vec solpetsc = PETSC_NULL; - VecCreate(PETSC_COMM_SELF, &solpetsc); - VecSetSizes(solpetsc, PETSC_DECIDE, idxptr.size()-1); - VecSetFromOptions(solpetsc); - VecSet(solpetsc, 0.0); - VecAssemblyBegin(solpetsc); - VecAssemblyEnd(solpetsc); - - KSP myksp = PETSC_NULL; - KSP* ksp = &myksp; - KSPCreate(PETSC_COMM_SELF, ksp); - KSPSetOperators(*ksp, Apetsc, Apetsc); - KSPSetFromOptions(*ksp); - - PC pc; - KSPGetPC(*ksp,&pc); - PCSetType(pc,PCLU); - PCFactorSetMatSolverType(pc,MATSOLVERMUMPS); - - KSPSolve(*ksp, bpetsc, solpetsc); - - std::vector sol(b.size(), 0); - VecGetValues(solpetsc, b.size(), addr.data(), sol.data()); - return sol; -} - vec sl::solve(mat A, vec b, std::string soltype, bool diagscaling) { if (soltype != "lu" && soltype != "cholesky") diff --git a/src/expression/operation/sl.h b/src/expression/operation/sl.h index 5ec587c0..8c1ec5df 100755 --- a/src/expression/operation/sl.h +++ b/src/expression/operation/sl.h @@ -251,7 +251,6 @@ namespace sl expression array3x3(expression term11, expression term12, expression term13, expression term21, expression term22, expression term23, expression term31, expression term32, expression term33); - std::vector solve(std::vector idxptr, std::vector indices, std::vector data, std::vector b); // Direct resolution (with or without diagonal scaling): vec solve(mat A, vec b, std::string soltype = "lu", bool diagscaling = false); // Multi-rhs direct resolution: