Merge pull request ECP-copa#115 from DARMA-tasking/extended-lammps-su…

…pport

extended `LAMMPS` commands support

.gitignore

@@ -2,3 +2,9 @@
 build
 .cproject
 .project
+*.jpg
+*.pvtu
+*.vtu
+
+cabanaMD.err
+cabanaMD.out

input/in.lb

@@ -0,0 +1,41 @@
+# Example for exercising load balancing (not intended to be physical)
+
+units		lj
+atom_style	atomic
+
+newton off
+
+lattice		fcc 0.8442
+region		lo-box block 0 10 0 10 0 10
+region		mid-box block 5 10 5 10 5 10
+region		hi-box block 0 5 15 20 15 20
+region		hj-box block 15 20 15 20 0 5
+create_box	4 box
+
+mass		1 2.0
+mass		2 8.0
+mass		3 9.0
+mass		4 1.0
+
+create_atoms	1 region lo-box
+create_atoms	2 region mid-box
+create_atoms	3 region hi-box
+create_atoms	4 region hj-box
+
+velocity	1 create 10.4 87287 loop geom
+velocity	2 create 41.6 87287 loop geom
+velocity	3 create 20.8 87287 loop geom
+velocity	4 create 10.4 87287 loop geom
+
+pair_style	lj/cut 2.5
+pair_coeff	1 1 1.0 1.0 2.5
+
+neighbor	0.3 bin
+neigh_modify	every 20 one 50
+comm_modify	cutoff * 20
+fix		1 all nve
+thermo 10
+
+dump dmpvtk all vtk 10 dump%_*.vtu
+
+run		100

src/cabanamd.h

@@ -96,6 +96,11 @@ class CbnMD : public CabanaMD
 
     void dump_binary( int ) override;
     void check_correctness( int ) override;
+
+  private:
+    void print_summary( std::ofstream &out, int step, T_V_FLOAT T, T_F_FLOAT PE,
+                        T_V_FLOAT KE, double time, double rate,
+                        T_INT exchanged = -1 );
 };
 
 #include <cabanamd_impl.h>

src/cabanamd_impl.h

@@ -244,27 +244,8 @@ void CbnMD<t_System, t_Neighbor>::init( InputCL commandline )
         auto T = temp.compute( system );
         auto PE = pote.compute( system, force, neighbor ) / system->N;
         auto KE = kine.compute( system ) / system->N;
-        if ( !_print_lammps )
-        {
-#ifdef CabanaMD_ENABLE_LB
-            log( out, "\n", std::fixed, std::setprecision( 6 ),
-                 "#Timestep Temperature PotE ETot Time Atomsteps/s "
-                 "LBImbalance\n",
-                 step, " ", T, " ", PE, " ", PE + KE, " ",
-                 std::setprecision( 2 ), 0.0, " ", std::scientific, 0.0, " ",
-                 std::setprecision( 2 ), 0.0 );
-#else
-            log( out, "\n", std::fixed, std::setprecision( 6 ),
-                 "#Timestep Temperature PotE ETot Time Atomsteps/s\n", step,
-                 " ", T, " ", PE, " ", PE + KE, " ", std::setprecision( 2 ),
-                 0.0, " ", std::scientific, 0.0 );
-#endif
-        }
-        else
-        {
-            log( out, "\nStep Temp E_pair TotEng CPU\n", step, " ", T, " ", PE,
-                 " ", PE + KE, " ", 0.0 );
-        }
+
+        print_summary( out, step, T, PE, KE, 0.0, 0.0 );
     }
 
     if ( input->dumpbinaryflag )
@@ -313,6 +294,7 @@ void CbnMD<t_System, t_Neighbor>::run()
         integrate_timer.reset();
         integrator->initial_integrate( system );
         integrate_time += integrate_timer.seconds();
+        T_INT exchanged = -1;
 
         if ( step % input->comm_exchange_rate == 0 && step > 0 )
         {
@@ -328,7 +310,7 @@ void CbnMD<t_System, t_Neighbor>::run()
 
             // Exchange atoms across MPI ranks
             comm_timer.reset();
-            comm->exchange();
+            exchanged = comm->exchange();
             comm_time += comm_timer.seconds();
 
             // Sort atoms
@@ -389,31 +371,14 @@ void CbnMD<t_System, t_Neighbor>::run()
             auto T = temp.compute( system );
             auto PE = pote.compute( system, force, neighbor ) / system->N;
             auto KE = kine.compute( system ) / system->N;
+            double time = timer.seconds();
+            double rate =
+                1.0 * system->N * input->thermo_rate / ( time - last_time );
+
+            print_summary( out, step, T, PE, KE, time, rate, exchanged );
+
+            last_time = time;
 
-            if ( !_print_lammps )
-            {
-                double time = timer.seconds();
-                double rate =
-                    1.0 * system->N * input->thermo_rate / ( time - last_time );
-#ifdef CabanaMD_ENABLE_LB
-                log( out, std::fixed, std::setprecision( 6 ), step, " ", T, " ",
-                     PE, " ", PE + KE, " ", std::setprecision( 2 ), time, " ",
-                     std::scientific, rate, " ", std::setprecision( 2 ),
-                     lb->getImbalance() );
-#else
-                log( out, std::fixed, std::setprecision( 6 ), step, " ", T, " ",
-                     PE, " ", PE + KE, " ", std::setprecision( 2 ), time, " ",
-                     std::scientific, rate );
-#endif
-                last_time = time;
-            }
-            else
-            {
-                double time = timer.seconds();
-                log( out, std::fixed, std::setprecision( 6 ), "     ", step,
-                     " ", T, " ", PE, " ", PE + KE, " ", time );
-                last_time = time;
-            }
 #ifdef CabanaMD_ENABLE_LB
             double work = system->N_local + system->N_ghost;
             std::array<double, 6> vertices;
@@ -681,3 +646,37 @@ void CbnMD<t_System, t_Neighbor>::check_correctness( int step )
     }
     err.close();
 }
+
+template <class t_System, class t_Neighbor>
+void CbnMD<t_System, t_Neighbor>::print_summary( std::ofstream &out, int step,
+                                                 T_V_FLOAT T, T_F_FLOAT PE,
+                                                 T_V_FLOAT KE, double time,
+                                                 double rate, T_INT exchanged )
+{
+    if ( !_print_lammps )
+    {
+        if ( step == 0 )
+        {
+            log( out, "\n#Timestep Temperature PotE ETot Time Atomsteps/s "
+#ifdef CabanaMD_ENABLE_LB
+                      "LB-Imbalance MPI-Exchanged"
+#endif
+            );
+        }
+
+        log( out, step, "\t", std::fixed, std::setprecision( 6 ), T, "\t", PE,
+             "\t", PE + KE, "\t", std::setprecision( 2 ), time, "\t",
+             std::scientific, rate
+#ifdef CabanaMD_ENABLE_LB
+             ,
+             "\t", std::setprecision( 2 ), lb->getImbalance(), "\t",
+             ( exchanged != -1 ) ? std::to_string( exchanged ) : "-"
+#endif
+        );
+    }
+    else
+    {
+        log( out, "\nStep Temp E_pair TotEng CPU\n", step, " ", T, " ", PE, " ",
+             PE + KE, " ", time );
+    }
+}

src/comm_mpi.h

@@ -125,7 +125,7 @@ class Comm
     Comm( t_System *s, T_X_FLOAT comm_depth_ );
     void init();
     void create_domain_decomposition();
-    void exchange();
+    T_INT exchange();
     void exchange_halo();
     void update_halo();
     void update_force();

src/comm_mpi_impl.h

@@ -48,7 +48,7 @@
 
 #include <mpi.h>
 
-#include <algorithm>
+#include <output.h>
 
 template <class t_System>
 Comm<t_System>::Comm( t_System *s, T_X_FLOAT comm_depth_ )
@@ -189,7 +189,7 @@ void Comm<t_System>::reduce_min_float( T_FLOAT *vals, T_INT count )
 }
 
 template <class t_System>
-void Comm<t_System>::exchange()
+T_INT Comm<t_System>::exchange()
 {
     Kokkos::Profiling::pushRegion( "Comm::exchange" );
 
@@ -270,7 +270,11 @@ void Comm<t_System>::exchange()
     system->N_local = N_local;
     system->N_ghost = 0;
 
+    T_INT global_send = N_total_send;
+    reduce_int( &global_send, 1 );
+
     Kokkos::Profiling::popRegion();
+    return global_send;
 }
 
 template <class t_System>

src/inputFile.h

@@ -58,8 +58,9 @@
 #include <system.h>
 #include <types.h>
 
-#include <array>
 #include <fstream>
+#include <limits>
+#include <unordered_map>
 #include <vector>
 
 // Class replicating LAMMPS Random velocity initialization with GEOM option
@@ -162,17 +163,76 @@ class InputFile
     int lattice_style = LATTICE_FCC;
     double lattice_constant = 0.8442, lattice_offset_x = 0.0,
            lattice_offset_y = 0.0, lattice_offset_z = 0.0;
-    int lattice_nx, lattice_ny, lattice_nz;
-    std::array<int, 6> box = { 0, 40, 0, 40, 0, 40 };
+
+    struct Block
+    {
+        double xlo, xhi, ylo, yhi, zlo, zhi;
+    };
+    std::unordered_map<std::string, Block> regions;
+    std::unordered_map<std::string, int> regions_to_type;
+
+    bool in_region( T_FLOAT xtmp, T_FLOAT ytmp, T_FLOAT ztmp,
+                    const Block &block ) const
+    {
+        return ( ( xtmp >= lattice_constant * block.xlo ) &&
+                 ( ytmp >= lattice_constant * block.ylo ) &&
+                 ( ztmp >= lattice_constant * block.zlo ) &&
+                 ( xtmp < lattice_constant * block.xhi ) &&
+                 ( ytmp < lattice_constant * block.yhi ) &&
+                 ( ztmp < lattice_constant * block.zhi ) );
+    }
+
+    bool is_empty_region( const std::string &region_id ) const
+    {
+        return regions_to_type.count( region_id ) == 0;
+    }
+
+    bool in_any_region( T_FLOAT xtmp, T_FLOAT ytmp, T_FLOAT ztmp ) const
+    {
+        return std::any_of( regions.cbegin(), regions.cend(),
+                            [=]( auto &pair )
+                            {
+                                return in_region( xtmp, ytmp, ztmp,
+                                                  pair.second ) &&
+                                       !is_empty_region( pair.first );
+                            } );
+    }
+
+    std::vector<std::string> get_regions( T_FLOAT xtmp, T_FLOAT ytmp,
+                                          T_FLOAT ztmp ) const
+    {
+        std::vector<std::string> ret;
+        for ( auto &[region_id, block] : regions )
+        {
+            if ( in_region( xtmp, ytmp, ztmp, block ) &&
+                 !is_empty_region( region_id ) )
+            {
+                ret.emplace_back( region_id );
+            }
+        }
+
+        return ret;
+    }
+
+    T_X_FLOAT min_x = std::numeric_limits<T_X_FLOAT>::max();
+    T_X_FLOAT min_y = std::numeric_limits<T_X_FLOAT>::max();
+    T_X_FLOAT min_z = std::numeric_limits<T_X_FLOAT>::max();
+    T_X_FLOAT max_x = std::numeric_limits<T_X_FLOAT>::min();
+    T_X_FLOAT max_y = std::numeric_limits<T_X_FLOAT>::min();
+    T_X_FLOAT max_z = std::numeric_limits<T_X_FLOAT>::min();
 
     char *data_file;
     int data_file_type;
 
     std::string output_file;
     std::string error_file;
 
-    double temperature_target = 1.4;
-    int temperature_seed = 87287;
+    struct Velocity
+    {
+        double temp;
+        int seed;
+    };
+    std::unordered_map<int, Velocity> type_to_temperature;
 
     int integrator_type = INTEGRATOR_NVE;
     int nsteps = 100;
@@ -205,7 +265,7 @@ class InputFile
     bool read_data_flag = false;
     bool write_data_flag = false;
     bool write_vtk_flag = false;
-    int vtk_rate;
+    int vtk_rate = 0;
     std::string vtk_file;
 
     InputFile( InputCL cl, t_System *s );