Coagulation

src/fluid/boundary/axis.cpp

@@ -38,15 +38,25 @@ void Axis::SymmetrizeEx1Side(int jref, IdefixArray3D<real> Ex1) {
 
   auto Ex1Avg = this->Ex1Avg.deviceView();
 
-  idefix_for("Ex1_ini",0,data->np_tot[IDIR],
-      KOKKOS_LAMBDA(int i) {
-        Ex1Avg(i) = ZERO_F;
-      });
-
-  idefix_for("Ex1_Symmetrize",data->beg[KDIR],data->end[KDIR],0,data->np_tot[IDIR],
-    KOKKOS_LAMBDA(int k,int i) {
-      Kokkos::atomic_add(&Ex1Avg(i),  Ex1(k,jref,i));
+  // Deterministic reduction: one thread per i, fixed k loop order
+
+  const int kbeg = data->beg[KDIR];
+  const int kend = data->end[KDIR];
+
+  idefix_for("Ex1_SymmetrizeDet", 0, data->np_tot[IDIR],
+    KOKKOS_LAMBDA(int i) {
+      real sum = ZERO_F;
+      real c = ZERO_F;  // Kahan compensation
+      for(int k = kbeg; k < kend; k++) {
+        real y = Ex1(k, jref, i) - c;
+        real t = sum + y;
+        c = (t - sum) - y;
+        sum = t;
+      }
+      Ex1Avg(i) = sum;
     });
+
+
   if(needMPIExchange) {
     #ifdef WITH_MPI
       Kokkos::fence();
@@ -105,14 +115,20 @@ void Axis::RegularizeCurrentSide(int side) {
     IdefixArray1D<real> dx3 = data->dx[KDIR];
     IdefixArray1D<real> dx2 = data->dx[JDIR];
 
-    idefix_for("B_ini",0,data->np_tot[IDIR],
-          KOKKOS_LAMBDA(int i) {
-            BAvg(i) = ZERO_F;
-    });
-    idefix_for("Compute_Bcirculation",data->beg[KDIR],data->end[KDIR],0,data->np_tot[IDIR],
-        KOKKOS_LAMBDA(int k,int i) {
-          Kokkos::atomic_add(&BAvg(i), Vs(BX3s,k,jc,i)*dx3(k) ); // Compute the circulation of
-                                                                 // Bphi around the pole
+    const int kbeg = data->beg[KDIR];
+    const int kend = data->end[KDIR];
+    idefix_for("Compute_BcirculationDet", 0, data->np_tot[IDIR],
+      KOKKOS_LAMBDA(int i) {
+      real sum = ZERO_F;
+      real c = ZERO_F;  // Kahan compensation
+      for(int k = kbeg; k < kend; k++) {
+        real term = Vs(BX3s, k, jc, i) * dx3(k);
+        real y = term - c;
+        real t = sum + y;
+        c = (t - sum) - y;
+        sum = t;
+      }
+      BAvg(i) = sum;
     });
 
     if(needMPIExchange) {
@@ -202,19 +218,34 @@ void Axis::FixBx2sAxis(int side) {
       sign = -1;
     }
 
-    idefix_for("B_ini",0,data->np_tot[IDIR],0,2,
-          KOKKOS_LAMBDA(int i, int n) {
-            BAvg(i,n) = ZERO_F;
-    });
-    idefix_for("BHorizontal_compute",data->beg[KDIR],data->end[KDIR],0,data->np_tot[IDIR],
-        KOKKOS_LAMBDA(int k,int i) {
-          real Bthmid = sign*HALF_F*(Vs(BX2s,k,jaxe-1,i) + Vs(BX2s,k,jaxe+1,i));
-          real Bphimid = HALF_F*(Vs(BX3s,k,jin,i) + Vs(BX3s,k,jout,i));
-          //Bthmid = 0.0;
-          //Bphimid = 0.0;
+    const int kbeg = data->beg[KDIR];
+    const int kend = data->end[KDIR];
+    idefix_for("BHorizontal_computeDet", 0, data->np_tot[IDIR],
+    KOKKOS_LAMBDA(int i) {
+      real sumX = ZERO_F;
+      real cX   = ZERO_F;  // Kahan compensation for X
+      real sumY = ZERO_F;
+      real cY   = ZERO_F;  // Kahan compensation for Y
+
+      for(int k = kbeg; k < kend ; k++) {
+        real Bthmid  = sign*HALF_F*(Vs(BX2s,k,jaxe-1,i) + Vs(BX2s,k,jaxe+1,i));
+        real Bphimid =      HALF_F*(Vs(BX3s,k,jin,i)    + Vs(BX3s,k,jout,i));
+
+        real termX = Bthmid * cos(phi(k)) - Bphimid * sin(phi(k));
+        real yX = termX - cX;
+        real tX = sumX + yX;
+        cX = (tX - sumX) - yX;
+        sumX = tX;
+
+        real termY = Bthmid * sin(phi(k)) + Bphimid * cos(phi(k));
+        real yY = termY - cY;
+        real tY = sumY + yY;
+        cY = (tY - sumY) - yY;
+        sumY = tY;
+      }
 
-          Kokkos::atomic_add(&BAvg(i,IDIR), Bthmid * cos(phi(k)) - Bphimid * sin(phi(k)));
-          Kokkos::atomic_add(&BAvg(i,JDIR), Bthmid * sin(phi(k)) + Bphimid * cos(phi(k)));
+      BAvg(i,IDIR) = sumX;
+      BAvg(i,JDIR) = sumY;
     });
     if(needMPIExchange) {
       Kokkos::fence();

test.py

@@ -27,6 +27,12 @@
 sys.path.append(source_dir)
 from pytools.idfx_test_run import IdexPytestRunner
 
+# export it for sub-commands accessing ./pytools module
+if "PYTHONPATH" in os.environ:
+    os.environ["PYTHONPATH"] = source_dir + ":" + os.environ["PYTHONPATH"]
+else:
+    os.environ["PYTHONPATH"] = source_dir
+
 # should be global so it remember state (last build) and avoids
 # to build for each run if we just changed the ini file and run options.
 gblIdefixPytestRunner = IdexPytestRunner(__file__)

test/Dust/CoagBox/CMakeLists.txt

@@ -0,0 +1,3 @@
+enable_idefix_property(Idefix_PYTHON)
+add_idefix_source(coagulation.cpp)
+add_idefix_source(coagulation.hpp)

test/Dust/CoagBox/README.md

@@ -0,0 +1,47 @@
+# Testing coagulation in Idefix
+
+## Introduction 
+
+We solve the coagulation of dust grains using the code COALA (see details in [Lombart & Laibe, 2021, MNRAS](https://ui.adsabs.harvard.edu/abs/2021MNRAS.501.4298L/abstract)). The code is a high-order scheme based on Discontinuous Galerkin (DG) method to solve the Smoluchowski coagulation equation in its conservative form. In practice, the coagulation corresponds to a source term in Idefix that updates the dust densities (and not the momentum for the moment). In this source term, COALA provides a table that only depends on the mass grid (not on the time), thanks to the method used in COALA. The table contains the needed weights to evolve coagulation in time. Hence, we generate the table only once, before the simulation starts, then it is read during the run in the coagulation source term.
+
+## Requierments
+
+Make sure you have python>=3.10 (relative to COALA and these tests) and pybind11>=2.12 (relative to pydefix that is used here, see [Pydefix doc](https://idefix.readthedocs.io/latest/modules/pydefix.html)). You can create a virtual environement with the packages in `requirements.txt`, which is the one I am using right now.
+
+## The tests explained
+
+### Coagulation kernels
+
+We can use COALA with 4 different coagulation kernels, which define the nature of the interaction between grains. The precalculation table generated by COALA depends on the kernel used. The *constant kernel* and the *additive kernel* are the only with exact solutions. The *brownian kernel* and the *collision kernel* derive from the ballistic kernel and are known as the "physical kernels". The ballistic kernel depends on the cross-section of the collision and the differential velocity between grains. In the brownian kernel case, the analytical expression of the differential velocity from Brownian motion is used, in order to obtain a reference solution. In the collision kernel case, only the cross-section is used for precalculation, and the differential velocities have to be given by the hydro part, Idefix. 
+
+### Description of the tests and the files
+
+We prepare 4 input files for each coagulation kernel. The precalculation table is already written are and ready to be called by Idefix. The corresponding COALA solutions are also already computed for comparison. Hence the tests can be done without using COALA. We give here a brief description of the files:
+
+- `idefix-XXXKernel.ini`: The input files. Parameters needed for coagulation are `nSpecies`in the *dust* block and all the parameters in the *Coala* block. Do not change them!
+- `coagulation.hpp` and `coagulation.cpp`: the source files concerning the coagulation. This is where the precalculation tables are loaded, and where the coagulation source term is written. 
+- `setup.cpp`: the setup, calling the coagulation files. It also contains a function that check mass conservation.
+- `myInit.py`: To initialize the flow.
+- `testidefix.py`: To compute the standard tests (the plots and the errors), comparing Idefix to COALA's solutions. 
+- `testme.json`: To compute the reference tests, given a space of parameters to test.
+- `utils.py`: Common functions and classes, called by the `.py` files.
+- `computeCoalaTable.py`: To generate the precalculation table. Depends on `nSpecies`, `Q`, `kernel`, `kpol`, `massmin` and `massmax`.
+- `coala_coag.py`: The file calling COALA to generate the solution for comparison.
+
+The output files are stored in the `./data` directory. 
+
+>[!NOTE]
+> If you want to change any parameter, you need COALA to compute the new precalculation table and the solution for comparison (so don't do it!).
+
+>[!NOTE]
+> To use `testidefix.py` (but also `computeCoalaTable.py`, and `coala_coag.py`), you need to give the argument `-ini idefix-XXXKernel.ini` (replace by the input file you want to test).
+
+## How do you test ?
+### Standard tests
+-  First, initialise the IDEFIX_DIR environement to the root of the idefix directory, with `export IDEFIX_DIR=/path_to_idefix_root`;
+- Second, configure the code with `cmake $IDEFIX_DIR` and compile the code with `make -j`;
+- Thrid, do `export PYTHONPATH=/$MY_VIRTUAL_ENV/lib/python3.XX/site-packages:$PYTHONPATH` (replace XX by the python version in your $MY_VIRTUAL_ENV) and launch the code with `./idefix -i idefix-XXXKernel.ini`(replace `idefix-XXXKernel.ini` by the input file you want to test);
+- Finally, use `python testidefix.py -ini idefix-XXXKernel.ini` to test the run.
+
+### Reference tests
+- In the root of the idefix directory, use `./test.py -subdir /path_to_the_coagulation_test -v -x` (-v to enable verbose mode, -x to stop the test if one error is detected)

test/Dust/CoagBox/coagulation.cpp

@@ -0,0 +1,176 @@
+#include "idefix.hpp"
+#include "coagulation.hpp"
+#include "npy.hpp"
+
+Coagulation *CoagulationGlob;
+
+Coagulation::Coagulation(DataBlock &data, int kernel, real eps, real coagCFL) {
+    this->kernel = kernel;
+    this->eps = eps;
+    this->coagCFL = coagCFL;
+
+    //==============================================//
+    //==============================================//
+    // In this part, we load the coagulation tables //
+    //==============================================//
+    //==============================================//
+    const int nSpecies = data.dust.size();
+    //===== Load the mass grid =====// 
+    std::string filename_massgrid = "data/test-massgrid-kernel-"+std::to_string(kernel)+".npy";
+    std::vector<uint64_t> shape_massgrid;
+    bool fortran_order_massgrid;
+    std::vector<double> massgrid;
+    npy::LoadArrayFromNumpy(filename_massgrid, shape_massgrid, fortran_order_massgrid, massgrid);
+      //===== verification =====//
+    if(massgrid.size() != static_cast<size_t>(nSpecies+1)) {
+      IDEFIX_ERROR("Mass grid has an unexpected size");
+    }
+      //========================//
+
+    this->binWidth = IdefixArray1D<real>("binWidth", nSpecies);
+    IdefixHostArray1D<real> binWidthHost = Kokkos::create_mirror_view(this->binWidth);
+    for(int n = 0; n <nSpecies; n++) {
+      binWidthHost(n) = massgrid[n+1] - massgrid[n];
+    }
+    Kokkos::deep_copy(this->binWidth, binWidthHost);
+
+    //===== Loading the Flux Table =====//
+    std::string filename_coagTabFlux = "data/test-coagTabFlux-kernel-"+std::to_string(kernel)+".npy";
+    std::vector<uint64_t> shape_coagTabFlux;
+    bool fortran_order_coagTabFlux;
+    std::vector<double> coagTabFlux;
+    npy::LoadArrayFromNumpy(filename_coagTabFlux, shape_coagTabFlux, fortran_order_coagTabFlux, coagTabFlux);
+      //===== verification =====//
+    const size_t expectedSize = static_cast<size_t>(nSpecies)*static_cast<size_t>(nSpecies)*static_cast<size_t>(nSpecies);
+    if(coagTabFlux.size() != expectedSize) {
+      IDEFIX_ERROR("Coagulation flux table has an unexpected size");
+    }
+      //========================//
+
+    this->coagTabFlux = IdefixArray3D<real>("coagTabFlux", nSpecies, nSpecies, nSpecies);
+    IdefixHostArray3D<real> coagTabFluxHost = Kokkos::create_mirror_view(this->coagTabFlux);
+    for(int i=0; i<nSpecies; i++) {
+      for(int j=0; j<nSpecies; j++) {
+        for(int k=0; k<nSpecies; k++) {
+          int index = i*nSpecies*nSpecies + j*nSpecies + k;
+          coagTabFluxHost(i,j,k) = coagTabFlux[index];
+        }
+      }
+    }
+    Kokkos::deep_copy(this->coagTabFlux, coagTabFluxHost);
+}
+
+void CoagSourceTerm(DataBlock &dataRef, const real t, const real dtin) {
+  //==============================================//
+  /*
+  Coagulation source term : compute the coagulation product to finally upload the dust densities. 
+  This function is implemented via operator splitting : coagulation is integrated in a step.
+  Hence, enroll the function with `data.EnrollUserStepFirst(&CoagSourceTerm);`.
+  This version does not preserve bins momentum yet.
+  This version does not work on GPU yet.
+  */
+  //==============================================//
+  DataBlock *data = &dataRef;
+  data->PrimToCons();
+  IdefixArray1D<real> h = CoagulationGlob->binWidth;
+  IdefixArray3D<real> coagTabFlux = CoagulationGlob->coagTabFlux;
+  const int nSpecies = data->dust.size();
+  const int test_size = 20;  // Because we can't give nSpecies in the KOKKOS_LAMBDA
+
+  const real eps = CoagulationGlob->eps;
+  const real coagCFL = CoagulationGlob->coagCFL;
+  const bool physicalKernel = (CoagulationGlob->kernel==3);
+
+  idefix_for("CoagSourceTerm",0,data->np_tot[KDIR],0,data->np_tot[JDIR],0,data->np_tot[IDIR],
+            KOKKOS_LAMBDA (int k, int j, int i) {
+              real rho[test_size];
+              real g0[test_size];
+              real vx1[test_size];
+              real vx2[test_size];
+              real vx3[test_size];
+              real deltaV[test_size][test_size];
+              real Flux[test_size+1];
+
+              // Initial State
+              for(int n=0; n<nSpecies; n++) {
+                rho[n] = data->dust[n]->Uc(RHO,k,j,i);
+                vx1[n] = data->dust[n]->Uc(MX1,k,j,i)/rho[n];
+                vx2[n] = data->dust[n]->Uc(MX2,k,j,i)/rho[n];
+                vx3[n] = data->dust[n]->Uc(MX3,k,j,i)/rho[n];
+
+                /*
+                // To have a constant drift velocity that does not evolve in sub-cycling
+                // WARNING: dirty here! check that it corresponds to Compute_test_velocity in utils.py
+                real vx = 1e-2 + (1e-1 - 1e-2)*n/(nSpecies - 1);
+                real vy = 1.3*vx;
+                real vz = 1.7*vx;
+                vx1[n] = vx;
+                vx2[n] = vy;
+                vx3[n] = vz;
+                */
+              }
+
+              // Relative velocity
+              if(physicalKernel) {
+                for(int lp=0; lp<nSpecies; lp++) {
+                  for(int l=0; l<nSpecies; l++) {
+                    real dvx = vx1[lp] - vx1[l];
+                    real dvy = vx2[lp] - vx2[l];
+                    real dvz = vx3[lp] - vx3[l];
+                    deltaV[lp][l] = sqrt(dvx*dvx + dvy*dvy + dvz*dvz);
+                  }
+                }
+              }
+
+              //===== COAGULATION SUBCYCLING =====//
+              real time = 0.0;
+              while(time<dtin) {
+                Flux[0] = 0.0;  // No mass flux crossing the left boundary of the distribution
+                for(int n=0; n<nSpecies; n++) {
+                  g0[n] = rho[n]/h(n);
+                }
+                //===== FLUX CALCULATION =====//
+                for(int n=0; n<nSpecies; n++) {
+                  real flux = 0.0;
+                  for(int lp=0; lp <= n; lp++) {
+                    real g0_lp = g0[lp];
+                    for(int l=0; l<nSpecies; l++) {
+                      real g0_l = g0[l];
+                      real delta_v = physicalKernel ? deltaV[lp][l] : 1.0;
+                      flux += g0_lp*g0_l*coagTabFlux(n,lp,l)*delta_v;
+                    }
+                  }
+                  Flux[n+1] = flux;
+                }
+                //===== CFL CALCULATION =====//
+                real dt_sub = std::numeric_limits<double>::max();
+                for(int n=0; n<nSpecies; n++) {
+                  real dflux = Flux[n+1] - Flux[n];
+                  if(std::abs(dflux)>0.0 && rho[n]>eps*h(n)) {
+                    dt_sub = std::min(dt_sub, std::abs(rho[n]/dflux));
+                  }
+                }
+                real dt_coag = dt_sub*coagCFL;
+                dt_coag = std::min(dt_coag, dtin-time);  // Last substep
+                if(dt_coag<=0.0) {
+                  IDEFIX_ERROR("Coagulation time step is non-positive");
+                }
+                //===== EULER UPDATE =====//
+                for(int n=0; n<nSpecies; n++) {
+                  rho[n] = std::max(rho[n] - dt_coag*(Flux[n+1] - Flux[n]), eps*h(n));
+                }
+                time += dt_coag;
+              }
+              //===== UPDATE THE FINAL DENSITIES =====//
+              for(int n=0; n<nSpecies; n++) {
+                data->dust[n]->Uc(RHO,k,j,i) = rho[n];
+
+                // To keep a constant velocity, just for the test with kernel = 3
+                data->dust[n]->Uc(MX1,k,j,i) = rho[n]*vx1[n];
+                data->dust[n]->Uc(MX2,k,j,i) = rho[n]*vx2[n];
+                data->dust[n]->Uc(MX3,k,j,i) = rho[n]*vx3[n];
+              }
+            }
+  );
+  data->ConsToPrim();
+}

test/Dust/CoagBox/coagulation.hpp

@@ -0,0 +1,20 @@
+#include "idefix.hpp"
+#include "dataBlock/dataBlock.hpp"
+#include "fluid/fluid.hpp"
+
+class Coagulation {
+ public:
+  // Constructor: allocate and load all coagulation tables
+  Coagulation(DataBlock &data, int kernel, real eps, real coagCFL);
+
+  IdefixArray1D<real> binWidth;
+  IdefixArray3D<real> coagTabFlux;
+
+  int kernel;
+  real eps;
+  real coagCFL;
+};
+
+extern Coagulation *CoagulationGlob;
+
+void CoagSourceTerm(DataBlock &data, const real t, const real dtin);