diff --git a/CMakeLists.txt b/CMakeLists.txt index 2724a9ff..ec4907a8 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -84,13 +84,20 @@ option(USE_LEGACY_SORT "Enable Legacy Sort Implementation" OFF) option(VPIC_PRINT_MORE_DIGITS "Print more digits in VPIC timer info" OFF) -option(ENABLE_OPENSSL "Enable OpenSSL support for checksums" OFF) - option(DISABLE_DYNAMIC_RESIZING "Prevent particle arrays from dynamically resizing during a run" OFF) +option(USE_HDF5 "Enable HDF5 for use during IO. VPIC does not help you install HDF5" OFF) + +option(USE_OPENPMD "Enable OpenPMD for use during IO. VPIC does not help you install OpenPM" OFF) + +option(OUTPUT_CONVERT_GLOBAL_ID "Convert particle cell id to be global, such that it tells you a unique global location instead of a local offset" ON) + # option to set minimum number of particles set(SET_MIN_NUM_PARTICLES AUTO CACHE STRING "Select minimum number of particles to use, if using dynamic particle array resizing") +# TODO: better name for this? +set(PMD_MAX_IO_CHUNK AUTO CACHE STRING "Select the maxiumum IO write size to use when writing -- applies to particles only, and is specified as number of particles. currently only honored by OpenPMD backend") + #------------------------------------------------------------------------------# # Create include and link aggregates @@ -125,31 +132,22 @@ if(NOT SET_MIN_NUM_PARTICLES STREQUAL "AUTO") add_definitions(-DMIN_NP=${SET_MIN_NUM_PARTICLES}) endif() -#------------------------------------------------------------------------------# -# OpenSSL -#------------------------------------------------------------------------------# - -if(ENABLE_OPENSSL) - find_package(OpenSSL REQUIRED) - - include_directories(${OPENSSL_INCLUDE_DIR}) - string(REPLACE ";" " " string_libraries "${OPENSSL_LIBRARIES}") - set(VPIC_CXX_LIBRARIES "${VPIC_CXX_LIBRARIES} ${string_libraries}") -endif(ENABLE_OPENSSL) +if(NOT PMD_MAX_IO_CHUNK STREQUAL "AUTO") + add_definitions(-DPMD_MAX_IO_CHUNK=${PMD_MAX_IO_CHUNK}) +endif() find_package(Threads REQUIRED) #------------------------------------------------------------------------------# # Act on build options set in project.cmake #------------------------------------------------------------------------------# - #------------------------------------------------------------------------------# # Add options for building with the legacy particle sort implementation. #------------------------------------------------------------------------------# if(USE_LEGACY_SORT) add_definitions(-DVPIC_USE_LEGACY_SORT) - set(VPIC_CXX_FLAGS "${VPIC_CXX_FLAGS} -DVPIC_USE_LEGACY_SORT") + set(VPIC_CXX_FLAGS "${VPIC_CXX_FLAGS} -DVPIC_USE_LEGACY_SORT") endif(USE_LEGACY_SORT) #------------------------------------------------------------------------------# @@ -164,7 +162,7 @@ endif() if(USE_PTHREADS) add_definitions(-DVPIC_USE_PTHREADS) - set(VPIC_CXX_FLAGS "${VPIC_CXX_FLAGS} -DVPIC_USE_PTHREADS") + set(VPIC_CXX_FLAGS "${VPIC_CXX_FLAGS} -DVPIC_USE_PTHREADS") endif(USE_PTHREADS) if(USE_OPENMP) @@ -288,15 +286,23 @@ endif() # Miscellaneous options. #------------------------------------------------------------------------------# -if(ENABLE_OPENSSL) - add_definitions(-DENABLE_OPENSSL) -endif(ENABLE_OPENSSL) - if(VPIC_PRINT_MORE_DIGITS) add_definitions(-DVPIC_PRINT_MORE_DIGITS) set(VPIC_CXX_FLAGS "${VPIC_CXX_FLAGS} -DVPIC_PRINT_MORE_DIGITS") endif(VPIC_PRINT_MORE_DIGITS) +if(USE_HDF5) + # Enable HDF5, and the relevant defines + find_package(HDF5 REQUIRED) + if (NOT HDF5_IS_PARALLEL) + message(FATAL_ERROR "HDF5 Parallel support is required: ${HDF5_IS_PARALLEL}") + endif() + add_definitions(-DVPIC_ENABLE_HDF5) + string(REPLACE ";" " " string_libraries "${HDF5_C_LIBRARIES}") + set(VPIC_CXX_LIBRARIES "${VPIC_CXX_LIBRARIES} ${string_libraries}") + include_directories(${HDF5_INCLUDE_DIRS}) +endif(USE_HDF5) + #------------------------------------------------------------------------------# # Handle vpic compile script last. #------------------------------------------------------------------------------# @@ -312,48 +318,12 @@ if(ENABLE_COVERAGE_BUILD) set(VPIC_CXX_FLAGS "${VPIC_CXX_FLAGS} --coverage") endif(ENABLE_COVERAGE_BUILD) -# process Makefile.run.in to get a simple Makefile.run for a run. Points to -# local built exe wrapper, and has example deck/platform. -configure_file(${CMAKE_SOURCE_DIR}/sample/Makefile.run.in - ${CMAKE_BINARY_DIR}/bin/Makefile.run) - -# Append all defines to VPIC_DEFINES, so it can be seen during input deck building -get_directory_property(ALL_DEFINES DIRECTORY ${CMAKE_SOURCE_DIR} COMPILE_DEFINITIONS) -#string(REPLACE ";" " -D" EEK "${ALL_DEFINES}") -foreach(d ${ALL_DEFINES}) - set(VPIC_DEFINES "${VPIC_DEFINES} -D${d}") -endforeach() - -# install script -configure_file(${CMAKE_SOURCE_DIR}/bin/vpic.in - ${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/vpic-install) -install(FILES ${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/vpic-install - DESTINATION bin - RENAME vpic - PERMISSIONS - OWNER_READ OWNER_WRITE OWNER_EXECUTE - GROUP_READ GROUP_EXECUTE - WORLD_READ WORLD_EXECUTE - ) - install(FILES deck/main.cc deck/wrapper.cc DESTINATION share/vpic) install(FILES deck/wrapper.h DESTINATION include/vpic) install(DIRECTORY src/ DESTINATION include/vpic FILES_MATCHING PATTERN "*.h") -# local script -configure_file(${CMAKE_SOURCE_DIR}/bin/vpic-local.in - ${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/vpic) - -file(COPY ${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/vpic - DESTINATION ${CMAKE_BINARY_DIR}/bin - FILE_PERMISSIONS - OWNER_READ OWNER_WRITE OWNER_EXECUTE - GROUP_READ GROUP_EXECUTE - WORLD_READ WORLD_EXECUTE -) - #------------------------------------------------------------------------------# # Add library target #------------------------------------------------------------------------------# @@ -377,10 +347,49 @@ else() set(VPIC_SRC) install(TARGETS vpic LIBRARY DESTINATION lib ARCHIVE DESTINATION lib) endif() + target_include_directories(vpic INTERFACE ${CMAKE_SOURCE_DIR}/src) -target_link_libraries(vpic ${VPIC_EXPOSE} ${MPI_CXX_LIBRARIES} ${CMAKE_THREAD_LIBS_INIT} ${OPENSSL_LIBRARIES} ${CMAKE_DL_LIBS}) +target_link_libraries(vpic ${VPIC_EXPOSE} ${MPI_CXX_LIBRARIES} ${CMAKE_THREAD_LIBS_INIT} ${CMAKE_DL_LIBS} ${HDF5_C_LIBRARIES}) target_compile_options(vpic ${VPIC_EXPOSE} ${MPI_C_COMPILE_FLAGS}) +# get absolute paths to linked libraries, and their transitive dependencies +function(openpmdreclibs tgtname outname) + get_target_property(PC_PRIVATE_LIBS_TGT ${tgtname} INTERFACE_LINK_LIBRARIES) + foreach(PC_LIB IN LISTS PC_PRIVATE_LIBS_TGT) + if(TARGET ${PC_LIB}) + openpmdreclibs(${PC_LIB} ${outname}) + else() + if(PC_LIB) + string(APPEND ${outname} " ${PC_LIB}") + endif() + endif() + endforeach() + set(${outname} ${${outname}} PARENT_SCOPE) +endfunction() + +if(USE_OPENPMD) + # Enable openPMD, and the relevant defines + find_package(openPMD REQUIRED CONFIG COMPONENTS MPI) + target_link_libraries(vpic PRIVATE openPMD::openPMD) + target_compile_definitions(vpic PRIVATE "-DVPIC_ENABLE_OPENPMD") + + add_definitions(-DVPIC_ENABLE_OPENPMD) + + # legacy stuff for 2-phase compile + get_target_property(openPMD_LIBRARIES openPMD::openPMD LOCATION) + string(REPLACE ";" " " string_libraries "${openPMD_LIBRARIES}") + set(VPIC_CXX_LIBRARIES "${VPIC_CXX_LIBRARIES} ${string_libraries}") + get_target_property(openPMD_TYPE openPMD::openPMD TYPE) + if("${openPMD_TYPE}" STREQUAL "STATIC_LIBRARY") + openpmdreclibs(openPMD openPMD_TRANSITIVE_LIBS) + set(VPIC_CXX_LIBRARIES "${VPIC_CXX_LIBRARIES} ${openPMD_TRANSITIVE_LIBS}") + endif() +endif(USE_OPENPMD) + +if(OUTPUT_CONVERT_GLOBAL_ID) + add_definitions(-DOUTPUT_CONVERT_GLOBAL_ID) +endif(OUTPUT_CONVERT_GLOBAL_ID) + macro(build_a_vpic name deck) if(NOT EXISTS ${deck}) message(FATAL_ERROR "Could not find deck '${deck}'") @@ -438,6 +447,40 @@ if(ENABLE_PERFORMANCE_TESTS) include_directories(${CATCH_DIR}) add_subdirectory(test/performance) endif(ENABLE_PERFORMANCE_TESTS) -#~---------------------------------------------------------------------------~-# -# vim: set tabstop=2 shiftwidth=2 expandtab : -#~---------------------------------------------------------------------------~-# + +# process Makefile.run.in to get a simple Makefile.run for a run. Points to +# local built exe wrapper, and has example deck/platform. +configure_file(${CMAKE_SOURCE_DIR}/sample/Makefile.run.in + ${CMAKE_BINARY_DIR}/bin/Makefile.run) + +# Append all defines to VPIC_DEFINES, so it can be seen during input deck building +get_directory_property(ALL_DEFINES DIRECTORY ${CMAKE_SOURCE_DIR} COMPILE_DEFINITIONS) +#string(REPLACE ";" " -D" EEK "${ALL_DEFINES}") +foreach(d ${ALL_DEFINES}) + set(VPIC_DEFINES "${VPIC_DEFINES} -D${d}") +endforeach() + +# install script +configure_file(${CMAKE_SOURCE_DIR}/bin/vpic.in + ${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/vpic-install) +install(FILES ${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/vpic-install + DESTINATION bin + RENAME vpic + PERMISSIONS + OWNER_READ OWNER_WRITE OWNER_EXECUTE + GROUP_READ GROUP_EXECUTE + WORLD_READ WORLD_EXECUTE +) + +# Configure local script to generate bin/vpic +configure_file(${CMAKE_SOURCE_DIR}/bin/vpic-local.in + ${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/vpic) + +file(COPY ${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/vpic + DESTINATION ${CMAKE_BINARY_DIR}/bin + FILE_PERMISSIONS + OWNER_READ OWNER_WRITE OWNER_EXECUTE + GROUP_READ GROUP_EXECUTE + WORLD_READ WORLD_EXECUTE +) + diff --git a/bin/vpic-local.in b/bin/vpic-local.in index 3b425409..e372c91f 100644 --- a/bin/vpic-local.in +++ b/bin/vpic-local.in @@ -2,6 +2,6 @@ deck=`echo $1 | sed 's,\.cxx,,g;s,\.cc,,g;s,\.cpp,,g;s,.*\/,,g'` -echo "${CMAKE_CXX_COMPILER} ${CMAKE_CXX_FLAGS} -I. -I${CMAKE_SOURCE_DIR}/src ${VPIC_CXX_FLAGS} -DINPUT_DECK=$1 ${CMAKE_SOURCE_DIR}/deck/main.cc ${CMAKE_SOURCE_DIR}/deck/wrapper.cc -o $deck.${CMAKE_SYSTEM_NAME} -Wl,-rpath,${CMAKE_BINARY_DIR} -L${CMAKE_BINARY_DIR} -lvpic ${VPIC_CXX_LIBRARIES} -lpthread -ldl" +echo "${CMAKE_CXX_COMPILER} ${VPIC_DEFINES} ${CMAKE_CXX_FLAGS} -I. -I${CMAKE_SOURCE_DIR}/src ${VPIC_CXX_FLAGS} -DINPUT_DECK=$1 ${CMAKE_SOURCE_DIR}/deck/main.cc ${CMAKE_SOURCE_DIR}/deck/wrapper.cc -o $deck.${CMAKE_SYSTEM_NAME} -Wl,-rpath,${CMAKE_BINARY_DIR} -L${CMAKE_BINARY_DIR} -lvpic ${VPIC_CXX_LIBRARIES} -lpthread -ldl" -${CMAKE_CXX_COMPILER} ${CMAKE_CXX_FLAGS} -I. -I${CMAKE_SOURCE_DIR}/src ${VPIC_CXX_FLAGS} -DINPUT_DECK=$1 ${CMAKE_SOURCE_DIR}/deck/main.cc ${CMAKE_SOURCE_DIR}/deck/wrapper.cc -o $deck.${CMAKE_SYSTEM_NAME} -Wl,-rpath,${CMAKE_BINARY_DIR} -L${CMAKE_BINARY_DIR} -lvpic ${VPIC_CXX_LIBRARIES} -lpthread -ldl +${CMAKE_CXX_COMPILER} ${VPIC_DEFINES} ${CMAKE_CXX_FLAGS} -I. -I${CMAKE_SOURCE_DIR}/src ${VPIC_CXX_FLAGS} -DINPUT_DECK=$1 ${CMAKE_SOURCE_DIR}/deck/main.cc ${CMAKE_SOURCE_DIR}/deck/wrapper.cc -o $deck.${CMAKE_SYSTEM_NAME} -Wl,-rpath,${CMAKE_BINARY_DIR} -L${CMAKE_BINARY_DIR} -lvpic ${VPIC_CXX_LIBRARIES} -lpthread -ldl diff --git a/deck/main.cc b/deck/main.cc index f9f7fb1b..8a9d2352 100644 --- a/deck/main.cc +++ b/deck/main.cc @@ -98,6 +98,12 @@ int main(int argc, char** argv) } simulation = new vpic_simulation(); simulation->initialize( argc, argv ); + + // do post init setup to consume deck values + // which includes setting dump starts steps, as we didn't know it sooner + // TODO: make this use sane functions + simulation->dump_strategy->num_step = simulation->num_step; + REGISTER_OBJECT( &simulation, checkpt_main, restore_main, NULL ); } diff --git a/sample/harrisHDF5 b/sample/harrisHDF5 new file mode 100644 index 00000000..c6c326de --- /dev/null +++ b/sample/harrisHDF5 @@ -0,0 +1,450 @@ +// Magnetic reconnection in a Harris equilibrium thin current sheet +// +// This input deck reproduces the PIC simulations found in: +// William Daughton. "Nonlinear dynamics of thin current sheets." Phys. +// Plasmas. 9(9): 3668-3678. September 2002. +// +// This input deck was written by: +// Kevin J Bowers, Ph.D. +// Plasma Physics Group (X-1) +// Applied Physics Division +// Los Alamos National Lab +// August 2003 - original version +// October 2003 - heavily revised to utilize input deck syntactic sugar +// March/April 2004 - rewritten for domain decomposition V4PIC + +// If you want to use global variables (for example, to store the dump +// intervals for your diagnostics section), it must be done in the globals +// section. Variables declared the globals section will be preserved across +// restart dumps. For example, if the globals section is: +// begin_globals { +// double variable; +// } end_globals +// the double "variable" will be visible to other input deck sections as +// "global->variable". Note: Variables declared in the globals section are set +// to zero before the user's initialization block is executed. Up to 16K +// of global variables can be defined. + + +// Deck only works if VPIC was build with HDF support. Check for that: +#ifndef VPIC_ENABLE_HDF5 +#error "VPIC_ENABLE_HDF5" is required +#endif + +begin_globals { + double energies_interval; + double fields_interval; + double ehydro_interval; + double ihydro_interval; + double eparticle_interval; + double iparticle_interval; + double restart_interval; +}; + +begin_initialization { + // At this point, there is an empty grid and the random number generator is + // seeded with the rank. The grid, materials, species need to be defined. + // Then the initial non-zero fields need to be loaded at time level 0 and the + // particles (position and momentum both) need to be loaded at time level 0. + + + // Example of how to call / set dumping + //field_dump_flag.disableEMAT(); + + + double input_mass_ratio; + int input_seed; + + // Arguments can be passed from the command line to the input deck + if( num_cmdline_arguments!=3 ) { + // Set sensible defaults + input_mass_ratio = 1.0; + input_seed = 0; + + sim_log( "Defaulting to mass_ratio of " << input_mass_ratio << " and seed of " << input_seed ); + sim_log( "For Custom Usage: " << cmdline_argument[0] << " mass_ratio seed" ); + } + else { + input_mass_ratio = atof(cmdline_argument[1]); // Ion mass / electron mass + input_seed = atof(cmdline_argument[2]); // Ion mass / electron mass + sim_log( "Detected input mass_ratio of " << input_mass_ratio << " and seed of " << input_seed ); + } + seed_entropy( input_seed ); + + // Diagnostic messages can be passed written (usually to stderr) + sim_log( "Computing simulation parameters"); + + // Define the system of units for this problem (natural units) + double L = 1; // Length normalization (sheet thickness) + double ec = 1; // Charge normalization + double me = 1; // Mass normalization + double c = 1; // Speed of light + double eps0 = 1; // Permittivity of space + + // Physics parameters + double mi_me = input_mass_ratio; // Ion mass / electron mass + double rhoi_L = 1; // Ion thermal gyroradius / Sheet thickness + double Ti_Te = 1; // Ion temperature / electron temperature + double wpe_wce = 3; // Electron plasma freq / electron cycltron freq + double theta = 0; // Orientation of the simulation wrt current sheet + double taui = 100; // Simulation wci's to run + + // Numerical parameters + double Lx = 16*L; // How big should the box be in the x direction + double Ly = 16*L; // How big should the box be in the y direction + double Lz = 16*L; // How big should the box be in the z direction + double nx = 64; // Global resolution in the x direction + double ny = 64; // Global resolution in the y direction + double nz = 1; // Global resolution in the z direction + double nppc = 64; // Average number of macro particles per cell (both species combined!) + double cfl_req = 0.99; // How close to Courant should we try to run + double wpedt_max = 0.36; // How big a timestep is allowed if Courant is not too restrictive + double damp = 0.001; // Level of radiation damping + + // Derived quantities + double mi = me*mi_me; // Ion mass + double kTe = me*c*c/(2*wpe_wce*wpe_wce*(1+Ti_Te)); // Electron temperature + double kTi = kTe*Ti_Te; // Ion temperature + double vthe = sqrt(2*kTe/me); // Electron thermal velocity (B.D. convention) + double vthi = sqrt(2*kTi/mi); // Ion thermal velocity (B.D. convention) + double wci = vthi/(rhoi_L*L); // Ion cyclotron frequency + double wce = wci*mi_me; // Electron cyclotron frequency + double wpe = wce*wpe_wce; // Electron plasma frequency + double wpi = wpe/sqrt(mi_me); // Ion plasma frequency + double vdre = c*c*wce/(wpe*wpe*L*(1+Ti_Te)); // Electron drift velocity + double vdri = -Ti_Te*vdre; // Ion drift velocity + double b0 = me*wce/ec; // Asymptotic magnetic field strength + double n0 = me*eps0*wpe*wpe/(ec*ec); // Peak electron density (also peak ion density) + double Npe = 2*n0*Ly*Lz*L*tanh(0.5*Lx/L); // Number of physical electrons in box + double Npi = Npe; // Number of physical ions in box + double Ne = 0.5*nppc*nx*ny*nz; // Total macro electrons in box + Ne = trunc_granular(Ne,nproc()); // Make it divisible by number of processors + double Ni = Ne; // Total macro ions in box + double we = Npe/Ne; // Weight of a macro electron + double wi = Npi/Ni; // Weight of a macro ion + double gdri = 1/sqrt(1-vdri*vdri/(c*c)); // gamma of ion drift frame + double gdre = 1/sqrt(1-vdre*vdre/(c*c)); // gamma of electron drift frame + double udri = vdri*gdri; // 4-velocity of ion drift frame + double udre = vdre*gdre; // 4-velocity of electron drift frame + double uthi = sqrt(kTi/mi)/c; // Normalized ion thermal velocity (K.B. convention) + double uthe = sqrt(kTe/me)/c; // Normalized electron thermal velocity (K.B. convention) + double cs = cos(theta); + double sn = sin(theta); + + // Determine the timestep + double dg = courant_length(Lx,Ly,Lz,nx,ny,nz); // Courant length + double dt = cfl_req*dg/c; // Courant limited time step + if( wpe*dt>wpedt_max ) dt=wpedt_max/wpe; // Override time step if plasma frequency limited + + //////////////////////////////////////// + // Setup high level simulation parmeters + + num_step = int(0.2*taui/(wci*dt)); + status_interval = int(1./(wci*dt)); + field_interval = 1; //status_interval; + hydro_interval = 1; //status_interval; + sync_shared_interval = status_interval; + clean_div_e_interval = status_interval; + clean_div_b_interval = status_interval; + + global->energies_interval = status_interval; + global->fields_interval = 1; //status_interval; + global->ehydro_interval = 1; //status_interval; + global->ihydro_interval = 1; //status_interval; + global->eparticle_interval = status_interval; + global->iparticle_interval = status_interval; + global->restart_interval = status_interval; + + /////////////////////////// + // Setup the space and time + + // Setup basic grid parameters + define_units( c, eps0 ); + define_timestep( dt ); + + // Parition a periodic box among the processors sliced uniformly along y + define_periodic_grid( -0.5*Lx, 0, 0, // Low corner + 0.5*Lx, Ly, Lz, // High corner + nx, ny, nz, // Resolution + 1, nproc(), 1 ); // Topology + + // Override some of the boundary conditions to put a particle reflecting + // perfect electrical conductor on the -x and +x boundaries + set_domain_field_bc( BOUNDARY(-1,0,0), pec_fields ); + set_domain_field_bc( BOUNDARY( 1,0,0), pec_fields ); + set_domain_particle_bc( BOUNDARY(-1,0,0), reflect_particles ); + set_domain_particle_bc( BOUNDARY( 1,0,0), reflect_particles ); + + define_material( "vacuum", 1 ); + // Note: define_material defaults to isotropic materials with mu=1,sigma=0 + // Tensor electronic, magnetic and conductive materials are supported + // though. See "shapes" for how to define them and assign them to regions. + // Also, space is initially filled with the first material defined. + + // If you pass NULL to define field array, the standard field array will + // be used (if damp is not provided, no radiation damping will be used). + define_field_array( NULL, damp ); + + //////////////////// + // Setup the species + + // Allow 50% more local_particles in case of non-uniformity + // VPIC will pick the number of movers to use for each species + // Both species use out-of-place sorting + species_t * ion = define_species( "ion", ec, mi, 1.5*Ni/nproc(), -1, 40, 1 ); + species_t * electron = define_species( "electron", -ec, me, 1.5*Ne/nproc(), -1, 20, 1 ); + + /////////////////////////////////////////////////// + // Log diagnostic information about this simulation + + sim_log( "" ); + sim_log( "System of units" ); + sim_log( "L = " << L ); + sim_log( "ec = " << ec ); + sim_log( "me = " << me ); + sim_log( "c = " << c ); + sim_log( "eps0 = " << eps0 ); + sim_log( "" ); + sim_log( "Physics parameters" ); + sim_log( "rhoi/L = " << rhoi_L ); + sim_log( "Ti/Te = " << Ti_Te ); + sim_log( "wpe/wce = " << wpe_wce ); + sim_log( "mi/me = " << mi_me ); + sim_log( "theta = " << theta ); + sim_log( "taui = " << taui ); + sim_log( "" ); + sim_log( "Numerical parameters" ); + sim_log( "num_step = " << num_step ); + sim_log( "dt = " << dt ); + sim_log( "Lx = " << Lx << ", Lx/L = " << Lx/L ); + sim_log( "Ly = " << Ly << ", Ly/L = " << Ly/L ); + sim_log( "Lz = " << Lz << ", Lz/L = " << Lz/L ); + sim_log( "nx = " << nx << ", dx = " << Lx/nx << ", L/dx = " << L*nx/Lx ); + sim_log( "ny = " << ny << ", dy = " << Ly/ny << ", L/dy = " << L*ny/Ly ); + sim_log( "nz = " << nz << ", dz = " << Lz/nz << ", L/dz = " << L*nz/Lz ); + sim_log( "nppc = " << nppc ); + sim_log( "courant = " << c*dt/dg ); + sim_log( "damp = " << damp ); + sim_log( "" ); + sim_log( "Ion parameters" ); + sim_log( "qpi = " << ec << ", mi = " << mi << ", qpi/mi = " << ec/mi ); + sim_log( "vthi = " << vthi << ", vthi/c = " << vthi/c << ", kTi = " << kTi ); + sim_log( "vdri = " << vdri << ", vdri/c = " << vdri/c ); + sim_log( "wpi = " << wpi << ", wpi dt = " << wpi*dt << ", n0 = " << n0 ); + sim_log( "wci = " << wci << ", wci dt = " << wci*dt ); + sim_log( "rhoi = " << vthi/wci << ", L/rhoi = " << L/(vthi/wci) << ", dx/rhoi = " << (Lx/nx)/(vthi/wci) ); + sim_log( "debyei = " << vthi/wpi << ", L/debyei = " << L/(vthi/wpi) << ", dx/debyei = " << (Lx/nx)/(vthi/wpi) ); + sim_log( "Npi = " << Npi << ", Ni = " << Ni << ", Npi/Ni = " << Npi/Ni << ", wi = " << wi ); + sim_log( "" ); + sim_log( "Electron parameters" ); + sim_log( "qpe = " << -ec << ", me = " << me << ", qpe/me = " << -ec/me ); + sim_log( "vthe = " << vthe << ", vthe/c = " << vthe/c << ", kTe = " << kTe ); + sim_log( "vdre = " << vdre << ", vdre/c = " << vdre/c ); + sim_log( "wpe = " << wpe << ", wpe dt = " << wpe*dt << ", n0 = " << n0 ); + sim_log( "wce = " << wce << ", wce dt = " << wce*dt ); + sim_log( "rhoe = " << vthe/wce << ", L/rhoe = " << L/(vthe/wce) << ", dx/rhoe = " << (Lx/nx)/(vthe/wce) ); + sim_log( "debyee = " << vthe/wpe << ", L/debyee = " << L/(vthe/wpe) << ", dx/debyee = " << (Lx/nx)/(vthe/wpe) ); + sim_log( "Npe = " << Npe << ", Ne = " << Ne << ", Npe/Ne = " << Npe/Ne << ", we = " << we ); + sim_log( "" ); + sim_log( "Miscellaneous" ); + sim_log( "nptotal = " << Ni + Ne ); + sim_log( "nproc = " << nproc() ); + sim_log( "" ); + + //////////////////////////// + // Load fields and particles + + sim_log( "Loading fields" ); + + set_region_field( everywhere, 0, 0, 0, // Electric field + 0, -sn*b0*tanh(x/L), cs*b0*tanh(x/L) ); // Magnetic field + // Note: everywhere is a region that encompasses the entire simulation + // In general, regions are specied as logical equations (i.e. x>0 && x+y<2) + + sim_log( "Loading particles" ); + + double ymin = rank()*Ly/nproc(), ymax = (rank()+1)*Ly/nproc(); + + repeat( Ni/nproc() ) { + double x, y, z, ux, uy, uz, d0; + + // Pick an appropriately distributed random location for the pair + do { + x = L*atanh( uniform( rng(0), -1, 1 ) ); + } while( x<=-0.5*Lx || x>=0.5*Lx ); + y = uniform( rng(0), ymin, ymax ); + z = uniform( rng(0), 0, Lz ); + + // For the ion, pick an isothermal normalized momentum in the drift frame + // (this is a proper thermal equilibrium in the non-relativistic limit), + // boost it from the drift frame to the frame with the magnetic field + // along z and then rotate it into the lab frame. Then load the particle. + // Repeat the process for the electron. + + ux = normal( rng(0), 0, uthi ); + uy = normal( rng(0), 0, uthi ); + uz = normal( rng(0), 0, uthi ); + d0 = gdri*uy + sqrt(ux*ux+uy*uy+uz*uz+1)*udri; + uy = d0*cs - uz*sn; + uz = d0*sn + uz*cs; + inject_particle( ion, x, y, z, ux, uy, uz, wi, 0, 0 ); + + ux = normal( rng(0), 0, uthe ); + uy = normal( rng(0), 0, uthe ); + uz = normal( rng(0), 0, uthe ); + d0 = gdre*uy + sqrt(ux*ux+uy*uy+uz*uz+1)*udre; + uy = d0*cs - uz*sn; + uz = d0*sn + uz*cs; + inject_particle( electron, x, y, z, ux, uy, uz, we, 0, 0 ); + } + + // Upon completion of the initialization, the following occurs: + // - The synchronization error (tang E, norm B) is computed between domains + // and tang E / norm B are synchronized by averaging where discrepancies + // are encountered. + // - The initial divergence error of the magnetic field is computed and + // one pass of cleaning is done (for good measure) + // - The bound charge density necessary to give the simulation an initially + // clean divergence e is computed. + // - The particle momentum is uncentered from u_0 to u_{-1/2} + // - The user diagnostics are called on the initial state + // - The physics loop is started + // + // The physics loop consists of: + // - Advance particles from x_0,u_{-1/2} to x_1,u_{1/2} + // - User particle injection at x_{1-age}, u_{1/2} (use inject_particles) + // - User current injection (adjust field(x,y,z).jfx, jfy, jfz) + // - Advance B from B_0 to B_{1/2} + // - Advance E from E_0 to E_1 + // - User field injection to E_1 (adjust field(x,y,z).ex,ey,ez,cbx,cby,cbz) + // - Advance B from B_{1/2} to B_1 + // - (periodically) Divergence clean electric field + // - (periodically) Divergence clean magnetic field + // - (periodically) Synchronize shared tang e and norm b + // - Increment the time step + // - Call user diagnostics + // - (periodically) Print a status message + + // Explicitly enable HDF5 backend for IO dump + // WARNING: Call this after you have set `num_step` (for now.. soon fixed) + + enable_hdf5_dump(); +} + +begin_diagnostics { + +# define should_dump(x) (global->x##_interval>0 && remainder(step(),global->x##_interval)==0) + + if( step()==-10 ) { + // A grid dump contains all grid parameters, field boundary conditions, + // particle boundary conditions and domain connectivity information. This + // is stored in a binary format. Each rank makes a grid dump + dump_grid("grid"); + + // A materials dump contains all the materials parameters. This is in a + // text format. Only rank 0 makes the materials dump + dump_materials("materials"); + + // A species dump contains the physics parameters of a species. This is in + // a text format. Only rank 0 makes the species dump + dump_species("species"); + } + + // Energy dumps store all the energies in various directions of E and B + // and the total kinetic (not including rest mass) energies of each species + // species in a simple text format. By default, the energies are appended to + // the file. However, if a "0" is added to the dump_energies call, a new + // energies dump file will be created. The energies are in the units of the + // problem and are all time centered appropriately. Note: When restarting a + // simulation from a restart dump made at a prior time step to the last + // energies dump, the energies file will have a "hiccup" of intervening + // time levels. This "hiccup" will not occur if the simulation is aborted + // immediately following a restart dump. Energies dumps are in a text + // format and the layout is documented at the top of the file. Only rank 0 + // makes makes an energies dump. + if( should_dump(energies) ) dump_energies( "energies", step()==0 ? 0 : 1 ); + + // Field dumps store the raw electromagnetic fields, sources and material + // placement and a number of auxilliary fields. E, B and RHOB are + // timecentered, JF and TCA are half a step old. Material fields are static + // and the remaining fields (DIV E ERR, DIV B ERR and RHOF) are for + // debugging purposes. By default, field dump filenames are tagged with + // step(). However, if a "0" is added to the call, the filename will not be + // tagged. The JF that gets stored is accumulated with a charge-conserving + // algorithm. As a result, JF is not valid until at least one timestep has + // been completed. Field dumps are in a binary format. Each rank makes a + // field dump. + if( step()==-10 ) dump_fields("fields"); // Get first valid total J + if( should_dump(fields) ) dump_fields("fields"); + + // Hydro dumps store particle charge density, current density and + // stress-energy tensor. All these quantities are known at the time + // t = time(). All these quantities are accumulated trilinear + // node-centered. By default, species dump filenames are tagged with + // step(). However, if a "0" is added to the call, the filename will not + // be tagged. Note that the current density accumulated by this routine is + // purely diagnostic. It is not used by the simulation and it is not + // accumulated using a self-consistent charge-conserving method. Hydro dumps + // are in a binary format. Each rank makes a hydro dump. + if(should_dump(ehydro) ) dump_hydro("electron","ehydro"); + if( should_dump(ihydro) ) dump_hydro("ion", "ihydro"); + + // Particle dumps store the particle data for a given species. The data + // written is known at the time t = time(). By default, particle dumps + // are tagged with step(). However, if a "0" is added to the call, the + // filename will not be tagged. Particle dumps are in a binary format. + // Each rank makes a particle dump. + if( should_dump(eparticle) ) dump_particles("electron","eparticle"); + if( should_dump(iparticle) ) dump_particles("ion", "iparticle"); + + // A checkpt is made by calling checkpt( fbase, tag ) where fname is a string + // and tag is an integer. A typical usage is: + // checkpt( "checkpt", step() ). + // This will cause each process to write their simulation state to a file + // whose name is based on fbase, tag and the node's rank. For the above + // usage, if called on step 314 on a 4 process run, the four files: + // checkpt.314.0, checkpt.314.1, checkpt.314.2, checkpt.314.3 + // to be written. The simulation can then be restarted from this point by + // invoking the application with "--restore checkpt.314". checkpt must be + // the _VERY_ LAST_ diagnostic called. If not, diagnostics performed after + // the checkpt but before the next timestep will be missed on restore. + // Restart dumps are in a binary format unique to the each simulation. + + if( should_dump(restart) ) checkpt( "checkpt", step() ); + + // If you want to write a checkpt after a certain amount of simulation time, + // use uptime() in conjunction with checkpt. For example, this will cause + // the simulation state to be written after 7.5 hours of running to the + // same file every time (useful for dealing with quotas on big machines). + //if( uptime()>=27000 ) { + // checkpt( "timeout", 0 ); + // abort(0); + //} + +# undef should_dump + +} + +begin_particle_injection { + + // No particle injection for this simulation + +} + +begin_current_injection { + + // No current injection for this simulation + +} + +begin_field_injection { + + // No field injection for this simulation + +} + +begin_particle_collisions{ + + // No collisions for this simulation + +} diff --git a/sample/harrisOpenPMD b/sample/harrisOpenPMD new file mode 100644 index 00000000..697e0250 --- /dev/null +++ b/sample/harrisOpenPMD @@ -0,0 +1,455 @@ +// Magnetic reconnection in a Harris equilibrium thin current sheet +// +// This input deck reproduces the PIC simulations found in: +// William Daughton. "Nonlinear dynamics of thin current sheets." Phys. +// Plasmas. 9(9): 3668-3678. September 2002. +// +// This input deck was written by: +// Kevin J Bowers, Ph.D. +// Plasma Physics Group (X-1) +// Applied Physics Division +// Los Alamos National Lab +// August 2003 - original version +// October 2003 - heavily revised to utilize input deck syntactic sugar +// March/April 2004 - rewritten for domain decomposition V4PIC + +// If you want to use global variables (for example, to store the dump +// intervals for your diagnostics section), it must be done in the globals +// section. Variables declared the globals section will be preserved across +// restart dumps. For example, if the globals section is: +// begin_globals { +// double variable; +// } end_globals +// the double "variable" will be visible to other input deck sections as +// "global->variable". Note: Variables declared in the globals section are set +// to zero before the user's initialization block is executed. Up to 16K +// of global variables can be defined. + + +// Deck only works if VPIC was build with HDF support. Check for that: +#ifndef VPIC_ENABLE_OPENPMD +#error "VPIC_ENABLE_OPENPMD" is required +#endif + +begin_globals { + double energies_interval; + double fields_interval; + double ehydro_interval; + double ihydro_interval; + double eparticle_interval; + double iparticle_interval; + double restart_interval; +}; + +begin_initialization { + + enable_openpmd_dump(); + + // TODO: this should be done through a setter once we have a common options + // interface + //dump_strategy->file_type = ".bp"; + + // At this point, there is an empty grid and the random number generator is + // seeded with the rank. The grid, materials, species need to be defined. + // Then the initial non-zero fields need to be loaded at time level 0 and the + // particles (position and momentum both) need to be loaded at time level 0. + + // Example of how to call / set dumping + //field_dump_flag.disableEMAT(); + + double input_mass_ratio; + int input_seed; + + // Arguments can be passed from the command line to the input deck + if( num_cmdline_arguments!=3 ) { + // Set sensible defaults + input_mass_ratio = 1.0; + input_seed = 0; + + sim_log( "Defaulting to mass_ratio of " << input_mass_ratio << " and seed of " << input_seed ); + sim_log( "For Custom Usage: " << cmdline_argument[0] << " mass_ratio seed" ); + } + else { + input_mass_ratio = atof(cmdline_argument[1]); // Ion mass / electron mass + input_seed = atof(cmdline_argument[2]); // Ion mass / electron mass + sim_log( "Detected input mass_ratio of " << input_mass_ratio << " and seed of " << input_seed ); + } + seed_entropy( input_seed ); + + // Diagnostic messages can be passed written (usually to stderr) + sim_log( "Computing simulation parameters"); + + // Define the system of units for this problem (natural units) + double L = 1; // Length normalization (sheet thickness) + double ec = 1; // Charge normalization + double me = 1; // Mass normalization + double c = 1; // Speed of light + double eps0 = 1; // Permittivity of space + + // Physics parameters + double mi_me = input_mass_ratio; // Ion mass / electron mass + double rhoi_L = 1; // Ion thermal gyroradius / Sheet thickness + double Ti_Te = 1; // Ion temperature / electron temperature + double wpe_wce = 3; // Electron plasma freq / electron cycltron freq + double theta = 0; // Orientation of the simulation wrt current sheet + double taui = 100; // Simulation wci's to run + + // Numerical parameters + double Lx = 16*L; // How big should the box be in the x direction + double Ly = 16*L; // How big should the box be in the y direction + double Lz = 16*L; // How big should the box be in the z direction + double nx = 64; // Global resolution in the x direction + double ny = 64; // Global resolution in the y direction + double nz = 1; // Global resolution in the z direction + double nppc = 64; // Average number of macro particles per cell (both species combined!) + double cfl_req = 0.99; // How close to Courant should we try to run + double wpedt_max = 0.36; // How big a timestep is allowed if Courant is not too restrictive + double damp = 0.001; // Level of radiation damping + + // Derived quantities + double mi = me*mi_me; // Ion mass + double kTe = me*c*c/(2*wpe_wce*wpe_wce*(1+Ti_Te)); // Electron temperature + double kTi = kTe*Ti_Te; // Ion temperature + double vthe = sqrt(2*kTe/me); // Electron thermal velocity (B.D. convention) + double vthi = sqrt(2*kTi/mi); // Ion thermal velocity (B.D. convention) + double wci = vthi/(rhoi_L*L); // Ion cyclotron frequency + double wce = wci*mi_me; // Electron cyclotron frequency + double wpe = wce*wpe_wce; // Electron plasma frequency + double wpi = wpe/sqrt(mi_me); // Ion plasma frequency + double vdre = c*c*wce/(wpe*wpe*L*(1+Ti_Te)); // Electron drift velocity + double vdri = -Ti_Te*vdre; // Ion drift velocity + double b0 = me*wce/ec; // Asymptotic magnetic field strength + double n0 = me*eps0*wpe*wpe/(ec*ec); // Peak electron density (also peak ion density) + double Npe = 2*n0*Ly*Lz*L*tanh(0.5*Lx/L); // Number of physical electrons in box + double Npi = Npe; // Number of physical ions in box + double Ne = 0.5*nppc*nx*ny*nz; // Total macro electrons in box + Ne = trunc_granular(Ne,nproc()); // Make it divisible by number of processors + double Ni = Ne; // Total macro ions in box + double we = Npe/Ne; // Weight of a macro electron + double wi = Npi/Ni; // Weight of a macro ion + double gdri = 1/sqrt(1-vdri*vdri/(c*c)); // gamma of ion drift frame + double gdre = 1/sqrt(1-vdre*vdre/(c*c)); // gamma of electron drift frame + double udri = vdri*gdri; // 4-velocity of ion drift frame + double udre = vdre*gdre; // 4-velocity of electron drift frame + double uthi = sqrt(kTi/mi)/c; // Normalized ion thermal velocity (K.B. convention) + double uthe = sqrt(kTe/me)/c; // Normalized electron thermal velocity (K.B. convention) + double cs = cos(theta); + double sn = sin(theta); + + // Determine the timestep + double dg = courant_length(Lx,Ly,Lz,nx,ny,nz); // Courant length + double dt = cfl_req*dg/c; // Courant limited time step + if( wpe*dt>wpedt_max ) dt=wpedt_max/wpe; // Override time step if plasma frequency limited + + //////////////////////////////////////// + // Setup high level simulation parmeters + + num_step = int(0.2*taui/(wci*dt)); + status_interval = int(1./(wci*dt)); + field_interval = status_interval; + hydro_interval = status_interval; + sync_shared_interval = status_interval; + clean_div_e_interval = status_interval; + clean_div_b_interval = status_interval; + + global->energies_interval = status_interval; + global->fields_interval = status_interval; + global->ehydro_interval = status_interval; + global->ihydro_interval = status_interval; + global->eparticle_interval = status_interval; + global->iparticle_interval = status_interval; + global->restart_interval = status_interval; + + /////////////////////////// + // Setup the space and time + + // Setup basic grid parameters + define_units( c, eps0 ); + define_timestep( dt ); + + // Parition a periodic box among the processors sliced uniformly along y + define_periodic_grid( -0.5*Lx, 0, 0, // Low corner + 0.5*Lx, Ly, Lz, // High corner + nx, ny, nz, // Resolution + 1, nproc(), 1 ); // Topology + + // Override some of the boundary conditions to put a particle reflecting + // perfect electrical conductor on the -x and +x boundaries + set_domain_field_bc( BOUNDARY(-1,0,0), pec_fields ); + set_domain_field_bc( BOUNDARY( 1,0,0), pec_fields ); + set_domain_particle_bc( BOUNDARY(-1,0,0), reflect_particles ); + set_domain_particle_bc( BOUNDARY( 1,0,0), reflect_particles ); + + define_material( "vacuum", 1 ); + // Note: define_material defaults to isotropic materials with mu=1,sigma=0 + // Tensor electronic, magnetic and conductive materials are supported + // though. See "shapes" for how to define them and assign them to regions. + // Also, space is initially filled with the first material defined. + + // If you pass NULL to define field array, the standard field array will + // be used (if damp is not provided, no radiation damping will be used). + define_field_array( NULL, damp ); + + //////////////////// + // Setup the species + + // Allow 50% more local_particles in case of non-uniformity + // VPIC will pick the number of movers to use for each species + // Both species use out-of-place sorting + species_t * ion = define_species( "ion", ec, mi, 1.5*Ni/nproc(), -1, 40, 1 ); + species_t * electron = define_species( "electron", -ec, me, 1.5*Ne/nproc(), -1, 20, 1 ); + + /////////////////////////////////////////////////// + // Log diagnostic information about this simulation + + sim_log( "" ); + sim_log( "System of units" ); + sim_log( "L = " << L ); + sim_log( "ec = " << ec ); + sim_log( "me = " << me ); + sim_log( "c = " << c ); + sim_log( "eps0 = " << eps0 ); + sim_log( "" ); + sim_log( "Physics parameters" ); + sim_log( "rhoi/L = " << rhoi_L ); + sim_log( "Ti/Te = " << Ti_Te ); + sim_log( "wpe/wce = " << wpe_wce ); + sim_log( "mi/me = " << mi_me ); + sim_log( "theta = " << theta ); + sim_log( "taui = " << taui ); + sim_log( "" ); + sim_log( "Numerical parameters" ); + sim_log( "num_step = " << num_step ); + sim_log( "dt = " << dt ); + sim_log( "Lx = " << Lx << ", Lx/L = " << Lx/L ); + sim_log( "Ly = " << Ly << ", Ly/L = " << Ly/L ); + sim_log( "Lz = " << Lz << ", Lz/L = " << Lz/L ); + sim_log( "nx = " << nx << ", dx = " << Lx/nx << ", L/dx = " << L*nx/Lx ); + sim_log( "ny = " << ny << ", dy = " << Ly/ny << ", L/dy = " << L*ny/Ly ); + sim_log( "nz = " << nz << ", dz = " << Lz/nz << ", L/dz = " << L*nz/Lz ); + sim_log( "nppc = " << nppc ); + sim_log( "courant = " << c*dt/dg ); + sim_log( "damp = " << damp ); + sim_log( "" ); + sim_log( "Ion parameters" ); + sim_log( "qpi = " << ec << ", mi = " << mi << ", qpi/mi = " << ec/mi ); + sim_log( "vthi = " << vthi << ", vthi/c = " << vthi/c << ", kTi = " << kTi ); + sim_log( "vdri = " << vdri << ", vdri/c = " << vdri/c ); + sim_log( "wpi = " << wpi << ", wpi dt = " << wpi*dt << ", n0 = " << n0 ); + sim_log( "wci = " << wci << ", wci dt = " << wci*dt ); + sim_log( "rhoi = " << vthi/wci << ", L/rhoi = " << L/(vthi/wci) << ", dx/rhoi = " << (Lx/nx)/(vthi/wci) ); + sim_log( "debyei = " << vthi/wpi << ", L/debyei = " << L/(vthi/wpi) << ", dx/debyei = " << (Lx/nx)/(vthi/wpi) ); + sim_log( "Npi = " << Npi << ", Ni = " << Ni << ", Npi/Ni = " << Npi/Ni << ", wi = " << wi ); + sim_log( "" ); + sim_log( "Electron parameters" ); + sim_log( "qpe = " << -ec << ", me = " << me << ", qpe/me = " << -ec/me ); + sim_log( "vthe = " << vthe << ", vthe/c = " << vthe/c << ", kTe = " << kTe ); + sim_log( "vdre = " << vdre << ", vdre/c = " << vdre/c ); + sim_log( "wpe = " << wpe << ", wpe dt = " << wpe*dt << ", n0 = " << n0 ); + sim_log( "wce = " << wce << ", wce dt = " << wce*dt ); + sim_log( "rhoe = " << vthe/wce << ", L/rhoe = " << L/(vthe/wce) << ", dx/rhoe = " << (Lx/nx)/(vthe/wce) ); + sim_log( "debyee = " << vthe/wpe << ", L/debyee = " << L/(vthe/wpe) << ", dx/debyee = " << (Lx/nx)/(vthe/wpe) ); + sim_log( "Npe = " << Npe << ", Ne = " << Ne << ", Npe/Ne = " << Npe/Ne << ", we = " << we ); + sim_log( "" ); + sim_log( "Miscellaneous" ); + sim_log( "nptotal = " << Ni + Ne ); + sim_log( "nproc = " << nproc() ); + sim_log( "" ); + + //////////////////////////// + // Load fields and particles + + sim_log( "Loading fields" ); + + set_region_field( everywhere, 0, 0, 0, // Electric field + 0, -sn*b0*tanh(x/L), cs*b0*tanh(x/L) ); // Magnetic field + // Note: everywhere is a region that encompasses the entire simulation + // In general, regions are specied as logical equations (i.e. x>0 && x+y<2) + + sim_log( "Loading particles" ); + + double ymin = rank()*Ly/nproc(), ymax = (rank()+1)*Ly/nproc(); + + repeat( Ni/nproc() ) { + double x, y, z, ux, uy, uz, d0; + + // Pick an appropriately distributed random location for the pair + do { + x = L*atanh( uniform( rng(0), -1, 1 ) ); + } while( x<=-0.5*Lx || x>=0.5*Lx ); + y = uniform( rng(0), ymin, ymax ); + z = uniform( rng(0), 0, Lz ); + + // For the ion, pick an isothermal normalized momentum in the drift frame + // (this is a proper thermal equilibrium in the non-relativistic limit), + // boost it from the drift frame to the frame with the magnetic field + // along z and then rotate it into the lab frame. Then load the particle. + // Repeat the process for the electron. + + ux = normal( rng(0), 0, uthi ); + uy = normal( rng(0), 0, uthi ); + uz = normal( rng(0), 0, uthi ); + d0 = gdri*uy + sqrt(ux*ux+uy*uy+uz*uz+1)*udri; + uy = d0*cs - uz*sn; + uz = d0*sn + uz*cs; + inject_particle( ion, x, y, z, ux, uy, uz, wi, 0, 0 ); + + ux = normal( rng(0), 0, uthe ); + uy = normal( rng(0), 0, uthe ); + uz = normal( rng(0), 0, uthe ); + d0 = gdre*uy + sqrt(ux*ux+uy*uy+uz*uz+1)*udre; + uy = d0*cs - uz*sn; + uz = d0*sn + uz*cs; + inject_particle( electron, x, y, z, ux, uy, uz, we, 0, 0 ); + } + + // Upon completion of the initialization, the following occurs: + // - The synchronization error (tang E, norm B) is computed between domains + // and tang E / norm B are synchronized by averaging where discrepancies + // are encountered. + // - The initial divergence error of the magnetic field is computed and + // one pass of cleaning is done (for good measure) + // - The bound charge density necessary to give the simulation an initially + // clean divergence e is computed. + // - The particle momentum is uncentered from u_0 to u_{-1/2} + // - The user diagnostics are called on the initial state + // - The physics loop is started + // + // The physics loop consists of: + // - Advance particles from x_0,u_{-1/2} to x_1,u_{1/2} + // - User particle injection at x_{1-age}, u_{1/2} (use inject_particles) + // - User current injection (adjust field(x,y,z).jfx, jfy, jfz) + // - Advance B from B_0 to B_{1/2} + // - Advance E from E_0 to E_1 + // - User field injection to E_1 (adjust field(x,y,z).ex,ey,ez,cbx,cby,cbz) + // - Advance B from B_{1/2} to B_1 + // - (periodically) Divergence clean electric field + // - (periodically) Divergence clean magnetic field + // - (periodically) Synchronize shared tang e and norm b + // - Increment the time step + // - Call user diagnostics + // - (periodically) Print a status message +} + +begin_diagnostics { + +# define should_dump(x) (global->x##_interval>0 && remainder(step(),global->x##_interval)==0) + + if( step()==-10 ) { + // A grid dump contains all grid parameters, field boundary conditions, + // particle boundary conditions and domain connectivity information. This + // is stored in a binary format. Each rank makes a grid dump + dump_grid("grid"); + + // A materials dump contains all the materials parameters. This is in a + // text format. Only rank 0 makes the materials dump + dump_materials("materials"); + + // A species dump contains the physics parameters of a species. This is in + // a text format. Only rank 0 makes the species dump + dump_species("species"); + } + + // Energy dumps store all the energies in various directions of E and B + // and the total kinetic (not including rest mass) energies of each species + // species in a simple text format. By default, the energies are appended to + // the file. However, if a "0" is added to the dump_energies call, a new + // energies dump file will be created. The energies are in the units of the + // problem and are all time centered appropriately. Note: When restarting a + // simulation from a restart dump made at a prior time step to the last + // energies dump, the energies file will have a "hiccup" of intervening + // time levels. This "hiccup" will not occur if the simulation is aborted + // immediately following a restart dump. Energies dumps are in a text + // format and the layout is documented at the top of the file. Only rank 0 + // makes makes an energies dump. + if( should_dump(energies) ) dump_energies( "energies", step()==0 ? 0 : 1 ); + + // Field dumps store the raw electromagnetic fields, sources and material + // placement and a number of auxilliary fields. E, B and RHOB are + // timecentered, JF and TCA are half a step old. Material fields are static + // and the remaining fields (DIV E ERR, DIV B ERR and RHOF) are for + // debugging purposes. By default, field dump filenames are tagged with + // step(). However, if a "0" is added to the call, the filename will not be + // tagged. The JF that gets stored is accumulated with a charge-conserving + // algorithm. As a result, JF is not valid until at least one timestep has + // been completed. Field dumps are in a binary format. Each rank makes a + // field dump. + + // TODO: passing in the field extension as part of the name doesn't work for + // the other functions, as they use it to look up species + std::string openpm_field_name = "fields.h5"; + //std::string openpm_field_name = "fields.bp"; + if( step()==-10 ) dump_fields(openpm_field_name.c_str()); // Get first valid total J + if( should_dump(fields) ) dump_fields(openpm_field_name.c_str()); + + // Hydro dumps store particle charge density, current density and + // stress-energy tensor. All these quantities are known at the time + // t = time(). All these quantities are accumulated trilinear + // node-centered. By default, species dump filenames are tagged with + // step(). However, if a "0" is added to the call, the filename will not + // be tagged. Note that the current density accumulated by this routine is + // purely diagnostic. It is not used by the simulation and it is not + // accumulated using a self-consistent charge-conserving method. Hydro dumps + // are in a binary format. Each rank makes a hydro dump. + if(should_dump(ehydro) ) dump_hydro("electron","ehydro"); + //if( should_dump(ihydro) ) dump_hydro_hdf5("ion", "ihydro"); + + // Particle dumps store the particle data for a given species. The data + // written is known at the time t = time(). By default, particle dumps + // are tagged with step(). However, if a "0" is added to the call, the + // filename will not be tagged. Particle dumps are in a binary format. + // Each rank makes a particle dump. + if( should_dump(eparticle) ) dump_particles("electron","eparticle"); + //if( should_dump(iparticle) ) dump_particles_hdf5("ion", "iparticle"); + + // A checkpt is made by calling checkpt( fbase, tag ) where fname is a string + // and tag is an integer. A typical usage is: + // checkpt( "checkpt", step() ). + // This will cause each process to write their simulation state to a file + // whose name is based on fbase, tag and the node's rank. For the above + // usage, if called on step 314 on a 4 process run, the four files: + // checkpt.314.0, checkpt.314.1, checkpt.314.2, checkpt.314.3 + // to be written. The simulation can then be restarted from this point by + // invoking the application with "--restore checkpt.314". checkpt must be + // the _VERY_ LAST_ diagnostic called. If not, diagnostics performed after + // the checkpt but before the next timestep will be missed on restore. + // Restart dumps are in a binary format unique to the each simulation. + + if( should_dump(restart) ) checkpt( "checkpt", step() ); + + // If you want to write a checkpt after a certain amount of simulation time, + // use uptime() in conjunction with checkpt. For example, this will cause + // the simulation state to be written after 7.5 hours of running to the + // same file every time (useful for dealing with quotas on big machines). + //if( uptime()>=27000 ) { + // checkpt( "timeout", 0 ); + // abort(0); + //} + +# undef should_dump + +} + +begin_particle_injection { + + // No particle injection for this simulation + +} + +begin_current_injection { + + // No current injection for this simulation + +} + +begin_field_injection { + + // No field injection for this simulation + +} + +begin_particle_collisions{ + + // No collisions for this simulation + +} diff --git a/sample/read_openpmd.py b/sample/read_openpmd.py new file mode 100644 index 00000000..84b19009 --- /dev/null +++ b/sample/read_openpmd.py @@ -0,0 +1,39 @@ + +import openpmd_api as api + +# example: data handling +import numpy as np + +file_name = "./fields.h5" +series = api.Series( file_name, api.Access_Type.read_only) + +print(list(series.iterations)) + +from pprint import pprint +#pprint(vars(series)) +#pprint(vars(series.iterations)) + +i = series.iterations[1]; + +print("openPMD version: ", + series.openPMD) + +# record +cB = i.meshes["B"] + +# record components +cbx = cB["x"] + +x_data = cbx.load_chunk() + +series.flush() + +extent = cbx.shape + +print( + "First values in E_x " + "of shape: ", + extent) + + +print(x_data) diff --git a/src/field_advance/field_advance.h b/src/field_advance/field_advance.h index d1cee710..0d435e8a 100644 --- a/src/field_advance/field_advance.h +++ b/src/field_advance/field_advance.h @@ -13,7 +13,7 @@ // // This module implements the following the difference equations on a // superhexahedral domain decomposed Yee-mesh: -// +// // advance_b -> Finite Differenced Faraday // cB_new = cB_old - frac c dt curl E // @@ -32,7 +32,7 @@ // rapidly reduce RMS divergence error assuming divergences errors // are due to accumulation of numerical roundoff when integrating // Faraday. See clean_div.c for details. -// +// // div_clean_e -> Modified Marder pass on electric fields // E_new = E_old + drive D dt grad err_mul div ( epsr E_old - rho/eps0 ) // Since the total rho may not be known everywhere (for example in @@ -65,7 +65,7 @@ // fmatx,fmaty,fmatz are all on the "face // mesh". rhof,rhob,div_e_err,nmat are on the "nodes mesh". // div_b_err,cmat are on the "cell mesh". -// +// // Above, for "edge mesh" quantities, interior means that the // component is not a tangential field directly on the surface of the // domain. For "face mesh" quantities, interior means that the @@ -97,7 +97,7 @@ // ... // material_coefficients = new_material_coefficients(grid,material_list); // fields = new_fields(grid); -// +// // ... Set the initial field values and place materials ... // // synchronize_fields(fields,grid); @@ -107,7 +107,7 @@ // initial fields or errors in the source terms or different floating // point properties on different nodes cause the shared faces to have // different fields). -// +// // To advance the fields in a PIC simulation with TCA radation damping // and periodic divergence cleaning, the following sequence is // suggested: @@ -118,7 +118,7 @@ // if( should_clean_div_e ) { // ... adjust rho_f, rho_b and/or rho_c as necessary // do { -// rms_err = clean_div_e( fields, material_coefficients, grid ); +// rms_err = clean_div_e( fields, material_coefficients, grid ); // } while( rms_err_too_high ); // } // if( should_clean_div_b ) { diff --git a/src/grid/grid.h b/src/grid/grid.h index 7654fe94..b95a034e 100644 --- a/src/grid/grid.h +++ b/src/grid/grid.h @@ -1,4 +1,4 @@ -/* +/* * Written by: * Kevin J. Bowers, Ph.D. * Plasma Physics Group (X-1) @@ -46,7 +46,7 @@ enum grid_enums { // B_tang -> Symmetric | B_tang -> Anti-symmetric // E_norm -> Symmetric | E_norm -> Anti-symmetric (see note) // div B -> Symmetric | div B -> Anti-symmetric - // + // // Note: B_norm is tricky. For a symmetry plane, B_norm on the // boundary must be zero as there are no magnetic charges (a // non-zero B_norm would imply an infinitesimal layer of magnetic @@ -80,7 +80,7 @@ typedef struct grid { int64_t step; // Current timestep double t0; // Simulation time corresponding to step 0 - // Phase 2 grid data structures + // Phase 2 grid data structures float x0, y0, z0; // Min corner local domain (must be coherent) float x1, y1, z1; // Max corner local domain (must be coherent) int nx, ny, nz; // Local voxel mesh resolution. Voxels are @@ -99,6 +99,9 @@ typedef struct grid { // 0 ... nproc-1 ... comm boundary condition // <0 ... locally applied boundary condition + int gpx, gpy, gpz = -1; // Store global processor decomposition to let us figure + // out where we are in the global decomposition + // Phase 3 grid data structures // NOTE: VOXEL INDEXING LIMITS NUMBER OF VOXELS TO 2^31 (INCLUDING // GHOSTS) PER NODE. NEIGHBOR INDEXING FURTHER LIMITS TO @@ -135,6 +138,20 @@ typedef struct grid { #define VOXEL(x,y,z, nx,ny,nz) ((x) + ((nx)+2)*((y) + ((ny)+2)*(z))) +// TODO: make the asymmetry in how nx+2 is handled more obvious +#define UNVOXEL(rank, ix, iy, iz, nx, ny, nz) BEGIN_PRIMITIVE { \ + int _ix, _iy, _iz; \ + _ix = (rank); /* ix = ix+gpx*( iy+gpy*iz ) */ \ + _iy = _ix/int(nx); /* iy = iy+gpy*iz */ \ + _ix -= _iy*int(nx); /* ix = ix */ \ + _iz = _iy/int(ny); /* iz = iz */ \ + _iy -= _iz*int(ny); /* iy = iy */ \ + (ix) = _ix; \ + (iy) = _iy; \ + (iz) = _iz; \ + } END_PRIMITIVE + + // Advance the voxel mesh index (v) and corresponding voxel mesh // coordinates (x,y,z) in a region with min- and max-corners of // (xl,yl,zl) and (xh,yh,zh) of a (nx,ny,nz) resolution voxel mesh in @@ -147,7 +164,7 @@ typedef struct grid { // inner loops.) // // This is written with seeming extraneously if tests in order to get -// the compiler to generate branceless conditional move and add +// the compiler to generate branceless conditional move and add // instructions (none of the branches below are actual branches in // assembly). @@ -311,7 +328,7 @@ end_send_port( int i, // x port coord ([-1,0,1]) // ordering (e.g. inner loop increments x-index). // // jobs are indexed from 0 to n_job-1. jobs are _always_ have the -// number of voxels an integer multiple of the bundle size. If job +// number of voxels an integer multiple of the bundle size. If job // is set to n_job, this function will determine the parameters of // the final incomplete bundle. diff --git a/src/grid/partition.cc b/src/grid/partition.cc index 96664b78..ff9b09f4 100644 --- a/src/grid/partition.cc +++ b/src/grid/partition.cc @@ -1,4 +1,4 @@ -/* +/* * Written by: * Kevin J. Bowers, Ph.D. * Plasma Physics Group (X-1) @@ -39,7 +39,7 @@ partition_periodic_box( grid_t * g, int gnx, int gny, int gnz, int gpx, int gpy, int gpz ) { double f; - int rank, px, py, pz; + int rank, px, py, pz; // Make sure the grid can be setup @@ -55,6 +55,11 @@ partition_periodic_box( grid_t * g, // Setup basic variables RANK_TO_INDEX( world_rank, px,py,pz ); + // Capture global processor decomposition + g->gpx = gpx; + g->gpy = gpy; + g->gpz = gpz; + g->dx = (gx1-gx0)/(double)gnx; g->dy = (gy1-gy0)/(double)gny; g->dz = (gz1-gz0)/(double)gnz; @@ -96,7 +101,7 @@ partition_absorbing_box( grid_t * g, int gnx, int gny, int gnz, int gpx, int gpy, int gpz, int pbc ) { - int px, py, pz; + int px, py, pz; partition_periodic_box( g, gx0, gy0, gz0, @@ -108,30 +113,30 @@ partition_absorbing_box( grid_t * g, RANK_TO_INDEX( world_rank, px,py,pz ); - if( px==0 && gnx>1 ) { + if( px==0 && gnx>1 ) { set_fbc(g,BOUNDARY(-1,0,0),absorb_fields); set_pbc(g,BOUNDARY(-1,0,0),pbc); - } + } if( px==gpx-1 && gnx>1 ) { set_fbc(g,BOUNDARY( 1,0,0),absorb_fields); set_pbc(g,BOUNDARY( 1,0,0),pbc); } - if( py==0 && gny>1 ) { + if( py==0 && gny>1 ) { set_fbc(g,BOUNDARY(0,-1,0),absorb_fields); set_pbc(g,BOUNDARY(0,-1,0),pbc); - } + } if( py==gpy-1 && gny>1 ) { set_fbc(g,BOUNDARY(0, 1,0),absorb_fields); set_pbc(g,BOUNDARY(0, 1,0),pbc); } - if( pz==0 && gnz>1 ) { + if( pz==0 && gnz>1 ) { set_fbc(g,BOUNDARY(0,0,-1),absorb_fields); set_pbc(g,BOUNDARY(0,0,-1),pbc); - } + } if( pz==gpz-1 && gnz>1 ) { set_fbc(g,BOUNDARY(0,0, 1),absorb_fields); @@ -148,7 +153,7 @@ partition_metal_box( grid_t * g, double gx1, double gy1, double gz1, int gnx, int gny, int gnz, int gpx, int gpy, int gpz ) { - int px, py, pz; + int px, py, pz; partition_periodic_box( g, gx0, gy0, gz0, diff --git a/src/species_advance/species_advance.cc b/src/species_advance/species_advance.cc index 0e85a646..2ed53cbb 100644 --- a/src/species_advance/species_advance.cc +++ b/src/species_advance/species_advance.cc @@ -1,4 +1,4 @@ -/* +/* * Written by: * Kevin J. Bowers, Ph.D. * Plasma Physics Group (X-1) @@ -146,7 +146,7 @@ species( const char * name, sp->sort_out_of_place = sort_out_of_place; MALLOC_ALIGNED( sp->partition, g->nv+1, 128 ); - sp->g = g; + sp->g = g; /* id, next are set by append species */ diff --git a/src/species_advance/species_advance_aos.h b/src/species_advance/species_advance_aos.h index 3e1af9ad..47fa1a78 100644 --- a/src/species_advance/species_advance_aos.h +++ b/src/species_advance/species_advance_aos.h @@ -12,6 +12,8 @@ #ifndef _species_advance_aos_h_ #define _species_advance_aos_h_ +// TODO: should we restrict the direct include of this header? + typedef int32_t species_id; // Must be 32-bit wide for particle_injector_t // FIXME: Eventually particle_t (definitely) and their other formats diff --git a/src/util/io/FileIO.h b/src/util/io/FileIO.h index 0d8ed6da..74221451 100644 --- a/src/util/io/FileIO.h +++ b/src/util/io/FileIO.h @@ -13,6 +13,7 @@ #define FileIO_h #include +#include #include "FileIOData.h" /*! diff --git a/src/util/util_base.h b/src/util/util_base.h index bc9db329..4f2ada6c 100644 --- a/src/util/util_base.h +++ b/src/util/util_base.h @@ -1,4 +1,4 @@ -/* +/* * Written by: * Kevin J. Bowers, Ph.D. * Plasma Physics Group (X-1) @@ -21,7 +21,7 @@ #endif // C99 does requires some key macros of stdint to only be defined in -// C++ implementations if explicitly requested. +// C++ implementations if explicitly requested. #define __STDC_LIMIT_MACROS @@ -102,7 +102,7 @@ typedef struct collective collective_t; #ifndef RESTRICT #define RESTRICT __restrict -#endif +#endif // Normal pointers (e.g. a *) are in whatever address space the given // compile unit uses. However, sometimes it is necessary to declare @@ -154,7 +154,7 @@ typedef struct collective collective_t; // allow correct autogeneration when no alignment necessary ... sigh // ... -#define PAD(s,a) ( (a) - ( (s) & ( (a)-1 ) ) ) +#define PAD(s,a) ( (a) - ( (s) & ( (a)-1 ) ) ) // POW2_CEIL rounds "u" up to the nearest multiple of the power of two // "a". If u is a multiple of "a", its value is unchanged. "a" should @@ -344,7 +344,7 @@ void detect_old_style_arguments(int* pargc, char *** pargv); #define MALLOC(x,n) \ util_malloc( "MALLOC( "#x", "#n" (%lu bytes) ) at " \ __FILE__ "(" EXPAND_AND_STRINGIFY(__LINE__) ") failed", \ - &(x), (n)*sizeof(*(x)) ) + &(x), (n)*sizeof(*(x)) ) void util_malloc( const char * err_fmt, // Has exactly one %lu in it @@ -370,7 +370,7 @@ util_free( void * mem_ref ); #n" (%lu bytes), " \ #a" (%lu bytes) ) at " \ __FILE__ "(" EXPAND_AND_STRINGIFY(__LINE__) ") failed", \ - &(x), (n)*sizeof(*(x)), (a) ) + &(x), (n)*sizeof(*(x)), (a) ) void diff --git a/src/vpic/dump.cc b/src/vpic/dump.cc index 62505147..21804882 100644 --- a/src/vpic/dump.cc +++ b/src/vpic/dump.cc @@ -10,9 +10,10 @@ */ #include +#include -#include "vpic.h" #include "dumpmacros.h" +#include "vpic.h" #include "../util/io/FileUtils.h" /* -1 means no ranks talk */ @@ -22,6 +23,36 @@ // COMPATIBLE WITH EXISTING EXTERNAL 3RD PARTY VISUALIZATION SOFTWARE. // IN THE LONG RUN, THIS EXTERNAL SOFTWARE WILL NEED TO BE UPDATED. +std::array global_particle_index(int local_i, grid_t* grid, int rank) +{ + int ix, iy, iz, rx, ry, rz; + // Convert rank to local x/y/z + UNVOXEL(rank, rx, ry, rz, grid->gpx, grid->gpy, grid->gpz); + + // Calculate local ix/iy/iz + UNVOXEL(local_i, ix, iy, iz, grid->nx+2, grid->ny+2, grid->nz+2); + + // Account for the "first" ghost cell + ix = ix - 1; + iy = iy - 1; + iz = iz - 1; + + // Convert ix/iy/iz to global + int gix = ix + (grid->nx * (rx)); + int giy = iy + (grid->ny * (ry)); + int giz = iz + (grid->nz * (rz)); + + // calculate global grid sizes + int gnx = grid->nx * grid->gpx; + int gny = grid->ny * grid->gpy; + int gnz = grid->nz * grid->gpz; + + // TODO: find a better way to account for the hard coded ghosts in VOXEL + int global_i = VOXEL(gix, giy, giz, gnx-2, gny-2, gnz-2); + + return { global_i, gix, giy, giz }; +} + int vpic_simulation::dump_mkdir(const char * dname) { return FileUtils::makeDirectory(dname); } // dump_mkdir @@ -34,6 +65,64 @@ int vpic_simulation::dump_cwd(char * dname, size_t size) { * ASCII dump IO *****************************************************************************/ +void vpic_simulation::enable_binary_dump() { + dump_strategy = std::unique_ptr(new BinaryDump( rank(), nproc() )); +} + +#ifdef VPIC_ENABLE_HDF5 +void vpic_simulation::enable_hdf5_dump() { + std::cout << "Enabling HDF5 IO backend" << std::endl; + dump_strategy = std::unique_ptr(new HDF5Dump( rank(), nproc() )); + dump_strategy->num_step = num_step; +} +#endif + +#ifdef VPIC_ENABLE_OPENPMD +void vpic_simulation::enable_openpmd_dump() { + std::cout << "Enabling openPMD IO backend" << std::endl; + dump_strategy = std::unique_ptr(new OpenPMDDump( rank(), nproc() )); +} +#endif + +void vpic_simulation::dump_particles( const char *sp_name, + const char *fbase, + int ftag ) +{ + species_t* sp = find_species_name(sp_name, species_list); + dump_strategy->dump_particles( + fbase, + sp, + grid, + step(), + interpolator_array, + ftag + ); +} + +void vpic_simulation::dump_fields( const char *fbase, int ftag ) +{ + dump_strategy->dump_fields( + fbase, + step(),grid, + field_array, + ftag + ); +} + +void vpic_simulation::dump_hydro( const char *sp_name, const char *fbase, int ftag ) +{ + species_t * sp = find_species_name(sp_name, species_list); + dump_strategy->dump_hydro( + fbase, + step(), + hydro_array, + sp, + interpolator_array, + grid, + ftag + ); +} + void vpic_simulation::dump_energies( const char *fname, int append ) { @@ -115,26 +204,6 @@ vpic_simulation::dump_materials( const char *fname ) { * Binary dump IO *****************************************************************************/ -/* -enum dump_types { - grid_dump = 0, - field_dump = 1, - hydro_dump = 2, - particle_dump = 3, - restart_dump = 4 -}; -*/ - -// TODO: should this be an enum? -namespace dump_type { - const int grid_dump = 0; - const int field_dump = 1; - const int hydro_dump = 2; - const int particle_dump = 3; - const int restart_dump = 4; - const int history_dump = 5; -} // namespace - void vpic_simulation::dump_grid( const char *fbase ) { char fname[256]; @@ -149,14 +218,14 @@ vpic_simulation::dump_grid( const char *fbase ) { if( status==fail ) ERROR(( "Could not open \"%s\".", fname )); /* IMPORTANT: these values are written in WRITE_HEADER_V0 */ - nxout = grid->nx; - nyout = grid->ny; - nzout = grid->nz; - dxout = grid->dx; - dyout = grid->dy; - dzout = grid->dz; + size_t nxout = grid->nx; + size_t nyout = grid->ny; + size_t nzout = grid->nz; + float dxout = grid->dx; + float dyout = grid->dy; + float dzout = grid->dz; - WRITE_HEADER_V0( dump_type::grid_dump, -1, 0, fileIO ); + WRITE_HEADER_V0( dump_type::grid_dump, -1, 0, fileIO, step(), rank(), nproc()); dim[0] = 3; dim[1] = 3; @@ -178,148 +247,6 @@ vpic_simulation::dump_grid( const char *fbase ) { if( fileIO.close() ) ERROR(( "File close failed on dump grid!!!" )); } -void -vpic_simulation::dump_fields( const char *fbase, int ftag ) { - char fname[256]; - FileIO fileIO; - int dim[3]; - - if( !fbase ) ERROR(( "Invalid filename" )); - - if( rank()==0 ) MESSAGE(( "Dumping fields to \"%s\"", fbase )); - - if( ftag ) sprintf( fname, "%s.%li.%i", fbase, (long)step(), rank() ); - else sprintf( fname, "%s.%i", fbase, rank() ); - - FileIOStatus status = fileIO.open(fname, io_write); - if( status==fail ) ERROR(( "Could not open \"%s\".", fname )); - - /* IMPORTANT: these values are written in WRITE_HEADER_V0 */ - nxout = grid->nx; - nyout = grid->ny; - nzout = grid->nz; - dxout = grid->dx; - dyout = grid->dy; - dzout = grid->dz; - - WRITE_HEADER_V0( dump_type::field_dump, -1, 0, fileIO ); - - dim[0] = grid->nx+2; - dim[1] = grid->ny+2; - dim[2] = grid->nz+2; - WRITE_ARRAY_HEADER( field_array->f, 3, dim, fileIO ); - fileIO.write( field_array->f, dim[0]*dim[1]*dim[2] ); - if( fileIO.close() ) ERROR(( "File close failed on dump fields!!!" )); -} - -void -vpic_simulation::dump_hydro( const char *sp_name, - const char *fbase, - int ftag ) { - species_t *sp; - char fname[256]; - FileIO fileIO; - int dim[3]; - - sp = find_species_name( sp_name, species_list ); - if( !sp ) ERROR(( "Invalid species \"%s\"", sp_name )); - - clear_hydro_array( hydro_array ); - accumulate_hydro_p( hydro_array, sp, interpolator_array ); - synchronize_hydro_array( hydro_array ); - - if( !fbase ) ERROR(( "Invalid filename" )); - - if( rank()==0 ) - MESSAGE(("Dumping \"%s\" hydro fields to \"%s\"",sp->name,fbase)); - - if( ftag ) sprintf( fname, "%s.%li.%i", fbase, (long)step(), rank() ); - else sprintf( fname, "%s.%i", fbase, rank() ); - FileIOStatus status = fileIO.open(fname, io_write); - if( status==fail) ERROR(( "Could not open \"%s\".", fname )); - - /* IMPORTANT: these values are written in WRITE_HEADER_V0 */ - nxout = grid->nx; - nyout = grid->ny; - nzout = grid->nz; - dxout = grid->dx; - dyout = grid->dy; - dzout = grid->dz; - - WRITE_HEADER_V0( dump_type::hydro_dump,sp->id,sp->q/sp->m,fileIO); - - dim[0] = grid->nx+2; - dim[1] = grid->ny+2; - dim[2] = grid->nz+2; - WRITE_ARRAY_HEADER( hydro_array->h, 3, dim, fileIO ); - fileIO.write( hydro_array->h, dim[0]*dim[1]*dim[2] ); - if( fileIO.close() ) ERROR(( "File close failed on dump hydro!!!" )); -} - -void -vpic_simulation::dump_particles( const char *sp_name, - const char *fbase, - int ftag ) -{ - species_t *sp; - char fname[256]; - FileIO fileIO; - int dim[1], buf_start; - static particle_t * ALIGNED(128) p_buf = NULL; -# define PBUF_SIZE 32768 // 1MB of particles - - sp = find_species_name( sp_name, species_list ); - if( !sp ) ERROR(( "Invalid species name \"%s\".", sp_name )); - - if( !fbase ) ERROR(( "Invalid filename" )); - - if( !p_buf ) MALLOC_ALIGNED( p_buf, PBUF_SIZE, 128 ); - - if( rank()==0 ) - MESSAGE(("Dumping \"%s\" particles to \"%s\"",sp->name,fbase)); - - if( ftag ) sprintf( fname, "%s.%li.%i", fbase, (long)step(), rank() ); - else sprintf( fname, "%s.%i", fbase, rank() ); - FileIOStatus status = fileIO.open(fname, io_write); - if( status==fail ) ERROR(( "Could not open \"%s\"", fname )); - - /* IMPORTANT: these values are written in WRITE_HEADER_V0 */ - nxout = grid->nx; - nyout = grid->ny; - nzout = grid->nz; - dxout = grid->dx; - dyout = grid->dy; - dzout = grid->dz; - - WRITE_HEADER_V0( dump_type::particle_dump, sp->id, sp->q/sp->m, fileIO ); - - dim[0] = sp->np; - WRITE_ARRAY_HEADER( p_buf, 1, dim, fileIO ); - - // Copy a PBUF_SIZE hunk of the particle list into the particle - // buffer, timecenter it and write it out. This is done this way to - // guarantee the particle list unchanged while not requiring too - // much memory. - - // FIXME: WITH A PIPELINED CENTER_P, PBUF NOMINALLY SHOULD BE QUITE - // LARGE. - - particle_t * sp_p = sp->p; sp->p = p_buf; - int sp_np = sp->np; sp->np = 0; - int sp_max_np = sp->max_np; sp->max_np = PBUF_SIZE; - for( buf_start=0; buf_startnp = sp_np-buf_start; if( sp->np > PBUF_SIZE ) sp->np = PBUF_SIZE; - COPY( sp->p, &sp_p[buf_start], sp->np ); - center_p( sp, interpolator_array ); - fileIO.write( sp->p, sp->np ); - } - sp->p = sp_p; - sp->np = sp_np; - sp->max_np = sp_max_np; - - if( fileIO.close() ) ERROR(("File close failed on dump particles!!!")); -} - /*------------------------------------------------------------------------------ * New dump logic *---------------------------------------------------------------------------*/ @@ -516,6 +443,8 @@ vpic_simulation::global_header( const char * base, if( fileIO.close() ) ERROR(( "File close failed on global header!!!" )); } +// TODO: why is there field_dump and dump_fields? +// TODO: this could probably move into the dump_strategy void vpic_simulation::field_dump( DumpParameters & dumpParams ) { @@ -554,12 +483,12 @@ vpic_simulation::field_dump( DumpParameters & dumpParams ) { # define f(x,y,z) f[ VOXEL(x,y,z, grid->nx,grid->ny,grid->nz) ] /* IMPORTANT: these values are written in WRITE_HEADER_V0 */ - nxout = (grid->nx)/istride; - nyout = (grid->ny)/jstride; - nzout = (grid->nz)/kstride; - dxout = (grid->dx)*istride; - dyout = (grid->dy)*jstride; - dzout = (grid->dz)*kstride; + size_t nxout = (grid->nx)/istride; + size_t nyout = (grid->ny)/jstride; + size_t nzout = (grid->nz)/kstride; + float dxout = (grid->dx)*istride; + float dyout = (grid->dy)*jstride; + float dzout = (grid->dz)*kstride; /* Banded output will write data as a single block-array as opposed to * the Array-of-Structure format that is used for native storage. @@ -571,7 +500,7 @@ vpic_simulation::field_dump( DumpParameters & dumpParams ) { if(dumpParams.format == band) { - WRITE_HEADER_V0(dump_type::field_dump, -1, 0, fileIO); + WRITE_HEADER_V0(dump_type::field_dump, -1, 0, fileIO, step(), rank(), nproc()); dim[0] = nxout+2; dim[1] = nyout+2; @@ -632,7 +561,7 @@ vpic_simulation::field_dump( DumpParameters & dumpParams ) { } else { // band_interleave - WRITE_HEADER_V0(dump_type::field_dump, -1, 0, fileIO); + WRITE_HEADER_V0(dump_type::field_dump, -1, 0, fileIO, step(), rank(), nproc()); dim[0] = nxout+2; dim[1] = nyout+2; @@ -699,16 +628,13 @@ vpic_simulation::hydro_dump( const char * speciesname, int dim[3]; - /* define to do C-style indexing */ -# define hydro(x,y,z) hydro_array->h[VOXEL(x,y,z, grid->nx,grid->ny,grid->nz)] - /* IMPORTANT: these values are written in WRITE_HEADER_V0 */ - nxout = (grid->nx)/istride; - nyout = (grid->ny)/jstride; - nzout = (grid->nz)/kstride; - dxout = (grid->dx)*istride; - dyout = (grid->dy)*jstride; - dzout = (grid->dz)*kstride; + size_t nxout = (grid->nx)/istride; + size_t nyout = (grid->ny)/jstride; + size_t nzout = (grid->nz)/kstride; + float dxout = (grid->dx)*istride; + float dyout = (grid->dy)*jstride; + float dzout = (grid->dz)*kstride; /* Banded output will write data as a single block-array as opposed to * the Array-of-Structure format that is used for native storage. @@ -720,7 +646,7 @@ vpic_simulation::hydro_dump( const char * speciesname, */ if(dumpParams.format == band) { - WRITE_HEADER_V0(dump_type::hydro_dump, sp->id, sp->q/sp->m, fileIO); + WRITE_HEADER_V0(dump_type::hydro_dump, sp->id, sp->q/sp->m, fileIO, step(), rank(), nproc()); dim[0] = nxout+2; dim[1] = nyout+2; @@ -764,7 +690,7 @@ vpic_simulation::hydro_dump( const char * speciesname, } else { // band_interleave - WRITE_HEADER_V0(dump_type::hydro_dump, sp->id, sp->q/sp->m, fileIO); + WRITE_HEADER_V0(dump_type::hydro_dump, sp->id, sp->q/sp->m, fileIO, step(), rank(), nproc()); dim[0] = nxout; dim[1] = nyout; diff --git a/src/vpic/dump.h b/src/vpic/dump.h new file mode 100644 index 00000000..1d17ee8a --- /dev/null +++ b/src/vpic/dump.h @@ -0,0 +1,19 @@ +#ifndef dump_h +#define dump_h + +#include +#include "../grid/grid.h" + +// TODO: should this be an enum? +namespace dump_type { + const int grid_dump = 0; + const int field_dump = 1; + const int hydro_dump = 2; + const int particle_dump = 3; + const int restart_dump = 4; + const int history_dump = 5; +} // namespace + +// TODO: namesapce? +std::array global_particle_index(int local_i, grid_t* grid, int rank); +#endif diff --git a/src/vpic/dump_strategy.cc b/src/vpic/dump_strategy.cc new file mode 100644 index 00000000..e03ca36c --- /dev/null +++ b/src/vpic/dump_strategy.cc @@ -0,0 +1,161 @@ +//BinaryDump::BinaryDump(int _rank, int _nproc) : Dump_Strategy(_rank, _nproc) +//{ + //// empty +//} +#include "dump_strategy.h" + +void BinaryDump::dump_fields( + const char *fbase, + int step, + grid_t* grid, + field_array_t* field_array, + int ftag + ) +{ + char fname[256]; + FileIO fileIO; + int dim[3]; + + if( !fbase ) ERROR(( "Invalid filename" )); + + if( rank==0 ) MESSAGE(( "Dumping fields to \"%s\"", fbase )); + + if( ftag ) sprintf( fname, "%s.%li.%i", fbase, (long)step, rank ); + else sprintf( fname, "%s.%i", fbase, rank ); + + FileIOStatus status = fileIO.open(fname, io_write); + if( status==fail ) ERROR(( "Could not open \"%s\".", fname )); + + /* IMPORTANT: these values are written in WRITE_HEADER_V0 */ + size_t nxout = grid->nx; + size_t nyout = grid->ny; + size_t nzout = grid->nz; + float dxout = grid->dx; + float dyout = grid->dy; + float dzout = grid->dz; + + WRITE_HEADER_V0( dump_type::field_dump, -1, 0, fileIO, step , rank, nproc); + + dim[0] = grid->nx+2; + dim[1] = grid->ny+2; + dim[2] = grid->nz+2; + WRITE_ARRAY_HEADER( field_array->f, 3, dim, fileIO ); + fileIO.write( field_array->f, dim[0]*dim[1]*dim[2] ); + if( fileIO.close() ) ERROR(( "File close failed on dump fields!!!" )); +} + +void BinaryDump::dump_particles( + const char *fbase, + species_t* sp, + grid_t* grid, + int step, + interpolator_array_t* interpolator_array, + int ftag + ) +{ + char fname[256]; + FileIO fileIO; + int dim[1], buf_start; + static particle_t * ALIGNED(128) p_buf = NULL; + + // TODO: reconcile this with MAX_IO_CHUNK, and update Cmake option + // description to explain what backends use it +# define PBUF_SIZE 32768 // 1MB of particles + + if( !sp ) ERROR(( "Invalid species name \"%s\".", sp->name )); + + if( !fbase ) ERROR(( "Invalid filename" )); + + if( !p_buf ) MALLOC_ALIGNED( p_buf, PBUF_SIZE, 128 ); + + if( rank==0 ) + MESSAGE(("Dumping \"%s\" particles to \"%s\"",sp->name,fbase)); + + if( ftag ) sprintf( fname, "%s.%li.%i", fbase, (long)step, rank ); + else sprintf( fname, "%s.%i", fbase, rank ); + FileIOStatus status = fileIO.open(fname, io_write); + if( status==fail ) ERROR(( "Could not open \"%s\"", fname )); + + /* IMPORTANT: these values are written in WRITE_HEADER_V0 */ + size_t nxout = grid->nx; + size_t nyout = grid->ny; + size_t nzout = grid->nz; + float dxout = grid->dx; + float dyout = grid->dy; + float dzout = grid->dz; + + WRITE_HEADER_V0( dump_type::particle_dump, sp->id, sp->q/sp->m, fileIO, step, rank, nproc); + + dim[0] = sp->np; + WRITE_ARRAY_HEADER( p_buf, 1, dim, fileIO ); + + // Copy a PBUF_SIZE hunk of the particle list into the particle + // buffer, timecenter it and write it out. This is done this way to + // guarantee the particle list unchanged while not requiring too + // much memory. + + // FIXME: WITH A PIPELINED CENTER_P, PBUF NOMINALLY SHOULD BE QUITE + // LARGE. + + particle_t * sp_p = sp->p; sp->p = p_buf; + int sp_np = sp->np; sp->np = 0; + int sp_max_np = sp->max_np; sp->max_np = PBUF_SIZE; + for( buf_start=0; buf_startnp = sp_np-buf_start; if( sp->np > PBUF_SIZE ) sp->np = PBUF_SIZE; + COPY( sp->p, &sp_p[buf_start], sp->np ); + center_p( sp, interpolator_array ); + fileIO.write( sp->p, sp->np ); + } + sp->p = sp_p; + sp->np = sp_np; + sp->max_np = sp_max_np; + + if( fileIO.close() ) ERROR(("File close failed on dump particles!!!")); +} +void BinaryDump::dump_hydro( + const char *fbase, + int step, + hydro_array_t* hydro_array, + species_t* sp, + interpolator_array_t* interpolator_array, + grid_t* grid, + int ftag + ) +{ + char fname[256]; + FileIO fileIO; + int dim[3]; + + if( !sp ) ERROR(( "Invalid species \"%s\"", sp->name )); + + clear_hydro_array( hydro_array ); + accumulate_hydro_p( hydro_array, sp, interpolator_array ); + synchronize_hydro_array( hydro_array ); + + if( !fbase ) ERROR(( "Invalid filename" )); + + if( rank==0 ) + MESSAGE(("Dumping \"%s\" hydro fields to \"%s\"",sp->name,fbase)); + + if( ftag ) sprintf( fname, "%s.%li.%i", fbase, (long)step, rank ); + else sprintf( fname, "%s.%i", fbase, rank ); + FileIOStatus status = fileIO.open(fname, io_write); + if( status==fail) ERROR(( "Could not open \"%s\".", fname )); + + /* IMPORTANT: these values are written in WRITE_HEADER_V0 */ + size_t nxout = grid->nx; + size_t nyout = grid->ny; + size_t nzout = grid->nz; + float dxout = grid->dx; + float dyout = grid->dy; + float dzout = grid->dz; + + WRITE_HEADER_V0(dump_type::hydro_dump, sp->id, sp->q/sp->m, fileIO, step, rank, nproc); + + dim[0] = grid->nx+2; + dim[1] = grid->ny+2; + dim[2] = grid->nz+2; + WRITE_ARRAY_HEADER( hydro_array->h, 3, dim, fileIO ); + fileIO.write( hydro_array->h, dim[0]*dim[1]*dim[2] ); + if( fileIO.close() ) ERROR(( "File close failed on dump hydro!!!" )); +} diff --git a/src/vpic/dump_strategy.h b/src/vpic/dump_strategy.h new file mode 100644 index 00000000..85f70b1e --- /dev/null +++ b/src/vpic/dump_strategy.h @@ -0,0 +1,2047 @@ +#ifndef Dump_Strategy_h +#define Dump_Strategy_h + +#include +#include +#include + +//#define DUMP_INFO_DEBUG 1 +//#define H5_ASYNC 1 +#ifdef H5_ASYNC +#include "h5_vol_external_async_native.h" +#endif +//#define CHUNK_FLAG 1 + +//#define METADATA_COLL_WRITE 1 +//#define TRUE 1 + +#define HAS_FIELD_COMP 1 +#define HAS_PARTICLE_COMP 1 +#define HAS_HYDRO_COMP 1 + +//#define HAS_INDEPENDENT_IO 1 + +#include +#include // TODO: it would be good if this didn't have to know about MPI + +// TODO: should I drop the ./src here? +#include "../field_advance/field_advance.h" +#include "../sf_interface/sf_interface.h" +#include "../species_advance/species_advance.h" +#include "../util/io/FileIO.h" +#include "../util/io/FileUtils.h" +#include "../util/util_base.h" + +#include "dump.h" +#include "dumpmacros.h" + +#ifdef VPIC_ENABLE_HDF5 +#include "hdf5.h" // from the lib +#include "hdf5_header_info.h" // from vpic +#endif + +#ifdef VPIC_ENABLE_OPENPMD +#include +#endif + +//#define N_FILE_N_PROCESS 1 +//#define TEST_MPIIO 1 + +// TODO: delete this +#define _LOG_PREFIX \ + __FILE__ "(" EXPAND_AND_STRINGIFY(__LINE__) ")[" << rank << "]: " + +/* +#define io_log(x) do { \ +if( rank==0 ) { \ +std::cerr << _LOG_PREFIX << x << std::endl; \ +std::cerr.flush(); \ +} \ +} while(0) +*/ + +// Runtime inheritance is obviously not very "VPIC like", as we will [probably] +// incur a penalty for the vtable lookup, but given we're about to do IO this +// is very negligible. +class Dump_Strategy { +public: + int rank, nproc, num_step; + + Dump_Strategy(int _rank, int _nproc) : rank(_rank), nproc(_nproc) {} // empty + + virtual ~Dump_Strategy(){}; + + virtual void dump_fields(const char *fbase, int step, grid_t *grid, + field_array_t *field_array, int ftag) = 0; + virtual void dump_hydro(const char *fbase, int step, + hydro_array_t *hydro_array, species_t *sp, + interpolator_array_t *interpolator_array, + grid_t *grid, int ftag) = 0; + virtual void dump_particles(const char *fbase, species_t *sp, grid_t *grid, + int step, + interpolator_array_t *interpolator_array, + int ftag) = 0; +}; + +class BinaryDump : public Dump_Strategy { +public: + using Dump_Strategy::Dump_Strategy; // inherit constructor + // BinaryDump(int _rank, int _nproc ) : Dump_Strategy(_rank, _nproc ){ } // + // empty + + // TODO: now we pass rank and step, ftag has odd semanticds + void dump_fields(const char *fbase, int step, grid_t *grid, + field_array_t *field_array, int ftag); + void dump_hydro(const char *fbase, int step, hydro_array_t *hydro_array, + species_t *sp, interpolator_array_t *interpolator_array, + grid_t *grid, int ftag); + void dump_particles(const char *fbase, species_t *sp, grid_t *grid, int step, + interpolator_array_t *interpolator_array, int ftag); +}; + +#ifdef VPIC_ENABLE_HDF5 + +struct field_dump_flag_t { + bool ex = true, ey = true, ez = true, div_e_err = true; + bool cbx = true, cby = true, cbz = true, div_b_err = true; + bool tcax = true, tcay = true, tcaz = true, rhob = true; + bool jfx = true, jfy = true, jfz = true, rhof = true; + bool ematx = true, ematy = true, ematz = true, nmat = true; + bool fmatx = true, fmaty = true, fmatz = true, cmat = true; + void disableE() { ex = false, ey = false, ez = false, div_e_err = false; } + + void disableCB() { cbx = false, cby = false, cbz = false, div_b_err = false; } + + void disableTCA() { tcax = false, tcay = false, tcaz = false, rhob = false; } + + void disableJF() { jfx = false, jfy = false, jfz = false, rhof = false; } + + void disableEMAT() { + ematx = false, ematy = false, ematz = false, nmat = false; + } + + void disableFMAT() { + fmatx = false, fmaty = false, fmatz = false, cmat = false; + } + + void resetToDefaults() { + ex = true, ey = true, ez = true, div_e_err = true; + cbx = true, cby = true, cbz = true, div_b_err = true; + tcax = true, tcay = true, tcaz = true, rhob = true; + jfx = true, jfy = true, jfz = true, rhof = true; + ematx = true, ematy = true, ematz = true, nmat = true; + fmatx = true, fmaty = true, fmatz = true, cmat = true; + } + + bool enabledE() { return ex && ey && ez; } + + bool enabledCB() { return cbx && cby && cbz; } + + bool enabledTCA() { return tcax && tcay && tcaz; } + + bool enabledJF() { return jfx && jfy && jfz; } + + bool enabledEMAT() { return ematx && ematy && ematz; } + + bool enabledFMAT() { return fmatx && fmaty && fmatz; } +}; + +struct hydro_dump_flag_t { + bool jx = true, jy = true, jz = true, rho = true; + bool px = true, py = true, pz = true, ke = true; + bool txx = true, tyy = true, tzz = true; + bool tyz = true, tzx = true, txy = true; + + void disableJ() { jx = false, jy = false, jz = false, rho = false; } + + void disableP() { px = false, py = false, pz = false, ke = false; } + + void disableTD() // Stress diagonal + { + txx = false, tyy = false, tzz = false; + } + + void disableTOD() // Stress off-diagonal + { + tyz = false, tzx = false, txy = false; + } + void resetToDefaults() { + jx = true, jy = true, jz = true, rho = true; + px = true, py = true, pz = true, ke = true; + txx = true, tyy = true, tzz = true; + tyz = true, tzx = true, txy = true; + } + + bool enabledJ() { return jx && jy && jz; } + + bool enabledP() { return px && py && pz; } + + bool enabledTD() { return txx && tyy && tzz; } + + bool enabledTOD() { return tyz && tzx && txy; } +}; +class HDF5Dump : public Dump_Strategy { + std::unordered_map tframe_map; + +public: + using Dump_Strategy::Dump_Strategy; // inherit constructor + + // TODO: replace these with a common dump interface + // Declare vars to use + hydro_dump_flag_t hydro_dump_flag; + field_dump_flag_t field_dump_flag; + +#define DUMP_DIR_FORMAT "./%s" + + // TODO: naming a macro so close to existing functions AND data is not a good + // define to do C-style indexing +#define _hydro(x, y, z) \ + hydro_array->h[VOXEL(x, y, z, grid->nx, grid->ny, grid->nz)] + + /** + * @brief Dump field data to the HDf5 file + * Author: Bin Dong dbin@lbl.gov + * https://crd.lbl.gov/bin-dong + * Nov 2020 + * @param fbase + * @param step + * @param grid + * @param field_array + * @param ftag + */ + void dump_fields(const char *fbase, int step, grid_t *grid, + field_array_t *field_array, int ftag) { + size_t step_for_viou = step; + + int mpi_size, mpi_rank; + MPI_Comm_size(MPI_COMM_WORLD, &mpi_size); + MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank); + +#ifdef DUMP_INFO_DEBUG + printf("MPI rank = %d, size = %d \n", mpi_rank, mpi_size); + // printf("base dir for field: %s \n", fdParams.baseDir); + // printf("stride x y z = (%ld, %ld, %ld)\n", fdParams.stride_x, + // fdParams.stride_y, fdParams.stride_z); + printf("grid x, y z = (%d, %d, %d) \n", grid->nx, grid->ny, grid->nz); + printf("domain loc (x0, y0, z0) -> (x1, y1, z1) = (%f, %f, %f) -> (%f, %f, " + "%f) \n", + grid->x0, grid->y0, grid->z0, grid->x1, grid->y1, grid->z1); + // printf("global->topology_x, y, z = %f, %f, %f \n ", global->topology_x, + // global->topology_y, global->topology_z); + printf("grid -> sx, sy, sz = (%d, %d, %d), nv=%d \n", grid->sx, grid->sy, + grid->sz, grid->nv); +#endif + + char fname[256]; + char field_scratch[128]; + char subfield_scratch[128]; + + sprintf(field_scratch, DUMP_DIR_FORMAT, "field_hdf5"); + FileUtils::makeDirectory(field_scratch); + sprintf(subfield_scratch, "%s/T.%zu/", field_scratch, step_for_viou); + FileUtils::makeDirectory(subfield_scratch); + + sprintf(fname, "%s/%s_%zu.h5", subfield_scratch, "fields", step_for_viou); + double el1 = uptime(); + + // int file_exist(const char *filename) + //{ + // struct stat buffer; + // return (stat(filename, &buffer) == 0); + //} + + // https://support.hdfgroup.org/ftp/HDF5/current/src/unpacked/examples/h5_compound.c +#ifdef HAS_FIELD_COMP + if (!mpi_rank) + printf("Using Field Compund type !\n"); + hid_t field_comp_type_it = H5Tcreate(H5T_COMPOUND, sizeof(field_t)); + H5Tinsert(field_comp_type_it, "ex", HOFFSET(field_t, ex), H5T_NATIVE_FLOAT); + H5Tinsert(field_comp_type_it, "ey", HOFFSET(field_t, ey), H5T_NATIVE_FLOAT); + H5Tinsert(field_comp_type_it, "ez", HOFFSET(field_t, ez), H5T_NATIVE_FLOAT); + H5Tinsert(field_comp_type_it, "div_e_err", HOFFSET(field_t, div_e_err), + H5T_NATIVE_FLOAT); + + H5Tinsert(field_comp_type_it, "cbx", HOFFSET(field_t, cbx), + H5T_NATIVE_FLOAT); + H5Tinsert(field_comp_type_it, "cby", HOFFSET(field_t, cby), + H5T_NATIVE_FLOAT); + H5Tinsert(field_comp_type_it, "cbz", HOFFSET(field_t, cbz), + H5T_NATIVE_FLOAT); + H5Tinsert(field_comp_type_it, "div_b_err", HOFFSET(field_t, div_b_err), + H5T_NATIVE_FLOAT); + + H5Tinsert(field_comp_type_it, "tcax", HOFFSET(field_t, tcax), + H5T_NATIVE_FLOAT); + H5Tinsert(field_comp_type_it, "tcay", HOFFSET(field_t, tcay), + H5T_NATIVE_FLOAT); + H5Tinsert(field_comp_type_it, "tcaz", HOFFSET(field_t, tcaz), + H5T_NATIVE_FLOAT); + H5Tinsert(field_comp_type_it, "rhob", HOFFSET(field_t, rhob), + H5T_NATIVE_FLOAT); + + H5Tinsert(field_comp_type_it, "jfx", HOFFSET(field_t, jfx), + H5T_NATIVE_FLOAT); + H5Tinsert(field_comp_type_it, "jfy", HOFFSET(field_t, jfy), + H5T_NATIVE_FLOAT); + H5Tinsert(field_comp_type_it, "jfz", HOFFSET(field_t, jfz), + H5T_NATIVE_FLOAT); + H5Tinsert(field_comp_type_it, "rhof", HOFFSET(field_t, rhof), + H5T_NATIVE_FLOAT); + + H5Tinsert(field_comp_type_it, "ematx", HOFFSET(field_t, ematx), + H5T_NATIVE_SHORT); + H5Tinsert(field_comp_type_it, "ematy", HOFFSET(field_t, ematy), + H5T_NATIVE_SHORT); + H5Tinsert(field_comp_type_it, "ematz", HOFFSET(field_t, ematz), + H5T_NATIVE_SHORT); + H5Tinsert(field_comp_type_it, "nmat", HOFFSET(field_t, nmat), + H5T_NATIVE_SHORT); + + H5Tinsert(field_comp_type_it, "fmatx", HOFFSET(field_t, fmatx), + H5T_NATIVE_SHORT); + H5Tinsert(field_comp_type_it, "fmaty", HOFFSET(field_t, fmaty), + H5T_NATIVE_SHORT); + H5Tinsert(field_comp_type_it, "fmatz", HOFFSET(field_t, fmatz), + H5T_NATIVE_SHORT); + H5Tinsert(field_comp_type_it, "cmat", HOFFSET(field_t, cmat), + H5T_NATIVE_SHORT); +#endif + + // struct stat buffer; + // if((stat(fname, &buffer) == 0)){ + // file_exist_flag = 1; + // if(!mpi_rank) + // printf("Write original files /w HDF5! \n"); + // } + // file_exist_flag = 0; + + hid_t plist_id; + hid_t file_id; + plist_id = H5Pcreate(H5P_FILE_ACCESS); + H5Pset_fapl_mpio(plist_id, MPI_COMM_WORLD, MPI_INFO_NULL); + // H5Pset_alignment(plist_id, 4194304, 4194304); + /*if(H5Pset_libver_bounds(plist_id, H5F_LIBVER_LATEST, H5F_LIBVER_LATEST) < + 0){ exit(-1); + }*/ + +#ifdef METADATA_COLL_WRITE + if (!mpi_rank) + printf("Enable collective metadata write !\n"); + H5Pset_coll_metadata_write(plist_id, TRUE); +#endif + file_id = H5Fcreate(fname, H5F_ACC_TRUNC, H5P_DEFAULT, plist_id); + H5Pclose(plist_id); + + sprintf(fname, "Timestep_%zu", step_for_viou); + hid_t group_id; + group_id = H5Gcreate(file_id, fname, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + + el1 = uptime() - el1; + // io_log("TimeHDF5Open): " << el1 << " s"); //Easy to handle results for + // scripts + double el2 = uptime(); + + /* + // Create a variable list of field values to output. + size_t numvars = std::min(global->fdParams.output_vars.bitsum(), + total_field_variables); size_t * varlist = new size_t[numvars]; + + for(size_t i(0), c(0); ifdParams.output_vars.bitset(i)) varlist[c++] = i; + + printf("\nBEGIN_OUTPUT: numvars = %zd \n", numvars);*/ + +#define fpp(x, y, z) f[VOXEL(x, y, z, grid->nx, grid->ny, grid->nz)] + /* + typedef struct field { + float ex, ey, ez, div_e_err; // Electric field and div E error + float cbx, cby, cbz, div_b_err; // Magnetic field and div B error + float tcax, tcay, tcaz, rhob; // TCA fields and bound charge + density float jfx, jfy, jfz, rhof; // Free current and charge + density material_id ematx, ematy, ematz, nmat; // Material at edge + centers and nodes material_id fmatx, fmaty, fmatz, cmat; // Material at + face and cell centers } field_t;*/ + // Local voxel mesh resolution. Voxels are + // indexed FORTRAN style 0:nx+1,0:ny+1,0:nz+1 + // with voxels 1:nx,1:ny,1:nz being non-ghost + // voxels. + + float *temp_buf = + (float *)malloc(sizeof(float) * (grid->nx) * (grid->ny) * (grid->nz)); + hsize_t temp_buf_index; + hid_t dset_id; + // char *field_var_name[] = + // {"ex","ey","ez","div_e_err","cbx","cby","cbz","div_b_err","tcax","tcay","tcaz","rhob","jfx","jfy","jfz","rhof"}; + // Comment out for test only + + plist_id = H5Pcreate(H5P_DATASET_XFER); +#ifdef HAS_INDEPENDENT_IO + if (!mpi_rank) + printf("\n ###\n VPIC Independent I/O! \n ###\n"); + H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT); +#else + H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE); +#endif + + // H5Sselect_hyperslab(filespace, H5S_SELECT_SET, (hsize_t *) &offset, NULL, + // (hsize_t *) &numparticles, NULL); + + // global->topology_x + + hsize_t field_global_size[3], field_local_size[3], global_offset[3], + global_count[3]; + field_global_size[0] = (grid->nx * grid->gpx); + field_global_size[1] = (grid->ny * grid->gpy); + field_global_size[2] = (grid->nz * grid->gpz); + + field_local_size[0] = grid->nx; + field_local_size[1] = grid->ny; + field_local_size[2] = grid->nz; + + int gpx = grid->gpx; + int gpy = grid->gpy; + int gpz = grid->gpz; + + // Convert rank to local decomposition + int rx, ry, rz; + UNVOXEL(mpi_rank, rx, ry, rz, grid->gpx, grid->gpy, grid->gpz); + + int mpi_rank_x, mpi_rank_y, mpi_rank_z; + mpi_rank_x = rx; + mpi_rank_y = ry; + mpi_rank_z = rz; + + global_offset[0] = (grid->nx) * mpi_rank_x; + global_offset[1] = (grid->ny) * mpi_rank_y; + global_offset[2] = (grid->nz) * mpi_rank_z; + + global_count[0] = (grid->nx); + global_count[1] = (grid->ny); + global_count[2] = (grid->nz); + +#ifdef DUMP_INFO_DEBUG + if (mpi_rank < 4) { + printf("grid nx, ny nz = (%d, %d, %d) \n", grid->nx, grid->ny, grid->nz); + printf("mpi-rank = %d, rank index = (%d, %d, %d) \n", mpi_rank, + mpi_rank_x, mpi_rank_y, mpi_rank_z); + printf("global size = %llu %llu %llu \n", field_global_size[0], + field_global_size[1], field_global_size[2]); + printf("global_offset = %llu %llu %llu \n", global_offset[0], + global_offset[1], global_offset[2]); + printf("global_count = %llu %llu %llu \n", global_count[0], + global_count[1], global_count[2]); + fflush(stdout); + } +#endif + + hid_t filespace; //= H5Screate_simple(3, field_global_size, NULL); + hid_t memspace; // = H5Screate_simple(3, field_local_size, NULL); + // if(!file_exist_flag){ + filespace = H5Screate_simple(3, field_global_size, NULL); + //} + memspace = H5Screate_simple(3, field_local_size, NULL); + + hsize_t chunk_dims[3]; + chunk_dims[0] = 288; // grid->nx; //8 x 8 x 8 + chunk_dims[1] = 24; // grid->ny; // + chunk_dims[2] = 24; // grid->nz; + + hid_t dataspace_id; + hid_t dcpl_id = H5Pcreate(H5P_DATASET_CREATE); +#ifdef CHUNK_FLAG + H5Pset_chunk(dcpl_id, 3, chunk_dims); + if (!mpi_rank) + printf("Enable chunking !\n"); +#endif + +#define DUMP_FIELD_TO_HDF5(DSET_NAME, ATTRIBUTE_NAME, ELEMENT_TYPE) \ + { \ + dset_id = H5Dcreate(group_id, DSET_NAME, ELEMENT_TYPE, filespace, \ + H5P_DEFAULT, dcpl_id, H5P_DEFAULT); \ + temp_buf_index = 0; \ + for (size_t i(1); i < grid->nx + 1; i++) { \ + for (size_t j(1); j < grid->ny + 1; j++) { \ + for (size_t k(1); k < grid->nz + 1; k++) { \ + temp_buf[temp_buf_index] = \ + FIELD_ARRAY_NAME->fpp(i, j, k).ATTRIBUTE_NAME; \ + temp_buf_index = temp_buf_index + 1; \ + } \ + } \ + } \ + dataspace_id = H5Dget_space(dset_id); \ + H5Sselect_hyperslab(dataspace_id, H5S_SELECT_SET, global_offset, NULL, \ + global_count, NULL); \ + H5Dwrite(dset_id, ELEMENT_TYPE, memspace, dataspace_id, plist_id, \ + temp_buf); \ + H5Sclose(dataspace_id); \ + H5Dclose(dset_id); \ + } + /* + typedef struct field { + float ex, ey, ez, div_e_err; // Electric field and div E + error float cbx, cby, cbz, div_b_err; // Magnetic field and div + B error float tcax, tcay, tcaz, rhob; // TCA fields and bound + charge density float jfx, jfy, jfz, rhof; // Free current + and charge density material_id ematx, ematy, ematz, nmat; // Material + at edge centers and nodes material_id fmatx, fmaty, fmatz, cmat; // + Material at face and cell centers } field_t;*/ + +#ifdef HAS_FIELD_COMP + field_t *field_buf; + temp_buf_index = 0; + int global_index; + field_buf = (field_t *)malloc(sizeof(field_t) * (grid->nx) * (grid->ny) * + (grid->nz)); + for (size_t i(1); i < grid->nx + 1; i++) { + for (size_t j(1); j < grid->ny + 1; j++) { + for (size_t k(1); k < grid->nz + 1; k++) { + field_buf[temp_buf_index] = FIELD_ARRAY_NAME->fpp(i, j, k); + temp_buf_index++; + } + } + } + dset_id = H5Dcreate(group_id, "field", field_comp_type_it, filespace, + H5P_DEFAULT, dcpl_id, H5P_DEFAULT); + dataspace_id = H5Dget_space(dset_id); + H5Sselect_hyperslab(dataspace_id, H5S_SELECT_SET, global_offset, NULL, + global_count, NULL); + H5Dwrite(dset_id, field_comp_type_it, memspace, dataspace_id, plist_id, + field_buf); + free(field_buf); + H5Sclose(dataspace_id); + H5Dclose(dset_id); + H5Tclose(field_comp_type_it); +#else + + if (field_dump_flag.ex) + DUMP_FIELD_TO_HDF5("ex", ex, H5T_NATIVE_FLOAT); + if (field_dump_flag.ey) + DUMP_FIELD_TO_HDF5("ey", ey, H5T_NATIVE_FLOAT); + if (field_dump_flag.ez) + DUMP_FIELD_TO_HDF5("ez", ez, H5T_NATIVE_FLOAT); + if (field_dump_flag.div_e_err) + DUMP_FIELD_TO_HDF5("div_e_err", div_e_err, H5T_NATIVE_FLOAT); + + if (field_dump_flag.cbx) + DUMP_FIELD_TO_HDF5("cbx", cbx, H5T_NATIVE_FLOAT); + if (field_dump_flag.cby) + DUMP_FIELD_TO_HDF5("cby", cby, H5T_NATIVE_FLOAT); + if (field_dump_flag.cbz) + DUMP_FIELD_TO_HDF5("cbz", cbz, H5T_NATIVE_FLOAT); + if (field_dump_flag.div_b_err) + DUMP_FIELD_TO_HDF5("div_b_err", div_b_err, H5T_NATIVE_FLOAT); + + if (field_dump_flag.tcax) + DUMP_FIELD_TO_HDF5("tcax", tcax, H5T_NATIVE_FLOAT); + if (field_dump_flag.tcay) + DUMP_FIELD_TO_HDF5("tcay", tcay, H5T_NATIVE_FLOAT); + if (field_dump_flag.tcaz) + DUMP_FIELD_TO_HDF5("tcaz", tcaz, H5T_NATIVE_FLOAT); + if (field_dump_flag.rhob) + DUMP_FIELD_TO_HDF5("rhob", rhob, H5T_NATIVE_FLOAT); + + if (field_dump_flag.jfx) + DUMP_FIELD_TO_HDF5("jfx", jfx, H5T_NATIVE_FLOAT); + if (field_dump_flag.jfy) + DUMP_FIELD_TO_HDF5("jfy", jfy, H5T_NATIVE_FLOAT); + if (field_dump_flag.jfz) + DUMP_FIELD_TO_HDF5("jfz", jfz, H5T_NATIVE_FLOAT); + if (field_dump_flag.rhof) + DUMP_FIELD_TO_HDF5("rhof", rhof, H5T_NATIVE_FLOAT); + + // H5T_NATIVE_SHORT for material_id (typedef int16_t material_id) + if (field_dump_flag.ematx) + DUMP_FIELD_TO_HDF5("ematx", ematx, H5T_NATIVE_SHORT); + if (field_dump_flag.ematy) + DUMP_FIELD_TO_HDF5("ematy", ematy, H5T_NATIVE_SHORT); + if (field_dump_flag.ematz) + DUMP_FIELD_TO_HDF5("ematz", ematz, H5T_NATIVE_SHORT); + if (field_dump_flag.nmat) + DUMP_FIELD_TO_HDF5("nmat", nmat, H5T_NATIVE_SHORT); + + if (field_dump_flag.fmatx) + DUMP_FIELD_TO_HDF5("fmatx", fmatx, H5T_NATIVE_SHORT); + if (field_dump_flag.fmaty) + DUMP_FIELD_TO_HDF5("fmaty", fmaty, H5T_NATIVE_SHORT); + if (field_dump_flag.fmatz) + DUMP_FIELD_TO_HDF5("fmatz", fmatz, H5T_NATIVE_SHORT); + if (field_dump_flag.cmat) + DUMP_FIELD_TO_HDF5("cmat", cmat, H5T_NATIVE_SHORT); + +#endif + + H5D_mpio_actual_io_mode_t actual_io_mode; + H5Pget_mpio_actual_io_mode(plist_id, &actual_io_mode); + /* + + switch(actual_io_mode){ + case H5D_MPIO_NO_COLLECTIVE: + io_log("H5Pget_mpio_actual_io_mode: H5D_MPIO_NO_COLLECTIVE: "); + break; + case H5D_MPIO_CHUNK_INDEPENDENT: + io_log("H5Pget_mpio_actual_io_mode: H5D_MPIO_CHUNK_INDEPENDENT: "); + break; + case H5D_MPIO_CHUNK_COLLECTIVE: + io_log("H5Pget_mpio_actual_io_mode: H5D_MPIO_CHUNK_COLLECTIVE: "); + break; + case H5D_MPIO_CHUNK_MIXED: + io_log("H5Pget_mpio_actual_io_mode: H5D_MPIO_CHUNK_MIXED: "); + break; + case H5D_MPIO_CONTIGUOUS_COLLECTIVE: + io_log("H5Pget_mpio_actual_io_mode: H5D_MPIO_CONTIGUOUS_COLLECTIVE: "); + break; + default : + io_log("H5Pget_mpio_actual_io_mode: None returend: "); + break; + } + + H5D_mpio_actual_chunk_opt_mode_t actual_chunk_opt_mode; + H5Pget_mpio_actual_chunk_opt_mode(plist_id, &actual_chunk_opt_mode); + switch(actual_chunk_opt_mode){ + case H5D_MPIO_NO_CHUNK_OPTIMIZATION: + io_log("H5Pget_mpio_actual_chunk_opt_mode: +H5D_MPIO_NO_CHUNK_OPTIMIZATION: "); break; case H5D_MPIO_MULTI_CHUNK: + io_log("H5Pget_mpio_actual_chunk_opt_mode: H5D_MPIO_MULTI_CHUNK: "); + break; + // case H5D_MPIO_MULTI_CHUNK_NO_OPT: + // io_log("H5Pget_mpio_actual_chunk_opt_mode: +H5D_MPIO_MULTI_CHUNK_NO_OPT: "); + // break; + case H5D_MPIO_LINK_CHUNK: + io_log("H5Pget_mpio_actual_chunk_opt_mode: H5D_MPIO_LINK_CHUNK: "); + break; + default : + io_log("H5Pget_mpio_actual_chunk_opt_mode: None returend: "); + break; + } + + uint32_t local_no_collective_cause, global_no_collective_cause; + H5Pget_mpio_no_collective_cause(plist_id, &local_no_collective_cause, +&global_no_collective_cause); + + switch(local_no_collective_cause){ + case H5D_MPIO_COLLECTIVE: + io_log("local_no_collective_cause: H5D_MPIO_COLLECTIVE: "); + break; + case H5D_MPIO_SET_INDEPENDENT: + io_log("local_no_collective_cause: H5D_MPIO_SET_INDEPENDENT: "); + break; + case H5D_MPIO_DATA_TRANSFORMS: + io_log("local_no_collective_cause: H5D_MPIO_DATA_TRANSFORMS: "); + break; + //case H5D_MPIO_SET_MPIPOSIX: + // io_log("local_no_collective_cause: H5D_MPIO_SET_MPIPOSIX: "); + // break; + case H5D_MPIO_NOT_SIMPLE_OR_SCALAR_DATASPACES: + io_log("local_no_collective_cause: H5D_MPIO_NOT_SIMPLE_OR_SCALAR_DATASPACES: +"); break; + //case H5D_MPIO_POINT_SELECTIONS: + // io_log("local_no_collective_cause: H5D_MPIO_POINT_SELECTIONS: "); + // break; + // case H5D_MPIO_FILTERS: + // io_log("local_no_collective_cause: H5D_MPIO_FILTERS: "); + // break; + default : + io_log("local_no_collective_cause: None returend: "); + break; +} + + +switch(global_no_collective_cause){ + case H5D_MPIO_COLLECTIVE: + io_log("global_no_collective_cause: H5D_MPIO_COLLECTIVE: "); + break; + case H5D_MPIO_SET_INDEPENDENT: + io_log("global_no_collective_cause: H5D_MPIO_SET_INDEPENDENT: "); + break; + case H5D_MPIO_DATA_TRANSFORMS: + io_log("global_no_collective_cause: H5D_MPIO_DATA_TRANSFORMS: "); + break; + //case H5D_MPIO_SET_MPIPOSIX: + // io_log("global_no_collective_cause: H5D_MPIO_SET_MPIPOSIX: "); + // break; + case H5D_MPIO_NOT_SIMPLE_OR_SCALAR_DATASPACES: + io_log("global_no_collective_cause: +H5D_MPIO_NOT_SIMPLE_OR_SCALAR_DATASPACES: "); break; + //case H5D_MPIO_POINT_SELECTIONS: + // io_log("global_no_collective_cause: H5D_MPIO_POINT_SELECTIONS: "); + // break; + // case H5D_MPIO_FILTERS: + // io_log("global_no_collective_cause: H5D_MPIO_FILTERS: "); + // break; + default : + io_log("global_no_collective_cause: None returend: "); + break; +} +*/ + + el2 = uptime() - el2; + // io_log("TimeHDF5Write: " << el2 << " s"); + + double el3 = uptime(); + + /* + //Write metadata (geo original and geo dx/dy/dz) for ArrayUDF + float attr_data[2][3]; + attr_data[0][0] = grid->x0; + attr_data[0][1] = grid->y0; + attr_data[0][2] = grid->z0; + attr_data[1][0] = grid->dx; + attr_data[1][1] = grid->dy; + attr_data[1][2] = grid->dz; + hsize_t dims[2]; + dims[0] = 2; + dims[1] = 3; + if(!file_exist_flag){ + hid_t va_geo_dataspace_id = H5Screate_simple(2, dims, NULL); + hid_t va_geo_attribute_id = H5Acreate2(file_id, "VPIC-ArrayUDF-GEO", + H5T_IEEE_F32BE, va_geo_dataspace_id, H5P_DEFAULT, H5P_DEFAULT); + H5Awrite(va_geo_attribute_id, H5T_NATIVE_FLOAT, attr_data); + H5Sclose(va_geo_dataspace_id); + H5Aclose(va_geo_attribute_id); + } + */ + free(temp_buf); + // if(!file_exist_flag) + H5Sclose(filespace); + H5Sclose(memspace); + H5Pclose(plist_id); + H5Gclose(group_id); + H5Fclose(file_id); + + el3 = uptime() - el3; + // io_log("TimeHDF5Close: " << el3 << " s"); + + if (mpi_rank == 0) { + char const *output_xml_file = "./field_hdf5/hdf5_field.xdmf"; + char dimensions_3d[128]; + sprintf(dimensions_3d, "%lld %lld %lld", field_global_size[0], + field_global_size[1], field_global_size[2]); + char dimensions_4d[128]; + sprintf(dimensions_4d, "%lld %lld %lld %d", field_global_size[0], + field_global_size[1], field_global_size[2], 3); + char orignal[128]; + sprintf(orignal, "%f %f %f", grid->x0, grid->y0, grid->z0); + char dxdydz[128]; + sprintf(dxdydz, "%f %f %f", grid->dx, grid->dy, grid->dz); + + // TODO: remove or let the user set + int field_interval = 1; + + // TODO: remove this dependence on number of steps + // std::cout << "num_step " << num_step << std::endl; + + int nframes = num_step / field_interval + 1; + static int field_tframe = 0; + +#ifdef DUMP_INFO_DEBUG + printf(" meta file : %s \n", output_xml_file); + printf(" array dims per var: %s \n", dimensions_3d); + printf("array dims all vars: %s \n", dimensions_4d); + printf(" orignal: %s \n", orignal); + printf(" dxdydz: %s \n", dxdydz); + printf(" nframes: %d \n", nframes); + printf(" field_interval: %d \n", field_interval); + printf(" current step: %zd \n", step_for_viou); + printf(" current step: %zd \n", step_for_viou); + + // printf(" Simulation time: %f \n", grid->t0); + printf(" tframe: %d \n", field_tframe); +#endif + + // TODO: this footer dumping is more likely better done in a + // destructor, rather than hoping a multiple division works out + if (field_tframe >= 1) { + if (field_tframe == (nframes - 1)) { + invert_field_xml_item(output_xml_file, "fields", step_for_viou, + dimensions_4d, dimensions_3d, 1); + } else { + invert_field_xml_item(output_xml_file, "fields", step_for_viou, + dimensions_4d, dimensions_3d, 0); + } + } else { + create_file_with_header(output_xml_file, dimensions_3d, orignal, dxdydz, + nframes, field_interval); + if (field_tframe == (nframes - 1)) { + invert_field_xml_item(output_xml_file, "fields", step_for_viou, + dimensions_4d, dimensions_3d, 1); + } else { + invert_field_xml_item(output_xml_file, "fields", step_for_viou, + dimensions_4d, dimensions_3d, 0); + } + } + field_tframe++; + } + } + /** + * @brief dump_particles to the HDF5 file + * Author: Bin Dong dbin@lbl.gov + * https://crd.lbl.gov/bin-dong + * Nov 2020 + * @param fbase + * @param sp + * @param grid + * @param step + * @param interpolator_array + * @param ftag + */ + void dump_particles(const char *fbase, species_t *sp, grid_t *grid, int step, + interpolator_array_t *interpolator_array, int ftag) { + static int file_index = 0; + file_index++; + int mpi_rank; + MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank); + double dump_particles_uptime = uptime(); + time_t seconds = time(NULL); + // printf("Atrank = %d, file_index = %d, dump_particles_uptime = %f, + // epoch_seconds = %ld \n ", mpi_rank, file_index, dump_particles_uptime, + // seconds); + + size_t step_for_viou = step; + char fname[256]; + char group_name[256]; + char particle_scratch[128]; + char subparticle_scratch[128]; + + int np_local; + + float *Pf; + int *Pi; + + // get the total number of particles. in this example, output only electrons + // sp = species_list; + sprintf(particle_scratch, DUMP_DIR_FORMAT, "particle_hdf5"); + FileUtils::makeDirectory(particle_scratch); + sprintf(subparticle_scratch, "%s/T.%ld/", particle_scratch, step_for_viou); + FileUtils::makeDirectory(subparticle_scratch); + + // TODO: Allow the user to set this + int stride_particle_dump = 1; + + np_local = (sp->np + stride_particle_dump - 1) / stride_particle_dump; + + // make a copy of the part of particle data to be dumped + double ec1 = uptime(); + + int sp_np = sp->np; + int sp_max_np = sp->max_np; + particle_t *ALIGNED(128) p_buf = NULL; + if (!p_buf) + MALLOC_ALIGNED(p_buf, np_local, 128); + particle_t *sp_p = sp->p; + sp->p = p_buf; + sp->np = np_local; + sp->max_np = np_local; + + for (long long iptl = 0, i = 0; iptl < sp_np; + iptl += stride_particle_dump, ++i) { + COPY(&sp->p[i], &sp_p[iptl], 1); + } + + center_p(sp, interpolator_array); + + ec1 = uptime() - ec1; + + // if(!mpi_rank || mpi_rank == 2047 ) + // std::cout << "on mpi_rank: " << mpi_rank << ", time in copying + // particle data: " << ec1 << " s" << ", np_local = " << np_local << + // std::endl; + +#ifndef N_FILE_N_PROCESS + int np_local_max, np_local_min; + MPI_Reduce(&np_local, &np_local_max, 1, MPI_INT, MPI_MAX, 0, + MPI_COMM_WORLD); + MPI_Reduce(&np_local, &np_local_min, 1, MPI_INT, MPI_MIN, 0, + MPI_COMM_WORLD); + // io_log("on mpi_rank: " << mpi_rank << ", time in copying particle data: " + // << ec1 << " s" << ", np_local = " << np_local << ",np_local_max = " << + // np_local_max << ", local_min = "<< np_local_min); +#endif + + Pf = (float *)sp->p; + Pi = (int *)sp->p; + + // open HDF5 file in "particle/T./" subdirectory + // filename: eparticle.h5p +#ifndef N_FILE_N_PROCESS + sprintf(fname, "%s/%s_%ld.h5", subparticle_scratch, sp->name, + step_for_viou); +#else + sprintf(fname, "%s/%s_%ld_p%d.h5", subparticle_scratch, sp->name, + step_for_viou, mpi_rank); +#endif + + sprintf(group_name, "/Timestep_%ld", step_for_viou); + double el1 = uptime(); + + long long total_particles, offset; + long long numparticles = np_local; + +#ifndef N_FILE_N_PROCESS + MPI_Allreduce(&numparticles, &total_particles, 1, MPI_LONG_LONG, MPI_SUM, + MPI_COMM_WORLD); + MPI_Scan(&numparticles, &offset, 1, MPI_LONG_LONG, MPI_SUM, MPI_COMM_WORLD); + offset -= numparticles; +#else + total_particles = np_local; + offset = 0; +#endif + + hid_t file_plist_id = H5Pcreate(H5P_FILE_ACCESS); + +#ifndef N_FILE_N_PROCESS + H5Pset_fapl_mpio(file_plist_id, MPI_COMM_WORLD, MPI_INFO_NULL); +#endif + +#ifdef H5_ASYNC + if (!mpi_rank) + printf("Enable async on particle data"); + + assert(H5Pset_vol_async(file_plist_id)); +#endif + + hid_t file_id = H5Fcreate(fname, H5F_ACC_TRUNC, H5P_DEFAULT, file_plist_id); + // if(!mpi_rank ) + // io_log("++Particle H5Fcreate) "); + + hid_t group_id = + H5Gcreate(file_id, group_name, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + // if(!mpi_rank ) + // io_log("++Particle H5Gcreate) "); + +#ifdef HAS_PARTICLE_COMP + if (!mpi_rank) + printf("Using Partilce Compund type !\n"); + hid_t particle_comp_type_it = H5Tcreate(H5T_COMPOUND, sizeof(particle_t)); + H5Tinsert(particle_comp_type_it, "dx", HOFFSET(particle_t, dx), + H5T_NATIVE_FLOAT); + H5Tinsert(particle_comp_type_it, "dy", HOFFSET(particle_t, dy), + H5T_NATIVE_FLOAT); + H5Tinsert(particle_comp_type_it, "dz", HOFFSET(particle_t, dz), + H5T_NATIVE_FLOAT); + + H5Tinsert(particle_comp_type_it, "i", HOFFSET(particle_t, i), + H5T_NATIVE_INT); + + H5Tinsert(particle_comp_type_it, "ux", HOFFSET(particle_t, ux), + H5T_NATIVE_FLOAT); + H5Tinsert(particle_comp_type_it, "uy", HOFFSET(particle_t, uy), + H5T_NATIVE_FLOAT); + H5Tinsert(particle_comp_type_it, "uz", HOFFSET(particle_t, uz), + H5T_NATIVE_FLOAT); + H5Tinsert(particle_comp_type_it, "w", HOFFSET(particle_t, w), + H5T_NATIVE_FLOAT); +#endif + + hid_t filespace = H5Screate_simple(1, (hsize_t *)&total_particles, NULL); + H5Sselect_hyperslab(filespace, H5S_SELECT_SET, (hsize_t *)&offset, NULL, + (hsize_t *)&numparticles, NULL); + + // if(!mpi_rank ) + // io_log("++Particle H5Sselect_hyperslab) "); + + // plist_id = H5P_DEFAULT; + hid_t io_plist_id = H5Pcreate(H5P_DATASET_XFER); + +#ifndef N_FILE_N_PROCESS +#ifdef HAS_INDEPENDENT_IO + if (!mpi_rank) { + printf("\n ###\n VPIC Independent I/O! \n ###\n"); + } + H5Pset_dxpl_mpio(io_plist_id, H5FD_MPIO_INDEPENDENT); +#else + H5Pset_dxpl_mpio(io_plist_id, H5FD_MPIO_COLLECTIVE); +#endif +#endif + +#ifdef H5_ASYNC + H5Pset_dxpl_async(io_plist_id, true); +#endif + hsize_t linearspace_count_temp = numparticles; + hid_t linearspace = H5Screate_simple(1, &linearspace_count_temp, NULL); + + hsize_t memspace_count_temp; + hid_t memspace; +#ifdef HAS_PARTICLE_COMP + memspace_count_temp = numparticles; + memspace = H5Screate_simple(1, &memspace_count_temp, NULL); +#else + memspace_count_temp = numparticles * 8; + memspace = H5Screate_simple(1, &memspace_count_temp, NULL); + hsize_t memspace_start = 0, memspace_stride = 8, memspace_count = np_local; + H5Sselect_hyperslab(memspace, H5S_SELECT_SET, &memspace_start, + &memspace_stride, &memspace_count, NULL); +#endif + el1 = uptime() - el1; + // if(!mpi_rank || mpi_rank == 2047 ) + // io_log("Particle TimeHDF5Open): " << el1 << " s"); //Easy to handle + // results for scripts + double el2 = uptime(); + int ierr; + +#define WRITE_H5_FILE(group_id_p, data_buf_p, type_p, dname_p) \ + { \ + hid_t dset_id = H5Dcreate(group_id_p, dname_p, type_p, filespace, \ + H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); \ + H5Dwrite(dset_id, type_p, memspace, filespace, io_plist_id, data_buf_p); \ + H5Dclose(dset_id); \ + } + + // MPI_Info_set(info, "romio_cb_write", "disable"); +#define WRITE_MPI_FILE(dname_p, offset_p, data_buf_p, count_p, type_p) \ + { \ + MPI_File fh; \ + MPI_Status status; \ + sprintf(fname, "%s/%s_%ld_%s.h5", subparticle_scratch, sp->name, \ + step_for_viou, dname_p); \ + if (mpi_rank == 0) \ + printf("fname= %s \n", fname); \ + MPI_Info info; \ + MPI_Info_create(&info); \ + MPI_File_open(MPI_COMM_WORLD, fname, MPI_MODE_WRONLY | MPI_MODE_CREATE, \ + info, &fh); \ + MPI_File_write_at(fh, offset_p, data_buf_p, count_p, type_p, &status); \ + MPI_Info_free(&info); \ + MPI_File_close(&fh); \ + } + +#ifdef HAS_PARTICLE_COMP + hid_t dset_id = H5Dcreate(group_id, "particle", particle_comp_type_it, + filespace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + H5Dwrite(dset_id, particle_comp_type_it, memspace, filespace, io_plist_id, + sp->p); + H5Dclose(dset_id); +#else +#ifdef TEST_MPIIO + // Here we don't use the stripe but just for performance test + if (!mpi_rank) + printf("Test MPI-IO\n"); + WRITE_MPI_FILE("dX", offset * sizeof(float), Pf, numparticles, MPI_FLOAT); + WRITE_MPI_FILE("dY", offset * sizeof(float), Pf, numparticles, MPI_FLOAT); + WRITE_MPI_FILE("dZ", offset * sizeof(float), Pf, numparticles, MPI_FLOAT); + WRITE_MPI_FILE("i", offset * sizeof(int), Pf, numparticles, MPI_INT); + WRITE_MPI_FILE("ux", offset * sizeof(float), Pf, numparticles, MPI_FLOAT); + WRITE_MPI_FILE("uy", offset * sizeof(float), Pf, numparticles, MPI_FLOAT); + WRITE_MPI_FILE("uz", offset * sizeof(float), Pf, numparticles, MPI_FLOAT); + WRITE_MPI_FILE("q", offset * sizeof(float), Pf, numparticles, MPI_FLOAT); +#else +#ifndef N_FILE_N_PROCESS + if (!mpi_rank) + printf("Test HDF5-IO Single \n"); +#else + if (!mpi_rank) + printf("Test HDF5-IO N Files N Process\n"); +#endif + // if(!mpi_rank ) + // io_log("++Particle Starting to write ) "); + WRITE_H5_FILE(group_id, Pf, H5T_NATIVE_FLOAT, "dX") + WRITE_H5_FILE(group_id, Pf + 1, H5T_NATIVE_FLOAT, "dY") + WRITE_H5_FILE(group_id, Pf + 2, H5T_NATIVE_FLOAT, "dZ") + WRITE_H5_FILE(group_id, Pi + 3, H5T_NATIVE_INT, "i") + WRITE_H5_FILE(group_id, Pf + 4, H5T_NATIVE_FLOAT, "ux") + WRITE_H5_FILE(group_id, Pf + 5, H5T_NATIVE_FLOAT, "uy") + WRITE_H5_FILE(group_id, Pf + 6, H5T_NATIVE_FLOAT, "uz") + WRITE_H5_FILE(group_id, Pf + 7, H5T_NATIVE_FLOAT, "q") +#endif +#endif + el2 = uptime() - el2; + // io_log("Particle TimeHDF5Write: " << el2 << " s"); + + double el3 = uptime(); + H5Sclose(memspace); + H5Sclose(filespace); + H5Pclose(file_plist_id); + H5Pclose(io_plist_id); + H5Gclose(group_id); + + H5Fclose(file_id); + +#ifdef H5_ASYNC + H5VLasync_finalize(); +#endif + el3 = uptime() - el3; + // io_log("Particle TimeHDF5Close: " << el3 << " s"); + } + + /** + * @brief Dump hydro data to the HDf5 file + * Author: Bin Dong dbin@lbl.gov + * https://crd.lbl.gov/bin-dong + * Nov 2020 + * @param fbase + * @param step + * @param hydro_array + * @param sp + * @param interpolator_array + * @param grid + * @param ftag + */ + void dump_hydro(const char *fbase, int step, hydro_array_t *hydro_array, + species_t *sp, interpolator_array_t *interpolator_array, + grid_t *grid, int ftag) { + size_t step_for_viou = step; + +#define DUMP_HYDRO_TO_HDF5(DSET_NAME, ATTRIBUTE_NAME, ELEMENT_TYPE) \ + { \ + dset_id = H5Dcreate(group_id, DSET_NAME, ELEMENT_TYPE, filespace, \ + H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); \ + temp_buf_index = 0; \ + for (size_t i(1); i < grid->nx + 1; i++) { \ + for (size_t j(1); j < grid->ny + 1; j++) { \ + for (size_t k(1); k < grid->nz + 1; k++) { \ + temp_buf[temp_buf_index] = _hydro(i, j, k).ATTRIBUTE_NAME; \ + temp_buf_index = temp_buf_index + 1; \ + } \ + } \ + } \ + dataspace_id = H5Dget_space(dset_id); \ + H5Sselect_hyperslab(dataspace_id, H5S_SELECT_SET, global_offset, NULL, \ + global_count, NULL); \ + H5Dwrite(dset_id, ELEMENT_TYPE, memspace, dataspace_id, plist_id, \ + temp_buf); \ + H5Sclose(dataspace_id); \ + H5Dclose(dset_id); \ + } + //#define DUMP_INFO_DEBUG 1 + int mpi_size, mpi_rank; + MPI_Comm_size(MPI_COMM_WORLD, &mpi_size); + MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank); + + if (!sp) { + ERROR(("Invalid species")); + } + + clear_hydro_array(hydro_array); + accumulate_hydro_p(hydro_array, sp, interpolator_array); + synchronize_hydro_array(hydro_array); + + char hname[256]; + char hydro_scratch[128]; + char subhydro_scratch[128]; + + sprintf(hydro_scratch, "./%s", "hydro_hdf5"); + FileUtils::makeDirectory(hydro_scratch); + sprintf(subhydro_scratch, "%s/T.%zu/", hydro_scratch, step_for_viou); + FileUtils::makeDirectory(subhydro_scratch); + + sprintf(hname, "%s/hydro_%s_%zu.h5", subhydro_scratch, sp->name, + step_for_viou); + double el1 = uptime(); + hid_t plist_id = H5Pcreate(H5P_FILE_ACCESS); + + /* + if(H5Pset_libver_bounds(plist_id, H5F_LIBVER_LATEST, H5F_LIBVER_LATEST) < + 0){ exit(-1); + }*/ + // if((fid = H5Fcreate(FILENAME, H5F_ACC_TRUNC, H5P_DEFAULT, plist_id)) < 0) + // ERROR_RETURN; + + H5Pset_fapl_mpio(plist_id, MPI_COMM_WORLD, MPI_INFO_NULL); + hid_t file_id = H5Fcreate(hname, H5F_ACC_TRUNC, H5P_DEFAULT, plist_id); + H5Pclose(plist_id); + + sprintf(hname, "Timestep_%zu", step_for_viou); + hid_t group_id = + H5Gcreate(file_id, hname, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + + el1 = uptime() - el1; + // io_log("TimeHDF5Open: " << el1 << " s"); //Easy to handle results for + // scripts + double el2 = uptime(); + + // Create a variable list of field values to output. + // size_t numvars = std::min(global->fdParams.output_vars.bitsum(), + // total_field_variables); size_t *varlist = new size_t[numvars]; + + // for (size_t i(0), c(0); i < total_field_variables; i++) + // if (global->fdParams.output_vars.bitset(i)) + // varlist[c++] = i; + + // printf("\nBEGIN_OUTPUT: numvars = %zu \n", numvars); + + // typedef struct hydro { + // float jx, jy, jz, rho; // Current and charge density => , + // float px, py, pz, ke; // Momentum and K.E. density => , float txx, tyy, tzz; // Stress diagonal => + // , i==j float tyz, tzx, txy; // Stress off-diagonal => , i!=j float _pad[2]; // 16-byte align + //} hydro_t; +#ifdef HAS_HYDRO_COMP + // if(!mpi_rank) + // printf("Using Field Compund type !\n"); + hid_t hydro_comp_type_it = H5Tcreate(H5T_COMPOUND, sizeof(hydro_t)); + H5Tinsert(hydro_comp_type_it, "jx", HOFFSET(hydro_t, jx), H5T_NATIVE_FLOAT); + H5Tinsert(hydro_comp_type_it, "jy", HOFFSET(hydro_t, jy), H5T_NATIVE_FLOAT); + H5Tinsert(hydro_comp_type_it, "jz", HOFFSET(hydro_t, jz), H5T_NATIVE_FLOAT); + H5Tinsert(hydro_comp_type_it, "rho", HOFFSET(hydro_t, rho), + H5T_NATIVE_FLOAT); + + H5Tinsert(hydro_comp_type_it, "px", HOFFSET(hydro_t, px), H5T_NATIVE_FLOAT); + H5Tinsert(hydro_comp_type_it, "py", HOFFSET(hydro_t, py), H5T_NATIVE_FLOAT); + H5Tinsert(hydro_comp_type_it, "pz", HOFFSET(hydro_t, pz), H5T_NATIVE_FLOAT); + H5Tinsert(hydro_comp_type_it, "ke", HOFFSET(hydro_t, ke), H5T_NATIVE_FLOAT); + + H5Tinsert(hydro_comp_type_it, "txx", HOFFSET(hydro_t, txx), + H5T_NATIVE_FLOAT); + H5Tinsert(hydro_comp_type_it, "tyy", HOFFSET(hydro_t, tyy), + H5T_NATIVE_FLOAT); + H5Tinsert(hydro_comp_type_it, "tzz", HOFFSET(hydro_t, tzz), + H5T_NATIVE_FLOAT); + + H5Tinsert(hydro_comp_type_it, "tyz", HOFFSET(hydro_t, tyz), + H5T_NATIVE_FLOAT); + H5Tinsert(hydro_comp_type_it, "tzx", HOFFSET(hydro_t, tzx), + H5T_NATIVE_FLOAT); + H5Tinsert(hydro_comp_type_it, "txy", HOFFSET(hydro_t, txy), + H5T_NATIVE_FLOAT); + H5Tinsert(hydro_comp_type_it, "pad", HOFFSET(hydro_t, _pad), + H5T_NATIVE_DOUBLE); +#endif + // typedef struct hydro_array { + // hydro_t * ALIGNED(128) h; + // grid_t * g; + //} hydro_array_t; + + float *temp_buf = + (float *)malloc(sizeof(float) * (grid->nx) * (grid->ny) * (grid->nz)); + hsize_t temp_buf_index; + hid_t dset_id; + // char *field_var_name[] = + // {"ex","ey","ez","div_e_err","cbx","cby","cbz","div_b_err","tcax","tcay","tcaz","rhob","jfx","jfy","jfz","rhof"}; + plist_id = H5Pcreate(H5P_DATASET_XFER); + // Comment out for test only + H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE); + // H5Sselect_hyperslab(filespace, H5S_SELECT_SET, (hsize_t *) &offset, NULL, + // (hsize_t *) &numparticles, NULL); + + // global->topology_x + + hsize_t hydro_global_size[3], hydro_local_size[3], global_offset[3], + global_count[3]; + hydro_global_size[0] = (grid->nx * grid->gpx); + hydro_global_size[1] = (grid->ny * grid->gpy); + hydro_global_size[2] = (grid->nz * grid->gpz); + + hydro_local_size[0] = grid->nx; + hydro_local_size[1] = grid->ny; + hydro_local_size[2] = grid->nz; + + int mpi_rank_x, mpi_rank_y, mpi_rank_z; + UNVOXEL(mpi_rank, mpi_rank_x, mpi_rank_y, mpi_rank_z, grid->gpx, grid->gpy, + grid->gpz); + + global_offset[0] = (grid->nx) * mpi_rank_x; + global_offset[1] = (grid->ny) * mpi_rank_y; + global_offset[2] = (grid->nz) * mpi_rank_z; + + global_count[0] = (grid->nx); + global_count[1] = (grid->ny); + global_count[2] = (grid->nz); + +#ifdef DUMP_INFO_DEBUG + printf("global size = %llu %llu %llu \n", hydro_global_size[0], + hydro_global_size[1], hydro_global_size[2]); + printf("global_offset = %llu %llu %llu \n", global_offset[0], + global_offset[1], global_offset[2]); + printf("global_count = %llu %llu %llu \n", global_count[0], + global_count[1], global_count[2]); + printf("mpi-rank = %d, rank index = (%d, %d, %d) \n", mpi_rank, mpi_rank_x, + mpi_rank_y, mpi_rank_z); + fflush(stdout); +#endif + + hid_t filespace = H5Screate_simple(3, hydro_global_size, NULL); + hid_t memspace = H5Screate_simple(3, hydro_local_size, NULL); + hid_t dataspace_id; + + // typedef struct hydro { + // float jx, jy, jz, rho; // Current and charge density => , + // float px, py, pz, ke; // Momentum and K.E. density => , float txx, tyy, tzz; // Stress diagonal => + // , i==j float tyz, tzx, txy; // Stress off-diagonal => , i!=j float _pad[2]; // 16-byte align + //} hydro_t; + +#ifdef HAS_HYDRO_COMP + hydro_t *hydro_buf = (hydro_t *)malloc(sizeof(hydro_t) * (grid->nx) * + (grid->ny) * (grid->nz)); + temp_buf_index = 0; + for (size_t i(1); i < grid->nx + 1; i++) { + for (size_t j(1); j < grid->ny + 1; j++) { + for (size_t k(1); k < grid->nz + 1; k++) { + hydro_buf[temp_buf_index] = _hydro(i, j, k); + temp_buf_index = temp_buf_index + 1; + } + } + } + dset_id = H5Dcreate(group_id, "hydro", hydro_comp_type_it, filespace, + H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + + dataspace_id = H5Dget_space(dset_id); + H5Sselect_hyperslab(dataspace_id, H5S_SELECT_SET, global_offset, NULL, + global_count, NULL); + H5Dwrite(dset_id, hydro_comp_type_it, memspace, dataspace_id, plist_id, + hydro_buf); + free(hydro_buf); + H5Sclose(dataspace_id); + H5Dclose(dset_id); + H5Tclose(hydro_comp_type_it); +#else + + if (hydro_dump_flag.jx) + DUMP_HYDRO_TO_HDF5("jx", jx, H5T_NATIVE_FLOAT); + if (hydro_dump_flag.jy) + DUMP_HYDRO_TO_HDF5("jy", jy, H5T_NATIVE_FLOAT); + if (hydro_dump_flag.jz) + DUMP_HYDRO_TO_HDF5("jz", jz, H5T_NATIVE_FLOAT); + if (hydro_dump_flag.rho) + DUMP_HYDRO_TO_HDF5("rho", rho, H5T_NATIVE_FLOAT); + + if (hydro_dump_flag.px) + DUMP_HYDRO_TO_HDF5("px", px, H5T_NATIVE_FLOAT); + if (hydro_dump_flag.py) + DUMP_HYDRO_TO_HDF5("py", py, H5T_NATIVE_FLOAT); + if (hydro_dump_flag.pz) + DUMP_HYDRO_TO_HDF5("pz", pz, H5T_NATIVE_FLOAT); + if (hydro_dump_flag.ke) + DUMP_HYDRO_TO_HDF5("ke", ke, H5T_NATIVE_FLOAT); + + if (hydro_dump_flag.txx) + DUMP_HYDRO_TO_HDF5("txx", txx, H5T_NATIVE_FLOAT); + if (hydro_dump_flag.tyy) + DUMP_HYDRO_TO_HDF5("tyy", tyy, H5T_NATIVE_FLOAT); + if (hydro_dump_flag.tzz) + DUMP_HYDRO_TO_HDF5("tzz", tzz, H5T_NATIVE_FLOAT); + + if (hydro_dump_flag.tyz) + DUMP_HYDRO_TO_HDF5("tyz", tyz, H5T_NATIVE_FLOAT); + if (hydro_dump_flag.tzx) + DUMP_HYDRO_TO_HDF5("tzx", tzx, H5T_NATIVE_FLOAT); + if (hydro_dump_flag.txy) + DUMP_HYDRO_TO_HDF5("txy", txy, H5T_NATIVE_FLOAT); + + el2 = uptime() - el2; + // io_log("TimeHDF5Write: " << el2 << " s"); + +#endif + double el3 = uptime(); + + // Write metadata (geo original and geo dx/dy/dz) for ArrayUDF + /* + float attr_data[2][3]; + attr_data[0][0] = grid->x0; + attr_data[0][1] = grid->y0; + attr_data[0][2] = grid->z0; + attr_data[1][0] = grid->dx; + attr_data[1][1] = grid->dy; + attr_data[1][2] = grid->dz; + hsize_t dims[2]; + dims[0] = 2; + dims[1] = 3; + hid_t va_geo_dataspace_id = H5Screate_simple(2, dims, NULL); + hid_t va_geo_attribute_id = H5Acreate2(file_id, "VPIC-ArrayUDF-GEO", + H5T_IEEE_F32BE, va_geo_dataspace_id, H5P_DEFAULT, H5P_DEFAULT); + H5Awrite(va_geo_attribute_id, H5T_NATIVE_FLOAT, attr_data); + H5Sclose(va_geo_dataspace_id); + H5Aclose(va_geo_attribute_id);*/ + + free(temp_buf); + H5Sclose(filespace); + H5Sclose(memspace); + H5Pclose(plist_id); + H5Gclose(group_id); + H5Fclose(file_id); + + el3 = uptime() - el3; + // io_log("TimeHDF5Close: " << el3 << " s"); + + if (mpi_rank == 0) { + char output_xml_file[128]; + sprintf(output_xml_file, "./%s/%s%s%s", "hydro_hdf5", "hydro-", sp->name, + ".xdmf"); + char dimensions_3d[128]; + sprintf(dimensions_3d, "%lld %lld %lld", hydro_global_size[0], + hydro_global_size[1], hydro_global_size[2]); + char dimensions_4d[128]; + sprintf(dimensions_4d, "%lld %lld %lld %d", hydro_global_size[0], + hydro_global_size[1], hydro_global_size[2], 3); + char orignal[128]; + sprintf(orignal, "%f %f %f", grid->x0, grid->y0, grid->z0); + char dxdydz[128]; + sprintf(dxdydz, "%f %f %f", grid->dx, grid->dy, grid->dz); + + // TODO: remove or let user set + int hydro_interval = 1; + + // TODO: remove this dependence on number of steps + int nframes = num_step / hydro_interval + 1; + + const int tframe = tframe_map[sp->id]; + +#ifdef DUMP_INFO_DEBUG + printf(" meta file : %s \n", output_xml_file); + printf(" array dims per var: %s \n", dimensions_3d); + printf("array dims all vars: %s \n", dimensions_4d); + printf(" orignal: %s \n", orignal); + printf(" dxdydz: %s \n", dxdydz); + printf(" nframes: %d \n", nframes); + printf(" hydro_fields_interval: %d \n", hydro_interval); + printf(" current step: %zu \n", step_for_viou); + printf(" Simulation time: %f \n", grid->t0); + printf(" tframe: %d \n", tframe); +#endif + + // TODO: why doesnt this just use the cstr? + char speciesname_new[128]; + sprintf(speciesname_new, "hydro_%s", sp->name); + if (tframe >= 1) { + if (tframe == (nframes - 1)) { + invert_hydro_xml_item(output_xml_file, speciesname_new, step_for_viou, + dimensions_4d, dimensions_3d, 1); + } else { + invert_hydro_xml_item(output_xml_file, speciesname_new, step_for_viou, + dimensions_4d, dimensions_3d, 0); + } + } else { + create_file_with_header(output_xml_file, dimensions_3d, orignal, dxdydz, + nframes, hydro_interval); + if (tframe == (nframes - 1)) { + invert_hydro_xml_item(output_xml_file, speciesname_new, step_for_viou, + dimensions_4d, dimensions_3d, 1); + } else { + invert_hydro_xml_item(output_xml_file, speciesname_new, step_for_viou, + dimensions_4d, dimensions_3d, 0); + } + } + tframe_map[sp->id]++; + } + } +}; +#endif + +#ifdef VPIC_ENABLE_OPENPMD +class OpenPMDDump : public Dump_Strategy { +public: + // openPMD::Series* series; + using Dump_Strategy::Dump_Strategy; // inherit constructor + + // std::string file_type = ".h5"; + std::string file_type = ".bp"; + + void dump_fields(const char *fbase, int step, grid_t *grid, + field_array_t *field_array, int ftag) { + std::cout << "Writing openPMD data" << std::endl; + + std::string full_file_name = fbase + file_type; + + // if (series == nullptr) { + std::cout << "init series" << std::endl; + openPMD::Series series = openPMD::Series( + full_file_name, openPMD::AccessType::CREATE, MPI_COMM_WORLD); + //} + + std::cout << "Writing iteration " << step << std::endl; + auto i = series.iterations[step]; + // TODO: it would be nice to set these... + // series.setAuthor( "Axel Huebl "); + // series.setMachine( "Hall Probe 5000, Model 3"); + i.setAttribute("vacuum", true); + + auto cB = i.meshes["B"]; + auto E = i.meshes["E"]; + auto J = i.meshes["J"]; + auto Tca = i.meshes["Tca"]; + auto Emat = i.meshes["Emat"]; + auto Fmat = i.meshes["Fmat"]; + auto Rho = i.meshes["Rho"]; + auto DivErr = i.meshes["DivErr"]; + + // record components + auto Cbx = cB["x"]; + auto Cby = cB["y"]; + auto Cbz = cB["z"]; + + auto Ex = E["x"]; + auto Ey = E["y"]; + auto Ez = E["z"]; + + auto Jx = J["x"]; + auto Jy = J["y"]; + auto Jz = J["z"]; + + auto Tcax = Tca["x"]; + auto Tcay = Tca["y"]; + auto Tcaz = Tca["z"]; + + auto Ematx = Emat["x"]; + auto Ematy = Emat["y"]; + auto Ematz = Emat["z"]; + + auto Fmatx = Fmat["x"]; + auto Fmaty = Fmat["y"]; + auto Fmatz = Fmat["z"]; + + auto RhoB = Rho["B"]; + auto RhoF = Rho["F"]; + + auto DivEErr = DivErr["E"]; + auto DivBErr = DivErr["B"]; + + // TODO: set unitDimension so the anaylsis software knows what fields + // things are + // + // // TODO: add timers for the convert and for the write + + size_t gnx = (grid->nx * grid->gpx); + size_t gny = (grid->ny * grid->gpy); + size_t gnz = (grid->nz * grid->gpz); + openPMD::Extent global_extent = {gnx, gny, gnz}; + + openPMD::Datatype datatype = openPMD::determineDatatype(); + openPMD::Dataset dataset = openPMD::Dataset(datatype, global_extent); + + Cbx.resetDataset(dataset); + Cby.resetDataset(dataset); + Cbz.resetDataset(dataset); + + Ex.resetDataset(dataset); + Ey.resetDataset(dataset); + Ez.resetDataset(dataset); + + Jx.resetDataset(dataset); + Jy.resetDataset(dataset); + Jz.resetDataset(dataset); + + Tcax.resetDataset(dataset); + Tcay.resetDataset(dataset); + Tcaz.resetDataset(dataset); + + Ematx.resetDataset(dataset); + Ematy.resetDataset(dataset); + Ematz.resetDataset(dataset); + + Fmatx.resetDataset(dataset); + Fmaty.resetDataset(dataset); + Fmatz.resetDataset(dataset); + + RhoB.resetDataset(dataset); + RhoF.resetDataset(dataset); + + DivEErr.resetDataset(dataset); + DivBErr.resetDataset(dataset); + + // TODO: hoist this conversion code, as is it used elsewhere + // Convert rank to local x/y/z + int rx, ry, rz; + UNVOXEL(rank, rx, ry, rz, grid->gpx, grid->gpy, grid->gpz); + + size_t nx = grid->nx; + size_t ny = grid->ny; + size_t nz = grid->nz; + + // NOTE: this assumes a static mesh decomposition in nx/ny/nz + size_t global_offset_x = (nx)*rx; + size_t global_offset_y = (ny)*ry; + size_t global_offset_z = (nz)*rz; + + openPMD::Offset chunk_offset = {global_offset_x, global_offset_y, + global_offset_z}; + openPMD::Extent chunk_extent = {nx, ny, nz}; + + std::cout << "Local offset " + << " x: " << global_offset_x << " y: " << global_offset_y + << " z: " << global_offset_z << std::endl; + + // Store a local copy of the data which we pull out of the AoS + std::vector cbx_data; + std::vector cby_data; + std::vector cbz_data; + + std::vector ex_data; + std::vector ey_data; + std::vector ez_data; + + std::vector jx_data; + std::vector jy_data; + std::vector jz_data; + + std::vector tcax_data; + std::vector tcay_data; + std::vector tcaz_data; + + // TODO: these are material_id (ints not floats) + std::vector ematx_data; + std::vector ematy_data; + std::vector ematz_data; + + std::vector fmatx_data; + std::vector fmaty_data; + std::vector fmatz_data; + // end todo + + std::vector rhob_data; + std::vector rhof_data; + + std::vector divb_data; + std::vector dive_data; + + size_t nv = nx * ny * nz; + + // TODO: resize here will zero out the data which we don't need, we + // could change to a different semantic to avoid this + cbx_data.resize(nv); + cby_data.resize(nv); + cbz_data.resize(nv); + + ex_data.resize(nv); + ey_data.resize(nv); + ez_data.resize(nv); + + jx_data.resize(nv); + jy_data.resize(nv); + jz_data.resize(nv); + + tcax_data.resize(nv); + tcay_data.resize(nv); + tcaz_data.resize(nv); + + ematx_data.resize(nv); + ematy_data.resize(nv); + ematz_data.resize(nv); + + fmatx_data.resize(nv); + fmaty_data.resize(nv); + fmatz_data.resize(nv); + + rhob_data.resize(nv); + rhof_data.resize(nv); + + divb_data.resize(nv); + dive_data.resize(nv); + + // TODO: make this AoS to SoA conversion a function + + // We could do 1D here, but we don't really care about the ghosts, and we + // can thread over nz/ny (collapsed?) + // Go over non-ghosts and grab just that data into a dense array + for (size_t k = 1; k < grid->nz + 1; k++) { + for (size_t j = 1; j < grid->ny + 1; j++) { + for (size_t i = 1; i < grid->nx + 1; i++) { + int local_index = VOXEL(i - 1, j - 1, k - 1, grid->nx - 2, + grid->ny - 2, grid->nz - 2); + int global_index = VOXEL(i, j, k, grid->nx, grid->ny, grid->nz); + + cbx_data[local_index] = field_array->f[global_index].cbx; + cby_data[local_index] = field_array->f[global_index].cby; + cbz_data[local_index] = field_array->f[global_index].cbz; + + ex_data[local_index] = field_array->f[global_index].ex; + ey_data[local_index] = field_array->f[global_index].ey; + ez_data[local_index] = field_array->f[global_index].ez; + + jx_data[local_index] = field_array->f[global_index].jfx; + jy_data[local_index] = field_array->f[global_index].jfy; + jz_data[local_index] = field_array->f[global_index].jfz; + + tcax_data[local_index] = field_array->f[global_index].tcax; + tcay_data[local_index] = field_array->f[global_index].tcay; + tcaz_data[local_index] = field_array->f[global_index].tcaz; + + ematx_data[local_index] = field_array->f[global_index].ematx; + ematy_data[local_index] = field_array->f[global_index].ematy; + ematz_data[local_index] = field_array->f[global_index].ematz; + + fmatx_data[local_index] = field_array->f[global_index].fmatx; + fmaty_data[local_index] = field_array->f[global_index].fmaty; + fmatz_data[local_index] = field_array->f[global_index].fmatz; + + rhob_data[local_index] = field_array->f[global_index].rhob; + rhof_data[local_index] = field_array->f[global_index].rhof; + + dive_data[local_index] = field_array->f[global_index].div_e_err; + divb_data[local_index] = field_array->f[global_index].div_b_err; + } + } + } + + Cbx.storeChunk(cbx_data, chunk_offset, chunk_extent); + Cby.storeChunk(cby_data, chunk_offset, chunk_extent); + Cbz.storeChunk(cbz_data, chunk_offset, chunk_extent); + + Ex.storeChunk(ex_data, chunk_offset, chunk_extent); + Ey.storeChunk(ey_data, chunk_offset, chunk_extent); + Ez.storeChunk(ez_data, chunk_offset, chunk_extent); + + Jx.storeChunk(jx_data, chunk_offset, chunk_extent); + Jy.storeChunk(jy_data, chunk_offset, chunk_extent); + Jz.storeChunk(jz_data, chunk_offset, chunk_extent); + + Tcax.storeChunk(tcax_data, chunk_offset, chunk_extent); + Tcay.storeChunk(tcay_data, chunk_offset, chunk_extent); + Tcaz.storeChunk(tcaz_data, chunk_offset, chunk_extent); + + Ematx.storeChunk(ematx_data, chunk_offset, chunk_extent); + Ematy.storeChunk(ematy_data, chunk_offset, chunk_extent); + Ematz.storeChunk(ematz_data, chunk_offset, chunk_extent); + + Fmatx.storeChunk(fmatx_data, chunk_offset, chunk_extent); + Fmaty.storeChunk(fmaty_data, chunk_offset, chunk_extent); + Fmatz.storeChunk(fmatz_data, chunk_offset, chunk_extent); + + RhoB.storeChunk(rhob_data, chunk_offset, chunk_extent); + RhoF.storeChunk(rhof_data, chunk_offset, chunk_extent); + + DivEErr.storeChunk(dive_data, chunk_offset, chunk_extent); + DivBErr.storeChunk(divb_data, chunk_offset, chunk_extent); + + series.flush(); + } + + void dump_particles(const char *fbase, species_t *sp, grid_t *grid, int step, + interpolator_array_t *interpolator_array, int ftag) { + std::string full_file_name = fbase + file_type; + + std::cout << "writing particles to " << full_file_name << std::endl; + + // if (series == nullptr) { + openPMD::Series series = openPMD::Series( + full_file_name, openPMD::AccessType::CREATE, MPI_COMM_WORLD); + //} + + auto i = series.iterations[step]; + + // TODO: set these + i.setTime((float)step); + i.setDt(1.0); + i.setTimeUnitSI(1.0); + + auto &p = i.particles[sp->name]; + + const int np = sp->np; + + // TODO: this could be a function call as it's used elsewhere (in hdf5) + unsigned long long total_particles, offset; + unsigned long long numparticles = np; + MPI_Allreduce(&numparticles, &total_particles, 1, MPI_LONG_LONG, MPI_SUM, + MPI_COMM_WORLD); + MPI_Scan(&numparticles, &offset, 1, MPI_LONG_LONG, MPI_SUM, MPI_COMM_WORLD); + offset -= numparticles; + + openPMD::Extent global_extent = {total_particles}; + openPMD::Datatype datatype = openPMD::determineDatatype(); + openPMD::Dataset dataset = openPMD::Dataset(datatype, global_extent); + + auto px = p["position"]["x"]; + auto pxo = p["positionOffset"]["x"]; + + auto py = p["position"]["y"]; + auto pyo = p["positionOffset"]["y"]; + + auto pz = p["position"]["z"]; + auto pzo = p["positionOffset"]["z"]; + + auto ux = p["velocity"]["x"]; + auto uy = p["velocity"]["y"]; + auto uz = p["velocity"]["z"]; + + px.resetDataset(dataset); + pxo.resetDataset(dataset); + + py.resetDataset(dataset); + pyo.resetDataset(dataset); + + pz.resetDataset(dataset); + pzo.resetDataset(dataset); + + ux.resetDataset(dataset); + uy.resetDataset(dataset); + uz.resetDataset(dataset); + // convert data to SoA, allowing the user to chunk the operation + + // TODO: Add code the convert to global offsets +#ifndef PMD_MAX_IO_CHUNK // in particles +#define PMD_MAX_IO_CHUNK 16777216; // 512MB total write +#endif + const int max_chunk = PMD_MAX_IO_CHUNK; + + // Loop over all particles in chunks + for (int i = 0; i < np; i += max_chunk) { + // We have to be careful as the last chunk may not be full + // Find how many are left and do that many + size_t to_write = std::min(np - i, max_chunk); + + // Convert the chunk ready to write + std::vector x_pos; + std::vector x_off; + x_pos.resize(to_write); + x_off.resize(to_write); + + std::vector y_pos; + std::vector y_off; + y_pos.resize(to_write); + y_off.resize(to_write); + + std::vector z_pos; + std::vector z_off; + z_pos.resize(to_write); + z_off.resize(to_write); + + std::vector ux_pos; + ux_pos.resize(to_write); + + std::vector uy_pos; + uy_pos.resize(to_write); + + std::vector uz_pos; + uz_pos.resize(to_write); + + for (int j = 0; j < to_write; j++) { + // TODO: do I need to center the particles? + auto &particle = sp->p[i + j]; + + x_pos[j] = particle.dx; + y_pos[j] = particle.dy; + z_pos[j] = particle.dz; + + ux_pos[j] = particle.ux; + uy_pos[j] = particle.uy; + uz_pos[j] = particle.uz; + + std::array gi = global_particle_index(particle.i, grid, rank); + x_off[j] = (float)gi[1]; + y_off[j] = (float)gi[2]; + z_off[j] = (float)gi[3]; + } + + // Base offset plus i to account for chunks + auto o = openPMD::Offset{offset + i}; + auto e = openPMD::Extent{to_write}; + px.storeChunk(x_pos, o, e); + pxo.storeChunk(x_off, o, e); + + py.storeChunk(y_pos, o, e); + pyo.storeChunk(y_off, o, e); + + pz.storeChunk(z_pos, o, e); + pzo.storeChunk(z_off, o, e); + + ux.storeChunk(ux_pos, o, e); + uy.storeChunk(uy_pos, o, e); + uz.storeChunk(uz_pos, o, e); + + series.flush(); + } + } + void dump_hydro(const char *fbase, int step, hydro_array_t *hydro_array, + species_t *sp, interpolator_array_t *interpolator_array, + grid_t *grid, int ftag) { + std::string full_file_name = fbase + file_type; + + std::cout << "OpenPMD dumping hydro to " << full_file_name << std::endl; + + // if (series == nullptr) { + openPMD::Series series = openPMD::Series( + full_file_name, openPMD::AccessType::CREATE, MPI_COMM_WORLD); + //} + + auto i = series.iterations[step]; + + // TODO: set these + i.setTime((float)step); + i.setDt(1.0); + i.setTimeUnitSI(1.0); + + if (!sp) + ERROR(("Invalid species \"%s\"", sp->name)); + + // TODO: do we want each backend to have to explicitly call these + // manually? Or, as it is common, should we hoist it to the VPIC + // call-site + clear_hydro_array(hydro_array); + accumulate_hydro_p(hydro_array, sp, interpolator_array); + synchronize_hydro_array(hydro_array); + + if (!fbase) + ERROR(("Invalid filename")); + + if (rank == 0) + MESSAGE(("Dumping \"%s\" hydro fields to \"%s\"", sp->name, fbase)); + + // Write data + // float jx, jy, jz, rho; // Current and charge density => , + // float px, py, pz, ke; // Momentum and K.E. density => , float txx, tyy, tzz; // Stress diagonal => , i==j float tyz, tzx, txy; // Stress off-diagonal => , i!=j + auto J = i.meshes["J"]; + auto P = i.meshes["P"]; + auto T = i.meshes["T"]; + auto _Ke = i.meshes["Ke"]; + auto _Rho = i.meshes["Rho"]; + + auto Jx = J["x"]; + auto Jy = J["y"]; + auto Jz = J["z"]; + + auto Px = P["x"]; + auto Py = P["y"]; + auto Pz = P["z"]; + + auto Txx = T["xx"]; + auto Tyy = T["yy"]; + auto Tzz = T["zz"]; + auto Tyz = T["yz"]; + auto Tzx = T["zx"]; + auto Txy = T["xy"]; + + auto Rho = _Rho["rho"]; // TODO: bad name.. + auto Ke = _Ke["ke"]; // TODO: bad name.. + + size_t gnx = (grid->nx * grid->gpx); + size_t gny = (grid->ny * grid->gpy); + size_t gnz = (grid->nz * grid->gpz); + openPMD::Extent global_extent = {gnx, gny, gnz}; + + openPMD::Datatype datatype = openPMD::determineDatatype(); + openPMD::Dataset dataset = openPMD::Dataset(datatype, global_extent); + + Jx.resetDataset(dataset); + Jy.resetDataset(dataset); + Jz.resetDataset(dataset); + + Px.resetDataset(dataset); + Py.resetDataset(dataset); + Pz.resetDataset(dataset); + + Txx.resetDataset(dataset); + Tyy.resetDataset(dataset); + Tzz.resetDataset(dataset); + Tyz.resetDataset(dataset); + Tzx.resetDataset(dataset); + Txy.resetDataset(dataset); + + Rho.resetDataset(dataset); + Ke.resetDataset(dataset); + + // TODO: hoist this conversion code, as is it used elsewhere + // Convert rank to local x/y/z + int rx, ry, rz; + UNVOXEL(rank, rx, ry, rz, grid->gpx, grid->gpy, grid->gpz); + + size_t nx = grid->nx; + size_t ny = grid->ny; + size_t nz = grid->nz; + + // NOTE: this assumes a static mesh decomposition in nx/ny/nz + size_t global_offset_x = (nx)*rx; + size_t global_offset_y = (ny)*ry; + size_t global_offset_z = (nz)*rz; + + openPMD::Offset chunk_offset = {global_offset_x, global_offset_y, + global_offset_z}; + openPMD::Extent chunk_extent = {nx, ny, nz}; + + std::cout << "Local offset " + << " x: " << global_offset_x << " y: " << global_offset_y + << " z: " << global_offset_z << std::endl; + + std::vector jx_data; + std::vector jy_data; + std::vector jz_data; + + std::vector px_data; + std::vector py_data; + std::vector pz_data; + + std::vector txx_data; + std::vector tyy_data; + std::vector tzz_data; + std::vector tyz_data; + std::vector tzx_data; + std::vector txy_data; + + std::vector rho_data; + std::vector ke_data; + + size_t nv = nx * ny * nz; + + jx_data.resize(nv); + jy_data.resize(nv); + jz_data.resize(nv); + + px_data.resize(nv); + py_data.resize(nv); + pz_data.resize(nv); + + txx_data.resize(nv); + tyy_data.resize(nv); + tzz_data.resize(nv); + tyz_data.resize(nv); + tzx_data.resize(nv); + txy_data.resize(nv); + + rho_data.resize(nv); + ke_data.resize(nv); + + // Transpose AoS to SoAs + for (size_t k = 1; k < grid->nz + 1; k++) { + for (size_t j = 1; j < grid->ny + 1; j++) { + for (size_t i = 1; i < grid->nx + 1; i++) { + int local_index = VOXEL(i - 1, j - 1, k - 1, grid->nx - 2, + grid->ny - 2, grid->nz - 2); + int global_index = VOXEL(i, j, k, grid->nx, grid->ny, grid->nz); + + jx_data[local_index] = hydro_array->h[global_index].jx; + jy_data[local_index] = hydro_array->h[global_index].jy; + jz_data[local_index] = hydro_array->h[global_index].jz; + + px_data[local_index] = hydro_array->h[global_index].px; + py_data[local_index] = hydro_array->h[global_index].py; + pz_data[local_index] = hydro_array->h[global_index].pz; + + txx_data[local_index] = hydro_array->h[global_index].txx; + tyy_data[local_index] = hydro_array->h[global_index].tyy; + tzz_data[local_index] = hydro_array->h[global_index].tzz; + tyz_data[local_index] = hydro_array->h[global_index].tyz; + tzx_data[local_index] = hydro_array->h[global_index].tzx; + txy_data[local_index] = hydro_array->h[global_index].txy; + + rho_data[local_index] = hydro_array->h[global_index].rho; + ke_data[local_index] = hydro_array->h[global_index].ke; + } + } + } + + Jx.storeChunk(jx_data, chunk_offset, chunk_extent); + Jy.storeChunk(jy_data, chunk_offset, chunk_extent); + Jz.storeChunk(jz_data, chunk_offset, chunk_extent); + + Px.storeChunk(px_data, chunk_offset, chunk_extent); + Py.storeChunk(py_data, chunk_offset, chunk_extent); + Pz.storeChunk(pz_data, chunk_offset, chunk_extent); + + Txx.storeChunk(txx_data, chunk_offset, chunk_extent); + Tyy.storeChunk(tyy_data, chunk_offset, chunk_extent); + Tzz.storeChunk(tzz_data, chunk_offset, chunk_extent); + Tyz.storeChunk(tyz_data, chunk_offset, chunk_extent); + Tzx.storeChunk(tzx_data, chunk_offset, chunk_extent); + Txy.storeChunk(txy_data, chunk_offset, chunk_extent); + + Rho.storeChunk(rho_data, chunk_offset, chunk_extent); + Ke.storeChunk(ke_data, chunk_offset, chunk_extent); + + series.flush(); + } +}; +#endif + +/* + template + struct IODump : private Policy { + using Policy::dump_particles; + using Policy::dump_fields; + using Policy::dump_hydro; + }; + */ + +#endif diff --git a/src/vpic/dumpmacros.h b/src/vpic/dumpmacros.h index 9e46bf6b..bbb1b743 100644 --- a/src/vpic/dumpmacros.h +++ b/src/vpic/dumpmacros.h @@ -4,7 +4,7 @@ /* FIXME: WHEN THESE MACROS WERE HOISTED AND VARIOUS HACKS DONE TO THEM THEY BECAME _VERY_ _DANGEROUS. */ -#define WRITE_HEADER_V0(dump_type,sp_id,q_m,fileIO) do { \ +#define WRITE_HEADER_V0(dump_type,sp_id,q_m,fileIO,step,rank,nproc) do { \ /* Binary compatibility information */ \ WRITE( char, CHAR_BIT, fileIO ); \ WRITE( char, sizeof(short int), fileIO ); \ @@ -19,7 +19,7 @@ WRITE( int, 0 /* Version */, fileIO ); \ WRITE( int, dump_type, fileIO ); \ /* High level information */ \ - WRITE( int, step(), fileIO ); \ + WRITE( int, step, fileIO ); \ WRITE( int, nxout, fileIO ); \ WRITE( int, nyout, fileIO ); \ WRITE( int, nzout, fileIO ); \ @@ -33,21 +33,21 @@ WRITE( float, grid->cvac, fileIO ); \ WRITE( float, grid->eps0, fileIO ); \ WRITE( float, 0 /* damp */, fileIO ); \ - WRITE( int, rank(), fileIO ); \ - WRITE( int, nproc(), fileIO ); \ + WRITE( int, rank, fileIO ); \ + WRITE( int, nproc, fileIO ); \ /* Species parameters */ \ WRITE( int, sp_id, fileIO ); \ WRITE( float, q_m, fileIO ); \ } while(0) - + // Note dim _MUST_ be a pointer to an int - + #define WRITE_ARRAY_HEADER(p,ndim,dim,fileIO) do { \ WRITE( int, sizeof(p[0]), fileIO ); \ WRITE( int, ndim, fileIO ); \ fileIO.write( dim, ndim ); \ } while(0) - + // The WRITE macro copies the output "value" into a temporary variable // of the requested output "type" so that the write to the "file" // occurs from a known binary data type. For example, if grid.dx were @@ -60,12 +60,12 @@ // single precision write copies. However, specialty types could be // created so that the type cast __WRITE_tmp = (type)(value) // automatically does the underlying conversion in C++ - + #define WRITE(type,value,fileIO) do { \ type __WRITE_tmp = (type)(value); \ fileIO.write( &__WRITE_tmp, 1 ); \ } while(0) - + // Note: strlen does not include the terminating \0 #define WRITE_STRING(string,fileIO) do { \ int __WRITE_STRING_len = 0; \ @@ -74,103 +74,102 @@ if( __WRITE_STRING_len>0 ) \ fileIO.write( string, __WRITE_STRING_len ); \ } while(0) - + #define READ(type,value,fileIO) do { \ type __READ_tmp; \ fileIO.read(&__READ_tmp, 1 ); \ (value) = __READ_tmp; \ } while(0) -#define F_WRITE_HEADER_V0(dump_type,sp_id,q_m,fileIO) do { \ - /* Binary compatibility information */ \ - F_WRITE( char, CHAR_BIT, fileIO ); \ - F_WRITE( char, sizeof(short int), fileIO ); \ - F_WRITE( char, sizeof(int), fileIO ); \ - F_WRITE( char, sizeof(float), fileIO ); \ - F_WRITE( char, sizeof(double), fileIO ); \ - F_WRITE( short int, 0xcafe, fileIO ); \ - F_WRITE( int, 0xdeadbeef, fileIO ); \ - F_WRITE( float, 1.0, fileIO ); \ - F_WRITE( double, 1.0, fileIO ); \ - /* Dump type and header format version */ \ - F_WRITE( int, 0 /* Version */, fileIO ); \ - F_WRITE( int, dump_type, fileIO ); \ - /* High level information */ \ - F_WRITE( int, step(), fileIO ); \ - F_WRITE( int, imxstr-2, fileIO ); \ - F_WRITE( int, jmxstr-2, fileIO ); \ - F_WRITE( int, kmxstr-2, fileIO ); \ - F_WRITE( float, grid->dt, fileIO ); \ - F_WRITE( float, dxstr, fileIO ); \ - F_WRITE( float, dystr, fileIO ); \ - F_WRITE( float, dzstr, fileIO ); \ - F_WRITE( float, grid->x0, fileIO ); \ - F_WRITE( float, grid->y0, fileIO ); \ - F_WRITE( float, grid->z0, fileIO ); \ - F_WRITE( float, grid->cvac, fileIO ); \ - F_WRITE( float, grid->eps0, fileIO ); \ - F_WRITE( float, 0 /*damp*/, fileIO ); \ - F_WRITE( int, rank(), fileIO ); \ - F_WRITE( int, nproc(), fileIO ); \ - /* Species parameters */ \ - F_WRITE( int, sp_id, fileIO ); \ - F_WRITE( float, q_m, fileIO ); \ - } while(0) - -#define F_WRITE_HEADER_PAR(dump_type,sp_id,q_m,fileIO) do { \ - /* Binary compatibility information */ \ - F_WRITE( char, CHAR_BIT, fileIO ); \ - F_WRITE( char, sizeof(short int), fileIO ); \ - F_WRITE( char, sizeof(int), fileIO ); \ - F_WRITE( char, sizeof(float), fileIO ); \ - F_WRITE( char, sizeof(double), fileIO ); \ - F_WRITE( short int, 0xcafe, fileIO ); \ - F_WRITE( int, 0xdeadbeef, fileIO ); \ - F_WRITE( float, 1.0, fileIO ); \ - F_WRITE( double, 1.0, fileIO ); \ - /* Dump type and header format version */ \ - F_WRITE( int, 0 /* Version */, fileIO ); \ - F_WRITE( int, dump_type, fileIO ); \ - /* High level information */ \ - F_WRITE( int, step(), fileIO ); \ - F_WRITE( int, grid->nx, fileIO ); \ - F_WRITE( int, grid->ny, fileIO ); \ - F_WRITE( int, grid->nz, fileIO ); \ - F_WRITE( float, grid->dt, fileIO ); \ - F_WRITE( float, grid->dx, fileIO ); \ - F_WRITE( float, grid->dy, fileIO ); \ - F_WRITE( float, grid->dz, fileIO ); \ - F_WRITE( float, grid->x0, fileIO ); \ - F_WRITE( float, grid->y0, fileIO ); \ - F_WRITE( float, grid->z0, fileIO ); \ - F_WRITE( float, grid->cvac, fileIO ); \ - F_WRITE( float, grid->eps0, fileIO ); \ - F_WRITE( float, 0 /*damp*/, fileIO ); \ - F_WRITE( int, rank(), fileIO ); \ - F_WRITE( int, nproc(), fileIO ); \ - /* Species parameters */ \ - F_WRITE( int, sp_id, fileIO ); \ - F_WRITE( float, q_m, fileIO ); \ - } while(0) - +//#define F_WRITE_HEADER_V0(dump_type,sp_id,q_m,fileIO) do { \ + ///* Binary compatibility information */ \ + //F_WRITE( char, CHAR_BIT, fileIO ); \ + //F_WRITE( char, sizeof(short int), fileIO ); \ + //F_WRITE( char, sizeof(int), fileIO ); \ + //F_WRITE( char, sizeof(float), fileIO ); \ + //F_WRITE( char, sizeof(double), fileIO ); \ + //F_WRITE( short int, 0xcafe, fileIO ); \ + //F_WRITE( int, 0xdeadbeef, fileIO ); \ + //F_WRITE( float, 1.0, fileIO ); \ + //F_WRITE( double, 1.0, fileIO ); \ + ///* Dump type and header format version */ \ + //F_WRITE( int, 0 /* Version */, fileIO ); \ + //F_WRITE( int, dump_type, fileIO ); \ + ///* High level information */ \ + //F_WRITE( int, step(), fileIO ); \ + //F_WRITE( int, imxstr-2, fileIO ); \ + //F_WRITE( int, jmxstr-2, fileIO ); \ + //F_WRITE( int, kmxstr-2, fileIO ); \ + //F_WRITE( float, grid->dt, fileIO ); \ + //F_WRITE( float, dxstr, fileIO ); \ + //F_WRITE( float, dystr, fileIO ); \ + //F_WRITE( float, dzstr, fileIO ); \ + //F_WRITE( float, grid->x0, fileIO ); \ + //F_WRITE( float, grid->y0, fileIO ); \ + //F_WRITE( float, grid->z0, fileIO ); \ + //F_WRITE( float, grid->cvac, fileIO ); \ + //F_WRITE( float, grid->eps0, fileIO ); \ + //F_WRITE( float, 0 /*damp*/, fileIO ); \ + //F_WRITE( int, rank(), fileIO ); \ + //F_WRITE( int, nproc(), fileIO ); \ + ///* Species parameters */ \ + //F_WRITE( int, sp_id, fileIO ); \ + //F_WRITE( float, q_m, fileIO ); \ + //} while(0) + +//#define F_WRITE_HEADER_PAR(dump_type,sp_id,q_m,fileIO) do { \ + ///* Binary compatibility information */ \ + //F_WRITE( char, CHAR_BIT, fileIO ); \ + //F_WRITE( char, sizeof(short int), fileIO ); \ + //F_WRITE( char, sizeof(int), fileIO ); \ + //F_WRITE( char, sizeof(float), fileIO ); \ + //F_WRITE( char, sizeof(double), fileIO ); \ + //F_WRITE( short int, 0xcafe, fileIO ); \ + //F_WRITE( int, 0xdeadbeef, fileIO ); \ + //F_WRITE( float, 1.0, fileIO ); \ + //F_WRITE( double, 1.0, fileIO ); \ + ///* Dump type and header format version */ \ + //F_WRITE( int, 0 /* Version */, fileIO ); \ + //F_WRITE( int, dump_type, fileIO ); \ + ///* High level information */ \ + //F_WRITE( int, step(), fileIO ); \ + //F_WRITE( int, grid->nx, fileIO ); \ + //F_WRITE( int, grid->ny, fileIO ); \ + //F_WRITE( int, grid->nz, fileIO ); \ + //F_WRITE( float, grid->dt, fileIO ); \ + //F_WRITE( float, grid->dx, fileIO ); \ + //F_WRITE( float, grid->dy, fileIO ); \ + //F_WRITE( float, grid->dz, fileIO ); \ + //F_WRITE( float, grid->x0, fileIO ); \ + //F_WRITE( float, grid->y0, fileIO ); \ + //F_WRITE( float, grid->z0, fileIO ); \ + //F_WRITE( float, grid->cvac, fileIO ); \ + //F_WRITE( float, grid->eps0, fileIO ); \ + //F_WRITE( float, 0 /*damp*/, fileIO ); \ + //F_WRITE( int, rank(), fileIO ); \ + //F_WRITE( int, nproc(), fileIO ); \ + ///* Species parameters */ \ + //F_WRITE( int, sp_id, fileIO ); \ + //F_WRITE( float, q_m, fileIO ); \ + //} while(0) + // Note dim _MUST_ be a pointer to an int - -#define F_WRITE_ARRAY_HEADER(psiz,ndim,dim,fileIO) do { \ - F_WRITE( int, psiz, fileIO ); \ - F_WRITE( int, ndim, fileIO ); \ - fileIO.write( dim, ndim ); \ - } while(0) - -#define F_WRITE(type,value,fileIO) do { \ - type __F_WRITE_tmp = (type)(value); \ - fileIO.write( &__F_WRITE_tmp, 1 ); \ - } while(0) - -#define F_READ(type,value,fileIO) do { \ - type __F_READ_tmp; \ - fileIO.read( &__F_READ_tmp, 1 ); \ - (value) = __F_READ_tmp; \ - } while(0) +//#define F_WRITE_ARRAY_HEADER(psiz,ndim,dim,fileIO) do { \ + //F_WRITE( int, psiz, fileIO ); \ + //F_WRITE( int, ndim, fileIO ); \ + //fileIO.write( dim, ndim ); \ + //} while(0) + +//#define F_WRITE(type,value,fileIO) do { \ + //type __F_WRITE_tmp = (type)(value); \ + //fileIO.write( &__F_WRITE_tmp, 1 ); \ + //} while(0) + +//#define F_READ(type,value,fileIO) do { \ + //type __F_READ_tmp; \ + //fileIO.read( &__F_READ_tmp, 1 ); \ + //(value) = __F_READ_tmp; \ + //} while(0) #define ABORT(cond) if( cond ) ERROR(( #cond )) diff --git a/src/vpic/hdf5_header_info.h b/src/vpic/hdf5_header_info.h new file mode 100644 index 00000000..e3810612 --- /dev/null +++ b/src/vpic/hdf5_header_info.h @@ -0,0 +1,125 @@ +#ifndef VPIC_HDF5_HEAD_INFO +#define VPIC_HDF5_HEAD_INFO + +#define FIELD_ARRAY_NAME field_array + +namespace VPIC_HDF { + // XML header stuff + static const char *header = "\n\n\n\t\n"; + static const char *header_topology = "\t\t\n"; + static const char *header_geom = "\t\t\n"; + static const char *header_origin = "\t\t\t \n\t\t\t%s\n"; + static const char *header_dxdydz = "\t\t\t \n\t\t\t%s\n"; + static const char *footer_geom = "\t\t\n"; + static const char *grid_line = "\t\t \n \ + \t\t\t\n"; + static const char *footer = "\t\t\n\t\n\n"; + + static const char *main_body_head = "\t\t\t \n \ + \t\t\t\t \n \ + \t\t\t\t \n"; + static const char *main_body_foot = "\t\t\t\n"; + + static const char *main_body_attributeV = "\ + \t\t\t\t \n \ + \t\t\t\t\t \n \ + \t\t\t\t\t\t T.%d/%s_%d.h5:/Timestep_%d/%s \n \ + \t\t\t\t\t\t T.%d/%s_%d.h5:/Timestep_%d/%s \n \ + \t\t\t\t\t\t T.%d/%s_%d.h5:/Timestep_%d/%s \n \ + \t\t\t\t\t \n \ + \t\t\t\t \n "; + + static const char *main_body_attributeS = "\ + \t\t\t\t \n \ + \t\t\t\t\t\t T.%d/%s_%d.h5:/Timestep_%d/%s \n \ + \t\t\t\t \n "; + +} // end namespace + +#define create_file_with_header(xml_file_name, dimensions, orignal, dxdydz, nframes, fields_interval) \ + { \ + FILE *fp; \ + fp = fopen(xml_file_name, "w"); \ + fputs(VPIC_HDF::header, fp); \ + fprintf(fp, VPIC_HDF::header_topology, dimensions); \ + fputs(VPIC_HDF::header_geom, fp); \ + fprintf(fp, VPIC_HDF::header_origin, orignal); \ + fprintf(fp, VPIC_HDF::header_dxdydz, dxdydz); \ + fputs(VPIC_HDF::footer_geom, fp); \ + fprintf(fp, VPIC_HDF::grid_line, nframes); \ + int i; \ + for (i = 0; i < nframes; i++) \ + fprintf(fp, "%d ", i*fields_interval); \ + fputs(VPIC_HDF::grid_line_footer, fp); \ + fclose(fp); \ + } +#define write_main_body_attribute(fpp, main_body_attribute_p, attribute_name, dims_4d_p, dims_3d_p, file_name_pre_p, time_step_p, a1, a2, a3) \ + { \ + fprintf(fpp, main_body_attribute_p, attribute_name, dims_4d_p, \ + dims_3d_p, time_step_p, file_name_pre_p, time_step_p, time_step_p, a1, \ + dims_3d_p, time_step_p, file_name_pre_p, time_step_p, time_step_p, a2, \ + dims_3d_p, time_step_p, file_name_pre_p, time_step_p, time_step_p, a3); \ + } + +#define invert_field_xml_item(xml_file_name, speciesname_p, time_step, dims_4d, dims_3d, add_footer_flag) \ + { \ + FILE *fp; \ + fp = fopen(xml_file_name, "a"); \ + fprintf(fp, VPIC_HDF::main_body_head, time_step); \ + if (field_dump_flag.enabledE()) \ + write_main_body_attribute(fp, VPIC_HDF::main_body_attributeV, "E", dims_4d, dims_3d, speciesname_p, time_step, "ex", "ey", "ez"); \ + if (field_dump_flag.div_e_err) \ + fprintf(fp, VPIC_HDF::main_body_attributeS, "div_e_err", dims_3d, time_step, speciesname_p, time_step, time_step, "div_e_err"); \ + if (field_dump_flag.enabledCB()) \ + write_main_body_attribute(fp, VPIC_HDF::main_body_attributeV, "B", dims_4d, dims_3d, speciesname_p, time_step, "cbx", "cby", "cbz"); \ + if (field_dump_flag.div_b_err) \ + fprintf(fp, VPIC_HDF::main_body_attributeS, "div_b_err", dims_3d, time_step, speciesname_p, time_step, time_step, "div_b_err"); \ + if (field_dump_flag.enabledTCA()) \ + write_main_body_attribute(fp, VPIC_HDF::main_body_attributeV, "TCA", dims_4d, dims_3d, speciesname_p, time_step, "tcax", "tcay", "tcaz"); \ + if (field_dump_flag.rhob) \ + fprintf(fp, VPIC_HDF::main_body_attributeS, "rhob", dims_3d, time_step, speciesname_p, time_step, time_step, "rhob"); \ + if (field_dump_flag.enabledJF()) \ + write_main_body_attribute(fp, VPIC_HDF::main_body_attributeV, "JF", dims_4d, dims_3d, speciesname_p, time_step, "jfx", "jfy", "jfz"); \ + if (field_dump_flag.rhof) \ + fprintf(fp, VPIC_HDF::main_body_attributeS, "rhof", dims_3d, time_step, speciesname_p, time_step, time_step, "rhof"); \ + if (field_dump_flag.enabledEMAT()) \ + write_main_body_attribute(fp, VPIC_HDF::main_body_attributeV, "EMAT", dims_4d, dims_3d, speciesname_p, time_step, "ematx", "ematy", "ematz"); \ + if (field_dump_flag.nmat) \ + fprintf(fp, VPIC_HDF::main_body_attributeS, "nmat", dims_3d, time_step, speciesname_p, time_step, time_step, "nmat"); \ + if (field_dump_flag.enabledFMAT()) \ + write_main_body_attribute(fp, VPIC_HDF::main_body_attributeV, "FMAT", dims_4d, dims_3d, speciesname_p, time_step, "fmatx", "fmaty", "fmatz"); \ + if (field_dump_flag.cmat) \ + fprintf(fp, VPIC_HDF::main_body_attributeS, "cmat", dims_3d, time_step, speciesname_p, time_step, time_step, "cmat"); \ + fprintf(fp, "%s", VPIC_HDF::main_body_foot); \ + if (add_footer_flag) \ + fputs(VPIC_HDF::footer, fp); \ + fclose(fp); \ + } +#define invert_hydro_xml_item(xml_file_name, speciesname_p, time_step, dims_4d, dims_3d, add_footer_flag) \ + { \ + FILE *fp; \ + fp = fopen(xml_file_name, "a"); \ + fprintf(fp, VPIC_HDF::main_body_head, time_step); \ + if (hydro_dump_flag.enabledJ()) \ + write_main_body_attribute(fp, VPIC_HDF::main_body_attributeV, "J", dims_4d, dims_3d, speciesname_p, time_step, "jx", "jy", "jz"); \ + if (hydro_dump_flag.rho) \ + fprintf(fp, VPIC_HDF::main_body_attributeS, "rho", dims_3d, time_step, speciesname_p, time_step, time_step, "rho"); \ + if (hydro_dump_flag.enabledP()) \ + write_main_body_attribute(fp, VPIC_HDF::main_body_attributeV, "P", dims_4d, dims_3d, speciesname_p, time_step, "px", "py", "pz"); \ + if (hydro_dump_flag.ke) \ + fprintf(fp, VPIC_HDF::main_body_attributeS, "ke", dims_3d, time_step, speciesname_p, time_step, time_step, "ke"); \ + if (hydro_dump_flag.enabledTD()) \ + write_main_body_attribute(fp, VPIC_HDF::main_body_attributeV, "TD", dims_4d, dims_3d, speciesname_p, time_step, "txx", "tyy", "tzz"); \ + if (hydro_dump_flag.enabledTOD()) \ + write_main_body_attribute(fp, VPIC_HDF::main_body_attributeV, "TOD", dims_4d, dims_3d, speciesname_p, time_step, "tyz", "tzx", "txy"); \ + fprintf(fp, "%s", VPIC_HDF::main_body_foot); \ + if (add_footer_flag) { \ + fputs(VPIC_HDF::footer, fp); \ + } \ + fclose(fp); \ + } + +#endif // VPIC_HDF5_HEAD_INFO diff --git a/src/vpic/vpic.cc b/src/vpic/vpic.cc index 4f150afc..cc081358 100644 --- a/src/vpic/vpic.cc +++ b/src/vpic/vpic.cc @@ -7,9 +7,10 @@ * March/April 2004 - Heavily revised and extended from earlier V4PIC versions * */ - + #include "vpic.h" - +#include "dump_strategy.h" + /* Note that, when a vpic_simulation is created (and thus registered with the checkpt service), it is created empty; none of the simulation objects on which it depends have been created yet. (These get created @@ -71,14 +72,18 @@ reanimate_vpic_simulation( vpic_simulation * vpic ) { } -vpic_simulation::vpic_simulation() { - CLEAR( this, 1 ); +vpic_simulation::vpic_simulation() +{ + // TODO: why is this a good idea? + // Is this just trying to 0 initialize everything? + // CLEAR( this, 1 ); - /* Set non-zero defaults */ - verbose = 1; - num_comm_round = 3; - num_div_e_round = 2; - num_div_b_round = 2; + // Now done in the class def / header + ///* Set non-zero defaults */ + //verbose = 1; + //num_comm_round = 3; + //num_div_e_round = 2; + //num_div_b_round = 2; #if defined(VPIC_USE_PTHREADS) // Pthreads case. int n_rng = serial.n_pipeline; @@ -99,7 +104,7 @@ vpic_simulation::vpic_simulation() { // if( n_rng(new BinaryDump( rank(), nproc() )); + enable_binary_dump(); + + // TODO: this this still makes sense now we have a dump strategy +//#ifdef VPIC_ENABLE_HDF5 + // Default init hdf5 dump flags + //field_interval = 1; + //hydro_interval = 1; + //field_dump_flag = field_dump_flag_t(); + //hydro_dump_flag = hydro_dump_flag_t(); +//#endif + + field_interval = 1; + hydro_interval = 1; } - + vpic_simulation::~vpic_simulation() { UNREGISTER_OBJECT( this ); delete_emitter_list( emitter_list ); @@ -123,4 +145,4 @@ vpic_simulation::~vpic_simulation() { delete_rng_pool( sync_entropy ); delete_rng_pool( entropy ); } - + diff --git a/src/vpic/vpic.h b/src/vpic/vpic.h index 9c22048e..2a4b1767 100644 --- a/src/vpic/vpic.h +++ b/src/vpic/vpic.h @@ -15,6 +15,7 @@ #include #include +#include // unique_ptr #include "../boundary/boundary.h" #include "../collision/collision.h" @@ -24,6 +25,7 @@ #include "../util/bitfield.h" #include "../util/checksum.h" #include "../util/system.h" +#include "dump_strategy.h" #ifndef USER_GLOBAL_SIZE #define USER_GLOBAL_SIZE 16384 @@ -121,56 +123,85 @@ class vpic_simulation { public: vpic_simulation(); ~vpic_simulation(); + vpic_simulation(const vpic_simulation&) = delete; + vpic_simulation& operator=(const vpic_simulation&) = delete; + void initialize( int argc, char **argv ); void modify( const char *fname ); int advance( void ); void finalize( void ); + // TODO: decide if I should collapse this to an enum + // An enum would stop these ifdefs being so leaky + void enable_binary_dump(); +#ifdef VPIC_ENABLE_HDF5 + void enable_hdf5_dump(); +#endif +#ifdef VPIC_ENABLE_OPENPMD + void enable_openpmd_dump(); +#endif + + // TODO: remake these protected + + // Very likely a user will forgot to delete this if they change the strategy, + // a smart ptr will save us from the small leak + // TODO: this does not survive the dump right now + std::unique_ptr dump_strategy; + + int num_step = 1; // Number of steps to take + protected: // Directly initialized by user - int verbose; // Should system be verbose - int num_step; // Number of steps to take - int num_comm_round; // Num comm round - int status_interval; // How often to print status messages - int clean_div_e_interval; // How often to clean div e - int num_div_e_round; // How many clean div e rounds per div e interval - int clean_div_b_interval; // How often to clean div b - int num_div_b_round; // How many clean div b rounds per div b interval - int sync_shared_interval; // How often to synchronize shared faces + int verbose = 1; // Should system be verbose + int num_comm_round = 3; // Num comm round + int status_interval = 0; // How often to print status messages + + int clean_div_e_interval = 0; // How often to clean div e + int num_div_e_round = 2; // How many clean div e rounds per div e interval + + int clean_div_b_interval = 0; // How often to clean div b + int num_div_b_round = 2; // How many clean div b rounds per div b interval + + int sync_shared_interval = 0; // How often to synchronize shared faces // FIXME: THESE INTERVALS SHOULDN'T BE PART OF vpic_simulation // THE BIG LIST FOLLOWING IT SHOULD BE CLEANED UP TOO - double quota; - int checkpt_interval; - int hydro_interval; - int field_interval; - int particle_interval; - - size_t nxout, nyout, nzout; - size_t px, py, pz; - float dxout, dyout, dzout; - - int ndfld; - int ndhyd; - int ndpar; - int ndhis; - int ndgrd; - int head_option; - int istride; - int jstride; - int kstride; - int stride_option; - int pstride; - int nprobe; + double quota = 0; + int checkpt_interval = 0; + int hydro_interval = 0; + int field_interval = 0; + int particle_interval = 0; + + // TODO: these can probably now be removed, as they should only be used by dump? + // TODO: check if any decks used them + //size_t nxout, nyout, nzout; + //float dxout, dyout, dzout; + + size_t px = 0; + size_t py = 0; + size_t pz = 0; + + int ndfld = 0; + int ndhyd = 0; + int ndpar = 0; + int ndhis = 0; + int ndgrd = 0; + int head_option = 0; + int istride = 0; + int jstride = 0; + int kstride = 0; + int stride_option = 0; + int pstride = 0; + int nprobe = 0; int ijkprobe[NVARHISMX][4]; float xyzprobe[NVARHISMX][3]; - int block_dump; - int stepdigit; - int rankdigit; - int ifenergies; + int block_dump = 0; + int stepdigit = 0; + int rankdigit = 0; + int ifenergies = 0; // Helper initialized by user @@ -180,21 +211,21 @@ class vpic_simulation { random numbers. Keeping the synchronous generators in sync is the generator users responsibility. */ - rng_pool_t * entropy; // Local entropy pool - rng_pool_t * sync_entropy; // Synchronous entropy pool - grid_t * grid; // define_*_grid et al - material_t * material_list; // define_material - field_array_t * field_array; // define_field_array - interpolator_array_t * interpolator_array; // define_interpolator_array - accumulator_array_t * accumulator_array; // define_accumulator_array - hydro_array_t * hydro_array; // define_hydro_array - species_t * species_list; // define_species / - // species helpers - particle_bc_t * particle_bc_list; // define_particle_bc / - // boundary helpers - emitter_t * emitter_list; // define_emitter / - // emitter helpers - collision_op_t * collision_op_list; // collision helpers + rng_pool_t * entropy = nullptr; // Local entropy pool + rng_pool_t * sync_entropy = nullptr; // Synchronous entropy pool + grid_t * grid = nullptr; // define_*_grid et al + material_t * material_list = nullptr; // define_material + field_array_t * field_array = nullptr; // define_field_array + interpolator_array_t * interpolator_array = nullptr; // define_interpolator_array + accumulator_array_t * accumulator_array = nullptr; // define_accumulator_array + hydro_array_t * hydro_array = nullptr; // define_hydro_array + species_t * species_list = nullptr; // define_species / + // species helpers + particle_bc_t * particle_bc_list = nullptr; // define_particle_bc / + // boundary helpers + emitter_t * emitter_list = nullptr; // define_emitter / + // emitter helpers + collision_op_t * collision_op_list = nullptr; // collision helpers // User defined checkpt preserved variables // Note: user_global is aliased with user_global_t (see deck_wrapper.cxx) @@ -223,7 +254,7 @@ class vpic_simulation { /////////////// // Dump helpers - int dump_mkdir(const char * dname); + static int dump_mkdir(const char * dname); int dump_cwd(char * dname, size_t size); // Text dumps @@ -233,12 +264,16 @@ class vpic_simulation { // Binary dumps void dump_grid( const char *fbase ); + void dump_fields( const char *fbase, int fname_tag = 1 ); + void dump_hydro( const char *sp_name, const char *fbase, int fname_tag = 1 ); + void dump_particles( const char *sp_name, const char *fbase, int fname_tag = 1 ); + // convenience functions for simlog output void create_field_list(char * strlist, DumpParameters & dumpParams); void create_hydro_list(char * strlist, DumpParameters & dumpParams); diff --git a/test/integrated/to_completion/CMakeLists.txt b/test/integrated/to_completion/CMakeLists.txt index 5baecf43..690e738c 100644 --- a/test/integrated/to_completion/CMakeLists.txt +++ b/test/integrated/to_completion/CMakeLists.txt @@ -58,11 +58,11 @@ add_test(${generate_restore} ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} # Run using the restore file # TODO: caps? -set(perform_restore "perform_${RESTART_BINARY}") -build_a_vpic(${perform_restore} ${CMAKE_CURRENT_SOURCE_DIR}/${RESTART_DECK}.deck) -add_test(${perform_restore} ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} - ${MPIEXEC_NUMPROC} ${MPIEXEC_PREFLAGS} ${perform_restore} - ${MPIEXEC_POSTFLAGS} ${RESTART_ARGS}) +#set(perform_restore "perform_${RESTART_BINARY}") +#build_a_vpic(${perform_restore} ${CMAKE_CURRENT_SOURCE_DIR}/${RESTART_DECK}.deck) +#add_test(${perform_restore} ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} +#${MPIEXEC_NUMPROC} ${MPIEXEC_PREFLAGS} ${perform_restore} +#${MPIEXEC_POSTFLAGS} ${RESTART_ARGS}) # TODO: re-enable modify test #list(APPEND MODIFY_BINARY restore-modify) @@ -76,4 +76,4 @@ add_test(${perform_restore} ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} set(RESTORE_LABEL "restore_group") set_tests_properties(${perform_restore} PROPERTIES DEPENDS ${generate_restore}) #set_property(TEST ${generate_restore} PROPERTY FIXTURES_SETUP ${RESTORE_LABEL}) -#set_property(TEST ${perform_restore} PROPERTY FIXTURES_REQUIRED ${RESTORE_LABEL}) \ No newline at end of file +#set_property(TEST ${perform_restore} PROPERTY FIXTURES_REQUIRED ${RESTORE_LABEL}) diff --git a/test/unit/energy_comparison/3d_test.cc b/test/unit/energy_comparison/3d_test.cc index 4c44736b..d84e4627 100644 --- a/test/unit/energy_comparison/3d_test.cc +++ b/test/unit/energy_comparison/3d_test.cc @@ -312,7 +312,7 @@ TEST_CASE( "Check if Weibel gives correct energy (within tol)", "[energy]" ) ofs.close(); // Init and run sim - vpic_simulation simulation = vpic_simulation(); + vpic_simulation simulation; // TODO: We should do this in a safer manner simulation.initialize( 0, NULL ); diff --git a/test/unit/energy_comparison/weibel_driver.cc b/test/unit/energy_comparison/weibel_driver.cc index eb4702db..cbc64e2c 100644 --- a/test/unit/energy_comparison/weibel_driver.cc +++ b/test/unit/energy_comparison/weibel_driver.cc @@ -310,7 +310,7 @@ TEST_CASE( "Check if Weibel gives correct energy (within tol)", "[energy]" ) ofs.close(); // Init and run sim - vpic_simulation simulation = vpic_simulation(); + vpic_simulation simulation; // TODO: We should do this in a safer manner simulation.initialize( 0, NULL ); diff --git a/test/unit/grid_heating/gridHeatingTestElec.cxx b/test/unit/grid_heating/gridHeatingTestElec.cxx index dc0e1ca5..9e804588 100644 --- a/test/unit/grid_heating/gridHeatingTestElec.cxx +++ b/test/unit/grid_heating/gridHeatingTestElec.cxx @@ -249,7 +249,7 @@ begin_initialization { TEST_CASE( "Check if Weibel gives correct energy (within tol)", "[energy]" ) { // Init and run sim - vpic_simulation simulation = vpic_simulation(); + vpic_simulation simulation; // TODO: We should do this in a safer manner simulation.initialize( 0, NULL );