diff --git a/.gitignore b/.gitignore index 40e1f1a..a228424 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,9 @@ build .cproject .project +*.jpg +*.pvtu +*.vtu + +cabanaMD.err +cabanaMD.out diff --git a/input/in.lb b/input/in.lb new file mode 100644 index 0000000..9d84391 --- /dev/null +++ b/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 diff --git a/src/cabanamd.h b/src/cabanamd.h index d5d6222..63be658 100644 --- a/src/cabanamd.h +++ b/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 diff --git a/src/cabanamd_impl.h b/src/cabanamd_impl.h index 8b047ab..b49b951 100644 --- a/src/cabanamd_impl.h +++ b/src/cabanamd_impl.h @@ -244,27 +244,8 @@ void CbnMD::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::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::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::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 vertices; @@ -681,3 +646,37 @@ void CbnMD::check_correctness( int step ) } err.close(); } + +template +void CbnMD::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 ); + } +} diff --git a/src/comm_mpi.h b/src/comm_mpi.h index 13e4156..00db288 100644 --- a/src/comm_mpi.h +++ b/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(); diff --git a/src/comm_mpi_impl.h b/src/comm_mpi_impl.h index d72a1a0..163ea6b 100644 --- a/src/comm_mpi_impl.h +++ b/src/comm_mpi_impl.h @@ -48,7 +48,7 @@ #include -#include +#include template Comm::Comm( t_System *s, T_X_FLOAT comm_depth_ ) @@ -189,7 +189,7 @@ void Comm::reduce_min_float( T_FLOAT *vals, T_INT count ) } template -void Comm::exchange() +T_INT Comm::exchange() { Kokkos::Profiling::pushRegion( "Comm::exchange" ); @@ -270,7 +270,11 @@ void Comm::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 diff --git a/src/inputFile.h b/src/inputFile.h index d70cd54..08ab4c2 100644 --- a/src/inputFile.h +++ b/src/inputFile.h @@ -58,8 +58,9 @@ #include #include -#include #include +#include +#include #include // Class replicating LAMMPS Random velocity initialization with GEOM option @@ -162,8 +163,63 @@ 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 box = { 0, 40, 0, 40, 0, 40 }; + + struct Block + { + double xlo, xhi, ylo, yhi, zlo, zhi; + }; + std::unordered_map regions; + std::unordered_map 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 ®ion_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 get_regions( T_FLOAT xtmp, T_FLOAT ytmp, + T_FLOAT ztmp ) const + { + std::vector 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::max(); + T_X_FLOAT min_y = std::numeric_limits::max(); + T_X_FLOAT min_z = std::numeric_limits::max(); + T_X_FLOAT max_x = std::numeric_limits::min(); + T_X_FLOAT max_y = std::numeric_limits::min(); + T_X_FLOAT max_z = std::numeric_limits::min(); char *data_file; int data_file_type; @@ -171,8 +227,12 @@ class InputFile 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 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 ); diff --git a/src/inputFile_impl.h b/src/inputFile_impl.h index e5803f4..b62b338 100644 --- a/src/inputFile_impl.h +++ b/src/inputFile_impl.h @@ -243,19 +243,22 @@ void InputFile::check_lammps_command( std::string line, if ( words.at( 2 ).compare( "block" ) == 0 ) { known = true; - int box[6]; - box[0] = std::stoi( words.at( 3 ) ); - box[1] = std::stoi( words.at( 4 ) ); - box[2] = std::stoi( words.at( 5 ) ); - box[3] = std::stoi( words.at( 6 ) ); - box[4] = std::stoi( words.at( 7 ) ); - box[5] = std::stoi( words.at( 8 ) ); - if ( ( box[0] != 0 ) || ( box[2] != 0 ) || ( box[4] != 0 ) ) - log_err( err, "LAMMPS-Command: region only allows for boxes " - "with 0,0,0 offset in CabanaMD" ); - lattice_nx = box[1]; - lattice_ny = box[3]; - lattice_nz = box[5]; + auto region_id = words.at( 1 ); + Block block; + block.xlo = std::stod( words.at( 3 ) ); + block.xhi = std::stod( words.at( 4 ) ); + block.ylo = std::stod( words.at( 5 ) ); + block.yhi = std::stod( words.at( 6 ) ); + block.zlo = std::stod( words.at( 7 ) ); + block.zhi = std::stod( words.at( 8 ) ); + regions[region_id] = block; + + min_x = std::min( min_x, lattice_constant * block.xlo ); + min_y = std::min( min_y, lattice_constant * block.ylo ); + min_z = std::min( min_z, lattice_constant * block.zlo ); + max_x = std::max( max_x, lattice_constant * block.xhi ); + max_y = std::max( max_y, lattice_constant * block.yhi ); + max_z = std::max( max_z, lattice_constant * block.zhi ); } else { @@ -270,7 +273,33 @@ void InputFile::check_lammps_command( std::string line, } if ( keyword.compare( "create_atoms" ) == 0 ) { + // supported versions: + // create_atoms TYPE-ID box + // create_atoms TYPE-ID region REGION-ID known = true; + if ( words.at( 2 ) == "region" ) + { + auto type = std::stoi( words.at( 1 ) ); + auto region_id = words.at( 3 ); + if ( regions.find( region_id ) == regions.end() ) + { + log_err( err, "LAMMPS-Command: region '", region_id, + "' is not defined" ); + } + regions_to_type[region_id] = type; + } + else if ( words.at( 2 ) == "box" ) + { + auto type = std::stoi( words.at( 1 ) ); + auto region_id = regions.begin()->first; + regions_to_type[region_id] = type; + } + else + { + log_err( err, + "LAMMPS-Command: 'create_atoms' command only supports " + "'region' option in CabanaMD" ); + } } if ( keyword.compare( "mass" ) == 0 ) { @@ -367,18 +396,17 @@ void InputFile::check_lammps_command( std::string line, if ( keyword.compare( "velocity" ) == 0 ) { known = true; - if ( words.at( 1 ).compare( "all" ) != 0 ) - { - log_err( err, "LAMMPS-Command: 'velocity' command can only be " - "applied to 'all' in CabanaMD" ); - } if ( words.at( 2 ).compare( "create" ) != 0 ) { log_err( err, "LAMMPS-Command: 'velocity' command can only be used " "with option 'create' in CabanaMD" ); } - temperature_target = std::stod( words.at( 3 ) ); - temperature_seed = std::stoi( words.at( 4 ) ); + auto atom_type = + ( words.at( 1 ) == "all" ) ? 1 : std::stoi( words.at( 1 ) ); + auto temperature_target = std::stod( words.at( 3 ) ); + auto temperature_seed = std::stoi( words.at( 4 ) ); + type_to_temperature[atom_type] = { temperature_target, + temperature_seed }; } if ( keyword.compare( "neighbor" ) == 0 ) { @@ -485,6 +513,19 @@ void InputFile::check_lammps_command( std::string line, "commandline --force-iteration" ); } } + if ( keyword.compare( "group" ) == 0 ) + { + // supported version: + // group GROUP-ID region REGION-ID + known = true; + [[maybe_unused]] auto group_id = words.at( 1 ); + if ( words.at( 2 ) != "region" ) + { + log_err( err, "LAMMPS-Command: 'group' command can only be " + "used with 'region' in CabanaMD" ); + } + [[maybe_unused]] auto region_id = words.at( 3 ); + } if ( !known ) { @@ -507,10 +548,7 @@ void InputFile::create_lattice( Comm *comm ) Kokkos::deep_copy( h_mass, s.mass ); // Create the mesh. - T_X_FLOAT max_x = lattice_constant * lattice_nx; - T_X_FLOAT max_y = lattice_constant * lattice_ny; - T_X_FLOAT max_z = lattice_constant * lattice_nz; - std::array global_low = { 0.0, 0.0, 0.0 }; + std::array global_low = { min_x, min_y, min_z }; std::array global_high = { max_x, max_y, max_z }; if ( commandline.vacuum ) { @@ -533,15 +571,15 @@ void InputFile::create_lattice( Comm *comm ) T_INT iy_start = local_mesh_lo_y / lattice_constant - 0.5; T_INT iz_start = local_mesh_lo_z / lattice_constant - 0.5; T_INT ix_end = - std::max( std::min( static_cast( lattice_nx ), + std::max( std::min( max_x / lattice_constant, local_mesh_hi_x / lattice_constant + 0.5 ), static_cast( ix_start ) ); T_INT iy_end = - std::max( std::min( static_cast( lattice_ny ), + std::max( std::min( max_y / lattice_constant, local_mesh_hi_y / lattice_constant + 0.5 ), static_cast( iy_start ) ); T_INT iz_end = - std::max( std::min( static_cast( lattice_nz ), + std::max( std::min( max_z / lattice_constant, local_mesh_hi_z / lattice_constant + 0.5 ), static_cast( iz_start ) ); if ( ix_start == ix_end ) @@ -651,7 +689,6 @@ void InputFile::create_lattice( Comm *comm ) } T_INT n = 0; - for ( T_INT iz = iz_start; iz <= iz_end; iz++ ) { for ( T_INT iy = iy_start; iy <= iy_end; iy++ ) @@ -674,7 +711,10 @@ void InputFile::create_lattice( Comm *comm ) ( ztmp < local_mesh_hi_z ) && ( xtmp < max_x ) && ( ytmp < max_y ) && ( ztmp < max_z ) ) { - n++; + if ( in_any_region( xtmp, ytmp, ztmp ) ) + { + n++; + } } } } @@ -715,12 +755,21 @@ void InputFile::create_lattice( Comm *comm ) ( ztmp < local_mesh_hi_z ) && ( xtmp < max_x ) && ( ytmp < max_y ) && ( ztmp < max_z ) ) { - h_x( n, 0 ) = xtmp; - h_x( n, 1 ) = ytmp; - h_x( n, 2 ) = ztmp; - h_type( n ) = rand() % s.ntypes; - h_id( n ) = n + 1; - n++; + if ( auto region_ids = + get_regions( xtmp, ytmp, ztmp ); + !region_ids.empty() ) + { + h_x( n, 0 ) = xtmp; + h_x( n, 1 ) = ytmp; + h_x( n, 2 ) = ztmp; + + auto rnd = rand() % region_ids.size(); + auto type_id = regions_to_type[region_ids[rnd]]; + + h_type( n ) = type_id - 1; + h_id( n ) = n + 1; + n++; + } } } } @@ -763,7 +812,7 @@ void InputFile::create_lattice( Comm *comm ) { LAMMPS_RandomVelocityGeom random; double x_i[3] = { h_x( i, 0 ), h_x( i, 1 ), h_x( i, 2 ) }; - random.reset( temperature_seed, x_i ); + random.reset( type_to_temperature[h_type( i ) + 1].seed, x_i ); T_FLOAT mass_i = h_mass( h_type( i ) ); T_FLOAT vx = random.uniform() - 0.5; @@ -802,13 +851,16 @@ void InputFile::create_lattice( Comm *comm ) Temperature temp( comm ); T_V_FLOAT T = temp.compute( system ); - T_V_FLOAT T_init_scale = sqrt( temperature_target / T ); - + auto T_init_scale = [=]( int i ) + { + auto target = type_to_temperature[h_type( i ) + 1].temp; + return sqrt( target / T ); + }; for ( T_INT i = 0; i < system->N_local; i++ ) { - h_v( i, 0 ) *= T_init_scale; - h_v( i, 1 ) *= T_init_scale; - h_v( i, 2 ) *= T_init_scale; + h_v( i, 0 ) *= T_init_scale( i ); + h_v( i, 1 ) *= T_init_scale( i ); + h_v( i, 2 ) *= T_init_scale( i ); } system->deep_copy( host_system ); diff --git a/src/vtk_writer.h b/src/vtk_writer.h index 88d54a6..aacaf2f 100644 --- a/src/vtk_writer.h +++ b/src/vtk_writer.h @@ -20,7 +20,7 @@ namespace { -std::string set_width( const int value, const unsigned width = 3 ) +std::string set_width( const int value, const unsigned width = 4 ) { std::ostringstream oss; oss << std::setw( width ) << std::setfill( '0' ) << value; diff --git a/utils/paraview_process_data.py b/utils/paraview_process_data.py new file mode 100644 index 0000000..902c151 --- /dev/null +++ b/utils/paraview_process_data.py @@ -0,0 +1,138 @@ +import os + +# get all the .pvtu files in current directory +# not very pythonic, but it gets the job done +dumps = [] +domain_act = [] +domain_lb = [] +dirname = os.getcwd() +for f in os.listdir(dirname): + if f.endswith(".pvtu"): + if f.startswith("dump"): + dumps.append(os.path.join(dirname, f)) + elif f.startswith("domain_act"): + domain_act.append(os.path.join(dirname, f)) + elif f.startswith("domain_lb"): + domain_lb.append(os.path.join(dirname, f)) + +# trace generated using paraview version 5.12.0 +#import paraview +#paraview.compatibility.major = 5 +#paraview.compatibility.minor = 12 + +#### import the simple module from the paraview +from paraview.simple import * +#### disable automatic camera reset on 'Show' +paraview.simple._DisableFirstRenderCameraReset() + +# create a new 'XML Partitioned Unstructured Grid Reader' +particles = XMLPartitionedUnstructuredGridReader( + registrationName='particles', FileName=dumps) + +# Properties modified on dump_10pvtu +particles.TimeArray = 'None' + +# get active view +renderView1 = GetActiveViewOrCreate('RenderView') + +renderView1.ApplyIsometricView() +renderView1.AxesGrid.Visibility = 1 + +# get display properties +particlesDisplay = GetDisplayProperties(particles, view=renderView1) + +# get animation scene +animationScene1 = GetAnimationScene() + +# update animation scene based on data timesteps +animationScene1.UpdateAnimationUsingDataTimeSteps() + +# create a new 'Glyph' +glyph1 = Glyph(registrationName='Glyph', Input=particles, + GlyphType='Arrow') + +# show data in view +glyph1Display = Show(glyph1, renderView1, 'GeometryRepresentation') + +# trace defaults for the display properties. +glyph1Display.Representation = 'Surface' + +# set scalar coloring +ColorBy(glyph1Display, ('POINTS', 'Velocity')) + +# rescale color and/or opacity maps used to include current data range +glyph1Display.RescaleTransferFunctionToDataRange(True, False) + +# show color bar/color legend +glyph1Display.SetScalarBarVisibility(renderView1, True) + +# get color transfer function/color map for 'Velocity' +velocityLUT = GetColorTransferFunction('Velocity') + +# get opacity transfer function/opacity map for 'Velocity' +velocityPWF = GetOpacityTransferFunction('Velocity') + +# get 2D transfer function for 'Velocity' +velocityTF2D = GetTransferFunction2D('Velocity') + +# Properties modified on glyph1 +glyph1.GlyphMode = 'All Points' +glyph1.GlyphType = 'Sphere' + +glyph1.ScaleArray = ['POINTS', 'No scale array'] +glyph1.ScaleFactor = 1.0 + +# Rescale transfer function +velocityLUT.RescaleTransferFunction(-1.42767, 1.41958) + +# Rescale transfer function +velocityPWF.RescaleTransferFunction(-1.42767, 1.41958) + +# create a new 'XML Partitioned Unstructured Grid Reader' +domain_actual = XMLPartitionedUnstructuredGridReader( + registrationName='domain_actual', FileName=domain_act) + +# Properties modified on domain_act_10pvtu +domain_actual.TimeArray = 'None' + +# show data in view +domain_actualDisplay = Show(domain_actual, renderView1, 'UnstructuredGridRepresentation') + +# trace defaults for the display properties. +domain_actualDisplay.Representation = 'Surface' + +# Rescale transfer function +velocityLUT.RescaleTransferFunction(-1.42767, 1.41958) + +# Rescale transfer function +velocityPWF.RescaleTransferFunction(-1.42767, 1.41958) + +# change representation type +domain_actualDisplay.SetRepresentationType('Wireframe') + +# set scalar coloring +ColorBy(domain_actualDisplay, ('CELLS', 'rank')) + +# rescale color and/or opacity maps used to include current data range +domain_actualDisplay.RescaleTransferFunctionToDataRange(True, False) + +# show color bar/color legend +domain_actualDisplay.SetScalarBarVisibility(renderView1, True) + +# get color transfer function/color map for 'rank' +rankLUT = GetColorTransferFunction('rank') + +# get opacity transfer function/opacity map for 'rank' +rankPWF = GetOpacityTransferFunction('rank') + +# get 2D transfer function for 'rank' +rankTF2D = GetTransferFunction2D('rank') + +# Properties modified on domain_act_10pvtuDisplay +domain_actualDisplay.LineWidth = 2.0 + +# reset view to fit data bounds +renderView1.ResetCamera(False, 0.9) + +# update the view to ensure updated data information +renderView1.Update()