Skip to content
Snippets Groups Projects
Commit 2563d289 authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

new argument to ddm solve() to specify if the subdomain matrices don't change... 

new argument to ddm solve() to specify if the subdomain matrices don't change when changing the physical/artificial commutator.

If matrices don't change, always reuse the factorization in solve() - which eliminates 2 factorizations per solve per subdomain.

This should eventually be auto-detected, by checking if the bilinear terms in the volume formulations depend on the physical/artificial commutator.

This patch also changes the handling of verbosity for "inner" solve operations, by setting to verbose-1 (instead of the harcoded value 2)
parent df11a7be
No related branches found
No related tags found
1 merge request!9new argument to ddm solve() to specify if the subdomain matrices don't change...
Pipeline #8985 passed
Stage: gmshfem
Stage: build
Show whitespace changes
Inline Side-by-side
examples/maxwell/main.cpp
+ 1
− 1
View file @ 2563d289
...@@ -348,7 +348,7 @@ int main(int argc, char **argv) ...@@ -348,7 +348,7 @@ int main(int argc, char **argv)
} }
formulation.pre(); formulation.pre();
// Solve the DDM formulation // Solve the DDM formulation
formulation.solve("gmres"); formulation.solve("gmres", 1e-4, 100, false); // FIXME: last argument should be autodetected!
//save field E //save field E
for(int i = 0; i < nDom; ++i) { for(int i = 0; i < nDom; ++i) {
save(E(i), omega(i), "E" + std::to_string(i)); save(E(i), omega(i), "E" + std::to_string(i));
......
src/problem/Formulation.cpp
+ 96
− 73
View file @ 2563d289
...@@ -244,6 +244,9 @@ namespace gmshddm ...@@ -244,6 +244,9 @@ namespace gmshddm
template< class T_Scalar > template< class T_Scalar >
gmshfem::common::Timer Formulation< T_Scalar >::pre() gmshfem::common::Timer Formulation< T_Scalar >::pre()
{ {
int outerVerbosity = gmshfem::common::Options::instance()->verbose;
int innerVerbosity = outerVerbosity - 1;
gmshfem::common::Timer time; gmshfem::common::Timer time;
time.tick(); time.tick();
...@@ -293,8 +296,7 @@ namespace gmshddm ...@@ -293,8 +296,7 @@ namespace gmshddm
_physicalCommutator = true; _physicalCommutator = true;
_artificialCommutator = false; _artificialCommutator = false;
int verbose = gmshfem::common::Options::instance()->verbose; gmshfem::common::Options::instance()->verbose = innerVerbosity;
gmshfem::common::Options::instance()->verbose = 2;
// Pre-process the problem for all subdomains to get the dofs from the neighbours // Pre-process the problem for all subdomains to get the dofs from the neighbours
for(auto idom = 0ULL; idom < _volume.size(); ++idom) { for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
...@@ -612,24 +614,13 @@ namespace gmshddm ...@@ -612,24 +614,13 @@ namespace gmshddm
if(mpi::isItMySubdomain(idom)) { if(mpi::isItMySubdomain(idom)) {
for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) { for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) {
it->second->assemble(); it->second->assemble();
for(auto itTerm = it->second->begin(); itTerm != it->second->end(); ++itTerm) {
if((*itTerm)->isBilinear()) {
(*itTerm)->desactivate();
}
}
}
}
}
for(auto idom = 0ULL; idom < _surface.size(); ++idom) {
if(mpi::isItMySubdomain(idom)) {
for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) {
it->second->solve(); it->second->solve();
} }
} }
} }
gmshfem::common::Options::instance()->verbose = verbose; gmshfem::common::Options::instance()->verbose = outerVerbosity;
// get the local interface unknowns for each processor // get the local interface unknowns for each processor
for(auto idom = 0ULL; idom < _surface.size(); ++idom) { for(auto idom = 0ULL; idom < _surface.size(); ++idom) {
if(mpi::isItMySubdomain(idom)) { if(mpi::isItMySubdomain(idom)) {
...@@ -671,11 +662,18 @@ namespace gmshddm ...@@ -671,11 +662,18 @@ namespace gmshddm
} }
template< class T_Scalar > template< class T_Scalar >
gmshfem::common::Timer Formulation< T_Scalar >::solve(const std::string &solver, const double tolerance, const int iterMax) gmshfem::common::Timer Formulation< T_Scalar >::solve(const std::string &solver, const double tolerance, const int iterMax, const bool sameMatrixWithArtificialAndPhysicalSources)
{ {
// TODO: sameMatrixWithArtificialAndPhysicalSources could be automatically
// detected by examining if bilinear terms of volume formulations use the
// physical/artifical communicator
gmshfem::common::Timer time; gmshfem::common::Timer time;
time.tick(); time.tick();
int outerVerbosity = gmshfem::common::Options::instance()->verbose;
int innerVerbosity = outerVerbosity - 1;
unsigned int MPI_Rank = mpi::getMPIRank(); unsigned int MPI_Rank = mpi::getMPIRank();
if(MPI_Rank == 0) { if(MPI_Rank == 0) {
...@@ -705,33 +703,64 @@ namespace gmshddm ...@@ -705,33 +703,64 @@ namespace gmshddm
MatCreateShell(PETSC_COMM_WORLD, locSize, locSize, globSize, globSize, this, &A); MatCreateShell(PETSC_COMM_WORLD, locSize, locSize, globSize, globSize, this, &A);
MatShellSetOperation(A, MATOP_MULT, (void (*)(void)) & MatVectProduct< T_Scalar >); MatShellSetOperation(A, MATOP_MULT, (void (*)(void)) & MatVectProduct< T_Scalar >);
for(auto idom = 0ULL; idom < _volume.size(); ++idom) { _physicalCommutator = false;
if(mpi::isItMySubdomain(idom)) { _artificialCommutator = true;
_volume[idom]->removeSystem();
_volume[idom]->initSystem(); gmshfem::common::Options::instance()->verbose = innerVerbosity;
auto desactivateBilinear = [=](std::size_t idom)
{
for(auto it = _volume[idom]->begin(); it != _volume[idom]->end(); ++it) {
if((*it)->isBilinear()) {
(*it)->desactivate();
} }
} }
for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) {
for(auto itTerm = it->second->begin(); itTerm != it->second->end(); ++itTerm) {
if((*itTerm)->isBilinear()) {
(*itTerm)->desactivate();
}
}
}
};
_physicalCommutator = false; auto activateBilinear = [=](std::size_t idom)
_artificialCommutator = true; {
for(auto it = _volume[idom]->begin(); it != _volume[idom]->end(); ++it) {
if((*it)->isBilinear()) {
(*it)->activate();
}
}
for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) {
for(auto itTerm = it->second->begin(); itTerm != it->second->end(); ++itTerm) {
if((*itTerm)->isBilinear()) {
(*itTerm)->activate();
}
}
}
};
int verbose = gmshfem::common::Options::instance()->verbose;
gmshfem::common::Options::instance()->verbose = 2;
// Assemble the volumic systems for each subdomain // Assemble the volumic systems for each subdomain
for(auto idom = 0ULL; idom < _volume.size(); ++idom) { for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
if(mpi::isItMySubdomain(idom)) { if(mpi::isItMySubdomain(idom)) {
_volume[idom]->setAttribute("ddm::physicalCommutator", _physicalCommutator); _volume[idom]->setAttribute("ddm::physicalCommutator", _physicalCommutator);
_volume[idom]->setAttribute("ddm::artificialCommutator", _artificialCommutator); _volume[idom]->setAttribute("ddm::artificialCommutator", _artificialCommutator);
_volume[idom]->pre(); if(sameMatrixWithArtificialAndPhysicalSources) {
desactivateBilinear(idom);
_volume[idom]->setRHSToZero();
_volume[idom]->assemble(); _volume[idom]->assemble();
for(auto it = _volume[idom]->begin(); it != _volume[idom]->end(); ++it) {
if((*it)->isBilinear()) {
(*it)->desactivate();
} }
else{
_volume[idom]->removeSystem();
_volume[idom]->initSystem();
_volume[idom]->pre();
_volume[idom]->assemble();
desactivateBilinear(idom);
} }
} }
} }
gmshfem::common::Options::instance()->verbose = verbose;
gmshfem::common::Options::instance()->verbose = outerVerbosity;
if(solver == "jacobi") { if(solver == "jacobi") {
if(MPI_Rank == 0) { if(MPI_Rank == 0) {
...@@ -779,58 +808,51 @@ namespace gmshddm ...@@ -779,58 +808,51 @@ namespace gmshddm
// ****************************************************** // ******************************************************
// Final solution // Final solution
// ****************************************************** // ******************************************************
gmshfem::common::Options::instance()->verbose = innerVerbosity;
_physicalCommutator = true;
_artificialCommutator = true;
for(auto idom = 0ULL; idom < _volume.size(); ++idom) { for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
if(mpi::isItMySubdomain(idom)) { if(mpi::isItMySubdomain(idom)) {
_volume[idom]->setAttribute("ddm::physicalCommutator", _physicalCommutator);
_volume[idom]->setAttribute("ddm::artificialCommutator", _artificialCommutator);
if(sameMatrixWithArtificialAndPhysicalSources) {
_volume[idom]->setRHSToZero();
}
else {
_volume[idom]->removeSystem(); _volume[idom]->removeSystem();
_volume[idom]->initSystem(); _volume[idom]->initSystem();
for(auto it = _volume[idom]->begin(); it != _volume[idom]->end(); ++it) { _volume[idom]->pre();
if((*it)->isBilinear()) {
(*it)->activate();
}
} }
} }
} }
_physicalCommutator = true;
_artificialCommutator = true;
const PetscScalar *petscArray; const PetscScalar *petscArray;
PetscInt petscArraySize; PetscInt petscArraySize;
VecGetArrayRead(Y, &petscArray); VecGetArrayRead(Y, &petscArray);
VecGetSize(Y, &petscArraySize); VecGetSize(Y, &petscArraySize);
verbose = gmshfem::common::Options::instance()->verbose;
gmshfem::common::Options::instance()->verbose = 2;
for(auto idom = 0ULL; idom < _surface.size(); ++idom) {
if(mpi::isItMySubdomain(idom)) {
_volume[idom]->setAttribute("ddm::physicalCommutator", _physicalCommutator);
_volume[idom]->setAttribute("ddm::artificialCommutator", _artificialCommutator);
_volume[idom]->pre();
for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) {
for(auto itTerm = it->second->begin(); itTerm != it->second->end(); ++itTerm) {
if((*itTerm)->isBilinear()) {
(*itTerm)->activate();
}
}
}
}
}
gmshfem::common::Options::instance()->verbose = verbose;
const T_Scalar *array = PetscInterface< T_Scalar, PetscScalar, PetscInt >::arrayInterface(petscArray, petscArraySize); const T_Scalar *array = PetscInterface< T_Scalar, PetscScalar, PetscInt >::arrayInterface(petscArray, petscArraySize);
_fillG(array); _fillG(array);
PetscInterface< T_Scalar, PetscScalar, PetscInt >::freeArray(array); PetscInterface< T_Scalar, PetscScalar, PetscInt >::freeArray(array);
VecRestoreArrayRead(Y, &petscArray); VecRestoreArrayRead(Y, &petscArray);
verbose = gmshfem::common::Options::instance()->verbose;
gmshfem::common::Options::instance()->verbose = 2;
// SolveVolumePDE: // SolveVolumePDE:
for(auto idom = 0ULL; idom < _volume.size(); ++idom) { for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
if(mpi::isItMySubdomain(idom)) { if(mpi::isItMySubdomain(idom)) {
if(!sameMatrixWithArtificialAndPhysicalSources) {
activateBilinear(idom);
}
_volume[idom]->assemble(); _volume[idom]->assemble();
_volume[idom]->solve(); _volume[idom]->solve(true);
if(sameMatrixWithArtificialAndPhysicalSources) {
activateBilinear(idom);
}
} }
} }
gmshfem::common::Options::instance()->verbose = verbose;
gmshfem::common::Options::instance()->verbose = outerVerbosity;
MatDestroy(&A); MatDestroy(&A);
VecDestroy(&X); VecDestroy(&X);
...@@ -869,14 +891,11 @@ namespace gmshddm ...@@ -869,14 +891,11 @@ namespace gmshddm
_physicalCommutator = false; _physicalCommutator = false;
_artificialCommutator = true; _artificialCommutator = true;
int verbose = gmshfem::common::Options::instance()->verbose;
gmshfem::common::Options::instance()->verbose = 2;
for(auto idom = 0ULL; idom < _volume.size(); ++idom) { for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
_volume[idom]->setAttribute("ddm::physicalCommutator", _physicalCommutator); _volume[idom]->setAttribute("ddm::physicalCommutator", _physicalCommutator);
_volume[idom]->setAttribute("ddm::artificialCommutator", _artificialCommutator); _volume[idom]->setAttribute("ddm::artificialCommutator", _artificialCommutator);
_volume[idom]->pre(); _volume[idom]->pre();
} }
gmshfem::common::Options::instance()->verbose = verbose;
_IA = true; _IA = true;
...@@ -936,6 +955,9 @@ namespace gmshddm ...@@ -936,6 +955,9 @@ namespace gmshddm
template< class T_Scalar > template< class T_Scalar >
static int MatVectProduct(Mat A, Vec X, Vec Y) static int MatVectProduct(Mat A, Vec X, Vec Y)
{ {
int outerVerbosity = gmshfem::common::Options::instance()->verbose;
int innerVerbosity = outerVerbosity - 1;
Formulation< T_Scalar > *formulation; Formulation< T_Scalar > *formulation;
MatShellGetContext(A, (void **)&formulation); MatShellGetContext(A, (void **)&formulation);
...@@ -948,8 +970,8 @@ namespace gmshddm ...@@ -948,8 +970,8 @@ namespace gmshddm
Formulation< T_Scalar >::template PetscInterface< T_Scalar, PetscScalar, PetscInt >::freeArray(array); Formulation< T_Scalar >::template PetscInterface< T_Scalar, PetscScalar, PetscInt >::freeArray(array);
VecRestoreArrayRead(X, &petscArray); VecRestoreArrayRead(X, &petscArray);
int verbose = gmshfem::common::Options::instance()->verbose; gmshfem::common::Options::instance()->verbose = innerVerbosity;
gmshfem::common::Options::instance()->verbose = 2;
// SolveVolumePDE // SolveVolumePDE
for(auto idom = 0ULL; idom < formulation->_volume.size(); ++idom) { for(auto idom = 0ULL; idom < formulation->_volume.size(); ++idom) {
if(mpi::isItMySubdomain(idom)) { if(mpi::isItMySubdomain(idom)) {
...@@ -982,7 +1004,8 @@ namespace gmshddm ...@@ -982,7 +1004,8 @@ namespace gmshddm
} }
} }
gmshfem::common::Options::instance()->verbose = verbose; gmshfem::common::Options::instance()->verbose = outerVerbosity;
gmshfem::algebra::Vector< T_Scalar > g_local; gmshfem::algebra::Vector< T_Scalar > g_local;
for(auto idom = 0ULL; idom < formulation->_surface.size(); ++idom) { for(auto idom = 0ULL; idom < formulation->_surface.size(); ++idom) {
if(mpi::isItMySubdomain(idom)) { if(mpi::isItMySubdomain(idom)) {
......
src/problem/Formulation.h
+ 7
− 7
View file @ 2563d289
...@@ -60,7 +60,7 @@ namespace gmshddm ...@@ -60,7 +60,7 @@ namespace gmshddm
gmshfem::common::Timer pre(); gmshfem::common::Timer pre();
gmshfem::common::Timer solve(const std::string &solver, const double tolerance = 1e-6, const int iterMax = 1000); gmshfem::common::Timer solve(const std::string &solver, const double tolerance = 1e-6, const int iterMax = 1000, const bool sameMatrixWithArtificialAndPhysicalSources = false);
gmshfem::algebra::MatrixCCS< T_Scalar > computeMatrix(); gmshfem::algebra::MatrixCCS< T_Scalar > computeMatrix();
const gmshfem::algebra::Vector< T_Scalar > &getRHS() const; const gmshfem::algebra::Vector< T_Scalar > &getRHS() const;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment