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. diff --git a/src/gravity/laplacian.cpp b/src/gravity/laplacian.cpp index 2300db82d..c049db524 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,12 +40,23 @@ 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; } + #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 @@ -521,7 +533,100 @@ 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; + + // 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/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 = (nprocsIDIR==1) ? + ighost + (i+ighost*(nxi-1))%nxi : i; + + 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; + + real Fl,Fr; + real dqm, dqp, dq; + + // 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; + //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 { + //Compute Fl + dqm = localVar(k,jo,iref) - localVar(k,jom1,iref); + 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); + + 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) { + const int iscrh = i - side*(ighost +nxi); + localVar(k,j,i) = scrh(k,j,iscrh); + }); + 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..743d11a42 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, @@ -90,6 +91,7 @@ class Laplacian { bool isTwoPi{false}; bool havePreconditioner{false}; // Use of preconditionner (or not) + IdefixArray3D sBArray; ///< Array use by shearingbox boundary conditions DataBlock *data; 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; diff --git a/test/HD/ShearingBox/idefix-fargo.ini b/test/HD/ShearingBox/idefix-fargo.ini index 9d0dd373e..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 511d7e15c..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/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..3508d6f91 --- /dev/null +++ b/test/SelfGravity/ShearingSheet/idefix.ini @@ -0,0 +1,46 @@ +[Grid] +X1-grid 1 -.5 60 u .5 +X2-grid 1 -.5 60 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 0.2 +# 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..b1b65fab8 --- /dev/null +++ b/test/SelfGravity/ShearingSheet/python/testidefix.py @@ -0,0 +1,88 @@ +#!/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 matplotlib.pyplot as plt +import argparse + +# 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() + + +Vbeg = readVTK("../data.0000.vtk") +Vend = readVTK("../data.0002.vtk") + +Ly = Vbeg.yl[-1]-Vbeg.yl[0] +xLeft = Vbeg.x[0] +xRight = Vbeg.x[-1] + +deltaT = Vend.t[0]-Vbeg.t[0] + +conf = inifix.load("../idefix.ini") + +S = conf["Hydro"]["shearingBox"] + +dy = Ly/Vbeg.ny + + +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 not args.noplot: + fig, ax = plt.subplots() + 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(Vend.data["phiP"][0,:,0]) +indexMaxEndRightPot = np.argmin(Vend.data["phiP"][-1,:,0]) + +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(errorPot/dy) +print("Error=",err) + +if(err<1): # total cummulated shift is less than on cell + 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..39e9f6fd2 --- /dev/null +++ b/test/SelfGravity/ShearingSheet/setup.cpp @@ -0,0 +1,99 @@ +#include "idefix.hpp" +#include "setup.hpp" + +static real gammaIdeal; +static real omega; +static real shear; + +//#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) { + 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); + + 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; + + 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.2 * 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(); +} diff --git a/test/SelfGravity/ShearingSheet/testme.py b/test/SelfGravity/ShearingSheet/testme.py new file mode 100755 index 000000000..9da97a788 --- /dev/null +++ b/test/SelfGravity/ShearingSheet/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=['3','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) 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..62a152e9f --- /dev/null +++ b/test/SelfGravity/ShearingWave/idefix.ini @@ -0,0 +1,46 @@ +[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 roe +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 diff --git a/test/SelfGravity/ShearingWave/python/testidefix.py b/test/SelfGravity/ShearingWave/python/testidefix.py new file mode 100755 index 000000000..82ae582a8 --- /dev/null +++ b/test/SelfGravity/ShearingWave/python/testidefix.py @@ -0,0 +1,116 @@ +#!/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=False, + 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) + + 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 + 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 + 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/Sigma0 +amp_exact=np.array(amp_exact) + + + +err=np.sum((amp_exact-amplitude)**2)/np.sum(amp_exact**2) +print("Error=",err) + + +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: + 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..3f88a4651 --- /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); + +} + +// 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 + // 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; + real Sigma0 = 1./40.; + real Q = 1.; + + 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++) { + 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;, + d.Vc(VX2,k,j,i) = shear*x;, + 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..7199fb5ca --- /dev/null +++ b/test/SelfGravity/ShearingWave/testme.py @@ -0,0 +1,44 @@ +#!/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-12 +def testMe(test): + test.configure() + test.compile() + inifiles=["idefix.ini"]#,"idefix-noSG.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() +test.noplot = False + +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.single=False + test.reconstruction=2 + test.mpi=False + testMe(test) + + test.mpi=True + testMe(test)