diff --git a/packages/mumps.cmake b/packages/mumps.cmake index 5d28887d7..aa441b243 100644 --- a/packages/mumps.cmake +++ b/packages/mumps.cmake @@ -1,111 +1,110 @@ #=============================================================================== # @file mumps.cmake # # @author Nicolas Richart # # @date creation: Mon Nov 21 2011 # @date last modification: Mon Jan 18 2016 # # @brief package description for mumps support # # @section LICENSE # # Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de # Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des # Solides) # # 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 . # #=============================================================================== - package_declare(Mumps EXTERNAL DESCRIPTION "Add Mumps support in akantu" SYSTEM ON third-party/cmake/mumps.cmake ) package_declare_sources(Mumps solver/sparse_solver_mumps.cc solver/sparse_solver_mumps.hh ) set(_mumps_float_type ${AKANTU_FLOAT_TYPE}) if(AKANTU_FLOAT_TYPE STREQUAL "float" OR AKANTU_FLOAT_TYPE STREQUAL "double") set(_mumps_components ${AKANTU_FLOAT_TYPE}) else() if(DEFINED AKANTU_FLOAT_TYPE) message(FATAL_ERROR "MUMPS doea not support floating point type \"${AKANTU_FLOAT_TYPE}\"") endif() endif() package_get_option_name(parallel _par_option) if(${_par_option}) package_set_find_package_extra_options(Mumps ARGS COMPONENTS "parallel" ${_mumps_components}) package_add_third_party_script_variable(Mumps MUMPS_TYPE "par") package_set_package_system_dependency(Mumps deb libmumps) package_set_package_system_dependency(Mumps deb-src libmumps-dev) else() package_set_find_package_extra_options(Mumps ARGS COMPONENTS "sequential" ${_mumps_components}) package_add_third_party_script_variable(Mumps MUMPS_TYPE "seq") package_set_package_system_dependency(Mumps deb libmumps-seq) package_set_package_system_dependency(Mumps deb-src libmumps-seq-dev) endif() package_use_system(Mumps _use_system) if(NOT _use_system) enable_language(Fortran) set(AKANTU_USE_MUMPS_VERSION "4.10.0" CACHE STRING "Default Mumps version to compile") mark_as_advanced(AKANTU_USE_MUMPS_VERSION) set_property(CACHE AKANTU_USE_MUMPS_VERSION PROPERTY STRINGS "4.9.2" "4.10.0" "5.0.0") package_get_option_name(MPI _mpi_option) if(${_mpi_option}) package_add_dependencies(Mumps ScaLAPACK MPI) endif() package_add_dependencies(Mumps Scotch BLAS) endif() package_declare_documentation(Mumps "This package enables the \\href{http://mumps.enseeiht.fr/}{MUMPS} parallel direct solver for sparce matrices." "This is necessary to solve static or implicit problems." "" "Under Ubuntu (14.04 LTS) the installation can be performed using the commands:" "" "\\begin{command}" " > sudo apt-get install libmumps-seq-dev # for sequential" " > sudo apt-get install libmumps-dev # for parallel" "\\end{command}" "" "Under Mac OS X the installation requires the following steps:" "\\begin{command}" " > sudo port install mumps" "\\end{command}" "" "If you activate the advanced option AKANTU\\_USE\\_THIRD\\_PARTY\\_MUMPS the make system of akantu can automatically compile MUMPS. For this you will have to download MUMPS from \\url{http://mumps.enseeiht.fr/} or \\url{http://graal.ens-lyon.fr/MUMPS} and place it in \\shellcode{/third-party}" ) package_declare_extra_files_to_package(MUMPS PROJECT third-party/MUMPS_4.10.0_make.inc.cmake third-party/MUMPS_5.0.0.patch third-party/MUMPS_4.10.0.patch third-party/MUMPS_4.9.2_make.inc.cmake third-party/cmake/mumps.cmake cmake/Modules/FindMumps.cmake ) diff --git a/src/solver/solver_vector_petsc.cc b/src/solver/solver_vector_petsc.cc index 949917c6e..2097bca51 100644 --- a/src/solver/solver_vector_petsc.cc +++ b/src/solver/solver_vector_petsc.cc @@ -1,85 +1,94 @@ /** * @file solver_vector_petsc.cc * * @author Nicolas Richart * * @date creation Tue Jan 01 2019 * * @brief A Documented file. * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * 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 . * */ /* -------------------------------------------------------------------------- */ #include "solver_vector_petsc.hh" #include "dof_manager_petsc.hh" #include "mpi_communicator_data.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ SolverVectorPETSc::SolverVectorPETSc(DOFManagerPETSc & dof_manager, const ID & id) : SolverVector(dof_manager, id), dof_manager(dof_manager) { auto && mpi_comm = dof_manager.getMPIComm(); PETSc_call(VecCreate, mpi_comm, &vector); PETSc_call(PetscObjectSetName, reinterpret_cast(vector), id.c_str()); } /* -------------------------------------------------------------------------- */ SolverVectorPETSc::SolverVectorPETSc(const SolverVectorPETSc & vector, const ID & id) : SolverVector(vector, id), dof_manager(vector.dof_manager) { PETSc_call(VecDuplicate, vector.vector, &this->vector); PETSc_call(PetscObjectSetName, reinterpret_cast(this->vector), id.c_str()); PETSc_call(VecCopy, vector.vector, this->vector); } +/* -------------------------------------------------------------------------- */ +SolverVectorPETSc::SolverVectorPETSc(Vec vec, + const ID & id = "solver_vector_petsc") { + PETSc_call(VecDuplicate, vec.vector, &vector); + PETSc_call(PetscObjectSetName, reinterpret_cast(this->vector), id.c_str()); + PETSc_call(VecCopy, vector.vector, this->vector); +} + + /* -------------------------------------------------------------------------- */ SolverVectorPETSc::~SolverVectorPETSc() { if (vector) { PETSc_call(VecDestroy, &vector); } } /* -------------------------------------------------------------------------- */ void SolverVectorPETSc::resize() { PETSc_call(VecSetSizes, vector, dof_manager.getPureLocalSystemSize(), dof_manager.getSystemSize()); auto & is_ltog_mapping = dof_manager.getISLocalToGlobalMapping(); PETSc_call(VecSetLocalToGlobalMapping, vector, is_ltog_mapping); } /* -------------------------------------------------------------------------- */ void SolverVectorPETSc::clear() { PETSc_call(VecSet, vector, 0.); applyModifications(); } /* -------------------------------------------------------------------------- */ void SolverVectorPETSc::applyModifications() { PETSc_call(VecAssemblyBegin, vector); PETSc_call(VecAssemblyEnd, vector); } } // namespace akantu diff --git a/src/solver/solver_vector_petsc.hh b/src/solver/solver_vector_petsc.hh index 88d172bcb..63fb6c786 100644 --- a/src/solver/solver_vector_petsc.hh +++ b/src/solver/solver_vector_petsc.hh @@ -1,70 +1,73 @@ /** * @file solver_vector_petsc.hh * * @author Nicolas Richart * * @date creation Tue Jan 01 2019 * * @brief A Documented file. * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * 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 . * */ /* -------------------------------------------------------------------------- */ #include "solver_vector.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLVER_VECTOR_PETSC_HH__ #define __AKANTU_SOLVER_VECTOR_PETSC_HH__ namespace akantu { class DOFManagerPETSc; } // namespace akantu namespace akantu { class SolverVectorPETSc : public SolverVector { public: SolverVectorPETSc(DOFManagerPETSc & dof_manager, const ID & id = "solver_vector_petsc"); SolverVectorPETSc(const SolverVectorPETSc & vector, const ID & id = "solver_vector_petsc"); + SolverVectorPETSc(Vec vec, + const ID & id = "solver_vector_petsc"); + ~SolverVectorPETSc() override; // resize the vector to the size of the problem void resize() override; void clear() override; protected: void applyModifications(); public: AKANTU_GET_MACRO_NOT_CONST(Vector, vector, auto &); AKANTU_GET_MACRO(Vector, vector, const auto &); protected: DOFManagerPETSc & dof_manager; Vec vector; }; } // namespace akantu #endif /* __AKANTU_SOLVER_VECTOR_PETSC_HH__ */ diff --git a/src/solver/sparse_matrix.hh b/src/solver/sparse_matrix.hh index 77b8bac0f..d0469aca6 100644 --- a/src/solver/sparse_matrix.hh +++ b/src/solver/sparse_matrix.hh @@ -1,161 +1,161 @@ /** * @file sparse_matrix.hh * * @author Nicolas Richart * * @date creation: Mon Dec 13 2010 * @date last modification: Tue Feb 20 2018 * * @brief sparse matrix storage class (distributed assembled matrix) * This is a COO format (Coordinate List) * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * 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 . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SPARSE_MATRIX_HH__ #define __AKANTU_SPARSE_MATRIX_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { class DOFManager; class TermsToAssemble; class SolverVector; } namespace akantu { class SparseMatrix { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: SparseMatrix(DOFManager & dof_manager, const MatrixType & matrix_type, const ID & id = "sparse_matrix"); SparseMatrix(const SparseMatrix & matrix, const ID & id = "sparse_matrix"); virtual ~SparseMatrix(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// remove the existing profile virtual void clearProfile(); /// set the matrix to 0 virtual void clear() = 0; /// add a non-zero element to the profile virtual UInt add(UInt i, UInt j) = 0; /// assemble a local matrix in the sparse one virtual void add(UInt i, UInt j, Real value) = 0; /// save the profil in a file using the MatrixMarket file format virtual void saveProfile(__attribute__((unused)) const std::string &) const { AKANTU_TO_IMPLEMENT(); } /// save the matrix in a file using the MatrixMarket file format virtual void saveMatrix(__attribute__((unused)) const std::string &) const { AKANTU_TO_IMPLEMENT(); }; /// multiply the matrix by a coefficient virtual void mul(Real alpha) = 0; /// add matrices virtual void add(const SparseMatrix & matrix, Real alpha = 1.); /// Equivalent of *gemv in blas virtual void matVecMul(const SolverVector & x, SolverVector & y, Real alpha = 1., Real beta = 0.) const = 0; /// modify the matrix to "remove" the blocked dof virtual void applyBoundary(Real block_val = 1.) = 0; /// operator *= SparseMatrix & operator*=(Real alpha) { this->mul(alpha); return *this; } protected: /// This is the revert of add B += \alpha * *this; virtual void addMeTo(SparseMatrix & B, Real alpha) const = 0; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the values at potition i, j virtual inline Real operator()(__attribute__((unused)) UInt i, __attribute__((unused)) UInt j) const { AKANTU_TO_IMPLEMENT(); } /// return the values at potition i, j virtual inline Real & operator()(__attribute__((unused)) UInt i, __attribute__((unused)) UInt j) { AKANTU_TO_IMPLEMENT(); } AKANTU_GET_MACRO(NbNonZero, nb_non_zero, UInt); UInt size() const { return size_; } AKANTU_GET_MACRO(MatrixType, matrix_type, const MatrixType &); virtual UInt getRelease() const = 0; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: ID id; /// Underlying dof manager DOFManager & _dof_manager; /// sparce matrix type MatrixType matrix_type; /// Size of the matrix UInt size_; /// number of processors UInt nb_proc; /// number of non zero element UInt nb_non_zero; }; -Array & operator*=(Array & vect, const SparseMatrix & mat); +//Array & operator*=(Array & vect, const SparseMatrix & mat); } // akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "sparse_matrix_inline_impl.cc" #endif /* __AKANTU_SPARSE_MATRIX_HH__ */ diff --git a/src/solver/sparse_matrix_petsc.cc b/src/solver/sparse_matrix_petsc.cc index 9d87a0210..5e601d44d 100644 --- a/src/solver/sparse_matrix_petsc.cc +++ b/src/solver/sparse_matrix_petsc.cc @@ -1,252 +1,253 @@ /** * @file sparse_matrix_petsc.cc * * @author Aurelia Isabel Cuba Ramos * * @date creation: Mon Dec 13 2010 * @date last modification: Sat Feb 03 2018 * * @brief Implementation of PETSc matrix class * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * 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 . * */ /* -------------------------------------------------------------------------- */ #include "sparse_matrix_petsc.hh" #include "dof_manager_petsc.hh" #include "mpi_communicator_data.hh" #include "solver_vector_petsc.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ SparseMatrixPETSc::SparseMatrixPETSc(DOFManagerPETSc & dof_manager, const MatrixType & sparse_matrix_type, const ID & id) : SparseMatrix(dof_manager, matrix_type, id), dof_manager(dof_manager) { AKANTU_DEBUG_IN(); auto mpi_comm = dof_manager.getMPIComm(); PETSc_call(MatCreate, mpi_comm, &mat); - PETSc_call(PetscObjectSetName, reinterpret_cast(mat), - id.c_str()); if (sparse_matrix_type == _symmetric) PETSc_call(MatSetOption, mat, MAT_SYMMETRIC, PETSC_TRUE); PETSc_call(MatSetOption, mat, MAT_ROW_ORIENTED, PETSC_TRUE); - PETSc_call(MatSetSizes, mat, PETSC_DECIDE, PETSC_DECIDE, size_, size_); + auto local_size = dof_manager.getPureLocalSystemSize(); + PETSc_call(MatSetSizes, mat, local_size, local_size, size_, size_); auto & is_ltog_mapping = dof_manager.getISLocalToGlobalMapping(); + PETSc_call(MatSetUp, mat); + PETSc_call(MatSetLocalToGlobalMapping, mat, is_ltog_mapping, is_ltog_mapping); + PETSc_call(PetscObjectSetName, reinterpret_cast(mat), + id.c_str()); + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SparseMatrixPETSc::SparseMatrixPETSc(const SparseMatrixPETSc & matrix, const ID & id) : SparseMatrix(matrix, id), dof_manager(matrix.dof_manager) { PETSc_call(MatDuplicate, matrix.mat, MAT_COPY_VALUES, &mat); PETSc_call(PetscObjectSetName, reinterpret_cast(mat), id.c_str()); } /* -------------------------------------------------------------------------- */ SparseMatrixPETSc::~SparseMatrixPETSc() { AKANTU_DEBUG_IN(); if (mat) PETSc_call(MatDestroy, &mat); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Method to save the nonzero pattern and the values stored at each position * @param filename name of the file in which the information will be stored */ void SparseMatrixPETSc::saveMatrix(const std::string & filename) const { AKANTU_DEBUG_IN(); auto mpi_comm = dof_manager.getMPIComm(); /// create Petsc viewer PetscViewer viewer; PETSc_call(PetscViewerASCIIOpen, mpi_comm, filename.c_str(), &viewer); - - /// set the format - // PETSc_call(PetscViewerSetFormat, viewer, PETSC_VIEWER_DEFAULT); - - /// save the matrix - /// @todo Write should be done in serial -> might cause problems + PETSc_call(PetscViewerPushFormat, viewer, PETSC_VIEWER_ASCII_MATRIXMARKET); PETSc_call(MatView, mat, viewer); - - /// destroy the viewer + PETSc_call(PetscViewerPopFormat, viewer); PETSc_call(PetscViewerDestroy, &viewer); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /// Equivalent of *gemv in blas void SparseMatrixPETSc::matVecMul(const SolverVector & _x, SolverVector & _y, Real alpha, Real beta) const { auto & x = dynamic_cast(_x).getVector(); auto & y = dynamic_cast(_y).getVector(); - PETSc_call(MatMult, mat, x, y); + // y = alpha A x + beta y + SolverVectorPETSc w(x); + // w = A x + PETSc_call(MatMult, mat, x, w); if (alpha != 1.) { - PETSc_call(VecScale, y, alpha); + // w = alpha w + PETSc_call(VecScale, w, alpha); } - PETSc_call(VecAYPX, y, beta, x); + // y = w + beta y + PETSc_call(VecAYPX, y, beta, w); } /* -------------------------------------------------------------------------- */ void SparseMatrixPETSc::addMeToImpl(SparseMatrixPETSc & B, Real alpha) const { PETSc_call(MatAXPY, B.mat, alpha, mat, SAME_NONZERO_PATTERN); B.release++; } /* -------------------------------------------------------------------------- */ /** * Method to add another PETSc matrix to this PETSc matrix * @param matrix PETSc matrix to be added * @param alpha the factor specifying how many times the matrix should be added */ void SparseMatrixPETSc::addMeTo(SparseMatrix & B, Real alpha) const { if (auto * B_petsc = dynamic_cast(&B)) { this->addMeToImpl(*B_petsc, alpha); } else { AKANTU_TO_IMPLEMENT(); // this->addMeToTemplated(*this, alpha); } } /* -------------------------------------------------------------------------- */ /** * MatSetValues() generally caches the values. The matrix is ready to * use only after MatAssemblyBegin() and MatAssemblyEnd() have been * called. (http://www.mcs.anl.gov/petsc/) */ void SparseMatrixPETSc::performAssembly() { this->beginAssembly(); this->endAssembly(); this->release++; } /* -------------------------------------------------------------------------- */ void SparseMatrixPETSc::beginAssembly() { PETSc_call(MatAssemblyBegin, mat, MAT_FINAL_ASSEMBLY); } /* -------------------------------------------------------------------------- */ void SparseMatrixPETSc::endAssembly() { PETSc_call(MatAssemblyEnd, mat, MAT_FINAL_ASSEMBLY); - PETSc_call(MatSetOption, mat, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE); + PETSc_call(MatSetOption, mat, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE); // PETSc_call(MatSetOption, mat, MAT_NEW_NONZERO_ALLOCATIONS, PETSC_TRUE); // PETSc_call(MatSetOption, mat, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE); } /* -------------------------------------------------------------------------- */ -void SparseMatrixPETSc::applyBoundary(Real /*block_val*/) { +void SparseMatrixPETSc::applyBoundary(Real block_val) { AKANTU_DEBUG_IN(); - // const auto & blocked_dofs = this->dof_manager.getGlobalBlockedDOFs(); - - // for (auto && data : arange(blocked_dofs)) { - // auto & dof = std::get<0>(data); - // auto & blocked = std::get<1>(data); - - // if (blocked) - // blocked_list.push_back(dof); - // } + const auto & blocked_dofs = this->dof_manager.getGlobalBlockedDOFs(); + std::vector rows; + for (auto && data : enumerate(blocked)) { + if (std::get<1>(data)) { + rows.push_back(std::get<0>(data)); + } + } - // PETSc_call(MatZeroRowsColumnsLocal, mat, blocked_list.size(), - // blocked_list.storage(), block_val, 0, 0); + PETSc_call(MatZeroRowsColumnsLocal, mat, roew.size(), + rows.storage(), block_val, nullptr, nullptr); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrixPETSc::mul(Real alpha) { PETSc_call(MatScale, mat, alpha); this->release++; } /* -------------------------------------------------------------------------- */ void SparseMatrixPETSc::clear() { PETSc_call(MatZeroEntries, mat); this->release++; } /* -------------------------------------------------------------------------- */ // void SparseMatrixPETSc::clearProfile() { // SparseMatrix::clearProfile(); // PETSc_call(MatSetOption, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE); // PETSc_call(MatSetOption, MAT_NEW_NONZERO_ALLOCATIONS, PETSC_TRUE); // PETSc_call(MatSetOption, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE); // PETSc_call(MatZeroEntries, mat); // //this->value_release++; // } /* -------------------------------------------------------------------------- */ UInt SparseMatrixPETSc::add(UInt i, UInt j) { PETSc_call(MatSetValue, mat, i, j, 0, ADD_VALUES); this->release++; return 0; } /* -------------------------------------------------------------------------- */ void SparseMatrixPETSc::add(UInt i, UInt j, Real val) { PETSc_call(MatSetValue, mat, i, j, val, ADD_VALUES); this->release++; } /* -------------------------------------------------------------------------- */ void SparseMatrixPETSc::addLocal(UInt i, UInt j) { PETSc_call(MatSetValueLocal, mat, i, j, 0, ADD_VALUES); this->release++; } /* -------------------------------------------------------------------------- */ void SparseMatrixPETSc::addLocal(UInt i, UInt j, Real val) { PETSc_call(MatSetValueLocal, mat, i, j, val, ADD_VALUES); this->release++; } /* -------------------------------------------------------------------------- */ void SparseMatrixPETSc::addLocal(Vector & rows, Vector & cols, Matrix & vals) { PETSc_call(MatSetValuesLocal, mat, rows.size(), rows.storage(), cols.size(), cols.storage(), vals.storage(), ADD_VALUES); this->release++; } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/test/test_model/test_model_solver/CMakeLists.txt b/test/test_model/test_model_solver/CMakeLists.txt index 2d365b9bd..4ae04a38a 100644 --- a/test/test_model/test_model_solver/CMakeLists.txt +++ b/test/test_model/test_model_solver/CMakeLists.txt @@ -1,43 +1,50 @@ #=============================================================================== # @file CMakeLists.txt # # @author Nicolas Richart # # @date creation: Fri Sep 03 2010 # @date last modification: Sat Apr 01 2017 # # @brief test for the common solvers interface of the models # # @section LICENSE # # Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # 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 . # #=============================================================================== register_test(test_dof_manager_default SOURCES test_dof_manager_default.cc PACKAGE mumps ) register_test(test_model_solver SOURCES test_model_solver.cc PACKAGE mumps ) register_test(test_model_solver_dynamic_explicit SOURCES test_model_solver_dynamic.cc PACKAGE core COMPILE_OPTIONS "EXPLICIT=true" ) register_test(test_model_solver_dynamic_implicit SOURCES test_model_solver_dynamic.cc PACKAGE mumps COMPILE_OPTIONS "EXPLICIT=false" ) + + +register_test(test_model_solver_dynamic_petsc + SOURCES test_model_solver_dynamic_petsc.cc + PACKAGE petsc + COMPILE_OPTIONS "EXPLICIT=false" + ) diff --git a/test/test_model/test_model_solver/test_model_solver_dynamic_petsc.cc b/test/test_model/test_model_solver/test_model_solver_dynamic_petsc.cc new file mode 100644 index 000000000..41f3b40c4 --- /dev/null +++ b/test/test_model/test_model_solver/test_model_solver_dynamic_petsc.cc @@ -0,0 +1,838 @@ +/** + * @file test_model_solver_dynamic.cc + * + * @author Nicolas Richart + * + * @date creation: Wed Apr 13 2016 + * @date last modification: Tue Feb 20 2018 + * + * @brief Test default dof manager + * + * @section LICENSE + * + * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + * 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 . + * + */ + +/* -------------------------------------------------------------------------- */ +#include "communicator.hh" +#include "element_group.hh" +#include "mesh.hh" +#include "mesh_accessor.hh" +#include "non_linear_solver.hh" +/* -------------------------------------------------------------------------- */ +#include "boundary_condition_functor.hh" +#include "mpi_communicator_data.hh" +/* -------------------------------------------------------------------------- */ +#include "dumpable_inline_impl.hh" +#include "dumper_element_partition.hh" +#include "dumper_iohelper_paraview.hh" +/* -------------------------------------------------------------------------- */ +#include +/* -------------------------------------------------------------------------- */ +#include +#include +#include +/* -------------------------------------------------------------------------- */ + +#ifndef EXPLICIT +#define EXPLICIT true +#endif + +template +void CHECK_ERR_CXX(func && func_, PetscErrorCode ierr) { + if (PetscUnlikely(ierr != 0)) { + const char * desc; + PetscErrorMessage(ierr, &desc, nullptr); + AKANTU_EXCEPTION("Error in PETSc call to \'" << func_ << "\': " << desc); + } +} + +using namespace akantu; + +static void genMesh(Mesh & mesh, UInt nb_nodes); + +class MyModel { +public: + MyModel(Real F, Mesh & mesh, bool lumped) + : nb_dofs(mesh.getNbNodes()), nb_elements(mesh.getNbElement(_segment_2)), + lumped(lumped), E(1.), A(1.), rho(1.), mesh(mesh), + displacement(nb_dofs, 1, "disp"), velocity(nb_dofs, 1, "velo"), + acceleration(nb_dofs, 1, "accel"), blocked(nb_dofs, 1, "blocked"), + forces(nb_dofs, 1, "force_ext"), + internal_forces(nb_dofs, 1, "force_int"), + stresses(nb_elements, 1, "stress"), strains(nb_elements, 1, "strain"), + initial_lengths(nb_elements, 1, "L0") { + + auto n_global = mesh.getNbGlobalNodes(); + int n_local = 0; + + std::vector nodes_global_ids(nb_dofs); + for (auto && data : enumerate(nodes_global_ids)) { + auto n = std::get<0>(data); + n_local += mesh.isLocalOrMasterNode(n); + std::get<1>(data) = mesh.getNodeGlobalId(n); + } + + mpi_comm = dynamic_cast( + mesh.getCommunicator().getCommunicatorData()) + .getMPICommunicator(); + + MeshAccessor mesh_accessor(mesh); + + ierr = ISLocalToGlobalMappingCreate( + mpi_comm, 1, mesh.getNbNodes(), nodes_global_ids.data(), + PETSC_COPY_VALUES, &petsc_local_to_global); + CHECK_ERR_CXX("ISLocalToGlobalMappingCreate", ierr); + + auto setName = [](auto && Obj, auto && name) { + PetscObjectSetName(reinterpret_cast(Obj), name); + }; + + ierr = VecCreate(mpi_comm, &rhs); + ierr = VecSetSizes(rhs, n_local, n_global); + ierr = VecSetFromOptions(rhs); + ierr = VecSetLocalToGlobalMapping(rhs, petsc_local_to_global); + setName(rhs, "rhs"); + + ierr = VecDuplicate(rhs, &x); + ierr = VecDuplicate(rhs, &x_save); + ierr = VecDuplicate(rhs, &dx); + ierr = VecDuplicate(rhs, &f_int); + ierr = VecDuplicate(rhs, &f_dirichlet); + setName(x, "x"); + setName(x_save, "x save"); + setName(dx, "dx"); + setName(f_int, "f_int"); + setName(f_dirichlet, "f_dirichlet"); + + ierr = MatCreate(mpi_comm, &M); + ierr = MatSetSizes(M, n_local, n_local, n_global, n_global); + ierr = MatSetFromOptions(M); + ierr = MatSetOption(M, MAT_SYMMETRIC, PETSC_TRUE); + ierr = MatSetOption(M, MAT_ROW_ORIENTED, PETSC_TRUE); + ierr = MatSetUp(M); + ierr = MatSetLocalToGlobalMapping(M, petsc_local_to_global, + petsc_local_to_global); + setName(M, "M"); + + assembleMass(); + + ierr = MatDuplicate(M, MAT_DO_NOT_COPY_VALUES, &K); + setName(K, "K"); + ierr = MatDuplicate(M, MAT_DO_NOT_COPY_VALUES, &J); + setName(J, "J"); + + ierr = SNESCreate(mpi_comm, &snes); + ierr = SNESSetFromOptions(snes); + ierr = SNESSetFunction(snes, rhs, MyModel::FormFunction, this); + ierr = SNESSetJacobian(snes, J, J, MyModel::FormJacobian, this); + + PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INDEX); + + displacement.set(0.); + velocity.set(0.); + acceleration.set(0.); + + forces.set(0.); + blocked.set(false); + blocked(0, 0) = true; + blocked(nb_dofs - 1, 0) = true; + displacement(0, 0) = 0; + displacement(nb_dofs - 1, 0) = 1; + + for (auto && data : + zip(make_view(this->mesh.getConnectivity(_segment_2), 2), + make_view(this->initial_lengths))) { + const auto & conn = std::get<0>(data); + auto & L = std::get<1>(data); + + auto p1 = this->mesh.getNodes()(conn(0), _x); + auto p2 = this->mesh.getNodes()(conn(1), _x); + + L = std::abs(p2 - p1); + } + } + + // static PetscErrorCode SNESMonitor(SNES snes,PetscInt its,PetscReal + // fnorm,void *ctx) { + // auto & _this = *reinterpret_cast(ctx); + // //SNESMonitorDefault(snes, its, fnorm, PETSC_VIEWER_STDOUT_WORLD); + // } + + static PetscErrorCode FormFunction(SNES /*snes*/, Vec /*dx*/, Vec /*f*/, + void * ctx) { + auto & _this = *reinterpret_cast(ctx); + _this.assembleResidual(); + return 0; + } + + static PetscErrorCode FormJacobian(SNES /*snes*/, Vec /*dx*/, Mat /*J*/, + Mat /*P*/, void * ctx) { + auto & _this = *reinterpret_cast(ctx); + _this.assembleJacobian(); + return 0; + } + + ~MyModel() { + ierr = MatDestroy(&M); + ierr = MatDestroy(&K); + ierr = MatDestroy(&J); + + ierr = VecDestroy(&rhs); + ierr = VecDestroy(&x); + ierr = VecDestroy(&dx); + ierr = VecDestroy(&x_save); + ierr = VecDestroy(&f_int); + + PetscFinalize(); + } + + void solveStep() { + std::cout << "solveStep" << std::endl; + copy(x_save, displacement); + + ierr = SNESSolve(snes, NULL, dx); + CHECK_ERR_CXX("SNESSolve", ierr); + + setSolutionToDisplacement(); + assembleResidual(); + } + + void applyBC() { + std::vector rows; + for (auto && data : enumerate(blocked)) { + if (std::get<1>(data)) { + rows.push_back(std::get<0>(data)); + } + } + + copy(x, displacement); + ierr = MatZeroRowsColumnsLocal(J, rows.size(), rows.data(), 1., x, + f_dirichlet); + VecView(f_dirichlet, PETSC_VIEWER_STDOUT_WORLD); + CHECK_ERR_CXX("MatZeroRowsColumnsLocal", ierr); + } + + void setSolutionToDisplacement() { + std::cout << "setSolutionToDisplacement" << std::endl; + ierr = VecWAXPY(x, 1, x_save, dx); + copy(displacement, x); + } + + void assembleJacobian() { + std::cout << "assembleJacobian" << std::endl; + setSolutionToDisplacement(); + + assembleStiffness(); + + ierr = MatZeroEntries(J); + CHECK_ERR_CXX("MatZeroEntries", ierr); + + ierr = MatAXPY(J, 1., K, SAME_NONZERO_PATTERN); + CHECK_ERR_CXX("MatAXPY", ierr); + + MatView(J, PETSC_VIEWER_STDOUT_WORLD); + applyBC(); + MatView(J, PETSC_VIEWER_STDOUT_WORLD); + } + + void assembleMass() { + std::cout << "assembleMass" << std::endl; + ierr = MatZeroEntries(M); + CHECK_ERR_CXX("MatZeroEntries", ierr); + + Array m_all_el(this->nb_elements, 4); + + Matrix m(2, 2); + m(0, 0) = m(1, 1) = 2; + m(0, 1) = m(1, 0) = 1; + + // under integrated + // m(0, 0) = m(1, 1) = 3./2.; + // m(0, 1) = m(1, 0) = 3./2.; + + // lumping the mass matrix + // m(0, 0) += m(0, 1); + // m(1, 1) += m(1, 0); + // m(0, 1) = m(1, 0) = 0; + + for (auto && data : + zip(make_view(this->mesh.getConnectivity(_segment_2), 2), + make_view(m_all_el, 2, 2))) { + const auto & conn = std::get<0>(data); + auto & m_el = std::get<1>(data); + UInt n1 = conn(0); + UInt n2 = conn(1); + + Real p1 = this->mesh.getNodes()(n1, _x); + Real p2 = this->mesh.getNodes()(n2, _x); + + Real L = std::abs(p2 - p1); + + m_el = m; + m_el *= rho * A * L / 6.; + + Vector conn_int(conn.size()); + for (auto && data : zip(conn_int, conn)) { + std::get<0>(data) = std::get<1>(data); + } + + ierr = MatSetValuesLocal(M, conn_int.size(), conn_int.storage(), + conn_int.size(), conn_int.storage(), m.storage(), + ADD_VALUES); + } + + ierr = MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY); + ierr = MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY); + ierr = MatSetOption(M, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE); + + PetscViewer viewer; + ierr = PetscViewerASCIIOpen(mpi_comm, "M.mtx", &viewer); + PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATRIXMARKET); + ierr = MatView(M, viewer); + PetscViewerPopFormat(viewer); + ierr = PetscViewerDestroy(&viewer); + // this->getDOFManager().assembleElementalMatricesToMatrix( + // "M", "disp", m_all_el, _segment_2); + + is_mass_assembled = true; + } + + // MatrixType getMatrixType(const ID &) { return _symmetric; } + + // void assembleMatrix(const ID & matrix_id) { + // if (matrix_id == "K") { + // if (not is_stiffness_assembled) + // this->assembleStiffness(); + // } else if (matrix_id == "M") { + // if (not is_mass_assembled) + // this->assembleMass(); + // } else if (matrix_id == "C") { + // // pass, no damping matrix + // } else { + // AKANTU_EXCEPTION("This solver does not know what to do with a matrix " + // << matrix_id); + // } + // } + + void assembleLumpedMatrix(const ID & matrix_id) { + std::cout << "assembleLumpedMatrix" << std::endl; + AKANTU_EXCEPTION("This solver does not know what to do with a matrix " + << matrix_id); + } + + void assembleStiffness() { + std::cout << "assembleStiffness" << std::endl; + // SparseMatrix & K = this->getDOFManager().getMatrix("K"); + // K.clear(); + ierr = MatZeroEntries(K); + CHECK_ERR_CXX("MatZeroEntries", ierr); + + Matrix k(2, 2); + k(0, 0) = k(1, 1) = 1; + k(0, 1) = k(1, 0) = -1; + + Array k_all_el(this->nb_elements, 4); + + auto k_it = k_all_el.begin(2, 2); + auto cit = this->mesh.getConnectivity(_segment_2).begin(2); + auto cend = this->mesh.getConnectivity(_segment_2).end(2); + + for (; cit != cend; ++cit, ++k_it) { + const auto & conn = *cit; + UInt n1 = conn(0); + UInt n2 = conn(1); + + Real p1 = this->mesh.getNodes()(n1, _x); + Real p2 = this->mesh.getNodes()(n2, _x); + + Real L = std::abs(p2 - p1); + + auto & k_el = *k_it; + k_el = k; + k_el *= E * A / L; + + Vector conn_int(conn.size()); + for (auto && data : zip(conn_int, conn)) { + std::get<0>(data) = std::get<1>(data); + } + + ierr = MatSetValuesLocal(K, conn_int.size(), conn_int.storage(), + conn_int.size(), conn_int.storage(), + k_el.storage(), ADD_VALUES); + } + + ierr = MatAssemblyBegin(K, MAT_FINAL_ASSEMBLY); + CHECK_ERR_CXX("MatAssemblyBegin", ierr); + + ierr = MatAssemblyEnd(K, MAT_FINAL_ASSEMBLY); + CHECK_ERR_CXX("MatAssemblyEnd", ierr); + + ierr = MatSetOption(K, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE); + CHECK_ERR_CXX("MatSetOption", ierr); + + PetscViewer viewer; + ierr = PetscViewerASCIIOpen(mpi_comm, "K.mtx", &viewer); + CHECK_ERR_CXX("PetscViewerASCIIOpen", ierr); + + PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATRIXMARKET); + ierr = MatView(K, viewer); + CHECK_ERR_CXX("MatView", ierr); + + PetscViewerPopFormat(viewer); + ierr = PetscViewerDestroy(&viewer); + CHECK_ERR_CXX("PetscViewerDestroy", ierr); + + // this->getDOFManager().assembleElementalMatricesToMatrix( + // "K", "disp", k_all_el, _segment_2); + + is_stiffness_assembled = true; + } + + void copy(Array & y, Vec x) { + std::cout << "copy <-" << std::endl; + const PetscScalar * x_local; + ierr = VecGetArrayRead(x, &x_local); + + for (auto && data : zip(y, range(x_local + 0, x_local + y.size()))) { + std::get<0>(data) = std::get<1>(data); + } + ierr = VecRestoreArrayRead(x, &x_local); + + // VecView(x, PETSC_VIEWER_STDOUT_WORLD); + // std::cout << y.getID() << " " << Vector(y.storage(), y.size()) + // << std::endl; + } + + void print(const Array & x) const { + std::cout << x.getID() << " " << Vector(x.storage(), x.size()) + << std::endl; + } + + void copy(Vec x, const Array & y) { + std::cout << "copy ->" << std::endl; + PetscScalar * x_local; + ierr = VecGetArray(x, &x_local); + + for (auto && data : zip(y, range(x_local + 0, x_local + y.size()))) { + std::get<1>(data) = std::get<0>(data); + } + ierr = VecRestoreArray(x, &x_local); + + // std::cout << y.getID() << " " << Vector(y.storage(), y.size()) + // << std::endl; + // VecView(x, PETSC_VIEWER_STDOUT_WORLD); + } + + void assembleResidual() { + std::cout << "assembleResidual" << std::endl; + // this->getDOFManager().assembleToResidual("disp", forces); + setSolutionToDisplacement(); + copy(rhs, forces); + // VecAXPY(rhs, -1., f_dirichlet); + + print(displacement); + this->assembleResidual(_not_ghost); + // this->synchronize(SynchronizationTag::_user_1); + + // this->getDOFManager().assembleToResidual("disp", internal_forces, -1.); + VecAXPY(rhs, 1., f_int); + + for (auto && data : enumerate(blocked)) { + if(std::get<1>(data)) { + VecSetValueLocal(rhs, std::get<0>(data), 0., INSERT_VALUES); + } + } + VecAssemblyBegin(rhs); + VecAssemblyEnd(rhs); + + VecView(rhs, PETSC_VIEWER_STDOUT_WORLD); + } + + void assembleResidual(const GhostType & ghost_type) { + std::cout << "assembleResidual" << std::endl; + VecZeroEntries(f_int); + + auto cit = this->mesh.getConnectivity(_segment_2, ghost_type).begin(2); + auto cend = this->mesh.getConnectivity(_segment_2, ghost_type).end(2); + + auto strain_it = this->strains.begin(); + auto stress_it = this->stresses.begin(); + auto L_it = this->initial_lengths.begin(); + + for (; cit != cend; ++cit, ++strain_it, ++stress_it, ++L_it) { + const auto & conn = *cit; + UInt n1 = conn(0); + UInt n2 = conn(1); + + Real u1 = this->displacement(n1, _x); + Real u2 = this->displacement(n2, _x); + + *strain_it = (u2 - u1) / *L_it; + *stress_it = E * *strain_it; + Real f_n = A * *stress_it; + + std::cout << n1 << "[" << u1 << "]" + << " <-> " << n2 << "[" << u2 << "]" + << " : " << f_n << std::endl; + + ierr = VecSetValueLocal(f_int, n1, -f_n, ADD_VALUES); + ierr = VecSetValueLocal(f_int, n2, f_n, ADD_VALUES); + } + + ierr = VecAssemblyBegin(f_int); + ierr = VecAssemblyEnd(f_int); + // this->getDOFManager().assembleElementalArrayLocalArray( + // forces_internal_el, internal_forces, _segment_2, ghost_type); + } + + Real getPotentialEnergy() { + std::cout << "getPotentialEnergy" << std::endl; + copy(x, displacement); + Vec Ax; + + ierr = VecDuplicate(x, &Ax); + ierr = MatMult(K, x, Ax); + PetscScalar res; + ierr = VecDot(x, Ax, &res); + + return res / 2.; + } + + Real getKineticEnergy() { + std::cout << "getKineticEnergy" << std::endl; + return 0; + } + + // Real getExternalWorkIncrement() { + // Real res = 0; + + // auto it = velocity.begin(); + // auto end = velocity.end(); + // auto if_it = internal_forces.begin(); + // auto ef_it = forces.begin(); + // auto b_it = blocked.begin(); + + // for (UInt node = 0; it != end; ++it, ++if_it, ++ef_it, ++b_it, ++node) { + // if (mesh.isLocalOrMasterNode(node)) + // res += (*b_it ? -*if_it : *ef_it) * *it; + // } + + // mesh.getCommunicator().allReduce(res, SynchronizerOperation::_sum); + + // return res * this->getTimeStep(); + // } + + // void predictor() {} + // void corrector() {} + + // /* ------------------------------------------------------------------------ + // */ UInt getNbData(const Array & elements, + // const SynchronizationTag &) const { + // return elements.size() * sizeof(Real); + // } + + // void packData(CommunicationBuffer & buffer, const Array & + // elements, + // const SynchronizationTag & tag) const { + // if (tag == SynchronizationTag::_user_1) { + // for (const auto & el : elements) { + // buffer << this->stresses(el.element); + // } + // } + // } + + // void unpackData(CommunicationBuffer & buffer, const Array & + // elements, + // const SynchronizationTag & tag) { + // if (tag == SynchronizationTag::_user_1) { + // auto cit = this->mesh.getConnectivity(_segment_2, _ghost).begin(2); + + // for (const auto & el : elements) { + // Real stress; + // buffer >> stress; + + // Real f = A * stress; + + // Vector conn = cit[el.element]; + // this->internal_forces(conn(0), _x) += -f; + // this->internal_forces(conn(1), _x) += f; + // } + // } + // } + + Real getExternalWorkIncrement() { + std::cout << "getExternalWorkIncrement" << std::endl; + return 0.; + } + + template void applyBC(Functor && func, const ID & group_id) { + auto & group = mesh.getElementGroup(group_id).getNodeGroup().getNodes(); + + auto blocked_dofs = make_view(blocked, 1).begin(); + auto disps = make_view(displacement, 1).begin(); + auto poss = make_view(mesh.getNodes(), 1).begin(); + for (auto && node : group) { + auto disp = Vector(disps[node]); + auto pos = Vector(poss[node]); + auto flags = Vector(blocked_dofs[node]); + func(node, flags, disp, pos); + } + } + + const Mesh & getMesh() const { return mesh; } + + UInt getSpatialDimension() const { return 1; } + + auto & getBlockedDOFs() { return blocked; } + + void setTimeStep(Real dt) { + std::cout << "setTimeStep" << std::endl; + this->dt = dt; + } + +private: + PetscErrorCode ierr{0}; + MPI_Comm mpi_comm; + ISLocalToGlobalMapping petsc_local_to_global; + + UInt nb_dofs; + UInt nb_elements; + + bool lumped; + + bool is_stiffness_assembled{false}; + bool is_mass_assembled{false}; + bool is_lumped_mass_assembled{false}; + + Mat K{nullptr}, J{nullptr}, M{nullptr}; + Vec rhs{nullptr}, x{nullptr}, x_save{nullptr}, dx{nullptr}, f_int{nullptr}, + f_dirichlet{nullptr}; + + SNES snes; + + Real dt{0}; + Array save_displacement; + +public: + Real E, A, rho; + + Mesh & mesh; + Array displacement; + Array velocity; + Array acceleration; + Array blocked; + Array forces; + Array internal_forces; + + Array stresses; + Array strains; + + Array initial_lengths; +}; + +/* -------------------------------------------------------------------------- */ +class Sinusoidal : public BC::Dirichlet::DirichletFunctor { +public: + Sinusoidal(MyModel & model, Real amplitude, Real pulse_width, Real t) + : model(model), A(amplitude), k(2 * M_PI / pulse_width), + t(t), v{std::sqrt(model.E / model.rho)} {} + + void operator()(UInt n, Vector & /*flags*/, Vector & disp, + const Vector & coord) const { + auto x = coord(_x); + model.velocity(n, _x) = k * v * A * sin(k * (x - v * t)); + disp(_x) = A * cos(k * (x - v * t)); + } + +private: + MyModel & model; + Real A{1.}; + Real k{2 * M_PI}; + Real t{1.}; + Real v{1.}; +}; + +/* -------------------------------------------------------------------------- */ +int main(int argc, char * argv[]) { + initialize(argc, argv); + + PetscInitialize(&argc, &argv, nullptr, nullptr); + + UInt prank = Communicator::getStaticCommunicator().whoAmI(); + UInt global_nb_nodes = 3; + UInt max_steps = 400; + Real time_step = 0.001; + Mesh mesh(1); + Real F = -9.81; + bool _explicit = EXPLICIT; + const Real pulse_width = 0.2; + const Real A = 0.01; + + if (prank == 0) + genMesh(mesh, global_nb_nodes); + + mesh.distribute(); + + // mesh.makePeriodic(_x); + + MyModel model(F, mesh, _explicit); + + // model.forces.clear(); + // model.blocked.clear(); + + // model.applyBC(Sinusoidal(model, A, pulse_width, 0.), "all"); + // model.applyBC(BC::Dirichlet::FlagOnly(_x), "border"); + + // if (!_explicit) { + // model.getNewSolver("dynamic", TimeStepSolverType::_dynamic, + // NonLinearSolverType::_newton_raphson); + // model.setIntegrationScheme("dynamic", "disp", + // IntegrationSchemeType::_trapezoidal_rule_2, + // IntegrationScheme::_displacement); + // } else { + // model.getNewSolver("dynamic", TimeStepSolverType::_dynamic_lumped, + // NonLinearSolverType::_lumped); + // model.setIntegrationScheme("dynamic", "disp", + // IntegrationSchemeType::_central_difference, + // IntegrationScheme::_acceleration); + // } + + model.setTimeStep(time_step); + + if (prank == 0) { + std::cout << std::scientific; + std::cout << std::setw(14) << "time" + << "," << std::setw(14) << "wext" + << "," << std::setw(14) << "epot" + << "," << std::setw(14) << "ekin" + << "," << std::setw(14) << "total" + << "," << std::setw(14) << "max_disp" + << "," << std::setw(14) << "min_disp" << std::endl; + } + Real wext = 0.; + + // model.getDOFManager().clearResidual(); + // model.assembleResidual(); + + Real epot = 0; // model.getPotentialEnergy(); + Real ekin = 0; // model.getKineticEnergy(); + Real einit = ekin + epot; + Real etot = ekin + epot - wext - einit; + + Real max_disp = 0., min_disp = 0.; + for (auto && disp : model.displacement) { + max_disp = std::max(max_disp, disp); + min_disp = std::min(min_disp, disp); + } + + if (prank == 0) { + std::cout << std::setw(14) << 0. << "," << std::setw(14) << wext << "," + << std::setw(14) << epot << "," << std::setw(14) << ekin << "," + << std::setw(14) << etot << "," << std::setw(14) << max_disp + << "," << std::setw(14) << min_disp << std::endl; + } + + // #if EXPLICIT == false + // NonLinearSolver & solver = + // model.getDOFManager().getNonLinearSolver("dynamic"); + + // solver.set("max_iterations", 20); + // #endif + + auto * dumper = new DumperParaview("dynamic", "./paraview"); + mesh.registerExternalDumper(*dumper, "dynamic", true); + mesh.addDumpMesh(mesh); + + mesh.addDumpFieldExternalToDumper("dynamic", "displacement", + model.displacement); + mesh.addDumpFieldExternalToDumper("dynamic", "velocity", model.velocity); + mesh.addDumpFieldExternalToDumper("dynamic", "forces", model.forces); + mesh.addDumpFieldExternalToDumper("dynamic", "acceleration", + model.acceleration); + + mesh.dump(); + max_steps = 1; + for (UInt i = 1; i < max_steps + 1; ++i) { + // model.applyBC(Sinusoidal(model, A, pulse_width, time_step * (i - 1)), + // "border"); + + model.solveStep(); + mesh.dump(); + + epot = model.getPotentialEnergy(); + ekin = model.getKineticEnergy(); + wext += model.getExternalWorkIncrement(); + etot = ekin + epot - wext - einit; + + Real max_disp = 0., min_disp = 0.; + for (auto && disp : model.displacement) { + max_disp = std::max(max_disp, disp); + min_disp = std::min(min_disp, disp); + } + + if (prank == 0) { + std::cout << std::setw(14) << time_step * i << "," << std::setw(14) + << wext << "," << std::setw(14) << epot << "," << std::setw(14) + << ekin << "," << std::setw(14) << etot << "," << std::setw(14) + << max_disp << "," << std::setw(14) << min_disp << std::endl; + } + } + + // output.close(); + // finalize(); + // PetscFinalize(); + + return EXIT_SUCCESS; +} + +/* -------------------------------------------------------------------------- */ +void genMesh(Mesh & mesh, UInt nb_nodes) { + MeshAccessor mesh_accessor(mesh); + Array & nodes = mesh_accessor.getNodes(); + Array & conn = mesh_accessor.getConnectivity(_segment_2); + + nodes.resize(nb_nodes); + + auto & all = mesh.createNodeGroup("all_nodes"); + + for (UInt n = 0; n < nb_nodes; ++n) { + nodes(n, _x) = n * (1. / (nb_nodes - 1)); + all.add(n); + } + + mesh.createElementGroupFromNodeGroup("all", "all_nodes"); + + conn.resize(nb_nodes - 1); + for (UInt n = 0; n < nb_nodes - 1; ++n) { + conn(n, 0) = n; + conn(n, 1) = n + 1; + } + + Array & conn_points = mesh_accessor.getConnectivity(_point_1); + conn_points.resize(2); + + conn_points(0, 0) = 0; + conn_points(1, 0) = nb_nodes - 1; + + auto & border = mesh.createElementGroup("border", 0); + border.add({_point_1, 0, _not_ghost}, true); + border.add({_point_1, 1, _not_ghost}, true); + + mesh_accessor.makeReady(); +}