diff --git a/Exec/ParticleFilterTest/DarkMatterParticleContainer.H b/Exec/ParticleFilterTest/DarkMatterParticleContainer.H index afa40bbf..638cab59 100644 --- a/Exec/ParticleFilterTest/DarkMatterParticleContainer.H +++ b/Exec/ParticleFilterTest/DarkMatterParticleContainer.H @@ -64,7 +64,10 @@ AMREX_GPU_HOST_DEVICE AMREX_INLINE void store_dm_particle_single (amrex::Particl amrex::GpuArray const& phi, amrex::GpuArray const& dxi, const amrex::Real& dt, const amrex::Real& a_prev, - const amrex::Real& a_cur, const int& do_move, const amrex::Real& radius_inner, const amrex::Real& radius_outer, int index); + const amrex::Real& a_cur, const int& do_move, + const amrex::Real& radius_inner, + const amrex::Real& radius_outer, + int index, const bool is_file_write=false); #endif /* _DarkMatterParticleContainer_H_ */ diff --git a/Exec/ParticleFilterTest/DarkMatterParticleContainer.cpp b/Exec/ParticleFilterTest/DarkMatterParticleContainer.cpp index 139f8829..c3979c3e 100644 --- a/Exec/ParticleFilterTest/DarkMatterParticleContainer.cpp +++ b/Exec/ParticleFilterTest/DarkMatterParticleContainer.cpp @@ -7,7 +7,7 @@ using namespace amrex; #include - +#include /// These are helper functions used when initializing from a morton-ordered /// binary particle file. @@ -461,19 +461,19 @@ DarkMatterParticleContainer::moveKickDrift (amrex::MultiFab& acceleration, Array4 accel= accel_fab.array(); int nc=AMREX_SPACEDIM; - int num_output=0; + int num_output=0; for(int i=0;i::SuperParticleType p = pstruct2[i]; int mask=shell_filter_test(ptile.getConstParticleTileData(),i); if(mask>0) { for(int j=0;jWritePlotFile(dir, name, real_comp_names_shell); + if(write_hdf5!=1){ + //ShellPC->WritePlotFile(dir, name, real_comp_names_shell); + } #ifdef AMREX_USE_HDF5 - else + else{ ShellPC->WritePlotFileHDF5(dir, name, real_comp_names_shell, compression); + } #endif } Print()<<"After write\t"<TotalNumberOfParticles()<<"\t"<::Su amrex::GpuArray const& phi, amrex::GpuArray const& dxi, const amrex::Real& dt, const amrex::Real& a_prev, - const amrex::Real& a_cur, const int& do_move, const Real& radius_inner, const Real& radius_outer, int index) + const amrex::Real& a_cur, const int& do_move, const Real& radius_inner, + const Real& radius_outer, int index, const bool is_file_write) { amrex::Real half_dt = 0.5 * dt; amrex::Real a_cur_inv = 1.0 / a_cur; amrex::Real dt_a_cur_inv = dt * a_cur_inv; const GpuArray center({AMREX_D_DECL((phi[0]-plo[0])*0.5,(phi[1]-plo[1])*0.5,(phi[2]-plo[2])*0.5)}); - if (do_move == 1) - { - p2.rdata(0)=p.rdata(0); - bool result=false; - Real lenx = phi[0]-plo[0]; - Real leny = phi[1]-plo[1]; - Real lenz = phi[2]-plo[2]; - int maxind[3]; - maxind[0] = floor((radius_outer+lenx*0.5)/lenx); - maxind[1] = floor((radius_outer+leny*0.5)/leny); - maxind[2] = floor((radius_outer+lenz*0.5)/lenz); - - Real xlen, ylen, zlen; - //printf("Value is %d\n", maxind); - int local_index=-1; - for(int idir=-maxind[0];idir<=maxind[0];idir++) - for(int jdir=-maxind[1];jdir<=maxind[1];jdir++) - for(int kdir=-maxind[2];kdir<=maxind[2];kdir++) - { - xlen = p.pos(0)+(idir)*(phi[0]-plo[0]) - center[0]; - ylen = p.pos(1)+(jdir)*(phi[1]-plo[1]) - center[1]; - zlen = p.pos(2)+(kdir)*(phi[2]-plo[2]) - center[2]; - Real mag = sqrt(xlen*xlen+ylen*ylen+zlen*zlen); - result=result? true : (mag>radius_inner && magradius_inner && magradius_inner && magradius_inner && mag shell_particles; + +void writeBinaryVTK(const std::string& filename, + const std::vector& shell_particles); +void writeBinarySimple(const std::string& filename, + const std::vector& shell_particles); diff --git a/Exec/ParticleFilterTest/LightConeParticle.cpp b/Exec/ParticleFilterTest/LightConeParticle.cpp new file mode 100644 index 00000000..4ae779e5 --- /dev/null +++ b/Exec/ParticleFilterTest/LightConeParticle.cpp @@ -0,0 +1,156 @@ +#include +#include +#include +#include +#include + +#include + +void SwapEnd(float& val) { + // Swap endianess if necessary + char* bytes = reinterpret_cast(&val); + std::swap(bytes[0], bytes[3]); + std::swap(bytes[1], bytes[2]); +} + +void writeBinaryVTK(const std::string& filename, const std::vector& particles) { + int rank, size; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + MPI_Comm_size(MPI_COMM_WORLD, &size); + + size_t local_num_particles = particles.size(); + size_t total_num_particles = 0; + + // Get total particles across all ranks + MPI_Reduce(&local_num_particles, &total_num_particles, 1, MPI_UNSIGNED_LONG, MPI_SUM, 0, MPI_COMM_WORLD); + + // Compute offset for this rank's data + size_t offset = 0; + MPI_Exscan(&local_num_particles, &offset, 1, MPI_UNSIGNED_LONG, MPI_SUM, MPI_COMM_WORLD); + + // Header handling + size_t header_size = 0; + + if (rank == 0) { + std::ofstream file(filename, std::ios::binary); + if (!file) { + std::cerr << "Error: Could not open file " << filename << "\n"; + MPI_Abort(MPI_COMM_WORLD, 1); + } + + // Write the header + file << "# vtk DataFile Version 2.0\n"; + file << "Particle Cloud Data\n"; + file << "BINARY\n"; + file << "DATASET POLYDATA\n"; + file << "POINTS " << total_num_particles << " float\n"; + + // Determine header size + file.seekp(0, std::ios::end); + header_size = file.tellp(); + file.close(); + } + + // Broadcast the header size to all ranks + MPI_Bcast(&header_size, 1, MPI_UNSIGNED_LONG, 0, MPI_COMM_WORLD); + + // Use MPI collective I/O for binary data + MPI_File mpi_file; + MPI_File_open(MPI_COMM_WORLD, filename.c_str(), + MPI_MODE_WRONLY | MPI_MODE_APPEND, MPI_INFO_NULL, &mpi_file); + + // Compute byte offset for this rank + size_t byte_offset = header_size + sizeof(float) * 3 * offset; + + // Prepare local data + std::vector data(3 * local_num_particles); + for (size_t i = 0; i < local_num_particles; ++i) { + data[3 * i] = particles[i].x; + data[3 * i + 1] = particles[i].y; + data[3 * i + 2] = particles[i].z; + + // Convert to big-endian if needed + SwapEnd(data[3 * i]); + SwapEnd(data[3 * i + 1]); + SwapEnd(data[3 * i + 2]); + } + + // Write particle data collectively + MPI_File_write_at_all(mpi_file, byte_offset, data.data(), data.size(), MPI_FLOAT, MPI_STATUS_IGNORE); + + MPI_File_close(&mpi_file); + + if (rank == 0) { + std::cout << "Successfully wrote VTK file: " << filename << "\n"; + } +} + +void writeBinarySimple(const std::string& filename, const std::vector& particles) { + int rank, size; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + MPI_Comm_size(MPI_COMM_WORLD, &size); + + size_t local_num_particles = particles.size(); + size_t total_num_particles = 0; + + // Get total particles across all ranks + MPI_Reduce(&local_num_particles, &total_num_particles, 1, MPI_UNSIGNED_LONG, MPI_SUM, 0, MPI_COMM_WORLD); + + // Compute offset for this rank's data + size_t offset = 0; + MPI_Exscan(&local_num_particles, &offset, 1, MPI_UNSIGNED_LONG, MPI_SUM, MPI_COMM_WORLD); + + // Header handling + size_t header_size = 0; + + if (rank == 0) { + std::ofstream file(filename, std::ios::binary); + if (!file) { + std::cerr << "Error: Could not open file " << filename << "\n"; + MPI_Abort(MPI_COMM_WORLD, 1); + } + file.seekp(0, std::ios::end); + header_size = file.tellp(); + file.close(); + } + + // Broadcast the header size to all ranks + MPI_Bcast(&header_size, 1, MPI_UNSIGNED_LONG, 0, MPI_COMM_WORLD); + + // Use MPI collective I/O for binary data + MPI_File mpi_file; + MPI_File_open(MPI_COMM_WORLD, filename.c_str(), + MPI_MODE_WRONLY | MPI_MODE_APPEND, MPI_INFO_NULL, &mpi_file); + + // Compute byte offset for this rank + size_t byte_offset = header_size + sizeof(float) * 6 * offset; + + // Prepare local data + std::vector data(6 * local_num_particles); + for (size_t i = 0; i < local_num_particles; ++i) { + data[6 * i] = particles[i].x; + data[6 * i + 1] = particles[i].y; + data[6 * i + 2] = particles[i].z; + data[6 * i + 3] = particles[i].vx; + data[6 * i + 4] = particles[i].vy; + data[6 * i + 5] = particles[i].vz; + + // Convert to big-endian if needed + SwapEnd(data[6 * i]); + SwapEnd(data[6 * i + 1]); + SwapEnd(data[6 * i + 2]); + SwapEnd(data[6 * i + 3]); + SwapEnd(data[6 * i + 4]); + SwapEnd(data[6 * i + 5]); + } + + // Write particle data collectively + MPI_File_write_at_all(mpi_file, byte_offset, data.data(), data.size(), MPI_FLOAT, MPI_STATUS_IGNORE); + + MPI_File_close(&mpi_file); + + if (rank == 0) { + std::cout << "Successfully wrote VTK file: " << filename << "\n"; + } +} + diff --git a/Exec/ParticleFilterTest/Make.package b/Exec/ParticleFilterTest/Make.package index 86de0a5f..7ac175cd 100644 --- a/Exec/ParticleFilterTest/Make.package +++ b/Exec/ParticleFilterTest/Make.package @@ -1,2 +1,4 @@ CEXE_sources += Prob.cpp CEXE_sources += Prob_error.cpp +CEXE_sources += LightConeParticle.cpp +CEXE_headers += LightConeParticle.H diff --git a/Exec/ParticleFilterTest/Nyx_advance.cpp b/Exec/ParticleFilterTest/Nyx_advance.cpp index e68203c1..9c820350 100644 --- a/Exec/ParticleFilterTest/Nyx_advance.cpp +++ b/Exec/ParticleFilterTest/Nyx_advance.cpp @@ -3,11 +3,14 @@ #include #include #include +#include using namespace amrex; using std::string; +std::vector shell_particles; + Real Nyx::advance (Real time, Real dt, @@ -240,6 +243,13 @@ Nyx::advance_hydro_plus_particles (Real time, MultiFab::RegionTag amrMoveKickDrift_tag("MoveKickDrift_" + std::to_string(level)); + shell_particles.clear(); + shell_particles.shrink_to_fit(); + + std::string filename=Concatenate("lightcone_", int(100*(1/a_old-1)), 7); + std::string filename_vtk=filename+".vtk"; + std::string filename_bin=filename+".bin"; + for (int lev = level; lev <= finest_level_to_advance; lev++) { // We need grav_n_grow grow cells to track boundary particles @@ -260,6 +270,9 @@ Nyx::advance_hydro_plus_particles (Real time, } Nyx::theActiveParticles()[i]->moveKickDrift(grav_vec_old, lev, time, dt, a_old, a_half, where_width, radius_inner, radius_outer); } + writeBinaryVTK(filename_vtk, shell_particles); + writeBinarySimple(filename_bin, shell_particles); + // Only need the coarsest virtual particles here. if (lev == level && level < finest_level) diff --git a/Exec/ParticleFilterTest/Read_BinarySimple.py b/Exec/ParticleFilterTest/Read_BinarySimple.py new file mode 100644 index 00000000..2c42a8cc --- /dev/null +++ b/Exec/ParticleFilterTest/Read_BinarySimple.py @@ -0,0 +1,52 @@ +import numpy as np + +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D + +def read_binary_points(filename): + with open(filename, 'rb') as file: + # Read the binary data directly, assuming big-endian float32 + data = np.fromfile(file, dtype='>f4') # '>f4' indicates big-endian float32 + + # Reshape the data to Nx3 (x, y, z for each point) + if len(data) % 6 != 0: + raise ValueError("Data size is not a multiple of 6. The file might be corrupted or improperly formatted.") + + points = data.reshape((-1, 6)) + + num_points = len(data) // 6 + print(f"Total number of points: {num_points}") + + return points, num_points + +def plot_points(x, y, z): + fig = plt.figure(figsize=(10, 7)) + ax = fig.add_subplot(111, projection='3d') + + # Scatter plot for points + ax.scatter(x, y, z, c='blue', marker='o', s=1, alpha=0.8) + + # Set labels for the axes + ax.set_xlabel('X Coordinate') + ax.set_ylabel('Y Coordinate') + ax.set_zlabel('Z Coordinate') + + ax.set_title('Lightcone shell') + plt.show() + + +# Example usage +filename = "lightcone_0000010.bin" +points, num_points = read_binary_points(filename) + +# Display the first few points +x = np.zeros(num_points); +y = np.zeros(num_points); +z = np.zeros(num_points); +for i, point in enumerate(points[:num_points]): # Adjust the slice as needed + x[i], y[i], z[i], vx, vy, vz = point + print(f"{x[i]:.15g}, {y[i]:.15g}, {z[i]:.15g}, {vx:.15g}, {vy:.15g}, {vz:.15g}") + +plot_points(x, y, z) + +