Shearing box self gravity

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.

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<LaplacianBoundaryType,3> leftBound,
@@ -39,12 +40,23 @@ Laplacian::Laplacian(DataBlock *datain, std::array<LaplacianBoundaryType,3> 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<real>("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<real> 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<int> (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

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<real> sBArray;    ///< Array use by shearingbox boundary conditions
 
   DataBlock *data;
 

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;

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

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

test/SelfGravity/ShearingSheet/definitions.hpp

@@ -0,0 +1,6 @@
+#define COMPONENTS 2
+#define DIMENSIONS 2
+
+#define GEOMETRY CARTESIAN
+//#define DEBUG_BICGSTAB
+//#define ISOTHERMAL

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

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)