From 51022cc32044c6b818ce1dc37445f35fd219c50b Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Thu, 26 Oct 2023 09:30:04 +0200 Subject: [PATCH] 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 --- .github/workflows/idefix-ci.yml | 4 + CHANGELOG.md | 4 + CMakeLists.txt | 2 +- src/gravity/laplacian.cpp | 49 +++--- src/gravity/selfGravity.cpp | 141 +++++------------- src/gravity/selfGravity.hpp | 12 +- src/mpi.cpp | 5 + src/mpi.hpp | 7 + src/utils/iterativesolver/iterativesolver.hpp | 27 ++-- .../SelfGravity/DustyCollapse/definitions.hpp | 5 + test/SelfGravity/DustyCollapse/idefix.ini | 53 +++++++ .../DustyCollapse/python/testidefix.py | 43 ++++++ test/SelfGravity/DustyCollapse/setup.cpp | 66 ++++++++ test/SelfGravity/DustyCollapse/testme.py | 43 ++++++ test/SelfGravity/JeansInstability/testme.py | 6 + test/SelfGravity/RandomSphere/testme.py | 7 + 16 files changed, 325 insertions(+), 149 deletions(-) create mode 100644 test/SelfGravity/DustyCollapse/definitions.hpp create mode 100644 test/SelfGravity/DustyCollapse/idefix.ini create mode 100644 test/SelfGravity/DustyCollapse/python/testidefix.py create mode 100644 test/SelfGravity/DustyCollapse/setup.cpp create mode 100755 test/SelfGravity/DustyCollapse/testme.py diff --git a/.github/workflows/idefix-ci.yml b/.github/workflows/idefix-ci.yml index 61e7d293..199626e0 100644 --- a/.github/workflows/idefix-ci.yml +++ b/.github/workflows/idefix-ci.yml @@ -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] diff --git a/CHANGELOG.md b/CHANGELOG.md index 7b8bdf42..739c1a09 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -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). diff --git a/CMakeLists.txt b/CMakeLists.txt index 7777bbde..f3fa42e0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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) diff --git a/src/gravity/laplacian.cpp b/src/gravity/laplacian.cpp index 8953eb4e..2300db82 100644 --- a/src/gravity/laplacian.cpp +++ b/src/gravity/laplacian.cpp @@ -17,6 +17,7 @@ Laplacian::Laplacian(DataBlock *datain, std::array leftBound, std::array rightBound, bool havePreconditionnerIn) { + idfx::pushRegion("Laplacian::Laplacian"); // Save the parents data objects this->data = datain; this->havePreconditioner = havePreconditionnerIn; @@ -40,10 +41,8 @@ Laplacian::Laplacian(DataBlock *datain, std::array 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 @@ -52,8 +51,19 @@ Laplacian::Laplacian(DataBlock *datain, std::array 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)) { @@ -72,20 +82,6 @@ Laplacian::Laplacian(DataBlock *datain, std::array left } #endif - // Init MPI stack when needed - #ifdef WITH_MPI - this->arr4D = IdefixArray4D ("WorkingArrayMpi", 1, this->np_tot[KDIR], - this->np_tot[JDIR], - this->np_tot[IDIR]); - - int ntarget = 0; - std::vector 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(); } @@ -96,6 +92,21 @@ Laplacian::Laplacian(DataBlock *datain, std::array left } // Initialise the Laplacian coefficients PreComputeLaplacian(); + + // Init MPI stack when needed + #ifdef WITH_MPI + this->arr4D = IdefixArray4D ("WorkingArrayMpi", 1, this->np_tot[KDIR], + this->np_tot[JDIR], + this->np_tot[IDIR]); + + int ntarget = 0; + std::vector mapVars; + mapVars.push_back(ntarget); + + this->mpi.Init(data->mygrid, mapVars, this->nghost.data(), this->np_int.data()); + #endif + + idfx::popRegion(); } void Laplacian::InitInternalGrid() { diff --git a/src/gravity/selfGravity.cpp b/src/gravity/selfGravity.cpp index 67c0f834..661b7f59 100644 --- a/src/gravity/selfGravity.cpp +++ b/src/gravity/selfGravity.cpp @@ -5,8 +5,7 @@ // Licensed under CeCILL 2.1 License, see COPYING for more information // *********************************************************************************** - - +#include #include #include @@ -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("SelfGravity","solver",0); @@ -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(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 ("Density", this->np_tot[KDIR], this->np_tot[JDIR], @@ -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; } @@ -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], @@ -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 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 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); @@ -330,11 +314,6 @@ void SelfGravity::InitSolver() { throw std::runtime_error(msg.str()); } - - #ifdef DEBUG_GRAVITY - WriteField(rhoFile, density); - #endif - idfx::popRegion(); } @@ -344,15 +323,15 @@ void SelfGravity::SubstractMeanDensity() { // Loading needed attributes IdefixArray3D density = this->density; - IdefixArray3D dV = laplacian.dV; + IdefixArray3D 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; @@ -386,41 +365,15 @@ 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 &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 &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"); @@ -428,17 +381,6 @@ void SelfGravity::SolvePoisson() { 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); @@ -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(); @@ -499,12 +428,12 @@ void SelfGravity::AddSelfGravityPotential(IdefixArray3D &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], diff --git a/src/gravity/selfGravity.hpp b/src/gravity/selfGravity.hpp index f38b2f5f..badd729c 100644 --- a/src/gravity/selfGravity.hpp +++ b/src/gravity/selfGravity.hpp @@ -8,6 +8,7 @@ #ifndef GRAVITY_SELFGRAVITY_HPP_ #define GRAVITY_SELFGRAVITY_HPP_ +#include #include #include "idefix.hpp" @@ -42,16 +43,7 @@ class SelfGravity { IterativeSolver *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 &, int = 0); - void WriteField(std::ofstream &, IdefixArray1D &, int = 0); - std::ofstream rhoFile; - std::ofstream potentialFile; - std::ofstream geometryFile; - #endif + std::unique_ptr laplacian; real currentError{0}; // last error of the iterative solver int nsteps{0}; // # of steps of the latest iteration diff --git a/src/mpi.cpp b/src/mpi.cpp index 9ca66b68..ee2fdb5c 100644 --- a/src/mpi.cpp +++ b/src/mpi.cpp @@ -45,6 +45,7 @@ void Mpi::ExchangeAll() { void Mpi::Init(Grid *grid, std::vector inputMap, int nghost[3], int nint[3], bool inputHaveVs) { + idfx::pushRegion("Mpi::Init"); this->mygrid = grid; // increase the number of instances @@ -202,10 +203,13 @@ void Mpi::Init(Grid *grid, std::vector inputMap, // say this instance is initialized. isInitialized = true; + + idfx::popRegion(); } // Destructor (clean up persistent communication channels) Mpi::~Mpi() { + idfx::pushRegion("Mpi::~Mpi"); if(isInitialized) { // Properly clean up the mess #ifdef MPI_PERSISTENT @@ -236,6 +240,7 @@ Mpi::~Mpi() { } isInitialized = false; } + idfx::popRegion(); } void Mpi::ExchangeX1(IdefixArray4D Vc, IdefixArray4D Vs) { diff --git a/src/mpi.hpp b/src/mpi.hpp index 03aceb03..a4250a17 100644 --- a/src/mpi.hpp +++ b/src/mpi.hpp @@ -183,6 +183,7 @@ class Buffer { class Mpi { public: + Mpi() = default; // MPI Exchange functions void ExchangeAll(); ///< Exchange boundary elements in all directions (todo) void ExchangeX1(IdefixArray4D inputVc, @@ -210,8 +211,14 @@ class Mpi { ~Mpi(); private: + // Because the MPI class initialise internal pointers, we do not allow copies of this class + // These lines should not be removed as they constitute a safeguard + Mpi(const Mpi&); + Mpi operator=(const Mpi&); + static int nInstances; // total number of mpi instances in the code int thisInstance; // unique number of the current instance + int nReferences; // # of references to this instance bool isInitialized{false}; DataBlock *data; // pointer to datablock object diff --git a/src/utils/iterativesolver/iterativesolver.hpp b/src/utils/iterativesolver/iterativesolver.hpp index d103eba6..e128f7ff 100644 --- a/src/utils/iterativesolver/iterativesolver.hpp +++ b/src/utils/iterativesolver/iterativesolver.hpp @@ -37,6 +37,7 @@ class IterativeSolver { int maxiter; // Maximum iteration allowed to achieve convergence bool convStatus; // Convergence status bool restart{false}; + static constexpr bool isVerbose{false}; // Whether the solver should be verbose while iterating std::array beg; std::array end; @@ -112,10 +113,10 @@ void IterativeSolver::TestErrorL1() { // Checking convergence if(currentError <= this->targetError) { this->convStatus = true; - #ifdef DEBUG_BICGSTAB - idfx::cout << "IterativeSolver:: Squared, normalized norm L1 is " << currentError - << " at convergence." << std::endl; - #endif + if constexpr(isVerbose) { + idfx::cout << "IterativeSolver:: Squared, normalized norm L1 is " << currentError + << " at convergence." << std::endl; + } } idfx::popRegion(); @@ -159,7 +160,7 @@ void IterativeSolver::TestErrorL2() { // Squared error this->currentError = sqrt(normL2Vector.v[0] / normL2Vector.v[1]); - //idfx::cout << "Error=" << this->currentError << std::endl; + if constexpr(isVerbose) idfx::cout << "L2 Error=" << this->currentError << std::endl; // Checking Nans if(std::isnan(currentError)) { @@ -170,10 +171,10 @@ void IterativeSolver::TestErrorL2() { // Checking convergence if(currentError <= this->targetError) { this->convStatus = true; - #ifdef DEBUG_BICGSTAB - idfx::cout << "IterativeSolver:: Squared, normalized norm L2 is " << currentError - << " at convergence." << std::endl; - #endif + if constexpr(isVerbose) { + idfx::cout << "IterativeSolver:: Squared, normalized norm L2 is " << currentError + << " at convergence." << std::endl; + } } idfx::popRegion(); @@ -235,10 +236,10 @@ void IterativeSolver::TestErrorLINF() { // Checking convergence if(currentError <= this->targetError) { this->convStatus = true; - #ifdef DEBUG_BICGSTAB - idfx::cout << "IterativeSolver:: Squared, normalized norm LINF is " << currentError - << " at convergence." << std::endl; - #endif + if constexpr(this->isVerbose) { + idfx::cout << "IterativeSolver:: Squared, normalized norm LINF is " << currentError + << " at convergence." << std::endl; + } } idfx::popRegion(); diff --git a/test/SelfGravity/DustyCollapse/definitions.hpp b/test/SelfGravity/DustyCollapse/definitions.hpp new file mode 100644 index 00000000..77affc7e --- /dev/null +++ b/test/SelfGravity/DustyCollapse/definitions.hpp @@ -0,0 +1,5 @@ +#define COMPONENTS 1 +#define DIMENSIONS 1 +#define ISOTHERMAL + +#define GEOMETRY SPHERICAL diff --git a/test/SelfGravity/DustyCollapse/idefix.ini b/test/SelfGravity/DustyCollapse/idefix.ini new file mode 100644 index 00000000..37976368 --- /dev/null +++ b/test/SelfGravity/DustyCollapse/idefix.ini @@ -0,0 +1,53 @@ +[Grid] +X1-grid 1 0.1 4096 l 3e4 +X2-grid 1 0 64 u 3.14159265358979 +X3-grid 1 0 64 u 6.28318530717958 + +[TimeIntegrator] +CFL 0.2 +tstop 1 +nstages 1 +fixed_dt 1 + +[Hydro] +solver hllc +csiso constant 10.0 + +[Gravity] +potential selfgravity central +Mcentral 1e2 +gravCst 6.67430e-11 + +[SelfGravity] +solver PBICGSTAB +skip 100 +targetError 1e-7 +maxIter 10000 +boundary-X1-beg nullgrad +boundary-X1-end nullpot +boundary-X2-beg axis +boundary-X2-end axis +boundary-X3-beg periodic +boundary-X3-end periodic + +[Boundary] +X1-beg outflow +X1-end outflow +X2-beg periodic +X2-end periodic +X3-beg periodic +X3-end periodic + +[Output] +dmp 1 +vtk 1 +uservar phiP + +[Dust] +nSpecies 1 +drag tau 1.0 +drag_feedback yes + +[Setup] +Mtot 1e16 # particle cloud mass (kg) +R0 1e4 # particle cloud radius (m) diff --git a/test/SelfGravity/DustyCollapse/python/testidefix.py b/test/SelfGravity/DustyCollapse/python/testidefix.py new file mode 100644 index 00000000..370d89e7 --- /dev/null +++ b/test/SelfGravity/DustyCollapse/python/testidefix.py @@ -0,0 +1,43 @@ +import numpy as np +import os +import sys + +sys.path.append(os.getenv("IDEFIX_DIR")) +from pytools.vtk_io import readVTK + +G=6.6743e-11 +M=1e16 #mass of spheric clump +r0=1e4 #radius of spheric clump +rmax = 3e4 #right edge of simulation box +d=M/(4*np.pi*r0**3/3) #density of spheric clump +phi0 = -22.247666667 #gravitationnal potential at rmax + +def potential(r): + retv = np.empty_like(r) + interior = (r + +real MTOT; // total mass of spheric particle clump +real R0; // size of spheric particle clump + + +void ComputeUserVars(DataBlock & data, UserDefVariablesContainer &variables) { + // Make references to the user-defined arrays (variables is a container of IdefixHostArray3D) + // Note that the labels should match the variable names in the input file + IdefixHostArray3D phiP = variables["phiP"]; + Kokkos::deep_copy(phiP, data.gravity->phiP); +} + +Setup::Setup(Input &input, Grid &grid, DataBlock &data, Output &output) { + MTOT = input.Get("Setup", "Mtot", 0); + R0 = input.Get("Setup", "R0", 0); + + output.EnrollUserDefVariables(&ComputeUserVars); +} + + +real rho_profile (real r) { // density profile of spheric particle clump + real rho_part = MTOT*3.0/(4.0*M_PI*R0*R0*R0); + return (r < R0 ? rho_part : rho_part/10000); +} + + + +void Setup::InitFlow(DataBlock &data) { + DataBlockHost d(data); + d.SyncFromDevice(); + + IdefixHostArray1D r = d.x[IDIR]; + + // init flow + for(int k = 0; k < d.np_tot[KDIR] ; k++) { + for(int j = 0; j < d.np_tot[JDIR] ; j++) { + for(int i = 0; i < d.np_tot[IDIR] ; i++) { + d.Vc(RHO,k,j,i) = 1e-6; + EXPAND(\ + d.Vc(VX1,k,j,i) = ZERO_F; ,\ + d.Vc(VX2,k,j,i) = ZERO_F; ,\ + d.Vc(VX3,k,j,i) = ZERO_F; ) + + d.dustVc[0](RHO,k,j,i) = rho_profile(r(i)); + EXPAND(\ + d.dustVc[0](VX1,k,j,i) = ZERO_F; ,\ + d.dustVc[0](VX2,k,j,i) = ZERO_F; ,\ + d.dustVc[0](VX3,k,j,i) = ZERO_F; ) + + } + } + } + + + + d.SyncToDevice(); +} + +// Analyse data to produce an output +void MakeAnalysis(DataBlock & data) { + +} diff --git a/test/SelfGravity/DustyCollapse/testme.py b/test/SelfGravity/DustyCollapse/testme.py new file mode 100755 index 00000000..68903826 --- /dev/null +++ b/test/SelfGravity/DustyCollapse/testme.py @@ -0,0 +1,43 @@ +#!/usr/bin/env python3 + +""" + +@author: glesur +""" +import os +import sys +sys.path.append(os.getenv("IDEFIX_DIR")) + +import pytools.idfx_test as tst + +name="dump.0001.dmp" + +tolerance=1e-13 + +def testMe(test): + test.configure() + test.compile() + inifiles=["idefix.ini"] + + # loop on all the ini files for this test + for ini in inifiles: + test.run(inputFile=ini) + test.standardTest() + + +test=tst.idfxTest() + +if not test.all: + if(test.check): + test.standardTest() + else: + testMe(test) +else: + test.noplot = True + test.vectPot=False + test.single=False + test.reconstruction=2 + test.mpi=False + testMe(test) + test.mpi=True + testMe(test) diff --git a/test/SelfGravity/JeansInstability/testme.py b/test/SelfGravity/JeansInstability/testme.py index 435156d1..ff2f39a8 100755 --- a/test/SelfGravity/JeansInstability/testme.py +++ b/test/SelfGravity/JeansInstability/testme.py @@ -28,6 +28,10 @@ def testMe(test): test=tst.idfxTest() +# if no decomposition is specified, use that one +if not test.dec: + test.dec=['2'] + if not test.all: if(test.check): test.standardTest() @@ -40,3 +44,5 @@ def testMe(test): test.reconstruction=2 test.mpi=False testMe(test) + test.mpi=True + testMe(test) diff --git a/test/SelfGravity/RandomSphere/testme.py b/test/SelfGravity/RandomSphere/testme.py index beafdd77..38fc30d9 100755 --- a/test/SelfGravity/RandomSphere/testme.py +++ b/test/SelfGravity/RandomSphere/testme.py @@ -27,8 +27,15 @@ def testMe(test): test=tst.idfxTest() + +if not test.dec: + test.dec=['2','2','1'] + if not test.all: testMe(test) else: + test.mpi=False test.noplot=True testMe(test) + test.mpi=True + testMe(test)