From 496a6c30b148e5ccd7a405f5375c404de49f07be Mon Sep 17 00:00:00 2001 From: marc Date: Fri, 25 Oct 2024 10:37:45 +0200 Subject: [PATCH 01/28] Almost works --- src/gravity/laplacian.cpp | 101 ++++++++++++++++++++++++++++++++++-- src/gravity/laplacian.hpp | 1 + src/gravity/selfGravity.cpp | 6 +++ 3 files changed, 105 insertions(+), 3 deletions(-) diff --git a/src/gravity/laplacian.cpp b/src/gravity/laplacian.cpp index 2300db82d..f3b3644a8 100644 --- a/src/gravity/laplacian.cpp +++ b/src/gravity/laplacian.cpp @@ -12,6 +12,7 @@ #include "laplacian.hpp" #include "selfGravity.hpp" #include "dataBlock.hpp" +#include "fluid.hpp" Laplacian::Laplacian(DataBlock *datain, std::array leftBound, @@ -39,10 +40,12 @@ Laplacian::Laplacian(DataBlock *datain, std::array left this->lbound = leftBound; this->rbound = rightBound; - isPeriodic = true; + isPeriodic = true; // is this ever used ? this is not selfgravity.isPeriodic for(int dir = 0 ; dir < 3 ; dir++) { - if(lbound[dir] != LaplacianBoundaryType::periodic) isPeriodic = false; - if(rbound[dir] != LaplacianBoundaryType::periodic) isPeriodic = false; + if(lbound[dir] != LaplacianBoundaryType::periodic + && lbound[dir] != LaplacianBoundaryType::shearingbox) isPeriodic = false; + if(rbound[dir] != LaplacianBoundaryType::periodic + && lbound[dir] != LaplacianBoundaryType::shearingbox) isPeriodic = false; } #ifdef WITH_MPI @@ -521,7 +524,99 @@ void Laplacian::EnforceBoundary(int dir, BoundarySide side, LaplacianBoundaryTyp }); break; } + case shearingbox: { + const real S = data->hydro->sbS; + // Box size + const real Lx = data->mygrid->xend[IDIR] - data->mygrid->xbeg[IDIR]; + const real Ly = data->mygrid->xend[JDIR] - data->mygrid->xbeg[JDIR]; + + // total number of cells in y (active domain) + const int nxj = data->mygrid->np_int[JDIR]; + const real dy = Ly/nxj; + + // Compute offset in y modulo the box size + const real t = data->t; + + const int sign=2*side-1; + const real sbVelocity = sign*S*Lx; + real dL = std::fmod(sbVelocity*t,Ly); + + const int m = static_cast (std::floor(dL/dy+HALF_F)); + + //const real eps = dL - m *dy; + const real eps = dL/dy - m; + + const int sideShift =0;// (side +1)%2; + + idefix_for("BoundaryShearingBox", kbeg, kend, jbeg, jend, ibeg, iend, + KOKKOS_LAMBDA (int k, int j, int i) { + int iref, jref, kref; + // This hack takes care of cases where we have more ghost zones than active zones + + if(dir==IDIR) + iref = ighost + (i+ighost*(nxi-1))%nxi; + else + IDEFIX_ERROR( + "Laplacian:: Shearing box boundary condition should only be used in IDIR" + ); + + //localVar(k,j,i) = localVar(k,j,iref)+j; // this works + + // il faut prendre en compte les ghost zones d'une façon ou d'une autre + // ce n'est pas le cas ici... + + const int jo = jghost + ((j-m-jghost+sideShift)%nxj+nxj)%nxj; + const int jop1 = jghost + ((jo+1-jghost)%nxj+nxj)%nxj; + const int jom1 = jghost + ((jo-1-jghost)%nxj+nxj)%nxj; + + const int jom2 = jghost + ((jo-2-jghost)%nxj+nxj)%nxj; + const int jop2 = jghost + ((jo+2-jghost)%nxj+nxj)%nxj; + + real Fl,Fr; + real dqm, dqp, dq; +/* + if (eps <= 0) { + localVar(k,j,i) = (1+eps) * localVar(k,jo,iref) + - eps*localVar(k,jop1,iref); + } else { + localVar(k,j,i) = ((1-eps) * localVar(k,jom1,iref) + + eps*localVar(k,jo,iref)); + }*/ + + // the follwing, like for the fluxes in boundary.hpp is a second- + // order recontruction in eps + if(eps>=ZERO_F) { + // Compute Fl + dqm = localVar(k,jom1,iref) - localVar(k,jom2,iref); + dqp = localVar(k,jo,iref) - localVar(k,jom1,iref); + dq = dqm+dqp; //(dqp*dqm > ZERO_F ? TWO_F*dqp*dqm/(dqp + dqm) : ZERO_F); + + Fl = localVar(k,jom1,iref) + 0.5*dq*(1.0-eps); + //Compute Fr + dqm=dqp; + dqp = localVar(k,jop1,iref) - localVar(k,jo,iref); + dq = dqm+dqp; //(dqp*dqm > ZERO_F ? TWO_F*dqp*dqm/(dqp + dqm) : ZERO_F); + + Fr = localVar(k,jo,iref) + 0.5*dq*(1.0-eps); + } else { + //Compute Fl + dqm = localVar(k,jo,iref) - localVar(k,jom1,iref); + dqp = localVar(k,jop1,iref) - localVar(k,jo,iref); + dq = dqm+dqp; //(dqp*dqm > ZERO_F ? TWO_F*dqp*dqm/(dqp + dqm) : ZERO_F); + + Fl = localVar(k,jo,iref) - 0.5*dq*(1.0+eps); + // Compute Fr + dqm=dqp; + dqp = localVar(k,jop2,iref) - localVar(k,jop1,iref); + dq = dqm+dqp; //(dqp*dqm > ZERO_F ? TWO_F*dqp*dqm/(dqp + dqm) : ZERO_F); + + Fr = localVar(k,jop1,iref) - 0.5*dq*(1.0+eps); + } + localVar(k,j,i) = localVar(k,jo,iref) - eps*(Fr - Fl); + }); + break; + } case userdef: { if(this->haveUserDefBoundary) { // Warning: unlike hydro userdef boundary functions, the selfGravity diff --git a/src/gravity/laplacian.hpp b/src/gravity/laplacian.hpp index 285724399..b70b79a80 100644 --- a/src/gravity/laplacian.hpp +++ b/src/gravity/laplacian.hpp @@ -21,6 +21,7 @@ class Laplacian { // Types of boundary which can be treated enum LaplacianBoundaryType {internalgrav, periodic, + shearingbox, nullgrad, nullpot, userdef, diff --git a/src/gravity/selfGravity.cpp b/src/gravity/selfGravity.cpp index 43f764713..4618b57b4 100644 --- a/src/gravity/selfGravity.cpp +++ b/src/gravity/selfGravity.cpp @@ -55,6 +55,9 @@ void SelfGravity::Init(Input &input, DataBlock *datain) { this->isPeriodic = false; } else if(boundary.compare("periodic") == 0) { this->lbound[dir] = Laplacian::LaplacianBoundaryType::periodic; + } else if(boundary.compare("shearingbox") == 0) { + this->lbound[dir] = Laplacian::LaplacianBoundaryType::shearingbox; + // this is considered periodic so that mean density is substracted } else if(boundary.compare("nullgrad") == 0) { this->lbound[dir] = Laplacian::LaplacianBoundaryType::nullgrad; this->isPeriodic = false; @@ -90,6 +93,9 @@ void SelfGravity::Init(Input &input, DataBlock *datain) { this->isPeriodic = false; } else if(boundary.compare("periodic") == 0) { this->rbound[dir] = Laplacian::LaplacianBoundaryType::periodic; + } else if(boundary.compare("shearingbox") == 0) { + this->rbound[dir] = Laplacian::LaplacianBoundaryType::shearingbox; + // this is considered periodic so that mean density is substracted } else if(boundary.compare("nullgrad") == 0) { this->rbound[dir] = Laplacian::LaplacianBoundaryType::nullgrad; this->isPeriodic = false; From 62bce3ab3b95c2bc74e8ffd416cf3de18bf724a0 Mon Sep 17 00:00:00 2001 From: marc Date: Fri, 15 Nov 2024 10:29:22 +0100 Subject: [PATCH 02/28] test enforce periodic From 091aed2a74374b54b64372d936d536e7488a3dca Mon Sep 17 00:00:00 2001 From: marc Date: Fri, 15 Nov 2024 11:34:02 +0100 Subject: [PATCH 03/28] Tiref test --- src/gravity/laplacian.cpp | 45 ++++++++++++++++++++++++++++----------- src/gravity/laplacian.hpp | 2 ++ 2 files changed, 35 insertions(+), 12 deletions(-) diff --git a/src/gravity/laplacian.cpp b/src/gravity/laplacian.cpp index f3b3644a8..dc66a8564 100644 --- a/src/gravity/laplacian.cpp +++ b/src/gravity/laplacian.cpp @@ -48,6 +48,15 @@ Laplacian::Laplacian(DataBlock *datain, std::array left && lbound[dir] != LaplacianBoundaryType::shearingbox) isPeriodic = false; } + #if GEOMETRY == CARTESIAN + if(lbound[IDIR] == shearingbox || rbound[IDIR] == shearingbox) { + sBArray = IdefixArray3D("ShearingPotentialBoxArray", + data->np_tot[KDIR], + data->np_tot[JDIR], + data->nghost[IDIR]); + } + #endif // CARTESIAN + #ifdef WITH_MPI if(lbound[IDIR] == origin) { // create communicator for spherical radius @@ -525,6 +534,8 @@ void Laplacian::EnforceBoundary(int dir, BoundarySide side, LaplacianBoundaryTyp break; } case shearingbox: { + IdefixArray3D scrh = sBArray; + const real S = data->hydro->sbS; // Box size @@ -549,17 +560,20 @@ void Laplacian::EnforceBoundary(int dir, BoundarySide side, LaplacianBoundaryTyp const int sideShift =0;// (side +1)%2; + const int nxiglob = data->mygrid->np_int[IDIR]; + idefix_for("BoundaryShearingBox", kbeg, kend, jbeg, jend, ibeg, iend, KOKKOS_LAMBDA (int k, int j, int i) { int iref, jref, kref; // This hack takes care of cases where we have more ghost zones than active zones if(dir==IDIR) - iref = ighost + (i+ighost*(nxi-1))%nxi; + iref = i - side*(ighost +nxi); //ighost + (i+ighost*(nxi-1))%nxi;//; else IDEFIX_ERROR( "Laplacian:: Shearing box boundary condition should only be used in IDIR" ); + //idfx::cout<=ZERO_F) { // Compute Fl - dqm = localVar(k,jom1,iref) - localVar(k,jom2,iref); - dqp = localVar(k,jo,iref) - localVar(k,jom1,iref); + dqm = localVar(k,jom1,i) - localVar(k,jom2,i); + dqp = localVar(k,jo,i) - localVar(k,jom1,i); dq = dqm+dqp; //(dqp*dqm > ZERO_F ? TWO_F*dqp*dqm/(dqp + dqm) : ZERO_F); - Fl = localVar(k,jom1,iref) + 0.5*dq*(1.0-eps); + Fl = localVar(k,jom1,i) + 0.5*dq*(1.0-eps); //Compute Fr dqm=dqp; - dqp = localVar(k,jop1,iref) - localVar(k,jo,iref); + dqp = localVar(k,jop1,i) - localVar(k,jo,i); dq = dqm+dqp; //(dqp*dqm > ZERO_F ? TWO_F*dqp*dqm/(dqp + dqm) : ZERO_F); - Fr = localVar(k,jo,iref) + 0.5*dq*(1.0-eps); + Fr = localVar(k,jo,i) + 0.5*dq*(1.0-eps); } else { //Compute Fl - dqm = localVar(k,jo,iref) - localVar(k,jom1,iref); - dqp = localVar(k,jop1,iref) - localVar(k,jo,iref); + dqm = localVar(k,jo,i) - localVar(k,jom1,i); + dqp = localVar(k,jop1,i) - localVar(k,jo,i); dq = dqm+dqp; //(dqp*dqm > ZERO_F ? TWO_F*dqp*dqm/(dqp + dqm) : ZERO_F); - Fl = localVar(k,jo,iref) - 0.5*dq*(1.0+eps); + Fl = localVar(k,jo,i) - 0.5*dq*(1.0+eps); // Compute Fr dqm=dqp; - dqp = localVar(k,jop2,iref) - localVar(k,jop1,iref); + dqp = localVar(k,jop2,i) - localVar(k,jop1,i); dq = dqm+dqp; //(dqp*dqm > ZERO_F ? TWO_F*dqp*dqm/(dqp + dqm) : ZERO_F); - Fr = localVar(k,jop1,iref) - 0.5*dq*(1.0+eps); + Fr = localVar(k,jop1,i) - 0.5*dq*(1.0+eps); } - localVar(k,j,i) = localVar(k,jo,iref) - eps*(Fr - Fl); + scrh(k,j,iref) = localVar(k,jo,i) - eps*(Fr - Fl); }); + + + idefix_for("BoundaryShearingBoxCopy", kbeg, kend, jbeg, jend, ibeg, iend, + KOKKOS_LAMBDA ( int k, int j, int i) { + int iref = i - side*(ighost +nxi); //ighost + (i+ighost*(nxi-1))%nxi;//; + localVar(k,j,i) = scrh(k,j,iref); + }); break; } case userdef: { diff --git a/src/gravity/laplacian.hpp b/src/gravity/laplacian.hpp index b70b79a80..b569311eb 100644 --- a/src/gravity/laplacian.hpp +++ b/src/gravity/laplacian.hpp @@ -90,6 +90,8 @@ class Laplacian { bool isTwoPi{false}; bool havePreconditioner{false}; // Use of preconditionner (or not) + + IdefixArray3D sBArray; ///< Array use by shearingbox boundary conditions DataBlock *data; From 716e565b64428ee7dae210ece90788f418aa0004 Mon Sep 17 00:00:00 2001 From: marc Date: Mon, 18 Nov 2024 10:12:20 +0100 Subject: [PATCH 04/28] working as intended --- src/gravity/laplacian.cpp | 6 +++--- src/gravity/laplacian.hpp | 3 +-- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/src/gravity/laplacian.cpp b/src/gravity/laplacian.cpp index dc66a8564..8183595ed 100644 --- a/src/gravity/laplacian.cpp +++ b/src/gravity/laplacian.cpp @@ -564,11 +564,11 @@ void Laplacian::EnforceBoundary(int dir, BoundarySide side, LaplacianBoundaryTyp idefix_for("BoundaryShearingBox", kbeg, kend, jbeg, jend, ibeg, iend, KOKKOS_LAMBDA (int k, int j, int i) { - int iref, jref, kref; + int iref; // This hack takes care of cases where we have more ghost zones than active zones if(dir==IDIR) - iref = i - side*(ighost +nxi); //ighost + (i+ighost*(nxi-1))%nxi;//; + iref = i - side*(ighost +nxi); //ighost + (i+ighost*(nxi-1))%nxi;//; else IDEFIX_ERROR( "Laplacian:: Shearing box boundary condition should only be used in IDIR" @@ -633,7 +633,7 @@ void Laplacian::EnforceBoundary(int dir, BoundarySide side, LaplacianBoundaryTyp idefix_for("BoundaryShearingBoxCopy", kbeg, kend, jbeg, jend, ibeg, iend, KOKKOS_LAMBDA ( int k, int j, int i) { - int iref = i - side*(ighost +nxi); //ighost + (i+ighost*(nxi-1))%nxi;//; + int iref = i - side*(ighost +nxi); //ighost + (i+ighost*(nxi-1))%nxi;//; localVar(k,j,i) = scrh(k,j,iref); }); break; diff --git a/src/gravity/laplacian.hpp b/src/gravity/laplacian.hpp index b569311eb..743d11a42 100644 --- a/src/gravity/laplacian.hpp +++ b/src/gravity/laplacian.hpp @@ -90,9 +90,8 @@ class Laplacian { bool isTwoPi{false}; bool havePreconditioner{false}; // Use of preconditionner (or not) - - IdefixArray3D sBArray; ///< Array use by shearingbox boundary conditions + IdefixArray3D sBArray; ///< Array use by shearingbox boundary conditions DataBlock *data; From e34431c55860abe2423d141cc4aae53885d96c00 Mon Sep 17 00:00:00 2001 From: marc Date: Mon, 18 Nov 2024 14:04:18 +0100 Subject: [PATCH 05/28] No idefix err in idefix for --- src/gravity/laplacian.cpp | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/src/gravity/laplacian.cpp b/src/gravity/laplacian.cpp index 8183595ed..7187c9c8a 100644 --- a/src/gravity/laplacian.cpp +++ b/src/gravity/laplacian.cpp @@ -534,6 +534,11 @@ void Laplacian::EnforceBoundary(int dir, BoundarySide side, LaplacianBoundaryTyp break; } case shearingbox: { + if(dir!=IDIR) + IDEFIX_ERROR( + "Laplacian:: Shearing box boundary condition should only be used in IDIR" + ); + IdefixArray3D scrh = sBArray; const real S = data->hydro->sbS; @@ -564,15 +569,8 @@ void Laplacian::EnforceBoundary(int dir, BoundarySide side, LaplacianBoundaryTyp idefix_for("BoundaryShearingBox", kbeg, kend, jbeg, jend, ibeg, iend, KOKKOS_LAMBDA (int k, int j, int i) { - int iref; - // This hack takes care of cases where we have more ghost zones than active zones + int iref = i - side*(ighost +nxi); //ighost + (i+ighost*(nxi-1))%nxi;//; - if(dir==IDIR) - iref = i - side*(ighost +nxi); //ighost + (i+ighost*(nxi-1))%nxi;//; - else - IDEFIX_ERROR( - "Laplacian:: Shearing box boundary condition should only be used in IDIR" - ); //idfx::cout< Date: Fri, 21 Feb 2025 11:25:43 +0100 Subject: [PATCH 06/28] Fixed IDIR index shift --- src/gravity/laplacian.cpp | 49 +++++--- src/kokkos | 2 +- test/HD/ShearingBox/idefix-fargo.ini | 2 +- test/HD/ShearingBox/idefix.ini | 2 +- test/HD/ShearingBox/python/testidefix.py | 4 +- test/SelfGravity/ShearingSheet/analysis.cpp | 17 +++ test/SelfGravity/ShearingSheet/analysis.hpp | 20 ++++ .../SelfGravity/ShearingSheet/definitions.hpp | 6 + test/SelfGravity/ShearingSheet/idefix.ini | 47 ++++++++ .../ShearingSheet/python/testidefix.py | 93 ++++++++++++++ test/SelfGravity/ShearingSheet/setup.cpp | 113 ++++++++++++++++++ test/SelfGravity/ShearingSheet/testme.py | 40 +++++++ 12 files changed, 374 insertions(+), 21 deletions(-) create mode 100644 test/SelfGravity/ShearingSheet/analysis.cpp create mode 100644 test/SelfGravity/ShearingSheet/analysis.hpp create mode 100644 test/SelfGravity/ShearingSheet/definitions.hpp create mode 100644 test/SelfGravity/ShearingSheet/idefix.ini create mode 100755 test/SelfGravity/ShearingSheet/python/testidefix.py create mode 100644 test/SelfGravity/ShearingSheet/setup.cpp create mode 100755 test/SelfGravity/ShearingSheet/testme.py diff --git a/src/gravity/laplacian.cpp b/src/gravity/laplacian.cpp index 7187c9c8a..14018bcf0 100644 --- a/src/gravity/laplacian.cpp +++ b/src/gravity/laplacian.cpp @@ -534,6 +534,7 @@ void Laplacian::EnforceBoundary(int dir, BoundarySide side, LaplacianBoundaryTyp break; } case shearingbox: { + idfx::cout<<"This is SG shearing box boundary " <mygrid->np_int[IDIR]; + //const int nxiglob = data->mygrid->np_int[IDIR]; idefix_for("BoundaryShearingBox", kbeg, kend, jbeg, jend, ibeg, iend, KOKKOS_LAMBDA (int k, int j, int i) { - int iref = i - side*(ighost +nxi); //ighost + (i+ighost*(nxi-1))%nxi;//; + int iscrh = i - side*(ighost +nxi); + int iref; + if(data->mygrid->nproc[dir]==1) { + // no MPI domain decomposition: we look at the other side of the datablock + iref = ighost + (i+ighost*(nxi-1))%nxi; + } else { + // we need to copy the data to this side (done by MPI), then shift it + iref = i; + } - //idfx::cout<=ZERO_F) { // Compute Fl - dqm = localVar(k,jom1,i) - localVar(k,jom2,i); - dqp = localVar(k,jo,i) - localVar(k,jom1,i); + dqm = localVar(k,jom1,iref) - localVar(k,jom2,iref); + dqp = localVar(k,jo,iref) - localVar(k,jom1,iref); dq = dqm+dqp; //(dqp*dqm > ZERO_F ? TWO_F*dqp*dqm/(dqp + dqm) : ZERO_F); + //dq = (dqp*dqm > ZERO_F ? TWO_F*dqp*dqm/(dqp + dqm) : ZERO_F); - Fl = localVar(k,jom1,i) + 0.5*dq*(1.0-eps); + Fl = localVar(k,jom1,iref) + 0.5*dq*(1.0-eps); //Compute Fr dqm=dqp; - dqp = localVar(k,jop1,i) - localVar(k,jo,i); + dqp = localVar(k,jop1,iref) - localVar(k,jo,iref); dq = dqm+dqp; //(dqp*dqm > ZERO_F ? TWO_F*dqp*dqm/(dqp + dqm) : ZERO_F); + //dq = (dqp*dqm > ZERO_F ? TWO_F*dqp*dqm/(dqp + dqm) : ZERO_F); - Fr = localVar(k,jo,i) + 0.5*dq*(1.0-eps); + Fr = localVar(k,jo,iref) + 0.5*dq*(1.0-eps); } else { //Compute Fl - dqm = localVar(k,jo,i) - localVar(k,jom1,i); - dqp = localVar(k,jop1,i) - localVar(k,jo,i); + dqm = localVar(k,jo,iref) - localVar(k,jom1,iref); + dqp = localVar(k,jop1,iref) - localVar(k,jo,iref); dq = dqm+dqp; //(dqp*dqm > ZERO_F ? TWO_F*dqp*dqm/(dqp + dqm) : ZERO_F); + //dq = (dqp*dqm > ZERO_F ? TWO_F*dqp*dqm/(dqp + dqm) : ZERO_F); - Fl = localVar(k,jo,i) - 0.5*dq*(1.0+eps); + Fl = localVar(k,jo,iref) - 0.5*dq*(1.0+eps); // Compute Fr dqm=dqp; - dqp = localVar(k,jop2,i) - localVar(k,jop1,i); + dqp = localVar(k,jop2,iref) - localVar(k,jop1,iref); dq = dqm+dqp; //(dqp*dqm > ZERO_F ? TWO_F*dqp*dqm/(dqp + dqm) : ZERO_F); + //dq = (dqp*dqm > ZERO_F ? TWO_F*dqp*dqm/(dqp + dqm) : ZERO_F); - Fr = localVar(k,jop1,i) - 0.5*dq*(1.0+eps); + Fr = localVar(k,jop1,iref) - 0.5*dq*(1.0+eps); } - scrh(k,j,iref) = localVar(k,jo,i) - eps*(Fr - Fl); + //idfx::cout< +#include +#include + +Analysis::Analysis() { + } + +/* **************************************************************** */ + +void Analysis::PerformAnalysis(DataBlock &data) { + idfx::pushRegion("Analysis::PerformAnalysis"); + data.DumpToFile("test"); + idfx::popRegion(); +} diff --git a/test/SelfGravity/ShearingSheet/analysis.hpp b/test/SelfGravity/ShearingSheet/analysis.hpp new file mode 100644 index 000000000..f9e412b34 --- /dev/null +++ b/test/SelfGravity/ShearingSheet/analysis.hpp @@ -0,0 +1,20 @@ +#ifndef ANALYSIS_HPP_ +#define ANALYSIS_HPP_ + +#include "dataBlock.hpp" +#include "dataBlockHost.hpp" +#include "grid.hpp" +#include "idefix.hpp" +#include "input.hpp" +#include "output.hpp" +#include +#include + +class Analysis { +public: + // Constructor from Setup arguments + Analysis(); + void PerformAnalysis(DataBlock &); +}; + +#endif // ANALYSIS_HPP__ diff --git a/test/SelfGravity/ShearingSheet/definitions.hpp b/test/SelfGravity/ShearingSheet/definitions.hpp new file mode 100644 index 000000000..d7e581052 --- /dev/null +++ b/test/SelfGravity/ShearingSheet/definitions.hpp @@ -0,0 +1,6 @@ +#define COMPONENTS 2 +#define DIMENSIONS 2 + +#define GEOMETRY CARTESIAN +//#define DEBUG_BICGSTAB +//#define ISOTHERMAL diff --git a/test/SelfGravity/ShearingSheet/idefix.ini b/test/SelfGravity/ShearingSheet/idefix.ini new file mode 100644 index 000000000..c3e785650 --- /dev/null +++ b/test/SelfGravity/ShearingSheet/idefix.ini @@ -0,0 +1,47 @@ +[Grid] +X1-grid 1 -.5 15 u .5 +X2-grid 1 -.5 15 u .5 +X3-grid 1 -0.5 1 u 0.5 + +[TimeIntegrator] +CFL 0.3 +tstop 0.201 +first_dt 1e-4 +nstages 2 + +[Hydro] +solver hllc +rotation 1.0 +shearingBox -1.5 + +[Gravity] +potential selfgravity +bodyForce userdef +gravCst 0.0318309886 # 1/ pi that way Q-1 = Sigma 0 + +[SelfGravity] +solver PBICGSTAB +skip 1 +targetError 1e-4 +maxIter 10000 +boundary-X1-beg shearingbox +boundary-X1-end shearingbox +boundary-X2-beg periodic +boundary-X2-end periodic +boundary-X3-beg periodic +boundary-X3-end periodic + +[Boundary] +X1-beg shearingbox +X1-end shearingbox +X2-beg periodic +X2-end periodic +X3-beg periodic +X3-end periodic + +[Output] +vtk 0.1 +dmp 2.0 +analysis 0.1 +# log 100 +uservar phiP diff --git a/test/SelfGravity/ShearingSheet/python/testidefix.py b/test/SelfGravity/ShearingSheet/python/testidefix.py new file mode 100755 index 000000000..e3db07c3f --- /dev/null +++ b/test/SelfGravity/ShearingSheet/python/testidefix.py @@ -0,0 +1,93 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed Feb 19 15:42:19 2025 + +@author: vdmba +""" +import os +import sys +import numpy as np +import inifix +sys.path.append(os.getenv("IDEFIX_DIR")) +from pytools.idfx_io import readIdfxFile +from pytools.vtk_io import readVTK +import matplotlib.pyplot as plt +import argparse + +# check that the density the potential minima is sheared +# at the correct rate in the boundaries + + +parser = argparse.ArgumentParser() +parser.add_argument("-noplot", + default=True, + help="disable plotting", + action="store_true") + + +args, unknown=parser.parse_known_args() + + +Vbeg = readVTK("../data.0000.vtk") +Vend = readVTK("../data.0002.vtk") + +Iend=readIdfxFile('../test.2.0.idfx') # Last Idfx +Iref=readIdfxFile('../test.0.0.idfx') # Ref vtk (initial state) + +Ly = Vbeg.yl[-1]-Vbeg.yl[0] +xLeft = Vbeg.xl[0] +xRight = Vbeg.xl[-1] + +deltaT = Vend.t[0]-Vbeg.t[0] + +conf = inifix.load("../idefix.ini") + +S = conf["Hydro"]["shearingBox"] + +RHO = 0 +nghost = 2 +dy = Ly/Vbeg.ny + + +indexMaxInitLeftRHO = np.argmax(Iref.data["Vc"][RHO,0,nghost:-nghost,0]) +indexMaxInitRightRHO = np.argmax(Iref.data["Vc"][RHO,0,nghost:-nghost,-1]) + +posMaxInitLeftRHO = Vbeg.y[indexMaxInitLeftRHO] +posMaxInitRightRHO = Vbeg.y[indexMaxInitRightRHO] + +#breakpoint() +if args.noplot: + fig, ax = plt.subplots() + ax.plot(Vbeg.y,Iend.data["Pot"][RHO,0,nghost:-nghost,0]/Iend.data["Pot"][RHO,0,nghost:-nghost,0].max()) + ax.plot(Vbeg.y,Iend.data["Pot"][RHO,0,nghost:-nghost,-1]/Iend.data["Pot"][RHO,0,nghost:-nghost,-1].max()) + plt.show() + +# potential is not computed at 0th step: using the position of the initial density maximum +posMaxInitLeftPot = posMaxInitLeftRHO +posMaxInitRightPot = posMaxInitRightRHO + +indexMaxEndLeftPot = np.argmin(Iend.data["Pot"][RHO,0,nghost:-nghost,0]) +indexMaxEndRightPot = np.argmin(Iend.data["Pot"][RHO,0,nghost:-nghost,-1]) + +posMaxEndLeftPot = Vbeg.y[indexMaxEndLeftPot] +posMaxEndRightPot = Vbeg.y[indexMaxEndRightPot] + +deltaPosLeftPot = posMaxEndLeftPot - posMaxInitLeftPot +deltaPosRightPot = posMaxEndRightPot - posMaxInitRightPot + +errorLeftPot = abs(deltaPosLeftPot - S*xLeft*deltaT) +errorRightPot = abs(deltaPosRightPot - S*xRight*deltaT) + +errorPot = (errorLeftPot**2 + errorRightPot**2)**.5 + + +err=abs(1-errorPot/dy) +print("Error=",err) + +if(err<0.05): + print("SUCCESS") + sys.exit(0) +else: + print("FAILED") + sys.exit(1) diff --git a/test/SelfGravity/ShearingSheet/setup.cpp b/test/SelfGravity/ShearingSheet/setup.cpp new file mode 100644 index 000000000..816ed957b --- /dev/null +++ b/test/SelfGravity/ShearingSheet/setup.cpp @@ -0,0 +1,113 @@ +#include "idefix.hpp" +#include "setup.hpp" +#include "analysis.hpp" + +static real gammaIdeal; +static real omega; +static real shear; +Analysis *analysis; + + +//#define STRATIFIED + +void BodyForce(DataBlock &data, const real t, IdefixArray4D &force) { + idfx::pushRegion("BodyForce"); + IdefixArray1D x = data.x[IDIR]; + IdefixArray1D z = data.x[KDIR]; + + // GPUS cannot capture static variables + real omegaLocal=omega; + real shearLocal =shear; + + idefix_for("BodyForce", + data.beg[KDIR] , data.end[KDIR], + data.beg[JDIR] , data.end[JDIR], + data.beg[IDIR] , data.end[IDIR], + KOKKOS_LAMBDA (int k, int j, int i) { + EXPAND( + force(IDIR,k,j,i) = -2.0*omegaLocal*shearLocal*x(i);, + force(JDIR,k,j,i) = ZERO_F;, + + #ifdef STRATIFIED + force(KDIR,k,j,i) = - omegaLocal*omegaLocal*z(k); + #else + force(KDIR,k,j,i) = ZERO_F; + #endif + ); + }); + + + idfx::popRegion(); +} + +void AnalysisFunction(DataBlock &data) { + analysis->PerformAnalysis(data); +} + +void ComputeUserVars(DataBlock &data, UserDefVariablesContainer &variables) { + Kokkos::deep_copy(variables["phiP"], data.gravity->phiP); +} +// Initialisation routine. Can be used to allocate +// Arrays or variables which are used later on +Setup::Setup(Input &input, Grid &grid, DataBlock &data, Output &output) { + gammaIdeal=data.hydro->eos->GetGamma(); + + // Get rotation rate along vertical axis + omega=input.Get("Hydro","rotation",0); + shear=input.Get("Hydro","shearingBox",0); + + // Add our userstep to the timeintegrator + data.gravity->EnrollBodyForce(BodyForce); + + analysis = new Analysis(); + output.EnrollAnalysis(&AnalysisFunction); + output.EnrollUserDefVariables(&ComputeUserVars); + +} + +// This routine initialize the flow +// Note that data is on the device. +// One can therefore define locally +// a datahost and sync it, if needed +void Setup::InitFlow(DataBlock &data) { + // Create a host copy + DataBlockHost d(data); + real x,y,z; + + real cs2 = gammaIdeal*omega*omega; + + 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++) { + EXPAND( + x=d.x[IDIR](i);, + y=d.x[JDIR](j);, + z=d.x[KDIR](k); + ); + d.Vc(RHO,k,j,i) = 2 + cos(2.0*M_PI*y) + 0.1 * sin(2.0*M_PI*x); + +#ifdef STRATIFIED + d.Vc(RHO,k,j,i) *=exp(-z*z/(2.0)); +#endif +#ifndef ISOTHERMAL + d.Vc(PRS,k,j,i) = d.Vc(RHO,k,j,i)*cs2/gammaIdeal; +#endif + EXPAND( + d.Vc(VX1,k,j,i) = 1e-5*sin(2.0*M_PI*(y+4*z));, + d.Vc(VX2,k,j,i) = shear*x;, + d.Vc(VX3,k,j,i) = 0.0; + ); + + } + } + } + + // Send it all, if needed + d.SyncToDevice(); +} + + +// Analyse data to produce an output + +void MakeAnalysis(DataBlock & data) { +} diff --git a/test/SelfGravity/ShearingSheet/testme.py b/test/SelfGravity/ShearingSheet/testme.py new file mode 100755 index 000000000..f6d371b88 --- /dev/null +++ b/test/SelfGravity/ShearingSheet/testme.py @@ -0,0 +1,40 @@ +#!/usr/bin/env python3 + +""" + +@author: glesur +""" +import os +import sys +sys.path.append(os.getenv("IDEFIX_DIR")) + +import pytools.idfx_test as tst +tolerance=1e-15 +def testMe(test): + test.configure() + test.compile() + inifiles=["idefix.ini"] + + for ini in inifiles: + test.run(inputFile=ini) + if test.init and not test.mpi: + test.makeReference(filename="dump.0001.dmp") + test.standardTest() + test.nonRegressionTest(filename="dump.0001.dmp",tolerance=tolerance) + + +test=tst.idfxTest() +if not test.dec: + test.dec=['2','1','2'] + +if not test.all: + if(test.check): + test.checkOnly(filename="dump.0001.dmp",tolerance=tolerance) + else: + testMe(test) +else: + test.noplot = True + test.single=False + test.reconstruction=2 + test.mpi=False + testMe(test) From 24c161a3b03cfab2789a293628b206b155867b86 Mon Sep 17 00:00:00 2001 From: marc Date: Tue, 25 Feb 2025 12:57:37 +0100 Subject: [PATCH 07/28] Test suite --- test/SelfGravity/ShearingSheet/idefix.ini | 6 ++-- .../ShearingSheet/python/testidefix.py | 31 ++++++++----------- test/SelfGravity/ShearingSheet/setup.cpp | 2 +- test/SelfGravity/ShearingSheet/testme.py | 5 ++- 4 files changed, 21 insertions(+), 23 deletions(-) diff --git a/test/SelfGravity/ShearingSheet/idefix.ini b/test/SelfGravity/ShearingSheet/idefix.ini index c3e785650..c12370a65 100644 --- a/test/SelfGravity/ShearingSheet/idefix.ini +++ b/test/SelfGravity/ShearingSheet/idefix.ini @@ -1,6 +1,6 @@ [Grid] -X1-grid 1 -.5 15 u .5 -X2-grid 1 -.5 15 u .5 +X1-grid 1 -.5 60 u .5 +X2-grid 1 -.5 60 u .5 X3-grid 1 -0.5 1 u 0.5 [TimeIntegrator] @@ -41,7 +41,7 @@ X3-end periodic [Output] vtk 0.1 -dmp 2.0 +dmp 0.2 analysis 0.1 # log 100 uservar phiP diff --git a/test/SelfGravity/ShearingSheet/python/testidefix.py b/test/SelfGravity/ShearingSheet/python/testidefix.py index e3db07c3f..725907ab0 100755 --- a/test/SelfGravity/ShearingSheet/python/testidefix.py +++ b/test/SelfGravity/ShearingSheet/python/testidefix.py @@ -10,12 +10,11 @@ import numpy as np import inifix sys.path.append(os.getenv("IDEFIX_DIR")) -from pytools.idfx_io import readIdfxFile from pytools.vtk_io import readVTK import matplotlib.pyplot as plt import argparse -# check that the density the potential minima is sheared +# check that the potential minima is sheared # at the correct rate in the boundaries @@ -32,12 +31,9 @@ Vbeg = readVTK("../data.0000.vtk") Vend = readVTK("../data.0002.vtk") -Iend=readIdfxFile('../test.2.0.idfx') # Last Idfx -Iref=readIdfxFile('../test.0.0.idfx') # Ref vtk (initial state) - Ly = Vbeg.yl[-1]-Vbeg.yl[0] -xLeft = Vbeg.xl[0] -xRight = Vbeg.xl[-1] +xLeft = Vbeg.x[0] +xRight = Vbeg.x[-1] deltaT = Vend.t[0]-Vbeg.t[0] @@ -45,30 +41,28 @@ S = conf["Hydro"]["shearingBox"] -RHO = 0 -nghost = 2 dy = Ly/Vbeg.ny -indexMaxInitLeftRHO = np.argmax(Iref.data["Vc"][RHO,0,nghost:-nghost,0]) -indexMaxInitRightRHO = np.argmax(Iref.data["Vc"][RHO,0,nghost:-nghost,-1]) +indexMaxInitLeftRHO = np.argmax(Vbeg.data["RHO"][0,:,0]) +indexMaxInitRightRHO = np.argmax(Vbeg.data["RHO"][-1,:,0]) posMaxInitLeftRHO = Vbeg.y[indexMaxInitLeftRHO] posMaxInitRightRHO = Vbeg.y[indexMaxInitRightRHO] #breakpoint() -if args.noplot: +if not args.noplot: fig, ax = plt.subplots() - ax.plot(Vbeg.y,Iend.data["Pot"][RHO,0,nghost:-nghost,0]/Iend.data["Pot"][RHO,0,nghost:-nghost,0].max()) - ax.plot(Vbeg.y,Iend.data["Pot"][RHO,0,nghost:-nghost,-1]/Iend.data["Pot"][RHO,0,nghost:-nghost,-1].max()) + ax.plot(Vend.y,Vend.data["phiP"][0,:,0]/Vend.data["phiP"][0,:,:].max(),".") + ax.plot(Vend.y,Vend.data["phiP"][-1,:,0]/Vend.data["phiP"][-1,:,0].max(),".") plt.show() # potential is not computed at 0th step: using the position of the initial density maximum posMaxInitLeftPot = posMaxInitLeftRHO posMaxInitRightPot = posMaxInitRightRHO -indexMaxEndLeftPot = np.argmin(Iend.data["Pot"][RHO,0,nghost:-nghost,0]) -indexMaxEndRightPot = np.argmin(Iend.data["Pot"][RHO,0,nghost:-nghost,-1]) +indexMaxEndLeftPot = np.argmin(Vend.data["phiP"][0,:,0]) +indexMaxEndRightPot = np.argmin(Vend.data["phiP"][-1,:,0]) posMaxEndLeftPot = Vbeg.y[indexMaxEndLeftPot] posMaxEndRightPot = Vbeg.y[indexMaxEndRightPot] @@ -82,10 +76,11 @@ errorPot = (errorLeftPot**2 + errorRightPot**2)**.5 -err=abs(1-errorPot/dy) + +err=abs(errorPot/dy) print("Error=",err) -if(err<0.05): +if(err<1): # total cummulated shift is less than on cell print("SUCCESS") sys.exit(0) else: diff --git a/test/SelfGravity/ShearingSheet/setup.cpp b/test/SelfGravity/ShearingSheet/setup.cpp index 816ed957b..8123c4c5a 100644 --- a/test/SelfGravity/ShearingSheet/setup.cpp +++ b/test/SelfGravity/ShearingSheet/setup.cpp @@ -84,7 +84,7 @@ void Setup::InitFlow(DataBlock &data) { y=d.x[JDIR](j);, z=d.x[KDIR](k); ); - d.Vc(RHO,k,j,i) = 2 + cos(2.0*M_PI*y) + 0.1 * sin(2.0*M_PI*x); + d.Vc(RHO,k,j,i) = 2 + cos(2.0*M_PI*y);//+ 0.2 * sin(2.0*M_PI*x); #ifdef STRATIFIED d.Vc(RHO,k,j,i) *=exp(-z*z/(2.0)); diff --git a/test/SelfGravity/ShearingSheet/testme.py b/test/SelfGravity/ShearingSheet/testme.py index f6d371b88..9da97a788 100755 --- a/test/SelfGravity/ShearingSheet/testme.py +++ b/test/SelfGravity/ShearingSheet/testme.py @@ -25,7 +25,7 @@ def testMe(test): test=tst.idfxTest() if not test.dec: - test.dec=['2','1','2'] + test.dec=['3','1'] if not test.all: if(test.check): @@ -38,3 +38,6 @@ def testMe(test): test.reconstruction=2 test.mpi=False testMe(test) + + test.mpi=True + testMe(test) From 985d454de8569c0984421b2c82ab6284654680ba Mon Sep 17 00:00:00 2001 From: marc Date: Tue, 25 Feb 2025 13:03:05 +0100 Subject: [PATCH 08/28] Cleanup --- src/gravity/laplacian.cpp | 110 +++++++----------- .../ShearingSheet/python/testidefix.py | 2 +- 2 files changed, 42 insertions(+), 70 deletions(-) diff --git a/src/gravity/laplacian.cpp b/src/gravity/laplacian.cpp index 14018bcf0..f03222f5a 100644 --- a/src/gravity/laplacian.cpp +++ b/src/gravity/laplacian.cpp @@ -534,7 +534,6 @@ void Laplacian::EnforceBoundary(int dir, BoundarySide side, LaplacianBoundaryTyp break; } case shearingbox: { - idfx::cout<<"This is SG shearing box boundary " < (std::floor(dL/dy+HALF_F)); - //const real eps = dL - m *dy; const real eps = dL/dy - m; - const int sideShift =0;// (side +1)%2; - - //const int nxiglob = data->mygrid->np_int[IDIR]; - idefix_for("BoundaryShearingBox", kbeg, kend, jbeg, jend, ibeg, iend, KOKKOS_LAMBDA (int k, int j, int i) { - int iscrh = i - side*(ighost +nxi); - int iref; - if(data->mygrid->nproc[dir]==1) { - // no MPI domain decomposition: we look at the other side of the datablock - iref = ighost + (i+ighost*(nxi-1))%nxi; - } else { - // we need to copy the data to this side (done by MPI), then shift it - iref = i; - } + const int iscrh = i - side*(ighost +nxi); - //idfx::cout<mygrid->nproc[dir]==1) ? + ighost + (i+ighost*(nxi-1))%nxi : i; - //localVar(k,j,i) = localVar(k,j,iref)+j; // this works - - // il faut prendre en compte les ghost zones d'une façon ou d'une autre - // ce n'est pas le cas ici... - - const int jo = jghost + ((j-m-jghost+sideShift)%nxj+nxj)%nxj; + const int jo = jghost + ((j-m-jghost)%nxj+nxj)%nxj; const int jop1 = jghost + ((jo+1-jghost)%nxj+nxj)%nxj; const int jom1 = jghost + ((jo-1-jghost)%nxj+nxj)%nxj; const int jom2 = jghost + ((jo-2-jghost)%nxj+nxj)%nxj; const int jop2 = jghost + ((jo+2-jghost)%nxj+nxj)%nxj; - //if (j==0) - // idfx::cout<< t <<"\t"<=ZERO_F) { + // Compute Fl + dqm = localVar(k,jom1,iref) - localVar(k,jom2,iref); + dqp = localVar(k,jo,iref) - localVar(k,jom1,iref); + dq = dqm+dqp; + //dq = (dqp*dqm > ZERO_F ? TWO_F*dqp*dqm/(dqp + dqm) : ZERO_F); + + Fl = localVar(k,jom1,iref) + 0.5*dq*(1.0-eps); + //Compute Fr + dqm=dqp; + dqp = localVar(k,jop1,iref) - localVar(k,jo,iref); + dq = dqm+dqp; + //dq = (dqp*dqm > ZERO_F ? TWO_F*dqp*dqm/(dqp + dqm) : ZERO_F); + + Fr = localVar(k,jo,iref) + 0.5*dq*(1.0-eps); } else { - localVar(k,j,i) = ((1-eps) * localVar(k,jom1,iref) - + eps*localVar(k,jo,iref)); - }*/ - - // the follwing, like for the fluxes in boundary.hpp is a second- - // order recontruction in eps - if(eps>=ZERO_F) { - // Compute Fl - dqm = localVar(k,jom1,iref) - localVar(k,jom2,iref); - dqp = localVar(k,jo,iref) - localVar(k,jom1,iref); - dq = dqm+dqp; //(dqp*dqm > ZERO_F ? TWO_F*dqp*dqm/(dqp + dqm) : ZERO_F); - //dq = (dqp*dqm > ZERO_F ? TWO_F*dqp*dqm/(dqp + dqm) : ZERO_F); - - Fl = localVar(k,jom1,iref) + 0.5*dq*(1.0-eps); - //Compute Fr - dqm=dqp; - dqp = localVar(k,jop1,iref) - localVar(k,jo,iref); - dq = dqm+dqp; //(dqp*dqm > ZERO_F ? TWO_F*dqp*dqm/(dqp + dqm) : ZERO_F); - //dq = (dqp*dqm > ZERO_F ? TWO_F*dqp*dqm/(dqp + dqm) : ZERO_F); - - Fr = localVar(k,jo,iref) + 0.5*dq*(1.0-eps); - } else { - //Compute Fl - dqm = localVar(k,jo,iref) - localVar(k,jom1,iref); - dqp = localVar(k,jop1,iref) - localVar(k,jo,iref); - dq = dqm+dqp; //(dqp*dqm > ZERO_F ? TWO_F*dqp*dqm/(dqp + dqm) : ZERO_F); - //dq = (dqp*dqm > ZERO_F ? TWO_F*dqp*dqm/(dqp + dqm) : ZERO_F); - - Fl = localVar(k,jo,iref) - 0.5*dq*(1.0+eps); - // Compute Fr - dqm=dqp; - dqp = localVar(k,jop2,iref) - localVar(k,jop1,iref); - dq = dqm+dqp; //(dqp*dqm > ZERO_F ? TWO_F*dqp*dqm/(dqp + dqm) : ZERO_F); - //dq = (dqp*dqm > ZERO_F ? TWO_F*dqp*dqm/(dqp + dqm) : ZERO_F); - - Fr = localVar(k,jop1,iref) - 0.5*dq*(1.0+eps); - } - //idfx::cout< ZERO_F ? TWO_F*dqp*dqm/(dqp + dqm) : ZERO_F); + + Fl = localVar(k,jo,iref) - 0.5*dq*(1.0+eps); + // Compute Fr + dqm=dqp; + dqp = localVar(k,jop2,iref) - localVar(k,jop1,iref); + dq = dqm+dqp; + //dq = (dqp*dqm > ZERO_F ? TWO_F*dqp*dqm/(dqp + dqm) : ZERO_F); + + Fr = localVar(k,jop1,iref) - 0.5*dq*(1.0+eps); + } + scrh(k,j,iscrh) = localVar(k,jo,iref) - eps*(Fr - Fl); }); idefix_for("BoundaryShearingBoxCopy", kbeg, kend, jbeg, jend, ibeg, iend, KOKKOS_LAMBDA ( int k, int j, int i) { - int iscrh = i - side*(ighost +nxi); //ighost + (i+ighost*(nxi-1))%nxi;//; + const int iscrh = i - side*(ighost +nxi); localVar(k,j,i) = scrh(k,j,iscrh); }); break; diff --git a/test/SelfGravity/ShearingSheet/python/testidefix.py b/test/SelfGravity/ShearingSheet/python/testidefix.py index 725907ab0..b1b65fab8 100755 --- a/test/SelfGravity/ShearingSheet/python/testidefix.py +++ b/test/SelfGravity/ShearingSheet/python/testidefix.py @@ -3,7 +3,7 @@ """ Created on Wed Feb 19 15:42:19 2025 -@author: vdmba +@author: vdbma """ import os import sys From d5fdcdf552a67c36806279e89f77ef35cc422623 Mon Sep 17 00:00:00 2001 From: marc Date: Tue, 25 Feb 2025 13:32:21 +0100 Subject: [PATCH 09/28] Doc --- doc/source/modules/selfGravity.rst | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/doc/source/modules/selfGravity.rst b/doc/source/modules/selfGravity.rst index 6254caa8d..ed20fbdfa 100644 --- a/doc/source/modules/selfGravity.rst +++ b/doc/source/modules/selfGravity.rst @@ -87,15 +87,18 @@ The boundary conditions can be following | origin | | In spherical coordinates, artificially extends the grid used to compute the potential close to R=0. | | | | should only be used in X1-beg direction. | +-----------------------+------------------------------------------------------------------------------------------------------------------+ +| shearingbox | Shearingbox boundary conditions. | ++-----------------------+------------------------------------------------------------------------------------------------------------------+ | userdef | |User-defined boundary conditions. The boundary condition function should be enrolled in the setup constructor | | | | (see :ref:`userdefBoundaries`). | +-----------------------+------------------------------------------------------------------------------------------------------------------+ .. note:: The method in fully periodic setups requires the removal of the mean gas density - before solving Poisson equation. This is done automatically if all of the self-gravity boundaries are set to ``periodic``. - Hence, make sure to specify all self-gravity boundary conditions as periodic for such setups, otherwise the solver will + before solving Poisson equation. This is done automatically if all of the self-gravity boundaries are set to ``periodic`` or ``shearingbox``. + Hence, make sure to specify all self-gravity boundary conditions as periodic (or shearingbox) for such setups, otherwise the solver will fail to converge. + Note that in these cases, if the initial density is uniform, the self-gravity solver will fail as this correspond to a vaccum initial state. The selfGravity module in *Idefix* is fully parallelised. This means that one can have a MPI domain decomposition in any spatial direction either on CPU or GPU. From 839bc0305852c17fce2f9b390482ee1ae576621f Mon Sep 17 00:00:00 2001 From: marc Date: Tue, 25 Feb 2025 15:02:29 +0100 Subject: [PATCH 10/28] Removed old idfx analysis --- test/SelfGravity/ShearingSheet/analysis.cpp | 17 ----------------- test/SelfGravity/ShearingSheet/analysis.hpp | 20 -------------------- test/SelfGravity/ShearingSheet/idefix.ini | 7 +++---- test/SelfGravity/ShearingSheet/setup.cpp | 14 -------------- 4 files changed, 3 insertions(+), 55 deletions(-) delete mode 100644 test/SelfGravity/ShearingSheet/analysis.cpp delete mode 100644 test/SelfGravity/ShearingSheet/analysis.hpp diff --git a/test/SelfGravity/ShearingSheet/analysis.cpp b/test/SelfGravity/ShearingSheet/analysis.cpp deleted file mode 100644 index 743373ce1..000000000 --- a/test/SelfGravity/ShearingSheet/analysis.cpp +++ /dev/null @@ -1,17 +0,0 @@ -#include "analysis.hpp" -#include "fluid.hpp" -#include "idefix.hpp" -#include -#include -#include - -Analysis::Analysis() { - } - -/* **************************************************************** */ - -void Analysis::PerformAnalysis(DataBlock &data) { - idfx::pushRegion("Analysis::PerformAnalysis"); - data.DumpToFile("test"); - idfx::popRegion(); -} diff --git a/test/SelfGravity/ShearingSheet/analysis.hpp b/test/SelfGravity/ShearingSheet/analysis.hpp deleted file mode 100644 index f9e412b34..000000000 --- a/test/SelfGravity/ShearingSheet/analysis.hpp +++ /dev/null @@ -1,20 +0,0 @@ -#ifndef ANALYSIS_HPP_ -#define ANALYSIS_HPP_ - -#include "dataBlock.hpp" -#include "dataBlockHost.hpp" -#include "grid.hpp" -#include "idefix.hpp" -#include "input.hpp" -#include "output.hpp" -#include -#include - -class Analysis { -public: - // Constructor from Setup arguments - Analysis(); - void PerformAnalysis(DataBlock &); -}; - -#endif // ANALYSIS_HPP__ diff --git a/test/SelfGravity/ShearingSheet/idefix.ini b/test/SelfGravity/ShearingSheet/idefix.ini index c12370a65..3508d6f91 100644 --- a/test/SelfGravity/ShearingSheet/idefix.ini +++ b/test/SelfGravity/ShearingSheet/idefix.ini @@ -40,8 +40,7 @@ X3-beg periodic X3-end periodic [Output] -vtk 0.1 -dmp 0.2 -analysis 0.1 +vtk 0.1 +dmp 0.2 # log 100 -uservar phiP +uservar phiP diff --git a/test/SelfGravity/ShearingSheet/setup.cpp b/test/SelfGravity/ShearingSheet/setup.cpp index 8123c4c5a..3ed35acec 100644 --- a/test/SelfGravity/ShearingSheet/setup.cpp +++ b/test/SelfGravity/ShearingSheet/setup.cpp @@ -1,12 +1,9 @@ #include "idefix.hpp" #include "setup.hpp" -#include "analysis.hpp" static real gammaIdeal; static real omega; static real shear; -Analysis *analysis; - //#define STRATIFIED @@ -40,9 +37,6 @@ void BodyForce(DataBlock &data, const real t, IdefixArray4D &force) { idfx::popRegion(); } -void AnalysisFunction(DataBlock &data) { - analysis->PerformAnalysis(data); -} void ComputeUserVars(DataBlock &data, UserDefVariablesContainer &variables) { Kokkos::deep_copy(variables["phiP"], data.gravity->phiP); @@ -59,8 +53,6 @@ Setup::Setup(Input &input, Grid &grid, DataBlock &data, Output &output) { // Add our userstep to the timeintegrator data.gravity->EnrollBodyForce(BodyForce); - analysis = new Analysis(); - output.EnrollAnalysis(&AnalysisFunction); output.EnrollUserDefVariables(&ComputeUserVars); } @@ -105,9 +97,3 @@ void Setup::InitFlow(DataBlock &data) { // Send it all, if needed d.SyncToDevice(); } - - -// Analyse data to produce an output - -void MakeAnalysis(DataBlock & data) { -} From f5d585d25ecce39a479e5fc402709a01694c8a9c Mon Sep 17 00:00:00 2001 From: marc Date: Tue, 25 Feb 2025 15:08:33 +0100 Subject: [PATCH 11/28] lambda capture --- src/gravity/laplacian.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/gravity/laplacian.cpp b/src/gravity/laplacian.cpp index f03222f5a..c049db524 100644 --- a/src/gravity/laplacian.cpp +++ b/src/gravity/laplacian.cpp @@ -562,13 +562,15 @@ void Laplacian::EnforceBoundary(int dir, BoundarySide side, LaplacianBoundaryTyp const real eps = dL/dy - m; + const int nprocsIDIR = data->mygrid->nproc[dir]; + idefix_for("BoundaryShearingBox", kbeg, kend, jbeg, jend, ibeg, iend, KOKKOS_LAMBDA (int k, int j, int i) { const int iscrh = i - side*(ighost +nxi); // no MPI domain decomposition: we look at the other side of the datablock // MPI: we need to copy the data to this side (done by MPI), then shift it - const int iref = (data->mygrid->nproc[dir]==1) ? + const int iref = (nprocsIDIR==1) ? ighost + (i+ighost*(nxi-1))%nxi : i; const int jo = jghost + ((j-m-jghost)%nxj+nxj)%nxj; From e64cabca0f648e0bf3118db5e75b604b8724ac52 Mon Sep 17 00:00:00 2001 From: marc Date: Mon, 17 Mar 2025 15:13:38 +0100 Subject: [PATCH 12/28] Revert kokkos version change --- src/kokkos | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/kokkos b/src/kokkos index 6ecdf605e..09e775bfc 160000 --- a/src/kokkos +++ b/src/kokkos @@ -1 +1 @@ -Subproject commit 6ecdf605e0f7639adec599d25cf0e206d7b8f9f5 +Subproject commit 09e775bfc585840bb9ab1156cbd8d7d1c9e0cc6d From 416c989ef1d481b390eaef7219705c425c65adc2 Mon Sep 17 00:00:00 2001 From: marc Date: Mon, 17 Mar 2025 15:15:31 +0100 Subject: [PATCH 13/28] Changed reference version --- test/HD/ShearingBox/idefix-fargo.ini | 2 +- test/HD/ShearingBox/idefix.ini | 2 +- test/HD/ShearingBox/python/testidefix.py | 4 ++-- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/test/HD/ShearingBox/idefix-fargo.ini b/test/HD/ShearingBox/idefix-fargo.ini index be1e56f98..168c9ee2e 100644 --- a/test/HD/ShearingBox/idefix-fargo.ini +++ b/test/HD/ShearingBox/idefix-fargo.ini @@ -29,6 +29,6 @@ X3-beg periodic X3-end periodic [Output] -vtk 0.1 +# vtk 0.1 dmp 2.0 analysis 0.1 diff --git a/test/HD/ShearingBox/idefix.ini b/test/HD/ShearingBox/idefix.ini index 3b02100a7..a70674b8f 100644 --- a/test/HD/ShearingBox/idefix.ini +++ b/test/HD/ShearingBox/idefix.ini @@ -26,6 +26,6 @@ X3-beg periodic X3-end periodic [Output] -vtk 0.1 +# vtk 0.1 dmp 2.0 analysis 0.1 diff --git a/test/HD/ShearingBox/python/testidefix.py b/test/HD/ShearingBox/python/testidefix.py index 3584c951b..1a0c2509a 100755 --- a/test/HD/ShearingBox/python/testidefix.py +++ b/test/HD/ShearingBox/python/testidefix.py @@ -1,9 +1,9 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- """ -Created on Wed Feb 19 15:42:19 2025 +Created on Mon Jun 21 15:42:19 2021 -@author: vdmba +@author: lesurg """ import sys import numpy as np From c516a63d684f03d0c57e5569a2bba7b4cc89e84a Mon Sep 17 00:00:00 2001 From: marc Date: Tue, 8 Apr 2025 19:11:14 +0200 Subject: [PATCH 14/28] Passing with no SG, failing with SG --- test/SelfGravity/ShearingWave/definitions.hpp | 6 + test/SelfGravity/ShearingWave/idefix.ini | 47 ++++++++ .../ShearingWave/python/testidefix.py | 106 +++++++++++++++++ test/SelfGravity/ShearingWave/setup.cpp | 111 ++++++++++++++++++ test/SelfGravity/ShearingWave/testme.py | 43 +++++++ 5 files changed, 313 insertions(+) create mode 100644 test/SelfGravity/ShearingWave/definitions.hpp create mode 100644 test/SelfGravity/ShearingWave/idefix.ini create mode 100755 test/SelfGravity/ShearingWave/python/testidefix.py create mode 100644 test/SelfGravity/ShearingWave/setup.cpp create mode 100755 test/SelfGravity/ShearingWave/testme.py diff --git a/test/SelfGravity/ShearingWave/definitions.hpp b/test/SelfGravity/ShearingWave/definitions.hpp new file mode 100644 index 000000000..548c27fa3 --- /dev/null +++ b/test/SelfGravity/ShearingWave/definitions.hpp @@ -0,0 +1,6 @@ +#define COMPONENTS 2 +#define DIMENSIONS 2 + +#define GEOMETRY CARTESIAN +//#define DEBUG_BICGSTAB +#define ISOTHERMAL diff --git a/test/SelfGravity/ShearingWave/idefix.ini b/test/SelfGravity/ShearingWave/idefix.ini new file mode 100644 index 000000000..ad3c1dfa8 --- /dev/null +++ b/test/SelfGravity/ShearingWave/idefix.ini @@ -0,0 +1,47 @@ +[Grid] +X1-grid 1 -.5 128 u .5 +X2-grid 1 -.5 128 u .5 +X3-grid 1 -0.5 1 u 0.5 + +[TimeIntegrator] +CFL 0.3 +tstop 10 +first_dt 1e-4 +nstages 2 + +[Hydro] +solver hllc +rotation 1. +shearingBox -1.5 +csiso constant 0.07853981633974483 + +[Gravity] +potential selfgravity +bodyForce userdef +gravCst 1 + +[SelfGravity] +solver PBICGSTAB +skip 1 +targetError 1e-6 +maxIter 10000 +boundary-X1-beg shearingbox +boundary-X1-end shearingbox +boundary-X2-beg periodic +boundary-X2-end periodic +boundary-X3-beg periodic +boundary-X3-end periodic + +[Boundary] +X1-beg shearingbox +X1-end shearingbox +X2-beg periodic +X2-end periodic +X3-beg periodic +X3-end periodic + +[Output] +vtk 0.05 +dmp 0.2 +log 100 +# uservar phiP diff --git a/test/SelfGravity/ShearingWave/python/testidefix.py b/test/SelfGravity/ShearingWave/python/testidefix.py new file mode 100755 index 000000000..8aff88c06 --- /dev/null +++ b/test/SelfGravity/ShearingWave/python/testidefix.py @@ -0,0 +1,106 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed Feb 19 15:42:19 2025 + +@author: vdbma +""" +import os +import sys +import numpy as np +import inifix +sys.path.append(os.getenv("IDEFIX_DIR")) +from pytools.vtk_io import readVTK +import argparse +import scipy.integrate as si + +# check that the potential minima is sheared +# at the correct rate in the boundaries + + +parser = argparse.ArgumentParser() +parser.add_argument("-noplot", + default=True, + help="disable plotting", + action="store_true") + + +args, unknown=parser.parse_known_args() + +conf = inifix.load('../idefix.ini') + +# unit dependant stuff +Omega = conf["Hydro"]["rotation"] +G = conf["Gravity"]["gravCst"] +# user parameters +q = -conf["Hydro"]["shearingBox"] +isSG = "potential" in conf["Gravity"] + + +kappa = 2*Omega*Omega*(1-q) +sigma0 = 0.0005 +Sigma0 = 1/40 +kx0 = -4*np.pi +ky = 2*np.pi +Q = 1 + + +xi0 = (2-q)*Omega/Sigma0 +xi1 = -sigma0*xi0 +tmax = conf["TimeIntegrator"]["tstop"] + +cs = np.pi*G*Sigma0*Q/kappa + + +def f2s(t,y): + kx = kx0 + q*Omega*ky*t + k = np.sqrt(kx*kx+ky*ky) + + b = -2*q*Omega*kx*ky/k/k + c = 2*(2-q)*Omega*Omega+k**2*cs**2 - 2*np.pi *G *Sigma0 *k*isSG - 2*q*(2-q)*Omega*Omega*ky*ky/k/k + d = 2*Omega*(1-q*ky*ky/k/k)*Sigma0*xi1 + + sp,s = y + f1 = -b*sp-c*s-d + f2 = sp + return np.array([f1,f2]) + + +y0 = [0,sigma0] + + +amplitude = [] +amp_exact = [] +time = [] + +sol = si.solve_ivp(f2s,y0=y0,t_span = [0,tmax],t_eval=np.linspace(0,tmax,2048)) + + +for i in range(200): + V = readVTK(f"../data.{i:04}.vtk",geometry="cartesian") + t=V.t[0] + xx,yy = np.meshgrid(V.x,V.y,indexing="ij") + time.append(t) + + mode = np.exp(1j*((kx0+q*Omega*t*ky)*xx + ky*yy)) + amplitude.append(np.sum(V.data["RHO"][:,:,0]*mode.conjugate())/V.nx/V.ny*2) # 2 parce que int cos**2 = 1/2 + i_t = np.argmin(np.abs(t-sol.t)) + amp_exact.append(sol.y[1][i_t]) + + +time = np.array(time) +amplitude=np.array(amplitude).real +amp_exact=np.array(amp_exact) + + + + +err=np.abs((amp_exact*Sigma0-amplitude)/Sigma0).max() +print("Error=",err) + +if(err<1e-3): + print("SUCCESS") + sys.exit(0) +else: + print("FAILED") + sys.exit(1) diff --git a/test/SelfGravity/ShearingWave/setup.cpp b/test/SelfGravity/ShearingWave/setup.cpp new file mode 100644 index 000000000..7399dd491 --- /dev/null +++ b/test/SelfGravity/ShearingWave/setup.cpp @@ -0,0 +1,111 @@ +#include "idefix.hpp" +#include "setup.hpp" + +#ifndef ISOTHERMAL +static real gammaIdeal; +#endif +static real omega; +static real shear; +static real Gnewt; + +//#define STRATIFIED + +void BodyForce(DataBlock &data, const real t, IdefixArray4D &force) { + idfx::pushRegion("BodyForce"); + IdefixArray1D x = data.x[IDIR]; + IdefixArray1D z = data.x[KDIR]; + + // GPUS cannot capture static variables + real omegaLocal=omega; + real shearLocal =shear; + + idefix_for("BodyForce", + data.beg[KDIR] , data.end[KDIR], + data.beg[JDIR] , data.end[JDIR], + data.beg[IDIR] , data.end[IDIR], + KOKKOS_LAMBDA (int k, int j, int i) { + EXPAND( + force(IDIR,k,j,i) = -2.0*omegaLocal*shearLocal*x(i);, + force(JDIR,k,j,i) = ZERO_F;, + + #ifdef STRATIFIED + force(KDIR,k,j,i) = - omegaLocal*omegaLocal*z(k); + #else + force(KDIR,k,j,i) = ZERO_F; + #endif + ); + }); + + + idfx::popRegion(); +} + + +void ComputeUserVars(DataBlock &data, UserDefVariablesContainer &variables) { + Kokkos::deep_copy(variables["phiP"], data.gravity->phiP); +} +// Initialisation routine. Can be used to allocate +// Arrays or variables which are used later on +Setup::Setup(Input &input, Grid &grid, DataBlock &data, Output &output) { + #ifndef ISOTHERMAL + gammaIdeal=data.hydro->eos->GetGamma(); + #endif + + // Get rotation rate along vertical axis + omega=input.Get("Hydro","rotation",0); + shear=input.Get("Hydro","shearingBox",0); + Gnewt=input.Get("Gravity","gravCst",0); + + // Add our userstep to the timeintegrator + data.gravity->EnrollBodyForce(BodyForce); + + //output.EnrollUserDefVariables(&ComputeUserVars); + +} + +// This routine initialize the flow +// Note that data is on the device. +// One can therefore define locally +// a datahost and sync it, if needed +void Setup::InitFlow(DataBlock &data) { + // Create a host copy + DataBlockHost d(data); + real x,y,z = 0; + + // setup from Paardekooper 2012 + real kx = -4.*M_PI; + real ky = 2.*M_PI; + real sigma_init = 0.0005; + real Sigma0 = 1./40.; + real Q = 1.; + + real kappa = omega; + + real cs = M_PI*Sigma0*Q/kappa*Gnewt; // G = 1 + std::cout<< "cs = " << cs << std::endl; + + 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++) { + EXPAND( + x=d.x[IDIR](i);, + y=d.x[JDIR](j);, + z=d.x[KDIR](k); + ); + d.Vc(RHO,k,j,i) = Sigma0*(1.+ sigma_init*cos(kx*x+ky*y)); +#ifndef ISOTHERMAL + d.Vc(PRS,k,j,i) = d.Vc(RHO,k,j,i)*cs*cs; +#endif + EXPAND( + d.Vc(VX1,k,j,i) = 0.0;, //sigma_init/Sigma0*cos(kx*x+ky*y);, + d.Vc(VX2,k,j,i) = shear*x;,// + sigma_init/Sigma0*cos(kx*x+ky*y);, + d.Vc(VX3,k,j,i) = 0.0; + ); + + } + } + } + + // Send it all, if needed + d.SyncToDevice(); +} diff --git a/test/SelfGravity/ShearingWave/testme.py b/test/SelfGravity/ShearingWave/testme.py new file mode 100755 index 000000000..c4128efc8 --- /dev/null +++ b/test/SelfGravity/ShearingWave/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 +tolerance=1e-15 +def testMe(test): + test.configure() + test.compile() + inifiles=["idefix.ini"] + + for ini in inifiles: + test.run(inputFile=ini) + if test.init and not test.mpi: + test.makeReference(filename="dump.0001.dmp") + test.standardTest() + test.nonRegressionTest(filename="dump.0001.dmp",tolerance=tolerance) + + +test=tst.idfxTest() +if not test.dec: + test.dec=['8','1'] + +if not test.all: + if(test.check): + test.checkOnly(filename="dump.0001.dmp",tolerance=tolerance) + else: + testMe(test) +else: + test.noplot = True + test.single=False + test.reconstruction=2 + test.mpi=False + testMe(test) + + test.mpi=True + testMe(test) From 9a0e79f780f7b43a04daab8589547771c2710cb5 Mon Sep 17 00:00:00 2001 From: marc Date: Mon, 5 May 2025 17:56:12 +0200 Subject: [PATCH 15/28] Cleanup --- test/SelfGravity/ShearingWave/setup.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/SelfGravity/ShearingWave/setup.cpp b/test/SelfGravity/ShearingWave/setup.cpp index 7399dd491..3f88a4651 100644 --- a/test/SelfGravity/ShearingWave/setup.cpp +++ b/test/SelfGravity/ShearingWave/setup.cpp @@ -59,8 +59,6 @@ Setup::Setup(Input &input, Grid &grid, DataBlock &data, Output &output) { // Add our userstep to the timeintegrator data.gravity->EnrollBodyForce(BodyForce); - //output.EnrollUserDefVariables(&ComputeUserVars); - } // This routine initialize the flow @@ -73,6 +71,7 @@ void Setup::InitFlow(DataBlock &data) { real x,y,z = 0; // setup from Paardekooper 2012 + // but not quite because we are in 2D here not rho ~ delta(z) real kx = -4.*M_PI; real ky = 2.*M_PI; real sigma_init = 0.0005; @@ -82,6 +81,7 @@ void Setup::InitFlow(DataBlock &data) { real kappa = omega; real cs = M_PI*Sigma0*Q/kappa*Gnewt; // G = 1 + cs = 0.07853981633974483; std::cout<< "cs = " << cs << std::endl; for(int k = 0; k < d.np_tot[KDIR] ; k++) { @@ -97,8 +97,8 @@ void Setup::InitFlow(DataBlock &data) { d.Vc(PRS,k,j,i) = d.Vc(RHO,k,j,i)*cs*cs; #endif EXPAND( - d.Vc(VX1,k,j,i) = 0.0;, //sigma_init/Sigma0*cos(kx*x+ky*y);, - d.Vc(VX2,k,j,i) = shear*x;,// + sigma_init/Sigma0*cos(kx*x+ky*y);, + d.Vc(VX1,k,j,i) = 0.0;, + d.Vc(VX2,k,j,i) = shear*x;, d.Vc(VX3,k,j,i) = 0.0; ); From 7ef44e811d1b96316fc7a2c6ecc693968357fde0 Mon Sep 17 00:00:00 2001 From: marc Date: Mon, 5 May 2025 17:56:42 +0200 Subject: [PATCH 16/28] Correct 2D potential --- .../ShearingWave/python/testidefix.py | 22 ++++++++++++++----- test/SelfGravity/ShearingWave/testme.py | 7 +++--- 2 files changed, 20 insertions(+), 9 deletions(-) diff --git a/test/SelfGravity/ShearingWave/python/testidefix.py b/test/SelfGravity/ShearingWave/python/testidefix.py index 8aff88c06..82ae582a8 100755 --- a/test/SelfGravity/ShearingWave/python/testidefix.py +++ b/test/SelfGravity/ShearingWave/python/testidefix.py @@ -20,7 +20,7 @@ parser = argparse.ArgumentParser() parser.add_argument("-noplot", - default=True, + default=False, help="disable plotting", action="store_true") @@ -56,8 +56,10 @@ def f2s(t,y): kx = kx0 + q*Omega*ky*t k = np.sqrt(kx*kx+ky*ky) + pot = -4*np.pi *G *Sigma0 *isSG + b = -2*q*Omega*kx*ky/k/k - c = 2*(2-q)*Omega*Omega+k**2*cs**2 - 2*np.pi *G *Sigma0 *k*isSG - 2*q*(2-q)*Omega*Omega*ky*ky/k/k + c = 2*(2-q)*Omega*Omega+k**2*cs**2 + pot - 2*q*(2-q)*Omega*Omega*ky*ky/k/k d = 2*Omega*(1-q*ky*ky/k/k)*Sigma0*xi1 sp,s = y @@ -89,16 +91,24 @@ def f2s(t,y): time = np.array(time) -amplitude=np.array(amplitude).real +amplitude=np.array(amplitude).real/Sigma0 amp_exact=np.array(amp_exact) - -err=np.abs((amp_exact*Sigma0-amplitude)/Sigma0).max() +err=np.sum((amp_exact-amplitude)**2)/np.sum(amp_exact**2) print("Error=",err) -if(err<1e-3): + +if not args.noplot: + import matplotlib.pyplot as plt + fig, ax = plt.subplots() + ax.plot(sol.t, sol.y[1],color="orange",ls="-",label="linear") + ax.plot(time,amplitude,".",label="vtk") + ax.legend() + plt.show() + +if(err<.1): print("SUCCESS") sys.exit(0) else: diff --git a/test/SelfGravity/ShearingWave/testme.py b/test/SelfGravity/ShearingWave/testme.py index c4128efc8..7199fb5ca 100755 --- a/test/SelfGravity/ShearingWave/testme.py +++ b/test/SelfGravity/ShearingWave/testme.py @@ -9,11 +9,11 @@ sys.path.append(os.getenv("IDEFIX_DIR")) import pytools.idfx_test as tst -tolerance=1e-15 +tolerance=1e-12 def testMe(test): test.configure() test.compile() - inifiles=["idefix.ini"] + inifiles=["idefix.ini"]#,"idefix-noSG.ini"] for ini in inifiles: test.run(inputFile=ini) @@ -24,6 +24,8 @@ def testMe(test): test=tst.idfxTest() +test.noplot = False + if not test.dec: test.dec=['8','1'] @@ -33,7 +35,6 @@ def testMe(test): else: testMe(test) else: - test.noplot = True test.single=False test.reconstruction=2 test.mpi=False From b1d1242e7d91a756c1c823e7eaacc89df648d777 Mon Sep 17 00:00:00 2001 From: marc Date: Mon, 5 May 2025 17:57:01 +0200 Subject: [PATCH 17/28] Use Roe --- test/SelfGravity/ShearingWave/idefix.ini | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/test/SelfGravity/ShearingWave/idefix.ini b/test/SelfGravity/ShearingWave/idefix.ini index ad3c1dfa8..62a152e9f 100644 --- a/test/SelfGravity/ShearingWave/idefix.ini +++ b/test/SelfGravity/ShearingWave/idefix.ini @@ -10,7 +10,7 @@ first_dt 1e-4 nstages 2 [Hydro] -solver hllc +solver roe rotation 1. shearingBox -1.5 csiso constant 0.07853981633974483 @@ -44,4 +44,3 @@ X3-end periodic vtk 0.05 dmp 0.2 log 100 -# uservar phiP From cc097223c62e8f80a59bf197a29e3eee554a6a23 Mon Sep 17 00:00:00 2001 From: marc Date: Fri, 15 Nov 2024 10:29:22 +0100 Subject: [PATCH 18/28] test enforce periodic From 2c5fd24a339f6bc2faf049b35c5c9b28463c8de1 Mon Sep 17 00:00:00 2001 From: marc Date: Fri, 15 Nov 2024 11:34:02 +0100 Subject: [PATCH 19/28] Tiref test --- src/gravity/laplacian.cpp | 2 ++ src/gravity/laplacian.hpp | 2 ++ 2 files changed, 4 insertions(+) diff --git a/src/gravity/laplacian.cpp b/src/gravity/laplacian.cpp index c049db524..4970a66de 100644 --- a/src/gravity/laplacian.cpp +++ b/src/gravity/laplacian.cpp @@ -564,6 +564,8 @@ void Laplacian::EnforceBoundary(int dir, BoundarySide side, LaplacianBoundaryTyp const int nprocsIDIR = data->mygrid->nproc[dir]; + const int nxiglob = data->mygrid->np_int[IDIR]; + idefix_for("BoundaryShearingBox", kbeg, kend, jbeg, jend, ibeg, iend, KOKKOS_LAMBDA (int k, int j, int i) { const int iscrh = i - side*(ighost +nxi); diff --git a/src/gravity/laplacian.hpp b/src/gravity/laplacian.hpp index 743d11a42..ae65f6c98 100644 --- a/src/gravity/laplacian.hpp +++ b/src/gravity/laplacian.hpp @@ -90,6 +90,8 @@ class Laplacian { bool isTwoPi{false}; bool havePreconditioner{false}; // Use of preconditionner (or not) + + IdefixArray3D sBArray; ///< Array use by shearingbox boundary conditions IdefixArray3D sBArray; ///< Array use by shearingbox boundary conditions From 4405b32373b16631ca9a8443d0a235e2ca90c35c Mon Sep 17 00:00:00 2001 From: marc Date: Mon, 18 Nov 2024 10:12:20 +0100 Subject: [PATCH 20/28] working as intended --- src/gravity/laplacian.hpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/gravity/laplacian.hpp b/src/gravity/laplacian.hpp index ae65f6c98..743d11a42 100644 --- a/src/gravity/laplacian.hpp +++ b/src/gravity/laplacian.hpp @@ -90,8 +90,6 @@ class Laplacian { bool isTwoPi{false}; bool havePreconditioner{false}; // Use of preconditionner (or not) - - IdefixArray3D sBArray; ///< Array use by shearingbox boundary conditions IdefixArray3D sBArray; ///< Array use by shearingbox boundary conditions From a7909d6d6c37ebb03d483300f00718207b464914 Mon Sep 17 00:00:00 2001 From: marc Date: Fri, 21 Feb 2025 11:25:43 +0100 Subject: [PATCH 21/28] Fixed IDIR index shift --- src/gravity/laplacian.cpp | 6 +++++- src/kokkos | 2 +- test/HD/ShearingBox/python/testidefix.py | 4 ++-- test/SelfGravity/ShearingSheet/analysis.cpp | 17 ++++++++++++++++ test/SelfGravity/ShearingSheet/analysis.hpp | 20 +++++++++++++++++++ .../ShearingSheet/python/testidefix.py | 4 ++++ test/SelfGravity/ShearingSheet/setup.cpp | 18 +++++++++++++++++ 7 files changed, 67 insertions(+), 4 deletions(-) create mode 100644 test/SelfGravity/ShearingSheet/analysis.cpp create mode 100644 test/SelfGravity/ShearingSheet/analysis.hpp diff --git a/src/gravity/laplacian.cpp b/src/gravity/laplacian.cpp index 4970a66de..b628ba0bd 100644 --- a/src/gravity/laplacian.cpp +++ b/src/gravity/laplacian.cpp @@ -534,6 +534,7 @@ void Laplacian::EnforceBoundary(int dir, BoundarySide side, LaplacianBoundaryTyp break; } case shearingbox: { + idfx::cout<<"This is SG shearing box boundary " <mygrid->nproc[dir]; - const int nxiglob = data->mygrid->np_int[IDIR]; + //const int nxiglob = data->mygrid->np_int[IDIR]; idefix_for("BoundaryShearingBox", kbeg, kend, jbeg, jend, ibeg, iend, KOKKOS_LAMBDA (int k, int j, int i) { @@ -582,6 +583,9 @@ void Laplacian::EnforceBoundary(int dir, BoundarySide side, LaplacianBoundaryTyp const int jom2 = jghost + ((jo-2-jghost)%nxj+nxj)%nxj; const int jop2 = jghost + ((jo+2-jghost)%nxj+nxj)%nxj; + //if (j==0) + // idfx::cout<< t <<"\t"< +#include +#include + +Analysis::Analysis() { + } + +/* **************************************************************** */ + +void Analysis::PerformAnalysis(DataBlock &data) { + idfx::pushRegion("Analysis::PerformAnalysis"); + data.DumpToFile("test"); + idfx::popRegion(); +} diff --git a/test/SelfGravity/ShearingSheet/analysis.hpp b/test/SelfGravity/ShearingSheet/analysis.hpp new file mode 100644 index 000000000..f9e412b34 --- /dev/null +++ b/test/SelfGravity/ShearingSheet/analysis.hpp @@ -0,0 +1,20 @@ +#ifndef ANALYSIS_HPP_ +#define ANALYSIS_HPP_ + +#include "dataBlock.hpp" +#include "dataBlockHost.hpp" +#include "grid.hpp" +#include "idefix.hpp" +#include "input.hpp" +#include "output.hpp" +#include +#include + +class Analysis { +public: + // Constructor from Setup arguments + Analysis(); + void PerformAnalysis(DataBlock &); +}; + +#endif // ANALYSIS_HPP__ diff --git a/test/SelfGravity/ShearingSheet/python/testidefix.py b/test/SelfGravity/ShearingSheet/python/testidefix.py index b1b65fab8..bce149b4c 100755 --- a/test/SelfGravity/ShearingSheet/python/testidefix.py +++ b/test/SelfGravity/ShearingSheet/python/testidefix.py @@ -3,7 +3,11 @@ """ Created on Wed Feb 19 15:42:19 2025 +<<<<<<< HEAD @author: vdbma +======= +@author: vdmba +>>>>>>> 4b686185 (Fixed IDIR index shift) """ import os import sys diff --git a/test/SelfGravity/ShearingSheet/setup.cpp b/test/SelfGravity/ShearingSheet/setup.cpp index 3ed35acec..f01c7efaa 100644 --- a/test/SelfGravity/ShearingSheet/setup.cpp +++ b/test/SelfGravity/ShearingSheet/setup.cpp @@ -1,9 +1,18 @@ #include "idefix.hpp" #include "setup.hpp" +<<<<<<< HEAD +======= +#include "analysis.hpp" +>>>>>>> 4b686185 (Fixed IDIR index shift) static real gammaIdeal; static real omega; static real shear; +<<<<<<< HEAD +======= +Analysis *analysis; + +>>>>>>> 4b686185 (Fixed IDIR index shift) //#define STRATIFIED @@ -97,3 +106,12 @@ void Setup::InitFlow(DataBlock &data) { // Send it all, if needed d.SyncToDevice(); } +<<<<<<< HEAD +======= + + +// Analyse data to produce an output + +void MakeAnalysis(DataBlock & data) { +} +>>>>>>> 4b686185 (Fixed IDIR index shift) From e58c4b0f54a6367f9c7794058c40c48d40af83aa Mon Sep 17 00:00:00 2001 From: marc Date: Tue, 25 Feb 2025 13:03:05 +0100 Subject: [PATCH 22/28] Cleanup --- src/gravity/laplacian.cpp | 8 -------- 1 file changed, 8 deletions(-) diff --git a/src/gravity/laplacian.cpp b/src/gravity/laplacian.cpp index b628ba0bd..40d862376 100644 --- a/src/gravity/laplacian.cpp +++ b/src/gravity/laplacian.cpp @@ -534,7 +534,6 @@ void Laplacian::EnforceBoundary(int dir, BoundarySide side, LaplacianBoundaryTyp break; } case shearingbox: { - idfx::cout<<"This is SG shearing box boundary " <mygrid->nproc[dir]; - - //const int nxiglob = data->mygrid->np_int[IDIR]; - idefix_for("BoundaryShearingBox", kbeg, kend, jbeg, jend, ibeg, iend, KOKKOS_LAMBDA (int k, int j, int i) { const int iscrh = i - side*(ighost +nxi); @@ -583,9 +578,6 @@ void Laplacian::EnforceBoundary(int dir, BoundarySide side, LaplacianBoundaryTyp const int jom2 = jghost + ((jo-2-jghost)%nxj+nxj)%nxj; const int jop2 = jghost + ((jo+2-jghost)%nxj+nxj)%nxj; - //if (j==0) - // idfx::cout<< t <<"\t"< Date: Tue, 25 Feb 2025 13:18:10 +0100 Subject: [PATCH 23/28] typo --- test/SelfGravity/ShearingSheet/python/testidefix.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/test/SelfGravity/ShearingSheet/python/testidefix.py b/test/SelfGravity/ShearingSheet/python/testidefix.py index bce149b4c..b1b65fab8 100755 --- a/test/SelfGravity/ShearingSheet/python/testidefix.py +++ b/test/SelfGravity/ShearingSheet/python/testidefix.py @@ -3,11 +3,7 @@ """ Created on Wed Feb 19 15:42:19 2025 -<<<<<<< HEAD @author: vdbma -======= -@author: vdmba ->>>>>>> 4b686185 (Fixed IDIR index shift) """ import os import sys From a88aaccd5a620b0a5340b2e055d13fc5e0ea2538 Mon Sep 17 00:00:00 2001 From: marc Date: Tue, 25 Feb 2025 15:08:33 +0100 Subject: [PATCH 24/28] lambda capture --- src/gravity/laplacian.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/gravity/laplacian.cpp b/src/gravity/laplacian.cpp index 40d862376..c049db524 100644 --- a/src/gravity/laplacian.cpp +++ b/src/gravity/laplacian.cpp @@ -562,6 +562,8 @@ void Laplacian::EnforceBoundary(int dir, BoundarySide side, LaplacianBoundaryTyp const real eps = dL/dy - m; + const int nprocsIDIR = data->mygrid->nproc[dir]; + idefix_for("BoundaryShearingBox", kbeg, kend, jbeg, jend, ibeg, iend, KOKKOS_LAMBDA (int k, int j, int i) { const int iscrh = i - side*(ighost +nxi); From 72324dae198cffabb0557eba3760857a206f3143 Mon Sep 17 00:00:00 2001 From: marc Date: Mon, 17 Mar 2025 15:13:38 +0100 Subject: [PATCH 25/28] Revert kokkos version change --- src/kokkos | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/kokkos b/src/kokkos index 6ecdf605e..09e775bfc 160000 --- a/src/kokkos +++ b/src/kokkos @@ -1 +1 @@ -Subproject commit 6ecdf605e0f7639adec599d25cf0e206d7b8f9f5 +Subproject commit 09e775bfc585840bb9ab1156cbd8d7d1c9e0cc6d From a23af806f23d5c1d85e5f7ca6159e056d005ab7e Mon Sep 17 00:00:00 2001 From: marc Date: Mon, 17 Mar 2025 15:20:07 +0100 Subject: [PATCH 26/28] Revert accidental minor changes --- test/HD/ShearingBox/python/testidefix.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/HD/ShearingBox/python/testidefix.py b/test/HD/ShearingBox/python/testidefix.py index 3584c951b..1a0c2509a 100755 --- a/test/HD/ShearingBox/python/testidefix.py +++ b/test/HD/ShearingBox/python/testidefix.py @@ -1,9 +1,9 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- """ -Created on Wed Feb 19 15:42:19 2025 +Created on Mon Jun 21 15:42:19 2021 -@author: vdmba +@author: lesurg """ import sys import numpy as np From 21b02f2bbbe783ecff3a276955d2b44822107604 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Mon, 17 Mar 2025 22:09:42 -0400 Subject: [PATCH 27/28] properly initialise coordinates --- test/SelfGravity/ShearingSheet/setup.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/SelfGravity/ShearingSheet/setup.cpp b/test/SelfGravity/ShearingSheet/setup.cpp index f01c7efaa..8425355eb 100644 --- a/test/SelfGravity/ShearingSheet/setup.cpp +++ b/test/SelfGravity/ShearingSheet/setup.cpp @@ -73,7 +73,7 @@ Setup::Setup(Input &input, Grid &grid, DataBlock &data, Output &output) { void Setup::InitFlow(DataBlock &data) { // Create a host copy DataBlockHost d(data); - real x,y,z; + real x,y,z = 0; real cs2 = gammaIdeal*omega*omega; From 29cbeae351349bf453bb9c8d8d535433da45b180 Mon Sep 17 00:00:00 2001 From: marc Date: Mon, 16 Jun 2025 13:09:17 +0200 Subject: [PATCH 28/28] remove analysis --- test/SelfGravity/ShearingSheet/analysis.cpp | 17 ----------------- test/SelfGravity/ShearingSheet/analysis.hpp | 20 -------------------- test/SelfGravity/ShearingSheet/setup.cpp | 18 ------------------ 3 files changed, 55 deletions(-) delete mode 100644 test/SelfGravity/ShearingSheet/analysis.cpp delete mode 100644 test/SelfGravity/ShearingSheet/analysis.hpp diff --git a/test/SelfGravity/ShearingSheet/analysis.cpp b/test/SelfGravity/ShearingSheet/analysis.cpp deleted file mode 100644 index 743373ce1..000000000 --- a/test/SelfGravity/ShearingSheet/analysis.cpp +++ /dev/null @@ -1,17 +0,0 @@ -#include "analysis.hpp" -#include "fluid.hpp" -#include "idefix.hpp" -#include -#include -#include - -Analysis::Analysis() { - } - -/* **************************************************************** */ - -void Analysis::PerformAnalysis(DataBlock &data) { - idfx::pushRegion("Analysis::PerformAnalysis"); - data.DumpToFile("test"); - idfx::popRegion(); -} diff --git a/test/SelfGravity/ShearingSheet/analysis.hpp b/test/SelfGravity/ShearingSheet/analysis.hpp deleted file mode 100644 index f9e412b34..000000000 --- a/test/SelfGravity/ShearingSheet/analysis.hpp +++ /dev/null @@ -1,20 +0,0 @@ -#ifndef ANALYSIS_HPP_ -#define ANALYSIS_HPP_ - -#include "dataBlock.hpp" -#include "dataBlockHost.hpp" -#include "grid.hpp" -#include "idefix.hpp" -#include "input.hpp" -#include "output.hpp" -#include -#include - -class Analysis { -public: - // Constructor from Setup arguments - Analysis(); - void PerformAnalysis(DataBlock &); -}; - -#endif // ANALYSIS_HPP__ diff --git a/test/SelfGravity/ShearingSheet/setup.cpp b/test/SelfGravity/ShearingSheet/setup.cpp index 8425355eb..39e9f6fd2 100644 --- a/test/SelfGravity/ShearingSheet/setup.cpp +++ b/test/SelfGravity/ShearingSheet/setup.cpp @@ -1,18 +1,9 @@ #include "idefix.hpp" #include "setup.hpp" -<<<<<<< HEAD -======= -#include "analysis.hpp" ->>>>>>> 4b686185 (Fixed IDIR index shift) static real gammaIdeal; static real omega; static real shear; -<<<<<<< HEAD -======= -Analysis *analysis; - ->>>>>>> 4b686185 (Fixed IDIR index shift) //#define STRATIFIED @@ -106,12 +97,3 @@ void Setup::InitFlow(DataBlock &data) { // Send it all, if needed d.SyncToDevice(); } -<<<<<<< HEAD -======= - - -// Analyse data to produce an output - -void MakeAnalysis(DataBlock & data) { -} ->>>>>>> 4b686185 (Fixed IDIR index shift)