From 9d432d687f30ec65cdff0db2b68fee43cac69e15 Mon Sep 17 00:00:00 2001 From: Debojyoti Ghosh Date: Wed, 13 Dec 2023 18:14:34 -0800 Subject: [PATCH 1/2] initial implementation of a generalized ERF particle interface; retained tracer and hydro particles as specific types; added uniform distribution of particles as an initialization option --- Source/ERF.H | 30 +- Source/ERF.cpp | 14 +- Source/IO/Plotfile.cpp | 39 +-- Source/Particles/ERFPC.H | 173 ++++++++++ Source/Particles/ERFPCEvolve.cpp | 258 ++++++++++++++ Source/Particles/ERFPCInitializations.cpp | 391 ++++++++++++++++++++++ Source/Particles/ERFTracers.cpp | 86 +++++ Source/Particles/HydroPC.H | 68 ---- Source/Particles/HydroPC.cpp | 221 ------------ Source/Particles/Make.package | 8 +- Source/Particles/ParticleData.H | 251 +++++++++----- Source/Particles/TracerPC.H | 68 ---- Source/Particles/TracerPC.cpp | 261 --------------- Source/TimeIntegration/ERF_Advance.cpp | 2 +- 14 files changed, 1117 insertions(+), 753 deletions(-) create mode 100644 Source/Particles/ERFPC.H create mode 100644 Source/Particles/ERFPCEvolve.cpp create mode 100644 Source/Particles/ERFPCInitializations.cpp create mode 100644 Source/Particles/ERFTracers.cpp delete mode 100644 Source/Particles/HydroPC.H delete mode 100644 Source/Particles/HydroPC.cpp delete mode 100644 Source/Particles/TracerPC.H delete mode 100644 Source/Particles/TracerPC.cpp diff --git a/Source/ERF.H b/Source/ERF.H index 6b96dc255..8c8eaef19 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -26,10 +26,6 @@ #include #endif -#ifdef ERF_USE_PARTICLES -#include "ParticleData.H" -#endif - #include "prob_common.H" #include @@ -44,6 +40,10 @@ #include #include +#ifdef ERF_USE_PARTICLES +#include "ParticleData.H" +#endif + #include "Microphysics.H" #ifdef ERF_USE_RRTMGP @@ -660,9 +660,6 @@ private: ,"qt", "qp", "qv", "qc", "qi", "qrain", "qsnow", "qgraup" #ifdef ERF_COMPUTE_ERROR ,"xvel_err", "yvel_err", "zvel_err", "pp_err" -#endif -#ifdef ERF_USE_PARTICLES - ,"tracer_particle_count", "hydro_particle_count" #endif }; @@ -670,8 +667,25 @@ private: static SolverChoice solverChoice; #ifdef ERF_USE_PARTICLES - // Particle info + // Particle container with all particle species static ParticleData particleData; + + // variables and functions for tracers particles + bool m_use_tracer_particles; /*!< tracer particles that advect with flow */ + bool m_use_hydro_particles; /*!< tracer particles that fall with gravity */ + + /*! Read tracer and hydro particles parameters */ + void readTracersParams(); + + /*! Initialize tracer and hydro particles */ + void initializeTracers ( amrex::ParGDBBase*, + const amrex::Vector>& ); + /*! Evolve tracers and hydro particles */ + void evolveTracers( int, + amrex::Real, + amrex::Vector>&, + const amrex::Vector>& ); + #endif static int verbose; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index e368e0b4b..b0a8c4ddf 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -537,7 +537,7 @@ ERF::InitData () } #ifdef ERF_USE_PARTICLES - particleData.init_particles((amrex::ParGDBBase*)GetParGDB(),z_phys_nd); + initializeTracers((amrex::ParGDBBase*)GetParGDB(),z_phys_nd); #endif } else { // Restart from a checkpoint @@ -945,10 +945,6 @@ ERF::ReadParameters () pp.query("fixed_fast_dt", fixed_fast_dt); pp.query("fixed_mri_dt_ratio", fixed_mri_dt_ratio); -#ifdef ERF_USE_PARTICLES - particleData.init_particle_params(); -#endif - // If this is set, it must be even if (fixed_mri_dt_ratio > 0 && (fixed_mri_dt_ratio%2 != 0) ) { @@ -1094,12 +1090,12 @@ ERF::ReadParameters () if (!iterate) SetIterateToFalse(); } -#ifdef ERF_USE_MULTIBLOCK - solverChoice.pp_prefix = pp_prefix; +#ifdef ERF_USE_PARTICLES + readTracersParams(); #endif -#ifdef ERF_USE_PARTICLES - particleData.init_particle_params(); +#ifdef ERF_USE_MULTIBLOCK + solverChoice.pp_prefix = pp_prefix; #endif solverChoice.init_params(max_level); diff --git a/Source/IO/Plotfile.cpp b/Source/IO/Plotfile.cpp index b8327e13d..0b06da74c 100644 --- a/Source/IO/Plotfile.cpp +++ b/Source/IO/Plotfile.cpp @@ -63,17 +63,20 @@ ERF::setPlotVariables (const std::string& pp_plot_var_names, Vector } for (int i = 0; i < derived_names.size(); ++i) { if ( containerHasElement(plot_var_names, derived_names[i]) ) { -#ifdef ERF_USE_PARTICLES - if ( (solverChoice.use_terrain || (derived_names[i] != "z_phys" && derived_names[i] != "detJ") ) && - (particleData.use_tracer_particles || (derived_names[i] != "tracer_particle_count") ) && - (particleData.use_hydro_particles || (derived_names[i] != "hydro_particle_count") ) ) { -#else if (solverChoice.use_terrain || (derived_names[i] != "z_phys" && derived_names[i] != "detJ") ) { -#endif tmp_plot_names.push_back(derived_names[i]); } } } +#ifdef ERF_USE_PARTICLES + const auto& particles_namelist( particleData.getNamesUnalloc() ); + for (auto it = particles_namelist.cbegin(); it != particles_namelist.cend(); ++it) { + std::string tmp( (*it)+"_count" ); + if (containerHasElement(plot_var_names, tmp) ) { + tmp_plot_names.push_back(tmp); + } + } +#endif // Check to see if we found all the requested variables for (const auto& plot_name : plot_var_names) { @@ -667,21 +670,15 @@ ERF::WritePlotFile (int which, Vector plot_var_names) } #ifdef ERF_USE_PARTICLES - if (containerHasElement(plot_var_names, "tracer_particle_count")) - { - MultiFab temp_dat(mf[lev].boxArray(), mf[lev].DistributionMap(), 1, 0); - temp_dat.setVal(0); - particleData.tracer_particles->Increment(temp_dat, lev); - MultiFab::Copy(mf[lev], temp_dat, 0, mf_comp, 1, 0); - mf_comp += 1; - } - if (containerHasElement(plot_var_names, "hydro_particle_count")) - { - MultiFab temp_dat(mf[lev].boxArray(), mf[lev].DistributionMap(), 1, 0); - temp_dat.setVal(0); - particleData.hydro_particles->Increment(temp_dat, lev); - MultiFab::Copy(mf[lev], temp_dat, 0, mf_comp, 1, 0); - mf_comp += 1; + const auto& particles_namelist( particleData.getNames() ); + for (ParticlesNamesVector::size_type i = 0; i < particles_namelist.size(); i++) { + if (containerHasElement(plot_var_names, std::string(particles_namelist[i]+"_count"))) { + MultiFab temp_dat(mf[lev].boxArray(), mf[lev].DistributionMap(), 1, 0); + temp_dat.setVal(0); + particleData[particles_namelist[i]]->Increment(temp_dat, lev); + MultiFab::Copy(mf[lev], temp_dat, 0, mf_comp, 1, 0); + mf_comp += 1; + } } #endif diff --git a/Source/Particles/ERFPC.H b/Source/Particles/ERFPC.H new file mode 100644 index 000000000..6a8be3542 --- /dev/null +++ b/Source/Particles/ERFPC.H @@ -0,0 +1,173 @@ +#ifndef ERF_PC_H_ +#define ERF_PC_H_ + +#ifdef ERF_USE_PARTICLES + +#include +#include + +struct ERFParticlesRealIdx +{ + enum { + vx = 0, + vy, + vz, + mass, + ncomps + }; +}; + +struct ERFParticlesIntIdx +{ + enum { + i = 0, + j, + k, + ncomps + }; +}; + +namespace ERFParticleInitializations +{ + /* list of particle initializations */ + const std::string init_default = "default"; + const std::string init_uniform = "uniform"; + +} + +namespace ERFParticleNames +{ + const std::string tracers = "tracer_particles"; + const std::string hydro = "hydro_particles"; +} + +struct ERFParticlesAssignor +{ + template + AMREX_GPU_HOST_DEVICE + amrex::IntVect operator() ( P const& p, + amrex::GpuArray const& plo, + amrex::GpuArray const& dxi, + const amrex::Box& domain ) const noexcept + { + /* TODO: What is this about? */ + amrex::IntVect iv( + AMREX_D_DECL( int(amrex::Math::floor((p.pos(0)-plo[0])*dxi[0])), + int(amrex::Math::floor((p.pos(1)-plo[1])*dxi[1])), + p.idata(ERFParticlesIntIdx::k) ) ); + iv[0] += domain.smallEnd()[0]; + iv[1] += domain.smallEnd()[1]; + return iv; + } +}; + +class ERFPC : public amrex::ParticleContainer< ERFParticlesRealIdx::ncomps, + ERFParticlesIntIdx::ncomps, + 0, + 0, + amrex::DefaultAllocator, + ERFParticlesAssignor > +{ + public: + + /*! Constructor */ + ERFPC ( amrex::ParGDBBase* a_gdb, + const std::string& a_name = "particles" ) + : amrex::ParticleContainer< ERFParticlesRealIdx::ncomps, + ERFParticlesIntIdx::ncomps, + 0, + 0, + amrex::DefaultAllocator, + ERFParticlesAssignor> (a_gdb) + { + BL_PROFILE("ERFPCPC::ERFPC()"); + m_name = a_name; + readInputs(); + } + + /*! Constructor */ + ERFPC ( const amrex::Geometry& a_geom, + const amrex::DistributionMapping& a_dmap, + const amrex::BoxArray& a_ba, + const std::string& a_name = "particles" ) + : amrex::ParticleContainer< ERFParticlesRealIdx::ncomps, + ERFParticlesIntIdx::ncomps, + 0, + 0, + amrex::DefaultAllocator, + ERFParticlesAssignor> ( a_geom, a_dmap, a_ba ) + { + BL_PROFILE("ERFPCPC::ERFPC()"); + m_name = a_name; + readInputs(); + } + + /*! Initialize particles in domain */ + virtual void InitializeParticles(const std::unique_ptr& a_ptr = nullptr); + + /*! Evolve particles for one time step */ + virtual void EvolveParticles( int, + amrex::Real, + amrex::Vector>&, + const amrex::Vector>& ); + + /*! Get real-type particle attribute names */ + virtual amrex::Vector varNames() const + { + BL_PROFILE("ERFPCPC::varNames()"); + return {AMREX_D_DECL("xvel","yvel","zvel"),"mass"}; + } + + /*! Specify if particles should advect with flow */ + inline void setAdvectWithFlow (bool a_flag) + { + BL_PROFILE("ERFPCPC::setAdvectWithFlow()"); + m_advect_w_flow = a_flag; + } + /*! Specify if particles fall under gravity */ + inline void setAdvectWithGravity (bool a_flag) + { + BL_PROFILE("ERFPCPC::setAdvectWithGravity()"); + m_advect_w_gravity = a_flag; + } + + protected: + + bool m_advect_w_flow; /*!< advect with flow velocity */ + bool m_advect_w_gravity; /*!< advect under gravitational force */ + + std::string m_name; /*!< name of this particle species */ + + std::string m_initialization_type; /*!< initial particle distribution type */ + int m_ppc_init; /*!< initial number of particles per cell */ + + /*! read inputs from file */ + virtual void readInputs(); + + /*! Default particle initialization */ + void initializeParticlesUniformDistribution(const std::unique_ptr& a_ptr = nullptr); + + /*! Uses midpoint method to advance particles using flow velocity. */ + void AdvectWithFlow( amrex::MultiFab*, + int, + amrex::Real, + const std::unique_ptr& ); + + /*! Uses midpoint method to advance particles falling under gravity. */ + void AdvectWithGravity( int, + amrex::Real, + const std::unique_ptr& ); + + private: + + /*! Default particle initialization */ + void initializeParticlesDefault(const std::unique_ptr& a_ptr = nullptr); + /*! Default initialization for tracer particles for WoA case (ref: AA) */ + void initializeParticlesDefaultTracersWoA(const std::unique_ptr& a_ptr=nullptr); + /*! Default initialization for hydro particles (ref: AA) */ + void initializeParticlesDefaultHydro(const std::unique_ptr& a_ptr = nullptr); + +}; + +#endif +#endif diff --git a/Source/Particles/ERFPCEvolve.cpp b/Source/Particles/ERFPCEvolve.cpp new file mode 100644 index 000000000..c46a759d5 --- /dev/null +++ b/Source/Particles/ERFPCEvolve.cpp @@ -0,0 +1,258 @@ +#include + +#ifdef ERF_USE_PARTICLES + +#include +#include +#include + +using namespace amrex; + +/*! Evolve particles for one time step */ +void ERFPC::EvolveParticles( int a_lev, + Real a_dt_lev, + Vector>& a_flow_vars, + const Vector>& a_z_phys_nd ) +{ + BL_PROFILE("ERFPCPC::EvolveParticles()"); + + if (m_advect_w_flow) { + MultiFab* flow_vel( &a_flow_vars[a_lev][Vars::xvel] ); + AdvectWithFlow( flow_vel, a_lev, a_dt_lev, a_z_phys_nd[a_lev] ); + } + + if (m_advect_w_gravity) { + AdvectWithGravity( a_lev, a_dt_lev, a_z_phys_nd[a_lev] ); + } + return; +} + +/*! Uses midpoint method to advance particles using flow velocity. */ +void ERFPC::AdvectWithFlow( MultiFab* a_umac, + int a_lev, + Real a_dt, + const std::unique_ptr& a_z_height ) +{ + BL_PROFILE("ERFPCPC::AdvectWithUmac()"); + AMREX_ASSERT(OK(a_lev, a_lev, a_umac[0].nGrow()-1)); + AMREX_ASSERT(a_lev >= 0 && a_lev < GetParticles().size()); + + AMREX_D_TERM(AMREX_ASSERT(a_umac[0].nGrow() >= 1);, + AMREX_ASSERT(a_umac[1].nGrow() >= 1);, + AMREX_ASSERT(a_umac[2].nGrow() >= 1);); + + AMREX_D_TERM(AMREX_ASSERT(!a_umac[0].contains_nan());, + AMREX_ASSERT(!a_umac[1].contains_nan());, + AMREX_ASSERT(!a_umac[2].contains_nan());); + + const auto strttime = amrex::second(); + const Geometry& geom = m_gdb->Geom(a_lev); + const Box& domain = geom.Domain(); + const auto plo = geom.ProbLoArray(); + const auto dxi = geom.InvCellSizeArray(); + const auto dx = geom.CellSizeArray(); + + Vector > raii_umac(AMREX_SPACEDIM); + Vector umac_pointer(AMREX_SPACEDIM); + if (OnSameGrids(a_lev, a_umac[0])) + { + for (int i = 0; i < AMREX_SPACEDIM; i++) { + umac_pointer[i] = &a_umac[i]; + } + } + else + { + for (int i = 0; i < AMREX_SPACEDIM; i++) + { + int ng = a_umac[i].nGrow(); + raii_umac[i] = std::make_unique + (amrex::convert(m_gdb->ParticleBoxArray(a_lev), IntVect::TheDimensionVector(i)), + m_gdb->ParticleDistributionMap(a_lev), a_umac[i].nComp(), ng); + umac_pointer[i] = raii_umac[i].get(); + umac_pointer[i]->ParallelCopy(a_umac[i],0,0,a_umac[i].nComp(),ng,ng); + } + } + + for (int ipass = 0; ipass < 2; ipass++) + { +#ifdef AMREX_USE_OMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (ParIterType pti(*this, a_lev); pti.isValid(); ++pti) + { + int grid = pti.index(); + auto& ptile = ParticlesAt(a_lev, pti); + auto& aos = ptile.GetArrayOfStructs(); + const int n = aos.numParticles(); + auto *p_pbox = aos().data(); + const FArrayBox* fab[AMREX_SPACEDIM] = { AMREX_D_DECL(&((*umac_pointer[0])[grid]), + &((*umac_pointer[1])[grid]), + &((*umac_pointer[2])[grid])) }; + + //array of these pointers to pass to the GPU + amrex::GpuArray, AMREX_SPACEDIM> + const umacarr {{AMREX_D_DECL((*fab[0]).array(), + (*fab[1]).array(), + (*fab[2]).array() )}}; + + bool use_terrain = (a_z_height != nullptr); + auto zheight = use_terrain ? (*a_z_height)[grid].array() : Array4{}; + + amrex::ParallelFor(n, + [=] AMREX_GPU_DEVICE (int i) + { + ParticleType& p = p_pbox[i]; + if (p.id() <= 0) { return; } + ParticleReal v[AMREX_SPACEDIM]; + mac_interpolate(p, plo, dxi, umacarr, v); + if (ipass == 0) + { + for (int dim=0; dim < AMREX_SPACEDIM; dim++) + { + p.rdata(dim) = p.pos(dim); + p.pos(dim) += static_cast(ParticleReal(0.5)*a_dt*v[dim]); + } + } + else + { + for (int dim=0; dim < AMREX_SPACEDIM; dim++) + { + p.pos(dim) = p.rdata(dim) + static_cast(a_dt*v[dim]); + p.rdata(dim) = v[dim]; + } + + // also update z-coordinate here + IntVect iv( + AMREX_D_DECL(int(amrex::Math::floor((p.pos(0)-plo[0])*dxi[0])), + int(amrex::Math::floor((p.pos(1)-plo[1])*dxi[1])), + p.idata(0))); + iv[0] += domain.smallEnd()[0]; + iv[1] += domain.smallEnd()[1]; + ParticleReal zlo, zhi; + if (use_terrain) { + zlo = zheight(iv[0], iv[1], iv[2]); + zhi = zheight(iv[0], iv[1], iv[2]+1); + } else { + zlo = iv[2] * dx[2]; + zhi = (iv[2]+1) * dx[2]; + } + if (p.pos(2) > zhi) { // need to be careful here + p.idata(0) += 1; + } else if (p.pos(2) <= zlo) { + p.idata(0) -= 1; + } + } + }); + } + } + + if (m_verbose > 1) + { + auto stoptime = amrex::second() - strttime; + +#ifdef AMREX_LAZY + Lazy::QueueReduction( [=] () mutable { +#endif + ParallelReduce::Max(stoptime, ParallelContext::IOProcessorNumberSub(), + ParallelContext::CommunicatorSub()); + + amrex::Print() << "ERFPC::AdvectWithFlow() time: " << stoptime << '\n'; +#ifdef AMREX_LAZY + }); +#endif + } +} + +void ERFPC::AdvectWithGravity( int a_lev, + Real a_dt, + const std::unique_ptr& a_z_height ) +{ + BL_PROFILE("ERFPC::AdvectWithGravity()"); + AMREX_ASSERT(a_lev >= 0 && a_lev < GetParticles().size()); + + const auto strttime = amrex::second(); + const Geometry& geom = m_gdb->Geom(a_lev); + const Box& domain = geom.Domain(); + const auto plo = geom.ProbLoArray(); + const auto dx = geom.CellSizeArray(); + const auto dxi = geom.InvCellSizeArray(); + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (ParIterType pti(*this, a_lev); pti.isValid(); ++pti) + { + int grid = pti.index(); + auto& ptile = ParticlesAt(a_lev, pti); + auto& aos = ptile.GetArrayOfStructs(); + const int n = aos.numParticles(); + auto *p_pbox = aos().data(); + + bool use_terrain = (a_z_height != nullptr); + auto zheight = use_terrain ? (*a_z_height)[grid].array() : Array4{}; + + amrex::ParallelFor(n, [=] AMREX_GPU_DEVICE (int i) + { + ParticleType& p = p_pbox[i]; + if (p.id() <= 0) { return; } + + ParticleReal v = p.rdata(ERFParticlesRealIdx::vz); + + // Define acceleration to be (gravity minus drag) where drag is defined + // such the particles will reach a terminal velocity of 5.0 (totally arbitrary) + ParticleReal terminal_vel = 5.0; + ParticleReal grav = CONST_GRAV; + ParticleReal drag = CONST_GRAV * (v * v) / (terminal_vel*terminal_vel); + + ParticleReal half_dt = 0.5 * a_dt; + + // Update the particle velocity over first half of step (a_dt/2) + p.rdata(ERFParticlesRealIdx::vz) -= (grav - drag) * half_dt; + + // Update the particle position over (a_dt) + p.pos(2) += static_cast(ParticleReal(0.5)*a_dt*p.rdata(ERFParticlesRealIdx::vz)); + + // Update the particle velocity over second half of step (a_dt/2) + p.rdata(ERFParticlesRealIdx::vz) -= (grav - drag) * half_dt; + + // also update z-coordinate here + IntVect iv( + AMREX_D_DECL(int(amrex::Math::floor((p.pos(0)-plo[0])*dxi[0])), + int(amrex::Math::floor((p.pos(1)-plo[1])*dxi[1])), + p.idata(0))); + iv[0] += domain.smallEnd()[0]; + iv[1] += domain.smallEnd()[1]; + ParticleReal zlo, zhi; + if (use_terrain) { + zlo = zheight(iv[0], iv[1], iv[2]); + zhi = zheight(iv[0], iv[1], iv[2]+1); + } else { + zlo = iv[2] * dx[2]; + zhi = (iv[2]+1) * dx[2]; + } + if (p.pos(2) > zhi) { // need to be careful here + p.idata(0) += 1; + } else if (p.pos(2) <= zlo) { + p.idata(0) -= 1; + } + }); + } + + if (m_verbose > 1) + { + auto stoptime = amrex::second() - strttime; + +#ifdef AMREX_LAZY + Lazy::QueueReduction( [=] () mutable { +#endif + ParallelReduce::Max(stoptime, ParallelContext::IOProcessorNumberSub(), + ParallelContext::CommunicatorSub()); + + amrex::Print() << "ERFPC::AdvectWithGravity() time: " << stoptime << '\n'; +#ifdef AMREX_LAZY + }); +#endif + } +} + +#endif diff --git a/Source/Particles/ERFPCInitializations.cpp b/Source/Particles/ERFPCInitializations.cpp new file mode 100644 index 000000000..3efc4a426 --- /dev/null +++ b/Source/Particles/ERFPCInitializations.cpp @@ -0,0 +1,391 @@ +#include +#include + +#ifdef ERF_USE_PARTICLES + +#include + +using namespace amrex; + +/*! Read inputs from file */ +void ERFPC::readInputs() +{ + BL_PROFILE("ERFPC::readInputs"); + + ParmParse pp(m_name); + + m_initialization_type = ERFParticleInitializations::init_default; + pp.query("initial_distribution_type", m_initialization_type); + + m_ppc_init = 1; + pp.query("inital_particles_per_cell", m_ppc_init); + + m_advect_w_flow = (m_name == ERFParticleNames::tracers ? true : false); + pp.query("advect_with_flow", m_advect_w_flow); + + m_advect_w_gravity = (m_name == ERFParticleNames::hydro ? true : false); + pp.query("advect_with_gravity", m_advect_w_gravity); + + return; +} + +/*! Initialize particles in domain */ +void ERFPC::InitializeParticles(const std::unique_ptr& a_height_ptr) +{ + BL_PROFILE("ERFPC::initializeParticles"); + + if (m_initialization_type == ERFParticleInitializations::init_default) { + initializeParticlesDefault( a_height_ptr ); + } else if (m_initialization_type == ERFParticleInitializations::init_uniform) { + initializeParticlesUniformDistribution( a_height_ptr ); + } else { + amrex::Print() << "Error: " << m_initialization_type + << " is not a valid initialization for " + << m_name << " particle species.\n"; + amrex::Error("See error message!"); + } + return; +} + +/*! Default particle initialization */ +void ERFPC::initializeParticlesDefault(const std::unique_ptr& a_height_ptr) +{ + BL_PROFILE(" ERFPC::initializeParticlesDefault"); + + if (m_name == ERFParticleNames::tracers) { + initializeParticlesDefaultTracersWoA( a_height_ptr ); + } else if (m_name == ERFParticleNames::hydro) { + initializeParticlesDefaultHydro( a_height_ptr ); + } +} + +/*! Default initialization for tracer particles for WoA case (ref: AA) */ +void ERFPC::initializeParticlesDefaultTracersWoA(const std::unique_ptr& a_height_ptr) +{ + BL_PROFILE("ERFPC::initializeParticlesDefaultTracersWoA"); + + const int lev = 0; + const Real* dx = Geom(lev).CellSize(); + const Real* plo = Geom(lev).ProbLo(); + + for(MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) + { + const Box& tile_box = mfi.tilebox(); + Gpu::HostVector host_particles; + + if (a_height_ptr) { + + const auto& height = (*a_height_ptr)[mfi]; + const FArrayBox* height_ptr = nullptr; +#ifdef AMREX_USE_GPU + std::unique_ptr hostfab; + if (height.arena()->isManaged() || height.arena()->isDevice()) { + hostfab = std::make_unique(height.box(), height.nComp(), + The_Pinned_Arena()); + Gpu::dtoh_memcpy_async(hostfab->dataPtr(), height.dataPtr(), + height.size()*sizeof(Real)); + Gpu::streamSynchronize(); + height_ptr = hostfab.get(); + } +#else + height_ptr = &height; +#endif + for (IntVect iv = tile_box.smallEnd(); iv <= tile_box.bigEnd(); tile_box.next(iv)) { + if (iv[0] == 3) { + Real r[3] = {0.5, 0.5, 0.5}; // this means place at cell center + Real v[3] = {0.0, 0.0, 0.0}; // with 0 initial velocity + + Real x = plo[0] + (iv[0] + r[0])*dx[0]; + Real y = plo[1] + (iv[1] + r[1])*dx[1]; + Real z = (*height_ptr)(iv) + + r[2] * ( (*height_ptr)(iv + IntVect(AMREX_D_DECL(0, 0, 1))) + - (*height_ptr)(iv) ); + + ParticleType p; + p.id() = ParticleType::NextID(); + p.cpu() = ParallelDescriptor::MyProc(); + p.pos(0) = x; + p.pos(1) = y; + p.pos(2) = z; + + p.rdata(ERFParticlesRealIdx::vx) = v[0]; + p.rdata(ERFParticlesRealIdx::vy) = v[1]; + p.rdata(ERFParticlesRealIdx::vz) = v[2]; + + p.idata(ERFParticlesIntIdx::i) = iv[0]; + p.idata(ERFParticlesIntIdx::j) = iv[1]; + p.idata(ERFParticlesIntIdx::k) = iv[2]; + + host_particles.push_back(p); + } + } + + } else { + + for (IntVect iv = tile_box.smallEnd(); iv <= tile_box.bigEnd(); tile_box.next(iv)) { + if (iv[0] == 3) { + Real r[3] = {0.5, 0.5, 0.5}; // this means place at cell center + Real v[3] = {0.0, 0.0, 0.0}; // with 0 initial velocity + + Real x = plo[0] + (iv[0] + r[0])*dx[0]; + Real y = plo[1] + (iv[1] + r[1])*dx[1]; + Real z = plo[2] + (iv[2] + r[2])*dx[2]; + + ParticleType p; + p.id() = ParticleType::NextID(); + p.cpu() = ParallelDescriptor::MyProc(); + p.pos(0) = x; + p.pos(1) = y; + p.pos(2) = z; + + p.rdata(ERFParticlesRealIdx::vx) = v[0]; + p.rdata(ERFParticlesRealIdx::vy) = v[1]; + p.rdata(ERFParticlesRealIdx::vz) = v[2]; + + p.idata(ERFParticlesIntIdx::i) = iv[0]; + p.idata(ERFParticlesIntIdx::j) = iv[1]; + p.idata(ERFParticlesIntIdx::k) = iv[2]; + + host_particles.push_back(p); + } + } + + } + + auto& particles = GetParticles(lev); + auto& particle_tile = particles[std::make_pair(mfi.index(), mfi.LocalTileIndex())]; + auto old_size = particle_tile.GetArrayOfStructs().size(); + auto new_size = old_size + host_particles.size(); + particle_tile.resize(new_size); + + Gpu::copy(Gpu::hostToDevice, + host_particles.begin(), + host_particles.end(), + particle_tile.GetArrayOfStructs().begin() + old_size); + } + + return; +} + +/*! Default initialization for hydro particles (ref: AA) */ +void ERFPC::initializeParticlesDefaultHydro(const std::unique_ptr& a_height_ptr) +{ + BL_PROFILE("ERFPC::initializeParticlesDefaultHydro"); + + const int lev = 0; + const Real* dx = Geom(lev).CellSize(); + const Real* plo = Geom(lev).ProbLo(); + + for(MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) + { + const Box& tile_box = mfi.tilebox(); + Gpu::HostVector host_particles; + + if (a_height_ptr) { + + const auto& height = (*a_height_ptr)[mfi]; + const FArrayBox* height_ptr = nullptr; +#ifdef AMREX_USE_GPU + std::unique_ptr hostfab; + if (height.arena()->isManaged() || height.arena()->isDevice()) { + hostfab = std::make_unique(height.box(), height.nComp(), + The_Pinned_Arena()); + Gpu::dtoh_memcpy_async(hostfab->dataPtr(), height.dataPtr(), + height.size()*sizeof(Real)); + Gpu::streamSynchronize(); + height_ptr = hostfab.get(); + } +#else + height_ptr = &height; +#endif + for (IntVect iv = tile_box.smallEnd(); iv <= tile_box.bigEnd(); tile_box.next(iv)) { + // This is a random choice to put them above the ground and let them fall + if (iv[2] == 13) { + Real r[3] = {0.5, 0.5, 0.5}; // this means place at cell center + Real v[3] = {0.0, 0.0, 0.0}; // with 0 initial velocity + + Real x = plo[0] + (iv[0] + r[0])*dx[0]; + Real y = plo[1] + (iv[1] + r[1])*dx[1]; + Real z = (*height_ptr)(iv) + + r[2]*( (*height_ptr)(iv + IntVect(AMREX_D_DECL(0, 0, 1))) + - (*height_ptr)(iv) ); + + ParticleType p; + p.id() = ParticleType::NextID(); + p.cpu() = ParallelDescriptor::MyProc(); + p.pos(0) = x; + p.pos(1) = y; + p.pos(2) = z; + + p.rdata(ERFParticlesRealIdx::vx) = v[0]; + p.rdata(ERFParticlesRealIdx::vy) = v[1]; + p.rdata(ERFParticlesRealIdx::vz) = v[2]; + + p.idata(ERFParticlesIntIdx::i) = iv[0]; + p.idata(ERFParticlesIntIdx::j) = iv[1]; + p.idata(ERFParticlesIntIdx::k) = iv[2]; + + host_particles.push_back(p); + } + } + + } else { + + for (IntVect iv = tile_box.smallEnd(); iv <= tile_box.bigEnd(); tile_box.next(iv)) { + // This is a random choice to put them above the ground and let them fall + if (iv[2] == 23) { + Real r[3] = {0.5, 0.5, 0.5}; // this means place at cell center + Real v[3] = {0.0, 0.0, 0.0}; // with 0 initial velocity + + Real x = plo[0] + (iv[0] + r[0])*dx[0]; + Real y = plo[1] + (iv[1] + r[1])*dx[1]; + Real z = plo[2] + (iv[2] + r[2])*dx[2]; + + ParticleType p; + p.id() = ParticleType::NextID(); + p.cpu() = ParallelDescriptor::MyProc(); + p.pos(0) = x; + p.pos(1) = y; + p.pos(2) = z; + + p.rdata(ERFParticlesRealIdx::vx) = v[0]; + p.rdata(ERFParticlesRealIdx::vy) = v[1]; + p.rdata(ERFParticlesRealIdx::vz) = v[2]; + + p.idata(ERFParticlesIntIdx::i) = iv[0]; + p.idata(ERFParticlesIntIdx::j) = iv[1]; + p.idata(ERFParticlesIntIdx::k) = iv[2]; + + host_particles.push_back(p); + } + } + + } + + auto& particles = GetParticles(lev); + auto& particle_tile = particles[std::make_pair(mfi.index(), mfi.LocalTileIndex())]; + auto old_size = particle_tile.GetArrayOfStructs().size(); + auto new_size = old_size + host_particles.size(); + particle_tile.resize(new_size); + + Gpu::copy(Gpu::hostToDevice, + host_particles.begin(), + host_particles.end(), + particle_tile.GetArrayOfStructs().begin() + old_size); + } +} + +/*! Uniform distribution: the number of particles per grid cell is specified + * by "initial_particles_per_cell", and they are randomly distributed. */ +void ERFPC::initializeParticlesUniformDistribution(const std::unique_ptr& a_height_ptr) +{ + BL_PROFILE("ERFPC::initializeParticlesUniformDistribution"); + + const int lev = 0; + const Real* dx = Geom(lev).CellSize(); + const Real* plo = Geom(lev).ProbLo(); + + for(MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) + { + const Box& tile_box = mfi.tilebox(); + Gpu::HostVector host_particles; + + if (a_height_ptr) { + + const auto& height = (*a_height_ptr)[mfi]; + const FArrayBox* height_ptr = nullptr; +#ifdef AMREX_USE_GPU + std::unique_ptr hostfab; + if (height.arena()->isManaged() || height.arena()->isDevice()) { + hostfab = std::make_unique(height.box(), height.nComp(), + The_Pinned_Arena()); + Gpu::dtoh_memcpy_async(hostfab->dataPtr(), height.dataPtr(), + height.size()*sizeof(Real)); + Gpu::streamSynchronize(); + height_ptr = hostfab.get(); + } +#else + height_ptr = &height; +#endif + for (IntVect iv = tile_box.smallEnd(); iv <= tile_box.bigEnd(); tile_box.next(iv)) { + std::vector rnd_reals(m_ppc_init*3); + amrex::FillRandom(rnd_reals.data(), rnd_reals.size()); + for (int n = 0; n < m_ppc_init; n++) { + Real r[3] = {rnd_reals[3*n], rnd_reals[3*n+1], rnd_reals[3*n+2]}; + Real v[3] = {0.0, 0.0, 0.0}; + + Real x = plo[0] + (iv[0] + r[0])*dx[0]; + Real y = plo[1] + (iv[1] + r[1])*dx[1]; + Real z = (*height_ptr)(iv) + + r[2] * ( (*height_ptr)(iv + IntVect(AMREX_D_DECL(0, 0, 1))) + - (*height_ptr)(iv) ); + + ParticleType p; + p.id() = ParticleType::NextID(); + p.cpu() = ParallelDescriptor::MyProc(); + p.pos(0) = x; + p.pos(1) = y; + p.pos(2) = z; + + p.rdata(ERFParticlesRealIdx::vx) = v[0]; + p.rdata(ERFParticlesRealIdx::vy) = v[1]; + p.rdata(ERFParticlesRealIdx::vz) = v[2]; + + p.idata(ERFParticlesIntIdx::i) = iv[0]; + p.idata(ERFParticlesIntIdx::j) = iv[1]; + p.idata(ERFParticlesIntIdx::k) = iv[2]; + + host_particles.push_back(p); + } + } + + } else { + + for (IntVect iv = tile_box.smallEnd(); iv <= tile_box.bigEnd(); tile_box.next(iv)) { + std::vector rnd_reals(m_ppc_init*3); + amrex::FillRandom(rnd_reals.data(), rnd_reals.size()); + for (int n = 0; n < m_ppc_init; n++) { + Real r[3] = {rnd_reals[3*n], rnd_reals[3*n+1], rnd_reals[3*n+2]}; + Real v[3] = {0.0, 0.0, 0.0}; + + Real x = plo[0] + (iv[0] + r[0])*dx[0]; + Real y = plo[1] + (iv[1] + r[1])*dx[1]; + Real z = plo[2] + (iv[2] + r[2])*dx[2]; + + ParticleType p; + p.id() = ParticleType::NextID(); + p.cpu() = ParallelDescriptor::MyProc(); + p.pos(0) = x; + p.pos(1) = y; + p.pos(2) = z; + + p.rdata(ERFParticlesRealIdx::vx) = v[0]; + p.rdata(ERFParticlesRealIdx::vy) = v[1]; + p.rdata(ERFParticlesRealIdx::vz) = v[2]; + + p.idata(ERFParticlesIntIdx::i) = iv[0]; + p.idata(ERFParticlesIntIdx::j) = iv[1]; + p.idata(ERFParticlesIntIdx::k) = iv[2]; + + host_particles.push_back(p); + } + } + + } + + auto& particles = GetParticles(lev); + auto& particle_tile = particles[std::make_pair(mfi.index(), mfi.LocalTileIndex())]; + auto old_size = particle_tile.GetArrayOfStructs().size(); + auto new_size = old_size + host_particles.size(); + particle_tile.resize(new_size); + + Gpu::copy(Gpu::hostToDevice, + host_particles.begin(), + host_particles.end(), + particle_tile.GetArrayOfStructs().begin() + old_size); + } + + return; +} + +#endif diff --git a/Source/Particles/ERFTracers.cpp b/Source/Particles/ERFTracers.cpp new file mode 100644 index 000000000..f320f4de5 --- /dev/null +++ b/Source/Particles/ERFTracers.cpp @@ -0,0 +1,86 @@ +#include +#include +#include + +#ifdef ERF_USE_PARTICLES + +using namespace amrex; + +/*! Read tracer and hydro particles parameters */ +void ERF::readTracersParams() +{ + ParmParse pp(pp_prefix); + + m_use_tracer_particles = 0; + m_use_hydro_particles = 0; + + pp.query(std::string("use_"+ERFParticleNames::tracers).c_str(), m_use_tracer_particles); + pp.query(std::string("use_"+ERFParticleNames::hydro).c_str(), m_use_hydro_particles); + + if (m_use_tracer_particles) { + particleData.addName(ERFParticleNames::tracers); + } + + if (m_use_hydro_particles) { + particleData.addName(ERFParticleNames::hydro); + } + return; +} + +/*! Initialize tracer and hydro particles */ +void ERF::initializeTracers( ParGDBBase* a_gdb, + const Vector>& a_z_phys_nd ) +{ + auto& namelist_unalloc( particleData.getNamesUnalloc() ); + + for (auto it = namelist_unalloc.begin(); it != namelist_unalloc.end(); ++it) { + + std::string species_name( *it ); + + if (species_name == ERFParticleNames::tracers) { + + AMREX_ASSERT(m_use_tracer_particles); + ERFPC* pc = new ERFPC(a_gdb, ERFParticleNames::tracers); + pc->InitializeParticles(a_z_phys_nd[0]); + amrex::Print() << "Initialized " << pc->TotalNumberOfParticles() << " tracer particles.\n"; + particleData.pushBack(ERFParticleNames::tracers, pc); + + } else if (species_name == ERFParticleNames::hydro) { + + AMREX_ASSERT(m_use_hydro_particles); + ERFPC* pc = new ERFPC(a_gdb, ERFParticleNames::hydro); + pc->InitializeParticles(a_z_phys_nd[0]); + amrex::Print() << "Initialized " << pc->TotalNumberOfParticles() << " hydro particles.\n"; + particleData.pushBack(ERFParticleNames::hydro, pc); + + } + } + + if (m_use_tracer_particles) namelist_unalloc.remove( ERFParticleNames::tracers ); + if (m_use_hydro_particles) namelist_unalloc.remove( ERFParticleNames::hydro ); + + return; +} + +/*! Evolve tracers and hydro particles for one time step*/ +void ERF::evolveTracers( int a_lev, + Real a_dt_lev, + Vector>& a_vars_new, + const Vector>& a_z_phys_nd ) +{ + if (m_use_tracer_particles) { + particleData[ERFParticleNames::tracers]->EvolveParticles( a_lev, + a_dt_lev, + a_vars_new, + a_z_phys_nd ); + } + if (m_use_hydro_particles) { + particleData[ERFParticleNames::hydro]->EvolveParticles( a_lev, + a_dt_lev, + a_vars_new, + a_z_phys_nd ); + } + return; +} + +#endif diff --git a/Source/Particles/HydroPC.H b/Source/Particles/HydroPC.H deleted file mode 100644 index ac8683837..000000000 --- a/Source/Particles/HydroPC.H +++ /dev/null @@ -1,68 +0,0 @@ -#ifndef HYDRO_PC_H_ -#define HYDRO_PC_H_ - -#include - -struct HydroRealIdx -{ - enum { - vx = 0, - vy, vz, - mass, - ncomps - }; -}; - -struct HydroIntIdx -{ - enum { - k = 0, - ncomps - }; -}; - -struct HydroAssignor -{ - template - AMREX_GPU_HOST_DEVICE - amrex::IntVect operator() (P const& p, - amrex::GpuArray const& plo, - amrex::GpuArray const& dxi, - const amrex::Box& domain) const noexcept - { - amrex::IntVect iv( - AMREX_D_DECL(int(amrex::Math::floor((p.pos(0)-plo[0])*dxi[0])), - int(amrex::Math::floor((p.pos(1)-plo[1])*dxi[1])), - p.idata(0))); - iv[0] += domain.smallEnd()[0]; - iv[1] += domain.smallEnd()[1]; - return iv; - } -}; - -class HydroPC - : public amrex::ParticleContainer -{ - -public: - - HydroPC (amrex::ParGDBBase* gdb) - : amrex::ParticleContainer(gdb) - {} - - HydroPC (const amrex::Geometry & geom, - const amrex::DistributionMapping & dmap, - const amrex::BoxArray & ba) - : amrex::ParticleContainer(geom, dmap, ba) - {} - - void InitParticles (); - void InitParticles (const amrex::MultiFab& a_z_height); - - void AdvectWithGravity (int level, amrex::Real dt, bool use_terrain, const amrex::MultiFab& a_z_height); -}; - -#endif diff --git a/Source/Particles/HydroPC.cpp b/Source/Particles/HydroPC.cpp deleted file mode 100644 index bda03159b..000000000 --- a/Source/Particles/HydroPC.cpp +++ /dev/null @@ -1,221 +0,0 @@ -#include "HydroPC.H" - -// #include -#include - -using namespace amrex; - -void -HydroPC:: -InitParticles () -{ - BL_PROFILE("HydroPC::InitParticles"); - - const int lev = 0; - const Real* dx = Geom(lev).CellSize(); - const Real* plo = Geom(lev).ProbLo(); - - for(MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) - { - const Box& tile_box = mfi.tilebox(); - - Gpu::HostVector host_particles; - for (IntVect iv = tile_box.smallEnd(); iv <= tile_box.bigEnd(); tile_box.next(iv)) { - if (iv[2] == 23) { // This is a random choice to put them above the ground and let them fall - Real r[3] = {0.5, 0.5, 0.5}; // this means place at cell center - Real v[3] = {0.0, 0.0, 0.0}; // with 0 initial velocity - - Real x = plo[0] + (iv[0] + r[0])*dx[0]; - Real y = plo[1] + (iv[1] + r[1])*dx[1]; - Real z = plo[2] + (iv[2] + r[2])*dx[2]; - - ParticleType p; - p.id() = ParticleType::NextID(); - p.cpu() = ParallelDescriptor::MyProc(); - p.pos(0) = x; - p.pos(1) = y; - p.pos(2) = z; - - p.rdata(HydroRealIdx::vx) = v[0]; - p.rdata(HydroRealIdx::vy) = v[1]; - p.rdata(HydroRealIdx::vz) = v[2]; - - p.idata(HydroIntIdx::k) = iv[2]; // particles carry their z-index - - host_particles.push_back(p); - } - } - - auto& particles = GetParticles(lev); - auto& particle_tile = particles[std::make_pair(mfi.index(), mfi.LocalTileIndex())]; - auto old_size = particle_tile.GetArrayOfStructs().size(); - auto new_size = old_size + host_particles.size(); - particle_tile.resize(new_size); - - Gpu::copy(Gpu::hostToDevice, - host_particles.begin(), - host_particles.end(), - particle_tile.GetArrayOfStructs().begin() + old_size); - } -} - -void -HydroPC:: -InitParticles (const MultiFab& a_z_height) -{ - BL_PROFILE("HydroPC::InitParticles"); - - const int lev = 0; - const Real* dx = Geom(lev).CellSize(); - const Real* plo = Geom(lev).ProbLo(); - - for(MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) - { - const Box& tile_box = mfi.tilebox(); - const auto& height = a_z_height[mfi]; - const FArrayBox* height_ptr = nullptr; -#ifdef AMREX_USE_GPU - std::unique_ptr hostfab; - if (height.arena()->isManaged() || height.arena()->isDevice()) { - hostfab = std::make_unique(height.box(), height.nComp(), - The_Pinned_Arena()); - Gpu::dtoh_memcpy_async(hostfab->dataPtr(), height.dataPtr(), - height.size()*sizeof(Real)); - Gpu::streamSynchronize(); - height_ptr = hostfab.get(); - } -#else - height_ptr = &height; -#endif - Gpu::HostVector host_particles; - for (IntVect iv = tile_box.smallEnd(); iv <= tile_box.bigEnd(); tile_box.next(iv)) { - if (iv[2] == 13) { // This is a random choice to put them above the ground and let them fall - Real r[3] = {0.5, 0.5, 0.5}; // this means place at cell center - Real v[3] = {0.0, 0.0, 0.0}; // with 0 initial velocity - - Real x = plo[0] + (iv[0] + r[0])*dx[0]; - Real y = plo[1] + (iv[1] + r[1])*dx[1]; - Real z = (*height_ptr)(iv) + r[2]*((*height_ptr)(iv + IntVect(AMREX_D_DECL(0, 0, 1))) - (*height_ptr)(iv)); - - ParticleType p; - p.id() = ParticleType::NextID(); - p.cpu() = ParallelDescriptor::MyProc(); - p.pos(0) = x; - p.pos(1) = y; - p.pos(2) = z; - - p.rdata(HydroRealIdx::vx) = v[0]; - p.rdata(HydroRealIdx::vy) = v[1]; - p.rdata(HydroRealIdx::vz) = v[2]; - - p.idata(HydroIntIdx::k) = iv[2]; // particles carry their z-index - - host_particles.push_back(p); - } - } - - auto& particles = GetParticles(lev); - auto& particle_tile = particles[std::make_pair(mfi.index(), mfi.LocalTileIndex())]; - auto old_size = particle_tile.GetArrayOfStructs().size(); - auto new_size = old_size + host_particles.size(); - particle_tile.resize(new_size); - - Gpu::copy(Gpu::hostToDevice, - host_particles.begin(), - host_particles.end(), - particle_tile.GetArrayOfStructs().begin() + old_size); - } -} - -/* - /brief Uses midpoint method to advance particles using umac. -*/ -void -HydroPC::AdvectWithGravity (int lev, Real dt, bool use_terrain, const MultiFab& a_z_height) -{ - BL_PROFILE("HydroPC::AdvectWithGravity()"); - AMREX_ASSERT(lev >= 0 && lev < GetParticles().size()); - - const auto strttime = amrex::second(); - const Geometry& geom = m_gdb->Geom(lev); - const Box& domain = geom.Domain(); - const auto plo = geom.ProbLoArray(); - const auto dx = geom.CellSizeArray(); - const auto dxi = geom.InvCellSizeArray(); - -#ifdef AMREX_USE_OMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for (ParIterType pti(*this, lev); pti.isValid(); ++pti) - { - int grid = pti.index(); - auto& ptile = ParticlesAt(lev, pti); - auto& aos = ptile.GetArrayOfStructs(); - const int n = aos.numParticles(); - auto *p_pbox = aos().data(); - - auto zheight = use_terrain ? a_z_height[grid].array() : Array4{}; - - amrex::ParallelFor(n, [=] AMREX_GPU_DEVICE (int i) - { - ParticleType& p = p_pbox[i]; - if (p.id() <= 0) { return; } - - ParticleReal v = p.rdata(HydroRealIdx::vz); - - // Define acceleration to be (gravity minus drag) where drag is defined - // such the particles will reach a terminal velocity of 5.0 (totally arbitrary) - ParticleReal terminal_vel = 5.0; - ParticleReal grav = CONST_GRAV; - ParticleReal drag = CONST_GRAV * (v * v) / (terminal_vel*terminal_vel); - - ParticleReal half_dt = 0.5 * dt; - - // Update the particle velocity over first half of step (dt/2) - p.rdata(HydroRealIdx::vz) -= (grav - drag) * half_dt; - - // Update the particle position over (dt) - p.pos(2) += static_cast(ParticleReal(0.5)*dt*p.rdata(HydroRealIdx::vz)); - - // Update the particle velocity over second half of step (dt/2) - p.rdata(HydroRealIdx::vz) -= (grav - drag) * half_dt; - - // also update z-coordinate here - IntVect iv( - AMREX_D_DECL(int(amrex::Math::floor((p.pos(0)-plo[0])*dxi[0])), - int(amrex::Math::floor((p.pos(1)-plo[1])*dxi[1])), - p.idata(0))); - iv[0] += domain.smallEnd()[0]; - iv[1] += domain.smallEnd()[1]; - ParticleReal zlo, zhi; - if (use_terrain) { - zlo = zheight(iv[0], iv[1], iv[2]); - zhi = zheight(iv[0], iv[1], iv[2]+1); - } else { - zlo = iv[2] * dx[2]; - zhi = (iv[2]+1) * dx[2]; - } - if (p.pos(2) > zhi) { // need to be careful here - p.idata(0) += 1; - } else if (p.pos(2) <= zlo) { - p.idata(0) -= 1; - } - }); - } - - if (m_verbose > 1) - { - auto stoptime = amrex::second() - strttime; - -#ifdef AMREX_LAZY - Lazy::QueueReduction( [=] () mutable { -#endif - ParallelReduce::Max(stoptime, ParallelContext::IOProcessorNumberSub(), - ParallelContext::CommunicatorSub()); - - amrex::Print() << "HydroParticleContainer::AdvectWithGravity() time: " << stoptime << '\n'; -#ifdef AMREX_LAZY - }); -#endif - } -} diff --git a/Source/Particles/Make.package b/Source/Particles/Make.package index 0bacbcdd5..d00e0269d 100644 --- a/Source/Particles/Make.package +++ b/Source/Particles/Make.package @@ -1,6 +1,6 @@ -CEXE_sources += HydroPC.cpp -CEXE_sources += TracerPC.cpp +CEXE_sources += ERFTracers.cpp +CEXE_sources += ERFPCInitializations.cpp +CEXE_sources += ERFPCEvolve.cpp CEXE_headers += ParticleData.H -CEXE_headers += TracerPC.H -CEXE_headers += HydroPC.H +CEXE_headers += ERFPC.H diff --git a/Source/Particles/ParticleData.H b/Source/Particles/ParticleData.H index 83e4ce518..229b7da70 100644 --- a/Source/Particles/ParticleData.H +++ b/Source/Particles/ParticleData.H @@ -1,134 +1,201 @@ #ifndef _PARTICLE_DATA_H_ #define _PARTICLE_DATA_H_ +#ifdef ERF_USE_PARTICLES + +#include +#include +#include #include #include -#include -#include #include #include #include -#include -#include +#include /** * Container holding many of the particle-related data and options */ -struct ParticleData { - public: - void init_particle_params() - { - amrex::ParmParse pp(pp_prefix); +typedef std::map ParticleSpeciesMap; +typedef std::vector ParticlesNamesVector; +typedef std::list ParticlesNamesList; + +class ParticleData +{ + public: + + /*! Constructor */ + ParticleData() + { + BL_PROFILE("ParticleData::ParticleData()"); + m_particle_species.clear(); + m_namelist.clear(); + m_namelist_unalloc.clear(); + } - use_tracer_particles = 0; - pp.query("use_tracer_particles", use_tracer_particles); + /*! Destructor */ + ~ParticleData() + { + BL_PROFILE("ParticleData::~ParticleData()"); + for (ParticlesNamesVector::size_type i = 0; i < m_namelist.size(); i++) { + auto particles( m_particle_species[m_namelist[i]] ); + delete particles; + } + m_particle_species.clear(); + m_namelist.clear(); + m_namelist_unalloc.clear(); + } - use_hydro_particles = 0; - pp.query("use_hydro_particles", use_hydro_particles); - } + /*! Write checkpoint files */ + void Checkpoint( const std::string& a_fname ) + { + BL_PROFILE("ParticleData::Checkpoint()"); + for (ParticlesNamesVector::size_type i = 0; i < m_namelist.size(); i++) { + auto name( m_namelist[i] ); + auto particles( m_particle_species[name] ); + particles->Checkpoint( a_fname, name, true, particles->varNames() ); + } + } - void init_particles(amrex::ParGDBBase* a_gdb, const amrex::Vector>& z_phys_nd) + /*! Read from restart file */ + void Restart( amrex::ParGDBBase* a_gdb, const std::string& a_fname ) + { + BL_PROFILE("ParticleData::Restart()"); + AMREX_ASSERT(isEmpty()); + for (auto it = m_namelist_unalloc.begin(); it != m_namelist_unalloc.end(); ++it) { + std::string species_name( *it ); + ERFPC* pc = new ERFPC( a_gdb, species_name ); + pc->Restart(a_fname, species_name ); + pushBack( species_name, pc ); + } + m_namelist_unalloc.clear(); + } - { - // Initialize tracer particles - if (use_tracer_particles) { - tracer_particles = std::make_unique(a_gdb); - if (z_phys_nd[0] != nullptr) { - tracer_particles->InitParticles(*z_phys_nd[0]); - } else { - tracer_particles->InitParticles(); + /*! Redistribute/rebalance particles data */ + inline void Redistribute() + { + BL_PROFILE("ParticleData::Redistribute()"); + for (ParticlesNamesVector::size_type i = 0; i < m_namelist.size(); i++) { + m_particle_species[m_namelist[i]]->Redistribute(); } - tracer_particles->Redistribute(); - amrex::Print() << "Initialized " << tracer_particles->TotalNumberOfParticles() << " tracer particles." << std::endl; } - // Initialize hydro particles - if (use_hydro_particles) { - hydro_particles = std::make_unique(a_gdb); - if (z_phys_nd[0] != nullptr) { - hydro_particles->InitParticles(*z_phys_nd[0]); + /*! Get species of a given name */ + inline ERFPC* GetSpecies( const std::string& a_name ) + { + BL_PROFILE("ParticleData::GetSpecies()"); + ParticleSpeciesMap::iterator it (m_particle_species.find(a_name)); + if (it == m_particle_species.end()) { + amrex::Print() << "ERROR: unable to fine particle species with name \"" + << a_name << "\"!"; + return nullptr; } else { - hydro_particles->InitParticles(); + return it->second; } - hydro_particles->Redistribute(); - amrex::Print() << "Initialized " << hydro_particles->TotalNumberOfParticles() << " hydro particles." << std::endl; } - } - void Checkpoint(const std::string& filename) - { - if (use_tracer_particles) { - tracer_particles->Checkpoint(filename, "tracer_particles", true, tracer_particle_varnames); - } - if (use_hydro_particles) { - hydro_particles->Checkpoint(filename, "hydro_particles", true, hydro_particle_varnames); - } - } - - void Restart(amrex::ParGDBBase* a_gdb, std::string& restart_file) - { - if (use_tracer_particles) { - // tracer_particles = std::make_unique(Geom(0), dmap[0], grids[0]); - tracer_particles = std::make_unique(a_gdb); - std::string tracer_file("tracer_particles"); - tracer_particles->Restart(restart_file, tracer_file); - } - if (use_hydro_particles) { - // hydro_particles = std::make_unique(Geom(0), dmap[0], grids[0]); - hydro_particles = std::make_unique(a_gdb); - std::string hydro_file("hydro_particles"); - hydro_particles->Restart(restart_file, hydro_file); + /*! accessor */ + inline ERFPC* operator[] ( const std::string& a_name ) + { + BL_PROFILE("ParticleData::operator[]"); + ParticleSpeciesMap::iterator it (m_particle_species.find(a_name)); + if (it == m_particle_species.end()) { + amrex::Print() << "ERROR: unable to fine particle species with name \"" + << a_name << "\"!"; + return nullptr; + } else { + return it->second; + } } - } - - void advance_particles(int lev, amrex::Real dt_lev, amrex::Vector>& vars_new, - const amrex::Vector>& z_phys_nd) - { - // Update tracer particles - if (use_tracer_particles) { - amrex::MultiFab* mac_vel = &vars_new[lev][Vars::xvel]; - bool use_terrain = (z_phys_nd[lev] != nullptr); - if (use_terrain) { - tracer_particles->AdvectWithUmac(mac_vel, lev, dt_lev, use_terrain, *z_phys_nd[lev]); + + /*! Get species of a given name (const version) */ + inline const ERFPC* GetSpecies( const std::string& a_name ) const + { + BL_PROFILE("ParticleData::GetSpecies()"); + ParticleSpeciesMap::const_iterator it (m_particle_species.find(a_name)); + if (it == m_particle_species.end()) { + amrex::Print() << "ERROR: unable to fine particle species with name \"" + << a_name << "\"!"; + return nullptr; } else { - amrex::MultiFab dummy; - tracer_particles->AdvectWithUmac(mac_vel, lev, dt_lev, use_terrain, dummy); + return it->second; } } - // Update hydro particles - if (use_hydro_particles) { - bool use_terrain = (z_phys_nd[lev] != nullptr); - if (use_terrain) { - hydro_particles->AdvectWithGravity(lev, dt_lev, use_terrain, *z_phys_nd[lev]); + /*! accessor */ + inline const ERFPC* operator[] ( const std::string& a_name ) const + { + BL_PROFILE("ParticleData::operator[]"); + ParticleSpeciesMap::const_iterator it (m_particle_species.find(a_name)); + if (it == m_particle_species.end()) { + amrex::Print() << "ERROR: unable to fine particle species with name \"" + << a_name << "\"!"; + return nullptr; } else { - amrex::MultiFab dummy; - hydro_particles->AdvectWithGravity(lev, dt_lev, use_terrain, dummy); + return it->second; } } - } - void Redistribute() - { - if (use_tracer_particles) { - tracer_particles->Redistribute(); + /*! Add a particle species to this container */ + inline void pushBack(const std::string& a_name, + ERFPC* const a_pc ) + { + BL_PROFILE("ParticleData::pushBack()"); + AMREX_ASSERT(!contains(a_name)); + m_particle_species[a_name] = a_pc; + m_namelist.push_back(a_name); } - if (use_hydro_particles) { - hydro_particles->Redistribute(); + + /*! Add a name; particle container will be initialized later */ + inline void addName(const std::string& a_name ) + { + BL_PROFILE("ParticleData::addName()"); + m_namelist_unalloc.push_back(a_name); } - } - std::string pp_prefix {"erf"}; + /*! Returns list of names of particle species */ + inline const ParticlesNamesVector& getNames() const + { + BL_PROFILE("ParticleData::getNames()"); + return m_namelist; + } - int use_tracer_particles; - std::unique_ptr tracer_particles; - amrex::Vector tracer_particle_varnames = {AMREX_D_DECL("xvel", "yvel", "zvel")}; + /*! Returns list of names of particle species that are unallocated */ + inline ParticlesNamesList& getNamesUnalloc() + { + BL_PROFILE("ParticleData::getNamesUnalloc()"); + return m_namelist_unalloc; + } - int use_hydro_particles; - std::unique_ptr hydro_particles; - amrex::Vector hydro_particle_varnames = {AMREX_D_DECL("xvel", "yvel", "zvel"), "mass"}; + /*! queries if containe has species of a certain name */ + inline bool contains( const std::string& a_name ) const + { + BL_PROFILE("ParticleData::contains()"); + ParticleSpeciesMap::const_iterator it (m_particle_species.find(a_name)); + return (it != m_particle_species.end()); + } + + /*! query if container is empty */ + inline bool isEmpty() const + { + BL_PROFILE("ParticleData::isEmpty()"); + return (m_particle_species.size() == 0); + } + + + private: + + /*! Vector of all particle species */ + ParticleSpeciesMap m_particle_species; + /*! Vector of particle species names */ + ParticlesNamesVector m_namelist; + /*! List with names of unallocated species */ + ParticlesNamesList m_namelist_unalloc; }; + +#endif #endif + diff --git a/Source/Particles/TracerPC.H b/Source/Particles/TracerPC.H deleted file mode 100644 index f3c0fc740..000000000 --- a/Source/Particles/TracerPC.H +++ /dev/null @@ -1,68 +0,0 @@ -#ifndef TRACER_PC_H_ -#define TRACER_PC_H_ - -#include - -struct TracerRealIdx -{ - enum { - vx = 0, - vy, vz, - ncomps - }; -}; - -struct TracerIntIdx -{ - enum { - k = 0, - ncomps - }; -}; - -struct TracerAssignor -{ - template - AMREX_GPU_HOST_DEVICE - amrex::IntVect operator() (P const& p, - amrex::GpuArray const& plo, - amrex::GpuArray const& dxi, - const amrex::Box& domain) const noexcept - { - amrex::IntVect iv( - AMREX_D_DECL(int(amrex::Math::floor((p.pos(0)-plo[0])*dxi[0])), - int(amrex::Math::floor((p.pos(1)-plo[1])*dxi[1])), - p.idata(0))); - iv[0] += domain.smallEnd()[0]; - iv[1] += domain.smallEnd()[1]; - return iv; - } -}; - -class TracerPC - : public amrex::ParticleContainer -{ - -public: - - TracerPC (amrex::ParGDBBase* gdb) - : amrex::ParticleContainer(gdb) - {} - - TracerPC (const amrex::Geometry & geom, - const amrex::DistributionMapping & dmap, - const amrex::BoxArray & ba) - : amrex::ParticleContainer(geom, dmap, ba) - {} - - void InitParticles (); - void InitParticles (const amrex::MultiFab& a_z_height); - - void AdvectWithUmac (amrex::MultiFab* umac, int level, amrex::Real dt, bool use_terrain, - amrex::MultiFab& a_z_height); -}; - -#endif diff --git a/Source/Particles/TracerPC.cpp b/Source/Particles/TracerPC.cpp deleted file mode 100644 index 8d2c20b1d..000000000 --- a/Source/Particles/TracerPC.cpp +++ /dev/null @@ -1,261 +0,0 @@ -#include "TracerPC.H" - -#include - -using namespace amrex; - -void -TracerPC:: -InitParticles () -{ - BL_PROFILE("TracerPC::InitParticles"); - - const int lev = 0; - const Real* dx = Geom(lev).CellSize(); - const Real* plo = Geom(lev).ProbLo(); - - for(MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) - { - const Box& tile_box = mfi.tilebox(); - Gpu::HostVector host_particles; - for (IntVect iv = tile_box.smallEnd(); iv <= tile_box.bigEnd(); tile_box.next(iv)) { - if (iv[0] == 3) { - Real r[3] = {0.5, 0.5, 0.5}; // this means place at cell center - Real v[3] = {0.0, 0.0, 0.0}; // with 0 initial velocity - - Real x = plo[0] + (iv[0] + r[0])*dx[0]; - Real y = plo[1] + (iv[1] + r[1])*dx[1]; - Real z = plo[2] + (iv[2] + r[2])*dx[2]; - - ParticleType p; - p.id() = ParticleType::NextID(); - p.cpu() = ParallelDescriptor::MyProc(); - p.pos(0) = x; - p.pos(1) = y; - p.pos(2) = z; - - p.rdata(TracerRealIdx::vx) = v[0]; - p.rdata(TracerRealIdx::vy) = v[1]; - p.rdata(TracerRealIdx::vz) = v[2]; - - p.idata(TracerIntIdx::k) = iv[2]; // particles carry their z-index - - host_particles.push_back(p); - } - } - - auto& particles = GetParticles(lev); - auto& particle_tile = particles[std::make_pair(mfi.index(), mfi.LocalTileIndex())]; - auto old_size = particle_tile.GetArrayOfStructs().size(); - auto new_size = old_size + host_particles.size(); - particle_tile.resize(new_size); - - Gpu::copy(Gpu::hostToDevice, - host_particles.begin(), - host_particles.end(), - particle_tile.GetArrayOfStructs().begin() + old_size); - } -} - -void -TracerPC:: -InitParticles (const MultiFab& a_z_height) -{ - BL_PROFILE("TracerPC::InitParticles"); - - const int lev = 0; - const Real* dx = Geom(lev).CellSize(); - const Real* plo = Geom(lev).ProbLo(); - - for(MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) - { - const Box& tile_box = mfi.tilebox(); - const auto& height = a_z_height[mfi]; - const FArrayBox* height_ptr = nullptr; -#ifdef AMREX_USE_GPU - std::unique_ptr hostfab; - if (height.arena()->isManaged() || height.arena()->isDevice()) { - hostfab = std::make_unique(height.box(), height.nComp(), - The_Pinned_Arena()); - Gpu::dtoh_memcpy_async(hostfab->dataPtr(), height.dataPtr(), - height.size()*sizeof(Real)); - Gpu::streamSynchronize(); - height_ptr = hostfab.get(); - } -#else - height_ptr = &height; -#endif - Gpu::HostVector host_particles; - for (IntVect iv = tile_box.smallEnd(); iv <= tile_box.bigEnd(); tile_box.next(iv)) { - if (iv[0] == 3) { - Real r[3] = {0.5, 0.5, 0.5}; // this means place at cell center - Real v[3] = {0.0, 0.0, 0.0}; // with 0 initial velocity - - Real x = plo[0] + (iv[0] + r[0])*dx[0]; - Real y = plo[1] + (iv[1] + r[1])*dx[1]; - Real z = (*height_ptr)(iv) + r[2]*((*height_ptr)(iv + IntVect(AMREX_D_DECL(0, 0, 1))) - (*height_ptr)(iv)); - - ParticleType p; - p.id() = ParticleType::NextID(); - p.cpu() = ParallelDescriptor::MyProc(); - p.pos(0) = x; - p.pos(1) = y; - p.pos(2) = z; - - p.rdata(TracerRealIdx::vx) = v[0]; - p.rdata(TracerRealIdx::vy) = v[1]; - p.rdata(TracerRealIdx::vz) = v[2]; - - p.idata(TracerIntIdx::k) = iv[2]; // particles carry their z-index - - host_particles.push_back(p); - } - } - - auto& particles = GetParticles(lev); - auto& particle_tile = particles[std::make_pair(mfi.index(), mfi.LocalTileIndex())]; - auto old_size = particle_tile.GetArrayOfStructs().size(); - auto new_size = old_size + host_particles.size(); - particle_tile.resize(new_size); - - Gpu::copy(Gpu::hostToDevice, - host_particles.begin(), - host_particles.end(), - particle_tile.GetArrayOfStructs().begin() + old_size); - } -} - -/* - /brief Uses midpoint method to advance particles using umac. -*/ -void -TracerPC::AdvectWithUmac (MultiFab* umac, int lev, Real dt, bool use_terrain, MultiFab& a_z_height) -{ - BL_PROFILE("TracerPC::AdvectWithUmac()"); - AMREX_ASSERT(OK(lev, lev, umac[0].nGrow()-1)); - AMREX_ASSERT(lev >= 0 && lev < GetParticles().size()); - - AMREX_D_TERM(AMREX_ASSERT(umac[0].nGrow() >= 1);, - AMREX_ASSERT(umac[1].nGrow() >= 1);, - AMREX_ASSERT(umac[2].nGrow() >= 1);); - - AMREX_D_TERM(AMREX_ASSERT(!umac[0].contains_nan());, - AMREX_ASSERT(!umac[1].contains_nan());, - AMREX_ASSERT(!umac[2].contains_nan());); - - const auto strttime = amrex::second(); - const Geometry& geom = m_gdb->Geom(lev); - const Box& domain = geom.Domain(); - const auto plo = geom.ProbLoArray(); - const auto dxi = geom.InvCellSizeArray(); - const auto dx = geom.CellSizeArray(); - - Vector > raii_umac(AMREX_SPACEDIM); - Vector umac_pointer(AMREX_SPACEDIM); - if (OnSameGrids(lev, umac[0])) - { - for (int i = 0; i < AMREX_SPACEDIM; i++) { - umac_pointer[i] = &umac[i]; - } - } - else - { - for (int i = 0; i < AMREX_SPACEDIM; i++) - { - int ng = umac[i].nGrow(); - raii_umac[i] = std::make_unique - (amrex::convert(m_gdb->ParticleBoxArray(lev), IntVect::TheDimensionVector(i)), - m_gdb->ParticleDistributionMap(lev), umac[i].nComp(), ng); - umac_pointer[i] = raii_umac[i].get(); - umac_pointer[i]->ParallelCopy(umac[i],0,0,umac[i].nComp(),ng,ng); - } - } - - for (int ipass = 0; ipass < 2; ipass++) - { -#ifdef AMREX_USE_OMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for (ParIterType pti(*this, lev); pti.isValid(); ++pti) - { - int grid = pti.index(); - auto& ptile = ParticlesAt(lev, pti); - auto& aos = ptile.GetArrayOfStructs(); - const int n = aos.numParticles(); - auto *p_pbox = aos().data(); - const FArrayBox* fab[AMREX_SPACEDIM] = { AMREX_D_DECL(&((*umac_pointer[0])[grid]), - &((*umac_pointer[1])[grid]), - &((*umac_pointer[2])[grid])) }; - - //array of these pointers to pass to the GPU - amrex::GpuArray, AMREX_SPACEDIM> - const umacarr {{AMREX_D_DECL((*fab[0]).array(), - (*fab[1]).array(), - (*fab[2]).array() )}}; - - auto zheight = use_terrain ? a_z_height[grid].array() : Array4{}; - - amrex::ParallelFor(n, - [=] AMREX_GPU_DEVICE (int i) - { - ParticleType& p = p_pbox[i]; - if (p.id() <= 0) { return; } - ParticleReal v[AMREX_SPACEDIM]; - mac_interpolate(p, plo, dxi, umacarr, v); - if (ipass == 0) - { - for (int dim=0; dim < AMREX_SPACEDIM; dim++) - { - p.rdata(dim) = p.pos(dim); - p.pos(dim) += static_cast(ParticleReal(0.5)*dt*v[dim]); - } - } - else - { - for (int dim=0; dim < AMREX_SPACEDIM; dim++) - { - p.pos(dim) = p.rdata(dim) + static_cast(dt*v[dim]); - p.rdata(dim) = v[dim]; - } - - // also update z-coordinate here - IntVect iv( - AMREX_D_DECL(int(amrex::Math::floor((p.pos(0)-plo[0])*dxi[0])), - int(amrex::Math::floor((p.pos(1)-plo[1])*dxi[1])), - p.idata(0))); - iv[0] += domain.smallEnd()[0]; - iv[1] += domain.smallEnd()[1]; - ParticleReal zlo, zhi; - if (use_terrain) { - zlo = zheight(iv[0], iv[1], iv[2]); - zhi = zheight(iv[0], iv[1], iv[2]+1); - } else { - zlo = iv[2] * dx[2]; - zhi = (iv[2]+1) * dx[2]; - } - if (p.pos(2) > zhi) { // need to be careful here - p.idata(0) += 1; - } else if (p.pos(2) <= zlo) { - p.idata(0) -= 1; - } - } - }); - } - } - - if (m_verbose > 1) - { - auto stoptime = amrex::second() - strttime; - -#ifdef AMREX_LAZY - Lazy::QueueReduction( [=] () mutable { -#endif - ParallelReduce::Max(stoptime, ParallelContext::IOProcessorNumberSub(), - ParallelContext::CommunicatorSub()); - - amrex::Print() << "TracerParticleContainer::AdvectWithUmac() time: " << stoptime << '\n'; -#ifdef AMREX_LAZY - }); -#endif - } -} diff --git a/Source/TimeIntegration/ERF_Advance.cpp b/Source/TimeIntegration/ERF_Advance.cpp index b7f987cbb..17391a0f5 100644 --- a/Source/TimeIntegration/ERF_Advance.cpp +++ b/Source/TimeIntegration/ERF_Advance.cpp @@ -134,6 +134,6 @@ ERF::Advance (int lev, Real time, Real dt_lev, int /*iteration*/, int /*ncycle*/ #endif #ifdef ERF_USE_PARTICLES - particleData.advance_particles(lev, dt_lev, vars_new, z_phys_nd); + evolveTracers( lev, dt_lev, vars_new, z_phys_nd ); #endif } From 5efefd43534d14f7af7fc5fa1ca892e75ea38e50 Mon Sep 17 00:00:00 2001 From: "David J. Gardner" Date: Thu, 4 Jan 2024 15:23:09 -0800 Subject: [PATCH 2/2] white space fixes --- Source/ERF.H | 8 ++--- Source/Particles/ERFPC.H | 40 +++++++++++------------ Source/Particles/ERFPCEvolve.cpp | 10 +++--- Source/Particles/ERFPCInitializations.cpp | 18 +++++----- Source/Particles/ERFTracers.cpp | 12 +++---- Source/Particles/ParticleData.H | 28 ++++++++-------- 6 files changed, 58 insertions(+), 58 deletions(-) diff --git a/Source/ERF.H b/Source/ERF.H index 7170ca0c9..c7dcb83fa 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -678,12 +678,12 @@ private: void readTracersParams(); /*! Initialize tracer and hydro particles */ - void initializeTracers ( amrex::ParGDBBase*, + void initializeTracers ( amrex::ParGDBBase*, const amrex::Vector>& ); /*! Evolve tracers and hydro particles */ - void evolveTracers( int, - amrex::Real, - amrex::Vector>&, + void evolveTracers( int, + amrex::Real, + amrex::Vector>&, const amrex::Vector>& ); #endif diff --git a/Source/Particles/ERFPC.H b/Source/Particles/ERFPC.H index 6a8be3542..2de337549 100644 --- a/Source/Particles/ERFPC.H +++ b/Source/Particles/ERFPC.H @@ -46,12 +46,12 @@ struct ERFParticlesAssignor template AMREX_GPU_HOST_DEVICE amrex::IntVect operator() ( P const& p, - amrex::GpuArray const& plo, + amrex::GpuArray const& plo, amrex::GpuArray const& dxi, const amrex::Box& domain ) const noexcept { /* TODO: What is this about? */ - amrex::IntVect iv( + amrex::IntVect iv( AMREX_D_DECL( int(amrex::Math::floor((p.pos(0)-plo[0])*dxi[0])), int(amrex::Math::floor((p.pos(1)-plo[1])*dxi[1])), p.idata(ERFParticlesIntIdx::k) ) ); @@ -72,13 +72,13 @@ class ERFPC : public amrex::ParticleContainer< ERFParticlesRealIdx::ncomps, /*! Constructor */ ERFPC ( amrex::ParGDBBase* a_gdb, - const std::string& a_name = "particles" ) + const std::string& a_name = "particles" ) : amrex::ParticleContainer< ERFParticlesRealIdx::ncomps, ERFParticlesIntIdx::ncomps, 0, 0, amrex::DefaultAllocator, - ERFParticlesAssignor> (a_gdb) + ERFParticlesAssignor> (a_gdb) { BL_PROFILE("ERFPCPC::ERFPC()"); m_name = a_name; @@ -89,7 +89,7 @@ class ERFPC : public amrex::ParticleContainer< ERFParticlesRealIdx::ncomps, ERFPC ( const amrex::Geometry& a_geom, const amrex::DistributionMapping& a_dmap, const amrex::BoxArray& a_ba, - const std::string& a_name = "particles" ) + const std::string& a_name = "particles" ) : amrex::ParticleContainer< ERFParticlesRealIdx::ncomps, ERFParticlesIntIdx::ncomps, 0, @@ -106,9 +106,9 @@ class ERFPC : public amrex::ParticleContainer< ERFParticlesRealIdx::ncomps, virtual void InitializeParticles(const std::unique_ptr& a_ptr = nullptr); /*! Evolve particles for one time step */ - virtual void EvolveParticles( int, - amrex::Real, - amrex::Vector>&, + virtual void EvolveParticles( int, + amrex::Real, + amrex::Vector>&, const amrex::Vector>& ); /*! Get real-type particle attribute names */ @@ -119,16 +119,16 @@ class ERFPC : public amrex::ParticleContainer< ERFParticlesRealIdx::ncomps, } /*! Specify if particles should advect with flow */ - inline void setAdvectWithFlow (bool a_flag) - { + inline void setAdvectWithFlow (bool a_flag) + { BL_PROFILE("ERFPCPC::setAdvectWithFlow()"); - m_advect_w_flow = a_flag; + m_advect_w_flow = a_flag; } /*! Specify if particles fall under gravity */ - inline void setAdvectWithGravity (bool a_flag) - { + inline void setAdvectWithGravity (bool a_flag) + { BL_PROFILE("ERFPCPC::setAdvectWithGravity()"); - m_advect_w_gravity = a_flag; + m_advect_w_gravity = a_flag; } protected: @@ -137,7 +137,7 @@ class ERFPC : public amrex::ParticleContainer< ERFParticlesRealIdx::ncomps, bool m_advect_w_gravity; /*!< advect under gravitational force */ std::string m_name; /*!< name of this particle species */ - + std::string m_initialization_type; /*!< initial particle distribution type */ int m_ppc_init; /*!< initial number of particles per cell */ @@ -148,14 +148,14 @@ class ERFPC : public amrex::ParticleContainer< ERFParticlesRealIdx::ncomps, void initializeParticlesUniformDistribution(const std::unique_ptr& a_ptr = nullptr); /*! Uses midpoint method to advance particles using flow velocity. */ - void AdvectWithFlow( amrex::MultiFab*, - int, - amrex::Real, + void AdvectWithFlow( amrex::MultiFab*, + int, + amrex::Real, const std::unique_ptr& ); /*! Uses midpoint method to advance particles falling under gravity. */ - void AdvectWithGravity( int, - amrex::Real, + void AdvectWithGravity( int, + amrex::Real, const std::unique_ptr& ); private: diff --git a/Source/Particles/ERFPCEvolve.cpp b/Source/Particles/ERFPCEvolve.cpp index c46a759d5..2d7270d2e 100644 --- a/Source/Particles/ERFPCEvolve.cpp +++ b/Source/Particles/ERFPCEvolve.cpp @@ -28,9 +28,9 @@ void ERFPC::EvolveParticles( int a_lev, } /*! Uses midpoint method to advance particles using flow velocity. */ -void ERFPC::AdvectWithFlow( MultiFab* a_umac, - int a_lev, - Real a_dt, +void ERFPC::AdvectWithFlow( MultiFab* a_umac, + int a_lev, + Real a_dt, const std::unique_ptr& a_z_height ) { BL_PROFILE("ERFPCPC::AdvectWithUmac()"); @@ -163,8 +163,8 @@ void ERFPC::AdvectWithFlow( MultiFab* a_umac, } } -void ERFPC::AdvectWithGravity( int a_lev, - Real a_dt, +void ERFPC::AdvectWithGravity( int a_lev, + Real a_dt, const std::unique_ptr& a_z_height ) { BL_PROFILE("ERFPC::AdvectWithGravity()"); diff --git a/Source/Particles/ERFPCInitializations.cpp b/Source/Particles/ERFPCInitializations.cpp index 3efc4a426..45626dcac 100644 --- a/Source/Particles/ERFPCInitializations.cpp +++ b/Source/Particles/ERFPCInitializations.cpp @@ -97,8 +97,8 @@ void ERFPC::initializeParticlesDefaultTracersWoA(const std::unique_ptr& a_height_ptr) { @@ -316,8 +316,8 @@ void ERFPC::initializeParticlesUniformDistribution(const std::unique_ptr>& a_z_phys_nd ) { if (m_use_tracer_particles) { - particleData[ERFParticleNames::tracers]->EvolveParticles( a_lev, - a_dt_lev, - a_vars_new, + particleData[ERFParticleNames::tracers]->EvolveParticles( a_lev, + a_dt_lev, + a_vars_new, a_z_phys_nd ); } if (m_use_hydro_particles) { - particleData[ERFParticleNames::hydro]->EvolveParticles( a_lev, - a_dt_lev, - a_vars_new, + particleData[ERFParticleNames::hydro]->EvolveParticles( a_lev, + a_dt_lev, + a_vars_new, a_z_phys_nd ); } return; diff --git a/Source/Particles/ParticleData.H b/Source/Particles/ParticleData.H index 229b7da70..df508d777 100644 --- a/Source/Particles/ParticleData.H +++ b/Source/Particles/ParticleData.H @@ -84,12 +84,12 @@ class ParticleData } /*! Get species of a given name */ - inline ERFPC* GetSpecies( const std::string& a_name ) + inline ERFPC* GetSpecies( const std::string& a_name ) { BL_PROFILE("ParticleData::GetSpecies()"); ParticleSpeciesMap::iterator it (m_particle_species.find(a_name)); if (it == m_particle_species.end()) { - amrex::Print() << "ERROR: unable to fine particle species with name \"" + amrex::Print() << "ERROR: unable to fine particle species with name \"" << a_name << "\"!"; return nullptr; } else { @@ -103,7 +103,7 @@ class ParticleData BL_PROFILE("ParticleData::operator[]"); ParticleSpeciesMap::iterator it (m_particle_species.find(a_name)); if (it == m_particle_species.end()) { - amrex::Print() << "ERROR: unable to fine particle species with name \"" + amrex::Print() << "ERROR: unable to fine particle species with name \"" << a_name << "\"!"; return nullptr; } else { @@ -117,7 +117,7 @@ class ParticleData BL_PROFILE("ParticleData::GetSpecies()"); ParticleSpeciesMap::const_iterator it (m_particle_species.find(a_name)); if (it == m_particle_species.end()) { - amrex::Print() << "ERROR: unable to fine particle species with name \"" + amrex::Print() << "ERROR: unable to fine particle species with name \"" << a_name << "\"!"; return nullptr; } else { @@ -126,12 +126,12 @@ class ParticleData } /*! accessor */ - inline const ERFPC* operator[] ( const std::string& a_name ) const + inline const ERFPC* operator[] ( const std::string& a_name ) const { BL_PROFILE("ParticleData::operator[]"); ParticleSpeciesMap::const_iterator it (m_particle_species.find(a_name)); if (it == m_particle_species.end()) { - amrex::Print() << "ERROR: unable to fine particle species with name \"" + amrex::Print() << "ERROR: unable to fine particle species with name \"" << a_name << "\"!"; return nullptr; } else { @@ -157,21 +157,21 @@ class ParticleData } /*! Returns list of names of particle species */ - inline const ParticlesNamesVector& getNames() const - { + inline const ParticlesNamesVector& getNames() const + { BL_PROFILE("ParticleData::getNames()"); - return m_namelist; + return m_namelist; } /*! Returns list of names of particle species that are unallocated */ - inline ParticlesNamesList& getNamesUnalloc() - { + inline ParticlesNamesList& getNamesUnalloc() + { BL_PROFILE("ParticleData::getNamesUnalloc()"); - return m_namelist_unalloc; + return m_namelist_unalloc; } /*! queries if containe has species of a certain name */ - inline bool contains( const std::string& a_name ) const + inline bool contains( const std::string& a_name ) const { BL_PROFILE("ParticleData::contains()"); ParticleSpeciesMap::const_iterator it (m_particle_species.find(a_name)); @@ -189,7 +189,7 @@ class ParticleData private: /*! Vector of all particle species */ - ParticleSpeciesMap m_particle_species; + ParticleSpeciesMap m_particle_species; /*! Vector of particle species names */ ParticlesNamesVector m_namelist; /*! List with names of unallocated species */