Select Git revision
GModelIO_Mesh.cpp
Forked from
gmsh / gmsh
Source project has a limited visibility.
Formulation.cpp 73.28 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 "Solvers.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>
#include <regex>
#ifdef HAVE_MPI
#include <mpi.h>
#endif
#ifdef HAVE_PETSC
#include <petsc.h>
#endif
namespace gmshddm
{
namespace problem
{
using gmshddm::common::s_printResources;
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(), _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(), _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 >
typename Formulation< T_Scalar >::sourceStatus Formulation< T_Scalar >::indexStatusInShots(unsigned id) const {
for (unsigned i = 0; i < _activeShots.size(); ++i) {
if (_numberedPhysicalSourceTerms[i].count(id) > 0) {
if (_activeShots[i]) {
return sourceStatus::ACTIVE;
} else {
return sourceStatus::INACTIVE;
}
}
}
return sourceStatus::NOT_IN_NUMBERED_TERMS;
}
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;
std::vector< MPI_Status > status;
// 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) {
gmshfem::msg::debug << "MPI_Irecv g rank=" << MPI_Rank << " source=" << mpi::rankOfSubdomain(it->first) << " tag=" << tagJI << " size=" << size << gmshfem::msg::endl;
MPI_Request req;
MPI_Irecv(&((*g_ji[currentSubdomain][currentInterface][currentField])[0]), size, MPI_C_DOUBLE_COMPLEX, mpi::rankOfSubdomain(it->first), tagJI, MPI_COMM_WORLD, &req);
requests.push_back(req);
}
}
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;
status.resize(requests.size());
MPI_Waitall(requests.size(), &requests[0], &status[0]);
for(auto s : status)
if(s.MPI_ERROR != MPI_SUCCESS)
gmshfem::msg::warning << "MPI_WaitAll g status " << s.MPI_ERROR << gmshfem::msg::endl;
gmshfem::msg::debug << "MPI_Wait g done on rank=" << MPI_Rank << gmshfem::msg::endl;
MPI_Barrier(MPI_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;
}
}
}
MPI_Barrier(MPI_COMM_WORLD);
}
}
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 >
void Formulation< T_Scalar >::numberedPhysicalSourceTerm(const unsigned int id, const unsigned int termTag)
{
_numberedPhysicalSourceTerms.at(id).insert(termTag);
}
template< class T_Scalar >
void Formulation< T_Scalar >::compareRHSexport()
{
Vec RHS;
Mat BlockRHS;
_buildPetscRHS(&RHS);
_buildPetscAllRHS(&BlockRHS);
{
PetscViewer viewer;
PetscViewerASCIIOpen(PETSC_COMM_WORLD, "output_vec.txt", &viewer);
PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_DENSE);
VecView(RHS, viewer);
PetscViewerDestroy(&viewer);
}
{
PetscViewer viewer;
PetscViewerASCIIOpen(PETSC_COMM_WORLD, "output_block.txt", &viewer);
PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_DENSE);
MatView(BlockRHS, viewer);
PetscViewerDestroy(&viewer);
}
}
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) {
unsigned termTag = (*it)->tag();
auto termStatus = indexStatusInShots(termTag);
if(_physicalSourceTerms.count(termTag) > 0 || termStatus != sourceStatus::NOT_IN_NUMBERED_TERMS) {
if(_physicalCommutator && termStatus != sourceStatus::INACTIVE) {
(*it)->activate();
}
else {
// Disable if all physical terms are disabled OR if it's a disabled shot
(*it)->deactivate();
}
}
else if(_artificialSourceTerms.count(termTag)) {
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< typename T_Scalar >
void Formulation< T_Scalar >::setShotNumber(unsigned N)
{
_activeShots.resize(N, false);
_numberedPhysicalSourceTerms.resize(N);
_multi_rhs.resize(N);
}
template< typename T_Scalar >
void Formulation< T_Scalar >::enableShot(unsigned i)
{
_activeShots.at(i) = true;
}
template< typename T_Scalar >
void Formulation< T_Scalar >::disableShot(unsigned i)
{
_activeShots.at(i) = false;
}
template< typename T_Scalar >
void Formulation< T_Scalar >::toggleShot(unsigned i)
{
_activeShots.at(i) = !_activeShots.at(i);
}
template< typename T_Scalar >
bool Formulation< T_Scalar >::isShotEnabled(unsigned i) const
{
return _activeShots.at(i);
}
template< class T_Scalar >
void Formulation< T_Scalar >::disableAllShots()
{
for (unsigned i = 0; i < _activeShots.size(); ++i)
_activeShots.at(i) = false;
}
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 >
void Formulation< T_Scalar >::addInterfaceField(field::InterfaceFieldInterface< T_Scalar > &interfaceField, field::InterfaceFieldInterface< T_Scalar > &interfaceFieldMapping)
{
_interfaceFields.insert(std::make_pair(&interfaceField, &interfaceFieldMapping));
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(bool mustAssemble)
{
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(MPI_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) {
if(it->second == nullptr) {
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("Unable to match the interface field '" + it->first->name() + "' with another interface field; give it manally");
}
}
}
}
// ******************************************************
// 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
_preProVolume();
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< MPI_Status > status;
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()) {
gmshfem::msg::debug << "MPI_Irecv dofs rank=" << MPI_Rank << " source=" << mpi::rankOfSubdomain(it->first) << " tag=" << tagJI << " size=" << dofs.size() << gmshfem::msg::endl;
MPI_Request req;
MPI_Irecv(&(*dofKeys[currentSubdomain][currentInterface][currentField])[0], dofs.size(), MPI_DOF, mpi::rankOfSubdomain(it->first), tagJI, MPI_COMM_WORLD, &req);
requests.push_back(req);
}
}
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;
status.resize(requests.size());
MPI_Waitall(requests.size(), &requests[0], &status[0]);
for(auto s : status)
if(s.MPI_ERROR != MPI_SUCCESS)
gmshfem::msg::warning << "MPI_WaitAll dofs status " << s.MPI_ERROR << gmshfem::msg::endl;
requests.clear();
status.clear();
gmshfem::msg::debug << "MPI_Wait dofs done on rank=" << MPI_Rank << gmshfem::msg::endl;
MPI_Barrier(MPI_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
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) {
gmshfem::msg::debug << "MPI_Irecv edges rank=" << MPI_Rank << " source=" << mpi::rankOfSubdomain(it->first) << " tag=" << tagJI << " size=" << size << gmshfem::msg::endl;
MPI_Request req;
MPI_Irecv(&edgeTagsNeighbor[currentSubdomain][currentInterface][currentField][0], size, MPI_UNSIGNED_LONG_LONG, mpi::rankOfSubdomain(it->first), tagJI, MPI_COMM_WORLD, &req);
requests.push_back(req);
}
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;
status.resize(requests.size());
MPI_Waitall(requests.size(), &requests[0], &status[0]);
for(auto s : status)
if(s.MPI_ERROR != MPI_SUCCESS)
gmshfem::msg::warning << "MPI_WaitAll edges status " << s.MPI_ERROR << gmshfem::msg::endl;
requests.clear();
status.clear();
gmshfem::msg::debug << "MPI_Wait edges done on rank=" << MPI_Rank << gmshfem::msg::endl;
MPI_Barrier(MPI_COMM_WORLD);
// Send the FaceTags vectors
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) {
gmshfem::msg::debug << "MPI_Irecv faces rank=" << MPI_Rank << " source=" << mpi::rankOfSubdomain(it->first) << " tag=" << tagJI << " size=" << size << gmshfem::msg::endl;
MPI_Request req;
MPI_Irecv(&faceTagsNeighbor[currentSubdomain][currentInterface][currentField][0], size, MPI_UNSIGNED_LONG_LONG, mpi::rankOfSubdomain(it->first), tagJI, MPI_COMM_WORLD, &req);
requests.push_back(req);
}
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;
status.resize(requests.size());
MPI_Waitall(requests.size(), &requests[0], &status[0]);
for(auto s : status)
if(s.MPI_ERROR != MPI_SUCCESS)
gmshfem::msg::warning << "MPI_WaitAll faces status " << s.MPI_ERROR << gmshfem::msg::endl;
requests.clear();
status.clear();
gmshfem::msg::debug << "MPI_Wait faces done on rank=" << MPI_Rank << gmshfem::msg::endl;
MPI_Barrier(MPI_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((*itField->first)(idom, it->first).multiplicity() > 1) { // To deal with compound fields
unsigned int multiplicity = 1;
while(firstKey >= GMSHFEM_DOF_FIELD_OFFSET || multiplicity < (*itField->first)(idom, it->first).multiplicity()) {
firstKey -= GMSHFEM_DOF_FIELD_OFFSET;
++multiplicity;
}
}
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 (firstKey = " << firstKey << ")" << 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(MPI_COMM_WORLD);
_infoNumberOfDofs();
if(mustAssemble) {
assemble();
time.tock();
if(MPI_Rank == 0) {
gmshfem::msg::info << "Done pre-processing and assembling in " << time << "s" << gmshfem::msg::endl;
}
return time;
}
else {
time.tock();
if(MPI_Rank == 0) {
gmshfem::msg::info << "Done pre-processing (no assembly) in " << time << "s" << gmshfem::msg::endl;
}
return time;
}
}
template< class T_Scalar >
gmshfem::common::Timer Formulation< T_Scalar >::assemble(int storeID)
{
gmshfem::common::Timer time;
time.tick();
const unsigned int MPI_Rank = mpi::getMPIRank();
_physicalCommutator = true;
_artificialCommutator = false;
togglePhysicalAndArtificialSourceTerms();
// Clear volume
for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
if(mpi::isItMySubdomain(idom)) {
// Clear formulation
_volume[idom]->setSystemToZero();
}
}
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 << " - Clearing " << it->second->name() << gmshfem::msg::endl;
it->second->setSystemToZero();
}
}
}
// Assemble locally on each processor
s_printResources("DDM pre-processing before volume assembly");
_assembleAllVolume();
MPI_Barrier(MPI_COMM_WORLD);
s_printResources("DDM pre-processing before volume solve");
// solve locally on each processor
_solveAllVolume();
MPI_Barrier(MPI_COMM_WORLD);
s_printResources("DDM pre-processing before surface assembly and solve");
_assembleAndSolveSurface();
MPI_Barrier(MPI_COMM_WORLD);
// get the local interface unknowns for each processor
_extractRHS();
if (storeID > -1) {
if (storeID >= _activeShots.size()) {
throw gmshfem::common::Exception("Invalid storeID argument to gmshddm::problem::Formulation::assemble()");
}
_multi_rhs.at(storeID) = _rhs; // Do a copy. Make moving an option ?
}
MPI_Barrier(MPI_COMM_WORLD);
s_printResources("DDM pre-processing end");
time.tock();
if(MPI_Rank == 0) {
gmshfem::msg::info << "Done assembling 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, const bool skipFinalSolutionComputation)
{
std::unique_ptr < AbstractIterativeSolver< T_Scalar >> solverObject;
// For backwards compatibility.
if(solver == "jacobi") {
solverObject = std::make_unique< JacobiIterativeSolver< T_Scalar > >();
}
else {
solverObject = std::make_unique< KSPIterativeSolver< T_Scalar > >(solver);
}
return solve(*solverObject, tolerance, iterMax, sameMatrixWithArtificialAndPhysicalSources, skipFinalSolutionComputation);
}
template< class T_Scalar >
gmshfem::common::Timer Formulation< T_Scalar >::solve(AbstractIterativeSolver<T_Scalar>& solver, const double tolerance, const int iterMax, const bool sameMatrixWithArtificialAndPhysicalSources, const bool skipFinalSolutionComputation)
{
// 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(MPI_COMM_WORLD);
gmshfem::common::Timer time;
time.tick();
if(MPI_Rank == 0) {
gmshfem::msg::info << "Solving " << _name << " using " << solver.solverName() << "..." << 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
_buildPetscRHS(&RHS);
VecGetSize(RHS, &globSize);
VecDuplicate(RHS, &X); // current iterate
VecDuplicate(RHS, &Y); // new iterate
MatCreateShell(MPI_COMM_WORLD, locSize, locSize, globSize, globSize, this, &A);
switch(solver.productType()) {
case ProductType::A:
MatShellSetOperation(A, MATOP_MULT, (void (*)(void)) & MatVectProductA< T_Scalar >);
break;
case ProductType::I_A:
MatShellSetOperation(A, MATOP_MULT, (void (*)(void)) & MatVectProductI_A< T_Scalar >);
break;
default:
throw gmshfem::common::Exception("Unknown product type.");
}
_physicalCommutator = false;
_artificialCommutator = true;
togglePhysicalAndArtificialSourceTerms();
// 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;
}
}
// Call the AbstractIterativeSolver object
solver.solve(A, X, Y, RHS, tolerance, iterMax);
_numberOfIterations = solver.absoluteResidual().size();
// ******************************************************
// Final solution
// ******************************************************
if(skipFinalSolutionComputation) {
time.tock();
return time;
}
_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;
}
}
_fromPetscToInterfaceFields(Y);
// 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(MPI_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;
}
template< class T_Scalar >
gmshfem::algebra::MatrixCCS< T_Scalar > Formulation< T_Scalar >::computeMatrix()
{
gmshfem::common::Timer time;
time.tick();
const int outerVerbosity = gmshfem::common::Options::instance()->verbose;
const int innerVerbosity = (outerVerbosity != 0 ? outerVerbosity - 1 : 0);
Mat A;
Vec X, Y, Y_seq;
VecScatter ctx;
// Create a distributed vector that acts as a RHS
PetscInt locSize = _rhs.size();
VecCreateMPI(MPI_COMM_WORLD, locSize, PETSC_DETERMINE, &Y);
VecDuplicate(Y, &X);
// Create a sequential full vector and a scatter object to synchronize the global vector of outputs,
// Needed to build the matrix
PetscInt globSize;
VecGetSize(Y, &globSize);
VecCreateSeq(MPI_COMM_SELF, globSize, &Y_seq);
VecScatterCreateToAll(Y, &ctx, &Y_seq);
if (mpi::getMPIRank() == 0)
gmshfem::msg::info << "Local size:" << locSize << " and global size:" << globSize << "." << gmshfem::msg::endl;
// Setup the A matrix
MatCreateShell(MPI_COMM_WORLD, locSize, locSize, globSize, globSize, this, &A);
MatShellSetOperation(A, MATOP_MULT, (void (*)(void)) & MatVectProductI_A< T_Scalar >);
_physicalCommutator = false;
_artificialCommutator = true;
togglePhysicalAndArtificialSourceTerms();
// Compute sub matrices then disable bilinear terms (to assemble only RHS)
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]->removeSystem();
_volume[idom]->initSystem();
_volume[idom]->pre();
_volume[idom]->assemble();
_deactivateBilinear(idom);
gmshfem::common::Options::instance()->verbose = outerVerbosity;
}
gmshfem::algebra::Vector< T_Scalar > vec(globSize);
std::vector< PetscInt > aiVector(globSize);
for(auto i = 0ULL; i < globSize; ++i) {
aiVector[i] = i;
}
std::vector< PetscScalar > aVector(globSize);
std::vector< unsigned long long > aj(globSize + 1, 0);
std::vector< unsigned long long > ai;
ai.reserve(globSize);
std::vector< T_Scalar > a;
a.reserve(globSize);
std::vector<PetscScalar> zeros(globSize, 0);
// Init RHS vector
VecSet(X, 0.0);
VecSetValue(X, 0, 1.0, INSERT_VALUES);
VecAssemblyBegin(X);
VecAssemblyEnd(X);
for(auto i = 0ULL; i < globSize; ++i) {
VecSet(X, 0.0);
VecSetValue(X, i, 1.0, INSERT_VALUES);
VecAssemblyBegin(X);
VecAssemblyEnd(X);
if(mpi::getMPIRank() == 0) {
gmshfem::msg::info << "Column number " << i << " over " << globSize << "." << gmshfem::msg::endl;
}
MatMult(A, X, Y);
// Scatter to rank 0
VecScatterBegin(ctx, Y, Y_seq, INSERT_VALUES, SCATTER_FORWARD);
VecScatterEnd(ctx, Y, Y_seq, INSERT_VALUES, SCATTER_FORWARD);
VecGetValues(Y_seq, globSize, &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() + globSize);
a.reserve(a.size() + globSize);
// Update RHS
VecSetValue(X, i, 0.0, INSERT_VALUES);
if(i + 1 < globSize)
VecSetValue(X, i + 1, 1.0, INSERT_VALUES);
VecAssemblyBegin(X);
VecAssemblyEnd(X);
}
for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
_activateBilinear(idom);
}
VecScatterDestroy(&ctx);
gmshfem::algebra::MatrixCCS< T_Scalar > mat(std::move(ai), std::move(aj), std::move(a));
time.tock();
if(mpi::getMPIRank() == 0) {
gmshfem::msg::info << "Computed the matrix in " << time.time() << "s." << gmshfem::msg::endl;
}
return mat;
}
template< class T_Scalar >
gmshfem::algebra::MatrixCCS< T_Scalar > Formulation< T_Scalar >::computeMatrixBlock()
{
gmshfem::common::Timer time;
time.tick();
const int outerVerbosity = gmshfem::common::Options::instance()->verbose;
const int innerVerbosity = (outerVerbosity != 0 ? outerVerbosity - 1 : 0);
Mat A;
Vec Y, Y_seq;
VecScatter ctx;
// Create a distributed vector that acts as a RHS
PetscInt locSize = _rhs.size();
VecCreateMPI(MPI_COMM_WORLD, locSize, PETSC_DETERMINE, &Y);
VecSetUp(Y);
// Create a sequential full vector and a scatter object to synchronize the global vector of outputs,
// Needed to build the matrix
PetscInt globSize;
VecGetSize(Y, &globSize);
// Create the identity matrix
Mat ID;
{
MatCreate(PETSC_COMM_WORLD, &ID);
MatSetSizes(ID, locSize, PETSC_DECIDE, globSize, globSize);
MatSetType(ID, MATMPIDENSE);
MatSetUp(ID);
MatZeroEntries(ID);
PetscInt start, end;
MatGetOwnershipRange(ID, &start, &end);
for(PetscInt i = start; i < end; ++i) {
PetscScalar val = 1.0;
MatSetValues(ID, 1, &i, 1, &i, &val, INSERT_VALUES);
}
MatAssemblyBegin(ID, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(ID, MAT_FINAL_ASSEMBLY);
}
VecCreateSeq(MPI_COMM_SELF, globSize, &Y_seq);
VecScatterCreateToAll(Y, &ctx, &Y_seq);
if (mpi::getMPIRank() == 0)
gmshfem::msg::info << "Local size:" << locSize << " and global size:" << globSize << "." << gmshfem::msg::endl;
// Setup the A matrix
MatCreateShell(MPI_COMM_WORLD, locSize, locSize, globSize, globSize, this, &A);
MatShellSetOperation(A, MATOP_MULT, (void (*)(void)) & MatVectProductI_A< T_Scalar >);
_physicalCommutator = false;
_artificialCommutator = true;
togglePhysicalAndArtificialSourceTerms();
// Compute sub matrices then disable bilinear terms (to assemble only RHS)
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]->removeSystem();
_volume[idom]->initSystem();
_volume[idom]->pre();
_volume[idom]->assemble();
_deactivateBilinear(idom);
gmshfem::common::Options::instance()->verbose = outerVerbosity;
}
Mat Yblock;
MatDuplicate(ID, MAT_DO_NOT_COPY_VALUES, &Yblock);
MatMatProductI_A< T_Scalar >(A, ID, Yblock);
std::vector< PetscInt > aiVector(globSize);
std::vector< PetscScalar > aVector(globSize);
std::vector< unsigned long long > aj(globSize + 1, 0);
std::vector< unsigned long long > ai;
ai.reserve(globSize);
std::vector< T_Scalar > a;
a.reserve(globSize);
for(auto i = 0ULL; i < globSize; ++i) {
aiVector[i] = i;
}
for(auto i = 0ULL; i < globSize; ++i) {
if(mpi::getMPIRank() == 0) {
gmshfem::msg::info << "Column number " << i << " over " << globSize << "." << gmshfem::msg::endl;
}
// Scatter to rank 0
MatGetColumnVector(Yblock, Y, i);
VecScatterBegin(ctx, Y, Y_seq, INSERT_VALUES, SCATTER_FORWARD);
VecScatterEnd(ctx, Y, Y_seq, INSERT_VALUES, SCATTER_FORWARD);
VecGetValues(Y_seq, globSize, &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() + globSize);
a.reserve(a.size() + globSize);
}
for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
_activateBilinear(idom);
}
VecScatterDestroy(&ctx);
gmshfem::algebra::MatrixCCS< T_Scalar > mat(std::move(ai), std::move(aj), std::move(a));
time.tock();
if(mpi::getMPIRank() == 0) {
gmshfem::msg::info << "Computed the matrix in " << time.time() << "s." << gmshfem::msg::endl;
}
return mat;
}
template< class T_Scalar >
void Formulation< T_Scalar >::setSolutionIntoInterfaceFields(const std::vector< T_Scalar > &solution)
{
_fillG(solution.data());
}
template< class T_Scalar >
const gmshfem::algebra::Vector< T_Scalar > &Formulation< T_Scalar >::getRHS() const
{
return _rhs;
}
template< class T_Scalar >
unsigned int Formulation< T_Scalar >::numberOfIterations() const
{
return _numberOfIterations;
}
template< class T_Scalar >
int MatVectProductA(Mat A, Vec X, Vec Y)
{
return MatVectProductImpl<T_Scalar>(A, X, Y, false);
}
template< class T_Scalar >
int MatVectProductI_A(Mat A, Vec X, Vec Y)
{
return MatVectProductImpl<T_Scalar>(A, X, Y, true);
}
template< class T_Scalar >
int MatMatProductI_A(Mat A, Mat X, Mat Y)
{
return MatMatProductImpl<T_Scalar>(A, X, Y, true);
}
template< class T_Scalar >
int MatMatProductA(Mat A, Mat X, Mat Y)
{
return MatMatProductImpl<T_Scalar>(A, X, Y, false);
}
template< class T_Scalar >
void Formulation<T_Scalar>::_fromPetscToInterfaceFields(Vec G)
{
const PetscScalar *petscArray;
PetscInt petscArraySize;
VecGetArrayRead(G, &petscArray);
VecGetLocalSize(G, &petscArraySize);
const T_Scalar *array = PetscInterface< T_Scalar, PetscScalar, PetscInt >::arrayInterface(petscArray, petscArraySize);
_fillG(array);
PetscInterface< T_Scalar, PetscScalar, PetscInt >::freeArray(array);
VecRestoreArrayRead(G, &petscArray);
}
template< class T_Scalar >
void Formulation< T_Scalar >::_preProVolume()
{
const int outerVerbosity = gmshfem::common::Options::instance()->verbose;
const int innerVerbosity = (outerVerbosity != 0 ? outerVerbosity - 1 : 0);
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;
}
}
}
template< class T_Scalar >
void Formulation< T_Scalar >::_infoNumberOfDofs() const
{
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;
}
}
}
}
template< class T_Scalar >
void Formulation< T_Scalar >::_assembleAllVolume()
{
const int outerVerbosity = gmshfem::common::Options::instance()->verbose;
const int innerVerbosity = (outerVerbosity != 0 ? outerVerbosity - 1 : 0);
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;
}
}
}
}
template< class T_Scalar >
void Formulation< T_Scalar >::_solveAllVolume()
{
const int outerVerbosity = gmshfem::common::Options::instance()->verbose;
const int innerVerbosity = (outerVerbosity != 0 ? outerVerbosity - 1 : 0);
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;
}
}
}
template< class T_Scalar >
void Formulation< T_Scalar >::_assembleAndSolveSurface()
{
const int outerVerbosity = gmshfem::common::Options::instance()->verbose;
const int innerVerbosity = (outerVerbosity != 0 ? outerVerbosity - 1 : 0);
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;
}
}
}
}
template< class T_Scalar >
void Formulation< T_Scalar >::_extractRHS()
{
_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);
}
}
}
}
}
}
template< class T_Scalar >
void Formulation< T_Scalar >::_activateBilinear(const unsigned long long idom)
{
for(auto formulation : *_volume[idom]) {
if(formulation->isBilinear()) {
formulation->activate();
}
}
for(auto formulation : _surface[idom]) {
for(auto term : *formulation.second) {
if(term->isBilinear()) {
term->activate();
}
}
}
}
template< class T_Scalar >
void Formulation< T_Scalar >::_deactivateBilinear(const unsigned long long idom)
{
for(auto formulation : *_volume[idom]) {
if(formulation->isBilinear()) {
formulation->deactivate();
}
}
for(auto formulation : _surface[idom]) {
for(auto term : *formulation.second) {
if(term->isBilinear()) {
term->deactivate();
}
}
}
}
template< class T_Scalar >
void Formulation< T_Scalar >::_buildPetscRHS(Vec *RHS)
{
auto locSize = _rhs.size();
VecCreateMPI(MPI_COMM_WORLD, locSize, PETSC_DETERMINE, RHS);
VecCopy(_rhs.getPetsc(), *RHS);
}
template< class T_Scalar >
void Formulation< T_Scalar >::_buildPetscAllRHS(Mat *RHS)
{
gmshfem::common::Timer time;
time.tick();
if(mpi::getMPIRank() == 0) {
gmshfem::msg::info << "Starting to build a block of RHS." << gmshfem::msg::endl;
}
auto locSize = _rhs.size();
auto nRHS = _activeShots.size();
MatCreateDense(MPI_COMM_WORLD, locSize, PETSC_DETERMINE, PETSC_DETERMINE, nRHS, NULL, RHS);
MatSetUp(*RHS);
std::vector< PetscInt > rowIndices(locSize);
PetscInt startRow, endRow;
MatGetOwnershipRange(*RHS, &startRow, &endRow);
for(PetscInt i = 0; i < locSize; i++) {
rowIndices[i] = startRow + i;
}
// For each RHS, copy the vector into the column of the block
for(PetscInt j = 0; j < nRHS; j++) {
Vec bj = _multi_rhs[j].getPetsc();
PetscScalar *colValues;
VecGetArray(bj, &colValues);
MatSetValues(*RHS, locSize, rowIndices.data(), 1, &j, colValues, INSERT_VALUES);
VecRestoreArray(bj, &colValues);
}
MatAssemblyBegin(*RHS, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(*RHS, MAT_FINAL_ASSEMBLY);
time.tock();
if(mpi::getMPIRank() == 0) {
gmshfem::msg::info << "Done assembling block of RHS in " << time << "s" << gmshfem::msg::endl;
}
}
template< class T_Scalar >
int MatVectProductImpl(Mat A, Vec X, Vec Y, bool IA)
{
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);
formulation->_fromPetscToInterfaceFields(X);
// 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(MPI_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(IA) {
VecAYPX(Y, -1., X);
}
PetscFunctionReturn(0);
}
template< class T_Scalar >
int MatMatProductImpl(Mat A, Mat X, Mat Y, bool IA)
{
// EXTREMLY UNOPTIMIZED
PetscFunctionBegin;
PetscInt n, p, n_local, p_local;
MatGetSize(X, &n, &p);
MatGetLocalSize(X, &n_local, &p_local);
Vec Xcol, Ycol;
VecCreateMPI(MPI_COMM_WORLD, n_local, n, &Xcol);
VecDuplicate(Xcol, &Ycol);
for (auto j = 0; j < p; ++j) {
MatGetColumnVector(X, Xcol, j);
MatVectProductImpl<T_Scalar>(A, Xcol, Ycol, IA);
// Copy to matrix
PetscInt low, high, i;
PetscScalar val;
// Get the ownership range for the current rank
VecGetOwnershipRange(Ycol, &low, &high);
// Loop over local rows and set the value for the specified column
for (i = low; i < high; ++i) {
VecGetValues(Ycol, 1, &i, &val);
MatSetValue(Y, i, j, val, INSERT_VALUES);
}
}
PetscFunctionReturn(0);
}
template class Formulation< std::complex< double > >;
template class Formulation< std::complex< float > >;
} // namespace problem
} // namespace gmshddm