diff --git a/src/solver/solver_vector_petsc.cc b/src/solver/solver_vector_petsc.cc index 32e0b49af..cff66630c 100644 --- a/src/solver/solver_vector_petsc.cc +++ b/src/solver/solver_vector_petsc.cc @@ -1,294 +1,294 @@ /** * Copyright (©) 2019-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * This file is part of Akantu * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see <http://www.gnu.org/licenses/>. */ /* -------------------------------------------------------------------------- */ #include "solver_vector_petsc.hh" #include "dof_manager_petsc.hh" #include "mpi_communicator_data.hh" /* -------------------------------------------------------------------------- */ #include <numeric> #include <petscvec.h> /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ SparseSolverVectorPETSc::SparseSolverVectorPETSc(DOFManagerPETSc & dof_manager, const ID & id) : SparseSolverVector(dof_manager, id) { auto && mpi_comm = dof_manager.getMPIComm(); PETSc_call(VecCreate, mpi_comm, &x); detail::PETScSetName(x, id); PETSc_call(VecSetFromOptions, x); auto local_system_size = dof_manager.getLocalSystemSize(); auto nb_local_dofs = dof_manager.getPureLocalSystemSize(); PETSc_call(VecSetSizes, x, nb_local_dofs, PETSC_DECIDE); VecType vec_type; PETSc_call(VecGetType, x, &vec_type); if (std::string(vec_type) == std::string(VECMPI)) { PetscInt lowest_gidx; PetscInt highest_gidx; PETSc_call(VecGetOwnershipRange, x, &lowest_gidx, &highest_gidx); std::vector<PetscInt> ghost_idx; for (auto && d : arange(local_system_size)) { int gidx = dof_manager.localToGlobalEquationNumber(d); if (gidx != -1) { if ((gidx < lowest_gidx) or (gidx >= highest_gidx)) { ghost_idx.push_back(gidx); } } } PETSc_call(VecMPISetGhost, x, ghost_idx.size(), ghost_idx.data()); } else { std::vector<int> idx(nb_local_dofs); std::iota(idx.begin(), idx.end(), 0); ISLocalToGlobalMapping is; PETSc_call(ISLocalToGlobalMappingCreate, PETSC_COMM_SELF, 1, idx.size(), idx.data(), PETSC_COPY_VALUES, &is); PETSc_call(VecSetLocalToGlobalMapping, x, is); PETSc_call(ISLocalToGlobalMappingDestroy, &is); } } /* -------------------------------------------------------------------------- */ SparseSolverVectorPETSc:: SparseSolverVectorPETSc( // NOLINT(bugprone-copy-constructor-init) const SparseSolverVectorPETSc & vector, const ID & id) : SparseSolverVector(vector, id) { if (vector.x != nullptr) { PETSc_call(VecDuplicate, vector.x, &x); PETSc_call(VecCopy, vector.x, x); detail::PETScSetName(x, id); } } /* -------------------------------------------------------------------------- */ void SparseSolverVectorPETSc::printself(std::ostream & stream, int indent) const { std::string space(indent, AKANTU_INDENT); stream << space << "SolverVectorPETSc [" << std::endl; stream << space << " + id: " << id << std::endl; PETSc_call(PetscViewerPushFormat, PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INDEX); PETSc_call(VecView, x, PETSC_VIEWER_STDOUT_WORLD); PETSc_call(PetscViewerPopFormat, PETSC_VIEWER_STDOUT_WORLD); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ SparseSolverVectorPETSc::SparseSolverVectorPETSc(Vec x, DOFManagerPETSc & dof_manager, const ID & id) : SparseSolverVector(dof_manager, id) { PETSc_call(VecDuplicate, x, &this->x); PETSc_call(VecCopy, x, this->x); detail::PETScSetName(x, id); } /* -------------------------------------------------------------------------- */ SparseSolverVectorPETSc::~SparseSolverVectorPETSc() { if (x != nullptr) { PETSc_call(VecDestroy, &x); } } /* -------------------------------------------------------------------------- */ void SparseSolverVectorPETSc::resize() { // the arrays are destroyed and recreated in the dof manager // resize is so not implemented } /* -------------------------------------------------------------------------- */ void SparseSolverVectorPETSc::set(Real val) { PETSc_call(VecSet, x, val); applyModifications(); } /* -------------------------------------------------------------------------- */ void SparseSolverVectorPETSc::applyModifications() { PETSc_call(VecAssemblyBegin, x); PETSc_call(VecAssemblyEnd, x); updateGhost(); } /* -------------------------------------------------------------------------- */ void SparseSolverVectorPETSc::updateGhost() { Vec x_ghosted{nullptr}; PETSc_call(VecGhostGetLocalForm, x, &x_ghosted); if (x_ghosted != nullptr) { PETSc_call(VecGhostUpdateBegin, x, INSERT_VALUES, SCATTER_FORWARD); PETSc_call(VecGhostUpdateEnd, x, INSERT_VALUES, SCATTER_FORWARD); } PETSc_call(VecGhostRestoreLocalForm, x, &x_ghosted); } /* -------------------------------------------------------------------------- */ void SparseSolverVectorPETSc::getValues(const Array<Int> & idx, Array<Real> & values) const { if (idx.empty()) { return; } ISLocalToGlobalMapping is_ltog_map; PETSc_call(VecGetLocalToGlobalMapping, x, &is_ltog_map); PetscInt n; Array<PetscInt> lidx(idx.size()); PETSc_call(ISGlobalToLocalMappingApply, is_ltog_map, IS_GTOLM_MASK, idx.size(), idx.data(), &n, lidx.data()); getValuesLocal(lidx, values); } /* -------------------------------------------------------------------------- */ void SparseSolverVectorPETSc::getValuesLocal(const Array<Int> & idx, Array<Real> & values) const { if (idx.empty()) { return; } Vec x_ghosted{nullptr}; PETSc_call(VecGhostGetLocalForm, x, &x_ghosted); // VecScatterBegin(scatter, x, x_local, INSERT_VALUES, SCATTER_FORWARD); // VecScatterEnd(scatter, x, x_local, INSERT_VALUES, SCATTER_FORWARD); if (x_ghosted == nullptr) { const PetscScalar * array; PETSc_call(VecGetArrayRead, x, &array); for (auto && data : zip(idx, make_view(values))) { auto i = std::get<0>(data); if (i != -1) { std::get<1>(data) = array[i]; } } PETSc_call(VecRestoreArrayRead, x, &array); return; } PETSc_call(VecSetOption, x_ghosted, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE); PETSc_call(VecGetValues, x_ghosted, idx.size(), idx.data(), values.data()); PETSc_call(VecGhostRestoreLocalForm, x, &x_ghosted); } /* -------------------------------------------------------------------------- */ void SparseSolverVectorPETSc::addValues(const Array<Int> & gidx, const Array<Real> & values, Real scale_factor) { auto to_add = values.data(); Array<Real> scaled_array; if (scale_factor != 1.) { scaled_array.copy(values, false); scaled_array *= scale_factor; to_add = scaled_array.data(); } PETSc_call(VecSetOption, x, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE); PETSc_call(VecSetValues, x, gidx.size(), gidx.data(), to_add, ADD_VALUES); applyModifications(); } /* -------------------------------------------------------------------------- */ void SparseSolverVectorPETSc::addValuesLocal(const Array<Int> & lidx, const Array<Real> & values, Real scale_factor) { Vec x_ghosted{nullptr}; PETSc_call(VecGhostGetLocalForm, x, &x_ghosted); if (x_ghosted == nullptr) { auto to_add = values.data(); Array<Real> scaled_array; if (scale_factor != 1.) { scaled_array.copy(values, false); scaled_array *= scale_factor; to_add = scaled_array.data(); } PETSc_call(VecSetOption, x, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE); PETSc_call(VecSetValuesLocal, x, lidx.size(), lidx.data(), to_add, ADD_VALUES); return; } PETSc_call(VecGhostRestoreLocalForm, x, &x_ghosted); ISLocalToGlobalMapping is_ltog_map; PETSc_call(VecGetLocalToGlobalMapping, x, &is_ltog_map); Array<Int> gidx(lidx.size()); PETSc_call(ISLocalToGlobalMappingApply, is_ltog_map, lidx.size(), lidx.data(), gidx.data()); addValues(gidx, values, scale_factor); } /* -------------------------------------------------------------------------- */ SparseSolverVectorPETSc::operator const Array<Real> &() const { const_cast<Array<Real> &>(cache).resize(local_size()); auto xl = internal::make_petsc_local_vector(x); auto cachep = internal::make_petsc_wraped_vector(this->cache); - PETSc_call(VecCopy, cachep, xl); + PETSc_call(VecCopy, xl, cachep); return cache; } /* -------------------------------------------------------------------------- */ SparseSolverVectorPETSc & SparseSolverVectorPETSc::operator=(const SparseSolverVectorPETSc & y) { if (size() != y.size()) { PETSc_call(VecDuplicate, y, &x); } PETSc_call(VecCopy, y.x, x); release_ = y.release_; return *this; } /* -------------------------------------------------------------------------- */ SparseSolverVector & SparseSolverVectorPETSc::copy(const SparseSolverVector & y) { const auto & y_ = aka::as_type<SparseSolverVectorPETSc>(y); return operator=(y_); } /* -------------------------------------------------------------------------- */ SparseSolverVector & SparseSolverVectorPETSc::operator+(const SparseSolverVector & y) { const auto & y_ = aka::as_type<SparseSolverVectorPETSc>(y); PETSc_call(VecAXPY, x, 1., y_.x); release_ = y_.release_; return *this; } bool SparseSolverVectorPETSc::isFinite() const { Real max, min; - PETSc_call(VecMax, x, PETSC_NULL, &max); - PETSc_call(VecMin, x, PETSC_NULL, &min); + PETSc_call(VecMax, x, PETSC_NULLPTR, &max); + PETSc_call(VecMin, x, PETSC_NULLPTR, &min); return std::isfinite(min) and std::isfinite(max); } } // namespace akantu diff --git a/src/solver/solver_vector_petsc.hh b/src/solver/solver_vector_petsc.hh index d8287cf37..4c82ee173 100644 --- a/src/solver/solver_vector_petsc.hh +++ b/src/solver/solver_vector_petsc.hh @@ -1,202 +1,203 @@ /** * Copyright (©) 2019-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * This file is part of Akantu * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see <http://www.gnu.org/licenses/>. */ /* -------------------------------------------------------------------------- */ #include "dof_manager_petsc.hh" #include "solver_vector.hh" /* -------------------------------------------------------------------------- */ #include <petscvec.h> /* -------------------------------------------------------------------------- */ #ifndef AKANTU_SOLVER_VECTOR_PETSC_HH_ #define AKANTU_SOLVER_VECTOR_PETSC_HH_ namespace akantu { class DOFManagerPETSc; } // namespace akantu namespace akantu { /* -------------------------------------------------------------------------- */ namespace internal { /* ------------------------------------------------------------------------ */ class PETScVector { public: virtual ~PETScVector() = default; operator Vec &() { return x; } operator const Vec &() const { return x; } Int size() const { PetscInt n; PETSc_call(VecGetSize, x, &n); return n; } Int local_size() const { PetscInt n; PETSc_call(VecGetLocalSize, x, &n); return n; } AKANTU_GET_MACRO_NOT_CONST(Vec, x, auto &); AKANTU_GET_MACRO(Vec, x, const auto &); protected: Vec x{nullptr}; }; } // namespace internal /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ class SparseSolverVectorPETSc : public SparseSolverVector, public internal::PETScVector { public: SparseSolverVectorPETSc(DOFManagerPETSc & dof_manager, const ID & id = "solver_vector_petsc"); SparseSolverVectorPETSc(const SparseSolverVectorPETSc & vector, const ID & id = "solver_vector_petsc"); SparseSolverVectorPETSc(Vec x, DOFManagerPETSc & dof_manager, const ID & id = "solver_vector_petsc"); ~SparseSolverVectorPETSc() override; // resize the vector to the size of the problem void resize() override; void set(Real val) override; operator const Array<Real> &() const override; SparseSolverVector & operator+(const SparseSolverVector & y) override; SparseSolverVector & copy(const SparseSolverVector & y) override; SparseSolverVectorPETSc & operator=(const SparseSolverVectorPETSc & y); /// get values using processors global indexes void getValues(const Array<Int> & idx, Array<Real> & values) const; /// get values using processors local indexes void getValuesLocal(const Array<Int> & idx, Array<Real> & values) const; /// adding values to the vector using the global indices void addValues(const Array<Int> & gidx, const Array<Real> & values, Real scale_factor = 1.); /// adding values to the vector using the local indices void addValuesLocal(const Array<Int> & lidx, const Array<Real> & values, Real scale_factor = 1.); Int size() const override { return internal::PETScVector::size(); } Int localSize() const override { return internal::PETScVector::local_size(); } void printself(std::ostream & stream, int indent = 0) const override; bool isDistributed() const override { return true; } bool isFinite() const override; protected: void applyModifications(); void updateGhost(); protected: // used for the conversion operator Array<Real> cache; }; /* -------------------------------------------------------------------------- */ namespace internal { /* ------------------------------------------------------------------------ */ template <class Array> class PETScWrapedVector : public PETScVector { public: PETScWrapedVector(Array && array) : array(array) { PETSc_call(VecCreateSeqWithArray, PETSC_COMM_SELF, 1, array.size(), array.data(), &x); } ~PETScWrapedVector() override { PETSc_call(VecDestroy, &x); } private: Array array; }; /* ------------------------------------------------------------------------ */ template <bool read_only> class PETScLocalVector : public PETScVector { public: PETScLocalVector(const Vec & g) : g(g) { + PETSc_call(VecCreateLocalVector, g, &x); PETSc_call(VecGetLocalVectorRead, g, x); } PETScLocalVector(const SparseSolverVectorPETSc & g) : PETScLocalVector(g.getVec()) {} ~PETScLocalVector() override { PETSc_call(VecRestoreLocalVectorRead, g, x); PETSc_call(VecDestroy, &x); } private: const Vec & g; }; template <> class PETScLocalVector<false> : public PETScVector { public: PETScLocalVector(Vec & g) : g(g) { PETSc_call(VecGetLocalVectorRead, g, x); } PETScLocalVector(SparseSolverVectorPETSc & g) : PETScLocalVector(g.getVec()) {} ~PETScLocalVector() override { PETSc_call(VecRestoreLocalVectorRead, g, x); PETSc_call(VecDestroy, &x); } private: Vec & g; }; /* ------------------------------------------------------------------------ */ template <class Array> decltype(auto) make_petsc_wraped_vector(Array && array) { return PETScWrapedVector<Array>(std::forward<Array>(array)); } template < typename V, std::enable_if_t<std::is_same<Vec, std::decay_t<V>>::value> * = nullptr> decltype(auto) make_petsc_local_vector(V && vec) { constexpr auto read_only = std::is_const<std::remove_reference_t<V>>::value; return PETScLocalVector<read_only>(vec); } template <typename V, std::enable_if_t<std::is_base_of< SparseSolverVector, std::decay_t<V>>::value> * = nullptr> decltype(auto) make_petsc_local_vector(V && vec) { constexpr auto read_only = std::is_const<std::remove_reference_t<V>>::value; return PETScLocalVector<read_only>( dynamic_cast< std::conditional_t<read_only, const SparseSolverVectorPETSc, SparseSolverVectorPETSc> &>(vec)); } } // namespace internal } // namespace akantu #endif /* AKANTU_SOLVER_VECTOR_PETSC_HH_ */ diff --git a/test/test_model/patch_tests/data/material.dat b/test/test_model/patch_tests/data/material.dat index 30e3f810e..7b94ffa66 100644 --- a/test/test_model/patch_tests/data/material.dat +++ b/test/test_model/patch_tests/data/material.dat @@ -1,6 +1,11 @@ material elastic [ name = steel rho = 7800 # density E = 2.1e11 # young's modulus nu = 0.0 # poisson's ratio +dof_manager_type = petsc +time_step_solver dynamic [ +name = tss +sparse_solver_type = petsc +] ] diff --git a/test/test_model/patch_tests/data/material_anisotropic_1.dat b/test/test_model/patch_tests/data/material_anisotropic_1.dat index 8d4afe6fc..b1059d47e 100644 --- a/test/test_model/patch_tests/data/material_anisotropic_1.dat +++ b/test/test_model/patch_tests/data/material_anisotropic_1.dat @@ -1,7 +1,12 @@ material elastic_anisotropic [ name = aluminum rho = 1.6465129043044597 #gramm/mol/Å C11 = 105.092023 #eV/Å^3 n1 = [-1] +dof_manager_type = petsc +time_step_solver dynamic [ +name = tss +sparse_solver_type = petsc +] ] diff --git a/test/test_model/patch_tests/data/material_anisotropic_2.dat b/test/test_model/patch_tests/data/material_anisotropic_2.dat index 50d2460d1..ab2665d72 100644 --- a/test/test_model/patch_tests/data/material_anisotropic_2.dat +++ b/test/test_model/patch_tests/data/material_anisotropic_2.dat @@ -1,13 +1,18 @@ material elastic_anisotropic [ name = aluminum rho = 1.6465129043044597 #gramm/mol/Å C11 = 105.092023 #eV/Å^3 C12 = 59.4637759 C22 = 105.092023 C33 = 30.6596356 n1 = [-1, 1] n2 = [ 1, 1] +dof_manager_type = petsc +time_step_solver dynamic [ +name = tss +sparse_solver_type = petsc +] ] diff --git a/test/test_model/patch_tests/data/material_anisotropic_3.dat b/test/test_model/patch_tests/data/material_anisotropic_3.dat index 550237ab4..b71d6a1a4 100644 --- a/test/test_model/patch_tests/data/material_anisotropic_3.dat +++ b/test/test_model/patch_tests/data/material_anisotropic_3.dat @@ -1,36 +1,41 @@ material elastic_anisotropic [ name = aluminum rho = 1.6465129043044597 #gramm/mol/Å C11 = 105.092023 #eV/Å^3 C12 = 59.4637759 C13 = 59.4637759 C14 = -9.47874286e-11 C15 = -9.34769857e-10 C16 = -3.49522561e-10 C22 = 105.092023 C23 = 59.4637759 C24 = -1.33279629e-10 C25 = 6.52290813e-11 C26 = -1.36100714e-10 C33 = 105.092023 C34 = 1.410383e-10 C35 = -3.02540691e-13 C36 = -8.72969656e-11 C44 = 30.6596356 C45 = 2.20340964e-11 C46 = -2.62259488e-10 C55 = 30.6596356 C56 = 8.56508319e-11 C66 = 30.6596356 n1 = [-1, 1, 0] n2 = [ 1, 1, 1] n3 = [ 1, 1, -2] # alpha = 1e-3 # viscous proportion useless unless manually activated in material_elastic_linear_anisotropic.cc +dof_manager_type = petsc +time_step_solver dynamic [ +name = tss +sparse_solver_type = petsc +] ] diff --git a/test/test_model/patch_tests/data/material_lumped.dat b/test/test_model/patch_tests/data/material_lumped.dat index 30e3f810e..7b94ffa66 100644 --- a/test/test_model/patch_tests/data/material_lumped.dat +++ b/test/test_model/patch_tests/data/material_lumped.dat @@ -1,6 +1,11 @@ material elastic [ name = steel rho = 7800 # density E = 2.1e11 # young's modulus nu = 0.0 # poisson's ratio +dof_manager_type = petsc +time_step_solver dynamic [ +name = tss +sparse_solver_type = petsc +] ]