diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..567609b --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +build/ diff --git a/CMakeLists.txt b/CMakeLists.txt index ead9106..bdd57b7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -3,7 +3,6 @@ project(ngspetsc) cmake_minimum_required(VERSION 3.8) -option (PETSC_WITH_HYPRE "Does PETSc have hypre enabled?" ON) option (PETSC_COMPLEX "Are we building the complex version of the interface?" OFF) set(CMAKE_MODULE_PATH "${CMAKE_MODULE_PATH}" "${CMAKE_CURRENT_SOURCE_DIR}/cmake/cmake_modules") @@ -32,12 +31,5 @@ else (NOT PETSC_COMPLEX) set (MODULE_NAME "ngs_petsc_complex") endif (NOT PETSC_COMPLEX) -if(NOT PETSC_WITH_HYPRE) - message(STATUS "Turning OFF Hypre Auxiliary space Preconditioners.") - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DPETSC_NO_HYPRE") -else(NOT PETSC_WITH_HYPRE) - message(STATUS "Turning ON Hypre Auxiliary space Preconditioners.") -endif(NOT PETSC_WITH_HYPRE) - add_subdirectory(src) add_subdirectory(python) diff --git a/README.rst b/README.rst index cb77cdd..db45cb4 100644 --- a/README.rst +++ b/README.rst @@ -11,8 +11,3 @@ override this with cmake -DPETSC_EXECUTABLE_RUNS=YES -If your PETSc installation has beed configured without hypre, you have to tell cmake: - - - cmake -DPETSC_WITH_HYPRE=OFF - diff --git a/src/petsc_interface.hpp b/src/petsc_interface.hpp index 8b8d81d..a588429 100644 --- a/src/petsc_interface.hpp +++ b/src/petsc_interface.hpp @@ -1,6 +1,11 @@ #include +// SZ +// consider removing this include from here #include "petsc.h" + +// SZ +// why you need python stuff in a C++ interface declaration? #include #include "typedefs.hpp" diff --git a/src/petsc_ksp.cpp b/src/petsc_ksp.cpp index ddb927f..929af4f 100644 --- a/src/petsc_ksp.cpp +++ b/src/petsc_ksp.cpp @@ -98,12 +98,8 @@ namespace ngs_petsc_interface KSPSetUp(GetKSP()); petsc_rhs = GetMatrix()->GetRowMap()->CreatePETScVector(); - // VecAssemblyBegin(petsc_rhs); - // VecAssemblyEnd(petsc_rhs); petsc_sol = GetMatrix()->GetColMap()->CreatePETScVector(); - // VecAssemblyBegin(petsc_sol); - // VecAssemblyEnd(petsc_sol); } diff --git a/src/petsc_ksp.hpp b/src/petsc_ksp.hpp index 0228ecc..f5b9fb9 100644 --- a/src/petsc_ksp.hpp +++ b/src/petsc_ksp.hpp @@ -1,6 +1,11 @@ #ifndef FILE_NGSPETSC_KSP_HPP #define FILE_NGSPETSC_KSP_HPP +// SZ +// Each specific hpp file should be able to resolve all the symbols defined in it +// So, at least for PETSc stuff, you can either include typedefs.hpp or forward declare the needed symbols (KSP, Mat, PC, etc) +// Since this is a possibly public header, you should avoid including petscksp.h + namespace ngs_petsc_interface { diff --git a/src/petsc_linalg.cpp b/src/petsc_linalg.cpp index 32fd1b0..6849cfa 100644 --- a/src/petsc_linalg.cpp +++ b/src/petsc_linalg.cpp @@ -125,6 +125,7 @@ namespace ngs_petsc_interface template PETScMat CreatePETScMatSeqBAIJFromSymmetric (shared_ptr> spmat, shared_ptr rss, shared_ptr css) { + // SZ: it seems you could interface to SBAIJ too static ngs::Timer t(string("CreatePETScMatSeqBAIJFromSymmetric::HEIGHT) + string(">>")); ngs::RegionTimer rt(t); // row map (map for a row) @@ -157,6 +158,8 @@ namespace ngs_petsc_interface } PETScMat petsc_mat; if (bh == 1) + // SZ: here and in all other similar calls: if you now the nonzero per row, then pass 0 (not -1) as + // total number of nonzeros { MatCreateSeqAIJ(PETSC_COMM_SELF, nrows, ncols, -1, &nzepr[1], &petsc_mat); } else { MatCreateSeqBAIJ(PETSC_COMM_SELF, bh, nrows, ncols, -1, &nzepr[1], &petsc_mat); } @@ -299,6 +302,13 @@ namespace ngs_petsc_interface } // CreatePETScMatSeqBAIJ + // SZ + // I would not use BAIJ by default if the block size > 2 + // AIJ is much more supported, BAIJ does not really deliver impressive performances improvements + // Some solvers (GAMG, BDDC, FIELDSPLIT) use block size information but this information + // is valid for AIJ too, if you use MatSetBlockSize(aijmat,bs); + // My suggestion here is to always default to AIJ and provide an extra argument to this + // function to instead use blocked formats. BTW, It seems SBAIJ support is missing, am I correct? PETScMat CreatePETScMatSeq (shared_ptr mat, shared_ptr rss, shared_ptr css) { if (auto spm = dynamic_pointer_cast>(mat)) @@ -758,6 +768,8 @@ namespace ngs_petsc_interface static ngs::Timer t("PETScMatrix::UpdateValues"); ngs::RegionTimer rt(t); // TODO: if we have converted the matrix from BAIJ to AIJ, is SetValuesBlocked inefficient?? + // SZ: answer: it shouldn't be inefficient: it expands by blocksize + // the blocked set of indices and then call MatSetValues auto parmat = dynamic_pointer_cast(ngs_mat); shared_ptr mat = (parmat == nullptr) ? ngs_mat : parmat->GetMatrix(); diff --git a/src/petsc_linalg.hpp b/src/petsc_linalg.hpp index 94f9112..acdc5c8 100644 --- a/src/petsc_linalg.hpp +++ b/src/petsc_linalg.hpp @@ -4,6 +4,8 @@ /** Matrices / Vectors **/ +// SZ +// See comments about header inclusions in petsc_ksp.hpp namespace ngs_petsc_interface { @@ -142,6 +144,7 @@ namespace ngs_petsc_interface Can be any kind of BaseMatrix **/ + // SZ An alternative is to have another Enum value in MAT_TYPE, namely SHELL, and you do not need a special class for it class FlatPETScMatrix : public PETScBaseMatrix { public: diff --git a/src/petsc_ops.hpp b/src/petsc_ops.hpp index 32d8b68..678351f 100644 --- a/src/petsc_ops.hpp +++ b/src/petsc_ops.hpp @@ -130,17 +130,22 @@ namespace ngs_petsc_interface void Ngs2PETSc (shared_ptr ngs_vec, ::Vec petsc_vec) { ngs_vec->Cumulate(); - VecAssemblyBegin(petsc_vec); auto fv = ngs_vec->FVDouble(); for (auto k : Range(loc)) buf[k] = fv(loc_inds[k]); VecSetValues(petsc_vec, loc, &glob_nums[0], &buf[0], INSERT_VALUES); + VecAssemblyBegin(petsc_vec); VecAssemblyEnd(petsc_vec); } // petsc -> ngsolve void PETSc2Ngs (shared_ptr ngs_vec, ::Vec petsc_vec) { ngs_vec->Distribute(); + // SZ + // This is odd, and rarely used from PETSc. Also, off-processor retrieval is disabled + // Can you use VecGetArrayRead/VecRestoreArrayRead semantics? + // In alternative, if Ngs2PETSc and PETSc2NGs needs off-process data, you can use + // the VecScatter object VecGetValues(petsc_vec, loc, &glob_nums[0], &buf[0]); auto fv = ngs_vec->FVDouble(); fv = 0.0; diff --git a/src/petsc_pc.cpp b/src/petsc_pc.cpp index 645b69d..ec986d4 100644 --- a/src/petsc_pc.cpp +++ b/src/petsc_pc.cpp @@ -189,7 +189,7 @@ namespace ngs_petsc_interface } -#ifdef PETSC_HAS_HYPRE +#ifdef PETSC_HAVE_HYPRE PETScHypreAuxiliarySpacePC :: PETScHypreAuxiliarySpacePC (shared_ptr _bfa, const ngs::Flags & _aflags, const string _aname) : PETSc2NGsPrecond(_bfa, _aflags, _aname) { @@ -455,6 +455,8 @@ namespace ngs_petsc_interface // Now set the PCs we already had before KSP* ksps; PetscInt n; + // SZ: be careful with this function, since those KSPs are different depending on the + // fieldsplit type, see the online documentation PCFieldSplitGetSubKSP(GetPETScPC(), &n, &ksps); for (auto k : Range(fields.Size())) { auto field = fields[k]; @@ -471,7 +473,7 @@ namespace ngs_petsc_interface ngs::RegisterPreconditioner registerPETSc2NGsPrecond("petsc_pc"); -#ifdef PETSC_HAS_HYPRE +#ifdef PETSC_HAVE_HYPRE ngs::RegisterPreconditioner registerPETScHypreAMS("petsc_pc_hypre_ams"); #endif @@ -513,7 +515,7 @@ namespace ngs_petsc_interface .def("Finalize", [](shared_ptr & pc) { pc->FinalizeLevel(); } ); -#ifdef PETSC_HAS_HYPRE +#ifdef PETSC_HAVE_HYPRE py::class_, PETSc2NGsPrecond> (m, "HypreAuxiliarySpacePrecond", "ADS/AMS from hypre package") .def ("SetGradientMatrix" , [](shared_ptr & pc, shared_ptr mat) diff --git a/src/petsc_pc.hpp b/src/petsc_pc.hpp index a49a6cb..883e87d 100644 --- a/src/petsc_pc.hpp +++ b/src/petsc_pc.hpp @@ -44,6 +44,8 @@ namespace ngs_petsc_interface // }; + // SZ + // Why do you inherit from FlatPETScMatrix? /** An NGSolve-BaseMatrix, wrapped to PETSc as a PC **/ class NGs2PETScPrecond : public PETScBasePrecond, public FlatPETScMatrix @@ -104,7 +106,7 @@ namespace ngs_petsc_interface Array> keep_alive; }; -#ifdef PETSC_HAS_HYPRE +#ifdef PETSC_HAVE_HYPRE class PETScHypreAuxiliarySpacePC : public PETSc2NGsPrecond { diff --git a/src/petsc_snes.cpp b/src/petsc_snes.cpp index a49d5ee..6613466 100644 --- a/src/petsc_snes.cpp +++ b/src/petsc_snes.cpp @@ -101,6 +101,7 @@ namespace ngs_petsc_interface { SetOptions (_opts, "", NULL); } // Create Vector to hold F(x) + // SZ: this not strictly needed, since Petsc will create a Vec for you if you pass NULL func_vec = jac_mat->GetRowMap()->CreatePETScVector(); // Set (non-linear) function evaluation, f = F(x) @@ -234,9 +235,18 @@ namespace ngs_petsc_interface { auto& self = *( (PETScSNES*) ctx); HeapReset hr(*self.use_lh); - + + // SZ: (do you really want users reading WTF in your code???) // actually, this happens if snes_mf or snes_mf_operator are set // if(A != self.jac_mat->GetPETScMat()) { throw Exception("A != JAC BUFFER, WTF?"); } + + // SZ: B is the correct choice here. However, I think to cover all possible matrix-free + // options, you need to add the following 4 lines (this will reset the shell matrix to its proper state) + // + //if (A && A != B) { + // ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); + // ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); + //} if(B != self.jac_mat->GetPETScMat()) { throw Exception("B != JAC BUFFER, WTF?"); } self.jac_mat->GetRowMap()->PETSc2NGs(*self.lin_vec, x); diff --git a/src/typedefs.hpp b/src/typedefs.hpp index a15eea5..64f018a 100644 --- a/src/typedefs.hpp +++ b/src/typedefs.hpp @@ -4,10 +4,6 @@ namespace ngs_petsc_interface { -#ifndef PETSC_NO_HYPRE -#define PETSC_HAS_HYPRE -#endif - namespace ngs = ngcomp; using ngs::Array; using ngs::Range; @@ -32,6 +28,8 @@ namespace ngs_petsc_interface #else static_assert( is_same::value, "Trying to compile the real interface with a complex PETSc installation, (set -DPETSC_COMPLEX=OFF)!"); #endif +//SZ +// Not sure if you support PETSc compiled with 64bit indices // static_assert( (is_same::value || is_same::value), "Need double or complex PETSc version!"); // static_assert( (is_same::value), "Not a double PETSc version!");