Skip to content

Commit

Permalink
Pulsar add external eb (#39)
Browse files Browse the repository at this point in the history 
* Adding pulsar input and reading pulsar parameters using Don's PR

* added external E and B fields and the appropriate input parameters

* EOL whitespace

* seeting up bounds for pulsar surface

* pushing changes for comparison

* particles injected below surface, dipole magnetic field inside and outside, momentum of all particles is zero, fixed function that applied vacuum EB fields on particles

* add particles based on analytical surface charge density

* using local Er and Er_cor to compute particle introduction

* fixing particle injection with local Er that includes Er_external

* changes to pulsar

* michel and li efields inside and outside star. Particle injection rho_GJ

* adding input parameters to toggle rho gj scaling, ninjection, maxden

* removing unwanted insidePulsarBound function

* Damp electric fields internal to the pulsar when compiling with PULSAR on.

* Separate out cartesian ijk -> spherical radius into a utility function so it can be used elsewhere.

* Add PulsarParm::Omega(t) function that returns omega(t) depending on ramp_omega_time input parameter. If ramp_omega_time < 0, PulsarParm::Omega(t) = omega_star.

* Ran the style commands suggested by Travis CI to remove trailing whitespace.

* Upsating GetPosition() based on recent changes in master branch

* Damp electric fields internal to the pulsar when compiling with PULSAR on.

* Undo line removal.

* Calling External E and B for pulsar from GetExternal functions

* a few fixes in the getexternalEB and addplasma to ensure test works

* fixed Ecell

* removing particles if inside R_Star

* recent input file

* eol fix

* Add functor for pulsar particle injection

* fix eol

* Add Functor file

* fix eol

* add updated input file

* push external-EB-on-grid to debug. particle injection off.

* inject fixed number of particles when injection is on

* adding file with external EM field on grid

* fix indent

* fixing for loop to include last dimension

* add notebook with scaling

* fixing E-corotating insdie. E-quad outside. B-dipole everywhere

* gather external pulsar fields, add num_ppc particles, damp Maxwell fields inside pulsar

* include particles only outside r_star

* fix EOL

* particle injection inner ring and particle removal in same location

* damp EB fields inside R_Star, intro and remove part inside R_star

* fix typo in cell cell z coordinate

* add fraction of particles through pcounts of cell.

* Revert "add fraction of particles through pcounts of cell."

This reverts commit f3daa96.

* ensuring that fractional particle injection is uniform in cell

* we dont need the p.id=-1 to remove particles 1-fractional part injection

* updated input file

* terminate ifdef PULSAR

* remove unnecessary call to add external EB on EB Multifab. we call it directly in Fieldgather instead

* clean commented out lines in GetExternalFields

* remove commented out code and old unused PulsarEField and PulsarBField

* remove unused function init

* removing AddPulsarPlasma functor which is not used

* remove commented out code and calls to add plasma functor

* adding missing ;

* FillBoundary for EB after dampingEB

* using gpuarray for staggering used to damp EB fields

* default values in fieldgather doShapeN for problo, probhi, and cur_time

* empty array

* pass by reference

* take user-defined cc to compute r

* Add Dons recommendation to change wt for fractional injection

* modifying particle wt, redefining insidePulsarBoundsCC, and removing unnecessary functions

* delete redundant function definition

* remove redundant function definitions

* remove unused variable

* remove unnecesary empty line

* adding user-inpu for modifying particle weight

* querying input for modiyfing particle weight

* remove header file decleration in Make package

* remove unnecesary header

* remove unused function call

* use input value of star center to allow for shifted centers in the input

* adding a separate file for gather pulsar external EB using grid resolution

* removing Pulsar related calls from the core gather functions

* call pulsar field gather separate from field gather

* deposit charge is modified and therefore not const

* reorder function definitions and add missingsemicolon

* fix typo of variable in pulsarparameters

* add the empty line back

* undo any changes in FieldGather so it is the same as WarpX development

* remove modifications made to fieldgather

* remove unnecessary blank line in at end of file

* remove tab

* EOL

* input parameters for damping EB and external EB to distinguish inside and outside

* adding a few more input parameters

* adding particle absortion radius as input

* remove unnecessary print statement

* add min max bounds for defining particle injection region

* EOL

* fix typo

* fix comments from the review

Co-authored-by: Donald E. Willcox <[email protected]>
Co-authored-by: Revathi Jambunathan <[email protected]>
  • Loading branch information
@RevathiJambunathan @dwillcox
3 people committed Jan 30, 2021
1 parent b3d1943 commit 035908f
Show file tree
Hide file tree
Showing 19 changed files with 3,142 additions and 441 deletions.
1,392 changes: 1,392 additions & 0 deletions Examples/Physics_applications/pulsar/Pulsar_scale_RstarPhysical.ipynb
View file

Large diffs are not rendered by default.

103 changes: 43 additions & 60 deletions Examples/Physics_applications/pulsar/inputs.corotating.3d.PM
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,8 +15,8 @@
#################################
max_step = 50000
amr.n_cell = 128 128 128
amr.max_grid_size = 64
amr.blocking_factor = 64
amr.max_grid_size = 128
amr.blocking_factor = 128
amr.max_level = 0
geometry.coord_sys = 0 # 0: Cartesian
geometry.is_periodic = 0 0 0 # Is periodic?
Expand All @@ -26,11 +26,12 @@ geometry.prob_hi = 180000 180000 180000
#################################
############ NUMERICS ###########
#################################
algo.maxwell_fdtd_solver = yee
#algo.maxwell_fdtd_solver = yee
warpx.verbose = 1
warpx.do_dive_cleaning = 0
warpx.use_filter = 1
warpx.cfl = .5
warpx.cfl = .95
warpx.do_pml = 1
my_constants.pi = 3.141592653589793
my_constants.dens = 5.544e6
my_constants.xc = 90000
Expand All @@ -40,20 +41,12 @@ my_constants.r_star = 12000
my_constants.omega = 6245.676
my_constants.c = 299792458.
my_constants.B_star = 8.0323e-6 # magnetic field of NS (T)
my_constants.dR = 1.4e3
my_constants.to = 2e-4
my_constants.dR = 1400
my_constants.to = 2.e-4
interpolation.nox = 3
interpolation.noy = 3
interpolation.noz = 3

#warpx.E_ext_grid_init_style = "parse_E_ext_grid_function"
#warpx.Ex_external_grid_function(x,y,z) = "((x-xc)*(x-xc)+(y-yc)*(y-yc) + (z-zc)*(z-zc))^0.5"
#warpx.Ey_external_grid_function(x,y,z) = "((x-xc)*(x-xc)+(y-yc)*(y-yc) + (z-zc)*(z-zc))^0.5"
#warpx.Ez_external_grid_function(x,y,z) = "((x-xc)*(x-xc)+(y-yc)*(y-yc) + (z-zc)*(z-zc))^0.5"





#################################
############ PLASMA #############
Expand All @@ -65,93 +58,83 @@ plasma_e.charge = -q_e
plasma_e.mass = m_e
plasma_e.injection_style = "NUniformPerCell"
plasma_e.profile = parse_density_function
plasma_e.density_function(x,y,z) = "( ((( (z-zc)*(z-zc) + (y-yc)*(y-yc) + (x-xc)*(x-xc) )^(0.5))<=(r_star+dR)) * ((( (z-zc)*(z-zc) + (y-yc)*(y-yc) + (x-xc)*(x-xc) )^(0.5))>=(r_star)) )*dens"
plasma_e.num_particles_per_cell_each_dim = 3 3 3
plasma_e.density_function(x,y,z) = "( ((( (z-zc)*(z-zc) + (y-yc)*(y-yc) + (x-xc)*(x-xc) )^(0.5))<=(r_star)) * ((( (z-zc)*(z-zc) + (y-yc)*(y-yc) + (x-xc)*(x-xc) )^(0.5))>=(r_star-dR)) )*dens"
plasma_e.num_particles_per_cell_each_dim = 4 4 4
plasma_e.momentum_distribution_type = parse_momentum_function

## Inject pairs rotating at omega(t) x r velocity
#plasma_e.momentum_function_ux(x,y,z) = "( (( (z-zc)*(z-zc) + (y-yc)*(y-yc) + (x-xc)*(x-xc) )^(0.5))>=r_star)*(-( (t<to)*(t/to)*omega + (t>=to)*omega )*(((x-xc)*(x-xc) + (y-yc)*(y-yc))^(0.5))/((1.0-( ( (t<to)*(t/to)*omega + (t>=to)*omega) *(((x-xc)*(x-xc) + (y-yc)*(y-yc))^(0.5))/c)^2)^(0.5)) * (y-yc)/(((x-xc)*(x-xc) + (y-yc)*(y-yc))^(0.5)))"
#plasma_e.momentum_function_uy(x,y,z) = "( (( (z-zc)*(z-zc) + (y-yc)*(y-yc) + (x-xc)*(x-xc) )^(0.5))>=r_star)*( ( (t<to)*(t/to)*omega + (t>=to)*omega ) *(((x-xc)*(x-xc) + (y-yc)*(y-yc))^(0.5))/((1.0-( ( (t<to)*t/to*omega + (t>=to)*omega )*(((x-xc)*(x-xc) + (y-yc)*(y-yc))^(0.5))/c)^2)^(0.5)) * (x-xc)/(((x-xc)*(x-xc) + (y-yc)*(y-yc))^(0.5)))"

## Inject pairs rotating at omega x r velocity
#plasma_e.momentum_function_ux(x,y,z) = "( (( (z-zc)*(z-zc) + (y-yc)*(y-yc) + (x-xc)*(x-xc) )^(0.5))<=r_star)*(-omega*(((x-xc)*(x-xc) + (y-yc)*(y-yc))^(0.5))/((1.0-(omega*(((x-xc)*(x-xc) + (y-yc)*(y-yc))^(0.5))/c)^2)^(0.5)) * (y-yc)/(((x-xc)*(x-xc) + (y-yc)*(y-yc))^(0.5)))"
#plasma_e.momentum_function_uy(x,y,z) = "( (( (z-zc)*(z-zc) + (y-yc)*(y-yc) + (x-xc)*(x-xc) )^(0.5))<=r_star)*(omega*(((x-xc)*(x-xc) + (y-yc)*(y-yc))^(0.5))/((1.0-(omega*(((x-xc)*(x-xc) + (y-yc)*(y-yc))^(0.5))/c)^2)^(0.5)) * (x-xc)/(((x-xc)*(x-xc) + (y-yc)*(y-yc))^(0.5)))"

## Inject stationary pairs
plasma_e.momentum_function_ux(x,y,z) = "0.0"
plasma_e.momentum_function_uy(x,y,z) = "0.0"

## uz is always 0 for the injection methods above
plasma_e.momentum_function_uz(x,y,z) = "0.0"

#plasma_e.momentum_distribution_type = "constant"
#plasma_e.ux = 0.0
#plasma_e.uy = 0.0
#plasma_e.uz = 0.0

plasma_e.do_continuous_injection = 0
plasma_e.density_min = 1E0

plasma_p.charge = q_e
plasma_p.mass = m_e
plasma_p.injection_style = "NUniformPerCell"
plasma_p.profile = parse_density_function
#plasma_p.density_function(x,y,z) = "( ((( (z-zc)*(z-zc) + (y-yc)*(y-yc) + (x-xc)*(x-xc) )^(0.5))>=r_star)* ((( (z-zc)*(z-zc) + (y-yc)*(y-yc) + (x-xc)*(x-xc) )^(0.5))<(r_star+dR)))*dens"
plasma_p.density_function(x,y,z) = "( ((( (z-zc)*(z-zc) + (y-yc)*(y-yc) + (x-xc)*(x-xc) )^(0.5))<=(r_star+dR)) * ((( (z-zc)*(z-zc) + (y-yc)*(y-yc) + (x-xc)*(x-xc) )^(0.5))>=(r_star)) )*dens"
plasma_p.num_particles_per_cell_each_dim = 3 3 3
plasma_p.density_function(x,y,z) = "( ((( (z-zc)*(z-zc) + (y-yc)*(y-yc) + (x-xc)*(x-xc) )^(0.5))<=(r_star)) * ((( (z-zc)*(z-zc) + (y-yc)*(y-yc) + (x-xc)*(x-xc) )^(0.5))>=(r_star-dR)) )*dens"
plasma_p.num_particles_per_cell_each_dim = 4 4 4
plasma_p.momentum_distribution_type = parse_momentum_function

## Inject pairs rotating at omega(t) x r velocity
#plasma_p.momentum_function_ux(x,y,z) = "( (( (z-zc)*(z-zc) + (y-yc)*(y-yc) + (x-xc)*(x-xc) )^(0.5))>=r_star)*(-( (t<to)*(t/to)*omega + (t>=to)*omega )*(((x-xc)*(x-xc) + (y-yc)*(y-yc))^(0.5))/((1.0-( ( (t<to)*(t/to)*omega + (t>=to)*omega) *(((x-xc)*(x-xc) + (y-yc)*(y-yc))^(0.5))/c)^2)^(0.5)) * (y-yc)/(((x-xc)*(x-xc) + (y-yc)*(y-yc))^(0.5)))"
#plasma_p.momentum_function_uy(x,y,z) = "( (( (z-zc)*(z-zc) + (y-yc)*(y-yc) + (x-xc)*(x-xc) )^(0.5))>=r_star)*( ( (t<to)*(t/to)*omega + (t>=to)*omega ) *(((x-xc)*(x-xc) + (y-yc)*(y-yc))^(0.5))/((1.0-( ( (t<to)*t/to*omega + (t>=to)*omega )*(((x-xc)*(x-xc) + (y-yc)*(y-yc))^(0.5))/c)^2)^(0.5)) * (x-xc)/(((x-xc)*(x-xc) + (y-yc)*(y-yc))^(0.5)))"

## Inject pairs rotating at omega x r velocity
#plasma_p.momentum_function_ux(x,y,z) = "( (( (z-zc)*(z-zc) + (y-yc)*(y-yc) + (x-xc)*(x-xc) )^(0.5))<=r_star)*(-omega*(((x-xc)*(x-xc) + (y-yc)*(y-yc))^(0.5))/((1.0-(omega*(((x-xc)*(x-xc) + (y-yc)*(y-yc))^(0.5))/c)^2)^(0.5)) * (y-yc)/(((x-xc)*(x-xc) + (y-yc)*(y-yc))^(0.5)))"
#plasma_p.momentum_function_uy(x,y,z) = "( (( (z-zc)*(z-zc) + (y-yc)*(y-yc) + (x-xc)*(x-xc) )^(0.5))<=r_star)*(omega*(((x-xc)*(x-xc) + (y-yc)*(y-yc))^(0.5))/((1.0-(omega*(((x-xc)*(x-xc) + (y-yc)*(y-yc))^(0.5))/c)^2)^(0.5)) * (x-xc)/(((x-xc)*(x-xc) + (y-yc)*(y-yc))^(0.5)))"

## Inject stationary pairs
plasma_p.momentum_function_ux(x,y,z) = "0.0"
plasma_p.momentum_function_uy(x,y,z) = "0.0"

## uz is always 0 for the injection methods above
plasma_p.momentum_function_uz(x,y,z) = "0.0"

#plasma_p.momentum_distribution_type = "constant"
#plasma_p.ux = 0.0
#plasma_p.uy = 0.0
#plasma_p.uz = 0.0

plasma_p.do_continuous_injection = 0
plasma_p.density_min = 1E0

diagnostics.diags_names = plt
diagnostics.diags_names = plt chk
plt.diag_type = Full
plt.period = 200
plt.period = 20
plt.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz part_per_cell rho divE

chk.format = checkpoint
chk.diag_type = Full
chk.period = 1000

#################################
######### PULSAR SETUP ##########
#################################
pulsar.pulsarType = "dead" # [dead/active]: sets particle injection type
pulsar.omega_star = 6245.676 # angular velocity of NS (rad/s)
<<<<<<< HEAD
pulsar.ramp_omega_time = 0.0 # time over which to ramp up NS angular velocity (s)
#pulsar.ramp_omega_time = 2e-4 # time over which to ramp up NS angular velocity (s)
=======
pulsar.ramp_omega_time = 0.e-5 # time over which to ramp up NS angular velocity (s)
>>>>>>> 2756eb86 (recent input file)
pulsar.ramp_omega_time = 0.0 # time over which to ramp up NS angular velocity (s)
# if ramp_omega_time < 0, omega = omega_star for any t
# consistency requires ramp_omega_time = my_constants.to
pulsar.center_star = 90000 90000 90000
pulsar.R_star = 12032 # radius of NS (m)
pulsar.B_star = 8.0323e-6 # magnetic field of NS (T)
pulsar.dR = 1075
pulsar.dR = 1400 # thickness on the surface of the pulsar
# consistency requires dR = my_constants.dR
pulsar.verbose = 0 # [0/1]: turn on verbosity for debugging print statements
pulsar.EB_external = 1 # [0/1]: to apply external E and B
pulsar.E_external_monopole = 0 # [0/1]
pulsar.EB_corotating_maxradius = 12032 # The radius where E-field changes from corotating
# inside the star to quadrapole outside.
# Default is R_star. (r<=cor_radius) i.e. the includes
# the r specified in the input
pulsar.max_ndens = 5.54e6 # max ndens == ndens used in initializing density
pulsar.Ninj_fraction = 0.1 # fraction of sigma injected
pulsar.Ninj_fraction = 0.2 # fraction of sigma injected
pulsar.ModifyParticleWeight = 1 # (0/1) the fractional injection is achieved
# by modifying particle weight if 1
# by modifying num_ppc if 0
pulsar.particle_inj_rmin = 10632 # default is Rstar-dR (consistent with density function above)
pulsar.particle_inj_rmax = 12032 # default is Rstar (consistent with density function)
pulsar.max_particle_absorption_radius = 12032 # Maximum radius within which particles are
# deleted every timestep.
# Default is Rstar
pulsar.rhoGJ_scale = 1e0 # scaling down of rho_GJ
pulsar.damp_EB_internal = 1 # Damp internal electric field
pulsar.damping_scale = 10.0 # Damping scale factor for internal electric field
pulsar.damp_EB_radius = 12032 # The radius within which EB should be damped.
# default is r_star (damping will include this radius)
pulsar.damping_scale = 1000.0 # Damping scale factor for internal electric field
pulsar.turnoffdeposition = 0 # (0/1) 0 for depositing (default)
# 1 for no deposition
pulsar.max_nodepos_radius = 0. # radius within which particle deposition for j/rho
# is turned off. (r<=max_radius)
pulsar.turnoff_plasmaEB_gather = 0 # (0/1) 0 is default. always gather
# 1 for no gather of self-consistent fields
pulsar.max_nogather_radius = 0. # radius within which self-consistent field gather
# is turned off
690 changes: 690 additions & 0 deletions Examples/Physics_applications/pulsar/pulsar_viz.ipynb
View file

Large diffs are not rendered by default.

31 changes: 11 additions & 20 deletions Source/Evolve/WarpXEvolve.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -156,33 +156,22 @@ WarpX::Evolve (int numsteps)
Bx = Bfield_fp[lev][0].get();
By = Bfield_fp[lev][1].get();
Bz = Bfield_fp[lev][2].get();
Gpu::ManagedVector<int> Ex_stag, Ey_stag, Ez_stag, Bx_stag, By_stag, Bz_stag;
Ex_stag.resize(3);
Ey_stag.resize(3);
Ez_stag.resize(3);
Bx_stag.resize(3);
By_stag.resize(3);
Bz_stag.resize(3);
amrex::IntVect ex_type = Ex->ixType().toIntVect();
amrex::IntVect ey_type = Ey->ixType().toIntVect();
amrex::IntVect ez_type = Ez->ixType().toIntVect();
amrex::IntVect bx_type = Bx->ixType().toIntVect();
amrex::IntVect by_type = By->ixType().toIntVect();
amrex::IntVect bz_type = Bz->ixType().toIntVect();
for (int idim = 0; idim < AMREX_SPACEDIM-1; ++idim) {
amrex::GpuArray<int, 3> Ex_stag, Ey_stag, Ez_stag, Bx_stag, By_stag, Bz_stag;
for (int idim = 0; idim < 3; ++idim)
{
Ex_stag[idim] = ex_type[idim];
Ey_stag[idim] = ey_type[idim];
Ez_stag[idim] = ez_type[idim];
Bx_stag[idim] = bx_type[idim];
By_stag[idim] = by_type[idim];
Bz_stag[idim] = bz_type[idim];
}
int const* const AMREX_RESTRICT Ex_stag_ptr = Ex_stag.data();
int const* const AMREX_RESTRICT Ey_stag_ptr = Ey_stag.data();
int const* const AMREX_RESTRICT Ez_stag_ptr = Ez_stag.data();
int const* const AMREX_RESTRICT Bx_stag_ptr = Bx_stag.data();
int const* const AMREX_RESTRICT By_stag_ptr = By_stag.data();
int const* const AMREX_RESTRICT Bz_stag_ptr = Bz_stag.data();
auto geom = Geom(lev).data();
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
Expand All @@ -200,15 +189,15 @@ WarpX::Evolve (int numsteps)
amrex::ParallelFor(tex, tey, tez,
[=] AMREX_GPU_DEVICE (int i, int j, int k)
{
PulsarParm::DampEField(i, j, k, geom, Exfab, Ex_stag_ptr);
PulsarParm::DampField(i, j, k, geom, Exfab, Ex_stag);
},
[=] AMREX_GPU_DEVICE (int i, int j, int k)
{
PulsarParm::DampEField(i, j, k, geom, Eyfab, Ey_stag_ptr);
PulsarParm::DampField(i, j, k, geom, Eyfab, Ey_stag);
},
[=] AMREX_GPU_DEVICE (int i, int j, int k)
{
PulsarParm::DampEField(i, j, k, geom, Ezfab, Ez_stag_ptr);
PulsarParm::DampField(i, j, k, geom, Ezfab, Ez_stag);
});
}
for ( MFIter mfi(*Bx, TilingIfNotGPU()); mfi.isValid(); ++mfi )
Expand All @@ -224,19 +213,21 @@ WarpX::Evolve (int numsteps)
amrex::ParallelFor(tex, tey, tez,
[=] AMREX_GPU_DEVICE (int i, int j, int k)
{
PulsarParm::DampEField(i, j, k, geom, Bxfab, Bx_stag_ptr);
PulsarParm::DampField(i, j, k, geom, Bxfab, Bx_stag);
},
[=] AMREX_GPU_DEVICE (int i, int j, int k)
{
PulsarParm::DampEField(i, j, k, geom, Byfab, By_stag_ptr);
PulsarParm::DampField(i, j, k, geom, Byfab, By_stag);
},
[=] AMREX_GPU_DEVICE (int i, int j, int k)
{
PulsarParm::DampEField(i, j, k, geom, Bzfab, Bz_stag_ptr);
PulsarParm::DampField(i, j, k, geom, Bzfab, Bz_stag);
});
}
}
}
FillBoundaryE(guard_cells.ng_FieldSolver, IntVect::TheZeroVector());
FillBoundaryB(guard_cells.ng_FieldSolver, IntVect::TheZeroVector());
#endif

if (warpx_py_beforeEsolve) warpx_py_beforeEsolve();
Expand Down
76 changes: 69 additions & 7 deletions Source/Initialization/InjectorPosition.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,10 @@
#include <AMReX_Dim3.H>
#include <AMReX_Utility.H>

#ifdef PULSAR
#include "Particles/PulsarParameters.H"
#endif

// struct whose getPositionUnitBox returns x, y and z for a particle with
// random distribution inside a unit cell.
struct InjectorPositionRandom
Expand Down Expand Up @@ -41,14 +45,31 @@ struct InjectorPositionRegular
amrex::RandomEngine const&) const noexcept
{
using namespace amrex;

#ifdef PULSAR
#if (defined WARPX_DIM_3D) || (defined WARPX_DIM_RZ)
int nx, ny, nz;
if (PulsarParm::ModifyParticleWtAtInjection == 1) {
nx = ref_fac*ppc.x;
ny = ref_fac*ppc.y;
nz = ref_fac*ppc.z;
} else if (PulsarParm::ModifyParticleWtAtInjection == 0) {
amrex::Real ninj_per_dim = std::cbrt(PulsarParm::Ninj_fraction);
nx = ref_fac*ppc.x*ninj_per_dim;
ny = ref_fac*ppc.y*ninj_per_dim;
nz = ref_fac*ppc.z*ninj_per_dim;
}
#else
amrex::Abort("pulsar sim does not support 2d yet");
#endif // endif for 3D
#else // else ifndef PULSAR
int const nx = ref_fac*ppc.x;
int const ny = ref_fac*ppc.y;
#if (defined WARPX_DIM_3D) || (defined WARPX_DIM_RZ)
int const nz = ref_fac*ppc.z;
#else
int const nz = 1;
#endif
#endif // 3D ifdef
#endif // PULSAR ifdef
int const ix_part = i_part / (ny*nz); // written this way backward compatibility
int const iz_part = (i_part-ix_part*(ny*nz)) / ny;
int const iy_part = (i_part-ix_part*(ny*nz)) - ny*iz_part;
Expand All @@ -58,6 +79,17 @@ struct InjectorPositionRegular
(0.5_rt + iz_part) / nz
};
}

#ifdef PULSAR
AMREX_GPU_HOST_DEVICE
amrex::XDim3
getppcInEachDim () const noexcept
{
using namespace amrex;
return XDim3{ppc.x, ppc.y, ppc.z};
}
#endif

private:
amrex::Dim3 ppc;
};
Expand Down Expand Up @@ -123,6 +155,25 @@ struct InjectorPosition
};
}

#ifdef PULSAR
AMREX_GPU_HOST_DEVICE
amrex::XDim3
getppcInEachDim () const noexcept
{
switch (type)
{
case Type::regular:
{
return object.regular.getppcInEachDim();
}
default:
{
return object.regular.getppcInEachDim();
}
};
}
#endif

// bool: whether position specified is within bounds.
AMREX_GPU_HOST_DEVICE
bool
Expand All @@ -132,15 +183,26 @@ struct InjectorPosition
y < ymax and y >= ymin and
z < zmax and z >= zmin);
}

// bool: whether position specified is withing pulsar injection region

#ifdef PULSAR
// bool: whether cell-center position specified is within pulsar injection region
// So that these cells can inject particles
// a buffer factor of 2*dR (2*dx) is used to ensure cells partially cut by sphere
// inject particles
AMREX_GPU_HOST_DEVICE
bool
insidePulsarBoundsCC (amrex::Real r, amrex::Real r_min, amrex::Real r_max, amrex::Real dR_factor) const noexcept
{
return (r<=(r_max + dR_factor) and r>=(r_min-dR_factor));
}

// bool: whether particle position specified is within pulsar injection region
// So that these cells can inject particles
AMREX_GPU_HOST_DEVICE
bool
insidePulsarBounds (amrex::Real r, amrex::Real R_star, amrex::Real dR_star) const noexcept
insidePulsarBounds (amrex::Real r, amrex::Real r_min, amrex::Real r_max) const noexcept
{
//return (r<=R_star and r>=(R_star-dR_star));
return (r<=(R_star+dR_star) and r>=(R_star-dR_star));
return (r<=(r_max) and r>=(r_min));
}
#endif

Expand Down
Loading

0 comments on commit 035908f

Please sign in to comment.