diff --git a/CMakeLists.txt b/CMakeLists.txt index 2724a9ff..2c912f52 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -84,10 +84,12 @@ 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 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") @@ -125,31 +127,18 @@ 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) - 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 +153,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 +277,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 +309,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 +338,45 @@ 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) + macro(build_a_vpic name deck) if(NOT EXISTS ${deck}) message(FATAL_ERROR "Could not find deck '${deck}'") @@ -438,6 +434,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/sample/harrisHDF5 b/sample/harrisHDF5 new file mode 100644 index 00000000..2b3b21bf --- /dev/null +++ b/sample/harrisHDF5 @@ -0,0 +1,445 @@ +// 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 = 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. + if( step()==-10 ) dump_fields_hdf5("fields"); // Get first valid total J + if( should_dump(fields) ) dump_fields_hdf5("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_hdf5("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_hdf5("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/harrisOpenPMD b/sample/harrisOpenPMD new file mode 100644 index 00000000..16abb432 --- /dev/null +++ b/sample/harrisOpenPMD @@ -0,0 +1,446 @@ +// 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 { + // 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. + + std::string openpm_field_name = "fields.h5"; + //std::string openpm_field_name = "fields.bp"; + if( step()==-10 ) dump_fields_openpmd(openpm_field_name.c_str()); // Get first valid total J + if( should_dump(fields) ) dump_fields_openpmd(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_hdf5("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_openpmd("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/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/vpic/dump.cc b/src/vpic/dump.cc index 62505147..61757331 100644 --- a/src/vpic/dump.cc +++ b/src/vpic/dump.cc @@ -10,11 +10,21 @@ */ #include +#include #include "vpic.h" #include "dumpmacros.h" #include "../util/io/FileUtils.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 + /* -1 means no ranks talk */ #define VERBOSE_rank -1 @@ -22,6 +32,10 @@ // COMPATIBLE WITH EXISTING EXTERNAL 3RD PARTY VISUALIZATION SOFTWARE. // IN THE LONG RUN, THIS EXTERNAL SOFTWARE WILL NEED TO BE UPDATED. +// TODO: this should live somewhere more sensible, but it's better than the +// global static it replaces +std::unordered_map tframe_map; + int vpic_simulation::dump_mkdir(const char * dname) { return FileUtils::makeDirectory(dname); } // dump_mkdir @@ -30,6 +44,39 @@ int vpic_simulation::dump_cwd(char * dname, size_t size) { return FileUtils::getCurrentWorkingDirectory(dname, size); } // dump_mkdir + +// TODO: move this somewhere more sensible +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 }; +} + + /***************************************************************************** * ASCII dump IO *****************************************************************************/ @@ -256,6 +303,1131 @@ vpic_simulation::dump_hydro( const char *sp_name, if( fileIO.close() ) ERROR(( "File close failed on dump hydro!!!" )); } +#ifdef VPIC_ENABLE_OPENPMD + +// TODO: remove this hack, and actually store the state properly +static openPMD::Series* series; + +void +vpic_simulation::dump_particles_openpmd( const char *sp_name, + const char *fbase, + int ftag ) +{ + + species_t *sp = find_species_name( sp_name, species_list ); + + if (series == nullptr) { + std::cout << "init series" << std::endl; + series = new openPMD::Series( + fbase, + 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]; + //openPMD::ParticleSpecies& 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"]; + + px.resetDataset(dataset); + pxo.resetDataset(dataset); + + // convert data to SoA, allowing the user to chunk the operation + const int max_chunk = 32768*8; // 1MB SoA + // 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.reserve(to_write); + x_off.reserve(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; + std::array gi = global_particle_index(particle.i, grid, rank()); + x_off[j] = (float)gi[1]; + } + + // 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); + } + + +} + +void vpic_simulation::dump_fields_openpmd(const char *fbase, int ftag) +{ + std::cout << "Writing openPMD data" << std::endl; + + if (series == nullptr) { + std::cout << "init series" << std::endl; + series = new openPMD::Series( + fbase, + openPMD::AccessType::CREATE, + MPI_COMM_WORLD + ); + } + + std::cout << "Writing itration " << 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"]; + + // 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"]; + + // TODO: set unitDimension so the anaylsis software knows what fields + // things are + + 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 = {gny, 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); + + // 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}; + + // 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; + + size_t nv = nx * ny * nz; + + cbx_data.reserve(nv); + cby_data.reserve(nv); + cbz_data.reserve(nv); + + ex_data.reserve(nv); + ey_data.reserve(nv); + ez_data.reserve(nv); + + jx_data.reserve(nv); + jy_data.reserve(nv); + jz_data.reserve(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; + } + } + } + + 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); + + series->flush(); +} +#endif + +#ifdef VPIC_ENABLE_HDF5 +#define DUMP_DIR_FORMAT "./%s" + +/* define to do C-style indexing */ +#define hydro(x, y, z) hydro_array->h[VOXEL(x, y, z, grid->nx, grid->ny, grid->nz)] + +// TODO: make function? +#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, main_body_head, time_step); \ + if (field_dump_flag.enabledE()) \ + write_main_body_attribute(fp, main_body_attributeV, "E", dims_4d, dims_3d, speciesname_p, time_step, "ex", "ey", "ez"); \ + if (field_dump_flag.div_e_err) \ + fprintf(fp, 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, main_body_attributeV, "B", dims_4d, dims_3d, speciesname_p, time_step, "cbx", "cby", "cbz"); \ + if (field_dump_flag.div_b_err) \ + fprintf(fp, 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, main_body_attributeV, "TCA", dims_4d, dims_3d, speciesname_p, time_step, "tcax", "tcay", "tcaz"); \ + if (field_dump_flag.rhob) \ + fprintf(fp, 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, main_body_attributeV, "JF", dims_4d, dims_3d, speciesname_p, time_step, "jfx", "jfy", "jfz"); \ + if (field_dump_flag.rhof) \ + fprintf(fp, 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, main_body_attributeV, "EMAT", dims_4d, dims_3d, speciesname_p, time_step, "ematx", "ematy", "ematz"); \ + if (field_dump_flag.nmat) \ + fprintf(fp, 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, main_body_attributeV, "FMAT", dims_4d, dims_3d, speciesname_p, time_step, "fmatx", "fmaty", "fmatz"); \ + if (field_dump_flag.cmat) \ + fprintf(fp, main_body_attributeS, "cmat", dims_3d, time_step, speciesname_p, time_step, time_step, "cmat"); \ + fprintf(fp, "%s", main_body_foot); \ + if (add_footer_flag) \ + fputs(footer, fp); \ + fclose(fp); \ + } +void +vpic_simulation::dump_fields_hdf5( const char *fbase, 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 + +#define DUMP_FIELD_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] = 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); \ + } + + char fname[256]; + char field_scratch[128]; + char subfield_scratch[128]; + + sprintf(field_scratch, DUMP_DIR_FORMAT, "field_hdf5"); + dump_mkdir(field_scratch); + sprintf(subfield_scratch, "%s/T.%zu/", field_scratch, step_for_viou); + dump_mkdir(subfield_scratch); + + sprintf(fname, "%s/%s_%zu.h5", subfield_scratch, "fields", step_for_viou); + double el1 = uptime(); + hid_t plist_id = H5Pcreate(H5P_FILE_ACCESS); + H5Pset_fapl_mpio(plist_id, MPI_COMM_WORLD, MPI_INFO_NULL); + hid_t 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 = H5Gcreate(file_id, fname, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + + el1 = uptime() - el1; + //sim_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"}; + 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 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); + + 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 + printf("global size = %d %d %d \n", field_global_size[0], field_global_size[1], field_global_size[2]); + printf("global_offset = %d %d %d \n", global_offset[0], global_offset[1], global_offset[2]); + printf("global_count = %d %d %d \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, field_global_size, NULL); + hid_t memspace = H5Screate_simple(3, field_local_size, NULL); + hid_t dataspace_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;*/ + + 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); + + el2 = uptime() - el2; + //sim_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; + 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; + //sim_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); + + 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: %lld \n", step_for_viou); + printf(" current step: %lld \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++; + } +} + +// TODO: fix this, it currently uses a static global and the logic only +// supports 1 species otherwise things get out of sync +void vpic_simulation::dump_hydro_hdf5( const char *speciesname, + const char *fbase, + 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); + + species_t *sp = find_species_name(speciesname, species_list); + if (!sp) + ERROR(("Invalid species name: %s", speciesname)); + + 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"); + dump_mkdir(hydro_scratch); + sprintf(subhydro_scratch, "%s/T.%zu/", hydro_scratch, step_for_viou); + dump_mkdir(subhydro_scratch); + + sprintf(hname, "%s/hydro_%s_%zu.h5", subhydro_scratch, speciesname, step_for_viou); + double el1 = uptime(); + hid_t plist_id = H5Pcreate(H5P_FILE_ACCESS); + 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; + //sim_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; + + //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; + RANK_TO_INDEX2(mpi_rank, mpi_rank_x, mpi_rank_y, mpi_rank_z); + + 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 = %d %d %d \n", hydro_global_size[0], hydro_global_size[1], hydro_global_size[2]); + printf("global_offset = %d %d %d \n", global_offset[0], global_offset[1], global_offset[2]); + printf("global_count = %d %d %d \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; + + 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; + //sim_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; + 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; + //sim_log("TimeHDF5Close: " << el3 << " s"); + + if (mpi_rank == 0) + { + char output_xml_file[128]; + sprintf(output_xml_file, "./%s/%s%s%s", "hydro_hdf5", "hydro-", speciesname, ".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); + + 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: %lld \n", step_for_viou); + printf(" Simulation time: %f \n", grid->t0); + printf(" tframe: %d \n", tframe); +#endif + + char speciesname_new[128]; + sprintf(speciesname_new, "hydro_%s", speciesname); + 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]++; + } +} + +// TODO": make the sp_name and speciesname varailbe naming consistent +void +vpic_simulation::dump_particles_hdf5( const char *sp_name, + const char *fbase, + int ftag ) +{ + size_t step_for_viou = step(); + char fname[256]; + char group_name[256]; + char particle_scratch[128]; + char subparticle_scratch[128]; + + int np_local; + species_t *sp; + + 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"); + dump_mkdir(particle_scratch); + sprintf(subparticle_scratch, "%s/T.%ld/", particle_scratch, step_for_viou); + dump_mkdir(subparticle_scratch); + + // TODO: Allow the user to set this + + int stride_particle_dump = 1; + while (sp) + { + 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; + int mpi_rank; + MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank); + + //std::cout << "on mpi_rank: " << mpi_rank << ", time in copying particle data: " << ec1 << " s" << ", np_local = " << np_local << std::endl; + //sim_log("on mpi_rank: " << mpi_rank << ", time in copying particle data: " << ec1 << " s" << ", np_local = " << np_local); + + Pf = (float *)sp->p; + Pi = (int *)sp->p; + + // open HDF5 file in "particle/T./" subdirectory + // filename: eparticle.h5p + sprintf(fname, "%s/%s_%ld.h5", subparticle_scratch, sp->name, step_for_viou); + sprintf(group_name, "/Timestep_%ld", step_for_viou); + double el1 = uptime(); + + hid_t plist_id = H5Pcreate(H5P_FILE_ACCESS); + H5Pset_fapl_mpio(plist_id, MPI_COMM_WORLD, MPI_INFO_NULL); + hid_t file_id = H5Fcreate(fname, H5F_ACC_TRUNC, H5P_DEFAULT, plist_id); + hid_t group_id = H5Gcreate(file_id, group_name, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + + H5Pclose(plist_id); + + long long total_particles, offset; + long long numparticles = np_local; + 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; + + hid_t filespace = H5Screate_simple(1, (hsize_t *)&total_particles, NULL); + + hsize_t memspace_count_temp = numparticles * 8; + hid_t memspace = H5Screate_simple(1, &memspace_count_temp, NULL); + + // The converted global_ids are stored compact, not strided + hsize_t linearspace_count_temp = numparticles; + hid_t linearspace = H5Screate_simple(1, &linearspace_count_temp, NULL); + + 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); + + 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); + + el1 = uptime() - el1; + //sim_log("Particle TimeHDF5Open): " << el1 << " s"); //Easy to handle results for scripts + + //double el2 = uptime(); + + // This point offset is silly, and loses the type safety (pf+1) + hid_t dset_id = H5Dcreate(group_id, "dX", H5T_NATIVE_FLOAT, filespace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + int ierr = H5Dwrite(dset_id, H5T_NATIVE_FLOAT, memspace, filespace, plist_id, Pf); + H5Dclose(dset_id); + + dset_id = H5Dcreate(group_id, "dY", H5T_NATIVE_FLOAT, filespace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + ierr = H5Dwrite(dset_id, H5T_NATIVE_FLOAT, memspace, filespace, plist_id, Pf + 1); + H5Dclose(dset_id); + + dset_id = H5Dcreate(group_id, "dZ", H5T_NATIVE_FLOAT, filespace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + ierr = H5Dwrite(dset_id, H5T_NATIVE_FLOAT, memspace, filespace, plist_id, Pf + 2); + H5Dclose(dset_id); + +#define OUTPUT_CONVERT_GLOBAL_ID 1 +#ifdef OUTPUT_CONVERT_GLOBAL_ID + // TODO: make a function out of this too, its used in openpmd + std::vector global_pi; + global_pi.reserve(numparticles); + // TODO: this could be parallel + for (int i = 0; i < numparticles; i++) + { + int local_i = sp->p[i].i; + + 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); + + //std::cout << rank() << " local i " << local_i << " becomes " << global_i << std::endl; + global_pi[i] = global_i; + } + +#undef UNVOXEL + dset_id = H5Dcreate(group_id, "i", H5T_NATIVE_INT, filespace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + ierr = H5Dwrite(dset_id, H5T_NATIVE_INT, linearspace, filespace, plist_id, global_pi.data()); + H5Dclose(dset_id); + +#else + dset_id = H5Dcreate(group_id, "i", H5T_NATIVE_INT, filespace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + ierr = H5Dwrite(dset_id, H5T_NATIVE_INT, memspace, filespace, plist_id, Pi + 3); + H5Dclose(dset_id); +#endif + + dset_id = H5Dcreate(group_id, "Ux", H5T_NATIVE_FLOAT, filespace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + ierr = H5Dwrite(dset_id, H5T_NATIVE_FLOAT, memspace, filespace, plist_id, Pf + 4); + H5Dclose(dset_id); + + dset_id = H5Dcreate(group_id, "Uy", H5T_NATIVE_FLOAT, filespace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + ierr = H5Dwrite(dset_id, H5T_NATIVE_FLOAT, memspace, filespace, plist_id, Pf + 5); + H5Dclose(dset_id); + + dset_id = H5Dcreate(group_id, "Uz", H5T_NATIVE_FLOAT, filespace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + ierr = H5Dwrite(dset_id, H5T_NATIVE_FLOAT, memspace, filespace, plist_id, Pf + 6); + H5Dclose(dset_id); + + dset_id = H5Dcreate(group_id, "q", H5T_NATIVE_FLOAT, filespace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + ierr = H5Dwrite(dset_id, H5T_NATIVE_FLOAT, memspace, filespace, plist_id, Pf + 7); + H5Dclose(dset_id); + + //el2 = uptime() - el2; + //sim_log("Particle TimeHDF5Write: " << el2 << " s"); + + double el3 = uptime(); + H5Sclose(memspace); + H5Sclose(filespace); + H5Pclose(plist_id); + H5Gclose(group_id); + H5Fclose(file_id); + el3 = uptime() - el3; + //sim_log("Particle TimeHDF5Close: " << el3 << " s"); + + sp->p = sp_p; + sp->np = sp_np; + sp->max_np = sp_max_np; + FREE_ALIGNED(p_buf); + + // Write metadata if step() == 0 + char meta_fname[256]; + + sprintf(meta_fname, "%s/grid_metadata_%s_%ld.h5", subparticle_scratch, sp->name, step_for_viou); + + double meta_el1 = uptime(); + + hid_t meta_plist_id = H5Pcreate(H5P_FILE_ACCESS); + H5Pset_fapl_mpio(meta_plist_id, MPI_COMM_WORLD, MPI_INFO_NULL); + hid_t meta_file_id = H5Fcreate(meta_fname, H5F_ACC_TRUNC, H5P_DEFAULT, meta_plist_id); + hid_t meta_group_id = H5Gcreate(meta_file_id, group_name, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + H5Pclose(meta_plist_id); + + long long meta_total_particles, meta_offset; + long long meta_numparticles = 1; + MPI_Allreduce(&meta_numparticles, &meta_total_particles, 1, MPI_LONG_LONG, MPI_SUM, MPI_COMM_WORLD); + MPI_Scan(&meta_numparticles, &meta_offset, 1, MPI_LONG_LONG, MPI_SUM, MPI_COMM_WORLD); + meta_offset -= meta_numparticles; + + hid_t meta_filespace = H5Screate_simple(1, (hsize_t *)&meta_total_particles, NULL); + hid_t meta_memspace = H5Screate_simple(1, (hsize_t *)&meta_numparticles, NULL); + meta_plist_id = H5Pcreate(H5P_DATASET_XFER); + H5Pset_dxpl_mpio(meta_plist_id, H5FD_MPIO_COLLECTIVE); + H5Sselect_hyperslab(meta_filespace, H5S_SELECT_SET, (hsize_t *)&meta_offset, NULL, (hsize_t *)&meta_numparticles, NULL); + meta_el1 = uptime() - meta_el1; + //sim_log("Metafile TimeHDF5Open): " << meta_el1 << " s"); //Easy to handle results for scripts + + double meta_el2 = uptime(); + + hid_t meta_dset_id = H5Dcreate(meta_group_id, "np_local", H5T_NATIVE_INT, meta_filespace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + ierr = H5Dwrite(meta_dset_id, H5T_NATIVE_INT, meta_memspace, meta_filespace, meta_plist_id, (int32_t *)&np_local); + H5Dclose(meta_dset_id); + //if (rank == 0) printf ("Written variable dX \n"); + + meta_dset_id = H5Dcreate(meta_group_id, "nx", H5T_NATIVE_INT, meta_filespace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + ierr = H5Dwrite(meta_dset_id, H5T_NATIVE_INT, meta_memspace, meta_filespace, meta_plist_id, &grid->nx); + H5Dclose(meta_dset_id); + //if (rank == 0) printf ("Written variable dY \n"); + + meta_dset_id = H5Dcreate(meta_group_id, "ny", H5T_NATIVE_INT, meta_filespace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + ierr = H5Dwrite(meta_dset_id, H5T_NATIVE_INT, meta_memspace, meta_filespace, meta_plist_id, &grid->ny); + H5Dclose(meta_dset_id); + //if (rank == 0) printf ("Written variable dZ \n"); + + meta_dset_id = H5Dcreate(meta_group_id, "nz", H5T_NATIVE_INT, meta_filespace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + ierr = H5Dwrite(meta_dset_id, H5T_NATIVE_INT, meta_memspace, meta_filespace, meta_plist_id, &grid->nz); + H5Dclose(meta_dset_id); + //if (rank == 0) printf ("Written variable i \n"); + + meta_dset_id = H5Dcreate(meta_group_id, "x0", H5T_NATIVE_FLOAT, meta_filespace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + ierr = H5Dwrite(meta_dset_id, H5T_NATIVE_FLOAT, meta_memspace, meta_filespace, meta_plist_id, &grid->x0); + H5Dclose(meta_dset_id); + + meta_dset_id = H5Dcreate(meta_group_id, "y0", H5T_NATIVE_FLOAT, meta_filespace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + ierr = H5Dwrite(meta_dset_id, H5T_NATIVE_FLOAT, meta_memspace, meta_filespace, meta_plist_id, &grid->y0); + H5Dclose(meta_dset_id); + + meta_dset_id = H5Dcreate(meta_group_id, "z0", H5T_NATIVE_FLOAT, meta_filespace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + ierr = H5Dwrite(meta_dset_id, H5T_NATIVE_FLOAT, meta_memspace, meta_filespace, meta_plist_id, &grid->z0); + H5Dclose(meta_dset_id); + + meta_dset_id = H5Dcreate(meta_group_id, "dx", H5T_NATIVE_FLOAT, meta_filespace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + ierr = H5Dwrite(meta_dset_id, H5T_NATIVE_FLOAT, meta_memspace, meta_filespace, meta_plist_id, &grid->dx); + H5Dclose(meta_dset_id); + + meta_dset_id = H5Dcreate(meta_group_id, "dy", H5T_NATIVE_FLOAT, meta_filespace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + ierr = H5Dwrite(meta_dset_id, H5T_NATIVE_FLOAT, meta_memspace, meta_filespace, meta_plist_id, &grid->dy); + H5Dclose(meta_dset_id); + + meta_dset_id = H5Dcreate(meta_group_id, "dz", H5T_NATIVE_FLOAT, meta_filespace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + ierr = H5Dwrite(meta_dset_id, H5T_NATIVE_FLOAT, meta_memspace, meta_filespace, meta_plist_id, &grid->dz); + H5Dclose(meta_dset_id); + + meta_el2 = uptime() - meta_el2; + //sim_log("Metafile TimeHDF5Write: " << meta_el2 << " s"); + double meta_el3 = uptime(); + H5Sclose(meta_memspace); + H5Sclose(meta_filespace); + H5Pclose(meta_plist_id); + H5Gclose(meta_group_id); + H5Fclose(meta_file_id); + meta_el3 = uptime() - meta_el3; + //sim_log("Metafile TimeHDF5Close: " << meta_el3 << " s"); + + sp = sp->next; + } +} +#endif + void vpic_simulation::dump_particles( const char *sp_name, const char *fbase, @@ -699,9 +1871,6 @@ 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; diff --git a/src/vpic/hdf5_header_info.h b/src/vpic/hdf5_header_info.h new file mode 100644 index 00000000..baed8f7d --- /dev/null +++ b/src/vpic/hdf5_header_info.h @@ -0,0 +1,122 @@ +#ifndef VPIC_HDF5_HEAD_INFO +#define VPIC_HDF5_HEAD_INFO + +#define FIELD_ARRAY_NAME field_array + +// XML header stuff +const char *header = "\n\n\n\t\n"; +const char *header_topology = "\t\t\n"; +const char *header_geom = "\t\t\n"; +const char *header_origin = "\t\t\t \n\t\t\t%s\n"; +const char *header_dxdydz = "\t\t\t \n\t\t\t%s\n"; +const char *footer_geom = "\t\t\n"; +const char *grid_line = "\t\t \n \ +\t\t\t\n"; +const char *footer = "\t\t\n\t\n\n"; + +const char *main_body_head = "\t\t\t \n \ +\t\t\t\t \n \ +\t\t\t\t \n"; +const char *main_body_foot = "\t\t\t\n"; + +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 "; + +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 "; + +#define create_file_with_header(xml_file_name, dimensions, orignal, dxdydz, nframes, fields_interval) \ + { \ + FILE *fp; \ + fp = fopen(xml_file_name, "w"); \ + fputs(header, fp); \ + fprintf(fp, header_topology, dimensions); \ + fputs(header_geom, fp); \ + fprintf(fp, header_origin, orignal); \ + fprintf(fp, header_dxdydz, dxdydz); \ + fputs(footer_geom, fp); \ + fprintf(fp, grid_line, nframes); \ + int i; \ + for (i = 0; i < nframes; i++) \ + fprintf(fp, "%d ", i*fields_interval); \ + fputs(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, main_body_head, time_step); \ + if (field_dump_flag.enabledE()) \ + write_main_body_attribute(fp, main_body_attributeV, "E", dims_4d, dims_3d, speciesname_p, time_step, "ex", "ey", "ez"); \ + if (field_dump_flag.div_e_err) \ + fprintf(fp, 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, main_body_attributeV, "B", dims_4d, dims_3d, speciesname_p, time_step, "cbx", "cby", "cbz"); \ + if (field_dump_flag.div_b_err) \ + fprintf(fp, 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, main_body_attributeV, "TCA", dims_4d, dims_3d, speciesname_p, time_step, "tcax", "tcay", "tcaz"); \ + if (field_dump_flag.rhob) \ + fprintf(fp, 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, main_body_attributeV, "JF", dims_4d, dims_3d, speciesname_p, time_step, "jfx", "jfy", "jfz"); \ + if (field_dump_flag.rhof) \ + fprintf(fp, 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, main_body_attributeV, "EMAT", dims_4d, dims_3d, speciesname_p, time_step, "ematx", "ematy", "ematz"); \ + if (field_dump_flag.nmat) \ + fprintf(fp, 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, main_body_attributeV, "FMAT", dims_4d, dims_3d, speciesname_p, time_step, "fmatx", "fmaty", "fmatz"); \ + if (field_dump_flag.cmat) \ + fprintf(fp, main_body_attributeS, "cmat", dims_3d, time_step, speciesname_p, time_step, time_step, "cmat"); \ + fprintf(fp, "%s", main_body_foot); \ + if (add_footer_flag) \ + fputs(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, main_body_head, time_step); \ + if (hydro_dump_flag.enabledJ()) \ + write_main_body_attribute(fp, main_body_attributeV, "J", dims_4d, dims_3d, speciesname_p, time_step, "jx", "jy", "jz"); \ + if (hydro_dump_flag.rho) \ + fprintf(fp, 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, main_body_attributeV, "P", dims_4d, dims_3d, speciesname_p, time_step, "px", "py", "pz"); \ + if (hydro_dump_flag.ke) \ + fprintf(fp, 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, 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, main_body_attributeV, "TOD", dims_4d, dims_3d, speciesname_p, time_step, "tyz", "tzx", "txy"); \ + fprintf(fp, "%s", main_body_foot); \ + if (add_footer_flag) \ + fputs(footer, fp); \ + fclose(fp); \ + } + + +#endif // VPIC_HDF5_HEAD_INFO diff --git a/src/vpic/vpic.cc b/src/vpic/vpic.cc index 4f150afc..bfd18767 100644 --- a/src/vpic/vpic.cc +++ b/src/vpic/vpic.cc @@ -7,9 +7,9 @@ * March/April 2004 - Heavily revised and extended from earlier V4PIC versions * */ - + #include "vpic.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 @@ -99,7 +99,7 @@ vpic_simulation::vpic_simulation() { // if( n_rng