Skip to content

Commit

Permalink
BUG: Fix self gravity (#202)
Browse files Browse the repository at this point in the history 
* fix Domain decomposition+Self gravity
* ensure that tests now validate MPI+Self-gravity
* add P. Segretin's test on dust+SelfGravity
* Update patch version to v2.0.01
  • Loading branch information
@glesur
glesur authored Oct 26, 2023
1 parent a5d51f5 commit 51022cc
Show file tree
Hide file tree
Showing 16 changed files with 325 additions and 149 deletions.
4 changes: 4 additions & 0 deletions .github/workflows/idefix-ci.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -186,6 +186,10 @@ jobs:
run: |
cd $IDEFIX_DIR/test/SelfGravity/UniformCollapse
./testme.py -all $TESTME_OPTIONS
- name: Dusty spherical collapse
run: |
cd $IDEFIX_DIR/test/SelfGravity/DustyCollapse
./testme.py -all $TESTME_OPTIONS
Planet:
needs: [Linter, ShocksHydro, ParabolicHydro, ShocksMHD, ParabolicMHD]
Expand Down
4 changes: 4 additions & 0 deletions CHANGELOG.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,10 @@ All notable changes to this project will be documented in this file.
The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).

## [2.0.1] 2023-10-26
### Changed
- fixed a bug in the self-gravity module introduced in #194 that resulted in instabilities when used in combination with MPI (#202).

## [2.0.0] 2023-10-23
### Changed
- Reorganisation of class instances embedded in datablocks to allow for multi-fluid problems using a generalised template class "Fluid" that replaces "Hydro". Most instances (like data.hydro.Vc) have been replaced by pointers (e.g data.hydro->Vc) (!307).
Expand Down
2 changes: 1 addition & 1 deletion CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@ set (CMAKE_CXX_STANDARD 17)

set(Idefix_VERSION_MAJOR 2)
set(Idefix_VERSION_MINOR 0)
set(Idefix_VERSION_PATCH 00)
set(Idefix_VERSION_PATCH 01)

project (idefix VERSION 1.0.0)
option(Idefix_MHD "enable MHD" OFF)
Expand Down
49 changes: 30 additions & 19 deletions src/gravity/laplacian.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@
Laplacian::Laplacian(DataBlock *datain, std::array<LaplacianBoundaryType,3> leftBound,
std::array<LaplacianBoundaryType,3> rightBound,
bool havePreconditionnerIn) {
idfx::pushRegion("Laplacian::Laplacian");
// Save the parents data objects
this->data = datain;
this->havePreconditioner = havePreconditionnerIn;
Expand All @@ -40,10 +41,8 @@ Laplacian::Laplacian(DataBlock *datain, std::array<LaplacianBoundaryType,3> left

isPeriodic = true;
for(int dir = 0 ; dir < 3 ; dir++) {
if(lbound[dir] != LaplacianBoundaryType::periodic &&
lbound[dir] != LaplacianBoundaryType::internalgrav) isPeriodic = false;
if(rbound[dir] != LaplacianBoundaryType::periodic &&
rbound[dir] != LaplacianBoundaryType::internalgrav) isPeriodic = false;
if(lbound[dir] != LaplacianBoundaryType::periodic) isPeriodic = false;
if(rbound[dir] != LaplacianBoundaryType::periodic) isPeriodic = false;
}

#ifdef WITH_MPI
Expand All @@ -52,8 +51,19 @@ Laplacian::Laplacian(DataBlock *datain, std::array<LaplacianBoundaryType,3> left
int remainDims[3] = {false, true, true};
MPI_SAFE_CALL(MPI_Cart_sub(data->mygrid->CartComm, remainDims, &originComm));
}
#endif

// Update internal boundaries in case of domain decomposition
for(int dir = 0 ; dir < DIMENSIONS ; dir++) {
if(data->mygrid->nproc[dir]>1) {
if(this->data->lbound[dir]==internal ) {
this->lbound[dir] = Laplacian::internalgrav;
}
if(this->data->rbound[dir]==internal) {
this->rbound[dir] = Laplacian::internalgrav;
}
}
}
#endif

#if GEOMETRY == SPHERICAL
if ((this->rbound[JDIR]==axis) || (this->lbound[JDIR]==axis)) {
Expand All @@ -72,20 +82,6 @@ Laplacian::Laplacian(DataBlock *datain, std::array<LaplacianBoundaryType,3> left
}
#endif

// Init MPI stack when needed
#ifdef WITH_MPI
this->arr4D = IdefixArray4D<real> ("WorkingArrayMpi", 1, this->np_tot[KDIR],
this->np_tot[JDIR],
this->np_tot[IDIR]);

int ntarget = 0;
std::vector<int> mapVars;
mapVars.push_back(ntarget);

this->mpi.Init(data->mygrid, mapVars, this->nghost.data(), this->np_int.data());
#endif


if(this->lbound[IDIR] == origin) {
InitInternalGrid();
}
Expand All @@ -96,6 +92,21 @@ Laplacian::Laplacian(DataBlock *datain, std::array<LaplacianBoundaryType,3> left
}
// Initialise the Laplacian coefficients
PreComputeLaplacian();

// Init MPI stack when needed
#ifdef WITH_MPI
this->arr4D = IdefixArray4D<real> ("WorkingArrayMpi", 1, this->np_tot[KDIR],
this->np_tot[JDIR],
this->np_tot[IDIR]);

int ntarget = 0;
std::vector<int> mapVars;
mapVars.push_back(ntarget);

this->mpi.Init(data->mygrid, mapVars, this->nghost.data(), this->np_int.data());
#endif

idfx::popRegion();
}

void Laplacian::InitInternalGrid() {
Expand Down
141 changes: 35 additions & 106 deletions src/gravity/selfGravity.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,8 +5,7 @@
// Licensed under CeCILL 2.1 License, see COPYING for more information
// ***********************************************************************************



#include <memory>
#include <string>
#include <vector>

Expand Down Expand Up @@ -106,20 +105,6 @@ void SelfGravity::Init(Input &input, DataBlock *datain) {
}
}

// Update internal boundaries in case of domain decomposition
#ifdef WITH_MPI
for(int dir = 0 ; dir < DIMENSIONS ; dir++) {
if(data->mygrid->nproc[dir]>1) {
if(this->data->lbound[dir]==internal ) {
this->lbound[dir] = Laplacian::internalgrav;
}
if(this->data->rbound[dir]==internal) {
this->rbound[dir] = Laplacian::internalgrav;
}
}
}
#endif

// Update solver when provided
if(input.CheckEntry("SelfGravity","solver") >= 0) {
std::string strSolver = input.Get<std::string>("SelfGravity","solver",0);
Expand Down Expand Up @@ -162,29 +147,28 @@ void SelfGravity::Init(Input &input, DataBlock *datain) {
}

// Make the Laplacian operator
laplacian = Laplacian(data, lbound, rbound, this->havePreconditioner );
laplacian = std::make_unique<Laplacian>(data, lbound, rbound, this->havePreconditioner );

np_tot = laplacian.np_tot;
np_tot = laplacian->np_tot;

// Instantiate the bicgstab solver
if(solver == BICGSTAB || solver == PBICGSTAB) {
iterativeSolver = new Bicgstab(laplacian, targetError, maxiter,
laplacian.np_tot, laplacian.beg, laplacian.end);
iterativeSolver = new Bicgstab(*laplacian.get(), targetError, maxiter,
laplacian->np_tot, laplacian->beg, laplacian->end);
} else if(solver == CG || solver == PCG) {
iterativeSolver = new Cg(laplacian, targetError, maxiter,
laplacian.np_tot, laplacian.beg, laplacian.end);
iterativeSolver = new Cg(*laplacian.get(), targetError, maxiter,
laplacian->np_tot, laplacian->beg, laplacian->end);
} else if(solver == MINRES || solver == PMINRES) {
iterativeSolver = new Minres(laplacian,
iterativeSolver = new Minres(*laplacian.get(),
targetError, maxiter,
laplacian.np_tot, laplacian.beg, laplacian.end);
laplacian->np_tot, laplacian->beg, laplacian->end);
} else {
real step = laplacian.ComputeCFL();
iterativeSolver = new Jacobi(laplacian, targetError, maxiter, step,
laplacian.np_tot, laplacian.beg, laplacian.end);
real step = laplacian->ComputeCFL();
iterativeSolver = new Jacobi(*laplacian.get(), targetError, maxiter, step,
laplacian->np_tot, laplacian->beg, laplacian->end);
}



// Arrays initialisation
this->density = IdefixArray3D<real> ("Density", this->np_tot[KDIR],
this->np_tot[JDIR],
Expand Down Expand Up @@ -246,7 +230,7 @@ void SelfGravity::ShowConfig() {
}

if(this->lbound[IDIR] == Laplacian::LaplacianBoundaryType::origin) {
idfx::cout << "SelfGravity: using origin boundary with " << laplacian.loffset[IDIR]
idfx::cout << "SelfGravity: using origin boundary with " << laplacian->loffset[IDIR]
<< " additional radial points." << std::endl;
}

Expand All @@ -268,9 +252,9 @@ void SelfGravity::InitSolver() {

// Initialise the density field
// todo: check bounds
int ioffset = laplacian.loffset[IDIR];
int joffset = laplacian.loffset[JDIR];
int koffset = laplacian.loffset[KDIR];
int ioffset = laplacian->loffset[IDIR];
int joffset = laplacian->loffset[JDIR];
int koffset = laplacian->loffset[KDIR];

idefix_for("InitDensity", data->beg[KDIR], data->end[KDIR],
data->beg[JDIR], data->end[JDIR],
Expand Down Expand Up @@ -300,13 +284,13 @@ void SelfGravity::InitSolver() {
// divide density by preconditionner if we're doing the preconditionned version
if(havePreconditioner) {
int ibeg, iend, jbeg, jend, kbeg, kend;
ibeg = laplacian.beg[IDIR];
iend = laplacian.end[IDIR];
jbeg = laplacian.beg[JDIR];
jend = laplacian.end[JDIR];
kbeg = laplacian.beg[KDIR];
kend = laplacian.end[KDIR];
IdefixArray3D<real> P = laplacian.precond;
ibeg = laplacian->beg[IDIR];
iend = laplacian->end[IDIR];
jbeg = laplacian->beg[JDIR];
jend = laplacian->end[JDIR];
kbeg = laplacian->beg[KDIR];
kend = laplacian->end[KDIR];
IdefixArray3D<real> P = laplacian->precond;
idefix_for("Precond density", kbeg, kend, jbeg, jend, ibeg, iend,
KOKKOS_LAMBDA (int k, int j, int i) {
density(k, j, i) = density(k,j,i) / P(k,j,i);
Expand All @@ -330,11 +314,6 @@ void SelfGravity::InitSolver() {
throw std::runtime_error(msg.str());
}


#ifdef DEBUG_GRAVITY
WriteField(rhoFile, density);
#endif

idfx::popRegion();
}

Expand All @@ -344,15 +323,15 @@ void SelfGravity::SubstractMeanDensity() {

// Loading needed attributes
IdefixArray3D<real> density = this->density;
IdefixArray3D<real> dV = laplacian.dV;
IdefixArray3D<real> dV = laplacian->dV;

int ibeg, iend, jbeg, jend, kbeg, kend;
ibeg = laplacian.beg[IDIR];
iend = laplacian.end[IDIR];
jbeg = laplacian.beg[JDIR];
jend = laplacian.end[JDIR];
kbeg = laplacian.beg[KDIR];
kend = laplacian.end[KDIR];
ibeg = laplacian->beg[IDIR];
iend = laplacian->end[IDIR];
jbeg = laplacian->beg[JDIR];
jend = laplacian->end[JDIR];
kbeg = laplacian->beg[KDIR];
kend = laplacian->end[KDIR];

// Do the reduction on a vector
MyVector meanDensityVector;
Expand Down Expand Up @@ -386,59 +365,22 @@ void SelfGravity::SubstractMeanDensity() {
density(k, j, i) -= mean;
});

#ifdef DEBUG_GRAVITY
idfx::cout << "SelfGravity:: Average value " << mean << " substracted to density" << std::endl;
#endif

idfx::popRegion();
}

void SelfGravity::EnrollUserDefBoundary(Laplacian::UserDefBoundaryFunc myFunc) {
laplacian.EnrollUserDefBoundary(myFunc);
laplacian->EnrollUserDefBoundary(myFunc);
}




#ifdef DEBUG_GRAVITY
// This routine is for debugging purpose
void SelfGravity::WriteField(std::ofstream &stream, IdefixArray3D<real> &in, int index) {
stream.precision(15);
stream << std::scientific;// << index;
for(int i = 0 ; i < this->np_tot[IDIR] ; i++) {
stream << in(0,0,i) << "\t";
}
stream << std::endl;
}

void SelfGravity::WriteField(std::ofstream &stream, IdefixArray1D<real> &in, int index) {
stream.precision(15);
stream << std::scientific;// << index;
for(int i = 0 ; i < this->np_tot[IDIR] ; i++) {
stream << in(i) << "\t";
}
stream << std::endl;
}
#endif

void SelfGravity::SolvePoisson() {
idfx::pushRegion("SelfGravity::SolvePoisson");

Kokkos::Timer timer;

elapsedTime -= timer.seconds();

#ifdef DEBUG_GRAVITY
rhoFile.open("rho.dat",std::ios::trunc);
potentialFile.open("potential.dat",std::ios::trunc);
geometryFile.open("geometry.dat",std::ios::trunc);
WriteField(geometryFile,this->x[IDIR]);
WriteField(geometryFile,this->dx[IDIR]);
WriteField(geometryFile,this->A[IDIR]);
WriteField(geometryFile,this->dV);
geometryFile.close();
#endif

InitSolver(); // (Re)initialise the solver

this->nsteps = iterativeSolver->Solve(potential, density);
Expand Down Expand Up @@ -472,19 +414,6 @@ void SelfGravity::SolvePoisson() {

currentError = iterativeSolver->GetError();

#ifdef DEBUG_GRAVITY
WriteField(potentialFile,potential,n);
if(this->convStatus == true) {
idfx::cout << "SelfGravity:: Reached convergence after " << n << " iterations" << std::endl;
rhoFile.close();
potentialFile.close();
} else if(n == maxiter) {
idfx::cout << "SelfGravity:: Reached max iter" << std::endl;
rhoFile.close();
potentialFile.close();
IDEFIX_WARNING("SelfGravity:: Failed to converge before reaching max iter");
}
#endif

elapsedTime += timer.seconds();
idfx::popRegion();
Expand All @@ -499,12 +428,12 @@ void SelfGravity::AddSelfGravityPotential(IdefixArray3D<real> &phiP) {
real gravCst = this->data->gravity->gravCst;

// Updating ghost cells before to return potential
laplacian.SetBoundaries(potential);
laplacian->SetBoundaries(potential);

// Adding self-gravity contribution
int ioffset = laplacian.loffset[IDIR];
int joffset = laplacian.loffset[JDIR];
int koffset = laplacian.loffset[KDIR];
int ioffset = laplacian->loffset[IDIR];
int joffset = laplacian->loffset[JDIR];
int koffset = laplacian->loffset[KDIR];
idefix_for("AddSelfGravityPotential", 0, data->np_tot[KDIR],
0, data->np_tot[JDIR],
0, data->np_tot[IDIR],
Expand Down
12 changes: 2 additions & 10 deletions src/gravity/selfGravity.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@
#ifndef GRAVITY_SELFGRAVITY_HPP_
#define GRAVITY_SELFGRAVITY_HPP_

#include <memory>
#include <vector>

#include "idefix.hpp"
Expand Down Expand Up @@ -42,16 +43,7 @@ class SelfGravity {
IterativeSolver<Laplacian> *iterativeSolver;

// The linear operator involved in Poisson equation
Laplacian laplacian;

#ifdef DEBUG_GRAVITY
// Used to get fields usefull for debugging
void WriteField(std::ofstream &, IdefixArray3D<real> &, int = 0);
void WriteField(std::ofstream &, IdefixArray1D<real> &, int = 0);
std::ofstream rhoFile;
std::ofstream potentialFile;
std::ofstream geometryFile;
#endif
std::unique_ptr<Laplacian> laplacian;

real currentError{0}; // last error of the iterative solver
int nsteps{0}; // # of steps of the latest iteration
Expand Down
Loading

0 comments on commit 51022cc

Please sign in to comment.