added kernels to compute structure functions online, but postprocessi…

…ng hdf5 save files works better, and to compute sgs energy fluxes
cselab · Oct 18, 2020 · de9a319 · de9a319
1 parent 3272540
commit de9a319
Show file tree

Hide file tree

Showing 18 changed files with 431 additions and 107 deletions.
diff --git a/makefiles/make.daint b/makefiles/make.daint
@@ -2,7 +2,7 @@ CXX=CC
 LD=CC
 
 ACCFFT_ROOT=../accfft/build/
-accfft?=true
+accfft?=false
 
 ifeq "$(compiler)" "intel"
 #CPPFLAGS += -dynamic
@@ -19,4 +19,4 @@ ifeq "$(accfft)" "true"
 LIBS += $(CRAY_CUDATOOLKIT_POST_LINK_OPTS)
 endif
 
-#CPPFLAGS += -DCUP_ASYNC_DUMP
+CPPFLAGS += -DCUP_ASYNC_DUMP
diff --git a/source/Definitions.h b/source/Definitions.h
@@ -36,6 +36,8 @@
 
 CubismUP_3D_NAMESPACE_BEGIN
 
+//#define ENERGY_FLUX_SPECTRUM 1
+
 enum { FE_CHI = 0, FE_U, FE_V, FE_W, FE_P, FE_TMPU, FE_TMPV, FE_TMPW };
 struct FluidElement
 {

diff --git a/source/Simulation.cpp b/source/Simulation.cpp
@@ -164,7 +164,7 @@ void Simulation::_init(const bool restart)
   else
     _ic();
   MPI_Barrier(sim.app_comm);
-  _serialize("init");
+  //_serialize("init");
 
   assert(sim.obstacle_vector != nullptr);
   if (sim.rank == 0)
@@ -397,6 +397,7 @@ void Simulation::setupOperators()
       checkpointPostVelocity = new Checkpoint(sim, "PostVelocity"));
 
   //sim.pipeline.push_back(new HITfiltering(sim));
+  sim.pipeline.push_back(new StructureFunctions(sim));
 
   sim.pipeline.push_back(new Analysis(sim));
 

diff --git a/source/SimulationData.cpp b/source/SimulationData.cpp
@@ -75,7 +75,8 @@ SimulationData::SimulationData(MPI_Comm mpicomm, ArgumentParser &parser)
   bChannelFixedMassFlux = parser("-channelFixedMassFlux").asBool(false);
 
   bRungeKutta23 = parser("-RungeKutta23").asBool(false);
-  bAdvection3rdOrder = parser("-Advection3rdOrder").asBool(true);
+  //bAdvection3rdOrder = parser("-Advection3rdOrder").asBool(true);
+  bAdvection3rdOrder = parser("-Advection3rdOrder").asBool(false);
 
   uMax_forced = parser("-uMax_forced").asDouble(0.0);
 

diff --git a/source/main_RL_HIT.cpp b/source/main_RL_HIT.cpp
@@ -23,7 +23,6 @@
 #include <sstream>
 
 using Real = cubismup3d::Real;
-static constexpr Real rew_baseline = - 1000.0;
 #define GIVE_LOCL_STATEVARS
 #define GIVE_GRID_STATEVARS
 
@@ -108,7 +107,7 @@ inline void app_main(
   // If every grid point is an agent: probably will allocate too much memory
   // and crash because smarties allocates a trajectory for each point
   // If only one agent: sequences will be garbled together and cannot
-  // send clean Sequences.
+  // send clean Episodes to smarties' learning algorithms.
   // Also, rememebr that separate agents are thread safe!
   // let's say that each fluid block has one agent
   const int nAgentPerBlock = 1;
@@ -172,10 +171,11 @@ inline void app_main(
     const Real timeUpdateLES = tau_eta / LES_RL_FREQ_A;
     const Real timeSimulationMax = LES_RL_N_TSIM * tau_integral;
     const int maxNumUpdatesPerSim = timeSimulationMax / timeUpdateLES;
-    if (0) { // enable all dumping. //  && wrank != 1
+    if (bEvaluating) { // enable all dumping. //  && wrank != 1
       sim.sim.b3Ddump = true;  sim.sim.muteAll  = false;
       sim.sim.b2Ddump = false; sim.sim.saveFreq = 0;
-      sim.sim.verbose = true;  sim.sim.saveTime = timeUpdateLES;
+      //sim.sim.verbose = true;  sim.sim.saveTime = timeUpdateLES;
+      sim.sim.verbose = true;  sim.sim.saveTime = tau_integral;
     }
     printf("Reset simulation up to time=0 with SGS for eps:%f nu:%f Re:%f. "
            "Max %d action turns per simulation.\n", target.eps,

diff --git a/source/operators/AdvectionDiffusion.cpp b/source/operators/AdvectionDiffusion.cpp
@@ -314,7 +314,8 @@ struct Upwind3rd
   const Real invU = 1 / std::max(EPS, _sim.uMax_measured);
   Upwind3rd(const SimulationData& s) : _sim(s) {}
 
-  template<StepType step, int dir> Real diffx(LabMPI& L, const FluidBlock& o, const Real uAbs[3], const int ix, const int iy, const int iz) const {
+  template<StepType step, int dir>
+  inline Real diffx(LabMPI& L, const FluidBlock& o, const Real uAbs[3], const int ix, const int iy, const int iz) const {
     const Real ucc = inp<step,dir>(L,ix,iy,iz);
     const Real um1 = inp<step,dir>(L,ix-1,iy,iz), um2 = inp<step,dir>(L,ix-2,iy,iz);
     const Real up1 = inp<step,dir>(L,ix+1,iy,iz), up2 = inp<step,dir>(L,ix+2,iy,iz);
@@ -329,7 +330,8 @@ struct Upwind3rd
       return uAbs[0]>0? 2*up1 +3*ucc -6*um1 +um2 : -up2 +6*up1 -3*ucc -2*um1;
     #endif
   }
-  template<StepType step, int dir> Real diffy(LabMPI& L, const FluidBlock& o, const Real uAbs[3], const int ix, const int iy, const int iz) const {
+  template<StepType step, int dir>
+  inline Real diffy(LabMPI& L, const FluidBlock& o, const Real uAbs[3], const int ix, const int iy, const int iz) const {
     const Real ucc = inp<step,dir>(L,ix,iy,iz);
     const Real um1 = inp<step,dir>(L,ix,iy-1,iz), um2 = inp<step,dir>(L,ix,iy-2,iz);
     const Real up1 = inp<step,dir>(L,ix,iy+1,iz), up2 = inp<step,dir>(L,ix,iy+2,iz);
@@ -344,7 +346,8 @@ struct Upwind3rd
       return uAbs[1]>0? 2*up1 +3*ucc -6*um1 +um2 : -up2 +6*up1 -3*ucc -2*um1;
     #endif
   }
-  template<StepType step, int dir> Real diffz(LabMPI& L, const FluidBlock& o, const Real uAbs[3], const int ix, const int iy, const int iz) const {
+  template<StepType step, int dir>
+  inline Real diffz(LabMPI& L, const FluidBlock& o, const Real uAbs[3], const int ix, const int iy, const int iz) const {
     const Real ucc = inp<step,dir>(L,ix,iy,iz);
     const Real um1 = inp<step,dir>(L,ix,iy,iz-1), um2 = inp<step,dir>(L,ix,iy,iz-2);
     const Real up1 = inp<step,dir>(L,ix,iy,iz+1), up2 = inp<step,dir>(L,ix,iy,iz+2);
@@ -359,7 +362,8 @@ struct Upwind3rd
       return uAbs[2]>0? 2*up1 +3*ucc -6*um1 +um2 : -up2 +6*up1 -3*ucc -2*um1;
     #endif
   }
-  template<StepType step, int dir> Real   lap(LabMPI& L, const FluidBlock& o, const int ix, const int iy, const int iz) const {
+  template<StepType step, int dir>
+  inline Real   lap(LabMPI& L, const FluidBlock& o, const int ix, const int iy, const int iz) const {
     return  inp<step,dir>(L,ix+1,iy,iz) + inp<step,dir>(L,ix-1,iy,iz)
           + inp<step,dir>(L,ix,iy+1,iz) + inp<step,dir>(L,ix,iy-1,iz)
           + inp<step,dir>(L,ix,iy,iz+1) + inp<step,dir>(L,ix,iy,iz-1)
@@ -532,7 +536,7 @@ void AdvectionDiffusion::operator()(const double dt)
   else
   {
     if(sim.bRungeKutta23) {
-      sim.startProfiler("AdvDiff Kernel");
+      sim.startProfiler("AdvDiff23 Kernel");
       if(sim.bAdvection3rdOrder) {
         const KernelAdvectDiffuse<RK1, Upwind3rd> K1(sim);
         compute(K1);
@@ -545,10 +549,12 @@ void AdvectionDiffusion::operator()(const double dt)
         compute(K2);
       }
       sim.stopProfiler();
-      sim.startProfiler("AdvDiff copy");
-      const UpdateAndCorrectInflow U(sim);
-      U.operate<false>();
-      sim.stopProfiler();
+      if(not sim.bUseFourierBC) {
+        sim.startProfiler("AdvDiff copy");
+        const UpdateAndCorrectInflow U(sim);
+        U.operate<false>();
+        sim.stopProfiler();
+      }
     } else {
       sim.startProfiler("AdvDiff Kernel");
       if(sim.bAdvection3rdOrder) {

diff --git a/source/operators/HITfiltering.cpp b/source/operators/HITfiltering.cpp
@@ -143,8 +143,9 @@ void HITfiltering::operator()(const double dt)
 {
   if ( not (sim.timeAnalysis>0 && (sim.time+dt) >= sim.nextAnalysisTime) )
     return;
+  if (sim.muteAll) return;
 
-  sim.startProfiler("SGS Kernel");
+  sim.startProfiler("HITfiltering Kernel");
   const int BPDX = sim.grid->getBlocksPerDimension(0);
   const int BPDY = sim.grid->getBlocksPerDimension(1);
   const int BPDZ = sim.grid->getBlocksPerDimension(2);
@@ -246,4 +247,164 @@ void HITfiltering::operator()(const double dt)
   check("SGS");
 }
 
+StructureFunctions::StructureFunctions(SimulationData& s) : Operator(s)
+{
+  std::random_device rd;
+  gen.seed(rd());
+}
+
+inline Real periodic_distance(const Real x1, const Real x0, const Real extent)
+{
+  const Real dx = x1 - x0;
+  if      (dx >   extent/2) return dx - extent;
+  else if (dx <= -extent/2) return dx + extent;
+  else return dx;
+}
+inline Real periodic_distance(const std::array<Real,3> & p1,
+                     const std::array<Real,3> & p0,
+                     const std::array<Real,3> & extent)
+{
+  const Real dx = periodic_distance(p1[0], p0[0], extent[0]);
+  const Real dy = periodic_distance(p1[1], p0[1], extent[1]);
+  const Real dz = periodic_distance(p1[2], p0[2], extent[2]);
+  return std::sqrt(dx*dx + dy*dy + dz*dz);
+}
+
+std::array<double, 6> StructureFunctions::pick_ref_point()
+{
+  std::uniform_int_distribution<int> distrib_ranks(0, sim.nprocs-1);
+  std::uniform_int_distribution<size_t> distrib_block(0, vInfo.size()-1);
+  std::uniform_int_distribution<int> distrib_elem(0, CUP_BLOCK_SIZE-1);
+  int ref_rank = distrib_ranks(gen);
+  MPI_Bcast(&ref_rank, 1, MPI_INT, 0, sim.app_comm);
+  const size_t ref_bid  = distrib_block(gen);
+  const int ref_iz = distrib_elem(gen);
+  const int ref_iy = distrib_elem(gen);
+  const int ref_ix = distrib_elem(gen);
+  const BlockInfo & ref_info = vInfo[ref_bid];
+  FluidBlock & ref_block = * (FluidBlock*) ref_info.ptrBlock;
+  const FluidElement & ref_elem = ref_block(ref_ix, ref_iy, ref_iz);
+  const std::array<Real,3> ref_pos = ref_info.pos<Real>(ref_ix, ref_iy, ref_iz);
+  const std::array<Real,3> ref_vel = {ref_elem.u, ref_elem.v, ref_elem.w};
+  std::array<double, 6> ref = {0};
+  if(sim.rank == ref_rank) {
+    ref[0] = ref_pos[0];
+    ref[1] = ref_pos[1];
+    ref[2] = ref_pos[2];
+    ref[3] = ref_vel[0];
+    ref[4] = ref_vel[1];
+    ref[5] = ref_vel[2];
+  }
+  MPI_Allreduce(MPI_IN_PLACE, ref.data(), 6, MPI_DOUBLE, MPI_SUM, sim.app_comm);
+  if(sim.rank == ref_rank) {
+    assert(std::fabs(ref[0] - ref_pos[0]) < 1e-8);
+    assert(std::fabs(ref[1] - ref_pos[1]) < 1e-8);
+    assert(std::fabs(ref[2] - ref_pos[2]) < 1e-8);
+    assert(std::fabs(ref[3] - ref_vel[0]) < 1e-8);
+    assert(std::fabs(ref[4] - ref_vel[1]) < 1e-8);
+    assert(std::fabs(ref[5] - ref_vel[2]) < 1e-8);
+  }
+  return ref;
+}
+
+
+void StructureFunctions::operator()(const double dt)
+{
+  if (sim.muteAll) return;
+  if (computeInterval <= 0 or (sim.time+dt) < nextComputeTime)
+    return;
+  nextComputeTime += computeInterval;
+
+  sim.startProfiler("StructureFunctions Kernel");
+
+  auto ref = pick_ref_point();
+  const std::array<Real,3> ref_pos = {ref[0], ref[1], ref[2]};
+  const std::array<Real,3> ref_vel = {ref[3], ref[4], ref[5]};
+
+  static constexpr size_t oneD_ref_gridN = 32; // LES resolution
+  const Real delta_increments = sim.maxextent / oneD_ref_gridN;
+  static constexpr size_t n_shells = oneD_ref_gridN / 2;
+
+  unsigned long counts[n_shells] = {0};
+  double sum_S2[n_shells] = {0.0}, sum_S3[n_shells] = {0.0};
+  double sum_S4[n_shells] = {0.0}, sum_S6[n_shells] = {0.0};
+  double sum_A3[n_shells] = {0.0};
+
+  const auto get_shell_id = [=](const Real delta) {
+    // first shell goes from 0 to 1.5 delta_increments
+    // then from 1.5 to 2.5 and so on, so we can compute with 0:N / eta
+    // NOTE: 'average' radius is = 3/4 * (B^4 - A^4) / (B^3 - A^3)
+    // where B and A are external and internal radius of shell respectively
+    if (delta <= delta_increments * 1.5) return (size_t) 0;
+    const Real delta_nondim = delta / delta_increments; //should be at least 1.5
+    assert(delta_nondim >= 1.5);
+    // in [1.5, 2.5) return 1, in [2.5, 3.5) return 2 and so on:
+    const size_t shell_id = std::max(delta_nondim - 0.5, (Real) 1);
+    return shell_id;
+  };
+  #pragma omp parallel for schedule(static) reduction(+ : counts[:n_shells], \
+                                                          sum_S2[:n_shells], \
+                                                          sum_S3[:n_shells], \
+                                                          sum_S4[:n_shells], \
+                                                          sum_S6[:n_shells], \
+                                                          sum_A3[:n_shells])
+  for (size_t i = 0; i < vInfo.size(); ++i)
+  {
+    const BlockInfo & info = vInfo[i];
+    FluidBlock& block = * (FluidBlock*) info.ptrBlock;
+
+    for (int iz = 0; iz < CUP_BLOCK_SIZE; ++iz)
+    for (int iy = 0; iy < CUP_BLOCK_SIZE; ++iy)
+    for (int ix = 0; ix < CUP_BLOCK_SIZE; ++ix)
+    {
+      const std::array<Real,3> pos = info.pos<Real>(ix, iy, iz);
+      const Real delta = periodic_distance(pos, ref_pos, sim.extent);
+      const size_t shell_id = get_shell_id(delta);
+
+      if (shell_id >= n_shells) continue;
+      const Real rx = pos[0] - ref_pos[0];
+      const Real ry = pos[1] - ref_pos[1];
+      const Real rz = pos[2] - ref_pos[2];
+      const Real rnorm = std::max(std::sqrt(rx*rx + ry*ry + rz*rz),
+                                  std::numeric_limits<Real>::epsilon());
+      const Real ex = rx/rnorm, ey = ry/rnorm, ez = rz/rnorm;
+      const Real du = block(ix,iy,iz).u - ref_vel[0];
+      const Real dv = block(ix,iy,iz).v - ref_vel[1];
+      const Real dw = block(ix,iy,iz).w - ref_vel[2];
+      const Real deltaU = du*ex + dv*ey + dw*ez;
+      counts[shell_id] += 1;
+      sum_S2[shell_id] += std::pow(deltaU, 2);
+      sum_S3[shell_id] += std::pow(deltaU, 3);
+      sum_S4[shell_id] += std::pow(deltaU, 4);
+      sum_S6[shell_id] += std::pow(deltaU, 6);
+      sum_A3[shell_id] += std::pow(std::fabs(deltaU), 3);
+    }
+  }
+
+  MPI_Allreduce(MPI_IN_PLACE, sum_S2, n_shells, MPI_DOUBLE, MPI_SUM, sim.app_comm);
+  MPI_Allreduce(MPI_IN_PLACE, sum_S3, n_shells, MPI_DOUBLE, MPI_SUM, sim.app_comm);
+  MPI_Allreduce(MPI_IN_PLACE, sum_S4, n_shells, MPI_DOUBLE, MPI_SUM, sim.app_comm);
+  MPI_Allreduce(MPI_IN_PLACE, sum_S6, n_shells, MPI_DOUBLE, MPI_SUM, sim.app_comm);
+  MPI_Allreduce(MPI_IN_PLACE, sum_A3, n_shells, MPI_DOUBLE, MPI_SUM, sim.app_comm);
+  MPI_Allreduce(MPI_IN_PLACE, counts, n_shells, MPI_UNSIGNED_LONG, MPI_SUM, sim.app_comm);
+
+  if(sim.rank==0 and not sim.muteAll)
+  {
+    std::vector<double> buffer;
+    buffer.insert(buffer.end(), sum_S2, sum_S2 + n_shells);
+    buffer.insert(buffer.end(), sum_S3, sum_S3 + n_shells);
+    buffer.insert(buffer.end(), sum_S4, sum_S4 + n_shells);
+    buffer.insert(buffer.end(), sum_S6, sum_S6 + n_shells);
+    buffer.insert(buffer.end(), sum_A3, sum_A3 + n_shells);
+    buffer.insert(buffer.end(), counts, counts + n_shells); // to double
+    FILE * pFile = fopen ("structureFunctionsAnalysis.raw", "ab");
+    fwrite (buffer.data(), sizeof(double), buffer.size(), pFile);
+    fflush(pFile); fclose(pFile);
+  }
+
+  sim.stopProfiler();
+
+  check("SGS");
+}
+
 CubismUP_3D_NAMESPACE_END
diff --git a/source/operators/HITfiltering.h b/source/operators/HITfiltering.h
@@ -24,5 +24,21 @@ class HITfiltering : public Operator
   std::string getName() { return "HITfiltering"; }
 };
 
+class StructureFunctions : public Operator
+{
+  std::mt19937 gen;
+  const Real computeInterval = sim.timeAnalysis / 10;
+  Real nextComputeTime = 0;
+
+  std::array<double, 6> pick_ref_point();
+
+public:
+  StructureFunctions(SimulationData& s);
+
+  void operator()(const double dt);
+
+  std::string getName() { return "StructureFunctions"; }
+};
+
 CubismUP_3D_NAMESPACE_END
 #endif