Draft: Lagrangian MC particle pusher

src/athena.hpp

@@ -52,7 +52,8 @@ enum VariableIndex {IDN=0, IM1=1, IVX=1, IM2=2, IVY=2, IM3=3, IVZ=3, IEN=4, ITM=
 // array indices for components of magnetic field
 enum BFieldIndex {IBX=0, IBY=1, IBZ=2};
 // array indices for particle arrays
-enum ParticlesIndex {PGID=0, PTAG=1, IPX=0, IPVX=1, IPY=2, IPVY=3, IPZ=4, IPVZ=5};
+enum ParticlesIndex {PGID=0, PTAG=1, PLASTMOVE=2, PLASTLEVEL=3,
+                     IPX=0, IPY=1, IPZ=2, IPVX=3, IPVY=4, IPVZ=5};
 
 // integer constants to specify spatial reconstruction methods
 enum ReconstructionMethod {dc, plm, ppm4, ppmx, wenoz};

src/bvals/bvals.hpp

@@ -247,8 +247,9 @@ class ParticlesBoundaryValues {
   ParticlesBoundaryValues(particles::Particles *ppart, ParameterInput *pin);
   ~ParticlesBoundaryValues();
 
-  int nprtcl_send, nprtcl_recv;
+  int nprtcl_send, nprtcl_recv, nprtcl_destroy;
   DualArray1D<ParticleLocationData> sendlist;
+  DualArray1D<ParticleLocationData> destroylist;
 
   // Data needed to count number of messages and particles to send between ranks
   int nsends; // number of MPI sends to neighboring ranks on this rank

src/hydro/hydro.cpp

@@ -32,6 +32,7 @@ Hydro::Hydro(MeshBlockPack *ppack, ParameterInput *pin) :
     coarse_w0("cprim",1,1,1,1,1),
     u1("cons1",1,1,1,1,1),
     uflx("uflx",1,1,1,1,1),
+    uflxidnsaved("uflxidn",1,1,1,1),
     utest("utest",1,1,1,1,1),
     fofc("fofc",1,1,1,1) {
   // Total number of MeshBlocks on this rank to be used in array dimensioning
@@ -287,4 +288,22 @@ Hydro::~Hydro() {
   if (psrc != nullptr) {delete psrc;}
 }
 
+//----------------------------------------------------------------------------------------
+// SetSaveUFlxIdn:  set flag to save density flux over entire step
+
+void Hydro::SetSaveUFlxIdn() {
+  int nmb = std::max((pmy_pack->nmb_thispack), (pmy_pack->pmesh->nmb_maxperrank));
+  auto &indcs = pmy_pack->pmesh->mb_indcs;
+  int ncells1 = indcs.nx1 + 2*(indcs.ng);
+  int ncells2 = (indcs.nx2 > 1)? (indcs.nx2 + 2*(indcs.ng)) : 1;
+  int ncells3 = (indcs.nx3 > 1)? (indcs.nx3 + 2*(indcs.ng)) : 1;
+
+  // allocated saved arrays for density flux
+  Kokkos::realloc(uflxidnsaved.x1f, nmb, ncells3, ncells2, ncells1+1);
+  Kokkos::realloc(uflxidnsaved.x2f, nmb, ncells3, ncells2+1, ncells1);
+  Kokkos::realloc(uflxidnsaved.x3f, nmb, ncells3+1, ncells2, ncells1);
+
+  uflxidn_saved = true;
+}
+
 } // namespace hydro

src/hydro/hydro.hpp

@@ -41,6 +41,7 @@ struct HydroTaskIDs {
   TaskID sendf;
   TaskID recvf;
   TaskID expl;
+  TaskID savef;
   TaskID restu;
   TaskID sendu;
   TaskID recvu;
@@ -88,6 +89,10 @@ class Hydro {
   DvceFaceFld5D<Real> uflx;   // fluxes of conserved quantities on cell faces
   Real dtnew;
 
+  // following used to save step-to-step flux values
+  bool uflxidn_saved = false;
+  DvceFaceFld4D<Real> uflxidnsaved;
+
   // following used for FOFC
   DvceArray4D<bool> fofc;  // flag for each cell to indicate if FOFC is needed
   bool use_fofc = false;   // flag to enable FOFC
@@ -96,6 +101,7 @@ class Hydro {
   HydroTaskIDs id;
 
   // functions...
+  void SetSaveUFlxIdn();
   void AssembleHydroTasks(std::map<std::string, std::shared_ptr<TaskList>> tl);
   // ...in "before_stagen_tl" list
   TaskStatus InitRecv(Driver *d, int stage);
@@ -105,6 +111,7 @@ class Hydro {
   TaskStatus SendFlux(Driver *d, int stage);
   TaskStatus RecvFlux(Driver *d, int stage);
   TaskStatus ExpRKUpdate(Driver *d, int stage);
+  TaskStatus SaveFlux(Driver *d, int stage);
   TaskStatus RestrictU(Driver *d, int stage);
   TaskStatus SendU(Driver *d, int stage);
   TaskStatus RecvU(Driver *d, int stage);

src/hydro/hydro_fluxes.cpp

@@ -356,4 +356,63 @@ template void Hydro::CalculateFluxes<Hydro_RSolver::hllc_sr>(Driver *pdriver, in
 template void Hydro::CalculateFluxes<Hydro_RSolver::llf_gr>(Driver *pdriver, int stage);
 template void Hydro::CalculateFluxes<Hydro_RSolver::hlle_gr>(Driver *pdriver, int stage);
 
+//----------------------------------------------------------------------------------------
+//! \fn TaskStatus Hydro::SaveFlux
+//! \brief Save IDN fluxes across all stages for full time step
+
+TaskStatus Hydro::SaveFlux(Driver *pdrive, int stage) {
+  if (uflxidn_saved) {
+    auto &indcs = pmy_pack->pmesh->mb_indcs;
+    int is = indcs.is, ie = indcs.ie;
+    int js = indcs.js, je = indcs.je;
+    int ks = indcs.ks, ke = indcs.ke;
+    bool &multi_d = pmy_pack->pmesh->multi_d;
+    bool &three_d = pmy_pack->pmesh->three_d;
+
+    Real dt = pmy_pack->pmesh->dt;
+    int nmb1 = pmy_pack->nmb_thispack - 1;
+    auto flx1 = uflx.x1f;
+    auto flx2 = uflx.x2f;
+    auto flx3 = uflx.x3f;
+    auto flxidn1 = uflxidnsaved.x1f;
+    auto flxidn2 = uflxidnsaved.x2f;
+    auto flxidn3 = uflxidnsaved.x3f;
+    auto &mbsize = pmy_pack->pmb->mb_size;
+
+    Real dtfactor = dt / 2.0;
+
+    // GNW 2024-JUL-5: Warning .. be careful when using this different time integrators
+    par_for("flux_save",DevExeSpace(),0,nmb1,ks,ke+1,js,je+1,is,ie+1,
+    KOKKOS_LAMBDA(const int m, const int k, const int j, const int i) {
+      // save dF1/dx1
+      if (j <= je && k <= ke)  {
+        if (stage == 1) {
+          flxidn1(m,k,j,i) = flx1(m,IDN,k,j,i) / mbsize.d_view(m).dx1 * dtfactor;
+        } else {
+          flxidn1(m,k,j,i) += flx1(m,IDN,k,j,i) / mbsize.d_view(m).dx1 * dtfactor;
+        }
+      }
+
+      // save dF2/dx2
+      if (multi_d && i <= ie && k <= ke) {
+        if (stage == 1) {
+          flxidn2(m,k,j,i) = flx2(m,IDN,k,j,i) / mbsize.d_view(m).dx2 * dtfactor;
+        } else {
+          flxidn2(m,k,j,i) += flx2(m,IDN,k,j,i) / mbsize.d_view(m).dx2 * dtfactor;
+        }
+      }
+
+      // save dF3/dx3
+      if (three_d && i <= ie && j <= je) {
+        if (stage == 1) {
+          flxidn3(m,k,j,i) = flx3(m,IDN,k,j,i) / mbsize.d_view(m).dx3 * dtfactor;
+        } else {
+          flxidn3(m,k,j,i) += flx3(m,IDN,k,j,i) / mbsize.d_view(m).dx3 * dtfactor;
+        }
+      }
+    });
+  }
+  return TaskStatus::complete;
+}
+
 } // namespace hydro

src/hydro/hydro_tasks.cpp

@@ -53,7 +53,8 @@ void Hydro::AssembleHydroTasks(std::map<std::string, std::shared_ptr<TaskList>>
   id.sendf = tl["stagen"]->AddTask(&Hydro::SendFlux, this, id.flux);
   id.recvf = tl["stagen"]->AddTask(&Hydro::RecvFlux, this, id.sendf);
   id.expl  = tl["stagen"]->AddTask(&Hydro::ExpRKUpdate, this, id.recvf);
-  id.restu = tl["stagen"]->AddTask(&Hydro::RestrictU, this, id.expl);
+  id.savef = tl["stagen"]->AddTask(&Hydro::SaveFlux, this, id.expl);
+  id.restu = tl["stagen"]->AddTask(&Hydro::RestrictU, this, id.savef);
   id.sendu = tl["stagen"]->AddTask(&Hydro::SendU, this, id.restu);
   id.recvu = tl["stagen"]->AddTask(&Hydro::RecvU, this, id.sendu);
   id.bcs   = tl["stagen"]->AddTask(&Hydro::ApplyPhysicalBCs, this, id.recvu);

src/main.cpp

@@ -295,6 +295,7 @@ int main(int argc, char *argv[]) {
     pmesh->pgen = std::make_unique<ProblemGenerator>(pinput, pmesh, restartfile);
     restartfile.Close();
   }
+  pmesh->FinalizeParticleDataStructures(pinput);
 
   //--- Step 6. --------------------------------------------------------------------------
   // Construct Driver and Outputs. Actual outputs (including initial conditions) are made

src/mesh/mesh.cpp

@@ -613,7 +613,12 @@ void Mesh::AddCoordinatesAndPhysics(ParameterInput *pinput) {
     pmb_pack->AddCoordinates(pinput);
     pmb_pack->AddPhysics(pinput);
   }
+}
+
+//----------------------------------------------------------------------------------------
+// \fn Mesh::FinalizeParticleDataStructures
 
+void Mesh::FinalizeParticleDataStructures(ParameterInput *pinput) {
   // Determine total number of particles across all ranks
   particles::Particles *ppart = pmb_pack->ppart;
   if (ppart != nullptr) {

src/mesh/mesh.hpp

@@ -151,6 +151,7 @@ class Mesh {
   void WriteMeshStructure();
   void NewTimeStep(const Real tlim);
   void AddCoordinatesAndPhysics(ParameterInput *pinput);
+  void FinalizeParticleDataStructures(ParameterInput *pinput);
   BoundaryFlag GetBoundaryFlag(const std::string& input_string);
   std::string GetBoundaryString(BoundaryFlag input_flag);
 

src/mhd/mhd.cpp

@@ -40,6 +40,7 @@ MHD::MHD(MeshBlockPack *ppack, ParameterInput *pin) :
     efld("efld",1,1,1,1),
     wsaved("wsaved",1,1,1,1,1),
     bccsaved("bccsaved",1,1,1,1,1),
+    uflxidnsaved("uflxidn",1,1,1,1),
     e3x1("e3x1",1,1,1,1),
     e2x1("e2x1",1,1,1,1),
     e1x2("e1x2",1,1,1,1),
@@ -357,4 +358,22 @@ void MHD::SetSaveWBcc() {
   wbcc_saved = true;
 }
 
+//----------------------------------------------------------------------------------------
+// SetSaveUFlxIdn:  set flag to save density flux over entire step
+
+void MHD::SetSaveUFlxIdn() {
+  int nmb = std::max((pmy_pack->nmb_thispack), (pmy_pack->pmesh->nmb_maxperrank));
+  auto &indcs = pmy_pack->pmesh->mb_indcs;
+  int ncells1 = indcs.nx1 + 2*(indcs.ng);
+  int ncells2 = (indcs.nx2 > 1)? (indcs.nx2 + 2*(indcs.ng)) : 1;
+  int ncells3 = (indcs.nx3 > 1)? (indcs.nx3 + 2*(indcs.ng)) : 1;
+
+  // allocated saved arrays for density flux
+  Kokkos::realloc(uflxidnsaved.x1f, nmb, ncells3, ncells2, ncells1+1);
+  Kokkos::realloc(uflxidnsaved.x2f, nmb, ncells3, ncells2+1, ncells1);
+  Kokkos::realloc(uflxidnsaved.x3f, nmb, ncells3+1, ncells2, ncells1);
+
+  uflxidn_saved = true;
+}
+
 } // namespace mhd

src/mhd/mhd.hpp

@@ -48,6 +48,7 @@ struct MHDTaskIDs {
   TaskID sendf;
   TaskID recvf;
   TaskID expl;
+  TaskID savef;
   TaskID restu;
   TaskID sendu;
   TaskID recvu;
@@ -115,6 +116,10 @@ class MHD {
   DvceArray5D<Real> wsaved;
   DvceArray5D<Real> bccsaved;
 
+  // following used to save step-to-step flux values
+  bool uflxidn_saved = false;
+  DvceFaceFld4D<Real> uflxidnsaved;
+
   // following used for FOFC algorithm
   DvceArray4D<bool> fofc;  // flag for each cell to indicate if FOFC is needed
   bool use_fofc = false;   // flag to enable FOFC
@@ -124,6 +129,7 @@ class MHD {
 
   // functions...
   void SetSaveWBcc();
+  void SetSaveUFlxIdn();
   void AssembleMHDTasks(std::map<std::string, std::shared_ptr<TaskList>> tl);
   // ...in "before_timeintegrator" task list
   TaskStatus SaveMHDState(Driver *d, int stage);
@@ -135,6 +141,7 @@ class MHD {
   TaskStatus SendFlux(Driver *d, int stage);
   TaskStatus RecvFlux(Driver *d, int stage);
   TaskStatus ExpRKUpdate(Driver *d, int stage);
+  TaskStatus SaveFlux(Driver *d, int stage);
   TaskStatus RestrictU(Driver *d, int stage);
   TaskStatus SendU(Driver *d, int stage);
   TaskStatus RecvU(Driver *d, int stage);

src/mhd/mhd_fluxes.cpp

@@ -417,4 +417,63 @@ template void MHD::CalculateFluxes<MHD_RSolver::hlle_sr>(Driver *pdriver, int st
 template void MHD::CalculateFluxes<MHD_RSolver::llf_gr>(Driver *pdriver, int stage);
 template void MHD::CalculateFluxes<MHD_RSolver::hlle_gr>(Driver *pdriver, int stage);
 
+//----------------------------------------------------------------------------------------
+//! \fn TaskStatus MHD::SaveFlux
+//! \brief Save IDN fluxes across all stages for full time step
+
+TaskStatus MHD::SaveFlux(Driver *pdrive, int stage) {
+  if (uflxidn_saved) {
+    auto &indcs = pmy_pack->pmesh->mb_indcs;
+    int is = indcs.is, ie = indcs.ie;
+    int js = indcs.js, je = indcs.je;
+    int ks = indcs.ks, ke = indcs.ke;
+    bool &multi_d = pmy_pack->pmesh->multi_d;
+    bool &three_d = pmy_pack->pmesh->three_d;
+
+    Real dt = pmy_pack->pmesh->dt;
+    int nmb1 = pmy_pack->nmb_thispack - 1;
+    auto flx1 = uflx.x1f;
+    auto flx2 = uflx.x2f;
+    auto flx3 = uflx.x3f;
+    auto flxidn1 = uflxidnsaved.x1f;
+    auto flxidn2 = uflxidnsaved.x2f;
+    auto flxidn3 = uflxidnsaved.x3f;
+    auto &mbsize = pmy_pack->pmb->mb_size;
+
+    Real dtfactor = dt / 2.0;
+
+    // GNW 2024-JUL-5: Warning .. be careful when using this different time integrators
+    par_for("flux_save",DevExeSpace(),0,nmb1,ks,ke+1,js,je+1,is,ie+1,
+    KOKKOS_LAMBDA(const int m, const int k, const int j, const int i) {
+      // save dF1/dx1
+      if (j <= je && k <= ke)  {
+        if (stage == 1) {
+          flxidn1(m,k,j,i) = flx1(m,IDN,k,j,i) / mbsize.d_view(m).dx1 * dtfactor;
+        } else {
+          flxidn1(m,k,j,i) += flx1(m,IDN,k,j,i) / mbsize.d_view(m).dx1 * dtfactor;
+        }
+      }
+
+      // save dF2/dx2
+      if (multi_d && i <= ie && k <= ke) {
+        if (stage == 1) {
+          flxidn2(m,k,j,i) = flx2(m,IDN,k,j,i) / mbsize.d_view(m).dx2 * dtfactor;
+        } else {
+          flxidn2(m,k,j,i) += flx2(m,IDN,k,j,i) / mbsize.d_view(m).dx2 * dtfactor;
+        }
+      }
+
+      // save dF3/dx3
+      if (three_d && i <= ie && j <= je) {
+        if (stage == 1) {
+          flxidn3(m,k,j,i) = flx3(m,IDN,k,j,i) / mbsize.d_view(m).dx3 * dtfactor;
+        } else {
+          flxidn3(m,k,j,i) += flx3(m,IDN,k,j,i) / mbsize.d_view(m).dx3 * dtfactor;
+        }
+      }
+    });
+  }
+  return TaskStatus::complete;
+}
+
 } // namespace mhd

src/mhd/mhd_tasks.cpp

@@ -45,7 +45,8 @@ void MHD::AssembleMHDTasks(std::map<std::string, std::shared_ptr<TaskList>> tl)
   id.sendf = tl["stagen"]->AddTask(&MHD::SendFlux, this, id.flux);
   id.recvf = tl["stagen"]->AddTask(&MHD::RecvFlux, this, id.sendf);
   id.expl  = tl["stagen"]->AddTask(&MHD::ExpRKUpdate, this, id.recvf);
-  id.restu = tl["stagen"]->AddTask(&MHD::RestrictU, this, id.expl);
+  id.savef = tl["stagen"]->AddTask(&MHD::SaveFlux, this, id.expl);
+  id.restu = tl["stagen"]->AddTask(&MHD::RestrictU, this, id.savef);
   id.sendu = tl["stagen"]->AddTask(&MHD::SendU, this, id.restu);
   id.recvu = tl["stagen"]->AddTask(&MHD::RecvU, this, id.sendu);
   id.efld  = tl["stagen"]->AddTask(&MHD::CornerE, this, id.recvu);

src/outputs/vtk_prtcl.cpp

@@ -182,15 +182,24 @@ void ParticleVTKOutput::WriteOutputFile(Mesh *pm, ParameterInput *pin) {
   }
 
   // Write Part 6: scalar particle data
-  // Write gid of points
+  bool have_written_pointdata_header = false;
+
   for (int n=0; n<(pm->pmb_pack->ppart->nidata); ++n) {
     std::stringstream msg;
+
+    if (!have_written_pointdata_header) {
+      have_written_pointdata_header = true;
+      msg << std::endl << std::endl << "POINT_DATA " << npout_total << std::endl;
+    }
+
     if (n == static_cast<int>(PGID)) {
-      msg << std::endl << std::endl << "POINT_DATA " << npout_total << std::endl
-          << "SCALARS gid float" << std::endl << "LOOKUP_TABLE default" << std::endl;
+      msg << std::endl << "SCALARS gid float" << std::endl
+          << "LOOKUP_TABLE default" << std::endl;
     } else if (n == static_cast<int>(PTAG)) {
-      msg << std::endl << std::endl << "POINT_DATA " << npout_total << std::endl
-          << "SCALARS ptag float" << std::endl << "LOOKUP_TABLE default" << std::endl;
+      msg << std::endl << "SCALARS ptag float" << std::endl
+          << "LOOKUP_TABLE default" << std::endl;
+    } else if (n == static_cast<int>(PLASTMOVE) || n == static_cast<int>(PLASTLEVEL)) {
+      continue;
     }
     if (global_variable::my_rank == 0) {
       partfile.Write_any_type_at(msg.str().c_str(),msg.str().size(),header_offset,"byte");