From 24e4a9a88368c7e9dc245cac8f5e5a2e5a89c174 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Tue, 19 Mar 2024 19:23:03 -0700 Subject: [PATCH] fix particle initialization so all particles are actually in the domain; also add the ability to specify where particles are created by defining a box in physical space (#1514) --- Exec/RegTests/ParticlesOverWoA/inputs | 6 +- Source/Particles/ERFPC.H | 7 +- Source/Particles/ERFPCEvolve.cpp | 2 +- Source/Particles/ERFPCInitializations.cpp | 213 ++++++++++++++++++---- Source/Particles/ParticleData.H | 8 +- 5 files changed, 195 insertions(+), 41 deletions(-) diff --git a/Exec/RegTests/ParticlesOverWoA/inputs b/Exec/RegTests/ParticlesOverWoA/inputs index 725789c3f..3633f1ba2 100644 --- a/Exec/RegTests/ParticlesOverWoA/inputs +++ b/Exec/RegTests/ParticlesOverWoA/inputs @@ -1,5 +1,5 @@ # ------------------ INPUTS TO MAIN PROGRAM ------------------- -max_step = 40 +max_step = 80 amrex.fpe_trap_invalid = 1 @@ -25,6 +25,10 @@ zhi.type = "SlipWall" # PARTICLES erf.use_tracer_particles = 1 +tracer_particles.initial_distribution_type = box +tracer_particles.particle_box_lo = 3.95 -1.0 -1.0 +tracer_particles.particle_box_hi = 4.00 2.0 1.0 +tracer_particles.place_randomly_in_cells = true # TIME STEP CONTROL erf.fixed_dt = 1E-3 diff --git a/Source/Particles/ERFPC.H b/Source/Particles/ERFPC.H index e6b854292..4574b8700 100644 --- a/Source/Particles/ERFPC.H +++ b/Source/Particles/ERFPC.H @@ -43,6 +43,7 @@ namespace ERFParticleInitializations { /* list of particle initializations */ const std::string init_default = "default"; + const std::string init_box = "box"; const std::string init_uniform = "uniform"; } @@ -215,12 +216,17 @@ class ERFPC : public amrex::ParticleContainer< ERFParticlesRealIdxAoS::ncomps, void initializeParticlesDefaultHydro (const std::unique_ptr& a_ptr = nullptr); /*! Default particle initialization */ void initializeParticlesUniformDistribution (const std::unique_ptr& a_ptr = nullptr); + void initializeParticlesUniformDistributionInBox (const std::unique_ptr& a_ptr, + const amrex::RealBox& particle_box); protected: bool m_advect_w_flow; /*!< advect with flow velocity */ bool m_advect_w_gravity; /*!< advect under gravitational force */ + amrex::RealBox particle_box; + bool place_randomly_in_cells; + std::string m_name; /*!< name of this particle species */ std::string m_initialization_type; /*!< initial particle distribution type */ @@ -233,7 +239,6 @@ class ERFPC : public amrex::ParticleContainer< ERFParticlesRealIdxAoS::ncomps, /*! Default particle initialization */ void initializeParticlesDefault (const std::unique_ptr& a_ptr = nullptr); - }; #endif diff --git a/Source/Particles/ERFPCEvolve.cpp b/Source/Particles/ERFPCEvolve.cpp index ba20d22bb..8d5fea7ee 100644 --- a/Source/Particles/ERFPCEvolve.cpp +++ b/Source/Particles/ERFPCEvolve.cpp @@ -64,7 +64,7 @@ void ERFPC::AdvectWithFlow ( MultiFab* a_umac, { for (int i = 0; i < AMREX_SPACEDIM; i++) { - int ng = a_umac[i].nGrow(); + IntVect ng = a_umac[i].nGrowVect(); raii_umac[i] = std::make_unique (convert(m_gdb->ParticleBoxArray(a_lev), IntVect::TheDimensionVector(i)), m_gdb->ParticleDistributionMap(a_lev), a_umac[i].nComp(), ng); diff --git a/Source/Particles/ERFPCInitializations.cpp b/Source/Particles/ERFPCInitializations.cpp index cc1dbb8e6..0ca781a8f 100644 --- a/Source/Particles/ERFPCInitializations.cpp +++ b/Source/Particles/ERFPCInitializations.cpp @@ -16,6 +16,30 @@ void ERFPC::readInputs () m_initialization_type = ERFParticleInitializations::init_default; pp.query("initial_distribution_type", m_initialization_type); + if (m_initialization_type == ERFParticleInitializations::init_box) + { + Vector particle_box_lo(AMREX_SPACEDIM); + Vector particle_box_hi(AMREX_SPACEDIM); + + // Defaults + for (int i = 0; i < AMREX_SPACEDIM; i++) { particle_box_lo[i] = Geom(0).ProbLo(i); } + for (int i = 0; i < AMREX_SPACEDIM; i++) { particle_box_hi[i] = Geom(0).ProbHi(i); } + + pp.queryAdd("particle_box_lo", particle_box_lo, AMREX_SPACEDIM); + AMREX_ASSERT(particle_box_lo.size() == AMREX_SPACEDIM); + + pp.queryAdd("particle_box_hi", particle_box_hi, AMREX_SPACEDIM); + AMREX_ASSERT(particle_box_hi.size() == AMREX_SPACEDIM); + + particle_box.setLo(particle_box_lo); + particle_box.setHi(particle_box_hi); + + // We default to placing the particles randomly within each cell, + // but can override this for regression testing + place_randomly_in_cells = true; + pp.query("place_randomly_in_cells", place_randomly_in_cells); + } + m_ppc_init = 1; pp.query("initial_particles_per_cell", m_ppc_init); @@ -35,6 +59,8 @@ void ERFPC::InitializeParticles (const std::unique_ptr& a_height_ptr) if (m_initialization_type == ERFParticleInitializations::init_default) { initializeParticlesDefault( a_height_ptr ); + } else if (m_initialization_type == ERFParticleInitializations::init_box) { + initializeParticlesUniformDistributionInBox( a_height_ptr , particle_box ); } else if (m_initialization_type == ERFParticleInitializations::init_uniform) { initializeParticlesUniformDistribution( a_height_ptr ); } else { @@ -76,7 +102,7 @@ void ERFPC::initializeParticlesDefaultTracersWoA (const std::unique_ptr& a_ Real x = plo[0] + (i + r[0])*dx[0]; Real y = plo[1] + (j + r[1])*dx[1]; - Real z = height_arr(i,j,k) - + r[2] * ( height_arr(i,j,k+1) - - height_arr(i,j,k) ); + Real zlo = 0.25 * (height_arr(i,j ,k ) + height_arr(i+1,j ,k ) + + height_arr(i,j+1,k ) + height_arr(i+1,j+1,k ) ); + Real zhi = 0.25 * (height_arr(i,j ,k+1) + height_arr(i+1,j ,k+1) + + height_arr(i,j+1,k+1) + height_arr(i+1,j+1,k ) ); + Real z = zlo + r[2] * (zhi - zlo); auto& p = aos[offset_arr(i,j,k)]; p.id() = pid + offset_arr(i,j,k); @@ -325,7 +355,15 @@ void ERFPC::initializeParticlesDefaultHydro (const std::unique_ptr& a_ * by "initial_particles_per_cell", and they are randomly distributed. */ void ERFPC::initializeParticlesUniformDistribution (const std::unique_ptr& a_height_ptr) { - BL_PROFILE("ERFPC::initializeParticlesUniformDistribution"); + initializeParticlesUniformDistributionInBox ( a_height_ptr, Geom(0).ProbDomain() ); +} + +/*! Uniform distribution: the number of particles per grid cell is specified + * by "initial_particles_per_cell", and they are randomly distributed. */ +void ERFPC::initializeParticlesUniformDistributionInBox (const std::unique_ptr& a_height_ptr, + const RealBox& particle_init_domain) +{ + BL_PROFILE("ERFPC::initializeParticlesUniformDistributionInBox"); const int lev = 0; const auto dx = Geom(lev).CellSizeArray(); @@ -340,10 +378,32 @@ void ERFPC::initializeParticlesUniformDistribution (const std::unique_ptr(pid + np) < LastParticleID, "Error: overflow on particle id numbers!" ); - if (a_height_ptr) { + if (a_height_ptr && place_randomly_in_cells) { + + const auto height_arr = (*a_height_ptr)[mfi].array(); - const auto height_arr = (*a_height_ptr)[mfi].array(); ParallelForRNG(tile_box, [=] AMREX_GPU_DEVICE (int i, int j, int k, const RandomEngine& rnd_engine) noexcept { int start = offset_arr(i,j,k); - for (int n = start; n < start+particles_per_cell; n++) { + for (int n = start; n < start+num_particles_arr(i,j,k); n++) { Real r[3] = {Random(rnd_engine), Random(rnd_engine), Random(rnd_engine)}; Real v[3] = {0.0, 0.0, 0.0}; Real x = plo[0] + (i + r[0])*dx[0]; Real y = plo[1] + (j + r[1])*dx[1]; - Real z = height_arr(i,j,k) - + r[2] * ( height_arr(i,j,k+1) - - height_arr(i,j,k) ); + + Real sx[] = { amrex::Real(1.) - r[0], r[0]}; + Real sy[] = { amrex::Real(1.) - r[1], r[1]}; + + Real height_at_pxy_lo = 0.; + for (int ii = 0; ii < 2; ++ii) { + for (int jj = 0; jj < 2; ++jj) { + height_at_pxy_lo += sx[ii] * sy[jj] * height_arr(i+ii,j+jj,k); + } + } + Real height_at_pxy_hi = 0.; + for (int ii = 0; ii < 2; ++ii) { + for (int jj = 0; jj < 2; ++jj) { + height_at_pxy_hi += sx[ii] * sy[jj] * height_arr(i+ii,j+jj,k+1); + } + } + + Real z = height_at_pxy_lo + r[2] * (height_at_pxy_hi - height_at_pxy_lo); auto& p = aos[n]; p.id() = pid + n; p.cpu() = my_proc; - p.pos(0) = x; - p.pos(1) = y; - p.pos(2) = z; + + p.pos(0) = x; p.pos(1) = y; p.pos(2) = z; p.idata(ERFParticlesIntIdxAoS::k) = k; - vx_ptr[n] = v[0]; - vy_ptr[n] = v[1]; - vz_ptr[n] = v[2]; + vx_ptr[n] = v[0]; vy_ptr[n] = v[1]; vz_ptr[n] = v[2]; mass_ptr[n] = 1.0e-6; } }); - } else { + } else if (a_height_ptr && !place_randomly_in_cells) { + + const auto height_arr = (*a_height_ptr)[mfi].array(); + + ParallelFor(tile_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + int start = offset_arr(i,j,k); + for (int n = start; n < start+num_particles_arr(i,j,k); n++) { + Real r[3] = {0.3, 0.7, 0.25}; + Real v[3] = {0.0, 0.0, 0.0}; + + Real x = plo[0] + (i + r[0])*dx[0]; + Real y = plo[1] + (j + r[1])*dx[1]; + + Real sx[] = { amrex::Real(1.) - r[0], r[0]}; + Real sy[] = { amrex::Real(1.) - r[1], r[1]}; + + Real height_at_pxy_lo = 0.; + for (int ii = 0; ii < 2; ++ii) { + for (int jj = 0; jj < 2; ++jj) { + height_at_pxy_lo += sx[ii] * sy[jj] * height_arr(i+ii,j+jj,k); + } + } + Real height_at_pxy_hi = 0.; + for (int ii = 0; ii < 2; ++ii) { + for (int jj = 0; jj < 2; ++jj) { + height_at_pxy_hi += sx[ii] * sy[jj] * height_arr(i+ii,j+jj,k+1); + } + } + + Real z = height_at_pxy_lo + r[2] * (height_at_pxy_hi - height_at_pxy_lo); + + auto& p = aos[n]; + p.id() = pid + n; + p.cpu() = my_proc; + + p.pos(0) = x; p.pos(1) = y; p.pos(2) = z; + + p.idata(ERFParticlesIntIdxAoS::k) = k; + + vx_ptr[n] = v[0]; vy_ptr[n] = v[1]; vz_ptr[n] = v[2]; + + mass_ptr[n] = 1.0e-6; + } + }); + + } else if (!a_height_ptr && place_randomly_in_cells) { ParallelForRNG(tile_box, [=] AMREX_GPU_DEVICE (int i, int j, int k, const RandomEngine& rnd_engine) noexcept { int start = offset_arr(i,j,k); - for (int n = start; n < start+particles_per_cell; n++) { + for (int n = start; n < start+num_particles_arr(i,j,k); n++) { Real r[3] = {Random(rnd_engine), Random(rnd_engine), Random(rnd_engine)}; Real v[3] = {0.0, 0.0, 0.0}; @@ -436,15 +557,39 @@ void ERFPC::initializeParticlesUniformDistribution (const std::unique_ptr