diff options
35 files changed, 1036 insertions, 993 deletions
diff --git a/Exec/Langmuir/ParticleProb.cpp b/Exec/Langmuir/ParticleProb.cpp deleted file mode 100644 index b14d97bf8..000000000 --- a/Exec/Langmuir/ParticleProb.cpp +++ /dev/null @@ -1,118 +0,0 @@ - -// -// Each problem must have its own version of PhysicalParticleContainer::InitData() -// to initialize the particle data. It must also initialize charge and mass. -// - -#include <cmath> - -#include <AMReX_BLProfiler.H> - -#include <ParticleContainer.H> -#include <WarpXConst.H> - -using namespace amrex; - -void -PhysicalParticleContainer::InitData() -{ - BL_PROFILE("PPC::InitData()"); - - if (species_id == 0) { - charge = -PhysConst::q_e; - mass = PhysConst::m_e; - } else if (species_id == 1) { - charge = PhysConst::q_e; - mass = PhysConst::m_e; - } else { - amrex::Abort("PhysicalParticleContainer::InitData(): species_id must be 0 or 1"); - } - - const int lev = 0; - - const Geometry& geom = Geom(lev); - const Real* dx = geom.CellSize(); - - Real weight, ux, uy, uz; - Real particle_xmin, particle_xmax, particle_ymin, particle_ymax, particle_zmin, particle_zmax; - int n_part_per_cell; - { - ParmParse pp("langmuirwave"); - n_part_per_cell = 1; - pp.query("num_particles_per_cell", n_part_per_cell); - weight = 1.e25; - pp.query("n_e", weight); - #if BL_SPACEDIM==3 - weight *= dx[0]*dx[1]*dx[2]/n_part_per_cell; - #elif BL_SPACEDIM==2 - weight *= dx[0]*dx[1]/n_part_per_cell; - #endif - - particle_xmin = particle_ymin = particle_zmin = -2.e-5; - particle_xmax = particle_ymax = particle_zmax = 2.e-5; - pp.query("particle_xmin", particle_xmin); - pp.query("particle_xmax", particle_xmax); - pp.query("particle_ymin", particle_ymin); - pp.query("particle_ymax", particle_ymax); - pp.query("particle_zmin", particle_zmin); - pp.query("particle_zmax", particle_zmax); - - ux = 0.; - uy = 0.; - uz = 0.; - if (species_id == 0) { // electrons - pp.query("ux", ux); - pp.query("uy", uy); - pp.query("uz", uz); - } - - ux *= PhysConst::c; - uy *= PhysConst::c; - uz *= PhysConst::c; - } - - std::array<Real,PIdx::nattribs> attribs; - attribs.fill(0.0); - attribs[PIdx::w ] = weight; - attribs[PIdx::ux] = ux; - attribs[PIdx::uy] = uy; - attribs[PIdx::uz] = uz; - - for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) - { - const Box& tile_box = mfi.tilebox(); - RealBox tile_real_box { tile_box, dx, geom.ProbLo() }; - - const int grid_id = mfi.index(); - const int tile_id = mfi.LocalTileIndex(); - - const IntVect& boxlo = tile_box.smallEnd(); - for (IntVect iv = boxlo; iv <= tile_box.bigEnd(); tile_box.next(iv)) - { - for (int i_part=0; i_part<n_part_per_cell;i_part++) - { - Real particle_shift = (0.5+i_part)/n_part_per_cell; -#if (BL_SPACEDIM == 3) - Real x = tile_real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift)*dx[0]; - Real y = tile_real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1]; - Real z = tile_real_box.lo(2) + (iv[2]-boxlo[2] + particle_shift)*dx[2]; -#elif (BL_SPACEDIM == 2) - Real x = tile_real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift)*dx[0]; - Real y = 0.0; - Real z = tile_real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1]; -#endif - - if (x >= particle_xmax || x < particle_xmin || - y >= particle_ymax || y < particle_ymin || - z >= particle_zmax || z < particle_zmin ) - { - continue; - } - else - { - AddOneParticle(lev, grid_id, tile_id, x, y, z, attribs); - } - } - } - } -} diff --git a/Exec/Langmuir/inputs b/Exec/Langmuir/inputs index f4bf84979..292cd19de 100644 --- a/Exec/Langmuir/inputs +++ b/Exec/Langmuir/inputs @@ -35,13 +35,22 @@ algo.particle_pusher = 0 warpx.cfl = 1.0 particles.nspecies = 1 +particles.species_names = electrons -langmuirwave.particle_xmin = -20.e-6 -langmuirwave.particle_xmax = 0.e-6 -langmuirwave.particle_ymin = -20.e-6 -langmuirwave.particle_ymax = 20.e-6 -langmuirwave.particle_zmin = -20.e-6 -langmuirwave.particle_zmax = 20.e-6 -langmuirwave.num_particles_per_cell = 10 -langmuirwave.n_e = 1.e25 # number of electrons per m^3 -langmuirwave.ux = 0.01 # ux = gamma*beta_x +electrons.charge = -q_e +electrons.mass = m_e +electrons.injection_style = "NDiagPerCell" + +electrons.xmin = -20.e-6 +electrons.xmax = 0.e-6 +electrons.ymin = -20.e-6 +electrons.ymax = 20.e-6 +electrons.zmin = -20.e-6 +electrons.zmax = 20.e-6 + +electrons.num_particles_per_cell = 10 +electrons.profile = constant +electrons.density = 1.e25 # number of electrons per m^3 + +electrons.momentum_distribution_type = "constant" +electrons.ux = 0.01 # ux = gamma*beta_x diff --git a/Exec/Langmuir/inputs.lb b/Exec/Langmuir/inputs.lb index 5e33e0dde..5652735b6 100644 --- a/Exec/Langmuir/inputs.lb +++ b/Exec/Langmuir/inputs.lb @@ -41,13 +41,22 @@ algo.particle_pusher = 0 warpx.cfl = 1.0 particles.nspecies = 1 +particles.species_names = electrons -langmuirwave.particle_xmin = -20.e-6 -langmuirwave.particle_xmax = 0.e-6 -langmuirwave.particle_ymin = -20.e-6 -langmuirwave.particle_ymax = 20.e-6 -langmuirwave.particle_zmin = -20.e-6 -langmuirwave.particle_zmax = 20.e-6 -langmuirwave.num_particles_per_cell = 10 -langmuirwave.n_e = 1.e25 # number of electrons per m^3 -langmuirwave.ux = 1.0 # ux = gamma*beta_x +electrons.charge = -q_e +electrons.mass = m_e +electrons.injection_style = "NDiagPerCell" + +electrons.xmin = -20.e-6 +electrons.xmax = 0.e-6 +electrons.ymin = -20.e-6 +electrons.ymax = 20.e-6 +electrons.zmin = -20.e-6 +electrons.zmax = 20.e-6 + +electrons.num_particles_per_cell = 10 +electrons.profile = constant +electrons.density = 1.e25 # number of electrons per m^3 + +electrons.momentum_distribution_type = "constant" +electrons.ux = 1.0 # ux = gamma*beta_x diff --git a/Exec/Langmuir/inputs.multi.rt b/Exec/Langmuir/inputs.multi.rt index 6fecda4fb..0f0a8ff55 100644 --- a/Exec/Langmuir/inputs.multi.rt +++ b/Exec/Langmuir/inputs.multi.rt @@ -38,15 +38,44 @@ warpx.cfl = 1.0 # Particles particles.nspecies = 2 +particles.species_names = electrons positrons -langmuirwave.particle_xmin = -20.e-6 -langmuirwave.particle_xmax = 0.e-6 -langmuirwave.particle_ymin = -20.e-6 -langmuirwave.particle_ymax = 20.e-6 -langmuirwave.particle_zmin = -20.e-6 -langmuirwave.particle_zmax = 20.e-6 -langmuirwave.num_particles_per_cell = 10 -langmuirwave.n_e = 1.e25 # number of electrons per m^3 -langmuirwave.ux = 0.01 # ux = gamma*beta_x -langmuirwave.uy = 0. # uy = gamma*beta_y -langmuirwave.uz = 0. # uz = gamma*beta_z +electrons.charge = -q_e +electrons.mass = m_e +electrons.injection_style = "NDiagPerCell" + +electrons.xmin = -20.e-6 +electrons.xmax = 0.e-6 +electrons.ymin = -20.e-6 +electrons.ymax = 20.e-6 +electrons.zmin = -20.e-6 +electrons.zmax = 20.e-6 + +electrons.num_particles_per_cell = 10 +electrons.profile = constant +electrons.density = 1.e25 # number of electrons per m^3 + +electrons.momentum_distribution_type = "constant" +electrons.ux = 0.01 # ux = gamma*beta_x +electrons.uy = 0. # uy = gamma*beta_y +electrons.uz = 0. # uz = gamma*beta_z + +positrons.charge = q_e +positrons.mass = m_e +positrons.injection_style = "NDiagPerCell" + +positrons.xmin = -20.e-6 +positrons.xmax = 0.e-6 +positrons.ymin = -20.e-6 +positrons.ymax = 20.e-6 +positrons.zmin = -20.e-6 +positrons.zmax = 20.e-6 + +positrons.num_particles_per_cell = 10 +positrons.profile = constant +positrons.density = 1.e25 # number of positrons per m^3 + +positrons.momentum_distribution_type = "constant" +positrons.ux = 0. # ux = gamma*beta_x +positrons.uy = 0. # uy = gamma*beta_y +positrons.uz = 0. # uz = gamma*beta_z diff --git a/Exec/Langmuir/inputs.nolb b/Exec/Langmuir/inputs.nolb index 1851bcf38..6ca55b3c7 100644 --- a/Exec/Langmuir/inputs.nolb +++ b/Exec/Langmuir/inputs.nolb @@ -41,13 +41,22 @@ algo.particle_pusher = 0 warpx.cfl = 1.0 particles.nspecies = 1 +particles.species_names = electrons -langmuirwave.particle_xmin = -20.e-6 -langmuirwave.particle_xmax = 0.e-6 -langmuirwave.particle_ymin = -20.e-6 -langmuirwave.particle_ymax = 20.e-6 -langmuirwave.particle_zmin = -20.e-6 -langmuirwave.particle_zmax = 20.e-6 -langmuirwave.num_particles_per_cell = 10 -langmuirwave.n_e = 1.e25 # number of electrons per m^3 -langmuirwave.ux = 1.0 # ux = gamma*beta_x +electrons.charge = -q_e +electrons.mass = m_e +electrons.injection_style = "NDiagPerCell" + +electrons.xmin = -20.e-6 +electrons.xmax = 0.e-6 +electrons.ymin = -20.e-6 +electrons.ymax = 20.e-6 +electrons.zmin = -20.e-6 +electrons.zmax = 20.e-6 + +electrons.num_particles_per_cell = 10 +electrons.profile = constant +electrons.density = 1.e25 # number of electrons per m^3 + +electrons.momentum_distribution_type = "constant" +electrons.ux = 1.0 # ux = gamma*beta_x diff --git a/Exec/Langmuir/inputs.rt b/Exec/Langmuir/inputs.rt index 2999aa4b9..3ade966f1 100644 --- a/Exec/Langmuir/inputs.rt +++ b/Exec/Langmuir/inputs.rt @@ -36,14 +36,23 @@ interpolation.noz = 1 # CFL warpx.cfl = 1.0 -langmuirwave.particle_xmin = -20.e-6 -langmuirwave.particle_xmax = 20.e-6 -langmuirwave.particle_ymin = -20.e-6 -langmuirwave.particle_ymax = 20.e-6 -langmuirwave.particle_zmin = -20.e-6 -langmuirwave.particle_zmax = 20.e-6 -langmuirwave.num_particles_per_cell = 10 -langmuirwave.n_e = 1.e25 # number of electrons per m^3 -langmuirwave.ux = 0. # ux = gamma*beta_x -langmuirwave.uy = 0. # uy = gamma*beta_y -langmuirwave.uz = 0. # uz = gamma*beta_z +particles.nspecies = 1 +particles.species_names = electrons + +electrons.charge = -q_e +electrons.mass = m_e +electrons.injection_style = "NDiagPerCell" + +electrons.xmin = -20.e-6 +electrons.xmax = 20.e-6 +electrons.ymin = -20.e-6 +electrons.ymax = 20.e-6 +electrons.zmin = -20.e-6 +electrons.zmax = 20.e-6 + +electrons.num_particles_per_cell = 10 +electrons.profile = constant +electrons.density = 1.e25 # number of electrons per m^3 + +electrons.momentum_distribution_type = "constant" + diff --git a/Exec/charged_beam/ParticleProb.cpp b/Exec/charged_beam/ParticleProb.cpp deleted file mode 100644 index 398ddf55b..000000000 --- a/Exec/charged_beam/ParticleProb.cpp +++ /dev/null @@ -1,110 +0,0 @@ - -// -// Each problem must have its own version of PhysicalParticleContainer::InitData() -// to initialize the particle data. It must also initialize charge and mass. -// - -#include <cmath> - -#include <AMReX_BLProfiler.H> - -#include <ParticleContainer.H> -#include <WarpXConst.H> - -using namespace amrex; - -void -PhysicalParticleContainer::InitData() -{ - BL_PROFILE("PPC::InitData()"); - - charge = -PhysConst::q_e; - mass = PhysConst::m_e; - - const int lev = 0; - - const Geometry& geom = Geom(lev); - const Real* dx = geom.CellSize(); - - Real weight, ux, uy, uz; - Real particle_xmin, particle_xmax, particle_ymin, particle_ymax, particle_zmin, particle_zmax; - int n_part_per_cell; - { - ParmParse pp("chargedbeam"); - n_part_per_cell = 1; - pp.query("num_particles_per_cell", n_part_per_cell); - weight = 1.e25; - pp.query("n_e", weight); - #if BL_SPACEDIM==3 - weight *= dx[0]*dx[1]*dx[2]/n_part_per_cell; - #elif BL_SPACEDIM==2 - weight *= dx[0]*dx[1]/n_part_per_cell; - #endif - - particle_xmin = particle_ymin = particle_zmin = -2.e-5; - particle_xmax = particle_ymax = particle_zmax = 2.e-5; - pp.query("particle_xmin", particle_xmin); - pp.query("particle_xmax", particle_xmax); - pp.query("particle_ymin", particle_ymin); - pp.query("particle_ymax", particle_ymax); - pp.query("particle_zmin", particle_zmin); - pp.query("particle_zmax", particle_zmax); - - ux = 0.; - uy = 0.; - uz = 0.; - pp.query("ux", ux); - pp.query("uy", uy); - pp.query("uz", uz); - - Real gamma = 1./std::sqrt(1.0 - ux*ux - uy*uy - uz*uz); - ux *= PhysConst::c*gamma; - uy *= PhysConst::c*gamma; - uz *= PhysConst::c*gamma; - } - - std::array<Real,PIdx::nattribs> attribs; - attribs.fill(0.0); - attribs[PIdx::w ] = weight; - attribs[PIdx::ux] = ux; - attribs[PIdx::uy] = uy; - attribs[PIdx::uz] = uz; - - for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) - { - const Box& tile_box = mfi.tilebox(); - RealBox tile_real_box { tile_box, dx, geom.ProbLo() }; - - const int grid_id = mfi.index(); - const int tile_id = mfi.LocalTileIndex(); - - const auto& boxlo = tile_box.smallEnd(); - for (IntVect iv = tile_box.smallEnd(); iv <= tile_box.bigEnd(); tile_box.next(iv)) - { - for (int i_part=0; i_part<n_part_per_cell;i_part++) - { - Real particle_shift = (0.5+i_part)/n_part_per_cell; -#if (BL_SPACEDIM == 3) - Real x = tile_real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift)*dx[0]; - Real y = tile_real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1]; - Real z = tile_real_box.lo(2) + (iv[2]-boxlo[2] + particle_shift)*dx[2]; -#elif (BL_SPACEDIM == 2) - Real x = tile_real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift)*dx[0]; - Real y = 0.0; - Real z = tile_real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1]; -#endif - - if (x >= particle_xmax || x < particle_xmin || - y >= particle_ymax || y < particle_ymin || - z >= particle_zmax || z < particle_zmin ) - { - continue; - } - else - { - AddOneParticle(lev, grid_id, tile_id, x, y, z, attribs); - } - } - } - } -} diff --git a/Exec/charged_beam/inputs b/Exec/charged_beam/inputs index a02198517..48b796d97 100644 --- a/Exec/charged_beam/inputs +++ b/Exec/charged_beam/inputs @@ -31,12 +31,23 @@ algo.particle_pusher = 0 # CFL warpx.cfl = 1.0 -chargedbeam.particle_xmin = -20.e-6 -chargedbeam.particle_xmax = 0.e-6 -chargedbeam.particle_ymin = -20.e-6 -chargedbeam.particle_ymax = 20.e-6 -chargedbeam.particle_zmin = -20.e-6 -chargedbeam.particle_zmax = 20.e-6 -chargedbeam.num_particles_per_cell = 10 -chargedbeam.n_e = 1.e25 # number of electrons per m^3 -chargedbeam.ux = 0.01 # ux = gamma*beta_x +particles.nspecies = 1 +particles.species_names = electrons + +electrons.charge = -q_e +electrons.mass = m_e +electrons.injection_style = "NDiagPerCell" + +electrons.xmin = -20.e-6 +electrons.xmax = 0.e-6 +electrons.ymin = -20.e-6 +electrons.ymax = 20.e-6 +electrons.zmin = -20.e-6 +electrons.zmax = 20.e-6 + +electrons.num_particles_per_cell = 10 +electrons.profile = constant +electrons.density = 1.e25 # number of electrons per m^3 + +electrons.momentum_distribution_type = "constant" +electrons.ux = 0.01 # ux = gamma*beta_x diff --git a/Exec/laser_acceleration/ParticleProb.cpp b/Exec/laser_acceleration/ParticleProb.cpp deleted file mode 100644 index 031e1e9cd..000000000 --- a/Exec/laser_acceleration/ParticleProb.cpp +++ /dev/null @@ -1,111 +0,0 @@ - -// -// Each problem must have its own version of PhysicalParticleContainer::InitData() -// to initialize the particle data. It must also initialize charge and mass. -// - -#include <cmath> - -#include <AMReX_BLProfiler.H> - -#include <ParticleContainer.H> -#include <WarpXConst.H> - -using namespace amrex; - -void -PhysicalParticleContainer::InitData() -{ - BL_PROFILE("SPC::InitData()"); - - // species_id 0 : the electrons of the plasma - // Note: the ions of the plasma are implicitly motionless, and so are not part of the simulation - if (species_id == 0) { - charge = -PhysConst::q_e; - mass = PhysConst::m_e; - } else { - amrex::Abort("PhysicalParticleContainer::InitData(): species_id must be 0 or 1"); - } - - const int lev = 0; - - const Geometry& geom = Geom(lev); - const Real* dx = geom.CellSize(); - - Real weight, gamma, uz; - Real particle_xmin, particle_xmax, particle_ymin, particle_ymax, particle_zmin, particle_zmax; - int n_part_per_cell; - { - ParmParse pp("lwfa"); - - // Calculate the particle weight - n_part_per_cell = 1; - pp.query("num_particles_per_cell", n_part_per_cell); - weight = 1.e22; - if (species_id == 0) { - pp.query("plasma_density", weight); - } - #if BL_SPACEDIM==3 - weight *= dx[0]*dx[1]*dx[2]/n_part_per_cell; - #elif BL_SPACEDIM==2 - weight *= dx[0]*dx[1]/n_part_per_cell; - #endif - - // Calculate the limits between which the particles will be initialized - particle_xmin = particle_ymin = particle_zmin = -2.e-5; - particle_xmax = particle_ymax = particle_zmax = 2.e-5; - if (species_id == 0) { - pp.query("plasma_xmin", particle_xmin); - pp.query("plasma_xmax", particle_xmax); - pp.query("plasma_ymin", particle_ymin); - pp.query("plasma_ymax", particle_ymax); - pp.query("plasma_zmin", particle_zmin); - pp.query("plasma_zmax", particle_zmax); - } - uz = 0.; - uz *= PhysConst::c; - } - - std::array<Real,PIdx::nattribs> attribs; - attribs.fill(0.0); - attribs[PIdx::w ] = weight; - attribs[PIdx::uz] = uz; - - for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) - { - const Box& tile_box = mfi.tilebox(); - RealBox tile_real_box { tile_box, dx, geom.ProbLo() }; - - const int grid_id = mfi.index(); - const int tile_id = mfi.LocalTileIndex(); - - const auto& boxlo = tile_box.smallEnd(); - for (IntVect iv = tile_box.smallEnd(); iv <= tile_box.bigEnd(); tile_box.next(iv)) - { - for (int i_part=0; i_part<n_part_per_cell;i_part++) - { - Real particle_shift = (0.5+i_part)/n_part_per_cell; -#if (BL_SPACEDIM == 3) - Real x = tile_real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift)*dx[0]; - Real y = tile_real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1]; - Real z = tile_real_box.lo(2) + (iv[2]-boxlo[2] + particle_shift)*dx[2]; -#elif (BL_SPACEDIM == 2) - Real x = tile_real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift)*dx[0]; - Real y = 0.0; - Real z = tile_real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1]; -#endif - - if (x >= particle_xmax || x < particle_xmin || - y >= particle_ymax || y < particle_ymin || - z >= particle_zmax || z < particle_zmin ) - { - continue; - } - else - { - AddOneParticle(lev, grid_id, tile_id, x, y, z, attribs); - } - } - } - } -} diff --git a/Exec/laser_acceleration/inputs b/Exec/laser_acceleration/inputs index 098b24582..09b60801a 100644 --- a/Exec/laser_acceleration/inputs +++ b/Exec/laser_acceleration/inputs @@ -31,17 +31,23 @@ algo.particle_pusher = 0 # CFL warpx.cfl = 1.0 particles.nspecies = 1 +particles.species_names = electrons -lwfa.num_particles_per_cell = 4 +electrons.charge = -q_e +electrons.mass = m_e +electrons.injection_style = "NDiagPerCell" -# The plasma -lwfa.plasma_xmin = -20.e-6 -lwfa.plasma_xmax = 20.e-6 -lwfa.plasma_ymin = -20.e-6 -lwfa.plasma_ymax = 20.e-6 -lwfa.plasma_zmin = 0.e-6 -lwfa.plasma_zmax = 20.e-6 -lwfa.plasma_density = 1e23 # number of electrons per m^3 +electrons.xmin = -20.e-6 +electrons.xmax = 20.e-6 +electrons.ymin = -20.e-6 +electrons.ymax = 20.e-6 +electrons.zmin = 0.e-6 + +electrons.num_particles_per_cell = 4 +electrons.profile = constant +electrons.density = 1.e23 # number of electrons per m^3 + +electrons.momentum_distribution_type = "constant" # Moving window warpx.do_moving_window = 1 diff --git a/Exec/laser_injection/ParticleProb.cpp b/Exec/laser_injection/ParticleProb.cpp deleted file mode 100644 index 35376c4a2..000000000 --- a/Exec/laser_injection/ParticleProb.cpp +++ /dev/null @@ -1,99 +0,0 @@ - -// -// Each problem must have its own version of PhysicalParticleContainer::InitData() -// to initialize the particle data. It must also initialize charge and mass. -// - -#include <cmath> - -#include <AMReX_BLProfiler.H> - -#include <ParticleContainer.H> -#include <WarpXConst.H> -#include <random> - -using namespace amrex; - -void -PhysicalParticleContainer::InitData() -{ - BL_PROFILE("PPC::InitData()"); - - charge = -PhysConst::q_e; - mass = PhysConst::m_e; - - const int lev = 0; - - const Geometry& geom = Geom(lev); - const Real* dx = geom.CellSize(); - - Real weight, u_th; - int n_part_per_cell; - { - ParmParse pp("uniform_plasma"); - n_part_per_cell = 1; - pp.query("num_particles_per_cell", n_part_per_cell); - weight = 1.e25; - pp.query("n_e", weight); - #if BL_SPACEDIM==3 - weight *= dx[0]*dx[1]*dx[2]/n_part_per_cell; - #elif BL_SPACEDIM==2 - weight *= dx[0]*dx[1]/n_part_per_cell; - #endif - - u_th = 0.; - pp.query("u_th", u_th); - u_th *= PhysConst::c; - } - - std::array<Real,PIdx::nattribs> attribs; - attribs.fill(0.0); - attribs[PIdx::w ] = weight; - - // Initialize random generator for normal distribution - std::default_random_engine generator; - std::normal_distribution<Real> velocity_distribution(0.0, u_th); - std::uniform_real_distribution<Real> position_distribution(0.0,1.0); - - for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) - { - const Box& tile_box = mfi.tilebox(); - RealBox tile_real_box { tile_box, dx, geom.ProbLo() }; - - const int grid_id = mfi.index(); - const int tile_id = mfi.LocalTileIndex(); - - const auto& boxlo = tile_box.smallEnd(); - for (IntVect iv = tile_box.smallEnd(); iv <= tile_box.bigEnd(); tile_box.next(iv)) - { - for (int i_part=0; i_part<n_part_per_cell;i_part++) - { - // Randomly generate the speed according to a normal distribution - Real ux = velocity_distribution(generator); - Real uy = velocity_distribution(generator); - Real uz = velocity_distribution(generator); - - // Randomly generate the positions (uniformly inside each cell) - Real particle_shift_x = position_distribution(generator); - Real particle_shift_y = position_distribution(generator); - Real particle_shift_z = position_distribution(generator); - -#if (BL_SPACEDIM == 3) - Real x = tile_real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift_x)*dx[0]; - Real y = tile_real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift_y)*dx[1]; - Real z = tile_real_box.lo(2) + (iv[2]-boxlo[2] + particle_shift_z)*dx[2]; -#elif (BL_SPACEDIM == 2) - Real x = tile_real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift_x)*dx[0]; - Real y = 0.0; - Real z = tile_real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift_z)*dx[1]; -#endif - - attribs[PIdx::ux] = ux; - attribs[PIdx::uy] = uy; - attribs[PIdx::uz] = uz; - - AddOneParticle(lev, grid_id, tile_id, x, y, z, attribs); - } - } - } -} diff --git a/Exec/plasma_acceleration/CustomDensityProb.cpp b/Exec/plasma_acceleration/CustomDensityProb.cpp new file mode 100644 index 000000000..3efcb13c5 --- /dev/null +++ b/Exec/plasma_acceleration/CustomDensityProb.cpp @@ -0,0 +1,12 @@ +#include <PlasmaInjector.H> + +#include <iostream> + +using namespace amrex; + +/// +/// This "custom" density profile just does constant +/// +Real CustomDensityProfile::getDensity(Real x, Real y, Real z) const { + return params[0]; +} diff --git a/Exec/plasma_acceleration/ParticleProb.cpp b/Exec/plasma_acceleration/ParticleProb.cpp deleted file mode 100644 index 7a538abdb..000000000 --- a/Exec/plasma_acceleration/ParticleProb.cpp +++ /dev/null @@ -1,126 +0,0 @@ - -// -// Each problem must have its own version of PhysicalParticleContainer::InitData() -// to initialize the particle data. It must also initialize charge and mass. -// - -#include <cmath> - -#include <AMReX_BLProfiler.H> - -#include <ParticleContainer.H> -#include <WarpXConst.H> - -using namespace amrex; - -void -PhysicalParticleContainer::InitData() -{ - BL_PROFILE("SPC::InitData()"); - - // species_id 0 : the beam - // species_id 1 : the electrons of the plasma - // Note: the ions of the plasma are implicitly motionless, and so are not part of the simulation - if (species_id == 0 or species_id == 1) { - charge = -PhysConst::q_e; - mass = PhysConst::m_e; - } else { - amrex::Abort("PhysicalParticleContainer::InitData(): species_id must be 0 or 1"); - } - - const int lev = 0; - - const Geometry& geom = Geom(lev); - const Real* dx = geom.CellSize(); - - Real weight, gamma, uz; - Real particle_xmin, particle_xmax, particle_ymin, particle_ymax, particle_zmin, particle_zmax; - int n_part_per_cell; - { - ParmParse pp("pwfa"); - - // Calculate the particle weight - n_part_per_cell = 1; - pp.query("num_particles_per_cell", n_part_per_cell); - weight = 1.e22; - if (species_id == 0) { - pp.query("beam_density", weight); - } else if (species_id == 1) { - pp.query("plasma_density", weight); - } - #if BL_SPACEDIM==3 - weight *= dx[0]*dx[1]*dx[2]/n_part_per_cell; - #elif BL_SPACEDIM==2 - weight *= dx[0]*dx[1]/n_part_per_cell; - #endif - - // Calculate the limits between which the particles will be initialized - particle_xmin = particle_ymin = particle_zmin = -2.e-5; - particle_xmax = particle_ymax = particle_zmax = 2.e-5; - if (species_id == 0) { - pp.query("beam_xmin", particle_xmin); - pp.query("beam_xmax", particle_xmax); - pp.query("beam_ymin", particle_ymin); - pp.query("beam_ymax", particle_ymax); - pp.query("beam_zmin", particle_zmin); - pp.query("beam_zmax", particle_zmax); - } else if (species_id == 1) { - pp.query("plasma_xmin", particle_xmin); - pp.query("plasma_xmax", particle_xmax); - pp.query("plasma_ymin", particle_ymin); - pp.query("plasma_ymax", particle_ymax); - pp.query("plasma_zmin", particle_zmin); - pp.query("plasma_zmax", particle_zmax); - } - uz = 0.; - if (species_id == 0) { // relativistic beam - gamma = 1.; - pp.query("beam_gamma", gamma); - uz = std::sqrt( gamma*gamma - 1 ); - } - uz *= PhysConst::c; - } - - std::array<Real,PIdx::nattribs> attribs; - attribs.fill(0.0); - attribs[PIdx::w ] = weight; - attribs[PIdx::uz] = uz; - - for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) - { - const Box& tile_box = mfi.tilebox(); - RealBox tile_real_box { tile_box, dx, geom.ProbLo() }; - - const int grid_id = mfi.index(); - const int tile_id = mfi.LocalTileIndex(); - - const auto& boxlo = tile_box.smallEnd(); - for (IntVect iv = tile_box.smallEnd(); iv <= tile_box.bigEnd(); tile_box.next(iv)) - { - for (int i_part=0; i_part<n_part_per_cell;i_part++) - { - Real particle_shift = (0.5+i_part)/n_part_per_cell; -#if (BL_SPACEDIM == 3) - Real x = tile_real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift)*dx[0]; - Real y = tile_real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1]; - Real z = tile_real_box.lo(2) + (iv[2]-boxlo[2] + particle_shift)*dx[2]; -#elif (BL_SPACEDIM == 2) - Real x = tile_real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift)*dx[0]; - Real y = 0.0; - Real z = tile_real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1]; -#endif - - if (x >= particle_xmax || x < particle_xmin || - y >= particle_ymax || y < particle_ymin || - z >= particle_zmax || z < particle_zmin ) - { - continue; - } - else - { - AddOneParticle(lev, grid_id, tile_id, x, y, z, attribs); - } - } - } - } -} diff --git a/Exec/plasma_acceleration/inputs b/Exec/plasma_acceleration/inputs index 26ade8abf..2aa9ac536 100644 --- a/Exec/plasma_acceleration/inputs +++ b/Exec/plasma_acceleration/inputs @@ -30,27 +30,57 @@ algo.particle_pusher = 0 # CFL warpx.cfl = 1.0 + +# Information about the particle species particles.nspecies = 2 +particles.species_names = beam plasma + +# +# The beam species information +# + +beam.charge = -q_e +beam.mass = m_e +beam.injection_style = "NDiagPerCell" + +beam.xmin = -20.e-6 +beam.xmax = 20.e-6 +beam.ymin = -20.e-6 +beam.ymax = 20.e-6 +beam.zmin = -150.e-6 +beam.zmax = -100.e-6 + +beam.profile = "constant" +beam.density = 1e23 # number of particles per m^3 +beam.num_particles_per_cell = 4 + +beam.momentum_distribution_type = "constant" +beam.ux = 0.0 +beam.uy = 0.0 +beam.uz = 1e9 + +# +# The plasma species information +# + +plasma.charge = -q_e +plasma.mass = m_e +plasma.injection_style = "NDiagPerCell" + +plasma.xmin = -200.e-6 +plasma.xmax = 200.e-6 +plasma.ymin = -200.e-6 +plasma.ymax = 200.e-6 +plasma.zmin = 0.e-6 + +plasma.profile = "custom" +plasma.custom_profile_params = 1e22 # number of particles per m^3 +plasma.num_particles_per_cell = 4 -pwfa.num_particles_per_cell = 4 - -# The beam -pwfa.beam_xmin = -20.e-6 -pwfa.beam_xmax = 20.e-6 -pwfa.beam_ymin = -20.e-6 -pwfa.beam_ymax = 20.e-6 -pwfa.beam_zmin = -150.e-6 -pwfa.beam_zmax = -100.e-6 -pwfa.beam_density = 1e23 # number of electrons per m^3 -pwfa.beam_gamma = 1e9 -# The plasma -pwfa.plasma_xmin = -200.e-6 -pwfa.plasma_xmax = 200.e-6 -pwfa.plasma_ymin = -200.e-6 -pwfa.plasma_ymax = 200.e-6 -pwfa.plasma_zmin = 0.e-6 -pwfa.plasma_zmax = 200.e-6 -pwfa.plasma_density = 1e22 # number of electrons per m^3 +plasma.momentum_distribution_type = "constant" +plasma.ux = 0.0 +plasma.uy = 0.0 +plasma.uz = 0.0 # Moving window warpx.do_moving_window = 1 @@ -61,5 +91,3 @@ warpx.moving_window_v = 1.0 # in units of the speed of light warpx.do_plasma_injection = 1 warpx.num_injected_species = 1 warpx.injected_plasma_species = 1 0 -warpx.injected_plasma_density = 1e22 1e22 -warpx.injected_plasma_ppc = 4 4 diff --git a/Exec/test_particle_injection/ParticleProb.cpp b/Exec/test_particle_injection/ParticleProb.cpp deleted file mode 100644 index 7cdee1d46..000000000 --- a/Exec/test_particle_injection/ParticleProb.cpp +++ /dev/null @@ -1,109 +0,0 @@ - -// -// Each problem must have its own version of PhysicalParticleContainer::InitData() -// to initialize the particle data. It must also initialize charge and mass. -// - -#include <cmath> - -#include <AMReX_BLProfiler.H> - -#include <ParticleContainer.H> -#include <WarpXConst.H> - -using namespace amrex; - -void -PhysicalParticleContainer::InitData() -{ - BL_PROFILE("PPC::InitData()"); - - charge = -PhysConst::q_e; - mass = PhysConst::m_e; - - const int lev = 0; - - const Geometry& geom = Geom(lev); - const Real* dx = geom.CellSize(); - - Real weight, ux, uy, uz; - Real particle_xmin, particle_xmax, particle_ymin, particle_ymax, particle_zmin, particle_zmax; - int n_part_per_cell; - { - ParmParse pp("langmuirwave"); - n_part_per_cell = 1; - pp.query("num_particles_per_cell", n_part_per_cell); - weight = 1.e25; - pp.query("n_e", weight); - #if BL_SPACEDIM==3 - weight *= dx[0]*dx[1]*dx[2]/n_part_per_cell; - #elif BL_SPACEDIM==2 - weight *= dx[0]*dx[1]/n_part_per_cell; - #endif - - particle_xmin = particle_ymin = particle_zmin = -2.e-5; - particle_xmax = particle_ymax = particle_zmax = 2.e-5; - pp.query("particle_xmin", particle_xmin); - pp.query("particle_xmax", particle_xmax); - pp.query("particle_ymin", particle_ymin); - pp.query("particle_ymax", particle_ymax); - pp.query("particle_zmin", particle_zmin); - pp.query("particle_zmax", particle_zmax); - - ux = 0.; - uy = 0.; - uz = 0.; - pp.query("ux", ux); - pp.query("uy", uy); - pp.query("uz", uz); - - ux *= PhysConst::c; - uy *= PhysConst::c; - uz *= PhysConst::c; - } - - std::array<Real,PIdx::nattribs> attribs; - attribs.fill(0.0); - attribs[PIdx::w ] = weight; - attribs[PIdx::ux] = ux; - attribs[PIdx::uy] = uy; - attribs[PIdx::uz] = uz; - - for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) - { - const Box& tile_box = mfi.tilebox(); - RealBox tile_real_box { tile_box, dx, geom.ProbLo() }; - - const int grid_id = mfi.index(); - const int tile_id = mfi.LocalTileIndex(); - - const IntVect& boxlo = tile_box.smallEnd(); - for (IntVect iv = boxlo; iv <= tile_box.bigEnd(); tile_box.next(iv)) - { - for (int i_part=0; i_part<n_part_per_cell;i_part++) - { - Real particle_shift = (0.5+i_part)/n_part_per_cell; -#if (BL_SPACEDIM == 3) - Real x = tile_real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift)*dx[0]; - Real y = tile_real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1]; - Real z = tile_real_box.lo(2) + (iv[2]-boxlo[2] + particle_shift)*dx[2]; -#elif (BL_SPACEDIM == 2) - Real x = tile_real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift)*dx[0]; - Real y = 0.0; - Real z = tile_real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1]; -#endif - - if (x >= particle_xmax || x < particle_xmin || - y >= particle_ymax || y < particle_ymin || - z >= particle_zmax || z < particle_zmin ) - { - continue; - } - else - { - AddOneParticle(lev, grid_id, tile_id, x, y, z, attribs); - } - } - } - } -} diff --git a/Exec/test_particle_injection/inputs b/Exec/test_particle_injection/inputs index 26e88f6c6..6681cbb3c 100644 --- a/Exec/test_particle_injection/inputs +++ b/Exec/test_particle_injection/inputs @@ -31,12 +31,23 @@ algo.particle_pusher = 0 # CFL warpx.cfl = 1.0 -langmuirwave.particle_xmin = -20.e-6 -langmuirwave.particle_xmax = 0.e-6 -langmuirwave.particle_ymin = -20.e-6 -langmuirwave.particle_ymax = 20.e-6 -langmuirwave.particle_zmin = -20.e-6 -langmuirwave.particle_zmax = 20.e-6 -langmuirwave.num_particles_per_cell = 10 -langmuirwave.n_e = 1.e25 # number of electrons per m^3 -langmuirwave.ux = 0.01 # ux = gamma*beta_x +particles.nspecies = 1 +particles.species_names = electrons + +electrons.charge = -q_e +electrons.mass = m_e +electrons.injection_style = "NDiagPerCell" + +electrons.xmin = -20.e-6 +electrons.xmax = 0.e-6 +electrons.ymin = -20.e-6 +electrons.ymax = 20.e-6 +electrons.zmin = -20.e-6 +electrons.zmax = 20.e-6 + +electrons.num_particles_per_cell = 10 +electrons.profile = constant +electrons.density = 1.e25 # number of electrons per m^3 + +electrons.momentum_distribution_type = "constant" +electrons.ux = 0.01 # ux = gamma*beta_x diff --git a/Exec/uniform_plasma/ParticleProb.cpp b/Exec/uniform_plasma/ParticleProb.cpp deleted file mode 100644 index cef618f61..000000000 --- a/Exec/uniform_plasma/ParticleProb.cpp +++ /dev/null @@ -1,114 +0,0 @@ - -// -// Each problem must have its own version of PhysicalParticleContainer::InitData() -// to initialize the particle data. It must also initialize charge and mass. -// - -#include <cmath> - -#include <AMReX_BLProfiler.H> - -#include <ParticleContainer.H> -#include <WarpXConst.H> -#include <random> - -using namespace amrex; - -void -PhysicalParticleContainer::InitData() -{ - BL_PROFILE("PPC::InitData()"); - - // species_id 0 : electrons - // species_id 1 : Hydrogen ions - // Note: the ions of the plasma are implicitly motionless, and so are not part of the simulation - if (species_id == 0) { - charge = -PhysConst::q_e; - mass = PhysConst::m_e; - } else { - charge = PhysConst::q_e; - mass = PhysConst::m_p; - } - - const int lev = 0; - - const Geometry& geom = Geom(lev); - const Real* dx = geom.CellSize(); - - Real weight, u_th, ux_m, uy_m, uz_m; - Real particle_shift_x, particle_shift_y, particle_shift_z; - int n_part_per_cell; - { - ParmParse pp("uniform_plasma"); - n_part_per_cell = 1; - pp.query("num_particles_per_cell", n_part_per_cell); - weight = 1.e25; - pp.query("n_e", weight); - #if BL_SPACEDIM==3 - weight *= dx[0]*dx[1]*dx[2]/n_part_per_cell; - #elif BL_SPACEDIM==2 - weight *= dx[0]*dx[1]/n_part_per_cell; - #endif - - u_th = 0.; - ux_m = 0.; - uy_m = 0.; - uz_m = 0.; - pp.query("u_th", u_th); - pp.query("ux_m", ux_m); - pp.query("uy_m", uy_m); - pp.query("uz_m", uz_m); - u_th *= PhysConst::c; - ux_m *= PhysConst::c; - uy_m *= PhysConst::c; - uz_m *= PhysConst::c; - } - - std::array<Real,PIdx::nattribs> attribs; - attribs.fill(0.0); - attribs[PIdx::w ] = weight; - - // Initialize random generator for normal distribution - std::default_random_engine generator; - std::normal_distribution<Real> velocity_distribution(0.0, u_th); - - for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) - { - const Box& tile_box = mfi.tilebox(); - RealBox tile_real_box { tile_box, dx, geom.ProbLo() }; - - const int grid_id = mfi.index(); - const int tile_id = mfi.LocalTileIndex(); - - const auto& boxlo = tile_box.smallEnd(); - for (IntVect iv = tile_box.smallEnd(); iv <= tile_box.bigEnd(); tile_box.next(iv)) - { - for (int i_part=0; i_part<n_part_per_cell;i_part++) - { - // Randomly generate the speed according to a normal distribution - Real ux_th = velocity_distribution(generator); - Real uy_th = velocity_distribution(generator); - Real uz_th = velocity_distribution(generator); - - // Randomly generate the positions (uniformly inside each cell) - Real particle_shift = (0.5+i_part)/n_part_per_cell; - -#if (BL_SPACEDIM == 3) - Real x = tile_real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift)*dx[0]; - Real y = tile_real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1]; - Real z = tile_real_box.lo(2) + (iv[2]-boxlo[2] + particle_shift)*dx[2]; -#elif (BL_SPACEDIM == 2) - Real x = tile_real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift)*dx[0]; - Real y = 0.0; - Real z = tile_real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1]; -#endif - - attribs[PIdx::ux] = ux_m + ux_th; - attribs[PIdx::uy] = uy_m + uy_th; - attribs[PIdx::uz] = uz_m + uz_th; - - AddOneParticle(lev, grid_id, tile_id, x, y, z, attribs); - } - } - } -} diff --git a/Exec/uniform_plasma/inputs b/Exec/uniform_plasma/inputs index 5e87c689e..419e8909a 100644 --- a/Exec/uniform_plasma/inputs +++ b/Exec/uniform_plasma/inputs @@ -32,7 +32,26 @@ algo.particle_pusher = 0 # CFL warpx.cfl = 1.0 -uniform_plasma.num_particles_per_cell = 40 -uniform_plasma.n_e = 1.e25 # number of electrons per m^3 -uniform_plasma.u_th = 0.01 # uth the std of the (unitless) momentum -uniform_plasma.uz_m = 10. # Mean momentum along x (unitless) +particles.nspecies = 2 +particles.species_names = electrons H_ions + +electrons.charge = -q_e +electrons.mass = m_e +electrons.injection_style = "NDiagPerCell" +electrons.num_particles_per_cell = 40 +electrons.profile = constant +electrons.density = 1.e25 # number of electrons per m^3 +electrons.momentum_distribution_type = "gaussian" +electrons.u_th = 0.01 # uth the std of the (unitless) momentum +electrons.uz_m = 10. # Mean momentum along z (unitless) + +H_ions.charge = q_e +H_ions.mass = m_p +H_ions.injection_style = "NDiagPerCell" +H_ions.num_particles_per_cell = 40 +H_ions.profile = constant +H_ions.density = 1.e25 # number of H_ions per m^3 +H_ions.momentum_distribution_type = "gaussian" +H_ions.u_th = 0.01 # uth the std of the (unitless) momentum +H_ions.uz_m = 10. # Mean momentum along z (unitless) + diff --git a/Exec/uniform_plasma/inputs.rt b/Exec/uniform_plasma/inputs.rt index 752cab814..bf6bb127a 100644 --- a/Exec/uniform_plasma/inputs.rt +++ b/Exec/uniform_plasma/inputs.rt @@ -31,6 +31,15 @@ warpx.cfl = 1.0 amr.plot_int = 10 -uniform_plasma.num_particles_per_cell = 10 -uniform_plasma.n_e = 1.e25 # number of electrons per m^3 -uniform_plasma.u_th = 0.01 # uth the std of the (unitless) momentum +particles.nspecies = 1 +particles.species_names = electrons + +electrons.charge = -q_e +electrons.mass = m_e +electrons.injection_style = "NDiagPerCell" +electrons.num_particles_per_cell = 10 +electrons.profile = constant +electrons.density = 1.e25 # number of electrons per m^3 +electrons.momentum_distribution_type = "gaussian" +electrons.u_th = 0.01 # uth the std of the (unitless) momentum + diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index 02de34fda..b8a15f9ca 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -65,7 +65,7 @@ useOMP = 0 numthreads = 2 compileTest = 0 doVis = 0 -runtime_params = langmuirwave.ux=0.01 langmuirwave.particle_xmax=0.e-6 +runtime_params = electrons.ux=0.01 electrons.xmax=0.e-6 [Langmuir_x] buildDir = Exec/Langmuir @@ -78,7 +78,7 @@ useOMP = 0 numthreads = 2 compileTest = 0 doVis = 0 -runtime_params = langmuirwave.ux=0.01 langmuirwave.particle_xmax=0.e-6 +runtime_params = electrons.ux=0.01 electrons.xmax=0.e-6 [Langmuir_y] buildDir = Exec/Langmuir @@ -91,7 +91,7 @@ useOMP = 0 numthreads = 2 compileTest = 0 doVis = 0 -runtime_params = langmuirwave.uy=0.01 langmuirwave.particle_ymax=0.e-6 +runtime_params = electrons.uy=0.01 electrons.ymax=0.e-6 [Langmuir_z] buildDir = Exec/Langmuir @@ -104,7 +104,7 @@ useOMP = 0 numthreads = 2 compileTest = 0 doVis = 0 -runtime_params = langmuirwave.uz=0.01 langmuirwave.particle_zmax=0.e-6 +runtime_params = electrons.uz=0.01 electrons.zmax=0.e-6 [Langmuir_multi] buildDir = Exec/Langmuir diff --git a/Source/CustomDensityProb.cpp b/Source/CustomDensityProb.cpp new file mode 100644 index 000000000..1df3d75ad --- /dev/null +++ b/Source/CustomDensityProb.cpp @@ -0,0 +1,13 @@ +#include <PlasmaInjector.H> + +#include <iostream> + +#include <AMReX.H> + +amrex::Real CustomDensityProfile::getDensity(amrex::Real x, + amrex::Real y, + amrex::Real z) const +{ + amrex::Abort("If running with a custom density profile, you must supply a CustomDensityProb.cpp file"); + return 0.0; +} diff --git a/Source/Make.package b/Source/Make.package index fae5bdc91..23968c898 100644 --- a/Source/Make.package +++ b/Source/Make.package @@ -6,7 +6,7 @@ CEXE_sources += WarpXWrappers.cpp CEXE_sources += WarpX.cpp WarpXInitData.cpp WarpXEvolve.cpp WarpXIO.cpp WarpXProb.cpp WarpXRegrid.cpp -CEXE_sources += ParticleProb.cpp ParticleIO.cpp +CEXE_sources += ParticleIO.cpp CEXE_sources += ParticleContainer.cpp WarpXParticleContainer.cpp PhysicalParticleContainer.cpp LaserParticleContainer.cpp CEXE_headers += WarpXWrappers.h @@ -15,6 +15,9 @@ CEXE_headers += WarpX.H WarpX_f.H WarpXConst.H CEXE_headers += ParticleContainer.H WarpXParticleContainer.H PhysicalParticleContainer.H LaserParticleContainer.H +CEXE_headers += PlasmaInjector.H +CEXE_sources += PlasmaInjector.cpp CustomDensityProb.cpp + F90EXE_sources += WarpX_f.F90 WarpX_picsar.F90 WarpX_laser.F90 ifeq ($(USE_OPENBC_POISSON),TRUE) diff --git a/Source/ParticleContainer.H b/Source/ParticleContainer.H index 8d5c4b7e0..20946e6e5 100644 --- a/Source/ParticleContainer.H +++ b/Source/ParticleContainer.H @@ -73,6 +73,10 @@ public: int nSpecies() {return nspecies;} +protected: + + std::vector<std::string> species_names; + private: // physical particles (+ laser) @@ -82,5 +86,6 @@ private: // runtime parameters int nspecies = 1; // physical particles only. If WarpX::use_laser, nspecies+1 == allcontainers.size(). + }; #endif /*WARPX_ParticleContainer_H_*/ diff --git a/Source/ParticleContainer.cpp b/Source/ParticleContainer.cpp index b27534feb..07cc8de09 100644 --- a/Source/ParticleContainer.cpp +++ b/Source/ParticleContainer.cpp @@ -13,7 +13,7 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) int n = WarpX::use_laser ? nspecies+1 : nspecies; allcontainers.resize(n); for (int i = 0; i < nspecies; ++i) { - allcontainers[i].reset(new PhysicalParticleContainer(amr_core, i)); + allcontainers[i].reset(new PhysicalParticleContainer(amr_core, i, species_names[i])); } if (WarpX::use_laser) { allcontainers[n-1].reset(new LaserParticleContainer(amr_core,n-1)); @@ -31,6 +31,11 @@ MultiParticleContainer::ReadParameters () pp.query("nspecies", nspecies); BL_ASSERT(nspecies >= 0); + if (nspecies > 0) { + pp.getarr("species_names", species_names); + BL_ASSERT(species_names.size() == nspecies); + } + initialized = true; } } diff --git a/Source/ParticleProb.cpp b/Source/ParticleProb.cpp deleted file mode 100644 index 4abec6ea8..000000000 --- a/Source/ParticleProb.cpp +++ /dev/null @@ -1,21 +0,0 @@ - -// -// Each problem must have its own version of PhysicalParticleContainer::InitData() -// to initialize the particle data. It must also initialize charge and mass. -// - -#include <cmath> - -#include <AMReX_BLProfiler.H> - -#include <ParticleContainer.H> -#include <WarpXConst.H> - -using namespace amrex; - -void -PhysicalParticleContainer::InitData() -{ - static_assert(false, - "Each problem must have its own version of PhysicalParticleContainer::InitData()"); -} diff --git a/Source/PhysicalParticleContainer.H b/Source/PhysicalParticleContainer.H index c6d13185e..40a8a6068 100644 --- a/Source/PhysicalParticleContainer.H +++ b/Source/PhysicalParticleContainer.H @@ -3,13 +3,16 @@ #include <map> +#include <PlasmaInjector.H> #include <WarpXParticleContainer.H> class PhysicalParticleContainer : public WarpXParticleContainer { public: - PhysicalParticleContainer (amrex::AmrCore* amr_core, int ispecies); + PhysicalParticleContainer (amrex::AmrCore* amr_core, + int ispecies, + const std::string& name); virtual ~PhysicalParticleContainer () {} virtual void AllocData () override; @@ -17,15 +20,39 @@ public: virtual void InitData () override; virtual void FieldGather(int lev, - const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, const amrex::MultiFab& Ez, - const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz) override; + const amrex::MultiFab& Ex, + const amrex::MultiFab& Ey, + const amrex::MultiFab& Ez, + const amrex::MultiFab& Bx, + const amrex::MultiFab& By, + const amrex::MultiFab& Bz) override; virtual void Evolve (int lev, - const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, const amrex::MultiFab& Ez, - const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz, - amrex::MultiFab& jx, amrex::MultiFab& jy, amrex::MultiFab& jz, amrex::Real t, amrex::Real dt) override; + const amrex::MultiFab& Ex, + const amrex::MultiFab& Ey, + const amrex::MultiFab& Ez, + const amrex::MultiFab& Bx, + const amrex::MultiFab& By, + const amrex::MultiFab& Bz, + amrex::MultiFab& jx, + amrex::MultiFab& jy, + amrex::MultiFab& jz, + amrex::Real t, + amrex::Real dt) override; virtual void PostRestart () override {} + + // Inject particles in Box 'part_box' + void AddParticles (int lev, amrex::Box part_box = amrex::Box()); + +protected: + + void AddNDiagPerCell (int lev, amrex::Box part_box = amrex::Box()); + void AddNRandomUniformPerCell (int lev, amrex::Box part_box = amrex::Box()); + void AddNRandomNormal (int lev, amrex::Box part_box = amrex::Box()); + + std::string species_name; + std::unique_ptr<PlasmaInjector> plasma_injector; }; #endif diff --git a/Source/PhysicalParticleContainer.cpp b/Source/PhysicalParticleContainer.cpp index 9e2532403..d499f7dfb 100644 --- a/Source/PhysicalParticleContainer.cpp +++ b/Source/PhysicalParticleContainer.cpp @@ -1,4 +1,5 @@ #include <limits> +#include <sstream> #include <ParticleContainer.H> #include <WarpX_f.H> @@ -6,9 +7,16 @@ using namespace amrex; -PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int ispecies) - : WarpXParticleContainer(amr_core, ispecies) +PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int ispecies, + const std::string& name) + : WarpXParticleContainer(amr_core, ispecies), + species_name(name) { + plasma_injector.reset(new PlasmaInjector(species_id, species_name)); +} + +void PhysicalParticleContainer::InitData() { + AddParticles(0); // Note - only one level right now } void @@ -21,6 +29,152 @@ PhysicalParticleContainer::AllocData () } void +PhysicalParticleContainer::AddNRandomUniformPerCell (int lev, Box part_box) { + BL_PROFILE("PhysicalParticleContainer::AddNRandomPerCell()"); + + charge = plasma_injector->getCharge(); + mass = plasma_injector->getMass(); + + const Geometry& geom = Geom(lev); + const Real* dx = geom.CellSize(); + + if (!part_box.ok()) part_box = geom.Domain(); + + Real scale_fac; + int n_part_per_cell = plasma_injector->numParticlesPerCell(); + +#if BL_SPACEDIM==3 + scale_fac = dx[0]*dx[1]*dx[2]/n_part_per_cell; +#elif BL_SPACEDIM==2 + scale_fac = dx[0]*dx[1]/n_part_per_cell; +#endif + + std::array<Real,PIdx::nattribs> attribs; + attribs.fill(0.0); + + // Initialize random generator for normal distribution + std::default_random_engine generator; + std::uniform_real_distribution<Real> position_distribution(0.0,1.0); + for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) + { + const Box& tile_box = mfi.tilebox(); + const Box& intersectBox = tile_box & part_box; + if (!intersectBox.ok()) continue; + + RealBox tile_real_box { intersectBox, dx, geom.ProbLo() }; + + const int grid_id = mfi.index(); + const int tile_id = mfi.LocalTileIndex(); + + const auto& boxlo = intersectBox.smallEnd(); + for (IntVect iv = intersectBox.smallEnd(); iv <= intersectBox.bigEnd(); intersectBox.next(iv)) + { + for (int i_part=0; i_part<n_part_per_cell;i_part++) + { + // Randomly generate the positions (uniformly inside each cell) + Real particle_shift_x = position_distribution(generator); + Real particle_shift_y = position_distribution(generator); + Real particle_shift_z = position_distribution(generator); + +#if (BL_SPACEDIM == 3) + Real x = tile_real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift_x)*dx[0]; + Real y = tile_real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift_y)*dx[1]; + Real z = tile_real_box.lo(2) + (iv[2]-boxlo[2] + particle_shift_z)*dx[2]; +#elif (BL_SPACEDIM == 2) + Real x = tile_real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift_x)*dx[0]; + Real y = 0.0; + Real z = tile_real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift_z)*dx[1]; +#endif + + if (plasma_injector->insideBounds(x, y, z)) { + Real weight; + std::array<Real, 3> u; + weight = plasma_injector->getDensity(x, y, z) * scale_fac; + plasma_injector->getMomentum(u); + attribs[PIdx::w ] = weight; + attribs[PIdx::ux] = u[0]; + attribs[PIdx::uy] = u[1]; + attribs[PIdx::uz] = u[2]; + AddOneParticle(lev, grid_id, tile_id, x, y, z, attribs); + } + } + } + } +} + +void +PhysicalParticleContainer::AddNRandomNormal (int lev, Box part_box) { + BL_PROFILE("PhysicalParticleContainer::AddNRandomNormal()"); + amrex::Abort("PhysicalParticleContainer::AddNRandomNormal() not implemented yet."); +} + +void +PhysicalParticleContainer::AddNDiagPerCell (int lev, Box part_box) { + BL_PROFILE("PhysicalParticleContainer::AddNDiagPerCell()"); + + charge = plasma_injector->getCharge(); + mass = plasma_injector->getMass(); + + const Geometry& geom = Geom(lev); + const Real* dx = geom.CellSize(); + + if (!part_box.ok()) part_box = geom.Domain(); + + Real scale_fac; + int n_part_per_cell = plasma_injector->numParticlesPerCell(); + +#if BL_SPACEDIM==3 + scale_fac = dx[0]*dx[1]*dx[2]/n_part_per_cell; +#elif BL_SPACEDIM==2 + scale_fac = dx[0]*dx[1]/n_part_per_cell; +#endif + + std::array<Real,PIdx::nattribs> attribs; + attribs.fill(0.0); + for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) { + const Box& tile_box = mfi.tilebox(); + const Box& intersectBox = tile_box & part_box; + if (!intersectBox.ok()) continue; + + RealBox tile_real_box { intersectBox, dx, geom.ProbLo() }; + + const int grid_id = mfi.index(); + const int tile_id = mfi.LocalTileIndex(); + + const auto& boxlo = intersectBox.smallEnd(); + for (IntVect iv = intersectBox.smallEnd(); iv <= intersectBox.bigEnd(); intersectBox.next(iv)) + { + for (int i_part=0; i_part<n_part_per_cell;i_part++) + { + Real particle_shift = (0.5+i_part)/n_part_per_cell; +#if (BL_SPACEDIM == 3) + Real x = tile_real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift)*dx[0]; + Real y = tile_real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1]; + Real z = tile_real_box.lo(2) + (iv[2]-boxlo[2] + particle_shift)*dx[2]; +#elif (BL_SPACEDIM == 2) + Real x = tile_real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift)*dx[0]; + Real y = 0.0; + Real z = tile_real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1]; +#endif + + if (plasma_injector->insideBounds(x, y, z)) { + Real weight; + std::array<Real, 3> u; + weight = plasma_injector->getDensity(x, y, z) * scale_fac; + plasma_injector->getMomentum(u); + attribs[PIdx::w ] = weight; + attribs[PIdx::ux] = u[0]; + attribs[PIdx::uy] = u[1]; + attribs[PIdx::uz] = u[2]; + AddOneParticle(lev, grid_id, tile_id, x, y, z, attribs); + } + } + } + } +} + + +void PhysicalParticleContainer::FieldGather (int lev, const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez, const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz) @@ -128,6 +282,17 @@ PhysicalParticleContainer::FieldGather (int lev, } void +PhysicalParticleContainer::AddParticles (int lev, Box part_box) +{ + if (plasma_injector->injection_style == "ndiagpercell") + AddNDiagPerCell(lev, part_box); + else if (plasma_injector->injection_style == "nrandomuniformpercell") + AddNRandomUniformPerCell(lev, part_box); + else if (plasma_injector->injection_style == "nrandomnormal") + AddNRandomNormal(lev, part_box); +} + +void PhysicalParticleContainer::Evolve (int lev, const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez, const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz, diff --git a/Source/PlasmaInjector.H b/Source/PlasmaInjector.H new file mode 100644 index 000000000..c346f79da --- /dev/null +++ b/Source/PlasmaInjector.H @@ -0,0 +1,161 @@ +#ifndef PLASMA_INJECTOR_H_ +#define PLASMA_INJECTOR_H_ + +#include <random> +#include <array> + +#include "AMReX_REAL.H" +#include "AMReX_Array.H" +#include "AMReX_ParmParse.H" + +/// +/// PlasmaDensityProfile describes how the charge density +/// is set in particle initialization. Subclasses must define a +/// getDensity function that describes the charge density as a +/// function of x, y, and z. +/// +class PlasmaDensityProfile +{ +public: + virtual ~PlasmaDensityProfile() {}; + virtual amrex::Real getDensity(amrex::Real x, + amrex::Real y, + amrex::Real z) const = 0; +protected: + std::string _species_name; +}; + +/// +/// This describes a constant density distribution. +/// +class ConstantDensityProfile : public PlasmaDensityProfile +{ +public: + ConstantDensityProfile(amrex::Real _density); + virtual amrex::Real getDensity(amrex::Real x, + amrex::Real y, + amrex::Real z) const override; + +private: + amrex::Real _density; +}; + +/// +/// This describes a custom density distribution. Users can supply +/// in their problem directory. +/// +/// +class CustomDensityProfile : public PlasmaDensityProfile +{ +public: + CustomDensityProfile(const std::string& species_name); + virtual amrex::Real getDensity(amrex::Real x, + amrex::Real y, + amrex::Real z) const override; +private: + amrex::Array<amrex::Real> params; +}; + +/// +/// PlasmaMomentumDistribution describes how the particle momenta +/// are set. Subclasses must define a getMomentum method that fills +/// a u with the 3 components of the particle momentum +/// +class PlasmaMomentumDistribution +{ +public: + using vec3 = std::array<amrex::Real, 3>; + virtual ~PlasmaMomentumDistribution() {}; + virtual void getMomentum(vec3& u) = 0; +}; + +/// +/// This is a constant momentum distribution - all particles will +/// have the same ux, uy, and uz +/// +class ConstantMomentumDistribution : public PlasmaMomentumDistribution +{ +public: + ConstantMomentumDistribution(amrex::Real ux, + amrex::Real uy, + amrex::Real uz); + virtual void getMomentum(vec3& u) override; + +private: + amrex::Real _ux; + amrex::Real _uy; + amrex::Real _uz; +}; + +/// +/// This is a Gaussian Random momentum distribution. +/// Particles will get random momenta, drawn from a normal. +/// ux_m, ux_y, and ux_z describe the mean components in the x, y, and z +/// directions, while u_th is the standard deviation of the random +/// component. +/// +class GaussianRandomMomentumDistribution : public PlasmaMomentumDistribution +{ +public: + GaussianRandomMomentumDistribution(amrex::Real ux_m, + amrex::Real uy_m, + amrex::Real uz_m, + amrex::Real u_th); + virtual void getMomentum(vec3& u) override; +private: + amrex::Real _ux_m; + amrex::Real _uy_m; + amrex::Real _uz_m; + amrex::Real _u_th; + + // xxxxx These random number generators are not thread safe. + std::normal_distribution<amrex::Real> momentum_distribution; + std::default_random_engine generator; +}; + +/// +/// The PlasmaInjector class parses and stores information about the plasma +/// type used in the particle container. This information is used to create the +/// particles on initialization and whenever the window moves. +/// +class PlasmaInjector +{ + +public: + + using vec3 = std::array<amrex::Real, 3>; + + PlasmaInjector(int ispecies, const std::string& name); + + amrex::Real getDensity(amrex::Real x, amrex::Real y, amrex::Real z); + + bool insideBounds(amrex::Real x, amrex::Real y, amrex::Real z); + + int numParticlesPerCell() {return num_particles_per_cell;}; + + void getMomentum(vec3& u); + + amrex::Real getCharge() {return charge;} + amrex::Real getMass() {return mass;} + + std::string injection_style; + +protected: + + amrex::Real mass, charge; + + amrex::Real xmin, xmax; + amrex::Real ymin, ymax; + amrex::Real zmin, zmax; + + amrex::Real density; + int num_particles_per_cell; + + int species_id; + std::string species_name; + + std::unique_ptr<PlasmaDensityProfile> rho_prof; + std::unique_ptr<PlasmaMomentumDistribution> mom_dist; +}; + +#endif diff --git a/Source/PlasmaInjector.cpp b/Source/PlasmaInjector.cpp new file mode 100644 index 000000000..027bb5b57 --- /dev/null +++ b/Source/PlasmaInjector.cpp @@ -0,0 +1,208 @@ +#include "PlasmaInjector.H" + +#include <sstream> + +#include <WarpXConst.H> +#include <AMReX.H> + +using namespace amrex; + +namespace { + void StringParseAbortMessage(const std::string& var, + const std::string& name) { + std::stringstream stringstream; + std::string string; + stringstream << var << " string '" << name << "' not recognized."; + string = stringstream.str(); + amrex::Abort(string.c_str()); + } + + Real parseChargeName(const std::string& name) { + if (name == "q_e") { + return PhysConst::q_e; + } else { + StringParseAbortMessage("Charge", name); + return 0.0; + } + } + + Real parseChargeString(const std::string& name) { + if(name.substr(0, 1) == "-") + return -1.0 * parseChargeName(name.substr(1, name.size() - 1)); + return parseChargeName(name); + } + + Real parseMassString(const std::string& name) { + if (name == "m_e") { + return PhysConst::m_e; + } else if (name == "m_p"){ + return PhysConst::m_p; + } else { + StringParseAbortMessage("Mass", name); + return 0.0; + } + } +} + +ConstantDensityProfile::ConstantDensityProfile(Real density) + : _density(density) +{} + +Real ConstantDensityProfile::getDensity(Real x, Real y, Real z) const +{ + return _density; +} + +CustomDensityProfile::CustomDensityProfile(const std::string& species_name) +{ + ParmParse pp(species_name); + pp.getarr("custom_profile_params", params); +} + +ConstantMomentumDistribution::ConstantMomentumDistribution(Real ux, + Real uy, + Real uz) + : _ux(ux), _uy(uy), _uz(uz) +{} + +void ConstantMomentumDistribution::getMomentum(vec3& u) { + u[0] = _ux; + u[1] = _uy; + u[2] = _uz; +} + +GaussianRandomMomentumDistribution::GaussianRandomMomentumDistribution(Real ux_m, + Real uy_m, + Real uz_m, + Real u_th) + : _ux_m(ux_m), _uy_m(uy_m), _uz_m(uz_m), _u_th(u_th), + momentum_distribution(0.0, u_th) +{ +} + +void GaussianRandomMomentumDistribution::getMomentum(vec3& u) { + Real ux_th = momentum_distribution(generator); + Real uy_th = momentum_distribution(generator); + Real uz_th = momentum_distribution(generator); + u[0] = _ux_m + ux_th; + u[1] = _uy_m + uy_th; + u[2] = _uz_m + uz_th; +} + +PlasmaInjector::PlasmaInjector(int ispecies, const std::string& name) + : species_id(ispecies), species_name(name) +{ + ParmParse pp(species_name); + + // parse charge and mass + std::string charge_s; + pp.get("charge", charge_s); + std::transform(charge_s.begin(), + charge_s.end(), + charge_s.begin(), + ::tolower); + charge = parseChargeString(charge_s); + + std::string mass_s; + pp.get("mass", mass_s); + std::transform(mass_s.begin(), + mass_s.end(), + mass_s.begin(), + ::tolower); + mass = parseMassString(mass_s); + + // parse plasma boundaries + xmin = std::numeric_limits<amrex::Real>::lowest(); + ymin = std::numeric_limits<amrex::Real>::lowest(); + zmin = std::numeric_limits<amrex::Real>::lowest(); + + xmax = std::numeric_limits<amrex::Real>::max(); + ymax = std::numeric_limits<amrex::Real>::max(); + zmax = std::numeric_limits<amrex::Real>::max(); + + pp.query("xmin", xmin); + pp.query("ymin", ymin); + pp.query("zmin", zmin); + pp.query("xmax", xmax); + pp.query("ymax", ymax); + pp.query("zmax", zmax); + + // parse density information + std::string rho_prof_s; + pp.get("profile", rho_prof_s); + std::transform(rho_prof_s.begin(), + rho_prof_s.end(), + rho_prof_s.begin(), + ::tolower); + if (rho_prof_s == "constant") { + Real density; + pp.get("density", density); + rho_prof.reset(new ConstantDensityProfile(density)); + } else if (rho_prof_s == "custom") { + rho_prof.reset(new CustomDensityProfile(species_name)); + } else { + StringParseAbortMessage("Density profile type", rho_prof_s); + } + pp.get("num_particles_per_cell", num_particles_per_cell); + + // parse momentum information + std::string mom_dist_s; + pp.get("momentum_distribution_type", mom_dist_s); + std::transform(mom_dist_s.begin(), + mom_dist_s.end(), + mom_dist_s.begin(), + ::tolower); + if (mom_dist_s == "constant") { + Real ux = 0.; + Real uy = 0.; + Real uz = 0.; + pp.query("ux", ux); + pp.query("uy", uy); + pp.query("uz", uz); + mom_dist.reset(new ConstantMomentumDistribution(ux, uy, uz)); + } else if (mom_dist_s == "gaussian") { + Real ux_m = 0.; + Real uy_m = 0.; + Real uz_m = 0.; + Real u_th = 0.; + pp.query("ux_m", ux_m); + pp.query("uy_m", uy_m); + pp.query("uz_m", uz_m); + pp.query("u_th", u_th); + mom_dist.reset(new GaussianRandomMomentumDistribution(ux_m, uy_m, uz_m, u_th)); + } else { + StringParseAbortMessage("Momentum distribution type", mom_dist_s); + } + + // get injection style + pp.get("injection_style", injection_style); + std::transform(injection_style.begin(), + injection_style.end(), + injection_style.begin(), + ::tolower); + if (injection_style == "nrandomnormal" or + injection_style == "nrandomuniformpercell" or + injection_style == "ndiagpercell") { + return; + } else { + StringParseAbortMessage("Injection style", injection_style); + } +} + +void PlasmaInjector::getMomentum(vec3& u) { + mom_dist->getMomentum(u); + u[0] *= PhysConst::c; + u[1] *= PhysConst::c; + u[2] *= PhysConst::c; +} + +bool PlasmaInjector::insideBounds(Real x, Real y, Real z) { + if (x >= xmax || x < xmin || + y >= ymax || y < ymin || + z >= zmax || z < zmin ) return false; + return true; +} + +Real PlasmaInjector::getDensity(Real x, Real y, Real z) { + return rho_prof->getDensity(x, y, z); +} diff --git a/Source/WarpX.H b/Source/WarpX.H index 47cc271c8..fad31e1bf 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -36,6 +36,15 @@ public: MultiParticleContainer& GetPartContainer () { return *mypc; } +#if (BL_SPACEDIM == 3) + void fillSlice(amrex::Real z_coord) const; + + void sampleAtPoints(const amrex::Array<amrex::Real>& x, + const amrex::Array<amrex::Real>& y, + const amrex::Array<amrex::Real>& z, + amrex::Array<amrex::Array<amrex::Real> >& result) const; +#endif + static void FillBoundary (amrex::MultiFab& mf, const amrex::Geometry& geom, const amrex::IntVect& nodalflag); static void shiftMF(amrex::MultiFab& mf, const amrex::Geometry& geom, int num_shift, @@ -113,6 +122,8 @@ protected: //! Delete level data. Called by AmrCore::regrid. virtual void ClearLevel (int lev) override; + amrex::Box getIndexBox(const amrex::RealBox& real_box) const; + private: // Singleton is used when the code is run from python @@ -157,6 +168,7 @@ private: // Particle container std::unique_ptr<MultiParticleContainer> mypc; + // Fields: First array for level, second for direction amrex::Array<amrex::Array< std::unique_ptr<amrex::MultiFab> > > current; amrex::Array<amrex::Array< std::unique_ptr<amrex::MultiFab> > > Efield; @@ -171,9 +183,7 @@ private: // Plasma injection parameters int do_plasma_injection = 0; int num_injected_species = -1; - amrex::Array<int> injected_plasma_ppc; amrex::Array<int> injected_plasma_species; - amrex::Array<amrex::Real> injected_plasma_density; // Other runtime parameters int verbose = 1; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 052b5ed09..13c46ec39 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -165,15 +165,9 @@ WarpX::ReadParameters () pp.query("do_plasma_injection", do_plasma_injection); if (do_plasma_injection) { pp.get("num_injected_species", num_injected_species); - injected_plasma_ppc.resize(num_injected_species); - pp.getarr("injected_plasma_ppc", injected_plasma_ppc, - 0, num_injected_species); injected_plasma_species.resize(num_injected_species); pp.getarr("injected_plasma_species", injected_plasma_species, 0, num_injected_species); - injected_plasma_density.resize(num_injected_species); - pp.getarr("injected_plasma_density", injected_plasma_density, - 0, num_injected_species); } pp.query("use_laser", use_laser); @@ -312,3 +306,157 @@ WarpX::Copy(MultiFab& dstmf, int dcomp, int ncomp, const MultiFab& srcmf, int sc dfab.SetBoxType(dst_typ); } } + +Box +WarpX::getIndexBox(const RealBox& real_box) const { + BL_ASSERT(max_level == 0); + + IntVect slice_lo, slice_hi; + + D_TERM(slice_lo[0]=floor((real_box.lo(0) - geom[0].ProbLo(0))/geom[0].CellSize(0));, + slice_lo[1]=floor((real_box.lo(1) - geom[0].ProbLo(1))/geom[0].CellSize(1));, + slice_lo[2]=floor((real_box.lo(2) - geom[0].ProbLo(2))/geom[0].CellSize(2));); + + D_TERM(slice_hi[0]=floor((real_box.hi(0) - geom[0].ProbLo(0))/geom[0].CellSize(0));, + slice_hi[1]=floor((real_box.hi(1) - geom[0].ProbLo(1))/geom[0].CellSize(1));, + slice_hi[2]=floor((real_box.hi(2) - geom[0].ProbLo(2))/geom[0].CellSize(2));); + + return Box(slice_lo, slice_hi) & geom[0].Domain(); +} + +#if (BL_SPACEDIM == 3) +void +WarpX::fillSlice(Real z_coord) const { + BL_ASSERT(max_level == 0); + + // Get our slice and convert to index space + RealBox real_slice = geom[0].ProbDomain(); + real_slice.setLo(2, z_coord); + real_slice.setHi(2, z_coord); + Box slice_box = getIndexBox(real_slice); + + // define the multifab that stores slice + BoxArray ba = Efield[0][0]->boxArray(); + const IndexType correct_typ(IntVect::TheZeroVector()); + ba.convert(correct_typ); + const DistributionMapping& dm = dmap[0]; + std::vector< std::pair<int, Box> > isects; + ba.intersections(slice_box, isects, false, 0); + Array<Box> boxes; + Array<int> procs; + for (int i = 0; i < isects.size(); ++i) { + procs.push_back(dm[isects[i].first]); + boxes.push_back(isects[i].second); + } + procs.push_back(ParallelDescriptor::MyProc()); + BoxArray slice_ba(&boxes[0], boxes.size()); + DistributionMapping slice_dmap(procs); + MultiFab slice(slice_ba, slice_dmap, 6, 0); + + const MultiFab* mfs[6]; + mfs[0] = Efield[0][0].get(); + mfs[1] = Efield[0][1].get(); + mfs[2] = Efield[0][2].get(); + mfs[3] = Bfield[0][0].get(); + mfs[4] = Bfield[0][1].get(); + mfs[5] = Bfield[0][2].get(); + + IntVect flags[6]; + flags[0] = WarpX::Ex_nodal_flag; + flags[1] = WarpX::Ey_nodal_flag; + flags[2] = WarpX::Ez_nodal_flag; + flags[3] = WarpX::Bx_nodal_flag; + flags[4] = WarpX::By_nodal_flag; + flags[5] = WarpX::Bz_nodal_flag; + + // Fill the slice with sampled data + const Real* dx = geom[0].CellSize(); + const Geometry& g = geom[0]; + for (MFIter mfi(slice); mfi.isValid(); ++mfi) { + int slice_gid = mfi.index(); + Box grid = slice_ba[slice_gid]; + int xlo = grid.smallEnd(0), ylo = grid.smallEnd(1), zlo = grid.smallEnd(2); + int xhi = grid.bigEnd(0), yhi = grid.bigEnd(1), zhi = grid.bigEnd(2); + + IntVect iv = grid.smallEnd(); + ba.intersections(Box(iv, iv), isects, true, 0); + BL_ASSERT(!isects.empty()); + int full_gid = isects[0].first; + + for (int k = zlo; k <= zhi; k++) { + for (int j = ylo; j <= yhi; j++) { + for (int i = xlo; i <= xhi; i++) { + for (int comp = 0; comp < 6; comp++) { + Real x = g.ProbLo(0) + i*dx[0]; + Real y = g.ProbLo(1) + j*dx[1]; + Real z = z_coord; + + D_TERM(iv[0]=floor((x - g.ProbLo(0) + 0.5*g.CellSize(0)*flags[comp][0])/g.CellSize(0));, + iv[1]=floor((y - g.ProbLo(1) + 0.5*g.CellSize(1)*flags[comp][1])/g.CellSize(1));, + iv[2]=floor((z - g.ProbLo(2) + 0.5*g.CellSize(2)*flags[comp][2])/g.CellSize(2));); + + slice[slice_gid](IntVect(i, j, k), comp) = (*mfs[comp])[full_gid](iv); + } + } + } + } + } +} + +void WarpX::sampleAtPoints(const Array<Real>& x, + const Array<Real>& y, + const Array<Real>& z, + Array<Array<Real> >& result) const { + BL_ASSERT((x.size() == y.size()) and (y.size() == z.size())); + BL_ASSERT(max_level == 0); + + const MultiFab* mfs[6]; + mfs[0] = Efield[0][0].get(); + mfs[1] = Efield[0][1].get(); + mfs[2] = Efield[0][2].get(); + mfs[3] = Bfield[0][0].get(); + mfs[4] = Bfield[0][1].get(); + mfs[5] = Bfield[0][2].get(); + + IntVect flags[6]; + flags[0] = WarpX::Ex_nodal_flag; + flags[1] = WarpX::Ey_nodal_flag; + flags[2] = WarpX::Ez_nodal_flag; + flags[3] = WarpX::Bx_nodal_flag; + flags[4] = WarpX::By_nodal_flag; + flags[5] = WarpX::Bz_nodal_flag; + + const unsigned npoints = x.size(); + const int ncomps = 6; + result.resize(ncomps); + for (int i = 0; i < ncomps; i++) { + result[i].resize(npoints, 0); + } + + BoxArray ba = Efield[0][0]->boxArray(); + const IndexType correct_typ(IntVect::TheZeroVector()); + ba.convert(correct_typ); + std::vector< std::pair<int, Box> > isects; + + IntVect iv; + const Geometry& g = geom[0]; + for (int i = 0; i < npoints; ++i) { + for (int comp = 0; comp < ncomps; comp++) { + D_TERM(iv[0]=floor((x[i] - g.ProbLo(0) + 0.5*g.CellSize(0)*flags[comp][0])/g.CellSize(0));, + iv[1]=floor((y[i] - g.ProbLo(1) + 0.5*g.CellSize(1)*flags[comp][1])/g.CellSize(1));, + iv[2]=floor((z[i] - g.ProbLo(2) + 0.5*g.CellSize(2)*flags[comp][2])/g.CellSize(2));); + + ba.intersections(Box(iv, iv), isects, true, 0); + const int grid = isects[0].first; + const int who = dmap[0][grid]; + if (who == ParallelDescriptor::MyProc()) { + result[comp][i] = (*mfs[comp])[grid](iv); + } + } + } + + for (int i = 0; i < ncomps; i++) { + ParallelDescriptor::ReduceRealSum(result[i].dataPtr(), result[i].size()); + } +} +#endif diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index 16b957dec..29336e66b 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/WarpXEvolve.cpp @@ -284,20 +284,11 @@ WarpX::InjectPlasma (int num_shift, int dir) const Real* dx = geom[lev].CellSize(); - for (int i = 0; i < num_injected_species; ++i) - { - int ppc = injected_plasma_ppc[i]; - Real density = injected_plasma_density[i]; -#if BL_SPACEDIM==3 - Real weight = density * dx[0]*dx[1]*dx[2]/ppc; -#elif BL_SPACEDIM==2 - Real weight = density * dx[0]*dx[1]/ppc; -#endif - + for (int i = 0; i < num_injected_species; ++i) { int ispecies = injected_plasma_species[i]; WarpXParticleContainer& pc = mypc->GetParticleContainer(ispecies); - - pc.AddParticles(lev, particleBox, weight, ppc); + auto& ppc = dynamic_cast<PhysicalParticleContainer&>(pc); + ppc.AddParticles(lev, particleBox); } } } diff --git a/Source/WarpXParticleContainer.H b/Source/WarpXParticleContainer.H index ac0ebfd8c..12152c286 100644 --- a/Source/WarpXParticleContainer.H +++ b/Source/WarpXParticleContainer.H @@ -69,9 +69,6 @@ public: amrex::Real x, amrex::Real y, amrex::Real z, const std::array<amrex::Real,PIdx::nattribs>& attribs); - // Inject 'ppc' particles per cell in Box 'part_box' - void AddParticles (int lev, const amrex::Box& part_box, amrex::Real weight, int ppc); - void ReadHeader (std::istream& is); void WriteHeader (std::ostream& os) const; diff --git a/Source/WarpXParticleContainer.cpp b/Source/WarpXParticleContainer.cpp index 2965824ff..31e4702ed 100644 --- a/Source/WarpXParticleContainer.cpp +++ b/Source/WarpXParticleContainer.cpp @@ -66,49 +66,6 @@ WarpXParticleContainer::AddOneParticle (ParticleTileType& particle_tile, } void -WarpXParticleContainer::AddParticles (int lev, const Box& part_box, Real weight, int ppc) -{ - const Geometry& geom = Geom(lev); - const Real* dx = geom.CellSize(); - - std::array<Real,PIdx::nattribs> attribs; - attribs.fill(0.0); - attribs[PIdx::w] = weight; - - for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) - { - const Box& tile_box = mfi.tilebox(); - const Box& intersectBox = tile_box & part_box; - if (intersectBox.ok()) - { - RealBox real_box { intersectBox, dx, geom.ProbLo() }; - - const int grid_id = mfi.index(); - const int tile_id = mfi.LocalTileIndex(); - - const IntVect& boxlo = tile_box.smallEnd(); - for (IntVect iv = boxlo; iv <= tile_box.bigEnd(); tile_box.next(iv)) - { - for (int i_part=0; i_part < ppc; i_part++) - { - Real particle_shift = (0.5+i_part)/ppc; -#if (BL_SPACEDIM == 3) - Real x = real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift)*dx[0]; - Real y = real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1]; - Real z = real_box.lo(2) + (iv[2]-boxlo[2] + particle_shift)*dx[2]; -#elif (BL_SPACEDIM == 2) - Real x = real_box.lo(0) + (iv[0]-boxlo[0] + particle_shift)*dx[0]; - Real y = 0.0; - Real z = real_box.lo(1) + (iv[1]-boxlo[1] + particle_shift)*dx[1]; -#endif - AddOneParticle(lev, grid_id, tile_id, x, y, z, attribs); - } - } - } - } -} - -void WarpXParticleContainer::AddNParticles (int n, const Real* x, const Real* y, const Real* z, const Real* vx, const Real* vy, const Real* vz, int nattr, const Real* attr, int uniqueparticles) diff --git a/Source/WarpXRegrid.cpp b/Source/WarpXRegrid.cpp index 32508b9ba..d3eced1e6 100644 --- a/Source/WarpXRegrid.cpp +++ b/Source/WarpXRegrid.cpp @@ -217,7 +217,7 @@ WarpX::getCostCountDM (const Array<long>& cost, const BoxArray& ba) { DistributionMapping res; int nprocs = ParallelDescriptor::NProcs(); - const int factor = 1.5; // A process can get up to 'factor' times of the average number of boxes. + const Real factor = 1.5; // A process can get up to 'factor' times of the average number of boxes. int nmax = (cost.size()+nprocs-1) / nprocs * factor; Real eff; res.KnapSackProcessorMap(cost, nprocs, &eff, true, nmax); |