[hvpic-k-devel] Dump infrastructure for (almost) all fields

CMakeLists.txt

@@ -74,10 +74,14 @@ option(VPIC_ENABLE_HDF5 "Enable HDF5 for use during IO. VPIC does not help you i
 
 option(VPIC_ENABLE_VARIABLE_CHARGE "Enable variable charge for particles" OFF)
 
+option(VPIC_ENABLE_EXTERNAL_FORCE "Enable external force for particles" OFF)
+
 option(VPIC_ENABLE_GPU_AWARE_MPI "Enable GPU-aware MPI. Allows GPUs to communicate with each other directly rather than going through the CPU. If VPIC_ENABLE_GPU_AWARE_MPI is enabled but not supported in MPI, the code may crash." OFF)
 
 option(VPIC_ENABLE_HALO_EXCHANGE "Enable halo exchange code path for remote communication. Temporary option until halo exchange code is accurately tested." ON)
 
+option(VPIC_SHAPE_QS "Use quadratic-sum particle shape" OFF)
+
 option(HYB_USE_SEPARATE_PE "Use separate electron energy equation instead of EoS." OFF)
 
 option(HYB_USE_STATIC_E "Use electrostatic hybrid field solver." OFF)
@@ -371,6 +375,12 @@ if (VPIC_ENABLE_VARIABLE_CHARGE)
    message("--     VPIC: Enabled variable charge for particles")
 endif(VPIC_ENABLE_VARIABLE_CHARGE)
 
+if (VPIC_ENABLE_EXTERNAL_FORCE)
+   add_definitions(-DEXTERNAL_FORCE)
+   set(VPIC_CPPFLAGS "${VPIC_CPPFLAGS} -DEXTERNAL_FORCE")
+   message("--     VPIC: Enabled external force for particles")
+endif(VPIC_ENABLE_EXTERNAL_FORCE)
+
 if (VPIC_ENABLE_GPU_AWARE_MPI)
   add_definitions(-DVPIC_ENABLE_GPU_AWARE_MPI)
   message("--     VPIC: Enabled GPU-Aware MPI")
@@ -381,6 +391,16 @@ if (VPIC_ENABLE_HALO_EXCHANGE)
   message("--     VPIC: Enabled Halo exchange path")
 endif(VPIC_ENABLE_HALO_EXCHANGE)
 
+if (VPIC_SHAPE_QS)
+  add_definitions(-DSHAPE_QS)
+  message("--     VPIC: Using quadratic-sum particle shape")
+  set(VPIC_CXX_FLAGS "${VPIC_CXX_FLAGS} -DSHAPE_QS")
+else()
+  add_definitions(-DSHAPE_NGP)
+  message("--     VPIC: Using nearest-grid-point particle shape")
+  set(VPIC_CXX_FLAGS "${VPIC_CXX_FLAGS} -DSHAPE_NGP")
+endif(VPIC_SHAPE_QS)
+
 if (HYB_USE_SEPARATE_PE)
    add_definitions(-DHYB_USE_SEPARATE_PE)
    set(VPIC_CPPFLAGS "${VPIC_CPPFLAGS} -DHYB_USE_SEPARATE_PE")

examples/iaw/NO-SMOOTH/iaw.cxx

@@ -456,7 +456,7 @@ sim_log( "Loading fields" );
    * Convenience functions for simlog output
    *------------------------------------------------------------------------*/
 
-  char varlist[512];
+  char varlist[1024]; // accommodate "allvars" for fields
   create_field_list(varlist, global->fdParams);
 
   sim_log ( "Fields variable list: " << varlist );

examples/iaw/NO-SMOOTH/translateIAW.f90

@@ -511,7 +511,12 @@ program translate
               bz(idxstart(n,1):idxstop(n,1), idxstart(n,2):idxstop(n,2), idxstart(n,3):idxstop(n,3)) = &
                  buffer(2:nc(1)-1,2:nc(2)-1,2:nc(3)-1)
 
-              read(10)buffer   ! skip div_b error
+              read(10)buffer   ! skip pe
+
+              read(10)buffer   ! skip magnetic0
+              read(10)buffer   ! skip
+              read(10)buffer   ! skip
+              read(10)buffer   ! skip te0
 
               read(10)buffer   ! skip tca
               read(10)buffer   ! skip 

examples/iaw/iaw.cxx

@@ -458,7 +458,7 @@ sim_log( "Loading fields" );
    * Convenience functions for simlog output
    *------------------------------------------------------------------------*/
 
-  char varlist[512];
+  char varlist[1024]; // accommodate "allvars" for fields
   create_field_list(varlist, global->fdParams);
 
   sim_log ( "Fields variable list: " << varlist );

examples/iaw/translateIAW.f90

@@ -511,7 +511,12 @@ program translate
               bz(idxstart(n,1):idxstop(n,1), idxstart(n,2):idxstop(n,2), idxstart(n,3):idxstop(n,3)) = &
                  buffer(2:nc(1)-1,2:nc(2)-1,2:nc(3)-1)
 
-              read(10)buffer   ! skip div_b error
+              read(10)buffer   ! skip pe
+
+              read(10)buffer   ! skip magnetic0
+              read(10)buffer   ! skip
+              read(10)buffer   ! skip
+              read(10)buffer   ! skip te0
 
               read(10)buffer   ! skip tca
               read(10)buffer   ! skip 

examples/islands/COLD-ION/islands_hyb.cxx

@@ -586,7 +586,7 @@ sim_log( "Loading fields" );
    * Convenience functions for simlog output
    *------------------------------------------------------------------------*/
 
-  char varlist[512];
+  char varlist[1024]; // accommodate "allvars" for fields
   create_field_list(varlist, global->fdParams);
 
   sim_log ( "Fields variable list: " << varlist );

examples/islands/islands_hyb.cxx

@@ -586,7 +586,7 @@ sim_log( "Loading fields" );
    * Convenience functions for simlog output
    *------------------------------------------------------------------------*/
 
-  char varlist[512];
+  char varlist[1024]; // accommodate "allvars" for fields
   create_field_list(varlist, global->fdParams);
 
   sim_log ( "Fields variable list: " << varlist );

examples/pcai/pcai.cxx

@@ -462,7 +462,7 @@ sim_log( "Loading fields" );
    * Convenience functions for simlog output
    *------------------------------------------------------------------------*/
 
-  char varlist[512];
+  char varlist[1024]; // accommodate "allvars" for fields
   create_field_list(varlist, global->fdParams);
 
   sim_log ( "Fields variable list: " << varlist );

examples/pcai/translate_pcai.f90

@@ -373,14 +373,14 @@ program translate
   fnames(5) = 'data/By'
   fnames(6) = 'data/Bz'
 
-!  fnames(7) = 'data/Uix'
-!  fnames(8) = 'data/Uiy'
-!  fnames(9) = 'data/Uiz'
+  fnames(7) = 'data/Uix'
+  fnames(8) = 'data/Uiy'
+  fnames(9) = 'data/Uiz'
   fnames(10) = 'data/ni'
-  fnames(11) = 'data/Uix'
-  fnames(12) = 'data/Uiy'
-  fnames(13) = 'data/Uiz'
-  fnames(14) = 'data/niold'
+  !fnames(11) = 'data/Uix'
+  !fnames(12) = 'data/Uiy'
+  !fnames(13) = 'data/Uiz'
+  !fnames(14) = 'data/niold'
   fnames(15) = 'data/aniso'
   !fnames(16) = 'data/Pi-xy1'
 
@@ -513,11 +513,17 @@ program translate
               bz(idxstart(n,1):idxstop(n,1), idxstart(n,2):idxstop(n,2), idxstart(n,3):idxstop(n,3)) = &
                  buffer(2:nc(1)-1,2:nc(2)-1,2:nc(3)-1)
 
-              read(10)buffer   ! skip div_b error
-              !read(10)buffer   ! skip tca
-              !read(10)buffer   ! skip 
-              !read(10)buffer   ! skip 
-              !read(10)buffer   ! skip rhob
+              read(10)buffer   ! skip pe
+
+              read(10)buffer   ! skip magnetic0
+              read(10)buffer   ! skip
+              read(10)buffer   ! skip
+              read(10)buffer   ! skip te0
+
+              read(10)buffer   ! skip tca
+              read(10)buffer   ! skip
+              read(10)buffer   ! skip
+              read(10)buffer   ! skip rhob
 
 
               read(10)buffer
@@ -537,23 +543,6 @@ program translate
                  buffer(2:nc(1)-1,2:nc(2)-1,2:nc(3)-1)
 
 
-              read(10)buffer
-              ux(idxstart(n,1):idxstop(n,1), idxstart(n,2):idxstop(n,2), idxstart(n,3):idxstop(n,3)) = &
-                 buffer(2:nc(1)-1,2:nc(2)-1,2:nc(3)-1)
-
-              read(10)buffer
-              uy(idxstart(n,1):idxstop(n,1), idxstart(n,2):idxstop(n,2), idxstart(n,3):idxstop(n,3)) = &
-                 buffer(2:nc(1)-1,2:nc(2)-1,2:nc(3)-1)
-
-              read(10)buffer
-              uz(idxstart(n,1):idxstop(n,1), idxstart(n,2):idxstop(n,2), idxstart(n,3):idxstop(n,3)) = &
-                 buffer(2:nc(1)-1,2:nc(2)-1,2:nc(3)-1)
-
-              read(10)buffer
-              ne(idxstart(n,1):idxstop(n,1), idxstart(n,2):idxstop(n,2), idxstart(n,3):idxstop(n,3)) = &
-                 buffer(2:nc(1)-1,2:nc(2)-1,2:nc(3)-1)
-
-
               close(10)
 
 
@@ -629,14 +618,14 @@ program translate
      call MPI_FILE_WRITE_AT_ALL(fh(4), offset, bx, ht%nx*ht%ny*ht%nz, MPI_REAL4, status, ierror)
      call MPI_FILE_WRITE_AT_ALL(fh(5), offset, by, ht%nx*ht%ny*ht%nz, MPI_REAL4, status, ierror)
      call MPI_FILE_WRITE_AT_ALL(fh(6), offset, bz, ht%nx*ht%ny*ht%nz, MPI_REAL4, status, ierror)
-!     call MPI_FILE_WRITE_AT_ALL(fh(7), offset, jx, ht%nx*ht%ny*ht%nz, MPI_REAL4, status, ierror)
-!     call MPI_FILE_WRITE_AT_ALL(fh(8), offset, jy, ht%nx*ht%ny*ht%nz, MPI_REAL4, status, ierror)
-!     call MPI_FILE_WRITE_AT_ALL(fh(9), offset, jz, ht%nx*ht%ny*ht%nz, MPI_REAL4, status, ierror)
+     call MPI_FILE_WRITE_AT_ALL(fh(7), offset, jx, ht%nx*ht%ny*ht%nz, MPI_REAL4, status, ierror)
+     call MPI_FILE_WRITE_AT_ALL(fh(8), offset, jy, ht%nx*ht%ny*ht%nz, MPI_REAL4, status, ierror)
+     call MPI_FILE_WRITE_AT_ALL(fh(9), offset, jz, ht%nx*ht%ny*ht%nz, MPI_REAL4, status, ierror)
      call MPI_FILE_WRITE_AT_ALL(fh(10), offset, rho, ht%nx*ht%ny*ht%nz, MPI_REAL4, status, ierror)
-     call MPI_FILE_WRITE_AT_ALL(fh(11), offset, ux, ht%nx*ht%ny*ht%nz, MPI_REAL4, status, ierror)
-     call MPI_FILE_WRITE_AT_ALL(fh(12), offset, uy, ht%nx*ht%ny*ht%nz, MPI_REAL4, status, ierror)
-     call MPI_FILE_WRITE_AT_ALL(fh(13), offset, uz, ht%nx*ht%ny*ht%nz, MPI_REAL4, status, ierror)
-     call MPI_FILE_WRITE_AT_ALL(fh(14), offset, ne, ht%nx*ht%ny*ht%nz, MPI_REAL4, status, ierror)
+!     call MPI_FILE_WRITE_AT_ALL(fh(11), offset, ux, ht%nx*ht%ny*ht%nz, MPI_REAL4, status, ierror)
+!     call MPI_FILE_WRITE_AT_ALL(fh(12), offset, uy, ht%nx*ht%ny*ht%nz, MPI_REAL4, status, ierror)
+!     call MPI_FILE_WRITE_AT_ALL(fh(13), offset, uz, ht%nx*ht%ny*ht%nz, MPI_REAL4, status, ierror)
+!     call MPI_FILE_WRITE_AT_ALL(fh(14), offset, ne, ht%nx*ht%ny*ht%nz, MPI_REAL4, status, ierror)
      call MPI_FILE_WRITE_AT_ALL(fh(15), offset, aniso, ht%nx*ht%ny*ht%nz, MPI_REAL4, status, ierror)
 
 

examples/shock/shock-hyb.cxx

@@ -542,7 +542,7 @@ sim_log( "Loading fields" );
    * Convenience functions for simlog output
    *------------------------------------------------------------------------*/
 
-  char varlist[512];
+  char varlist[1024]; // accommodate "allvars" for fields
   create_field_list(varlist, global->fdParams);
 
   sim_log ( "Fields variable list: " << varlist );

src/field_advance/field_advance.cc

@@ -169,6 +169,16 @@ field_array_t::copy_to_host() {
       host_field[i].sy = k_field(i, field_var::sy);
       host_field[i].sz = k_field(i, field_var::sz);
       host_field[i].se = k_field(i, field_var::se);
+
+#ifdef EXTERNAL_FORCE
+      host_field[i].Ex0   = k_field(i, field_var::Ex0);
+      host_field[i].Ey0   = k_field(i, field_var::Ey0);
+      host_field[i].Ez0   = k_field(i, field_var::Ez0);
+
+      host_field[i].Gx0   = k_field(i, field_var::Gx0);
+      host_field[i].Gy0   = k_field(i, field_var::Gy0);
+      host_field[i].Gz0   = k_field(i, field_var::Gz0);
+#endif
 
       host_field[i].ematx = k_field_edge(i, field_edge_var::ematx);
       host_field[i].ematy = k_field_edge(i, field_edge_var::ematy);
@@ -254,6 +264,16 @@ field_array_t::copy_to_device() {
       k_field(i, field_var::sy) = host_field[i].sy;
       k_field(i, field_var::sz) = host_field[i].sz;
       k_field(i, field_var::se) = host_field[i].se;
+
+#ifdef EXTERNAL_FORCE
+      k_field(i, field_var::Ex0)   = host_field[i].Ex0;
+      k_field(i, field_var::Ey0)   = host_field[i].Ey0;
+      k_field(i, field_var::Ez0)   = host_field[i].Ez0;
+
+      k_field(i, field_var::Gx0)   = host_field[i].Gx0;
+      k_field(i, field_var::Gy0)   = host_field[i].Gy0;
+      k_field(i, field_var::Gz0)   = host_field[i].Gz0;
+#endif
 
       k_field_edge(i, field_edge_var::ematx) = host_field[i].ematx;
       k_field_edge(i, field_edge_var::ematy) = host_field[i].ematy;

src/field_advance/field_advance.h

@@ -165,6 +165,10 @@ typedef struct field {
   float pex, pey, pez,div_b_err ;        // pressure etc
   float ux, uy, uz, ue;                  // Electron bulk flow velocity
   float sx, sy, sz, se;                  // Electron momentum and energy sources
+  #ifdef EXTERNAL_FORCE
+  float Ex0, Ey0, Ez0, _pad1; // External electric field (mover only, exclude from Ohm's law)
+  float Gx0, Gy0, Gz0, _pad2; // External gravity field (mover only, exclude from Ohm's law)
+  #endif
   material_id ematx, ematy, ematz, nmat; // Material at edge centers and nodes
   material_id fmatx, fmaty, fmatz, cmat; // Material at face and cell centers
 

src/field_advance/standard/hyb_advance_b.cc

@@ -85,12 +85,25 @@ hyb_advance_b(field_array_t * RESTRICT fa,
   Kokkos::Profiling::pushRegion("HybyridAdvanceB::Smooth_Ion_Moments");
   //Only smooth on the first subcycle
   if (isub==0) {
+#ifndef SHAPE_NGP
+    //fix edge rho/currents for local BCs
+    //adds adjacent ghost rho/current to cell and sets BC ghosts to 0
+    k_hyb_local_adjust_jf(fa, fa->g);
+    //fix edge rho/currents everywhere else
+    //sends ghosts across ranks and adds ghosts to adjacent cells
+    k_begin_remote_edge_hyb_jf(fa, fa->g, *(fa->fb) );
+    k_end_remote_edge_hyb_jf  (fa, fa->g, *(fa->fb) );
+#else
+    // NGP accumulates only on live cells; ghosts empty.
+    // So no need to fix edge rho/currents
+#endif
 
-    //smooth moments
+    //refresh all ghosts
     Kokkos::Profiling::pushRegion("HybyridAdvanceB::Smooth_Ion_Moments::Exchange_JF");
     k_begin_remote_ghost_hyb_jf(fa, fa->g, *(fa->fb) );
     k_end_remote_ghost_hyb_jf  (fa, fa->g, *(fa->fb) );
     Kokkos::Profiling::popRegion();
+    //fix ghosts for local BCs
     Kokkos::Profiling::pushRegion("HybyridAdvanceB::Smooth_Ion_Moments::Apply_Local_Ghost_JF");
     k_hyb_local_ghost_jf  (fa, fa->g);
     Kokkos::Profiling::popRegion();

src/field_advance/standard/local.cc

@@ -1755,6 +1755,9 @@ void k_local_adjust_rhob(field_array_t* fa, const grid_t* g) {
     adjust_rhob<ZXY>(fa, g, 0, 0, 1);
 }
 
+/*****************************************************************************
+ * Hybrid local ghosts
+ *****************************************************************************/
 
 //Hybrid local B field
 
@@ -2107,3 +2110,65 @@ k_hyb_local_ghost_lapl_b( field_array_t      * RESTRICT f,
 
 #undef APPLYLOCAL
 
+/*****************************************************************************
+ * Hybrid local adjusts
+ *****************************************************************************/
+
+//Hybrid adjust jf,rhof
+
+template<typename T> void adjust_hyb_local_jf(int i, int j, int k,
+                                              const int nx, const int ny, const int nz,
+                                              field_array_t* RESTRICT f, const grid_t* g) {}
+
+#define ADJUSTLOCAL(x_,y_,z_)                                                       \
+  int bc = g->bc[BOUNDARY(i,j,k)];                                                  \
+  k_field_t k_field = f->k_f_d;                                                     \
+  Kokkos::MDRangePolicy<Kokkos::Rank<2> > x_##_face({1,1},{n##y_+1,n##z_+1});       \
+  if(bc < 0 || bc >= world_size) {                                                  \
+    int x_ = (i+j+k)<0 ? 1 : n##x_;  /* adjust edges */                             \
+    switch(bc) {                                                                    \
+    case anti_symmetric_fields:                                                     \
+    case symmetric_fields:                                                          \
+      Kokkos::parallel_for("adjust_hyb_local_jf: anti_symmetric_fields",            \
+                           x_##_face, KOKKOS_LAMBDA(const int y_, const int z_) {   \
+          /*  add ghost cell to adjacent local BC cell, then set ghost to 0 */      \
+          k_field(VOXEL(x,y,z,nx,ny,nz), field_var::jfx)  += k_field(VOXEL(x+i,y+j,z+k,nx,ny,nz), field_var::jfx);  \
+          k_field(VOXEL(x,y,z,nx,ny,nz), field_var::jfy)  += k_field(VOXEL(x+i,y+j,z+k,nx,ny,nz), field_var::jfy);  \
+          k_field(VOXEL(x,y,z,nx,ny,nz), field_var::jfz)  += k_field(VOXEL(x+i,y+j,z+k,nx,ny,nz), field_var::jfz);  \
+          k_field(VOXEL(x,y,z,nx,ny,nz), field_var::rhof) += k_field(VOXEL(x+i,y+j,z+k,nx,ny,nz), field_var::rhof); \
+          k_field(VOXEL(x+i,y+j,z+k,nx,ny,nz), field_var::jfx)  = 0;    \
+          k_field(VOXEL(x+i,y+j,z+k,nx,ny,nz), field_var::jfy)  = 0;    \
+          k_field(VOXEL(x+i,y+j,z+k,nx,ny,nz), field_var::jfz)  = 0;    \
+          k_field(VOXEL(x+i,y+j,z+k,nx,ny,nz), field_var::rhof) = 0;    \
+        });                                                             \
+      break;                                                            \
+    default:                                                            \
+      ERROR(("Bad boundary condition encountered."));                   \
+      break;                                                            \
+    }                                                                   \
+  }
+
+template<> void adjust_hyb_local_jf<XYZ>(int i, int j, int k,
+                                         const int nx, const int ny, const int nz,
+				         field_array_t* RESTRICT f, const grid_t* g) { ADJUSTLOCAL(x,y,z); }
+template<> void adjust_hyb_local_jf<YZX>(int i, int j, int k,
+                                         const int nx, const int ny, const int nz,
+				         field_array_t* RESTRICT f, const grid_t* g) { ADJUSTLOCAL(y,z,x); }
+template<> void adjust_hyb_local_jf<ZXY>(int i, int j, int k,
+                                         const int nx, const int ny, const int nz,
+				         field_array_t* RESTRICT f, const grid_t* g) { ADJUSTLOCAL(z,x,y); }
+
+void
+k_hyb_local_adjust_jf( field_array_t * RESTRICT f,
+                        const grid_t *          g ) {
+    const int nx = g->nx, ny = g->ny, nz = g->nz;
+
+    adjust_hyb_local_jf<XYZ>(-1,0,0,nx,ny,nz,f,g);
+    adjust_hyb_local_jf<YZX>(0,-1,0,nx,ny,nz,f,g);
+    adjust_hyb_local_jf<ZXY>(0,0,-1,nx,ny,nz,f,g);
+    adjust_hyb_local_jf<XYZ>(1,0,0, nx,ny,nz,f,g);
+    adjust_hyb_local_jf<YZX>(0,1,0, nx,ny,nz,f,g);
+    adjust_hyb_local_jf<ZXY>(0,0,1, nx,ny,nz,f,g);
+}
+
+#undef ADJUSTLOCAL