Skip to content

Commit 51022cc

Browse files
glesurglesur
authored
BUG: Fix self gravity (#202)
* 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
1 parent a5d51f5 commit 51022cc

16 files changed

Lines changed: 325 additions & 149 deletions

File tree

‎.github/workflows/idefix-ci.yml‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -186,6 +186,10 @@ jobs:
186186
run: |
187187
cd $IDEFIX_DIR/test/SelfGravity/UniformCollapse
188188
./testme.py -all $TESTME_OPTIONS
189+
- name: Dusty spherical collapse
190+
run: |
191+
cd $IDEFIX_DIR/test/SelfGravity/DustyCollapse
192+
./testme.py -all $TESTME_OPTIONS
189193
190194
Planet:
191195
needs: [Linter, ShocksHydro, ParabolicHydro, ShocksMHD, ParabolicMHD]

‎CHANGELOG.md‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,10 @@ All notable changes to this project will be documented in this file.
44
The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
55
and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
66

7+
## [2.0.1] 2023-10-26
8+
### Changed
9+
- fixed a bug in the self-gravity module introduced in #194 that resulted in instabilities when used in combination with MPI (#202).
10+
711
## [2.0.0] 2023-10-23
812
### Changed
913
- 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).

‎CMakeLists.txt‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,7 @@ set (CMAKE_CXX_STANDARD 17)
44

55
set(Idefix_VERSION_MAJOR 2)
66
set(Idefix_VERSION_MINOR 0)
7-
set(Idefix_VERSION_PATCH 00)
7+
set(Idefix_VERSION_PATCH 01)
88

99
project (idefix VERSION 1.0.0)
1010
option(Idefix_MHD "enable MHD" OFF)

‎src/gravity/laplacian.cpp‎

Lines changed: 30 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,7 @@
1717
Laplacian::Laplacian(DataBlock *datain, std::array<LaplacianBoundaryType,3> leftBound,
1818
std::array<LaplacianBoundaryType,3> rightBound,
1919
bool havePreconditionnerIn) {
20+
idfx::pushRegion("Laplacian::Laplacian");
2021
// Save the parents data objects
2122
this->data = datain;
2223
this->havePreconditioner = havePreconditionnerIn;
@@ -40,10 +41,8 @@ Laplacian::Laplacian(DataBlock *datain, std::array<LaplacianBoundaryType,3> left
4041

4142
isPeriodic = true;
4243
for(int dir = 0 ; dir < 3 ; dir++) {
43-
if(lbound[dir] != LaplacianBoundaryType::periodic &&
44-
lbound[dir] != LaplacianBoundaryType::internalgrav) isPeriodic = false;
45-
if(rbound[dir] != LaplacianBoundaryType::periodic &&
46-
rbound[dir] != LaplacianBoundaryType::internalgrav) isPeriodic = false;
44+
if(lbound[dir] != LaplacianBoundaryType::periodic) isPeriodic = false;
45+
if(rbound[dir] != LaplacianBoundaryType::periodic) isPeriodic = false;
4746
}
4847

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

55+
// Update internal boundaries in case of domain decomposition
56+
for(int dir = 0 ; dir < DIMENSIONS ; dir++) {
57+
if(data->mygrid->nproc[dir]>1) {
58+
if(this->data->lbound[dir]==internal ) {
59+
this->lbound[dir] = Laplacian::internalgrav;
60+
}
61+
if(this->data->rbound[dir]==internal) {
62+
this->rbound[dir] = Laplacian::internalgrav;
63+
}
64+
}
65+
}
66+
#endif
5767

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

75-
// Init MPI stack when needed
76-
#ifdef WITH_MPI
77-
this->arr4D = IdefixArray4D<real> ("WorkingArrayMpi", 1, this->np_tot[KDIR],
78-
this->np_tot[JDIR],
79-
this->np_tot[IDIR]);
80-
81-
int ntarget = 0;
82-
std::vector<int> mapVars;
83-
mapVars.push_back(ntarget);
84-
85-
this->mpi.Init(data->mygrid, mapVars, this->nghost.data(), this->np_int.data());
86-
#endif
87-
88-
8985
if(this->lbound[IDIR] == origin) {
9086
InitInternalGrid();
9187
}
@@ -96,6 +92,21 @@ Laplacian::Laplacian(DataBlock *datain, std::array<LaplacianBoundaryType,3> left
9692
}
9793
// Initialise the Laplacian coefficients
9894
PreComputeLaplacian();
95+
96+
// Init MPI stack when needed
97+
#ifdef WITH_MPI
98+
this->arr4D = IdefixArray4D<real> ("WorkingArrayMpi", 1, this->np_tot[KDIR],
99+
this->np_tot[JDIR],
100+
this->np_tot[IDIR]);
101+
102+
int ntarget = 0;
103+
std::vector<int> mapVars;
104+
mapVars.push_back(ntarget);
105+
106+
this->mpi.Init(data->mygrid, mapVars, this->nghost.data(), this->np_int.data());
107+
#endif
108+
109+
idfx::popRegion();
99110
}
100111

101112
void Laplacian::InitInternalGrid() {

‎src/gravity/selfGravity.cpp‎

Lines changed: 35 additions & 106 deletions
Original file line numberDiff line numberDiff line change
@@ -5,8 +5,7 @@
55
// Licensed under CeCILL 2.1 License, see COPYING for more information
66
// ***********************************************************************************
77

8-
9-
8+
#include <memory>
109
#include <string>
1110
#include <vector>
1211

@@ -106,20 +105,6 @@ void SelfGravity::Init(Input &input, DataBlock *datain) {
106105
}
107106
}
108107

109-
// Update internal boundaries in case of domain decomposition
110-
#ifdef WITH_MPI
111-
for(int dir = 0 ; dir < DIMENSIONS ; dir++) {
112-
if(data->mygrid->nproc[dir]>1) {
113-
if(this->data->lbound[dir]==internal ) {
114-
this->lbound[dir] = Laplacian::internalgrav;
115-
}
116-
if(this->data->rbound[dir]==internal) {
117-
this->rbound[dir] = Laplacian::internalgrav;
118-
}
119-
}
120-
}
121-
#endif
122-
123108
// Update solver when provided
124109
if(input.CheckEntry("SelfGravity","solver") >= 0) {
125110
std::string strSolver = input.Get<std::string>("SelfGravity","solver",0);
@@ -162,29 +147,28 @@ void SelfGravity::Init(Input &input, DataBlock *datain) {
162147
}
163148

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

167-
np_tot = laplacian.np_tot;
152+
np_tot = laplacian->np_tot;
168153

169154
// Instantiate the bicgstab solver
170155
if(solver == BICGSTAB || solver == PBICGSTAB) {
171-
iterativeSolver = new Bicgstab(laplacian, targetError, maxiter,
172-
laplacian.np_tot, laplacian.beg, laplacian.end);
156+
iterativeSolver = new Bicgstab(*laplacian.get(), targetError, maxiter,
157+
laplacian->np_tot, laplacian->beg, laplacian->end);
173158
} else if(solver == CG || solver == PCG) {
174-
iterativeSolver = new Cg(laplacian, targetError, maxiter,
175-
laplacian.np_tot, laplacian.beg, laplacian.end);
159+
iterativeSolver = new Cg(*laplacian.get(), targetError, maxiter,
160+
laplacian->np_tot, laplacian->beg, laplacian->end);
176161
} else if(solver == MINRES || solver == PMINRES) {
177-
iterativeSolver = new Minres(laplacian,
162+
iterativeSolver = new Minres(*laplacian.get(),
178163
targetError, maxiter,
179-
laplacian.np_tot, laplacian.beg, laplacian.end);
164+
laplacian->np_tot, laplacian->beg, laplacian->end);
180165
} else {
181-
real step = laplacian.ComputeCFL();
182-
iterativeSolver = new Jacobi(laplacian, targetError, maxiter, step,
183-
laplacian.np_tot, laplacian.beg, laplacian.end);
166+
real step = laplacian->ComputeCFL();
167+
iterativeSolver = new Jacobi(*laplacian.get(), targetError, maxiter, step,
168+
laplacian->np_tot, laplacian->beg, laplacian->end);
184169
}
185170

186171

187-
188172
// Arrays initialisation
189173
this->density = IdefixArray3D<real> ("Density", this->np_tot[KDIR],
190174
this->np_tot[JDIR],
@@ -246,7 +230,7 @@ void SelfGravity::ShowConfig() {
246230
}
247231

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

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

269253
// Initialise the density field
270254
// todo: check bounds
271-
int ioffset = laplacian.loffset[IDIR];
272-
int joffset = laplacian.loffset[JDIR];
273-
int koffset = laplacian.loffset[KDIR];
255+
int ioffset = laplacian->loffset[IDIR];
256+
int joffset = laplacian->loffset[JDIR];
257+
int koffset = laplacian->loffset[KDIR];
274258

275259
idefix_for("InitDensity", data->beg[KDIR], data->end[KDIR],
276260
data->beg[JDIR], data->end[JDIR],
@@ -300,13 +284,13 @@ void SelfGravity::InitSolver() {
300284
// divide density by preconditionner if we're doing the preconditionned version
301285
if(havePreconditioner) {
302286
int ibeg, iend, jbeg, jend, kbeg, kend;
303-
ibeg = laplacian.beg[IDIR];
304-
iend = laplacian.end[IDIR];
305-
jbeg = laplacian.beg[JDIR];
306-
jend = laplacian.end[JDIR];
307-
kbeg = laplacian.beg[KDIR];
308-
kend = laplacian.end[KDIR];
309-
IdefixArray3D<real> P = laplacian.precond;
287+
ibeg = laplacian->beg[IDIR];
288+
iend = laplacian->end[IDIR];
289+
jbeg = laplacian->beg[JDIR];
290+
jend = laplacian->end[JDIR];
291+
kbeg = laplacian->beg[KDIR];
292+
kend = laplacian->end[KDIR];
293+
IdefixArray3D<real> P = laplacian->precond;
310294
idefix_for("Precond density", kbeg, kend, jbeg, jend, ibeg, iend,
311295
KOKKOS_LAMBDA (int k, int j, int i) {
312296
density(k, j, i) = density(k,j,i) / P(k,j,i);
@@ -330,11 +314,6 @@ void SelfGravity::InitSolver() {
330314
throw std::runtime_error(msg.str());
331315
}
332316

333-
334-
#ifdef DEBUG_GRAVITY
335-
WriteField(rhoFile, density);
336-
#endif
337-
338317
idfx::popRegion();
339318
}
340319

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

345324
// Loading needed attributes
346325
IdefixArray3D<real> density = this->density;
347-
IdefixArray3D<real> dV = laplacian.dV;
326+
IdefixArray3D<real> dV = laplacian->dV;
348327

349328
int ibeg, iend, jbeg, jend, kbeg, kend;
350-
ibeg = laplacian.beg[IDIR];
351-
iend = laplacian.end[IDIR];
352-
jbeg = laplacian.beg[JDIR];
353-
jend = laplacian.end[JDIR];
354-
kbeg = laplacian.beg[KDIR];
355-
kend = laplacian.end[KDIR];
329+
ibeg = laplacian->beg[IDIR];
330+
iend = laplacian->end[IDIR];
331+
jbeg = laplacian->beg[JDIR];
332+
jend = laplacian->end[JDIR];
333+
kbeg = laplacian->beg[KDIR];
334+
kend = laplacian->end[KDIR];
356335

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

389-
#ifdef DEBUG_GRAVITY
390-
idfx::cout << "SelfGravity:: Average value " << mean << " substracted to density" << std::endl;
391-
#endif
392-
393368
idfx::popRegion();
394369
}
395370

396371
void SelfGravity::EnrollUserDefBoundary(Laplacian::UserDefBoundaryFunc myFunc) {
397-
laplacian.EnrollUserDefBoundary(myFunc);
372+
laplacian->EnrollUserDefBoundary(myFunc);
398373
}
399374

400375

401376

402-
403-
#ifdef DEBUG_GRAVITY
404-
// This routine is for debugging purpose
405-
void SelfGravity::WriteField(std::ofstream &stream, IdefixArray3D<real> &in, int index) {
406-
stream.precision(15);
407-
stream << std::scientific;// << index;
408-
for(int i = 0 ; i < this->np_tot[IDIR] ; i++) {
409-
stream << in(0,0,i) << "\t";
410-
}
411-
stream << std::endl;
412-
}
413-
414-
void SelfGravity::WriteField(std::ofstream &stream, IdefixArray1D<real> &in, int index) {
415-
stream.precision(15);
416-
stream << std::scientific;// << index;
417-
for(int i = 0 ; i < this->np_tot[IDIR] ; i++) {
418-
stream << in(i) << "\t";
419-
}
420-
stream << std::endl;
421-
}
422-
#endif
423-
424377
void SelfGravity::SolvePoisson() {
425378
idfx::pushRegion("SelfGravity::SolvePoisson");
426379

427380
Kokkos::Timer timer;
428381

429382
elapsedTime -= timer.seconds();
430383

431-
#ifdef DEBUG_GRAVITY
432-
rhoFile.open("rho.dat",std::ios::trunc);
433-
potentialFile.open("potential.dat",std::ios::trunc);
434-
geometryFile.open("geometry.dat",std::ios::trunc);
435-
WriteField(geometryFile,this->x[IDIR]);
436-
WriteField(geometryFile,this->dx[IDIR]);
437-
WriteField(geometryFile,this->A[IDIR]);
438-
WriteField(geometryFile,this->dV);
439-
geometryFile.close();
440-
#endif
441-
442384
InitSolver(); // (Re)initialise the solver
443385

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

473415
currentError = iterativeSolver->GetError();
474416

475-
#ifdef DEBUG_GRAVITY
476-
WriteField(potentialFile,potential,n);
477-
if(this->convStatus == true) {
478-
idfx::cout << "SelfGravity:: Reached convergence after " << n << " iterations" << std::endl;
479-
rhoFile.close();
480-
potentialFile.close();
481-
} else if(n == maxiter) {
482-
idfx::cout << "SelfGravity:: Reached max iter" << std::endl;
483-
rhoFile.close();
484-
potentialFile.close();
485-
IDEFIX_WARNING("SelfGravity:: Failed to converge before reaching max iter");
486-
}
487-
#endif
488417

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

501430
// Updating ghost cells before to return potential
502-
laplacian.SetBoundaries(potential);
431+
laplacian->SetBoundaries(potential);
503432

504433
// Adding self-gravity contribution
505-
int ioffset = laplacian.loffset[IDIR];
506-
int joffset = laplacian.loffset[JDIR];
507-
int koffset = laplacian.loffset[KDIR];
434+
int ioffset = laplacian->loffset[IDIR];
435+
int joffset = laplacian->loffset[JDIR];
436+
int koffset = laplacian->loffset[KDIR];
508437
idefix_for("AddSelfGravityPotential", 0, data->np_tot[KDIR],
509438
0, data->np_tot[JDIR],
510439
0, data->np_tot[IDIR],

‎src/gravity/selfGravity.hpp‎

Lines changed: 2 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@
88
#ifndef GRAVITY_SELFGRAVITY_HPP_
99
#define GRAVITY_SELFGRAVITY_HPP_
1010

11+
#include <memory>
1112
#include <vector>
1213

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

4445
// The linear operator involved in Poisson equation
45-
Laplacian laplacian;
46-
47-
#ifdef DEBUG_GRAVITY
48-
// Used to get fields usefull for debugging
49-
void WriteField(std::ofstream &, IdefixArray3D<real> &, int = 0);
50-
void WriteField(std::ofstream &, IdefixArray1D<real> &, int = 0);
51-
std::ofstream rhoFile;
52-
std::ofstream potentialFile;
53-
std::ofstream geometryFile;
54-
#endif
46+
std::unique_ptr<Laplacian> laplacian;
5547

5648
real currentError{0}; // last error of the iterative solver
5749
int nsteps{0}; // # of steps of the latest iteration

0 commit comments

Comments
 (0)