Select Git revision
Cal_SmallFemTermOfFemEquation.h
Formulation.cpp 59.76 KiB
// GmshDDM - Copyright (C) 2019-2022, A. Royer, C. Geuzaine, Université de Liège
//
// See the LICENSE.txt file for license information. Please report all
// issues on https://gitlab.onelab.info/gmsh/ddm/issues
#include "Formulation.h"
#include "GmshDdm.h"
#include "MPIInterface.h"
#include "Options.h"
#include "gmshddmDefines.h"
#include "resources.h"
#include <gmshfem/CSVio.h>
#include <gmshfem/Dof.h>
#include <gmshfem/Message.h>
#include <gmshfem/Vector.h>
#ifdef HAVE_MPI
#include <mpi.h>
#endif
#ifdef HAVE_PETSC
#include <petsc.h>
#endif
namespace gmshddm
{
namespace problem
{
static void s_printResources(const std::string &step, bool first = false)
{
static gmshfem::common::Timer timeRef;
if(first) timeRef.tick();
gmshfem::common::Timer time = timeRef;
time.tock();
double cpuTime, peakRSS, currentRSS;
common::getResources(cpuTime, peakRSS, currentRSS);
if(mpi::getMPISize() > 1) {
double val[3] = {cpuTime, peakRSS, currentRSS};
double min[3] = {val[0], val[1], val[2]};
double max[3] = {val[0], val[1], val[2]};
MPI_Reduce(val, min, 3, MPI_DOUBLE, MPI_MIN, 0, MPI_COMM_WORLD);
MPI_Reduce(val, max, 3, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD);
if(mpi::getMPIRank() == 0) {
if(gmshddm::common::Options::instance()->memory) {
gmshfem::msg::info << "* Resources @ " << step << ": Wall = " << time << "s, "
<< "CPU = [" << min[0] << "s, " << max[0] << "s], "
<< "PeakRSS = [" << min[1] << "Mb, " << max[1] << "Mb], "
<< "CurrentRSS = [" << min[2] << "Mb, " << max[2] << "Mb]"
<< gmshfem::msg::endl;
}
}
MPI_Barrier(PETSC_COMM_WORLD);
}
else {
if(gmshddm::common::Options::instance()->memory) {
gmshfem::msg::info << "* Resources @ " << step << ": Wall = " << time << "s, "
<< "CPU = " << cpuTime << "s, "
<< "PeakRSS = " << peakRSS << "Mb, "
<< "CurrentRSS = " << currentRSS << "Mb"
<< gmshfem::msg::endl;
}
}
}
template< class T_Scalar >
Formulation< T_Scalar >::Formulation(const std::string &name, const std::vector< std::vector< unsigned int > > &topology) :
_name(name), _volume(), _surface(topology.size()), _interfaceFields(), _interfaceMapping(), _rhs(), _absoluteResidual(), _relativeResidual(), _numberOfIterations(0), _physicalCommutator(false), _artificialCommutator(false), _IA(false)
{
for(auto i = 0ULL; i < topology.size(); ++i) {
gmshfem::msg::debug << "Create volume formulation (" << i << ")." << gmshfem::msg::endl;
_volume.push_back(new gmshfem::problem::Formulation< T_Scalar >("volume_" + std::to_string(i)));
for(auto ii = 0ULL; ii < topology[i].size(); ++ii) {
const unsigned int j = topology[i][ii];
gmshfem::msg::debug << "Create surface formulation (" << i << ", " << j << ")." << gmshfem::msg::endl;
_surface[i].emplace(j, new gmshfem::problem::Formulation< T_Scalar >("surface_" + std::to_string(i) + "_" + std::to_string(j)));
}
}
}
template< class T_Scalar >
Formulation< T_Scalar >::Formulation(const std::string &name, const std::vector< gmshfem::problem::Formulation< T_Scalar > * > &volume, const std::vector< std::unordered_map< unsigned int, gmshfem::problem::Formulation< T_Scalar > * > > &surface) :
_name(name), _volume(volume), _surface(surface), _interfaceFields(), _interfaceMapping(), _rhs(), _absoluteResidual(), _relativeResidual(), _numberOfIterations(0), _physicalCommutator(false), _artificialCommutator(false), _IA(false)
{
}
template< class T_Scalar >
Formulation< T_Scalar >::~Formulation()
{
for(auto i = 0ULL; i < _volume.size(); ++i) {
delete _volume[i];
}
for(auto i = 0ULL; i < _surface.size(); ++i) {
for(auto it = _surface[i].begin(); it != _surface[i].end(); ++it) {
delete it->second;
}
}
}
template< class T_Scalar >
void Formulation< T_Scalar >::_runIterationOperations(const unsigned int iteration, const gmshfem::scalar::Precision< T_Scalar > rnorm)
{
_absoluteResidual.push_back(rnorm);
_relativeResidual.push_back(rnorm / _absoluteResidual[0]);
_numberOfIterations = iteration;
const unsigned int MPI_Rank = mpi::getMPIRank();
if(MPI_Rank == 0) {
gmshfem::msg::info << " - Iteration " << iteration << ": absolute residual = " << _absoluteResidual.back() << ", relative residual = " << _relativeResidual.back() << gmshfem::msg::endl;
}
s_printResources("DDM Iteration");
}
template< class T_Scalar >
gmshfem::algebra::Vector< T_Scalar > Formulation< T_Scalar >::_sortDofValues(const gmshfem::algebra::Vector< T_Scalar > &vector, const unsigned int fieldTag)
{
auto it = _interfaceMapping.find(fieldTag);
if(it == _interfaceMapping.end()) {
gmshfem::msg::debug << "There is any mapping for the field tagged '" << fieldTag << "'" << gmshfem::msg::endl;
return vector;
}
gmshfem::algebra::Vector< T_Scalar > newVector(vector.size());
#pragma omp parallel for
for(auto i = 0ULL; i < vector.size(); ++i) {
newVector[i] = vector[it->second[i]];
}
return newVector;
}
template< class T_Scalar >
void Formulation< T_Scalar >::_fillG(const T_Scalar *array)
{
const unsigned int MPI_Rank = mpi::getMPIRank();
const unsigned int MPI_Size = mpi::getMPISize();
if(MPI_Size == 1) {
unsigned int pos = 0;
for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) {
for(auto itField = _interfaceFields.begin(); itField != _interfaceFields.end(); ++itField) {
if((*itField->first).isDefined(idom, it->first)) {
const unsigned long long size = (*itField->first)(idom, it->first).numberOfUnknownDofs();
gmshfem::algebra::Vector< T_Scalar > g_ij(std::vector< T_Scalar >(&array[pos], &array[pos + size]));
(*itField->first)(idom, it->first).setUnknownVector(g_ij, gmshfem::dofs::RawOrder::Hash);
pos += size;
}
}
}
}
}
else {
std::vector< std::vector< std::vector< gmshfem::algebra::Vector< T_Scalar > * > > > g_ji;
std::vector< MPI_Request * > requests;
// Send and recv g_ij
for(auto recv : {true, false}) {
unsigned int pos = 0;
unsigned int currentSubdomain = 0;
unsigned int currentInterface = 0;
unsigned int currentField = 0;
for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
if(mpi::isItMySubdomain(idom)) {
if(recv) {
g_ji.push_back(std::vector< std::vector< gmshfem::algebra::Vector< T_Scalar > * > >());
}
for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) {
if(mpi::isItMySubdomain(it->first)) {
for(auto itField = _interfaceFields.begin(); itField != _interfaceFields.end(); ++itField) {
if((*itField->first).isDefined(idom, it->first)) {
const unsigned long long size = (*itField->first)(idom, it->first).numberOfUnknownDofs();
if(!recv) {
gmshfem::algebra::Vector< T_Scalar > g_ij(std::vector< T_Scalar >(&array[pos], &array[pos + size]));
(*itField->first)(idom, it->first).setUnknownVector(g_ij, gmshfem::dofs::RawOrder::Hash);
}
pos += size;
}
}
}
else {
if(recv) {
g_ji[currentSubdomain].push_back(std::vector< gmshfem::algebra::Vector< T_Scalar > * >());
}
for(auto itField = _interfaceFields.begin(); itField != _interfaceFields.end(); ++itField) {
if(!(*itField->first).isDefined(idom, it->first)) continue;
const unsigned int tagIJ = (*itField->first)(idom, it->first).tag();
const unsigned int tagJI = (*itField->second)(it->first, idom).tag();
const unsigned long long size = (*itField->first)(idom, it->first).numberOfUnknownDofs();
if(recv) {
g_ji[currentSubdomain][currentInterface].push_back(new gmshfem::algebra::Vector< T_Scalar >(size));
if(size) {
requests.push_back(new MPI_Request);
gmshfem::msg::debug << "MPI_Irecv g rank=" << MPI_Rank << " source=" << mpi::rankOfSubdomain(it->first) << " tag=" << tagJI << " size=" << size << gmshfem::msg::endl;
MPI_Irecv(&(*g_ji[currentSubdomain][currentInterface][currentField])[0], size, MPI_C_DOUBLE_COMPLEX, mpi::rankOfSubdomain(it->first), tagJI, MPI_COMM_WORLD, requests.back());
}
}
else {
gmshfem::algebra::Vector< T_Scalar > g_ij(std::vector< T_Scalar >(&array[pos], &array[pos + size]));
(*itField->first)(idom, it->first).setUnknownVector(g_ij, gmshfem::dofs::RawOrder::Hash);
if(size) {
gmshfem::msg::debug << "MPI_Send g rank=" << MPI_Rank << " dest=" << mpi::rankOfSubdomain(it->first) << " tag=" << tagIJ << " size=" << size << gmshfem::msg::endl;
MPI_Send(&g_ij[0], size, MPI_C_DOUBLE_COMPLEX, mpi::rankOfSubdomain(it->first), tagIJ, MPI_COMM_WORLD);
}
}
pos += size;
currentField += 1;
}
currentField = 0;
currentInterface += 1;
}
}
currentInterface = 0;
currentSubdomain += 1;
}
}
}
// wait until we receive the data
gmshfem::msg::debug << "MPI_Wait g starting on rank=" << MPI_Rank << gmshfem::msg::endl;
for(auto r : requests) {
MPI_Wait(r, MPI_STATUS_IGNORE);
delete r;
}
gmshfem::msg::debug << "MPI_Wait g done on rank=" << MPI_Rank << gmshfem::msg::endl;
MPI_Barrier(PETSC_COMM_WORLD);
{
// assign g_ji
unsigned int currentSubdomain = 0;
unsigned int currentInterface = 0;
unsigned int currentField = 0;
for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
if(mpi::isItMySubdomain(idom)) {
for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) {
if(mpi::isItMySubdomain(it->first)) continue;
for(auto itField = _interfaceFields.begin(); itField != _interfaceFields.end(); ++itField) {
if(!(*itField->first).isDefined(idom, it->first)) continue;
const unsigned int tagIJ = (*itField->first)(idom, it->first).tag();
(*itField->second)(it->first, idom).setUnknownVector(_sortDofValues(*g_ji[currentSubdomain][currentInterface][currentField], tagIJ), gmshfem::dofs::RawOrder::Hash);
delete g_ji[currentSubdomain][currentInterface][currentField];
currentField += 1;
}
currentField = 0;
currentInterface += 1;
}
currentInterface = 0;
currentSubdomain += 1;
}
}
}
}
}
template< class T_Scalar >
unsigned int Formulation< T_Scalar >::size() const
{
return _volume.size();
}
template< class T_Scalar >
unsigned int Formulation< T_Scalar >::size(unsigned int i) const
{
return _surface.at(i).size();
}
template< class T_Scalar >
const gmshfem::function::ScalarFunction< T_Scalar > Formulation< T_Scalar >::physicalSource(const gmshfem::function::ScalarFunction< T_Scalar > &f)
{
return commutator(f, &_physicalCommutator);
}
template< class T_Scalar >
const gmshfem::function::VectorFunction< T_Scalar > Formulation< T_Scalar >::physicalSource(const gmshfem::function::VectorFunction< T_Scalar > &f)
{
return commutator(f, &_physicalCommutator);
}
template< class T_Scalar >
void Formulation< T_Scalar >::physicalSourceTerm(const unsigned int termTag)
{
_physicalSourceTerms.insert(termTag);
}
template< class T_Scalar >
const gmshfem::function::ScalarFunction< T_Scalar > Formulation< T_Scalar >::artificialSource(const gmshfem::function::ScalarFunction< T_Scalar > &f)
{
return commutator(f, &_artificialCommutator);
}
template< class T_Scalar >
const gmshfem::function::VectorFunction< T_Scalar > Formulation< T_Scalar >::artificialSource(const gmshfem::function::VectorFunction< T_Scalar > &f)
{
return commutator(f, &_artificialCommutator);
}
template< class T_Scalar >
void Formulation< T_Scalar >::artificialSourceTerm(const unsigned int termTag)
{
_artificialSourceTerms.insert(termTag);
}
template< class T_Scalar >
void Formulation< T_Scalar >::togglePhysicalAndArtificialSourceTerms()
{
for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
if(mpi::isItMySubdomain(idom)) {
for(auto it = _volume[idom]->begin(); it != _volume[idom]->end(); ++it) {
if(_physicalSourceTerms.find((*it)->tag()) != _physicalSourceTerms.end()) {
if(_physicalCommutator) {
(*it)->activate();
}
else {
(*it)->deactivate();
}
}
else if(_artificialSourceTerms.find((*it)->tag()) != _artificialSourceTerms.end()) {
if(_artificialCommutator) {
(*it)->activate();
}
else {
(*it)->deactivate();
}
}
}
for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) {
for(auto itTerm = it->second->begin(); itTerm != it->second->end(); ++itTerm) {
if(_physicalSourceTerms.find((*itTerm)->tag()) != _physicalSourceTerms.end()) {
if(_physicalCommutator) {
(*itTerm)->activate();
}
else {
(*itTerm)->deactivate();
}
}
else if(_artificialSourceTerms.find((*itTerm)->tag()) != _artificialSourceTerms.end()) {
if(_artificialCommutator) {
(*itTerm)->activate();
}
else {
(*itTerm)->deactivate();
}
}
}
}
}
}
}
template< class T_Scalar >
void Formulation< T_Scalar >::addInterfaceField(field::InterfaceFieldInterface< T_Scalar > &interfaceField)
{
_interfaceFields.insert(std::make_pair(&interfaceField, nullptr));
gmshfem::msg::debug << "New InterfaceField named " << interfaceField.name() << " is added to the formulation '" << _name << "'." << gmshfem::msg::endl;
}
template< class T_Scalar >
field::InterfaceFieldInterface< T_Scalar > *Formulation< T_Scalar >::getInterfaceField(const std::string &name) const
{
for(auto it = _interfaceFields.begin(); it != _interfaceFields.end(); ++it) {
if(it->first->name() == name) {
return it->first;
}
}
return nullptr;
}
template< class T_Scalar >
gmshfem::problem::Formulation< T_Scalar > &Formulation< T_Scalar >::operator()(unsigned int i)
{
return *_volume[i];
}
template< class T_Scalar >
gmshfem::problem::Formulation< T_Scalar > &Formulation< T_Scalar >::operator()(unsigned int i, unsigned int j)
{
return *_surface[i].at(j);
}
template< class T_Scalar >
gmshfem::common::Timer Formulation< T_Scalar >::pre()
{
const int outerVerbosity = gmshfem::common::Options::instance()->verbose;
const int innerVerbosity = (outerVerbosity != 0 ? outerVerbosity - 1 : 0);
const unsigned int MPI_Rank = mpi::getMPIRank();
const unsigned int MPI_Size = mpi::getMPISize();
MPI_Barrier(PETSC_COMM_WORLD);
gmshfem::common::Timer time;
time.tick();
if(MPI_Rank == 0) {
gmshfem::msg::info << "Pre-processing " << _name << "..." << gmshfem::msg::endl;
}
if(_interfaceFields.size() == 0) {
gmshfem::msg::warning << "There is no interface field defined. Perhaps you have forgotten to declare them with 'addInterfaceField' function." << gmshfem::msg::endl;
}
s_printResources("DDM pre-processing begin", true);
// Build interface field mapping
if(MPI_Size != 1) {
const unsigned int nDom = _volume.size();
for(auto it = _interfaceFields.begin(); it != _interfaceFields.end(); ++it) {
for(auto i = 0U; i < nDom; ++i) {
for(auto j = 0U; j < nDom; ++j) {
if(!it->first->isDefined(i, j)) continue;
for(auto it2 = _interfaceFields.begin(); it2 != _interfaceFields.end(); ++it2) {
if(!it2->first->isDefined(j, i)) continue;
if((*it2->first)(j, i).domain() == (*it->first)(i, j).domain()) {
it->second = it2->first;
gmshfem::msg::debug << "Interface field '" << it->first->name() << "' is mapped with interface field '" << it2->first->name() << "'" << gmshfem::msg::endl;
break;
}
}
if(it->second != nullptr) break;
}
if(it->second != nullptr) break;
}
if(it->second == nullptr) {
throw gmshfem::common::Exception("Impossible to match the interface field '" + it->first->name() + "' with another interface field");
}
}
}
// ******************************************************
// Initialize Kyrlov solver (compute b of (I-A) g^n = b)
// ******************************************************
_physicalCommutator = true;
_artificialCommutator = false;
togglePhysicalAndArtificialSourceTerms();
// Pre-process the problem for all subdomains to get the dofs from the neighbours
for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
if(mpi::isItMySubdomain(idom)) {
gmshfem::common::Options::instance()->verbose = innerVerbosity;
_volume[idom]->pre();
gmshfem::common::Options::instance()->verbose = outerVerbosity;
}
}
if(MPI_Size == 1) {
for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) {
gmshfem::common::Options::instance()->verbose = innerVerbosity;
it->second->pre();
gmshfem::common::Options::instance()->verbose = outerVerbosity;
}
}
}
else {
// Create MPI structure
MPI_Aint struct_data[2] = {offsetof(gmshfem::dofs::RawDofKey, numType),
offsetof(gmshfem::dofs::RawDofKey, entity)};
int block_lengths[2] = {1, 1};
MPI_Datatype types[2] = {MPI_UNSIGNED_LONG_LONG, MPI_UNSIGNED_LONG_LONG};
// Commit the new datatype
MPI_Datatype MPI_DOF;
MPI_Type_create_struct(2, block_lengths, struct_data, types, &MPI_DOF);
MPI_Type_commit(&MPI_DOF);
// Send the RawDofKey vector
// If
//
// 1) the mesh on the interface is identical on both subdomains (same
// entities, same elements and same nodes, with the same tags, and
// stored in the same order)
//
// 2) for Dofs associated with entities other than nodes (edges and
// faces), the edges and face tags correspond on both subdomains
// (e.g. they have been prealably generated globally in Gmsh)
//
// then exchanging the RawDofKey vector is actually enough, provided
// that the Dofs are correctly sorted afterwards (see _interfaceMapping
// below)
//
// Since generating the edge and face tags beforehand is cumbersome
// (e.g. if the mesh was not generated with Gmsh), the code below also
// exchanges the edge and face tags, based only on the assumption that
// the meshes match (same entities, same elements and same nodes, with
// the same tags, and stored in the same order).
std::vector< MPI_Request * > requests;
std::vector< std::vector< std::vector< std::vector< gmshfem::dofs::RawDofKey > * > > > dofKeys;
for(auto recv : {true, false}) {
unsigned int currentSubdomain = 0;
unsigned int currentInterface = 0;
unsigned int currentField = 0;
for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
if(mpi::isItMySubdomain(idom)) {
if(recv) {
dofKeys.push_back(std::vector< std::vector< std::vector< gmshfem::dofs::RawDofKey > * > >());
}
for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) {
if(recv) {
gmshfem::common::Options::instance()->verbose = innerVerbosity;
it->second->pre();
gmshfem::common::Options::instance()->verbose = outerVerbosity;
}
if(mpi::isItMySubdomain(it->first)) continue;
if(recv) {
dofKeys[currentSubdomain].push_back(std::vector< std::vector< gmshfem::dofs::RawDofKey > * >());
}
for(auto itField = _interfaceFields.begin(); itField != _interfaceFields.end(); ++itField) {
if(!(*itField->first).isDefined(idom, it->first)) continue;
if(recv) {
dofKeys[currentSubdomain][currentInterface].push_back(new std::vector< gmshfem::dofs::RawDofKey >());
}
const unsigned int tagIJ = (*itField->first)(idom, it->first).tag();
const unsigned int tagJI = (*itField->second)(it->first, idom).tag();
std::vector< gmshfem::dofs::RawDof > dofs;
// get dofs
(*itField->first)(idom, it->first).getUnknownDofVector(dofs);
if(recv) {
dofKeys[currentSubdomain][currentInterface][currentField]->resize(dofs.size());
if(dofs.size()) {
requests.push_back(new MPI_Request);
gmshfem::msg::debug << "MPI_Irecv dofs rank=" << MPI_Rank << " source=" << mpi::rankOfSubdomain(it->first) << " tag=" << tagJI << " size=" << dofs.size() << gmshfem::msg::endl;
MPI_Irecv(&(*dofKeys[currentSubdomain][currentInterface][currentField])[0], dofs.size(), MPI_DOF, mpi::rankOfSubdomain(it->first), tagJI, MPI_COMM_WORLD, requests.back());
}
}
else {
std::vector< gmshfem::dofs::RawDofKey > myDofKeys(dofs.size());
#pragma omp parallel for
for(auto i = 0ULL; i < dofs.size(); ++i) {
myDofKeys[i].numType = dofs[i].numType;
myDofKeys[i].entity = dofs[i].entity;
dofs[i].numType -= GMSHFEM_DOF_FIELD_OFFSET * tagIJ;
dofs[i].numType += GMSHFEM_DOF_FIELD_OFFSET * tagJI;
}
// set dofs to the associated field
(*itField->second)(it->first, idom).setUnknownDofVector(dofs);
if(dofs.size()) {
gmshfem::msg::debug << "MPI_Send dofs rank=" << MPI_Rank << " dest=" << mpi::rankOfSubdomain(it->first) << " tag=" << tagIJ << " size=" << dofs.size() << gmshfem::msg::endl;
MPI_Send(&myDofKeys[0], dofs.size(), MPI_DOF, mpi::rankOfSubdomain(it->first), tagIJ, MPI_COMM_WORLD);
}
}
currentField += 1;
}
currentField = 0;
currentInterface += 1;
}
currentInterface = 0;
currentSubdomain += 1;
}
}
}
// wait until we receive the data
gmshfem::msg::debug << "MPI_Wait dofs starting on rank=" << MPI_Rank << gmshfem::msg::endl;
for(auto r : requests) {
MPI_Wait(r, MPI_STATUS_IGNORE);
delete r;
}
gmshfem::msg::debug << "MPI_Wait dofs done on rank=" << MPI_Rank << gmshfem::msg::endl;
MPI_Barrier(PETSC_COMM_WORLD);
// Generate the EdgeTags and FaceTags vectors
std::vector< std::vector< std::vector< std::vector< unsigned long long > > > > edgeTags, edgeTagsNeighbor, faceTags, faceTagsNeighbor;
{
// This is only useful for Dofs associated to edges, in case a global
// numbering of edges is not available. We might want to detect when
// this communication step is actually necessary.
unsigned int currentSubdomain = 0;
unsigned int currentInterface = 0;
unsigned int currentField = 0;
for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
if(mpi::isItMySubdomain(idom)) {
edgeTags.push_back(std::vector< std::vector< std::vector< unsigned long long > > >());
edgeTagsNeighbor.push_back(std::vector< std::vector< std::vector< unsigned long long > > >());
faceTags.push_back(std::vector< std::vector< std::vector< unsigned long long > > >());
faceTagsNeighbor.push_back(std::vector< std::vector< std::vector< unsigned long long > > >());
for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) {
if(mpi::isItMySubdomain(it->first)) continue;
edgeTags[currentSubdomain].push_back(std::vector< std::vector< unsigned long long > >());
edgeTagsNeighbor[currentSubdomain].push_back(std::vector< std::vector< unsigned long long > >());
faceTags[currentSubdomain].push_back(std::vector< std::vector< unsigned long long > >());
faceTagsNeighbor[currentSubdomain].push_back(std::vector< std::vector< unsigned long long > >());
for(auto itField = _interfaceFields.begin(); itField != _interfaceFields.end(); ++itField) {
if(!(*itField->first).isDefined(idom, it->first)) continue;
edgeTags[currentSubdomain][currentInterface].push_back(std::vector< unsigned long long >());
edgeTagsNeighbor[currentSubdomain][currentInterface].push_back(std::vector< unsigned long long >());
faceTags[currentSubdomain][currentInterface].push_back(std::vector< unsigned long long >());
faceTagsNeighbor[currentSubdomain][currentInterface].push_back(std::vector< unsigned long long >());
// build the interface map
const gmshfem::domain::Domain domain = (*itField->first)(idom, it->first).domain();
for(auto itEntity = domain.cbegin(); itEntity != domain.cend(); ++itEntity) {
// edge case
if((*itField->first)(idom, it->first).getFunctionSpace()->getNumberOfKeysByEdge()) {
std::vector< int > elementTypes, orientations;
std::vector< std::size_t > edgeNodes, localEdgeTags;
gmsh::model::mesh::getElementTypes(elementTypes, itEntity->first, itEntity->second);
for(auto iElementType = 0ULL; iElementType < elementTypes.size(); ++iElementType) {
gmsh::model::mesh::getElementEdgeNodes(elementTypes[iElementType], edgeNodes, itEntity->second, true);
gmsh::model::mesh::getEdges(edgeNodes, localEdgeTags, orientations);
edgeTags[currentSubdomain][currentInterface][currentField].reserve(edgeTags[currentSubdomain][currentInterface][currentField].size() + localEdgeTags.size());
for(auto i = 0ULL; i < localEdgeTags.size(); ++i) {
edgeTags[currentSubdomain][currentInterface][currentField].push_back(localEdgeTags[i]);
}
}
}
// face (triangle) case
if((*itField->first)(idom, it->first).getFunctionSpace()->getNumberOfKeysByTriangularFace()) {
std::vector< int > elementTypes, orientations;
std::vector< std::size_t > faceNodes, localFaceTags;
gmsh::model::mesh::getElementTypes(elementTypes, itEntity->first, itEntity->second);
for(auto iElementType = 0ULL; iElementType < elementTypes.size(); ++iElementType) {
gmsh::model::mesh::getElementFaceNodes(elementTypes[iElementType], 3, faceNodes, itEntity->second, true);
gmsh::model::mesh::getFaces(3, faceNodes, localFaceTags, orientations);
faceTags[currentSubdomain][currentInterface][currentField].reserve(faceTags[currentSubdomain][currentInterface][currentField].size() + localFaceTags.size());
for(auto i = 0ULL; i < localFaceTags.size(); ++i) {
faceTags[currentSubdomain][currentInterface][currentField].push_back(localFaceTags[i]);
}
}
}
// face (quadrangle) case
if((*itField->first)(idom, it->first).getFunctionSpace()->getNumberOfKeysByQuadrangularFace()) {
std::vector< int > elementTypes, orientations;
std::vector< std::size_t > faceNodes, localFaceTags;
gmsh::model::mesh::getElementTypes(elementTypes, itEntity->first, itEntity->second);
for(auto iElementType = 0ULL; iElementType < elementTypes.size(); ++iElementType) {
gmsh::model::mesh::getElementFaceNodes(elementTypes[iElementType], 4, faceNodes, itEntity->second, true);
gmsh::model::mesh::getFaces(4, faceNodes, localFaceTags, orientations);
faceTags[currentSubdomain][currentInterface][currentField].reserve(faceTags[currentSubdomain][currentInterface][currentField].size() + localFaceTags.size());
for(auto i = 0ULL; i < localFaceTags.size(); ++i) {
faceTags[currentSubdomain][currentInterface][currentField].push_back(localFaceTags[i]);
}
}
}
}
edgeTagsNeighbor[currentSubdomain][currentInterface][currentField].resize(edgeTags[currentSubdomain][currentInterface][currentField].size());
faceTagsNeighbor[currentSubdomain][currentInterface][currentField].resize(faceTags[currentSubdomain][currentInterface][currentField].size());
currentField += 1;
}
currentField = 0;
currentInterface += 1;
}
currentInterface = 0;
currentSubdomain += 1;
}
}
}
// Send the EdgeTags vectors
requests.clear();
for(auto recv : {true, false}) {
unsigned int currentSubdomain = 0;
unsigned int currentInterface = 0;
unsigned int currentField = 0;
for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
if(mpi::isItMySubdomain(idom)) {
for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) {
if(mpi::isItMySubdomain(it->first)) continue;
for(auto itField = _interfaceFields.begin(); itField != _interfaceFields.end(); ++itField) {
if(!(*itField->first).isDefined(idom, it->first)) continue;
auto size = edgeTags[currentSubdomain][currentInterface][currentField].size();
if(size) {
const unsigned int tagIJ = (*itField->first)(idom, it->first).tag();
const unsigned int tagJI = (*itField->second)(it->first, idom).tag();
if(recv) {
requests.push_back(new MPI_Request);
gmshfem::msg::debug << "MPI_Irecv edges rank=" << MPI_Rank << " source=" << mpi::rankOfSubdomain(it->first) << " tag=" << tagJI << " size=" << size << gmshfem::msg::endl;
MPI_Irecv(&edgeTagsNeighbor[currentSubdomain][currentInterface][currentField][0], size, MPI_UNSIGNED_LONG_LONG, mpi::rankOfSubdomain(it->first), tagJI, MPI_COMM_WORLD, requests.back());
}
else {
gmshfem::msg::debug << "MPI_Send edges rank=" << MPI_Rank << " dest=" << mpi::rankOfSubdomain(it->first) << " tag=" << tagIJ << " size=" << size << gmshfem::msg::endl;
MPI_Send(&edgeTags[currentSubdomain][currentInterface][currentField][0], size, MPI_UNSIGNED_LONG_LONG, mpi::rankOfSubdomain(it->first), tagIJ, MPI_COMM_WORLD);
}
}
currentField += 1;
}
currentField = 0;
currentInterface += 1;
}
currentInterface = 0;
currentSubdomain += 1;
}
}
}
// wait until we receive the data
gmshfem::msg::debug << "MPI_Wait edges starting on rank=" << MPI_Rank << gmshfem::msg::endl;
for(auto r : requests) {
MPI_Wait(r, MPI_STATUS_IGNORE);
delete r;
}
gmshfem::msg::debug << "MPI_Wait edges done on rank=" << MPI_Rank << gmshfem::msg::endl;
MPI_Barrier(PETSC_COMM_WORLD);
// Send the FaceTags vectors
requests.clear();
for(auto recv : {true, false}) {
unsigned int currentSubdomain = 0;
unsigned int currentInterface = 0;
unsigned int currentField = 0;
for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
if(mpi::isItMySubdomain(idom)) {
for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) {
if(mpi::isItMySubdomain(it->first)) continue;
for(auto itField = _interfaceFields.begin(); itField != _interfaceFields.end(); ++itField) {
if(!(*itField->first).isDefined(idom, it->first)) continue;
auto size = faceTags[currentSubdomain][currentInterface][currentField].size();
if(size) {
const unsigned int tagIJ = (*itField->first)(idom, it->first).tag();
const unsigned int tagJI = (*itField->second)(it->first, idom).tag();
if(recv) {
requests.push_back(new MPI_Request);
gmshfem::msg::debug << "MPI_Irecv faces rank=" << MPI_Rank << " source=" << mpi::rankOfSubdomain(it->first) << " tag=" << tagJI << " size=" << size << gmshfem::msg::endl;
MPI_Irecv(&faceTagsNeighbor[currentSubdomain][currentInterface][currentField][0], size, MPI_UNSIGNED_LONG_LONG, mpi::rankOfSubdomain(it->first), tagJI, MPI_COMM_WORLD, requests.back());
}
else {
gmshfem::msg::debug << "MPI_Send faces rank=" << MPI_Rank << " dest=" << mpi::rankOfSubdomain(it->first) << " tag=" << tagIJ << " size=" << size << gmshfem::msg::endl;
MPI_Send(&faceTags[currentSubdomain][currentInterface][currentField][0], size, MPI_UNSIGNED_LONG_LONG, mpi::rankOfSubdomain(it->first), tagIJ, MPI_COMM_WORLD);
}
}
currentField += 1;
}
currentField = 0;
currentInterface += 1;
}
currentInterface = 0;
currentSubdomain += 1;
}
}
}
// wait until we receive the data
gmshfem::msg::debug << "MPI_Wait faces starting on rank=" << MPI_Rank << gmshfem::msg::endl;
for(auto r : requests) {
MPI_Wait(r, MPI_STATUS_IGNORE);
delete r;
}
gmshfem::msg::debug << "MPI_Wait faces done on rank=" << MPI_Rank << gmshfem::msg::endl;
MPI_Barrier(PETSC_COMM_WORLD);
// Generate the interface mapping
{
unsigned int currentSubdomain = 0;
unsigned int currentInterface = 0;
unsigned int currentField = 0;
for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
if(mpi::isItMySubdomain(idom)) {
for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) {
if(mpi::isItMySubdomain(it->first)) continue;
for(auto itField = _interfaceFields.begin(); itField != _interfaceFields.end(); ++itField) {
if(!(*itField->first).isDefined(idom, it->first)) continue;
const unsigned int tagIJ = (*itField->first)(idom, it->first).tag();
const unsigned int tagJI = (*itField->second)(it->first, idom).tag();
const unsigned long long numberOfEdges = edgeTags[currentSubdomain][currentInterface][currentField].size();
const unsigned long long numberOfFaces = faceTags[currentSubdomain][currentInterface][currentField].size();
// build maps
std::unordered_map< unsigned long long, unsigned long long > edgeToEdge(numberOfEdges);
std::unordered_map< unsigned long long, unsigned long long > faceToFace(numberOfFaces);
for(auto i = 0U; i < numberOfEdges; ++i) {
edgeToEdge.insert(std::make_pair(edgeTags[currentSubdomain][currentInterface][currentField][i], edgeTagsNeighbor[currentSubdomain][currentInterface][currentField][i]));
}
for(auto i = 0U; i < numberOfFaces; ++i) {
faceToFace.insert(std::make_pair(faceTags[currentSubdomain][currentInterface][currentField][i], faceTagsNeighbor[currentSubdomain][currentInterface][currentField][i]));
}
std::vector< gmshfem::dofs::RawDof > myDofs;
(*itField->second)(it->first, idom).getUnknownDofVector(myDofs);
_interfaceMapping.insert(std::make_pair(tagIJ, std::vector< unsigned long long >(myDofs.size())));
const unsigned int numberOfNodeKey = (*itField->first)(idom, it->first).getFunctionSpace()->getNumberOfKeysByNode();
const unsigned int numberOfEdgeKey = (*itField->first)(idom, it->first).getFunctionSpace()->getNumberOfKeysByEdge();
const unsigned int numberOfFaceKey = (*itField->first)(idom, it->first).getFunctionSpace()->getNumberOfKeysByTriangularFace() + (*itField->first)(idom, it->first).getFunctionSpace()->getNumberOfKeysByQuadrangularFace();
#pragma omp parallel for
for(auto i = 0ULL; i < myDofs.size(); ++i) {
unsigned int firstKey = myDofs[i].numType - GMSHFEM_DOF_FIELD_OFFSET * tagJI;
if(numberOfNodeKey == 0) { // 1-Form as the firstKey for edges always starts with value '1', even if there is no node keys.
firstKey--;
}
if(firstKey < numberOfNodeKey) { // if the dof is assocated to a node
// nothing to do (note that we might want to handle meshes with different node tags on both sides)
}
else if(firstKey >= numberOfNodeKey && firstKey < numberOfNodeKey + numberOfEdgeKey) { // if the dof is assocated to a edge
myDofs[i].entity = edgeToEdge[myDofs[i].entity];
}
else if(firstKey >= numberOfNodeKey + numberOfEdgeKey && firstKey < numberOfNodeKey + numberOfEdgeKey + numberOfFaceKey) { // if the dof is assocated to a face
myDofs[i].entity = faceToFace[myDofs[i].entity];
}
else {
gmshfem::msg::error << "An unknown key is found" << gmshfem::msg::endl;
}
auto itFind = std::find_if(dofKeys[currentSubdomain][currentInterface][currentField]->begin(), dofKeys[currentSubdomain][currentInterface][currentField]->end(), [&](const gmshfem::dofs::RawDofKey &other) { return myDofs[i].numType == other.numType && myDofs[i].entity == other.entity; });
if(itFind == dofKeys[currentSubdomain][currentInterface][currentField]->end()) {
gmshfem::msg::error << "Unable to find the corresponding interface" << gmshfem::msg::endl;
}
_interfaceMapping[tagIJ][i] = std::distance(dofKeys[currentSubdomain][currentInterface][currentField]->begin(), itFind);
}
delete dofKeys[currentSubdomain][currentInterface][currentField];
currentField += 1;
}
currentField = 0;
currentInterface += 1;
}
currentInterface = 0;
currentSubdomain += 1;
}
}
}
// Free the MPI_DOF datatype
MPI_Type_free(&MPI_DOF);
}
MPI_Barrier(PETSC_COMM_WORLD);
for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
if(mpi::isItMySubdomain(idom)) {
gmshfem::msg::info << " - Subdomain " << idom << ":" << gmshfem::msg::endl;
gmshfem::msg::info << " * Subdomain size: " << _volume[idom]->getNumberOfUnknownDof() << " dofs" << gmshfem::msg::endl;
unsigned long long totalInterfaceSize = 0;
for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) {
totalInterfaceSize += it->second->getNumberOfUnknownDof();
}
gmshfem::msg::info << " * Interface size: " << totalInterfaceSize << " dofs:" << gmshfem::msg::endl;
for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) {
gmshfem::msg::info << " · Interface 'g(" << idom << ", " << it->first << ")' size: " << it->second->getNumberOfUnknownDof() << " dofs" << gmshfem::msg::endl;
}
}
}
s_printResources("DDM pre-processing before volume assembly");
// Assemble locally on each processor
for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
if(mpi::isItMySubdomain(idom)) {
gmshfem::msg::info << " - Assembling " << _volume[idom]->name() << gmshfem::msg::endl;
_volume[idom]->setAttribute("ddm::physicalCommutator", _physicalCommutator);
_volume[idom]->setAttribute("ddm::artificialCommutator", _artificialCommutator);
gmshfem::common::Options::instance()->verbose = innerVerbosity;
_volume[idom]->assemble();
gmshfem::common::Options::instance()->verbose = outerVerbosity;
if(gmshddm::common::Options::instance()->memory) {
gmshfem::msg::info << " - Subdomain " << idom << ": number of non zeros : " << _volume[idom]->getNumberOfNonZeros() << gmshfem::msg::endl;
}
}
}
MPI_Barrier(PETSC_COMM_WORLD);
s_printResources("DDM pre-processing before volume solve");
// solve locally on each processor
for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
if(mpi::isItMySubdomain(idom)) {
gmshfem::msg::info << " - Solving " << _volume[idom]->name() << gmshfem::msg::endl;
if(gmshddm::common::Options::instance()->memory) {
gmshfem::msg::info << " - Subdomain " << idom << ": estimated factorization memory : " << _volume[idom]->getEstimatedFactorizationMemoryUsage() << gmshfem::msg::endl;
}
gmshfem::common::Options::instance()->verbose = innerVerbosity;
_volume[idom]->solve();
gmshfem::common::Options::instance()->verbose = outerVerbosity;
}
}
MPI_Barrier(PETSC_COMM_WORLD);
s_printResources("DDM pre-processing before surface assembly and solve");
for(auto idom = 0ULL; idom < _surface.size(); ++idom) {
if(mpi::isItMySubdomain(idom)) {
for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) {
gmshfem::msg::info << " - Assembling and solving " << it->second->name() << gmshfem::msg::endl;
gmshfem::common::Options::instance()->verbose = innerVerbosity;
it->second->assemble();
it->second->solve();
gmshfem::common::Options::instance()->verbose = outerVerbosity;
}
}
}
MPI_Barrier(PETSC_COMM_WORLD);
// get the local interface unknowns for each processor
_rhs.clear();
for(auto idom = 0ULL; idom < _surface.size(); ++idom) {
if(mpi::isItMySubdomain(idom)) {
for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) {
for(auto itField = _interfaceFields.begin(); itField != _interfaceFields.end(); ++itField) {
if((*itField->first).isDefined(idom, it->first)) {
gmshfem::algebra::Vector< T_Scalar > g_ij;
(*itField->first)(idom, it->first).getUnknownVector(g_ij, gmshfem::dofs::RawOrder::Hash);
// concatenate the local unknowns g_ij of the current process and set it in the rhs for the iterative solver
_rhs.concatenate(g_ij);
}
}
}
}
}
MPI_Barrier(PETSC_COMM_WORLD);
s_printResources("DDM pre-processing end");
time.tock();
if(MPI_Rank == 0) {
gmshfem::msg::info << "Done pre-processing in " << time << "s" << gmshfem::msg::endl;
}
return time;
}
template< class T_Scalar >
gmshfem::common::Timer Formulation< T_Scalar >::solve(const std::string &solver, const double tolerance, const int iterMax, const bool sameMatrixWithArtificialAndPhysicalSources)
{
// sameMatrixWithArtificialAndPhysicalSources could be automatically set
// to true if bilinear terms of volume formulations do not use the
// physical/artifical commutator and if there are no non-homogeneous
// Dirichlet BCs
const int outerVerbosity = gmshfem::common::Options::instance()->verbose;
const int innerVerbosity = (outerVerbosity != 0 ? outerVerbosity - 1 : 0);
const unsigned int MPI_Rank = mpi::getMPIRank();
MPI_Barrier(PETSC_COMM_WORLD);
gmshfem::common::Timer time;
time.tick();
if(MPI_Rank == 0) {
gmshfem::msg::info << "Solving " << _name << " using " << solver << "..." << gmshfem::msg::endl;
}
s_printResources("DDM iterative solver begin");
// ******************************************************
// Iterative solver
// ******************************************************
// PETSc
Mat A;
Vec X, Y, RHS;
// local size per proc and global size
PetscInt locSize = _rhs.size();
PetscInt globSize;
// Create Petsc distributed vector - let Petsc compute the global vector size
VecCreateMPI(PETSC_COMM_WORLD, locSize, PETSC_DETERMINE, &RHS);
VecCopy(_rhs.getPetsc(), RHS);
VecGetSize(RHS, &globSize);
VecDuplicate(RHS, &X); // current iterate
VecDuplicate(RHS, &Y); // new iterate
MatCreateShell(PETSC_COMM_WORLD, locSize, locSize, globSize, globSize, this, &A);
MatShellSetOperation(A, MATOP_MULT, (void (*)(void)) & MatVectProduct< T_Scalar >);
_physicalCommutator = false;
_artificialCommutator = true;
togglePhysicalAndArtificialSourceTerms();
auto deactivateBilinear = [=](const unsigned long long idom) {
for(auto it = _volume[idom]->begin(); it != _volume[idom]->end(); ++it) {
if((*it)->isBilinear()) {
(*it)->deactivate();
}
}
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)->deactivate();
}
}
}
};
auto activateBilinear = [=](const unsigned long long idom) {
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();
}
}
}
};
// Assemble the volumic systems for each subdomain
for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
if(mpi::isItMySubdomain(idom)) {
_volume[idom]->setAttribute("ddm::physicalCommutator", _physicalCommutator);
_volume[idom]->setAttribute("ddm::artificialCommutator", _artificialCommutator);
gmshfem::common::Options::instance()->verbose = innerVerbosity;
if(sameMatrixWithArtificialAndPhysicalSources) {
deactivateBilinear(idom);
_volume[idom]->setRHSToZero();
_volume[idom]->assemble();
}
else {
_volume[idom]->removeSystem();
_volume[idom]->initSystem();
_volume[idom]->pre();
_volume[idom]->assemble();
deactivateBilinear(idom);
}
gmshfem::common::Options::instance()->verbose = outerVerbosity;
}
}
if(solver == "jacobi") {
_IA = false;
Vec W; // residual
VecDuplicate(RHS, &W);
VecSet(X, 0.);
PetscReal norm = 1;
int i = 0;
do { // jacobi
MatVectProduct< T_Scalar >(A, X, Y);
VecAYPX(Y, 1., RHS);
VecWAXPY(W, -1., X, Y);
VecCopy(Y, X);
VecNorm(W, NORM_2, &norm);
_runIterationOperations(i, norm);
i++;
} while(_relativeResidual.back() > tolerance && i < iterMax);
VecDestroy(&W);
}
else { // gmres, bcgs, bcgsl, ...
_IA = true;
KSP ksp;
PC pc;
KSPCreate(PETSC_COMM_WORLD, &ksp);
KSPSetFromOptions(ksp);
KSPSetOperators(ksp, A, A);
KSPSetType(ksp, solver.c_str());
KSPMonitorSet(ksp, PetscInterface< T_Scalar, PetscScalar, PetscInt >::kspPrint, this, PETSC_NULL);
KSPSetTolerances(ksp, tolerance, PETSC_DEFAULT, PETSC_DEFAULT, iterMax);
KSPGetPC(ksp, &pc);
PCSetType(pc, PCNONE);
if(solver == "gmres")
KSPGMRESSetRestart(ksp, iterMax);
else if(solver == "bcgsl")
KSPBCGSLSetEll(ksp, 6);
KSPSolve(ksp, RHS, Y);
KSPDestroy(&ksp);
}
// ******************************************************
// Final solution
// ******************************************************
_physicalCommutator = true;
_artificialCommutator = true;
togglePhysicalAndArtificialSourceTerms();
for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
if(mpi::isItMySubdomain(idom)) {
_volume[idom]->setAttribute("ddm::physicalCommutator", _physicalCommutator);
_volume[idom]->setAttribute("ddm::artificialCommutator", _artificialCommutator);
gmshfem::common::Options::instance()->verbose = innerVerbosity;
if(sameMatrixWithArtificialAndPhysicalSources) {
_volume[idom]->setRHSToZero();
}
else {
_volume[idom]->removeSystem();
_volume[idom]->initSystem();
_volume[idom]->pre();
}
gmshfem::common::Options::instance()->verbose = outerVerbosity;
}
}
const PetscScalar *petscArray;
PetscInt petscArraySize;
VecGetArrayRead(Y, &petscArray);
VecGetSize(Y, &petscArraySize);
const T_Scalar *array = PetscInterface< T_Scalar, PetscScalar, PetscInt >::arrayInterface(petscArray, petscArraySize);
_fillG(array);
PetscInterface< T_Scalar, PetscScalar, PetscInt >::freeArray(array);
VecRestoreArrayRead(Y, &petscArray);
// SolveVolumePDE:
for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
if(mpi::isItMySubdomain(idom)) {
gmshfem::common::Options::instance()->verbose = innerVerbosity;
if(!sameMatrixWithArtificialAndPhysicalSources) {
activateBilinear(idom);
}
_volume[idom]->assemble();
_volume[idom]->solve(true);
if(sameMatrixWithArtificialAndPhysicalSources) {
activateBilinear(idom);
}
gmshfem::common::Options::instance()->verbose = outerVerbosity;
}
}
MatDestroy(&A);
VecDestroy(&X);
VecDestroy(&Y);
VecDestroy(&RHS);
MPI_Barrier(PETSC_COMM_WORLD);
s_printResources("DDM iterative solver end");
time.tock();
if(MPI_Rank == 0) {
gmshfem::msg::info << "Done solving in " << time << "s" << gmshfem::msg::endl;
}
return time;
}
// TODO
template< class T_Scalar >
gmshfem::algebra::MatrixCCS< T_Scalar > Formulation< T_Scalar >::computeMatrix()
{
const int outerVerbosity = gmshfem::common::Options::instance()->verbose;
const int innerVerbosity = (outerVerbosity != 0 ? outerVerbosity - 1 : 0);
if(mpi::getMPISize() != 1) {
throw gmshfem::common::Exception("Formulation< T_Scalar >::computeMatrix() does not work on multiple processes");
}
Mat A;
Vec X, Y;
VecDuplicate(_rhs.getPetsc(), &X);
VecDuplicate(_rhs.getPetsc(), &Y);
MatCreateShell(PETSC_COMM_SELF, _rhs.size(), _rhs.size(), _rhs.size(), _rhs.size(), this, &A);
MatShellSetOperation(A, MATOP_MULT, (void (*)(void)) & MatVectProduct< T_Scalar >);
for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
_volume[idom]->removeSystem();
_volume[idom]->initSystem();
}
_physicalCommutator = false;
_artificialCommutator = true;
togglePhysicalAndArtificialSourceTerms();
for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
_volume[idom]->setAttribute("ddm::physicalCommutator", _physicalCommutator);
_volume[idom]->setAttribute("ddm::artificialCommutator", _artificialCommutator);
gmshfem::common::Options::instance()->verbose = innerVerbosity;
_volume[idom]->pre();
gmshfem::common::Options::instance()->verbose = outerVerbosity;
}
_IA = true;
gmshfem::algebra::Vector< T_Scalar > vec(_rhs.size());
std::vector< PetscInt > aiVector(_rhs.size());
for(auto i = 0ULL; i < _rhs.size(); ++i) {
aiVector[i] = i;
}
std::vector< PetscScalar > aVector(_rhs.size());
std::vector< unsigned long long > aj(_rhs.size() + 1, 0);
std::vector< unsigned long long > ai;
ai.reserve(_rhs.size());
std::vector< T_Scalar > a;
a.reserve(_rhs.size());
for(auto i = 0ULL; i < _rhs.size(); ++i) {
gmshfem::msg::info << "Column number " << i << " over " << _rhs.size() << "." << gmshfem::msg::endl;
vec[i] = 1.;
vec.removePetsc();
MatVectProduct< T_Scalar >(A, vec.getPetsc(), Y);
VecGetValues(Y, _rhs.size(), &aiVector[0], &aVector[0]);
aj[i + 1] = aj[i];
for(auto j = 0ULL; j < aVector.size(); ++j) {
if(aVector[j] != 0) {
aj[i + 1]++;
ai.push_back(j);
a.push_back(T_Scalar(aVector[j]));
}
}
ai.reserve(aj.size() + _rhs.size());
a.reserve(a.size() + _rhs.size());
vec[i] = 0.;
}
gmshfem::algebra::MatrixCCS< T_Scalar > mat(std::move(ai), std::move(aj), std::move(a));
return mat;
}
template< class T_Scalar >
const gmshfem::algebra::Vector< T_Scalar > &Formulation< T_Scalar >::getRHS() const
{
return _rhs;
}
template< class T_Scalar >
std::vector< gmshfem::scalar::Precision< T_Scalar > > Formulation< T_Scalar >::absoluteResidual() const
{
return _absoluteResidual;
}
template< class T_Scalar >
std::vector< gmshfem::scalar::Precision< T_Scalar > > Formulation< T_Scalar >::relativeResidual() const
{
return _relativeResidual;
}
template< class T_Scalar >
unsigned int Formulation< T_Scalar >::numberOfIterations() const
{
return _numberOfIterations;
}
template< class T_Scalar >
int MatVectProduct(Mat A, Vec X, Vec Y)
{
PetscFunctionBegin;
const int outerVerbosity = gmshfem::common::Options::instance()->verbose;
const int innerVerbosity = (outerVerbosity != 0 ? outerVerbosity - 1 : 0);
Formulation< T_Scalar > *formulation;
MatShellGetContext(A, (void **)&formulation);
const PetscScalar *petscArray;
PetscInt petscArraySize;
VecGetArrayRead(X, &petscArray);
VecGetSize(X, &petscArraySize);
const T_Scalar *array = Formulation< T_Scalar >::template PetscInterface< T_Scalar, PetscScalar, PetscInt >::arrayInterface(petscArray, petscArraySize);
formulation->_fillG(array);
Formulation< T_Scalar >::template PetscInterface< T_Scalar, PetscScalar, PetscInt >::freeArray(array);
VecRestoreArrayRead(X, &petscArray);
// SolveVolumePDE
for(auto idom = 0ULL; idom < formulation->_volume.size(); ++idom) {
if(mpi::isItMySubdomain(idom)) {
gmshfem::common::Options::instance()->verbose = innerVerbosity;
formulation->_volume[idom]->setRHSToZero();
formulation->_volume[idom]->assemble();
gmshfem::common::Options::instance()->verbose = outerVerbosity;
}
}
for(auto idom = 0ULL; idom < formulation->_volume.size(); ++idom) {
if(mpi::isItMySubdomain(idom)) {
gmshfem::common::Options::instance()->verbose = innerVerbosity;
formulation->_volume[idom]->solve(true);
gmshfem::common::Options::instance()->verbose = outerVerbosity;
}
}
// SolveSurfacePDE
for(auto idom = 0ULL; idom < formulation->_surface.size(); ++idom) {
if(mpi::isItMySubdomain(idom)) {
for(auto it = formulation->_surface[idom].begin(); it != formulation->_surface[idom].end(); ++it) {
gmshfem::common::Options::instance()->verbose = innerVerbosity;
it->second->setRHSToZero();
it->second->assemble();
gmshfem::common::Options::instance()->verbose = outerVerbosity;
}
}
}
for(auto idom = 0ULL; idom < formulation->_surface.size(); ++idom) {
if(mpi::isItMySubdomain(idom)) {
for(auto it = formulation->_surface[idom].begin(); it != formulation->_surface[idom].end(); ++it) {
gmshfem::common::Options::instance()->verbose = innerVerbosity;
it->second->solve(true);
gmshfem::common::Options::instance()->verbose = outerVerbosity;
}
}
}
MPI_Barrier(PETSC_COMM_WORLD);
gmshfem::algebra::Vector< T_Scalar > g_local;
for(auto idom = 0ULL; idom < formulation->_surface.size(); ++idom) {
if(mpi::isItMySubdomain(idom)) {
for(auto it = formulation->_surface[idom].begin(); it != formulation->_surface[idom].end(); ++it) {
for(auto itField = formulation->_interfaceFields.begin(); itField != formulation->_interfaceFields.end(); ++itField) {
if((*itField->first).isDefined(idom, it->first)) {
gmshfem::algebra::Vector< T_Scalar > g_ij;
// get the solved g_ij
(*itField->first)(idom, it->first).getUnknownVector(g_ij, gmshfem::dofs::RawOrder::Hash);
// concatenate to get one vector per process
g_local.concatenate(g_ij);
}
}
}
}
}
VecCopy(g_local.getPetsc(), Y);
if(formulation->_IA) {
VecAYPX(Y, -1., X);
}
PetscFunctionReturn(0);
}
template class Formulation< std::complex< double > >;
template class Formulation< std::complex< float > >;
} // namespace problem
} // namespace gmshddm