diff --git a/Examples/Physics_applications/pulsar/inputs.corotating.3d.PM b/Examples/Physics_applications/pulsar/inputs.corotating.3d.PM index d375825e8ea..1ae3af9cd2b 100644 --- a/Examples/Physics_applications/pulsar/inputs.corotating.3d.PM +++ b/Examples/Physics_applications/pulsar/inputs.corotating.3d.PM @@ -13,16 +13,18 @@ ################################# ####### GENERAL PARAMETERS ###### ################################# -max_step = 5 +max_step = 1 amr.n_cell = 128 128 128 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? -geometry.prob_lo = 0.0 0.0 0.0 +geometry.prob_lo = 0 0 0 geometry.prob_hi = 180000 180000 180000 +boundary.field_lo = pml pml pml +boundary.field_hi = pml pml pml + ################################# ############ NUMERICS ########### ################################# @@ -31,17 +33,16 @@ warpx.verbose = 1 warpx.do_dive_cleaning = 0 warpx.use_filter = 1 warpx.cfl = .95 -warpx.do_pml = 1 my_constants.pi = 3.141592653589793 my_constants.dens = 5.544e6 my_constants.xc = 90000 my_constants.yc = 90000 my_constants.zc = 90000 -my_constants.r_star = 12000 +my_constants.r_star = 12032 my_constants.omega = 6245.676 my_constants.c = 299792458. my_constants.B_star = 8.0323e-6 # magnetic field of NS (T) -my_constants.dR = 1400 +my_constants.dR = 1406.25 my_constants.to = 2.e-4 algo.particle_shape = 3 @@ -49,8 +50,8 @@ algo.particle_shape = 3 ################################# ############ PLASMA ############# ################################# -particles.nspecies = 2 -particles.species_names = plasma_e plasma_p +#particles.nspecies = 2 +#particles.species_names = plasma_e plasma_p plasma_e.charge = -q_e plasma_e.mass = m_e @@ -102,11 +103,11 @@ pulsar.ramp_omega_time = 0.0 # time over which to ramp up NS angular 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 = 1400 # thickness on the surface of the pulsar +pulsar.dR = 1406.25 # 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 = 0 # [0/1]: to apply external E and B -pulsar.E_external_monopole = 0 # [0/1] +pulsar.E_external_monopole = 1 # [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 @@ -135,11 +136,18 @@ pulsar.turnoff_plasmaEB_gather = 0 # (0/1) 0 is default. always gather pulsar.max_nogather_radius = 0. # radius within which self-consistent field gather pulsar.init_dipoleBfield = 1 # default is 1 pulsar.init_corotatingEfield = 1 # default -pulsar.corotatingE_maxradius = 12032 # default is Rstar +pulsar.corotatingE_maxradius = 16250.75 # default is Rstar pulsar.enforceDipoleB_maxradius = 12032 # default is Rstar pulsar.enforceCorotatingE = 1 # default 1 pulsar.enforceDipoleB = 1 # default 1 pulsar.singleParticleTest = 0 # default 0, if 1, then a particle pair will be introduced. # Single particle position/momentum will need to be specified -pulsar.DampBDipoleInRing = 1 +pulsar.DampBDipoleInRing = 0 pulsar.Bdamping_scale = 10000 +pulsar.LimitDipoleBInit = 0 # if 1 then Bdipole is iniialization only in a smaller region +pulsar.DipoleB_init_maxradius = 12032 # read if LimitBdipole is 1 +pulsar.use_theoreticalEB = 0 +pulsar.theory_max_rstar = 100000 +pulsar.EnforceTheoreticalEBInGrid = 1 +pulsar.AddExternalMonopoleOnly = 1 + diff --git a/Source/Diagnostics/ComputeDiagFunctors/EdotBFunctor.H b/Source/Diagnostics/ComputeDiagFunctors/EdotBFunctor.H index 5c445ae5f01..a12a8dd9da3 100644 --- a/Source/Diagnostics/ComputeDiagFunctors/EdotBFunctor.H +++ b/Source/Diagnostics/ComputeDiagFunctors/EdotBFunctor.H @@ -16,7 +16,7 @@ public: const int lev, const amrex::IntVect crse_ratio, int ncomp = 1); - virtual void operator() (amrex::MultiFab& mf_dst, int dcomp, const int /*i_buffer=0*/) const override; + virtual void operator() (amrex::MultiFab& mf_dst, int dcomp, const int /*i_buffer=0*/) const override; void ComputeEdotB(amrex::MultiFab& mf_dst, int dcomp) const; private: amrex::MultiFab const * const m_Ex_src = nullptr; diff --git a/Source/Diagnostics/ComputeDiagFunctors/EdotBFunctor.cpp b/Source/Diagnostics/ComputeDiagFunctors/EdotBFunctor.cpp index c91d129dae0..e4994d6e2fb 100644 --- a/Source/Diagnostics/ComputeDiagFunctors/EdotBFunctor.cpp +++ b/Source/Diagnostics/ComputeDiagFunctors/EdotBFunctor.cpp @@ -21,10 +21,6 @@ void EdotBFunctor::operator ()(amrex::MultiFab& mf_dst, int dcomp, const int /*i_buffer=0*/) const { using namespace amrex; - auto & warpx = WarpX::GetInstance(); - const auto dx = warpx.Geom(m_lev).CellSizeArray(); - const auto problo = warpx.Geom(m_lev).ProbLoArray(); - const auto probhi = warpx.Geom(m_lev).ProbHiArray(); const amrex::IntVect stag_dst = mf_dst.ixType().toIntVect(); // convert boxarray of source MultiFab to staggering of dst Multifab @@ -38,14 +34,14 @@ EdotBFunctor::operator ()(amrex::MultiFab& mf_dst, int dcomp, const int /*i_buff ComputeEdotB(mf_dst, dcomp); } else { const int ncomp = 1; - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( m_Ex_src->DistributionMap() == m_Ey_src->DistributionMap() and m_Ey_src->DistributionMap() == m_Ez_src->DistributionMap(), + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( m_Ex_src->DistributionMap() == m_Ey_src->DistributionMap() and m_Ey_src->DistributionMap() == m_Ez_src->DistributionMap(), " all sources must have the same Distribution map"); - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( m_Bx_src->DistributionMap() == m_By_src->DistributionMap() and m_By_src-> DistributionMap() == m_Bz_src->DistributionMap(), + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( m_Bx_src->DistributionMap() == m_By_src->DistributionMap() and m_By_src-> DistributionMap() == m_Bz_src->DistributionMap(), " all sources must have the same Distribution map"); amrex::MultiFab mf_tmp( ba_tmp, m_Ex_src->DistributionMap(), ncomp, 0); const int dcomp_tmp = 0; ComputeEdotB(mf_tmp, dcomp_tmp); - mf_dst.copy( mf_tmp, 0, dcomp, ncomp); + mf_dst.ParallelCopy( mf_tmp, 0, dcomp, ncomp); } } @@ -53,10 +49,6 @@ void EdotBFunctor::ComputeEdotB(amrex::MultiFab& mf_dst, int dcomp) const { using namespace amrex; - auto & warpx = WarpX::GetInstance(); - const auto dx = warpx.Geom(m_lev).CellSizeArray(); - const auto problo = warpx.Geom(m_lev).ProbLoArray(); - const auto probhi = warpx.Geom(m_lev).ProbHiArray(); const amrex::IntVect stag_Exsrc = m_Ex_src->ixType().toIntVect(); const amrex::IntVect stag_Eysrc = m_Ey_src->ixType().toIntVect(); const amrex::IntVect stag_Ezsrc = m_Ez_src->ixType().toIntVect(); diff --git a/Source/Diagnostics/ComputeDiagFunctors/PoyntingVectorFunctor.H b/Source/Diagnostics/ComputeDiagFunctors/PoyntingVectorFunctor.H index ee1a4438feb..acc581ba1e7 100644 --- a/Source/Diagnostics/ComputeDiagFunctors/PoyntingVectorFunctor.H +++ b/Source/Diagnostics/ComputeDiagFunctors/PoyntingVectorFunctor.H @@ -17,7 +17,7 @@ public: const amrex::IntVect crse_ratio, int vectorcomp, int ncomp = 1); - virtual void operator() (amrex::MultiFab& mf_dst, int dcomp, const int /*i_buffer=0*/) const override; + virtual void operator() (amrex::MultiFab& mf_dst, int dcomp, const int /*i_buffer=0*/) const override; void ComputePoyntingVector(amrex::MultiFab& mf_dst, int dcomp) const; private: diff --git a/Source/Diagnostics/ComputeDiagFunctors/PoyntingVectorFunctor.cpp b/Source/Diagnostics/ComputeDiagFunctors/PoyntingVectorFunctor.cpp index 603c6736bb7..09aca59eff1 100644 --- a/Source/Diagnostics/ComputeDiagFunctors/PoyntingVectorFunctor.cpp +++ b/Source/Diagnostics/ComputeDiagFunctors/PoyntingVectorFunctor.cpp @@ -25,10 +25,6 @@ PoyntingVectorFunctor::operator ()(amrex::MultiFab& mf_dst, int dcomp, const int /*i_buffer=0*/) const { using namespace amrex; - auto & warpx = WarpX::GetInstance(); - const auto dx = warpx.Geom(m_lev).CellSizeArray(); - const auto problo = warpx.Geom(m_lev).ProbLoArray(); - const auto probhi = warpx.Geom(m_lev).ProbHiArray(); const amrex::IntVect stag_dst = mf_dst.ixType().toIntVect(); // convert boxarray of source MultiFab to staggering of dst Multifab @@ -42,14 +38,14 @@ PoyntingVectorFunctor::operator ()(amrex::MultiFab& mf_dst, ComputePoyntingVector(mf_dst, dcomp); } else { const int ncomp = 1; - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( m_Ex_src->DistributionMap() == m_Ey_src->DistributionMap() and m_Ey_src->DistributionMap() == m_Ez_src->DistributionMap(), + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( m_Ex_src->DistributionMap() == m_Ey_src->DistributionMap() and m_Ey_src->DistributionMap() == m_Ez_src->DistributionMap(), " all sources must have the same Distribution map"); - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( m_Bx_src->DistributionMap() == m_By_src->DistributionMap() and m_By_src->DistributionMap() == m_Bz_src->DistributionMap(), + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( m_Bx_src->DistributionMap() == m_By_src->DistributionMap() and m_By_src->DistributionMap() == m_Bz_src->DistributionMap(), " all sources must have the same Distribution map"); amrex::MultiFab mf_tmp( ba_tmp, m_Ex_src->DistributionMap(), ncomp, 0 ); const int dcomp_tmp = 0; ComputePoyntingVector(mf_tmp, dcomp_tmp); - mf_dst.copy( mf_tmp, 0, dcomp, ncomp); + mf_dst.ParallelCopy( mf_tmp, 0, dcomp, ncomp); } } @@ -58,10 +54,6 @@ void PoyntingVectorFunctor::ComputePoyntingVector(amrex::MultiFab& mf_dst, int dcomp) const { using namespace amrex; - auto & warpx = WarpX::GetInstance(); - const auto dx = warpx.Geom(m_lev).CellSizeArray(); - const auto problo = warpx.Geom(m_lev).ProbLoArray(); - const auto probhi = warpx.Geom(m_lev).ProbHiArray(); const amrex::IntVect stag_Exsrc = m_Ex_src->ixType().toIntVect(); const amrex::IntVect stag_Eysrc = m_Ey_src->ixType().toIntVect(); const amrex::IntVect stag_Ezsrc = m_Ez_src->ixType().toIntVect(); @@ -90,7 +82,6 @@ PoyntingVectorFunctor::ComputePoyntingVector(amrex::MultiFab& mf_dst, int dcomp) cr[i] = m_crse_ratio[i]; } const int vectorcomp = m_vectorcomp; - amrex::Real mu0_inv = 1.0_rt/PhysConst::mu0; #ifdef AMREX_USE_OMP #pragma omp parallel if (Gpu::notInLaunchRegion()) diff --git a/Source/Diagnostics/ComputeDiagFunctors/SphericalComponentFunctor.H b/Source/Diagnostics/ComputeDiagFunctors/SphericalComponentFunctor.H index 01ebdd16dbe..c5c2a3df740 100644 --- a/Source/Diagnostics/ComputeDiagFunctors/SphericalComponentFunctor.H +++ b/Source/Diagnostics/ComputeDiagFunctors/SphericalComponentFunctor.H @@ -8,7 +8,7 @@ SphericalComponentFunctor final : public ComputeDiagFunctor { public: SphericalComponentFunctor ( const amrex::MultiFab * const mfx_src, - const amrex::MultiFab * const mfy_src, + const amrex::MultiFab * const mfy_src, const amrex::MultiFab * const mfz_src, const int lev, const amrex::IntVect crse_ratio, @@ -26,7 +26,7 @@ private: amrex::MultiFab const * const m_mfz_src = nullptr; int m_lev; int m_sphericalcomp; - int m_Efield; + int m_Efield; }; #endif diff --git a/Source/Diagnostics/ComputeDiagFunctors/SphericalComponentFunctor.cpp b/Source/Diagnostics/ComputeDiagFunctors/SphericalComponentFunctor.cpp index 23d87e91f52..fda23991952 100644 --- a/Source/Diagnostics/ComputeDiagFunctors/SphericalComponentFunctor.cpp +++ b/Source/Diagnostics/ComputeDiagFunctors/SphericalComponentFunctor.cpp @@ -12,7 +12,7 @@ #include SphericalComponentFunctor::SphericalComponentFunctor (amrex::MultiFab const * mfx_src, - amrex::MultiFab const * mfy_src, + amrex::MultiFab const * mfy_src, amrex::MultiFab const * mfz_src, int lev, amrex::IntVect crse_ratio, @@ -25,14 +25,10 @@ SphericalComponentFunctor::SphericalComponentFunctor (amrex::MultiFab const * mf {} -void +void SphericalComponentFunctor::operator ()(amrex::MultiFab& mf_dst, int dcomp, const int /*i_buffer=0*/) const { using namespace amrex; - auto & warpx = WarpX::GetInstance(); - const auto dx = warpx.Geom(m_lev).CellSizeArray(); - const auto problo = warpx.Geom(m_lev).ProbLoArray(); - const auto probhi = warpx.Geom(m_lev).ProbHiArray(); const amrex::IntVect stag_dst = mf_dst.ixType().toIntVect(); // convert boxarray of source MultiFab to staggering of dst Multifab @@ -46,12 +42,12 @@ SphericalComponentFunctor::operator ()(amrex::MultiFab& mf_dst, int dcomp, const ComputeSphericalFieldComponent(mf_dst, dcomp); } else { const int ncomp = 1; - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( m_mfx_src->DistributionMap() == m_mfy_src->DistributionMap() and m_mfy_src->DistributionMap() == m_mfz_src->DistributionMap(), + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( m_mfx_src->DistributionMap() == m_mfy_src->DistributionMap() and m_mfy_src->DistributionMap() == m_mfz_src->DistributionMap(), " all sources must have the same Distribution map"); amrex::MultiFab mf_tmp( ba_tmp, m_mfx_src->DistributionMap(), ncomp, 0); const int dcomp_tmp = 0; ComputeSphericalFieldComponent(mf_tmp, dcomp_tmp); - mf_dst.copy( mf_tmp, 0, dcomp, ncomp); + mf_dst.ParallelCopy( mf_tmp, 0, dcomp, ncomp); } } @@ -59,10 +55,10 @@ void SphericalComponentFunctor::ComputeSphericalFieldComponent( amrex::MultiFab& mf_dst, int dcomp) const { using namespace amrex; +#ifdef PULSAR auto & warpx = WarpX::GetInstance(); const auto dx = warpx.Geom(m_lev).CellSizeArray(); const auto problo = warpx.Geom(m_lev).ProbLoArray(); - const auto probhi = warpx.Geom(m_lev).ProbHiArray(); const amrex::IntVect stag_xsrc = m_mfx_src->ixType().toIntVect(); const amrex::IntVect stag_ysrc = m_mfy_src->ixType().toIntVect(); const amrex::IntVect stag_zsrc = m_mfz_src->ixType().toIntVect(); @@ -73,18 +69,16 @@ SphericalComponentFunctor::ComputeSphericalFieldComponent( amrex::MultiFab& mf_d amrex::GpuArray sfz; // staggering of source zfield amrex::GpuArray s_dst; amrex::GpuArray cr; - + amrex::GpuArray center_star_arr; for (int i=0; icompile(); + int integral_type = m_integral_type; // MFIter loop to interpolate fields to cell center and perform reduction #ifdef AMREX_USE_OMP @@ -190,11 +192,20 @@ public: amrex::ParallelDescriptor::ReduceRealSum(reduce_value); // If reduction operation is a sum, multiply the value by the cell volume so that the // result is the integral of the function over the simulation domain. + if (integral_type == 0) { #if (AMREX_SPACEDIM==2) reduce_value *= dx[0]*dx[1]; #else reduce_value *= dx[0]*dx[1]*dx[2]; #endif + } else { + // works only when dx = dy = dz +#if (AMREX_SPACEDIM==2) + reduce_value *= dx[0]; +#else + reduce_value *= dx[0]*dx[1]; +#endif + } } // Fill output array diff --git a/Source/Diagnostics/ReducedDiags/FieldReduction.cpp b/Source/Diagnostics/ReducedDiags/FieldReduction.cpp index 0ebc89a7d06..7695c5ca4cb 100644 --- a/Source/Diagnostics/ReducedDiags/FieldReduction.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldReduction.cpp @@ -61,6 +61,9 @@ FieldReduction::FieldReduction (std::string rd_name) std::string reduction_type_string; pp_rd_name.get("reduction_type", reduction_type_string); m_reduction_type = GetAlgorithmInteger (pp_rd_name, "reduction_type"); + if (m_reduction_type == 2) { + m_integral_type = GetAlgorithmInteger (pp_rd_name, "integration_type"); + } if (amrex::ParallelDescriptor::IOProcessor()) { diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 7076245c4f9..56ec8496b18 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -36,6 +36,7 @@ #include #include +#include #include #include #include @@ -119,23 +120,25 @@ WarpX::Evolve (int numsteps) } #ifdef PULSAR - if (PulsarParm::singleParticleTest == 1) { - if (PulsarParm::continuous_injection == 0) { - if (PulsarParm::injection_time - dt[0] <= cur_time && - cur_time <= PulsarParm::injection_time) { - amrex::Print() << " injecting particles \n"; + if (Pulsar::m_singleParticleTest == 1) { + if (Pulsar::m_continuous_injection == 0) { + if (Pulsar::m_injection_time <= cur_time && + cur_time <= Pulsar::m_injection_time + 0.5*dt[0]) { + amrex::Print() << " injecting particles \n"; mypc->PulsarParticleInjection(); mypc->Redistribute(); } } else { - if ( cur_time > PulsarParm::injection_time) { + if ( cur_time > Pulsar::m_injection_time) { mypc->PulsarParticleInjection(); mypc->Redistribute(); } } } else { - mypc->PulsarParticleInjection(); - mypc->Redistribute(); + if ( cur_time >= Pulsar::m_injection_time) { + mypc->PulsarParticleInjection(); + mypc->Redistribute(); + } } #endif // At the beginning, we have B^{n} and E^{n}. @@ -240,7 +243,7 @@ WarpX::Evolve (int numsteps) } #ifdef PULSAR - if (PulsarParm::damp_EB_internal) { + if (Pulsar::m_do_damp_EB_internal == 1 ) { MultiFab *Ex, *Ey, *Ez; MultiFab *Bx, *By, *Bz; for (int lev = 0; lev <= finest_level; ++lev) { @@ -257,6 +260,8 @@ WarpX::Evolve (int numsteps) amrex::IntVect by_type = By->ixType().toIntVect(); amrex::IntVect bz_type = Bz->ixType().toIntVect(); amrex::GpuArray Ex_stag, Ey_stag, Ez_stag, Bx_stag, By_stag, Bz_stag; + amrex::GpuArray center_star_data; + for (int idim = 0; idim < 3; ++idim) { Ex_stag[idim] = ex_type[idim]; @@ -265,8 +270,14 @@ WarpX::Evolve (int numsteps) Bx_stag[idim] = bx_type[idim]; By_stag[idim] = by_type[idim]; Bz_stag[idim] = bz_type[idim]; + center_star_data[idim] = Pulsar::m_center_star[idim]; } - auto geom = Geom(lev).data(); + amrex::Real max_EBdamping_radius_data = Pulsar::m_max_EBdamping_radius; + amrex::Real damping_scale_data = Pulsar::m_field_damping_scale; + amrex::Real Rstar_data = Pulsar::m_R_star; + const auto domain_xlo = WarpX::GetInstance().Geom(lev).ProbLoArray(); + const auto domain_xhi = WarpX::GetInstance().Geom(lev).ProbHiArray(); + const auto domain_dx = WarpX::GetInstance().Geom(lev).CellSizeArray(); #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif @@ -283,15 +294,24 @@ WarpX::Evolve (int numsteps) amrex::ParallelFor(tex, tey, tez, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - PulsarParm::DampField(i, j, k, geom, Exfab, Ex_stag); + Pulsar::DampField(i, j, k, Exfab, Ex_stag, center_star_data, + domain_xlo, domain_dx, + max_EBdamping_radius_data, + damping_scale_data, Rstar_data); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - PulsarParm::DampField(i, j, k, geom, Eyfab, Ey_stag); + Pulsar::DampField(i, j, k, Eyfab, Ey_stag, center_star_data, + domain_xlo, domain_dx, + max_EBdamping_radius_data, + damping_scale_data, Rstar_data); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - PulsarParm::DampField(i, j, k, geom, Ezfab, Ez_stag); + Pulsar::DampField(i, j, k, Ezfab, Ez_stag, center_star_data, + domain_xlo, domain_dx, + max_EBdamping_radius_data, + damping_scale_data, Rstar_data); }); } for ( MFIter mfi(*Bx, TilingIfNotGPU()); mfi.isValid(); ++mfi ) @@ -307,15 +327,24 @@ WarpX::Evolve (int numsteps) amrex::ParallelFor(tex, tey, tez, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - PulsarParm::DampField(i, j, k, geom, Bxfab, Bx_stag); + Pulsar::DampField(i, j, k, Bxfab, Bx_stag, center_star_data, + domain_xlo, domain_dx, + max_EBdamping_radius_data, + damping_scale_data, Rstar_data); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - PulsarParm::DampField(i, j, k, geom, Byfab, By_stag); + Pulsar::DampField(i, j, k, Byfab, By_stag, center_star_data, + domain_xlo, domain_dx, + max_EBdamping_radius_data, + damping_scale_data, Rstar_data); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - PulsarParm::DampField(i, j, k, geom, Bzfab, Bz_stag); + Pulsar::DampField(i, j, k, Bzfab, Bz_stag, center_star_data, + domain_xlo, domain_dx, + max_EBdamping_radius_data, + damping_scale_data, Rstar_data); }); } } diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 7425b670805..3579943ff9e 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -444,14 +444,14 @@ WarpX::PushPSATD () } #ifdef PULSAR amrex::Real a_dt = 0._rt; - if (PulsarParm::enforceDipoleB == 1) { - PulsarParm::ApplyDipoleBfield_BC( Bfield_fp[lev], lev, a_dt); + if (Pulsar::m_enforceDipoleB == 1) { + Pulsar::ApplyDipoleBfield_BC( Bfield_fp[lev], lev, a_dt); } - if (PulsarParm::enforceCorotatingE == 1) { - PulsarParm::ApplyCorotatingEfield_BC( Efield_fp[lev], lev, a_dt); + if (Pulsar::m_enforceCorotatingE == 1) { + Pulsar::ApplyCorotatingEfield_BC( Efield_fp[lev], lev, a_dt); } -#endif +#endif #endif } @@ -501,10 +501,10 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt, DtType a_dt_typ } #ifdef PULSAR - if (PulsarParm::enforceDipoleB == 1) { - PulsarParm::ApplyDipoleBfield_BC( Bfield_fp[lev], lev, a_dt); + if (Pulsar::m_enforceDipoleB == 1) { + Pulsar::ApplyDipoleBfield_BC( Bfield_fp[lev], lev, a_dt); } -#endif +#endif ApplyBfieldBoundary(lev, patch_type, a_dt_type); @@ -564,10 +564,10 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) } } #ifdef PULSAR - if (PulsarParm::enforceCorotatingE == 1) { - PulsarParm::ApplyCorotatingEfield_BC( Efield_fp[lev], lev, a_dt); + if (Pulsar::m_enforceCorotatingE == 1) { + Pulsar::ApplyCorotatingEfield_BC( Efield_fp[lev], lev, a_dt); } -#endif +#endif ApplyEfieldBoundary(lev, patch_type); } @@ -711,10 +711,10 @@ WarpX::MacroscopicEvolveE (int lev, PatchType patch_type, amrex::Real a_dt) { } } #ifdef PULSAR - if (PulsarParm::enforceCorotatingE == 1) { - PulsarParm::ApplyCorotatingEfield_BC( Efield_fp[lev], lev, a_dt); + if (Pulsar::m_enforceCorotatingE == 1) { + Pulsar::ApplyCorotatingEfield_BC( Efield_fp[lev], lev, a_dt); } -#endif +#endif ApplyEfieldBoundary(lev, patch_type); } diff --git a/Source/Initialization/InjectorPosition.H b/Source/Initialization/InjectorPosition.H index d574c365a19..5749a304d61 100644 --- a/Source/Initialization/InjectorPosition.H +++ b/Source/Initialization/InjectorPosition.H @@ -25,7 +25,12 @@ struct InjectorPositionRandom AMREX_GPU_HOST_DEVICE amrex::XDim3 getPositionUnitBox (int /*i_part*/, int /*ref_fac*/, - amrex::RandomEngine const& engine) const noexcept + amrex::RandomEngine const& engine +#ifdef PULSAR + , int ModifyParticleWtAtInjection = 0, + amrex::Real Ninj_fraction = 1. +#endif + ) const noexcept { return amrex::XDim3{amrex::Random(engine), amrex::Random(engine), amrex::Random(engine)}; } @@ -40,7 +45,12 @@ struct InjectorPositionRandomPlane AMREX_GPU_HOST_DEVICE amrex::XDim3 getPositionUnitBox (int /*i_part*/, int /*ref_fac*/, - amrex::RandomEngine const& engine) const noexcept + amrex::RandomEngine const& engine +#ifdef PULSAR + , int ModifyParticleWtAtInjection = 0, + amrex::Real Ninj_fraction = 1. +#endif + ) const noexcept { using namespace amrex::literals; if (dir == 0) return amrex::XDim3{0._rt, amrex::Random(engine), amrex::Random(engine)}; @@ -68,18 +78,23 @@ struct InjectorPositionRegular AMREX_GPU_HOST_DEVICE amrex::XDim3 getPositionUnitBox (int const i_part, int const ref_fac, - amrex::RandomEngine const&) const noexcept + amrex::RandomEngine const& +#ifdef PULSAR + , int ModifyParticleWtAtInjection = 0, + amrex::Real Ninj_fraction = 1. +#endif + ) const noexcept { using namespace amrex; #ifdef PULSAR #if (defined WARPX_DIM_3D) || (defined WARPX_DIM_RZ) int nx, ny, nz; - if (PulsarParm::ModifyParticleWtAtInjection == 1) { + if (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); + } else if (ModifyParticleWtAtInjection == 0) { + amrex::Real ninj_per_dim = std::cbrt(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; @@ -179,21 +194,41 @@ struct InjectorPosition AMREX_GPU_HOST_DEVICE amrex::XDim3 getPositionUnitBox (int const i_part, int const ref_fac, - amrex::RandomEngine const& engine) const noexcept + amrex::RandomEngine const& engine +#ifdef PULSAR + , int ModifyParticleWtAtInjection = 0, + amrex::Real Ninj_fraction = 1. +#endif + ) const noexcept { switch (type) { case Type::regular: { - return object.regular.getPositionUnitBox(i_part, ref_fac, engine); + return object.regular.getPositionUnitBox(i_part, ref_fac, engine +#ifdef PULSAR + , ModifyParticleWtAtInjection, + Ninj_fraction +#endif + ); } case Type::randomplane: { - return object.randomplane.getPositionUnitBox(i_part, ref_fac, engine); + return object.randomplane.getPositionUnitBox(i_part, ref_fac, engine +#ifdef PULSAR + , ModifyParticleWtAtInjection, + Ninj_fraction +#endif + ); } default: { - return object.random.getPositionUnitBox(i_part, ref_fac, engine); + return object.random.getPositionUnitBox(i_part, ref_fac, engine +#ifdef PULSAR + , ModifyParticleWtAtInjection, + Ninj_fraction +#endif + ); } }; } diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 730250b6afa..8c040695838 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -574,35 +574,35 @@ WarpX::InitLevelData (int lev, Real /*time*/) } } #ifdef PULSAR - if (PulsarParm::InitializeGrid_with_Pulsar_Bfield == 1) { + if (Pulsar::m_do_InitializeGrid_with_Pulsar_Bfield == 1) { bool Init_Bfield = true; - PulsarParm::InitializeExternalPulsarFieldsOnGrid (Bfield_fp[lev][0].get(), + Pulsar::InitializeExternalPulsarFieldsOnGrid (Bfield_fp[lev][0].get(), Bfield_fp[lev][1].get(), Bfield_fp[lev][2].get(), lev, Init_Bfield); if (lev > 0) { - PulsarParm::InitializeExternalPulsarFieldsOnGrid (Bfield_aux[lev][0].get(), + Pulsar::InitializeExternalPulsarFieldsOnGrid (Bfield_aux[lev][0].get(), Bfield_aux[lev][1].get(), Bfield_aux[lev][2].get(), lev, Init_Bfield); - PulsarParm::InitializeExternalPulsarFieldsOnGrid (Bfield_cp[lev][0].get(), + Pulsar::InitializeExternalPulsarFieldsOnGrid (Bfield_cp[lev][0].get(), Bfield_cp[lev][1].get(), Bfield_cp[lev][2].get(), lev, Init_Bfield); } } - if (PulsarParm::InitializeGrid_with_Pulsar_Efield == 1) { + if (Pulsar::m_do_InitializeGrid_with_Pulsar_Efield == 1) { bool Init_Bfield = false; - PulsarParm::InitializeExternalPulsarFieldsOnGrid (Efield_fp[lev][0].get(), + Pulsar::InitializeExternalPulsarFieldsOnGrid (Efield_fp[lev][0].get(), Efield_fp[lev][1].get(), Efield_fp[lev][2].get(), lev, Init_Bfield); if (lev > 0) { - PulsarParm::InitializeExternalPulsarFieldsOnGrid (Efield_aux[lev][0].get(), + Pulsar::InitializeExternalPulsarFieldsOnGrid (Efield_aux[lev][0].get(), Efield_aux[lev][1].get(), Efield_aux[lev][2].get(), lev, Init_Bfield); - PulsarParm::InitializeExternalPulsarFieldsOnGrid (Efield_cp[lev][0].get(), + Pulsar::InitializeExternalPulsarFieldsOnGrid (Efield_cp[lev][0].get(), Efield_cp[lev][1].get(), Efield_cp[lev][2].get(), lev, Init_Bfield); diff --git a/Source/Particles/Deposition/ChargeDeposition.H b/Source/Particles/Deposition/ChargeDeposition.H index 17b9b3231d3..8cf2d4749b6 100644 --- a/Source/Particles/Deposition/ChargeDeposition.H +++ b/Source/Particles/Deposition/ChargeDeposition.H @@ -83,6 +83,15 @@ void doChargeDepositionShapeN (const GetParticlePosition& GetPosition, constexpr int NODE = amrex::IndexType::NODE; constexpr int CELL = amrex::IndexType::CELL; +#ifdef PULSAR + int turnoff_deposition = Pulsar::m_turnoffdeposition; + amrex::GpuArray center_star_arr; + for(int i = 0; i < 3; ++i) { + center_star_arr[i] = Pulsar::m_center_star[i]; + } + amrex::Real max_nodepos_radius = Pulsar::m_max_nodepos_radius; +#endif + // Loop over particles and deposit into rho_fab #if defined(WARPX_USE_GPUCLOCK) amrex::Real* cost_real = nullptr; @@ -110,17 +119,17 @@ void doChargeDepositionShapeN (const GetParticlePosition& GetPosition, #ifdef PULSAR // temporarily adding this to check if turning off j-deposition // within a region that is corotating makes a difference - if (PulsarParm::turnoffdeposition == 1) { - const amrex::Real xc = PulsarParm::center_star[0]; - const amrex::Real yc = PulsarParm::center_star[1]; - const amrex::Real zc = PulsarParm::center_star[2]; + if (turnoff_deposition == 1) { + const amrex::Real xc = center_star_arr[0]; + const amrex::Real yc = center_star_arr[1]; + const amrex::Real zc = center_star_arr[2]; // find radius of the particle amrex::Real rp = std::sqrt( (xp-xc)*(xp-xc) + (yp-yc)*(yp-yc) + (zp-zc)*(zp-zc) ); // change the locally defined particle weight, wq, to 0 // within the region defined by r<=max_radius_nodepos - if (rp <= PulsarParm::max_nodepos_radius) { + if (rp <= max_nodepos_radius) { wq = 0.; } } diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H index 43a78cef84b..55e3d8c96df 100644 --- a/Source/Particles/Deposition/CurrentDeposition.H +++ b/Source/Particles/Deposition/CurrentDeposition.H @@ -111,6 +111,15 @@ void doDepositionShapeN(const GetParticlePosition& GetPosition, constexpr int NODE = amrex::IndexType::NODE; constexpr int CELL = amrex::IndexType::CELL; +#ifdef PULSAR + int turnoff_deposition = Pulsar::m_turnoffdeposition; + amrex::GpuArray center_star_arr; + for(int i = 0; i < 3; ++i) { + center_star_arr[i] = Pulsar::m_center_star[i]; + } + amrex::Real max_nodepos_radius = Pulsar::m_max_nodepos_radius; +#endif + // Loop over particles and deposit into jx_fab, jy_fab and jz_fab #if defined(WARPX_USE_GPUCLOCK) amrex::Real* cost_real = nullptr; @@ -142,17 +151,17 @@ void doDepositionShapeN(const GetParticlePosition& GetPosition, #ifdef PULSAR // temporarily adding this to check if turning off j-deposition // within a region that is corotating makes a difference - if (PulsarParm::turnoffdeposition == 1) { - const amrex::Real xc = PulsarParm::center_star[0]; - const amrex::Real yc = PulsarParm::center_star[1]; - const amrex::Real zc = PulsarParm::center_star[2]; + if (turnoff_deposition == 1) { + const amrex::Real xc = center_star_arr[0]; + const amrex::Real yc = center_star_arr[1]; + const amrex::Real zc = center_star_arr[2]; // find radius of the particle amrex::Real rp = std::sqrt( (xp-xc)*(xp-xc) + (yp-yc)*(yp-yc) + (zp-zc)*(zp-zc) ); // change the locally defined particle weight, wq, to 0 // within the region defined by r<=max_radius_nodepos - if (rp <= PulsarParm::max_nodepos_radius) { + if (rp <= max_nodepos_radius) { wq = 0.; } } @@ -421,6 +430,14 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, Real const clightsq = 1.0_rt / ( PhysConst::c * PhysConst::c ); +#ifdef PULSAR + int turnoff_deposition = Pulsar::m_turnoffdeposition; + amrex::GpuArray center_star_arr; + for(int i = 0; i < 3; ++i) { + center_star_arr[i] = Pulsar::m_center_star[i]; + } + amrex::Real max_nodepos_radius = Pulsar::m_max_nodepos_radius; +#endif // Loop over particles and deposit into Jx_arr, Jy_arr and Jz_arr #if defined(WARPX_USE_GPUCLOCK) amrex::Real* cost_real = nullptr; @@ -454,17 +471,17 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, #ifdef PULSAR // temporarily adding this to check if turning off j-deposition // within a region that is corotating makes a difference - if (PulsarParm::turnoffdeposition == 1) { - const amrex::Real xc = PulsarParm::center_star[0]; - const amrex::Real yc = PulsarParm::center_star[1]; - const amrex::Real zc = PulsarParm::center_star[2]; + if (turnoff_deposition == 1) { + const amrex::Real xc = center_star_arr[0]; + const amrex::Real yc = center_star_arr[1]; + const amrex::Real zc = center_star_arr[2]; // find radius of the particle amrex::Real rp = std::sqrt( (xp-xc)*(xp-xc) + (yp-yc)*(yp-yc) + (zp-zc)*(zp-zc) ); // change the locally defined particle weight, wq, to 0 // within the region defined by r<=max_radius_nodepos - if (rp <= PulsarParm::max_nodepos_radius) { + if (rp <= max_nodepos_radius) { wq = 0.; } } @@ -789,6 +806,15 @@ void doVayDepositionShapeN (const GetParticlePosition& GetPosition, amrex::Array4 const& jy_arr = jy_fab.array(); amrex::Array4 const& jz_arr = jz_fab.array(); +#ifdef PULSAR + int turnoff_deposition = Pulsar::m_turnoffdeposition; + amrex::GpuArray center_star_arr; + for(int i = 0; i < 3; ++i) { + center_star_arr[i] = Pulsar::m_center_star[i]; + } + amrex::Real max_nodepos_radius = Pulsar::m_max_nodepos_radius; +#endif + // Loop over particles and deposit (Dx,Dy,Dz) into jx_fab, jy_fab and jz_fab #if defined(WARPX_USE_GPUCLOCK) amrex::Real* cost_real = nullptr; @@ -819,17 +845,17 @@ void doVayDepositionShapeN (const GetParticlePosition& GetPosition, #ifdef PULSAR // temporarily adding this to check if turning off j-deposition // within a region that is corotating makes a difference - if (PulsarParm::turnoffdeposition == 1) { - const amrex::Real xc = PulsarParm::center_star[0]; - const amrex::Real yc = PulsarParm::center_star[1]; - const amrex::Real zc = PulsarParm::center_star[2]; + if (turnoff_deposition == 1) { + const amrex::Real xc = center_star_arr[0]; + const amrex::Real yc = center_star_arr[1]; + const amrex::Real zc = center_star_arr[2]; // find radius of the particle amrex::Real rp = std::sqrt( (xp-xc)*(xp-xc) + (yp-yc)*(yp-yc) + (zp-zc)*(zp-zc) ); // change the locally defined particle weight, wq, to 0 // within the region defined by r<=max_radius_nodepos - if (rp <= PulsarParm::max_nodepos_radius) { + if (rp <= max_nodepos_radius ) { wq = 0.; } } diff --git a/Source/Particles/Gather/FieldGather.H b/Source/Particles/Gather/FieldGather.H index 8401d752706..5ed82cbd5b5 100644 --- a/Source/Particles/Gather/FieldGather.H +++ b/Source/Particles/Gather/FieldGather.H @@ -416,7 +416,14 @@ void doGatherShapeN(const GetParticlePosition& getPosition, amrex::GpuArray dx_arr = {dx[0], dx[1], dx[2]}; amrex::GpuArray xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]}; - +#ifdef PULSAR + amrex::GpuArray center_star_arr; + for (int i = 0; i < 3; ++i) { + center_star_arr[i] = Pulsar::m_center_star[i]; + } + int do_turnoff_plasmaEB_gather_data = Pulsar::m_turnoff_plasmaEB_gather; + amrex::Real max_nogather_radius = Pulsar::m_max_nogather_radius; +#endif amrex::Array4 const& ex_arr = exfab->array(); amrex::Array4 const& ey_arr = eyfab->array(); amrex::Array4 const& ez_arr = ezfab->array(); @@ -447,16 +454,16 @@ void doGatherShapeN(const GetParticlePosition& getPosition, // of self-consistent fields within a region // makes a difference int dogather = 1; - if (PulsarParm::turnoff_plasmaEB_gather == 1) { - const amrex::Real xc = PulsarParm::center_star[0]; - const amrex::Real yc = PulsarParm::center_star[1]; - const amrex::Real zc = PulsarParm::center_star[2]; - // find radius of the particle + if (do_turnoff_plasmaEB_gather_data == 1) { + const amrex::Real xc = center_star_arr[0]; + const amrex::Real yc = center_star_arr[1]; + const amrex::Real zc = center_star_arr[2]; + // find radius of the particlestar[0]; amrex::Real rp = std::sqrt( (xp-xc)*(xp-xc) + (yp-yc)*(yp-yc) + (zp-zc)*(zp-zc) ); // set dogather to 0 if particle within user-defined nogather radius - if (rp <= PulsarParm::max_nogather_radius) { + if (rp <= max_nogather_radius) { dogather = 0; } } diff --git a/Source/Particles/Gather/GatherExternalPulsarFieldOnGrid.H b/Source/Particles/Gather/GatherExternalPulsarFieldOnGrid.H index c8e534bedb7..06e83b9ff53 100644 --- a/Source/Particles/Gather/GatherExternalPulsarFieldOnGrid.H +++ b/Source/Particles/Gather/GatherExternalPulsarFieldOnGrid.H @@ -29,534 +29,534 @@ * \param Field : 0 for Bfield, 1 for Efield * \param comp : x(ncomp=0), y(ncomp=1), or z(ncomp=2) component */ -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -amrex::Real GetExternalFieldAtCellLocation(int ii, int jj, int kk, - amrex::GpuArray const f_type, - amrex::GpuArray const domain_lo, - amrex::GpuArray const domain_hi, - amrex::GpuArray const dx, - amrex::Real cur_time, int Field, int comp) -{ - amrex::Real xcell, ycell, zcell; - amrex::Real r, theta, phi; - amrex::Real Fr, Ftheta, Fphi; - amrex::Real Fext_comp; - // compute cell coordinates, (x,y,z) corresponding to (ii,jj,kk) - // based on the staggered Bx location on the yee-grid - PulsarParm::ComputeCellCoordinates( ii, jj, kk, f_type, domain_lo, - dx, xcell, ycell, zcell ); - // Convert (x,y,z) to (r,theta,phi) - PulsarParm::ConvertCartesianToSphericalCoord( xcell, ycell, zcell, - domain_lo, domain_hi, - r, theta, phi); - if (Field == 0) { - // Compute External vacuum fields of the pulsar (Br, Btheta, Bphi) - PulsarParm::ExternalBFieldSpherical( r, theta, phi, cur_time, - Fr, Ftheta, Fphi); - } else if (Field == 1) { - // Compute External vacuum fields of the pulsar (Br, Btheta, Bphi) - PulsarParm::ExternalEFieldSpherical( r, theta, phi, cur_time, - Fr, Ftheta, Fphi); - } - - if (comp == 0) { - // Convert (Fr, Ftheta, Fphi) to Xcomponent - PulsarParm::ConvertSphericalToCartesianXComponent( Fr, Ftheta, Fphi, - r, theta, phi, Fext_comp); - } else if (comp == 1) { - // Convert (Fr, Ftheta, Fphi) to YComponent - PulsarParm::ConvertSphericalToCartesianYComponent( Fr, Ftheta, Fphi, - r, theta, phi, Fext_comp); - } else if (comp == 2) { - // Convert (Fr, Ftheta, Fphi) to ZComponent - PulsarParm::ConvertSphericalToCartesianZComponent( Fr, Ftheta, Fphi, - r, theta, phi, Fext_comp); - } - - return Fext_comp; -} - - -/** - * \brief Field gather for a single particle - * - * \tparam depos_order Particle shape order - * \tparam galerkin_interpolation Lower the order of the particle shape by - * this value (0/1) for the parallel field component - * \param xp, yp, zp Particle position coordinates - * \param Exp, Eyp, Ezp Electric field on particles. - * \param Bxp, Byp, Bzp Magnetic field on particles. - * \param ex_arr ey_arr ez_arr Array4 of the electric field, either full array or tile. - * \param bx_arr by_arr bz_arr Array4 of the magnetic field, either full array or tile. - * \param ex_type, ey_type, ez_type IndexType of the electric field - * \param bx_type, by_type, bz_type IndexType of the magnetic field - * \param dx 3D cell spacing - * \param xyzmin Physical lower bounds of domain in x, y, z. - * \param lo Index lower bounds of domain. - * \param n_rz_azimuthal_modes Number of azimuthal modes when using RZ geometry - */ -template -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void doPulsarFieldGatherShapeN (const amrex::ParticleReal xp, - const amrex::ParticleReal yp, - const amrex::ParticleReal zp, - amrex::ParticleReal& Exp, - amrex::ParticleReal& Eyp, - amrex::ParticleReal& Ezp, - amrex::ParticleReal& Bxp, - amrex::ParticleReal& Byp, - amrex::ParticleReal& Bzp, - amrex::Array4 const& ex_arr, - amrex::Array4 const& ey_arr, - amrex::Array4 const& ez_arr, - amrex::Array4 const& bx_arr, - amrex::Array4 const& by_arr, - amrex::Array4 const& bz_arr, - const amrex::IndexType ex_type, - const amrex::IndexType ey_type, - const amrex::IndexType ez_type, - const amrex::IndexType bx_type, - const amrex::IndexType by_type, - const amrex::IndexType bz_type, - const amrex::GpuArray& dx, - const amrex::GpuArray& xyzmin, - const amrex::Dim3& lo, - const long n_rz_azimuthal_modes, - amrex::GpuArray const& dom_lo = {}, - amrex::GpuArray const& dom_hi = {}, - amrex::Real cur_time = 0._rt ) -{ - using namespace amrex; - -#if defined(WARPX_DIM_XZ) - amrex::Abort("pulsar does not work in 2D"); -#endif - -#ifdef WARPX_DIM_RZ - amrex::Abort("pulsar does not work in RZ"); -#endif - - const amrex::Real dxi = 1.0_rt/dx[0]; - const amrex::Real dyi = 1.0_rt/dx[1]; - const amrex::Real dzi = 1.0_rt/dx[2]; - - const amrex::Real xmin = xyzmin[0]; - const amrex::Real ymin = xyzmin[1]; - const amrex::Real zmin = xyzmin[2]; - - constexpr int zdir = (AMREX_SPACEDIM - 1); - constexpr int NODE = amrex::IndexType::NODE; - constexpr int CELL = amrex::IndexType::CELL; - - // store staggering of ex-bz in GpuArray - amrex::GpuArray 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]; - } - // --- Compute shape factors - // x direction - // Get particle position - const amrex::Real x = (xp-xmin)*dxi; - - // j_[eb][xyz] leftmost grid point in x that the particle touches for the centering of each current - // sx_[eb][xyz] shape factor along x for the centering of each current - // There are only two possible centerings, node or cell centered, so at most only two shape factor - // arrays will be needed. - amrex::Real sx_node[depos_order + 1]; - amrex::Real sx_cell[depos_order + 1]; - amrex::Real sx_node_galerkin[depos_order + 1 - galerkin_interpolation] = {0._rt}; - amrex::Real sx_cell_galerkin[depos_order + 1 - galerkin_interpolation] = {0._rt}; - - int j_node = 0; - int j_cell = 0; - int j_node_v = 0; - int j_cell_v = 0; - Compute_shape_factor< depos_order > const compute_shape_factor; - Compute_shape_factor const compute_shape_factor_galerkin; - if ((ey_type[0] == NODE) || (ez_type[0] == NODE) || (bx_type[0] == NODE)) { - j_node = compute_shape_factor(sx_node, x); - } - if ((ey_type[0] == CELL) || (ez_type[0] == CELL) || (bx_type[0] == CELL)) { - j_cell = compute_shape_factor(sx_cell, x - 0.5_rt); - } - if ((ex_type[0] == NODE) || (by_type[0] == NODE) || (bz_type[0] == NODE)) { - j_node_v = compute_shape_factor_galerkin(sx_node_galerkin, x); - } - if ((ex_type[0] == CELL) || (by_type[0] == CELL) || (bz_type[0] == CELL)) { - j_cell_v = compute_shape_factor_galerkin(sx_cell_galerkin, x - 0.5_rt); - } - const amrex::Real (&sx_ex)[depos_order + 1 - galerkin_interpolation] = ((ex_type[0] == NODE) ? sx_node_galerkin : sx_cell_galerkin); - const amrex::Real (&sx_ey)[depos_order + 1 ] = ((ey_type[0] == NODE) ? sx_node : sx_cell ); - const amrex::Real (&sx_ez)[depos_order + 1 ] = ((ez_type[0] == NODE) ? sx_node : sx_cell ); - const amrex::Real (&sx_bx)[depos_order + 1 ] = ((bx_type[0] == NODE) ? sx_node : sx_cell ); - const amrex::Real (&sx_by)[depos_order + 1 - galerkin_interpolation] = ((by_type[0] == NODE) ? sx_node_galerkin : sx_cell_galerkin); - const amrex::Real (&sx_bz)[depos_order + 1 - galerkin_interpolation] = ((bz_type[0] == NODE) ? sx_node_galerkin : sx_cell_galerkin); - int const j_ex = ((ex_type[0] == NODE) ? j_node_v : j_cell_v); - int const j_ey = ((ey_type[0] == NODE) ? j_node : j_cell ); - int const j_ez = ((ez_type[0] == NODE) ? j_node : j_cell ); - int const j_bx = ((bx_type[0] == NODE) ? j_node : j_cell ); - int const j_by = ((by_type[0] == NODE) ? j_node_v : j_cell_v); - int const j_bz = ((bz_type[0] == NODE) ? j_node_v : j_cell_v); - - // y direction - const amrex::Real y = (yp-ymin)*dyi; - amrex::Real sy_node[depos_order + 1]; - amrex::Real sy_cell[depos_order + 1]; - amrex::Real sy_node_v[depos_order + 1 - galerkin_interpolation]; - amrex::Real sy_cell_v[depos_order + 1 - galerkin_interpolation]; - int k_node = 0; - int k_cell = 0; - int k_node_v = 0; - int k_cell_v = 0; - if ((ex_type[1] == NODE) || (ez_type[1] == NODE) || (by_type[1] == NODE)) { - k_node = compute_shape_factor(sy_node, y); - } - if ((ex_type[1] == CELL) || (ez_type[1] == CELL) || (by_type[1] == CELL)) { - k_cell = compute_shape_factor(sy_cell, y - 0.5_rt); - } - if ((ey_type[1] == NODE) || (bx_type[1] == NODE) || (bz_type[1] == NODE)) { - k_node_v = compute_shape_factor_galerkin(sy_node_v, y); - } - if ((ey_type[1] == CELL) || (bx_type[1] == CELL) || (bz_type[1] == CELL)) { - k_cell_v = compute_shape_factor_galerkin(sy_cell_v, y - 0.5_rt); - } - const amrex::Real (&sy_ex)[depos_order + 1 ] = ((ex_type[1] == NODE) ? sy_node : sy_cell ); - const amrex::Real (&sy_ey)[depos_order + 1 - galerkin_interpolation] = ((ey_type[1] == NODE) ? sy_node_v : sy_cell_v); - const amrex::Real (&sy_ez)[depos_order + 1 ] = ((ez_type[1] == NODE) ? sy_node : sy_cell ); - const amrex::Real (&sy_bx)[depos_order + 1 - galerkin_interpolation] = ((bx_type[1] == NODE) ? sy_node_v : sy_cell_v); - const amrex::Real (&sy_by)[depos_order + 1 ] = ((by_type[1] == NODE) ? sy_node : sy_cell ); - const amrex::Real (&sy_bz)[depos_order + 1 - galerkin_interpolation] = ((bz_type[1] == NODE) ? sy_node_v : sy_cell_v); - int const k_ex = ((ex_type[1] == NODE) ? k_node : k_cell ); - int const k_ey = ((ey_type[1] == NODE) ? k_node_v : k_cell_v); - int const k_ez = ((ez_type[1] == NODE) ? k_node : k_cell ); - int const k_bx = ((bx_type[1] == NODE) ? k_node_v : k_cell_v); - int const k_by = ((by_type[1] == NODE) ? k_node : k_cell ); - int const k_bz = ((bz_type[1] == NODE) ? k_node_v : k_cell_v); - - // z direction - const amrex::Real z = (zp-zmin)*dzi; - amrex::Real sz_node[depos_order + 1]; - amrex::Real sz_cell[depos_order + 1]; - amrex::Real sz_node_v[depos_order + 1 - galerkin_interpolation]; - amrex::Real sz_cell_v[depos_order + 1 - galerkin_interpolation]; - int l_node = 0; - int l_cell = 0; - int l_node_v = 0; - int l_cell_v = 0; - if ((ex_type[zdir] == NODE) || (ey_type[zdir] == NODE) || (bz_type[zdir] == NODE)) { - l_node = compute_shape_factor(sz_node, z); - } - if ((ex_type[zdir] == CELL) || (ey_type[zdir] == CELL) || (bz_type[zdir] == CELL)) { - l_cell = compute_shape_factor(sz_cell, z - 0.5_rt); - } - if ((ez_type[zdir] == NODE) || (bx_type[zdir] == NODE) || (by_type[zdir] == NODE)) { - l_node_v = compute_shape_factor_galerkin(sz_node_v, z); - } - if ((ez_type[zdir] == CELL) || (bx_type[zdir] == CELL) || (by_type[zdir] == CELL)) { - l_cell_v = compute_shape_factor_galerkin(sz_cell_v, z - 0.5_rt); - } - const amrex::Real (&sz_ex)[depos_order + 1 ] = ((ex_type[zdir] == NODE) ? sz_node : sz_cell ); - const amrex::Real (&sz_ey)[depos_order + 1 ] = ((ey_type[zdir] == NODE) ? sz_node : sz_cell ); - const amrex::Real (&sz_ez)[depos_order + 1 - galerkin_interpolation] = ((ez_type[zdir] == NODE) ? sz_node_v : sz_cell_v); - const amrex::Real (&sz_bx)[depos_order + 1 - galerkin_interpolation] = ((bx_type[zdir] == NODE) ? sz_node_v : sz_cell_v); - const amrex::Real (&sz_by)[depos_order + 1 - galerkin_interpolation] = ((by_type[zdir] == NODE) ? sz_node_v : sz_cell_v); - const amrex::Real (&sz_bz)[depos_order + 1 ] = ((bz_type[zdir] == NODE) ? sz_node : sz_cell ); - int const l_ex = ((ex_type[zdir] == NODE) ? l_node : l_cell ); - int const l_ey = ((ey_type[zdir] == NODE) ? l_node : l_cell ); - int const l_ez = ((ez_type[zdir] == NODE) ? l_node_v : l_cell_v); - int const l_bx = ((bx_type[zdir] == NODE) ? l_node_v : l_cell_v); - int const l_by = ((by_type[zdir] == NODE) ? l_node_v : l_cell_v); - int const l_bz = ((bz_type[zdir] == NODE) ? l_node : l_cell ); - - - // Each field is gathered in a separate block of - // AMREX_SPACEDIM nested loops because the deposition - // order can differ for each component of each field - // when galerkin_interpolation is set to 1 - - // Gather Pulsar external field using grid resolution on particle Exp - for (int iz=0; iz<=depos_order; iz++){ - for (int iy=0; iy<=depos_order; iy++){ - for (int ix=0; ix<= depos_order - galerkin_interpolation; ix++){ - // Add the external field contribution for pulsar - int ii = lo.x + j_ex + ix; - int jj = lo.y + k_ex + iy; - int kk = lo.z + l_ex + iz; - int Field = 1; // 1 for Efield, 0 for Bfield - int comp = 0; // 0 for xcomponent - amrex::Real ex_ext = GetExternalFieldAtCellLocation( ii, jj, kk, - Ex_stag, dom_lo, dom_hi, dx, - cur_time, Field, comp); - Exp += sx_ex[ix]*sy_ex[iy]*sz_ex[iz]*ex_ext; - } - } - } - // Gather Pulsar external field using grid resolution on particle Eyp - for (int iz=0; iz<=depos_order; iz++){ - for (int iy=0; iy<= depos_order - galerkin_interpolation; iy++){ - for (int ix=0; ix<=depos_order; ix++){ - // Add the external field contribution for pulsar - int ii = lo.x + j_ey + ix; - int jj = lo.y + k_ey + iy; - int kk = lo.z + l_ey + iz; - int Field = 1; // 1 for Efield, 0 for Bfield - int comp = 1; // 1 for ycomponent - amrex::Real ey_ext = GetExternalFieldAtCellLocation( ii, jj, kk, - Ey_stag, dom_lo, dom_hi, dx, - cur_time, Field, comp); - Eyp += sx_ey[ix]*sy_ey[iy]*sz_ey[iz]*ey_ext; - } - } - } - // Gather Pulsar external field using grid resolution on particle Ezp - for (int iz=0; iz<= depos_order - galerkin_interpolation; iz++){ - for (int iy=0; iy<=depos_order; iy++){ - for (int ix=0; ix<=depos_order; ix++){ - // Add the external field contribution for pulsar - int ii = lo.x + j_ez + ix; - int jj = lo.y + k_ez + iy; - int kk = lo.z + l_ez + iz; - int Field = 1; // 1 for Efield, 0 for Bfield - int comp = 2; // 2 for zcomponent - amrex::Real ez_ext = GetExternalFieldAtCellLocation( ii, jj, kk, - Ez_stag, dom_lo, dom_hi, dx, - cur_time, Field, comp); - Ezp += sx_ez[ix]*sy_ez[iy]*sz_ez[iz]*ez_ext; - } - } - } - // Gather Pulsar external field using grid resolution on particle Bzp - for (int iz=0; iz<=depos_order; iz++){ - for (int iy=0; iy<= depos_order - galerkin_interpolation; iy++){ - for (int ix=0; ix<= depos_order - galerkin_interpolation; ix++){ - // Add the external field contribution for pulsar - int ii = lo.x + j_bz + ix; - int jj = lo.y + k_bz + iy; - int kk = lo.z + l_bz + iz; - int Field = 0; // 1 for Efield, 0 for Bfield - int comp = 2; // 2 for zcomponent - amrex::Real bz_ext = GetExternalFieldAtCellLocation( ii, jj, kk, - Bz_stag, dom_lo, dom_hi, dx, - cur_time, Field, comp); - Bzp += sx_bz[ix]*sy_bz[iy]*sz_bz[iz]*bz_ext; - } - } - } - // Gather Pulsar external field using grid resolution on particle Byp - for (int iz=0; iz<= depos_order - galerkin_interpolation; iz++){ - for (int iy=0; iy<=depos_order; iy++){ - for (int ix=0; ix<= depos_order - galerkin_interpolation; ix++){ - // Add the external field contribution for pulsar - int ii = lo.x + j_by + ix; - int jj = lo.y + k_by + iy; - int kk = lo.z + l_by + iz; - int Field = 0; // 1 for Efield, 0 for Bfield - int comp = 1; // 1 for ycomponent - amrex::Real by_ext = GetExternalFieldAtCellLocation( ii, jj, kk, - By_stag, dom_lo, dom_hi, dx, - cur_time, Field, comp); - Byp += sx_by[ix]*sy_by[iy]*sz_by[iz]*by_ext; - } - } - } - // Gather Pulsar external field using grid resolution on particle Bxp - for (int iz=0; iz<= depos_order - galerkin_interpolation; iz++){ - for (int iy=0; iy<= depos_order - galerkin_interpolation; iy++){ - for (int ix=0; ix<=depos_order; ix++){ - // Add the external field contribution for pulsar - int ii = lo.x + j_bx + ix; - int jj = lo.y + k_bx + iy; - int kk = lo.z + l_bx + iz; - int Field = 0; // 1 for Efield, 0 for Bfield - int comp = 0; // 0 for xcomponent - amrex::Real bx_ext = GetExternalFieldAtCellLocation( ii, jj, kk, - Bx_stag, dom_lo, dom_hi, dx, - cur_time, Field, comp); - Bxp += sx_bx[ix]*sy_bx[iy]*sz_bx[iz]*bx_ext; - } - } - } -} - - - - -/** - * \brief Field gather for particles - * - * /param getPosition : A functor for returning the particle position. - * /param getExternalEField : A functor for assigning the external E field. - * /param getExternalBField : A functor for assigning the external B field. - * \param Exp, Eyp, Ezp : Pointer to array of electric field on particles. - * \param Bxp, Byp, Bzp : Pointer to array of magnetic field on particles. - * \param exfab eyfab ezfab : Array4 of the electric field, either full array or tile. - * \param ezfab bxfab bzfab : Array4 of the magnetic field, either full array or tile. - * \param np_to_gather : Number of particles for which field is gathered. - * \param dx : 3D cell size - * \param xyzmin : Physical lower bounds of domain. - * \param lo : Index lower bounds of domain. - * \param n_rz_azimuthal_modes : Number of azimuthal modes when using RZ geometry - */ -template -void doPulsarFieldGatherShapeN(const GetParticlePosition& getPosition, - const GetExternalEField& getExternalE, const GetExternalBField& getExternalB, - amrex::ParticleReal * const Exp, amrex::ParticleReal * const Eyp, - amrex::ParticleReal * const Ezp, amrex::ParticleReal * const Bxp, - amrex::ParticleReal * const Byp, amrex::ParticleReal * const Bzp, - amrex::FArrayBox const * const exfab, - amrex::FArrayBox const * const eyfab, - amrex::FArrayBox const * const ezfab, - amrex::FArrayBox const * const bxfab, - amrex::FArrayBox const * const byfab, - amrex::FArrayBox const * const bzfab, - const long np_to_gather, - const std::array& dx, - const std::array& xyzmin, - const amrex::Dim3 lo, - const long n_rz_azimuthal_modes, - amrex::GpuArray const& dom_lo = {}, - amrex::GpuArray const& dom_hi = {}, - amrex::Real cur_time = 0._rt ) -{ - - amrex::GpuArray dx_arr = {dx[0], dx[1], dx[2]}; - amrex::GpuArray xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]}; - - amrex::Array4 const& ex_arr = exfab->array(); - amrex::Array4 const& ey_arr = eyfab->array(); - amrex::Array4 const& ez_arr = ezfab->array(); - amrex::Array4 const& bx_arr = bxfab->array(); - amrex::Array4 const& by_arr = byfab->array(); - amrex::Array4 const& bz_arr = bzfab->array(); - - amrex::IndexType const ex_type = exfab->box().ixType(); - amrex::IndexType const ey_type = eyfab->box().ixType(); - amrex::IndexType const ez_type = ezfab->box().ixType(); - amrex::IndexType const bx_type = bxfab->box().ixType(); - amrex::IndexType const by_type = byfab->box().ixType(); - amrex::IndexType const bz_type = bzfab->box().ixType(); - - // Loop over particles and gather fields from - // {e,b}{x,y,z}_arr to {E,B}{xyz}p. - amrex::ParallelFor( - np_to_gather, - [=] AMREX_GPU_DEVICE (long ip) { - - amrex::ParticleReal xp, yp, zp; - getPosition(ip, xp, yp, zp); - getExternalE(ip, Exp[ip], Eyp[ip], Ezp[ip]); - getExternalB(ip, Bxp[ip], Byp[ip], Bzp[ip]); - - doPulsarFieldGatherShapeN( - xp, yp, zp, Exp[ip], Eyp[ip], Ezp[ip], Bxp[ip], Byp[ip], Bzp[ip], - ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, - ex_type, ey_type, ez_type, bx_type, by_type, bz_type, - dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes, - dom_lo, dom_hi, cur_time); - }); -} - -/** - * \brief Field gather for a single particle - * - * /param xp, yp, zp : Particle position coordinates - * \param Exp, Eyp, Ezp : Electric field on particles. - * \param Bxp, Byp, Bzp : Magnetic field on particles. - * \param ex_arr ey_arr ez_arr : Array4 of the electric field, either full array or tile. - * \param bx_arr by_arr bz_arr : Array4 of the magnetic field, either full array or tile. - * \param ex_type, ey_type, ez_type : IndexType of the electric field - * \param bx_type, by_type, bz_type : IndexType of the magnetic field - * \param dx : 3D cell spacing - * \param xyzmin : Physical lower bounds of domain in x, y, z. - * \param lo : Index lower bounds of domain. - * \param n_rz_azimuthal_modes : Number of azimuthal modes when using RZ geometry - * \param nox : order of the particle shape function - * \param galerkin_interpolation : whether to use lower order in v - */ -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void doPulsarFieldGatherShapeN (const amrex::ParticleReal xp, - const amrex::ParticleReal yp, - const amrex::ParticleReal zp, - amrex::ParticleReal& Exp, - amrex::ParticleReal& Eyp, - amrex::ParticleReal& Ezp, - amrex::ParticleReal& Bxp, - amrex::ParticleReal& Byp, - amrex::ParticleReal& Bzp, - amrex::Array4 const& ex_arr, - amrex::Array4 const& ey_arr, - amrex::Array4 const& ez_arr, - amrex::Array4 const& bx_arr, - amrex::Array4 const& by_arr, - amrex::Array4 const& bz_arr, - const amrex::IndexType ex_type, - const amrex::IndexType ey_type, - const amrex::IndexType ez_type, - const amrex::IndexType bx_type, - const amrex::IndexType by_type, - const amrex::IndexType bz_type, - const amrex::GpuArray& dx_arr, - const amrex::GpuArray& xyzmin_arr, - const amrex::Dim3& lo, - const int n_rz_azimuthal_modes, - const int nox, - const bool galerkin_interpolation, - amrex::GpuArray const& dom_lo = {}, - amrex::GpuArray const& dom_hi = {}, - amrex::Real cur_time = 0._rt ) -{ - if (galerkin_interpolation) { - if (nox == 1) { - doPulsarFieldGatherShapeN<1,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp, - ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, - ex_type, ey_type, ez_type, bx_type, by_type, bz_type, - dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes, - dom_lo, dom_hi, cur_time - ); - } else if (nox == 2) { - doPulsarFieldGatherShapeN<2,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp, - ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, - ex_type, ey_type, ez_type, bx_type, by_type, bz_type, - dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes, - dom_lo, dom_hi, cur_time - ); - } else if (nox == 3) { - doPulsarFieldGatherShapeN<3,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp, - ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, - ex_type, ey_type, ez_type, bx_type, by_type, bz_type, - dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes, - dom_lo, dom_hi, cur_time - ); - } - } else { - if (nox == 1) { - doPulsarFieldGatherShapeN<1,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp, - ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, - ex_type, ey_type, ez_type, bx_type, by_type, bz_type, - dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes, - dom_lo, dom_hi, cur_time - ); - } else if (nox == 2) { - doPulsarFieldGatherShapeN<2,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp, - ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, - ex_type, ey_type, ez_type, bx_type, by_type, bz_type, - dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes, - dom_lo, dom_hi, cur_time); - } else if (nox == 3) { - doPulsarFieldGatherShapeN<3,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp, - ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, - ex_type, ey_type, ez_type, bx_type, by_type, bz_type, - dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes, - dom_lo, dom_hi, cur_time); - } - } -} - +//AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +//amrex::Real GetExternalFieldAtCellLocation(int ii, int jj, int kk, +// amrex::GpuArray const f_type, +// amrex::GpuArray const domain_lo, +// amrex::GpuArray const domain_hi, +// amrex::GpuArray const dx, +// amrex::Real cur_time, int Field, int comp) +//{ +// amrex::Real xcell, ycell, zcell; +// amrex::Real r, theta, phi; +// amrex::Real Fr, Ftheta, Fphi; +// amrex::Real Fext_comp; +// // compute cell coordinates, (x,y,z) corresponding to (ii,jj,kk) +// // based on the staggered Bx location on the yee-grid +// PulsarParm::ComputeCellCoordinates( ii, jj, kk, f_type, domain_lo, +// dx, xcell, ycell, zcell ); +// // Convert (x,y,z) to (r,theta,phi) +// PulsarParm::ConvertCartesianToSphericalCoord( xcell, ycell, zcell, +// domain_lo, domain_hi, +// r, theta, phi); +// if (Field == 0) { +// // Compute External vacuum fields of the pulsar (Br, Btheta, Bphi) +// PulsarParm::ExternalBFieldSpherical( r, theta, phi, cur_time, +// Fr, Ftheta, Fphi); +// } else if (Field == 1) { +// // Compute External vacuum fields of the pulsar (Br, Btheta, Bphi) +// PulsarParm::ExternalEFieldSpherical( r, theta, phi, cur_time, +// Fr, Ftheta, Fphi); +// } +// +// if (comp == 0) { +// // Convert (Fr, Ftheta, Fphi) to Xcomponent +// PulsarParm::ConvertSphericalToCartesianXComponent( Fr, Ftheta, Fphi, +// r, theta, phi, Fext_comp); +// } else if (comp == 1) { +// // Convert (Fr, Ftheta, Fphi) to YComponent +// PulsarParm::ConvertSphericalToCartesianYComponent( Fr, Ftheta, Fphi, +// r, theta, phi, Fext_comp); +// } else if (comp == 2) { +// // Convert (Fr, Ftheta, Fphi) to ZComponent +// PulsarParm::ConvertSphericalToCartesianZComponent( Fr, Ftheta, Fphi, +// r, theta, phi, Fext_comp); +// } +// +// return Fext_comp; +//} +// +// +///** +// * \brief Field gather for a single particle +// * +// * \tparam depos_order Particle shape order +// * \tparam galerkin_interpolation Lower the order of the particle shape by +// * this value (0/1) for the parallel field component +// * \param xp, yp, zp Particle position coordinates +// * \param Exp, Eyp, Ezp Electric field on particles. +// * \param Bxp, Byp, Bzp Magnetic field on particles. +// * \param ex_arr ey_arr ez_arr Array4 of the electric field, either full array or tile. +// * \param bx_arr by_arr bz_arr Array4 of the magnetic field, either full array or tile. +// * \param ex_type, ey_type, ez_type IndexType of the electric field +// * \param bx_type, by_type, bz_type IndexType of the magnetic field +// * \param dx 3D cell spacing +// * \param xyzmin Physical lower bounds of domain in x, y, z. +// * \param lo Index lower bounds of domain. +// * \param n_rz_azimuthal_modes Number of azimuthal modes when using RZ geometry +// */ +//template +//AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +//void doPulsarFieldGatherShapeN (const amrex::ParticleReal xp, +// const amrex::ParticleReal yp, +// const amrex::ParticleReal zp, +// amrex::ParticleReal& Exp, +// amrex::ParticleReal& Eyp, +// amrex::ParticleReal& Ezp, +// amrex::ParticleReal& Bxp, +// amrex::ParticleReal& Byp, +// amrex::ParticleReal& Bzp, +// amrex::Array4 const& ex_arr, +// amrex::Array4 const& ey_arr, +// amrex::Array4 const& ez_arr, +// amrex::Array4 const& bx_arr, +// amrex::Array4 const& by_arr, +// amrex::Array4 const& bz_arr, +// const amrex::IndexType ex_type, +// const amrex::IndexType ey_type, +// const amrex::IndexType ez_type, +// const amrex::IndexType bx_type, +// const amrex::IndexType by_type, +// const amrex::IndexType bz_type, +// const amrex::GpuArray& dx, +// const amrex::GpuArray& xyzmin, +// const amrex::Dim3& lo, +// const long n_rz_azimuthal_modes, +// amrex::GpuArray const& dom_lo = {}, +// amrex::GpuArray const& dom_hi = {}, +// amrex::Real cur_time = 0._rt ) +//{ +// using namespace amrex; +// +//#if defined(WARPX_DIM_XZ) +// amrex::Abort("pulsar does not work in 2D"); +//#endif +// +//#ifdef WARPX_DIM_RZ +// amrex::Abort("pulsar does not work in RZ"); +//#endif +// +// const amrex::Real dxi = 1.0_rt/dx[0]; +// const amrex::Real dyi = 1.0_rt/dx[1]; +// const amrex::Real dzi = 1.0_rt/dx[2]; +// +// const amrex::Real xmin = xyzmin[0]; +// const amrex::Real ymin = xyzmin[1]; +// const amrex::Real zmin = xyzmin[2]; +// +// constexpr int zdir = (AMREX_SPACEDIM - 1); +// constexpr int NODE = amrex::IndexType::NODE; +// constexpr int CELL = amrex::IndexType::CELL; +// +// // store staggering of ex-bz in GpuArray +// amrex::GpuArray 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]; +// } +// // --- Compute shape factors +// // x direction +// // Get particle position +// const amrex::Real x = (xp-xmin)*dxi; +// +// // j_[eb][xyz] leftmost grid point in x that the particle touches for the centering of each current +// // sx_[eb][xyz] shape factor along x for the centering of each current +// // There are only two possible centerings, node or cell centered, so at most only two shape factor +// // arrays will be needed. +// amrex::Real sx_node[depos_order + 1]; +// amrex::Real sx_cell[depos_order + 1]; +// amrex::Real sx_node_galerkin[depos_order + 1 - galerkin_interpolation] = {0._rt}; +// amrex::Real sx_cell_galerkin[depos_order + 1 - galerkin_interpolation] = {0._rt}; +// +// int j_node = 0; +// int j_cell = 0; +// int j_node_v = 0; +// int j_cell_v = 0; +// Compute_shape_factor< depos_order > const compute_shape_factor; +// Compute_shape_factor const compute_shape_factor_galerkin; +// if ((ey_type[0] == NODE) || (ez_type[0] == NODE) || (bx_type[0] == NODE)) { +// j_node = compute_shape_factor(sx_node, x); +// } +// if ((ey_type[0] == CELL) || (ez_type[0] == CELL) || (bx_type[0] == CELL)) { +// j_cell = compute_shape_factor(sx_cell, x - 0.5_rt); +// } +// if ((ex_type[0] == NODE) || (by_type[0] == NODE) || (bz_type[0] == NODE)) { +// j_node_v = compute_shape_factor_galerkin(sx_node_galerkin, x); +// } +// if ((ex_type[0] == CELL) || (by_type[0] == CELL) || (bz_type[0] == CELL)) { +// j_cell_v = compute_shape_factor_galerkin(sx_cell_galerkin, x - 0.5_rt); +// } +// const amrex::Real (&sx_ex)[depos_order + 1 - galerkin_interpolation] = ((ex_type[0] == NODE) ? sx_node_galerkin : sx_cell_galerkin); +// const amrex::Real (&sx_ey)[depos_order + 1 ] = ((ey_type[0] == NODE) ? sx_node : sx_cell ); +// const amrex::Real (&sx_ez)[depos_order + 1 ] = ((ez_type[0] == NODE) ? sx_node : sx_cell ); +// const amrex::Real (&sx_bx)[depos_order + 1 ] = ((bx_type[0] == NODE) ? sx_node : sx_cell ); +// const amrex::Real (&sx_by)[depos_order + 1 - galerkin_interpolation] = ((by_type[0] == NODE) ? sx_node_galerkin : sx_cell_galerkin); +// const amrex::Real (&sx_bz)[depos_order + 1 - galerkin_interpolation] = ((bz_type[0] == NODE) ? sx_node_galerkin : sx_cell_galerkin); +// int const j_ex = ((ex_type[0] == NODE) ? j_node_v : j_cell_v); +// int const j_ey = ((ey_type[0] == NODE) ? j_node : j_cell ); +// int const j_ez = ((ez_type[0] == NODE) ? j_node : j_cell ); +// int const j_bx = ((bx_type[0] == NODE) ? j_node : j_cell ); +// int const j_by = ((by_type[0] == NODE) ? j_node_v : j_cell_v); +// int const j_bz = ((bz_type[0] == NODE) ? j_node_v : j_cell_v); +// +// // y direction +// const amrex::Real y = (yp-ymin)*dyi; +// amrex::Real sy_node[depos_order + 1]; +// amrex::Real sy_cell[depos_order + 1]; +// amrex::Real sy_node_v[depos_order + 1 - galerkin_interpolation]; +// amrex::Real sy_cell_v[depos_order + 1 - galerkin_interpolation]; +// int k_node = 0; +// int k_cell = 0; +// int k_node_v = 0; +// int k_cell_v = 0; +// if ((ex_type[1] == NODE) || (ez_type[1] == NODE) || (by_type[1] == NODE)) { +// k_node = compute_shape_factor(sy_node, y); +// } +// if ((ex_type[1] == CELL) || (ez_type[1] == CELL) || (by_type[1] == CELL)) { +// k_cell = compute_shape_factor(sy_cell, y - 0.5_rt); +// } +// if ((ey_type[1] == NODE) || (bx_type[1] == NODE) || (bz_type[1] == NODE)) { +// k_node_v = compute_shape_factor_galerkin(sy_node_v, y); +// } +// if ((ey_type[1] == CELL) || (bx_type[1] == CELL) || (bz_type[1] == CELL)) { +// k_cell_v = compute_shape_factor_galerkin(sy_cell_v, y - 0.5_rt); +// } +// const amrex::Real (&sy_ex)[depos_order + 1 ] = ((ex_type[1] == NODE) ? sy_node : sy_cell ); +// const amrex::Real (&sy_ey)[depos_order + 1 - galerkin_interpolation] = ((ey_type[1] == NODE) ? sy_node_v : sy_cell_v); +// const amrex::Real (&sy_ez)[depos_order + 1 ] = ((ez_type[1] == NODE) ? sy_node : sy_cell ); +// const amrex::Real (&sy_bx)[depos_order + 1 - galerkin_interpolation] = ((bx_type[1] == NODE) ? sy_node_v : sy_cell_v); +// const amrex::Real (&sy_by)[depos_order + 1 ] = ((by_type[1] == NODE) ? sy_node : sy_cell ); +// const amrex::Real (&sy_bz)[depos_order + 1 - galerkin_interpolation] = ((bz_type[1] == NODE) ? sy_node_v : sy_cell_v); +// int const k_ex = ((ex_type[1] == NODE) ? k_node : k_cell ); +// int const k_ey = ((ey_type[1] == NODE) ? k_node_v : k_cell_v); +// int const k_ez = ((ez_type[1] == NODE) ? k_node : k_cell ); +// int const k_bx = ((bx_type[1] == NODE) ? k_node_v : k_cell_v); +// int const k_by = ((by_type[1] == NODE) ? k_node : k_cell ); +// int const k_bz = ((bz_type[1] == NODE) ? k_node_v : k_cell_v); +// +// // z direction +// const amrex::Real z = (zp-zmin)*dzi; +// amrex::Real sz_node[depos_order + 1]; +// amrex::Real sz_cell[depos_order + 1]; +// amrex::Real sz_node_v[depos_order + 1 - galerkin_interpolation]; +// amrex::Real sz_cell_v[depos_order + 1 - galerkin_interpolation]; +// int l_node = 0; +// int l_cell = 0; +// int l_node_v = 0; +// int l_cell_v = 0; +// if ((ex_type[zdir] == NODE) || (ey_type[zdir] == NODE) || (bz_type[zdir] == NODE)) { +// l_node = compute_shape_factor(sz_node, z); +// } +// if ((ex_type[zdir] == CELL) || (ey_type[zdir] == CELL) || (bz_type[zdir] == CELL)) { +// l_cell = compute_shape_factor(sz_cell, z - 0.5_rt); +// } +// if ((ez_type[zdir] == NODE) || (bx_type[zdir] == NODE) || (by_type[zdir] == NODE)) { +// l_node_v = compute_shape_factor_galerkin(sz_node_v, z); +// } +// if ((ez_type[zdir] == CELL) || (bx_type[zdir] == CELL) || (by_type[zdir] == CELL)) { +// l_cell_v = compute_shape_factor_galerkin(sz_cell_v, z - 0.5_rt); +// } +// const amrex::Real (&sz_ex)[depos_order + 1 ] = ((ex_type[zdir] == NODE) ? sz_node : sz_cell ); +// const amrex::Real (&sz_ey)[depos_order + 1 ] = ((ey_type[zdir] == NODE) ? sz_node : sz_cell ); +// const amrex::Real (&sz_ez)[depos_order + 1 - galerkin_interpolation] = ((ez_type[zdir] == NODE) ? sz_node_v : sz_cell_v); +// const amrex::Real (&sz_bx)[depos_order + 1 - galerkin_interpolation] = ((bx_type[zdir] == NODE) ? sz_node_v : sz_cell_v); +// const amrex::Real (&sz_by)[depos_order + 1 - galerkin_interpolation] = ((by_type[zdir] == NODE) ? sz_node_v : sz_cell_v); +// const amrex::Real (&sz_bz)[depos_order + 1 ] = ((bz_type[zdir] == NODE) ? sz_node : sz_cell ); +// int const l_ex = ((ex_type[zdir] == NODE) ? l_node : l_cell ); +// int const l_ey = ((ey_type[zdir] == NODE) ? l_node : l_cell ); +// int const l_ez = ((ez_type[zdir] == NODE) ? l_node_v : l_cell_v); +// int const l_bx = ((bx_type[zdir] == NODE) ? l_node_v : l_cell_v); +// int const l_by = ((by_type[zdir] == NODE) ? l_node_v : l_cell_v); +// int const l_bz = ((bz_type[zdir] == NODE) ? l_node : l_cell ); +// +// +// // Each field is gathered in a separate block of +// // AMREX_SPACEDIM nested loops because the deposition +// // order can differ for each component of each field +// // when galerkin_interpolation is set to 1 +// +// // Gather Pulsar external field using grid resolution on particle Exp +// for (int iz=0; iz<=depos_order; iz++){ +// for (int iy=0; iy<=depos_order; iy++){ +// for (int ix=0; ix<= depos_order - galerkin_interpolation; ix++){ +// // Add the external field contribution for pulsar +// int ii = lo.x + j_ex + ix; +// int jj = lo.y + k_ex + iy; +// int kk = lo.z + l_ex + iz; +// int Field = 1; // 1 for Efield, 0 for Bfield +// int comp = 0; // 0 for xcomponent +// amrex::Real ex_ext = GetExternalFieldAtCellLocation( ii, jj, kk, +// Ex_stag, dom_lo, dom_hi, dx, +// cur_time, Field, comp); +// Exp += sx_ex[ix]*sy_ex[iy]*sz_ex[iz]*ex_ext; +// } +// } +// } +// // Gather Pulsar external field using grid resolution on particle Eyp +// for (int iz=0; iz<=depos_order; iz++){ +// for (int iy=0; iy<= depos_order - galerkin_interpolation; iy++){ +// for (int ix=0; ix<=depos_order; ix++){ +// // Add the external field contribution for pulsar +// int ii = lo.x + j_ey + ix; +// int jj = lo.y + k_ey + iy; +// int kk = lo.z + l_ey + iz; +// int Field = 1; // 1 for Efield, 0 for Bfield +// int comp = 1; // 1 for ycomponent +// amrex::Real ey_ext = GetExternalFieldAtCellLocation( ii, jj, kk, +// Ey_stag, dom_lo, dom_hi, dx, +// cur_time, Field, comp); +// Eyp += sx_ey[ix]*sy_ey[iy]*sz_ey[iz]*ey_ext; +// } +// } +// } +// // Gather Pulsar external field using grid resolution on particle Ezp +// for (int iz=0; iz<= depos_order - galerkin_interpolation; iz++){ +// for (int iy=0; iy<=depos_order; iy++){ +// for (int ix=0; ix<=depos_order; ix++){ +// // Add the external field contribution for pulsar +// int ii = lo.x + j_ez + ix; +// int jj = lo.y + k_ez + iy; +// int kk = lo.z + l_ez + iz; +// int Field = 1; // 1 for Efield, 0 for Bfield +// int comp = 2; // 2 for zcomponent +// amrex::Real ez_ext = GetExternalFieldAtCellLocation( ii, jj, kk, +// Ez_stag, dom_lo, dom_hi, dx, +// cur_time, Field, comp); +// Ezp += sx_ez[ix]*sy_ez[iy]*sz_ez[iz]*ez_ext; +// } +// } +// } +// // Gather Pulsar external field using grid resolution on particle Bzp +// for (int iz=0; iz<=depos_order; iz++){ +// for (int iy=0; iy<= depos_order - galerkin_interpolation; iy++){ +// for (int ix=0; ix<= depos_order - galerkin_interpolation; ix++){ +// // Add the external field contribution for pulsar +// int ii = lo.x + j_bz + ix; +// int jj = lo.y + k_bz + iy; +// int kk = lo.z + l_bz + iz; +// int Field = 0; // 1 for Efield, 0 for Bfield +// int comp = 2; // 2 for zcomponent +// amrex::Real bz_ext = GetExternalFieldAtCellLocation( ii, jj, kk, +// Bz_stag, dom_lo, dom_hi, dx, +// cur_time, Field, comp); +// Bzp += sx_bz[ix]*sy_bz[iy]*sz_bz[iz]*bz_ext; +// } +// } +// } +// // Gather Pulsar external field using grid resolution on particle Byp +// for (int iz=0; iz<= depos_order - galerkin_interpolation; iz++){ +// for (int iy=0; iy<=depos_order; iy++){ +// for (int ix=0; ix<= depos_order - galerkin_interpolation; ix++){ +// // Add the external field contribution for pulsar +// int ii = lo.x + j_by + ix; +// int jj = lo.y + k_by + iy; +// int kk = lo.z + l_by + iz; +// int Field = 0; // 1 for Efield, 0 for Bfield +// int comp = 1; // 1 for ycomponent +// amrex::Real by_ext = GetExternalFieldAtCellLocation( ii, jj, kk, +// By_stag, dom_lo, dom_hi, dx, +// cur_time, Field, comp); +// Byp += sx_by[ix]*sy_by[iy]*sz_by[iz]*by_ext; +// } +// } +// } +// // Gather Pulsar external field using grid resolution on particle Bxp +// for (int iz=0; iz<= depos_order - galerkin_interpolation; iz++){ +// for (int iy=0; iy<= depos_order - galerkin_interpolation; iy++){ +// for (int ix=0; ix<=depos_order; ix++){ +// // Add the external field contribution for pulsar +// int ii = lo.x + j_bx + ix; +// int jj = lo.y + k_bx + iy; +// int kk = lo.z + l_bx + iz; +// int Field = 0; // 1 for Efield, 0 for Bfield +// int comp = 0; // 0 for xcomponent +// amrex::Real bx_ext = GetExternalFieldAtCellLocation( ii, jj, kk, +// Bx_stag, dom_lo, dom_hi, dx, +// cur_time, Field, comp); +// Bxp += sx_bx[ix]*sy_bx[iy]*sz_bx[iz]*bx_ext; +// } +// } +// } +//} +// +// +// +// +///** +// * \brief Field gather for particles +// * +// * /param getPosition : A functor for returning the particle position. +// * /param getExternalEField : A functor for assigning the external E field. +// * /param getExternalBField : A functor for assigning the external B field. +// * \param Exp, Eyp, Ezp : Pointer to array of electric field on particles. +// * \param Bxp, Byp, Bzp : Pointer to array of magnetic field on particles. +// * \param exfab eyfab ezfab : Array4 of the electric field, either full array or tile. +// * \param ezfab bxfab bzfab : Array4 of the magnetic field, either full array or tile. +// * \param np_to_gather : Number of particles for which field is gathered. +// * \param dx : 3D cell size +// * \param xyzmin : Physical lower bounds of domain. +// * \param lo : Index lower bounds of domain. +// * \param n_rz_azimuthal_modes : Number of azimuthal modes when using RZ geometry +// */ +//template +//void doPulsarFieldGatherShapeN(const GetParticlePosition& getPosition, +// const GetExternalEField& getExternalE, const GetExternalBField& getExternalB, +// amrex::ParticleReal * const Exp, amrex::ParticleReal * const Eyp, +// amrex::ParticleReal * const Ezp, amrex::ParticleReal * const Bxp, +// amrex::ParticleReal * const Byp, amrex::ParticleReal * const Bzp, +// amrex::FArrayBox const * const exfab, +// amrex::FArrayBox const * const eyfab, +// amrex::FArrayBox const * const ezfab, +// amrex::FArrayBox const * const bxfab, +// amrex::FArrayBox const * const byfab, +// amrex::FArrayBox const * const bzfab, +// const long np_to_gather, +// const std::array& dx, +// const std::array& xyzmin, +// const amrex::Dim3 lo, +// const long n_rz_azimuthal_modes, +// amrex::GpuArray const& dom_lo = {}, +// amrex::GpuArray const& dom_hi = {}, +// amrex::Real cur_time = 0._rt ) +//{ +// +// amrex::GpuArray dx_arr = {dx[0], dx[1], dx[2]}; +// amrex::GpuArray xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]}; +// +// amrex::Array4 const& ex_arr = exfab->array(); +// amrex::Array4 const& ey_arr = eyfab->array(); +// amrex::Array4 const& ez_arr = ezfab->array(); +// amrex::Array4 const& bx_arr = bxfab->array(); +// amrex::Array4 const& by_arr = byfab->array(); +// amrex::Array4 const& bz_arr = bzfab->array(); +// +// amrex::IndexType const ex_type = exfab->box().ixType(); +// amrex::IndexType const ey_type = eyfab->box().ixType(); +// amrex::IndexType const ez_type = ezfab->box().ixType(); +// amrex::IndexType const bx_type = bxfab->box().ixType(); +// amrex::IndexType const by_type = byfab->box().ixType(); +// amrex::IndexType const bz_type = bzfab->box().ixType(); +// +// // Loop over particles and gather fields from +// // {e,b}{x,y,z}_arr to {E,B}{xyz}p. +// amrex::ParallelFor( +// np_to_gather, +// [=] AMREX_GPU_DEVICE (long ip) { +// +// amrex::ParticleReal xp, yp, zp; +// getPosition(ip, xp, yp, zp); +// getExternalE(ip, Exp[ip], Eyp[ip], Ezp[ip]); +// getExternalB(ip, Bxp[ip], Byp[ip], Bzp[ip]); +// +// doPulsarFieldGatherShapeN( +// xp, yp, zp, Exp[ip], Eyp[ip], Ezp[ip], Bxp[ip], Byp[ip], Bzp[ip], +// ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, +// ex_type, ey_type, ez_type, bx_type, by_type, bz_type, +// dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes, +// dom_lo, dom_hi, cur_time); +// }); +//} +// +///** +// * \brief Field gather for a single particle +// * +// * /param xp, yp, zp : Particle position coordinates +// * \param Exp, Eyp, Ezp : Electric field on particles. +// * \param Bxp, Byp, Bzp : Magnetic field on particles. +// * \param ex_arr ey_arr ez_arr : Array4 of the electric field, either full array or tile. +// * \param bx_arr by_arr bz_arr : Array4 of the magnetic field, either full array or tile. +// * \param ex_type, ey_type, ez_type : IndexType of the electric field +// * \param bx_type, by_type, bz_type : IndexType of the magnetic field +// * \param dx : 3D cell spacing +// * \param xyzmin : Physical lower bounds of domain in x, y, z. +// * \param lo : Index lower bounds of domain. +// * \param n_rz_azimuthal_modes : Number of azimuthal modes when using RZ geometry +// * \param nox : order of the particle shape function +// * \param galerkin_interpolation : whether to use lower order in v +// */ +//AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +//void doPulsarFieldGatherShapeN (const amrex::ParticleReal xp, +// const amrex::ParticleReal yp, +// const amrex::ParticleReal zp, +// amrex::ParticleReal& Exp, +// amrex::ParticleReal& Eyp, +// amrex::ParticleReal& Ezp, +// amrex::ParticleReal& Bxp, +// amrex::ParticleReal& Byp, +// amrex::ParticleReal& Bzp, +// amrex::Array4 const& ex_arr, +// amrex::Array4 const& ey_arr, +// amrex::Array4 const& ez_arr, +// amrex::Array4 const& bx_arr, +// amrex::Array4 const& by_arr, +// amrex::Array4 const& bz_arr, +// const amrex::IndexType ex_type, +// const amrex::IndexType ey_type, +// const amrex::IndexType ez_type, +// const amrex::IndexType bx_type, +// const amrex::IndexType by_type, +// const amrex::IndexType bz_type, +// const amrex::GpuArray& dx_arr, +// const amrex::GpuArray& xyzmin_arr, +// const amrex::Dim3& lo, +// const int n_rz_azimuthal_modes, +// const int nox, +// const bool galerkin_interpolation, +// amrex::GpuArray const& dom_lo = {}, +// amrex::GpuArray const& dom_hi = {}, +// amrex::Real cur_time = 0._rt ) +//{ +// if (galerkin_interpolation) { +// if (nox == 1) { +// doPulsarFieldGatherShapeN<1,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp, +// ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, +// ex_type, ey_type, ez_type, bx_type, by_type, bz_type, +// dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes, +// dom_lo, dom_hi, cur_time +// ); +// } else if (nox == 2) { +// doPulsarFieldGatherShapeN<2,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp, +// ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, +// ex_type, ey_type, ez_type, bx_type, by_type, bz_type, +// dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes, +// dom_lo, dom_hi, cur_time +// ); +// } else if (nox == 3) { +// doPulsarFieldGatherShapeN<3,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp, +// ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, +// ex_type, ey_type, ez_type, bx_type, by_type, bz_type, +// dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes, +// dom_lo, dom_hi, cur_time +// ); +// } +// } else { +// if (nox == 1) { +// doPulsarFieldGatherShapeN<1,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp, +// ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, +// ex_type, ey_type, ez_type, bx_type, by_type, bz_type, +// dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes, +// dom_lo, dom_hi, cur_time +// ); +// } else if (nox == 2) { +// doPulsarFieldGatherShapeN<2,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp, +// ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, +// ex_type, ey_type, ez_type, bx_type, by_type, bz_type, +// dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes, +// dom_lo, dom_hi, cur_time); +// } else if (nox == 3) { +// doPulsarFieldGatherShapeN<3,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp, +// ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, +// ex_type, ey_type, ez_type, bx_type, by_type, bz_type, +// dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes, +// dom_lo, dom_hi, cur_time); +// } +// } +//} +// #endif // ifdef pulsar #endif // GATHEREXTERNALPULSARFIELDONGRID_H_ diff --git a/Source/Particles/Gather/GetExternalFields.H b/Source/Particles/Gather/GetExternalFields.H index 03127957e75..0bec101ffc1 100644 --- a/Source/Particles/Gather/GetExternalFields.H +++ b/Source/Particles/Gather/GetExternalFields.H @@ -57,7 +57,7 @@ struct GetExternalField field_y += m_field_value[1]; field_z += m_field_value[2]; } - else if (m_type == Parser) + else if (m_type == ExternalFieldInitType::Parser) { amrex::ParticleReal x, y, z; m_get_position(i, x, y, z); diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index ec9c2615f93..87485910c47 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -690,7 +690,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) WARPX_PROFILE("PhysicalParticleContainer::AddPlasma()"); #ifdef PULSAR - amrex::Print() << " in add plasma" << PulsarParm::R_star << "\n"; + amrex::Print() << " in add plasma" << Pulsar::m_R_star << "\n"; #endif // If no part_realbox is provided, initialize particles in the whole domain @@ -716,8 +716,8 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) } #ifdef PULSAR // modifying particle weight to achieve fractional injection - if (PulsarParm::ModifyParticleWtAtInjection == 1) { - scale_fac = scale_fac*PulsarParm::Ninj_fraction; + if (Pulsar::m_ModifyParticleWtAtInjection == 1) { + scale_fac = scale_fac*Pulsar::m_Ninj_fraction; } #endif @@ -756,6 +756,17 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) const MultiFab& Ez_mf = WarpX::GetInstance().getEfield(lev,2); const MultiFab& rho_mf = WarpX::GetInstance().getrho_fp(lev); const Real dt = WarpX::GetInstance().getdt(0); + amrex::GpuArray center_star_arr; + for (int i = 0; i < 3; ++i) { + center_star_arr[i] = Pulsar::m_center_star[i]; + } + amrex::Real particle_inject_rmin = Pulsar::m_particle_inject_rmin; + amrex::Real particle_inject_rmax = Pulsar::m_particle_inject_rmax; + amrex::Real dR_star = Pulsar::m_dR_star; + int ModifyParticleWtAtInjection = Pulsar::m_ModifyParticleWtAtInjection; + amrex::Real Ninj_fraction = Pulsar::m_Ninj_fraction; + amrex::Real removeparticle_theta_min = Pulsar::m_removeparticle_theta_min; + amrex::Real removeparticle_theta_max = Pulsar::m_removeparticle_theta_max; #endif #ifdef WARPX_DIM_RZ @@ -856,35 +867,35 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) lo.z = applyBallisticCorrection(lo, inj_mom, gamma_boost, beta_boost, t); hi.z = applyBallisticCorrection(hi, inj_mom, gamma_boost, beta_boost, t); #ifdef PULSAR - amrex::Real xc = PulsarParm::center_star[0]; - amrex::Real yc = PulsarParm::center_star[1]; - amrex::Real zc = PulsarParm::center_star[2]; + amrex::Real xc = center_star_arr[0]; + amrex::Real yc = center_star_arr[1]; + amrex::Real zc = center_star_arr[2]; // Find cell-center amrex::Real x, y, z; x = overlap_corner[0] + i*dx[0] + 0.5*dx[0]; y = overlap_corner[1] + j*dx[1] + 0.5*dx[1]; z = overlap_corner[2] + k*dx[2] + 0.5*dx[2]; // radius of the cell-center - amrex::Real rad = std::sqrt( (x-xc)*(x-xc) + (y-yc)*(y-yc) + (z-zc)*(z-zc)); + amrex::Real rad = std::sqrt( (x-xc)*(x-xc) + (y-yc)*(y-yc) + (z-zc)*(z-zc)); // Adding buffer-factor to ensure all cells that intersect the ring // inject particles amrex::Real buffer_factor = 0.5; // is cell-center inside the pulsar ring - if (inj_pos->insidePulsarBoundsCC( rad, PulsarParm::particle_inject_rmin, - PulsarParm::particle_inject_rmax, - PulsarParm::dR_star*buffer_factor) ) + if (inj_pos->insidePulsarBoundsCC( rad, particle_inject_rmin, + particle_inject_rmax, + dR_star*buffer_factor) ) { auto index = overlap_box.index(iv); const amrex::XDim3 ppc_per_dim = inj_pos->getppcInEachDim(); - if (PulsarParm::ModifyParticleWtAtInjection == 1) { + if (ModifyParticleWtAtInjection == 1) { // instead of modiying number of particles, the weight is changed pcounts[index] = num_ppc; - } else if (PulsarParm::ModifyParticleWtAtInjection == 0) { + } else if (ModifyParticleWtAtInjection == 0) { // Modiying number of particles injected // (could lead to round-off errors) - pcounts[index] = int(ppc_per_dim.x*std::cbrt(PulsarParm::Ninj_fraction)) - * int(ppc_per_dim.y*std::cbrt(PulsarParm::Ninj_fraction)) - * int(ppc_per_dim.z*std::cbrt(PulsarParm::Ninj_fraction)); + pcounts[index] = int(ppc_per_dim.x*std::cbrt(Ninj_fraction)) + * int(ppc_per_dim.y*std::cbrt(Ninj_fraction)) + * int(ppc_per_dim.z*std::cbrt(Ninj_fraction)); } } #else @@ -1004,7 +1015,12 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) p.cpu() = cpuid; const XDim3 r = - inj_pos->getPositionUnitBox(i_part, lrrfac, engine); + inj_pos->getPositionUnitBox(i_part, lrrfac, engine +#ifdef PULSAR + , ModifyParticleWtAtInjection, + Ninj_fraction +#endif + ); auto pos = getCellCoords(overlap_corner, dx, r, iv); #if (AMREX_SPACEDIM == 3) @@ -1071,26 +1087,26 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) } #ifdef PULSAR //amrex::Print() << " old xyz : " << xb << " " << yb << " " << z0 << "\n"; - amrex::Real xc = PulsarParm::center_star[0]; - amrex::Real yc = PulsarParm::center_star[1]; - amrex::Real zc = PulsarParm::center_star[2]; + amrex::Real xc = center_star_arr[0]; + amrex::Real yc = center_star_arr[1]; + amrex::Real zc = center_star_arr[2]; amrex::Real rad = std::sqrt( (xb-xc)*(xb-xc) + (yb-yc)*(yb-yc) + (z0-zc)*(z0-zc)); - if (!inj_pos->insidePulsarBounds(rad,PulsarParm::particle_inject_rmin, - PulsarParm::particle_inject_rmax)) { + if (!inj_pos->insidePulsarBounds(rad,particle_inject_rmin, + particle_inject_rmax)) { p.id() = -1; continue; } // particle theta amrex::Real theta_p = 0.0; if (rad > 0.) theta_p = std::acos((z0-zc)/rad); - if (theta_p > PulsarParm::removeparticle_theta_min*MathConst::pi/180. and - theta_p < PulsarParm::removeparticle_theta_max*MathConst::pi/180.) { + if (theta_p > removeparticle_theta_min*MathConst::pi/180. and + theta_p < removeparticle_theta_max*MathConst::pi/180.) { p.id() = -1; } - - //// temporarily testing if on positrons emitted at eq - //// and only electrons emitted at poles - //if (q_pm < 0) { // for electrons inject at pole + + //// temporarily testing if on positrons emitted at eq + //// and only electrons emitted at poles + //if (q_pm < 0) { // for electrons inject at pole // if (theta_p >= 40*MathConst::pi/180. and // theta_p <= 90*MathConst::pi/180.) { // p.id() = -1; @@ -1099,17 +1115,17 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) // theta_p <= 140*MathConst::pi/180.) { // p.id() = -1; // } - //} else if (q_pm > 0 ) { + //} else if (q_pm > 0 ) { // if ( theta_p <= 75*MathConst::pi/180. and - // theta_p >= 0.*MathConst::pi/180.) { + // theta_p >= 0.*MathConst::pi/180.) { // p.id() = -1; // } // if ( theta_p <= 180*MathConst::pi/180. and - // theta_p >= 105.*MathConst::pi/180.) { + // theta_p >= 105.*MathConst::pi/180.) { // p.id() = -1; // } - //} - + //} + #endif u = inj_mom->getMomentum(pos.x, pos.y, z0, engine); @@ -1280,7 +1296,6 @@ PhysicalParticleContainer::AddPlasmaFlux (int lev, amrex::Real dt) #endif int Nmax_particles = 0; - int valid_particles_beforeAdd = TotalNumberOfParticles(); MFItInfo info; if (do_tiling && Gpu::notInLaunchRegion()) { @@ -2106,10 +2121,21 @@ PhysicalParticleContainer::PushP (int lev, Real dt, const std::array& dx = WarpX::CellSize(std::max(lev,0)); #ifdef PULSAR - const auto problo = WarpX::GetInstance().Geom(lev).ProbLoArray(); - const auto probhi = WarpX::GetInstance().Geom(lev).ProbHiArray(); - amrex::Real cur_time = WarpX::GetInstance().gett_new(lev); auto &warpx = WarpX::GetInstance(); + const auto problo = warpx.Geom(lev).ProbLoArray(); + const auto probhi = warpx.Geom(lev).ProbHiArray(); + amrex::Real cur_time = warpx.gett_new(lev); + amrex::Real omega_star_data = Pulsar::m_omega_star; + amrex::Real ramp_omega_time_data = Pulsar::m_omega_ramp_time; + amrex::Real Bstar_data = Pulsar::m_B_star; + amrex::Real Rstar_data = Pulsar::m_R_star; + amrex::Real dRstar_data = Pulsar::m_dR_star; + amrex::Real corotatingE_maxradius_data = Pulsar::m_corotatingE_maxradius; + int E_external_monopole_data = Pulsar::m_do_E_external_monopole; + int AddExternalMonopoleOnly = Pulsar::m_AddExternalMonopoleOnly; + int use_theoreticalEB = Pulsar::m_use_theoreticalEB; + amrex::Real theory_max_rstar = Pulsar::m_theory_max_rstar; + int singleParticleTest = Pulsar::m_singleParticleTest; // amrex::AllPrintToFile("PulsarParticle") << " step : " << warpx.getistep(0) << " time : " << cur_time << "\n"; #endif @@ -2146,6 +2172,12 @@ PhysicalParticleContainer::PushP (int lev, Real dt, amrex::GpuArray dx_arr = {dx[0], dx[1], dx[2]}; amrex::GpuArray xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]}; +#ifdef PULSAR + amrex::GpuArray center_star_arr; + for (int idim = 0; idim < 3; ++idim) { + center_star_arr[idim] = Pulsar::m_center_star[idim]; + } +#endif amrex::Array4 const& ex_arr = exfab.array(); amrex::Array4 const& ey_arr = eyfab.array(); @@ -2167,8 +2199,10 @@ PhysicalParticleContainer::PushP (int lev, Real dt, ParticleReal* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); amrex::Gpu::synchronize(); +#ifdef PULSAR amrex::Gpu::ManagedVector PulsarParticleDiag(55,0.0); amrex::Real * PulsarParticleDiagData = PulsarParticleDiag.data(); +#endif int* AMREX_RESTRICT ion_lev = nullptr; if (do_field_ionization) { @@ -2195,7 +2229,8 @@ PhysicalParticleContainer::PushP (int lev, Real dt, // auto& p = particles[ip]; // amrex::AllPrintToFile("PulsarParticle") << " part id : " << p.id() << " ip: " << ip << " q: " << q << " xp : " << xp << " yp " << yp << " zp " << zp << "\n"; amrex::Real r_p, theta_p, phi_p; - PulsarParm::ConvertCartesianToSphericalCoord( xp, yp, zp, problo, probhi, r_p, theta_p, phi_p); + Pulsar::ConvertCartesianToSphericalCoord( xp, yp, zp, center_star_arr, + r_p, theta_p, phi_p); // amrex::AllPrintToFile("PulsarParticle") << " rp : " << r_p << " thetap " << theta_p << " phip " << phi_p << "\n"; #endif amrex::ParticleReal Exp = 0._rt, Eyp = 0._rt, Ezp = 0._rt; @@ -2214,30 +2249,32 @@ PhysicalParticleContainer::PushP (int lev, Real dt, // ex_type, ey_type, ez_type, bx_type, by_type, bz_type, // dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes, // nox, galerkin_interpolation, problo, probhi, cur_time); - if (PulsarParm::AddExternalMonopoleOnly == 1) { + if (AddExternalMonopoleOnly == 1) { amrex::Real Erp, Ethetap, Ephip; amrex::Real r_p, theta_p, phi_p; - PulsarParm::ConvertCartesianToSphericalCoord( xp, yp, zp, - problo, probhi, r_p, - theta_p, phi_p); - if (r_p > PulsarParm::corotatingE_maxradius) { - PulsarParm::ExternalEMonopoleSpherical(r_p, theta_p, phi_p, cur_time, + Pulsar::ConvertCartesianToSphericalCoord( xp, yp, zp, center_star_arr, + r_p, theta_p, phi_p); + if (r_p > corotatingE_maxradius_data) { + Pulsar::ExternalEMonopoleSpherical(r_p, theta_p, phi_p, cur_time, + omega_star_data, + ramp_omega_time_data, + Bstar_data, Rstar_data, Erp, Ethetap, Ephip); amrex::Real Ex_monopole, Ey_monopole, Ez_monopole; - PulsarParm::ConvertSphericalToCartesianXComponent(Erp, Ethetap, + Pulsar::ConvertSphericalToCartesianXComponent(Erp, Ethetap, Ephip, r_p, theta_p, phi_p, Ex_monopole); - PulsarParm::ConvertSphericalToCartesianYComponent(Erp, Ethetap, + Pulsar::ConvertSphericalToCartesianYComponent(Erp, Ethetap, Ephip, r_p, theta_p, phi_p, Ey_monopole); - PulsarParm::ConvertSphericalToCartesianZComponent(Erp, Ethetap, + Pulsar::ConvertSphericalToCartesianZComponent(Erp, Ethetap, Ephip, r_p, theta_p, phi_p, Ez_monopole); Exp += Ex_monopole; Eyp += Ey_monopole; Ezp += Ez_monopole; } - } + } #endif } // Externally applied E-field in Cartesian co-ordinates @@ -2255,60 +2292,73 @@ PhysicalParticleContainer::PushP (int lev, Real dt, amrex::Real qp = q; if (ion_lev) { qp *= ion_lev[ip]; } #ifdef PULSAR -// amrex::AllPrintToFile("PulsarParticle") << "Exp " << Exp << " Eyp " << Eyp << " Ezp " << Ezp << " Bxp : " << Bxp << " Byp : " << Byp << " Bzp " << Bzp << "\n"; -// amrex::AllPrintToFile("PulsarParticle") << "uxp " << ux[ip] << " uyp " << uy[ip] << " uzp " << uz[ip] << "\n"; +// amrex::AllPrintToFile("PulsarParticle") << "Exp " << Exp << " Eyp " << Eyp << " Ezp " << Ezp << " Bxp : " << Bxp << " Byp : " << Byp << " Bzp " << Bzp << "\n"; +// amrex::AllPrintToFile("PulsarParticle") << "uxp " << ux[ip] << " uyp " << uy[ip] << " uzp " << uz[ip] << "\n"; amrex::Real Erp, Ethetap, Ephip; - PulsarParm::ExternalEFieldSpherical(r_p, theta_p, phi_p, cur_time, + Pulsar::ExternalEFieldSpherical(r_p, theta_p, phi_p, cur_time, + omega_star_data, + ramp_omega_time_data, + Bstar_data, Rstar_data, + corotatingE_maxradius_data, + E_external_monopole_data, Erp, Ethetap, Ephip); amrex::Real Er_comp, Etheta_comp, Ephi_comp; - PulsarParm::ConvertCartesianToSphericalRComponent( Exp, Eyp, Ezp, + Pulsar::ConvertCartesianToSphericalRComponent( Exp, Eyp, Ezp, theta_p, phi_p, Er_comp); - PulsarParm::ConvertCartesianToSphericalThetaComponent( Exp, Eyp, Ezp, + Pulsar::ConvertCartesianToSphericalThetaComponent( Exp, Eyp, Ezp, theta_p, phi_p, Etheta_comp); - PulsarParm::ConvertCartesianToSphericalPhiComponent( Exp, Eyp, Ezp, + Pulsar::ConvertCartesianToSphericalPhiComponent( Exp, Eyp, Ezp, theta_p, phi_p, Ephi_comp); amrex::Real Brp, Bthetap, Bphip; - PulsarParm::ExternalBFieldSpherical (r_p, theta_p, phi_p, cur_time, + Pulsar::ExternalBFieldSpherical (r_p, theta_p, phi_p, cur_time, + Bstar_data, Rstar_data, dRstar_data, Brp, Bthetap, Bphip); amrex::Real Br_comp, Btheta_comp, Bphi_comp; - PulsarParm::ConvertCartesianToSphericalRComponent( Bxp, Byp, Bzp, + Pulsar::ConvertCartesianToSphericalRComponent( Bxp, Byp, Bzp, theta_p, phi_p, Br_comp); - PulsarParm::ConvertCartesianToSphericalThetaComponent( Bxp, Byp, Bzp, + Pulsar::ConvertCartesianToSphericalThetaComponent( Bxp, Byp, Bzp, theta_p, phi_p, Btheta_comp); - PulsarParm::ConvertCartesianToSphericalPhiComponent( Bxp, Byp, Bzp, + Pulsar::ConvertCartesianToSphericalPhiComponent( Bxp, Byp, Bzp, theta_p, phi_p, Bphi_comp); -// amrex::AllPrintToFile("PulsarParticle") << " Er theory " << Erp << " Etheta_theory " << Ethetap << " Ephi_theory " << Ephip << "\n"; -// amrex::AllPrintToFile("PulsarParticle") << " Er sim " << Er_comp << " Etheta_comp " << Etheta_comp << " Ephi_comp " << Ephi_comp << "\n"; +// amrex::AllPrintToFile("PulsarParticle") << " Er theory " << Erp << " Etheta_theory " << Ethetap << " Ephi_theory " << Ephip << "\n"; +// amrex::AllPrintToFile("PulsarParticle") << " Er sim " << Er_comp << " Etheta_comp " << Etheta_comp << " Ephi_comp " << Ephi_comp << "\n"; // amrex::AllPrintToFile("PulsarParticle") << " Br theory " << Brp << " Btheta_theory " << Bthetap << " Bphi_theory " << Bphip << "\n"; // amrex::AllPrintToFile("PulsarParticle") << " Br sim " << Br_comp << " Btheta_comp " << Btheta_comp << " Bphi_comp " << Bphi_comp << "\n"; amrex::Real ur_p, utheta_p, uphi_p; - PulsarParm::ConvertCartesianToSphericalRComponent( ux[ip], uy[ip], uz[ip], + Pulsar::ConvertCartesianToSphericalRComponent( ux[ip], uy[ip], uz[ip], theta_p, phi_p, ur_p); - PulsarParm::ConvertCartesianToSphericalThetaComponent( ux[ip], uy[ip], uz[ip], + Pulsar::ConvertCartesianToSphericalThetaComponent( ux[ip], uy[ip], uz[ip], theta_p, phi_p, utheta_p); - PulsarParm::ConvertCartesianToSphericalPhiComponent( ux[ip], uy[ip], uz[ip], + Pulsar::ConvertCartesianToSphericalPhiComponent( ux[ip], uy[ip], uz[ip], theta_p, phi_p, uphi_p); amrex::Real Exp_theory, Eyp_theory, Ezp_theory; amrex::Real Bxp_theory, Byp_theory, Bzp_theory; - if (PulsarParm::use_theoreticalEB == 1) { - if (r_p < ( PulsarParm::theory_max_rstar) ) { - PulsarParm::ExternalEFieldSpherical(r_p, theta_p, phi_p, cur_time, + if (use_theoreticalEB == 1) { + if (r_p < ( theory_max_rstar) ) { + Pulsar::ExternalEFieldSpherical(r_p, theta_p, phi_p, cur_time, + omega_star_data, + ramp_omega_time_data, + Bstar_data, Rstar_data, + corotatingE_maxradius_data, + E_external_monopole_data, Erp, Ethetap, Ephip); - PulsarParm::ExternalBFieldSpherical(r_p, theta_p, phi_p, cur_time, + Pulsar::ExternalBFieldSpherical(r_p, theta_p, phi_p, cur_time, + Bstar_data, Rstar_data, + dRstar_data, Brp, Bthetap, Bphip); - PulsarParm::ConvertSphericalToCartesianXComponent(Erp, Ethetap, + Pulsar::ConvertSphericalToCartesianXComponent(Erp, Ethetap, Ephip, r_p, theta_p, phi_p, Exp_theory); - PulsarParm::ConvertSphericalToCartesianYComponent(Erp, Ethetap, Ephip, r_p, + Pulsar::ConvertSphericalToCartesianYComponent(Erp, Ethetap, Ephip, r_p, theta_p, phi_p, Eyp_theory); - PulsarParm::ConvertSphericalToCartesianZComponent(Erp, Ethetap, Ephip, r_p, + Pulsar::ConvertSphericalToCartesianZComponent(Erp, Ethetap, Ephip, r_p, theta_p, phi_p, Ezp_theory); - PulsarParm::ConvertSphericalToCartesianXComponent(Brp, Bthetap, Bphip, r_p, + Pulsar::ConvertSphericalToCartesianXComponent(Brp, Bthetap, Bphip, r_p, theta_p, phi_p, Bxp_theory); - PulsarParm::ConvertSphericalToCartesianYComponent(Brp, Bthetap, Bphip, r_p, + Pulsar::ConvertSphericalToCartesianYComponent(Brp, Bthetap, Bphip, r_p, theta_p, phi_p, Byp_theory); - PulsarParm::ConvertSphericalToCartesianZComponent(Brp, Bthetap, Bphip, r_p, + Pulsar::ConvertSphericalToCartesianZComponent(Brp, Bthetap, Bphip, r_p, theta_p, phi_p, Bzp_theory); Exp = Exp_theory; Eyp = Eyp_theory; @@ -2318,17 +2368,17 @@ PhysicalParticleContainer::PushP (int lev, Real dt, Bzp = Bzp_theory; } } - PulsarParm::ConvertCartesianToSphericalRComponent( Exp, Eyp, Ezp, - theta_p, phi_p, Er_comp); - PulsarParm::ConvertCartesianToSphericalThetaComponent( Exp, Eyp, Ezp, - theta_p, phi_p, Etheta_comp); - PulsarParm::ConvertCartesianToSphericalPhiComponent( Exp, Eyp, Ezp, - theta_p, phi_p, Ephi_comp); - PulsarParm::ConvertCartesianToSphericalRComponent( Bxp, Byp, Bzp, - theta_p, phi_p, Br_comp); - PulsarParm::ConvertCartesianToSphericalThetaComponent( Bxp, Byp, Bzp, - theta_p, phi_p, Btheta_comp); - PulsarParm::ConvertCartesianToSphericalPhiComponent( Bxp, Byp, Bzp, + Pulsar::ConvertCartesianToSphericalRComponent( Exp, Eyp, Ezp, + theta_p, phi_p, Er_comp); + Pulsar::ConvertCartesianToSphericalThetaComponent( Exp, Eyp, Ezp, + theta_p, phi_p, Etheta_comp); + Pulsar::ConvertCartesianToSphericalPhiComponent( Exp, Eyp, Ezp, + theta_p, phi_p, Ephi_comp); + Pulsar::ConvertCartesianToSphericalRComponent( Bxp, Byp, Bzp, + theta_p, phi_p, Br_comp); + Pulsar::ConvertCartesianToSphericalThetaComponent( Bxp, Byp, Bzp, + theta_p, phi_p, Btheta_comp); + Pulsar::ConvertCartesianToSphericalPhiComponent( Bxp, Byp, Bzp, theta_p, phi_p, Bphi_comp); amrex::Real uxip1 = ux[ip]; amrex::Real uyip1 = uy[ip]; @@ -2373,73 +2423,76 @@ PhysicalParticleContainer::PushP (int lev, Real dt, #endif UpdateMomentumBoris( ux[ip], uy[ip], uz[ip], Exp, Eyp, Ezp, Bxp, - Byp, Bzp, qp, m, dt, - PulsarParticleDiagData, ip, 0); + Byp, Bzp, qp, m, dt +#ifdef PULSAR + , PulsarParticleDiagData, ip, 0 +#endif + ); #ifdef PULSAR if (ip == 0) { amrex::Real Exp_theory, Eyp_theory, Ezp_theory; - PulsarParm::ConvertSphericalToCartesianXComponent(Erp, Ethetap, Ephip, r_p, - theta_p, phi_p, Exp_theory); - PulsarParm::ConvertSphericalToCartesianYComponent(Erp, Ethetap, Ephip, r_p, - theta_p, phi_p, Eyp_theory); - PulsarParm::ConvertSphericalToCartesianZComponent(Erp, Ethetap, Ephip, r_p, + Pulsar::ConvertSphericalToCartesianXComponent(Erp, Ethetap, Ephip, r_p, + theta_p, phi_p, Exp_theory); + Pulsar::ConvertSphericalToCartesianYComponent(Erp, Ethetap, Ephip, r_p, + theta_p, phi_p, Eyp_theory); + Pulsar::ConvertSphericalToCartesianZComponent(Erp, Ethetap, Ephip, r_p, theta_p, phi_p, Ezp_theory); amrex::Real Bxp_theory, Byp_theory, Bzp_theory; - PulsarParm::ConvertSphericalToCartesianXComponent(Brp, Bthetap, Bphip, r_p, - theta_p, phi_p, Bxp_theory); - PulsarParm::ConvertSphericalToCartesianYComponent(Brp, Bthetap, Bphip, r_p, - theta_p, phi_p, Byp_theory); - PulsarParm::ConvertSphericalToCartesianZComponent(Brp, Bthetap, Bphip, r_p, + Pulsar::ConvertSphericalToCartesianXComponent(Brp, Bthetap, Bphip, r_p, + theta_p, phi_p, Bxp_theory); + Pulsar::ConvertSphericalToCartesianYComponent(Brp, Bthetap, Bphip, r_p, + theta_p, phi_p, Byp_theory); + Pulsar::ConvertSphericalToCartesianZComponent(Brp, Bthetap, Bphip, r_p, theta_p, phi_p, Bzp_theory); UpdateMomentumBoris( uxip2, uyip2, uzip2, Exp_theory, Eyp_theory, Ezp_theory, Bxp_theory, Byp_theory, Bzp_theory, qp, m, dt, PulsarParticleDiagData, ip, 1); - PulsarParm::ConvertCartesianToSphericalRComponent( + Pulsar::ConvertCartesianToSphericalRComponent( PulsarParticleDiagData[31], PulsarParticleDiagData[32], PulsarParticleDiagData[33], theta_p, phi_p, PulsarParticleDiagData[37] ); - PulsarParm::ConvertCartesianToSphericalThetaComponent( + Pulsar::ConvertCartesianToSphericalThetaComponent( PulsarParticleDiagData[31], PulsarParticleDiagData[32], PulsarParticleDiagData[33], theta_p, phi_p, PulsarParticleDiagData[38] ); - PulsarParm::ConvertCartesianToSphericalPhiComponent( + Pulsar::ConvertCartesianToSphericalPhiComponent( PulsarParticleDiagData[31], PulsarParticleDiagData[32], PulsarParticleDiagData[33], theta_p, phi_p, PulsarParticleDiagData[39] ); - PulsarParm::ConvertCartesianToSphericalRComponent( + Pulsar::ConvertCartesianToSphericalRComponent( PulsarParticleDiagData[34], PulsarParticleDiagData[35], PulsarParticleDiagData[36], theta_p, phi_p, PulsarParticleDiagData[40] ); - PulsarParm::ConvertCartesianToSphericalThetaComponent( + Pulsar::ConvertCartesianToSphericalThetaComponent( PulsarParticleDiagData[34], PulsarParticleDiagData[35], PulsarParticleDiagData[36], theta_p, phi_p, PulsarParticleDiagData[41] ); - PulsarParm::ConvertCartesianToSphericalPhiComponent( + Pulsar::ConvertCartesianToSphericalPhiComponent( PulsarParticleDiagData[34], PulsarParticleDiagData[35], PulsarParticleDiagData[36], theta_p, phi_p, PulsarParticleDiagData[42] ); - PulsarParm::ConvertCartesianToSphericalRComponent( + Pulsar::ConvertCartesianToSphericalRComponent( PulsarParticleDiagData[43], PulsarParticleDiagData[44], PulsarParticleDiagData[45], theta_p, phi_p, PulsarParticleDiagData[49] ); - PulsarParm::ConvertCartesianToSphericalThetaComponent( + Pulsar::ConvertCartesianToSphericalThetaComponent( PulsarParticleDiagData[43], PulsarParticleDiagData[44], PulsarParticleDiagData[45], theta_p, phi_p, PulsarParticleDiagData[50] ); - PulsarParm::ConvertCartesianToSphericalPhiComponent( + Pulsar::ConvertCartesianToSphericalPhiComponent( PulsarParticleDiagData[43], PulsarParticleDiagData[44], PulsarParticleDiagData[45], theta_p, phi_p, PulsarParticleDiagData[51] ); - PulsarParm::ConvertCartesianToSphericalRComponent( + Pulsar::ConvertCartesianToSphericalRComponent( PulsarParticleDiagData[46], PulsarParticleDiagData[47], PulsarParticleDiagData[48], theta_p, phi_p, PulsarParticleDiagData[52] ); - PulsarParm::ConvertCartesianToSphericalThetaComponent( + Pulsar::ConvertCartesianToSphericalThetaComponent( PulsarParticleDiagData[46], PulsarParticleDiagData[47], PulsarParticleDiagData[48], theta_p, phi_p, PulsarParticleDiagData[53] ); - PulsarParm::ConvertCartesianToSphericalPhiComponent( + Pulsar::ConvertCartesianToSphericalPhiComponent( PulsarParticleDiagData[46], PulsarParticleDiagData[47], PulsarParticleDiagData[48], theta_p, phi_p, PulsarParticleDiagData[54] ); @@ -2465,7 +2518,8 @@ PhysicalParticleContainer::PushP (int lev, Real dt, } // amrex::AllPrintToFile("PulsarParticle") << " momentum update complete step : " << warpx.getistep(0) << " time : " << cur_time << "\n"; amrex::Gpu::synchronize(); - if (PulsarParm::singleParticleTest == 1) { +#ifdef PULSAR + if (Pulsar::m_singleParticleTest == 1) { const amrex::Real q = this->charge; if (q > 0) { amrex::AllPrintToFile("PulsarPositronDiagnostics") << " cur_time xp yp zp r_p theta_p phi_p ux uy uz ur utheta uphi Ex Ey Ez Bx By Bz Er Etheta Ephi Br Btheta Bphi Er_theory Etheta_theory Ephi_theory Br_theory Btheta_theory Bphi_theory qEx_force qEy_force qEz_force qvcrossB_x qvcrossB_y qvcrossB_z qE_r qE_theta qE_phi qvcrossB_r qvcrossB_theta qvcrossB_phi qEx_theory qEy_theory qEz_theory qvcrossB_x_theory qvcrossB_y_theory qvcrossB_z_theory qEr_theory qEtheta_theory qEphi_theory qvcrossB_r_theory qvcrossB_theta_theory qvcrossB_phi_theory\n"; @@ -2473,6 +2527,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, amrex::AllPrintToFile("PulsarElectronDiagnostics") << " cur_time xp yp zp r_p theta_p phi_p ux uy uz ur utheta uphi Ex Ey Ez Bx By Bz Er Etheta Ephi Br Btheta Bphi Er_theory Etheta_theory Ephi_theory Br_theory Btheta_theory Bphi_theory qEx_force qEy_force qEz_force qvcrossB_x qvcrossB_y qvcrossB_z qE_r qE_theta qE_phi qvcrossB_r qvcrossB_theta qvcrossB_phi qEx_theory qEy_theory qEz_theory qvcrossB_x_theory qvcrossB_y_theory qvcrossB_z_theory qEr_theory qEtheta_theory qEphi_theory qvcrossB_r_theory qvcrossB_theta_theory qvcrossB_phi_theory\n"; } } +#endif } void @@ -2738,9 +2793,20 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, const Dim3 lo = lbound(box); #ifdef PULSAR - const auto problo = WarpX::GetInstance().Geom(lev).ProbLoArray(); - const auto probhi = WarpX::GetInstance().Geom(lev).ProbHiArray(); auto &warpx = WarpX::GetInstance(); + const auto problo = warpx.Geom(lev).ProbLoArray(); + const auto probhi = warpx.Geom(lev).ProbHiArray(); + amrex::Real omega_star_data = Pulsar::m_omega_star; + amrex::Real ramp_omega_time_data = Pulsar::m_omega_ramp_time; + amrex::Real Bstar_data = Pulsar::m_B_star; + amrex::Real Rstar_data = Pulsar::m_R_star; + amrex::Real dRstar_data = Pulsar::m_dR_star; + amrex::Real corotatingE_maxradius_data = Pulsar::m_corotatingE_maxradius; + int E_external_monopole_data = Pulsar::m_do_E_external_monopole; + int AddExternalMonopoleOnly = Pulsar::m_AddExternalMonopoleOnly; + int use_theoreticalEB = Pulsar::m_use_theoreticalEB; + amrex::Real theory_max_rstar = Pulsar::m_theory_max_rstar; + int singleParticleTest = Pulsar::m_singleParticleTest; amrex::Gpu::synchronize(); amrex::Gpu::ManagedVector PulsarParticleDiag(55,0.0); amrex::Real * PulsarParticleDiagData = PulsarParticleDiag.data(); @@ -2820,6 +2886,10 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, #ifdef PULSAR auto& particles = pti.GetArrayOfStructs(); + amrex::GpuArray center_star_arr; + for (int idim = 0; idim < 3; ++idim) { + center_star_arr[idim] = Pulsar::m_center_star[idim]; + } #endif amrex::ParallelFor( np_to_push, [=] AMREX_GPU_DEVICE (long ip) @@ -2829,7 +2899,8 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, #ifdef PULSAR amrex::Real r_p, theta_p, phi_p; - PulsarParm::ConvertCartesianToSphericalCoord( xp, yp, zp, problo, probhi, r_p, theta_p, phi_p); + Pulsar::ConvertCartesianToSphericalCoord( xp, yp, zp, center_star_arr, + r_p, theta_p, phi_p); #endif if (save_previous_position) { x_old[ip] = xp; @@ -2855,23 +2926,25 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, // ex_type, ey_type, ez_type, bx_type, by_type, bz_type, // dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes, // nox, galerkin_interpolation, problo, probhi, cur_time); - if (PulsarParm::AddExternalMonopoleOnly == 1) { + if (AddExternalMonopoleOnly == 1) { amrex::Real Erp, Ethetap, Ephip; amrex::Real r_p, theta_p, phi_p; - PulsarParm::ConvertCartesianToSphericalCoord( xp, yp, zp, - problo, probhi, r_p, - theta_p, phi_p); - if (r_p > PulsarParm::corotatingE_maxradius) { - PulsarParm::ExternalEMonopoleSpherical(r_p, theta_p, phi_p, cur_time, + Pulsar::ConvertCartesianToSphericalCoord( xp, yp, zp, center_star_arr, + r_p, theta_p, phi_p); + if (r_p > corotatingE_maxradius_data) { + Pulsar::ExternalEMonopoleSpherical(r_p, theta_p, phi_p, cur_time, + omega_star_data, + ramp_omega_time_data, + Bstar_data, Rstar_data, Erp, Ethetap, Ephip); amrex::Real Ex_monopole, Ey_monopole, Ez_monopole; - PulsarParm::ConvertSphericalToCartesianXComponent(Erp, Ethetap, + Pulsar::ConvertSphericalToCartesianXComponent(Erp, Ethetap, Ephip, r_p, theta_p, phi_p, Ex_monopole); - PulsarParm::ConvertSphericalToCartesianYComponent(Erp, Ethetap, + Pulsar::ConvertSphericalToCartesianYComponent(Erp, Ethetap, Ephip, r_p, theta_p, phi_p, Ey_monopole); - PulsarParm::ConvertSphericalToCartesianZComponent(Erp, Ethetap, + Pulsar::ConvertSphericalToCartesianZComponent(Erp, Ethetap, Ephip, r_p, theta_p, phi_p, Ez_monopole); Exp += Ex_monopole; @@ -2889,13 +2962,19 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, scaleFields(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp); #ifdef PULSAR - + amrex::Real Erp, Ethetap, Ephip; - PulsarParm::ExternalEFieldSpherical(r_p, theta_p, phi_p, cur_time, - Erp, Ethetap, Ephip); + Pulsar::ExternalEFieldSpherical(r_p, theta_p, phi_p, cur_time, + omega_star_data, + ramp_omega_time_data, + Bstar_data, Rstar_data, + corotatingE_maxradius_data, + E_external_monopole_data, + Erp, Ethetap, Ephip); amrex::Real Brp, Bthetap, Bphip; - PulsarParm::ExternalBFieldSpherical (r_p, theta_p, phi_p, cur_time, - Brp, Bthetap, Bphip); + Pulsar::ExternalBFieldSpherical (r_p, theta_p, phi_p, cur_time, + Bstar_data, Rstar_data, dRstar_data, + Brp, Bthetap, Bphip); //amrex::Real Er_comp, Etheta_comp, Ephi_comp; //PulsarParm::ConvertCartesianToSphericalRComponent( Exp, Eyp, Ezp, // theta_p, phi_p, Er_comp); @@ -2919,25 +2998,31 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, // theta_p, phi_p, uphi_p); amrex::Real Exp_theory, Eyp_theory, Ezp_theory; amrex::Real Bxp_theory, Byp_theory, Bzp_theory; - if (PulsarParm::use_theoreticalEB == 1) { - if (r_p < PulsarParm::theory_max_rstar) { - PulsarParm::ExternalEFieldSpherical(r_p, theta_p, phi_p, cur_time, - Erp, Ethetap, Ephip); - PulsarParm::ConvertSphericalToCartesianXComponent(Erp, Ethetap, - Ephip, r_p, theta_p, - phi_p, Exp_theory); - PulsarParm::ConvertSphericalToCartesianYComponent(Erp, Ethetap, Ephip, r_p, - theta_p, phi_p, Eyp_theory); - PulsarParm::ConvertSphericalToCartesianZComponent(Erp, Ethetap, Ephip, r_p, - theta_p, phi_p, Ezp_theory); - PulsarParm::ExternalBFieldSpherical(r_p, theta_p, phi_p, cur_time, - Brp, Bthetap, Bphip); - PulsarParm::ConvertSphericalToCartesianXComponent(Brp, Bthetap, Bphip, r_p, - theta_p, phi_p, Bxp_theory); - PulsarParm::ConvertSphericalToCartesianYComponent(Brp, Bthetap, Bphip, r_p, - theta_p, phi_p, Byp_theory); - PulsarParm::ConvertSphericalToCartesianZComponent(Brp, Bthetap, Bphip, r_p, - theta_p, phi_p, Bzp_theory); + if (use_theoreticalEB == 1) { + if (r_p < theory_max_rstar) { + Pulsar::ExternalEFieldSpherical(r_p, theta_p, phi_p, cur_time, + omega_star_data, + ramp_omega_time_data, + Bstar_data, Rstar_data, + corotatingE_maxradius_data, + E_external_monopole_data, + Erp, Ethetap, Ephip); + Pulsar::ConvertSphericalToCartesianXComponent(Erp, Ethetap, + Ephip, r_p, theta_p, + phi_p, Exp_theory); + Pulsar::ConvertSphericalToCartesianYComponent(Erp, Ethetap, Ephip, r_p, + theta_p, phi_p, Eyp_theory); + Pulsar::ConvertSphericalToCartesianZComponent(Erp, Ethetap, Ephip, r_p, + theta_p, phi_p, Ezp_theory); + Pulsar::ExternalBFieldSpherical(r_p, theta_p, phi_p, cur_time, + Bstar_data, Rstar_data, dRstar_data, + Brp, Bthetap, Bphip); + Pulsar::ConvertSphericalToCartesianXComponent(Brp, Bthetap, Bphip, r_p, + theta_p, phi_p, Bxp_theory); + Pulsar::ConvertSphericalToCartesianYComponent(Brp, Bthetap, Bphip, r_p, + theta_p, phi_p, Byp_theory); + Pulsar::ConvertSphericalToCartesianZComponent(Brp, Bthetap, Bphip, r_p, + theta_p, phi_p, Bzp_theory); Exp = Exp_theory; Eyp = Eyp_theory; Ezp = Ezp_theory; @@ -2947,25 +3032,25 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, } } amrex::Real Er_comp, Etheta_comp, Ephi_comp; - PulsarParm::ConvertCartesianToSphericalRComponent( Exp, Eyp, Ezp, - theta_p, phi_p, Er_comp); - PulsarParm::ConvertCartesianToSphericalThetaComponent( Exp, Eyp, Ezp, - theta_p, phi_p, Etheta_comp); - PulsarParm::ConvertCartesianToSphericalPhiComponent( Exp, Eyp, Ezp, - theta_p, phi_p, Ephi_comp); + Pulsar::ConvertCartesianToSphericalRComponent( Exp, Eyp, Ezp, + theta_p, phi_p, Er_comp); + Pulsar::ConvertCartesianToSphericalThetaComponent( Exp, Eyp, Ezp, + theta_p, phi_p, Etheta_comp); + Pulsar::ConvertCartesianToSphericalPhiComponent( Exp, Eyp, Ezp, + theta_p, phi_p, Ephi_comp); amrex::Real Br_comp, Btheta_comp, Bphi_comp; - PulsarParm::ConvertCartesianToSphericalRComponent( Bxp, Byp, Bzp, - theta_p, phi_p, Br_comp); - PulsarParm::ConvertCartesianToSphericalThetaComponent( Bxp, Byp, Bzp, - theta_p, phi_p, Btheta_comp); - PulsarParm::ConvertCartesianToSphericalPhiComponent( Bxp, Byp, Bzp, - theta_p, phi_p, Bphi_comp); + Pulsar::ConvertCartesianToSphericalRComponent( Bxp, Byp, Bzp, + theta_p, phi_p, Br_comp); + Pulsar::ConvertCartesianToSphericalThetaComponent( Bxp, Byp, Bzp, + theta_p, phi_p, Btheta_comp); + Pulsar::ConvertCartesianToSphericalPhiComponent( Bxp, Byp, Bzp, + theta_p, phi_p, Bphi_comp); amrex::Real ur_p, utheta_p, uphi_p; - PulsarParm::ConvertCartesianToSphericalRComponent( ux[ip+offset], uy[ip+offset], uz[ip+offset], - theta_p, phi_p, ur_p); - PulsarParm::ConvertCartesianToSphericalThetaComponent( ux[ip+offset], uy[ip+offset], uz[ip+offset], - theta_p, phi_p, utheta_p); - PulsarParm::ConvertCartesianToSphericalPhiComponent( ux[ip+offset], uy[ip+offset], uz[ip+offset], + Pulsar::ConvertCartesianToSphericalRComponent( ux[ip+offset], uy[ip+offset], uz[ip+offset], + theta_p, phi_p, ur_p); + Pulsar::ConvertCartesianToSphericalThetaComponent( ux[ip+offset], uy[ip+offset], uz[ip+offset], + theta_p, phi_p, utheta_p); + Pulsar::ConvertCartesianToSphericalPhiComponent( ux[ip+offset], uy[ip+offset], uz[ip+offset], theta_p, phi_p, uphi_p); if (ip == 0) { PulsarParticleDiagData[0] = cur_time; @@ -3054,56 +3139,56 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, // Byp_theory); // PulsarParm::ConvertSphericalToCartesianZComponent(Brp, Bthetap, Bphip, r_p, theta_p, phi_p, // Bzp_theory); - PulsarParm::ConvertCartesianToSphericalRComponent( + Pulsar::ConvertCartesianToSphericalRComponent( PulsarParticleDiagData[31], PulsarParticleDiagData[32], PulsarParticleDiagData[33], theta_p, phi_p, PulsarParticleDiagData[37] ); - PulsarParm::ConvertCartesianToSphericalThetaComponent( + Pulsar::ConvertCartesianToSphericalThetaComponent( PulsarParticleDiagData[31], PulsarParticleDiagData[32], PulsarParticleDiagData[33], theta_p, phi_p, PulsarParticleDiagData[38] ); - PulsarParm::ConvertCartesianToSphericalPhiComponent( + Pulsar::ConvertCartesianToSphericalPhiComponent( PulsarParticleDiagData[31], PulsarParticleDiagData[32], PulsarParticleDiagData[33], theta_p, phi_p, PulsarParticleDiagData[39] ); - PulsarParm::ConvertCartesianToSphericalRComponent( + Pulsar::ConvertCartesianToSphericalRComponent( PulsarParticleDiagData[34], PulsarParticleDiagData[35], PulsarParticleDiagData[36], theta_p, phi_p, PulsarParticleDiagData[40] ); - PulsarParm::ConvertCartesianToSphericalThetaComponent( + Pulsar::ConvertCartesianToSphericalThetaComponent( PulsarParticleDiagData[34], PulsarParticleDiagData[35], PulsarParticleDiagData[36], theta_p, phi_p, PulsarParticleDiagData[41] ); - PulsarParm::ConvertCartesianToSphericalPhiComponent( + Pulsar::ConvertCartesianToSphericalPhiComponent( PulsarParticleDiagData[34], PulsarParticleDiagData[35], PulsarParticleDiagData[36], theta_p, phi_p, PulsarParticleDiagData[42] ); // theoretical forces - PulsarParm::ConvertCartesianToSphericalRComponent( + Pulsar::ConvertCartesianToSphericalRComponent( PulsarParticleDiagData[43], PulsarParticleDiagData[44], PulsarParticleDiagData[45], theta_p, phi_p, PulsarParticleDiagData[49] ); - PulsarParm::ConvertCartesianToSphericalThetaComponent( + Pulsar::ConvertCartesianToSphericalThetaComponent( PulsarParticleDiagData[43], PulsarParticleDiagData[44], PulsarParticleDiagData[45], theta_p, phi_p, PulsarParticleDiagData[50] ); - PulsarParm::ConvertCartesianToSphericalPhiComponent( + Pulsar::ConvertCartesianToSphericalPhiComponent( PulsarParticleDiagData[43], PulsarParticleDiagData[44], PulsarParticleDiagData[45], theta_p, phi_p, PulsarParticleDiagData[51] ); - PulsarParm::ConvertCartesianToSphericalRComponent( + Pulsar::ConvertCartesianToSphericalRComponent( PulsarParticleDiagData[46], PulsarParticleDiagData[47], PulsarParticleDiagData[48], theta_p, phi_p, PulsarParticleDiagData[52] ); - PulsarParm::ConvertCartesianToSphericalThetaComponent( + Pulsar::ConvertCartesianToSphericalThetaComponent( PulsarParticleDiagData[46], PulsarParticleDiagData[47], PulsarParticleDiagData[48], theta_p, phi_p, PulsarParticleDiagData[53] ); - PulsarParm::ConvertCartesianToSphericalPhiComponent( + Pulsar::ConvertCartesianToSphericalPhiComponent( PulsarParticleDiagData[46], PulsarParticleDiagData[47], PulsarParticleDiagData[48], theta_p, phi_p, PulsarParticleDiagData[54] ); - } + } #endif //#ifdef PULSAR @@ -3123,7 +3208,8 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, }); // amrex::AllPrintToFile("PulsarParticle") << " done with particle push " << cur_time << "\n"; amrex::Gpu::synchronize(); - if (PulsarParm::singleParticleTest == 1) { +#ifdef PULSAR + if (singleParticleTest == 1) { if ( q > 0) { // positrons amrex::AllPrintToFile("PulsarPositronDiagnostics") << " "; for (int i = 0; i < PulsarParticleDiag.size(); ++i) { @@ -3139,6 +3225,7 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, } amrex::Gpu::synchronize(); } +#endif } void @@ -3305,8 +3392,7 @@ PhysicalParticleContainer::getPairGenerationFilterFunc () #ifdef PULSAR void PhysicalParticleContainer::PulsarParticleInjection() { - - if (PulsarParm::singleParticleTest == 1) { + if (Pulsar::m_singleParticleTest == 1) { AddParticles(0); // Note - add on level 0 Redistribute(); // We then redistribute } else { @@ -3316,10 +3402,15 @@ void PhysicalParticleContainer::PulsarParticleInjection() { void PhysicalParticleContainer::PulsarParticleRemoval() { int lev = 0; - Gpu::DeviceScalar sumParticles(0); - Gpu::DeviceScalar sumWeight(0.0); - int sum_d ; + Gpu::DeviceScalar sumParticles(0); + Gpu::DeviceScalar sumWeight(0.0); + int sum_d ; const amrex::Real q = this->charge; + amrex::GpuArray center_star_arr; + for (int i_dim = 0; i_dim < 3; ++i_dim) { + center_star_arr[i_dim] = Pulsar::m_center_star[i_dim]; + } + amrex::Real max_particle_absorption_radius = Pulsar::m_max_particle_absorption_radius; // Remove Particles From inside sphere #ifdef _OPENMP #pragma omp parallel @@ -3333,9 +3424,9 @@ void PhysicalParticleContainer::PulsarParticleRemoval() { for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { const auto GetPosition = GetParticlePosition(pti); - Real xc = PulsarParm::center_star[0]; - Real yc = PulsarParm::center_star[1]; - Real zc = PulsarParm::center_star[2]; + Real xc = center_star_arr[0]; + Real yc = center_star_arr[1]; + Real zc = center_star_arr[2]; ParticleType* pp = pti.GetArrayOfStructs()().data(); auto& attribs = pti.GetAttribs(); auto& wp = attribs[PIdx::w]; @@ -3350,7 +3441,7 @@ void PhysicalParticleContainer::PulsarParticleRemoval() { Real r = std::sqrt((x-xc)*(x-xc) + (y-yc)*(y-yc) + (z-zc)*(z-zc)); - if (r <= (PulsarParm::max_particle_absorption_radius)) { + if (r <= (max_particle_absorption_radius)) { pp[i].id() = -1; // atomic add int const unity = 1; diff --git a/Source/Particles/PulsarParameters.H b/Source/Particles/PulsarParameters.H index dcc10855034..fe178ba7949 100644 --- a/Source/Particles/PulsarParameters.H +++ b/Source/Particles/PulsarParameters.H @@ -11,122 +11,83 @@ using namespace amrex; -namespace PulsarParm -{ - extern std::string pulsar_type; - - extern AMREX_GPU_DEVICE_MANAGED amrex::Real omega_star; - extern AMREX_GPU_DEVICE_MANAGED amrex::Real ramp_omega_time; - extern AMREX_GPU_DEVICE_MANAGED amrex::Real B_star; - extern AMREX_GPU_DEVICE_MANAGED amrex::Real R_star; - extern AMREX_GPU_DEVICE_MANAGED amrex::Real dR_star; - extern AMREX_GPU_DEVICE_MANAGED amrex::Real damping_scale; - extern AMREX_GPU_DEVICE_MANAGED int EB_external; - extern AMREX_GPU_DEVICE_MANAGED int E_external_monopole; - extern AMREX_GPU_DEVICE_MANAGED - amrex::GpuArray center_star; - extern AMREX_GPU_DEVICE_MANAGED amrex::Real max_ndens; - extern AMREX_GPU_DEVICE_MANAGED amrex::Real Ninj_fraction; - extern AMREX_GPU_DEVICE_MANAGED int ModifyParticleWtAtInjection; - extern AMREX_GPU_DEVICE_MANAGED amrex::Real rhoGJ_scale; - extern AMREX_GPU_DEVICE_MANAGED int damp_EB_internal; - extern AMREX_GPU_DEVICE_MANAGED amrex::Real max_EBcorotating_radius; - extern AMREX_GPU_DEVICE_MANAGED amrex::Real max_EBdamping_radius; - extern AMREX_GPU_DEVICE_MANAGED int turnoffdeposition; - extern AMREX_GPU_DEVICE_MANAGED amrex::Real max_nodepos_radius; - extern AMREX_GPU_DEVICE_MANAGED int turnoff_plasmaEB_gather; - extern AMREX_GPU_DEVICE_MANAGED amrex::Real max_nogather_radius; - extern AMREX_GPU_DEVICE_MANAGED int verbose; - extern AMREX_GPU_DEVICE_MANAGED amrex::Real max_particle_absorption_radius; - extern AMREX_GPU_DEVICE_MANAGED amrex::Real particle_inject_rmin; - extern AMREX_GPU_DEVICE_MANAGED amrex::Real particle_inject_rmax; - extern AMREX_GPU_DEVICE_MANAGED amrex::Real corotatingE_maxradius; - extern AMREX_GPU_DEVICE_MANAGED amrex::Real enforceDipoleB_maxradius; - extern AMREX_GPU_DEVICE_MANAGED amrex::Real InitializeGrid_with_Pulsar_Bfield; - extern AMREX_GPU_DEVICE_MANAGED amrex::Real InitializeGrid_with_Pulsar_Efield; - extern AMREX_GPU_DEVICE_MANAGED int enforceCorotatingE; - extern AMREX_GPU_DEVICE_MANAGED int enforceDipoleB; - extern AMREX_GPU_DEVICE_MANAGED int singleParticleTest; - extern AMREX_GPU_DEVICE_MANAGED amrex::Real Bdamping_scale; - extern AMREX_GPU_DEVICE_MANAGED int DampBDipoleInRing; - extern AMREX_GPU_DEVICE_MANAGED amrex::Real injection_time; - extern AMREX_GPU_DEVICE_MANAGED int continuous_injection; - extern AMREX_GPU_DEVICE_MANAGED amrex::Real removeparticle_theta_min; - extern AMREX_GPU_DEVICE_MANAGED amrex::Real removeparticle_theta_max; - extern AMREX_GPU_DEVICE_MANAGED int use_theoreticalEB; - extern AMREX_GPU_DEVICE_MANAGED amrex::Real theory_max_rstar; - extern AMREX_GPU_DEVICE_MANAGED int LimitDipoleBInit; - extern AMREX_GPU_DEVICE_MANAGED amrex::Real DipoleB_init_maxradius; - extern AMREX_GPU_DEVICE_MANAGED int AddExternalMonopoleOnly; - extern AMREX_GPU_DEVICE_MANAGED int AddMonopoleInsideRstarOnGrid; - extern AMREX_GPU_DEVICE_MANAGED int EnforceTheoreticalEBInGrid; - - void ReadParameters(); - void InitializeExternalPulsarFieldsOnGrid (amrex::MultiFab *mfx, + +class Pulsar { + +public: + //Constructor + Pulsar (); + + static void ReadParameters(); + static void InitializeExternalPulsarFieldsOnGrid (amrex::MultiFab *mfx, amrex::MultiFab *mfy, amrex::MultiFab *mfz, const int lev, const bool init_Bfield); - void ApplyCorotatingEfield_BC ( std::array< std::unique_ptr, 3> &Efield, + static void ApplyCorotatingEfield_BC ( std::array< std::unique_ptr, 3> &Efield, const int lev, const amrex::Real a_dt); - void ApplyDipoleBfield_BC ( std::array< std::unique_ptr, 3> &Bfield, + static void ApplyDipoleBfield_BC ( std::array< std::unique_ptr, 3> &Bfield, const int lev, const amrex::Real a_dt); AMREX_GPU_HOST_DEVICE AMREX_INLINE - amrex::Real Omega(const amrex::Real time) + static amrex::Real Omega(const amrex::Real omegastar, const amrex::Real cur_time, amrex::Real omega_ramp_time) { - amrex::Real omega = omega_star; - if (ramp_omega_time > 0.0 && time < ramp_omega_time) { - omega = omega_star * time / ramp_omega_time; + amrex::Real omega = omegastar; + if (omega_ramp_time > 0.0 && cur_time < omega_ramp_time) { + omega = omegastar * cur_time / omega_ramp_time; } return omega; } - AMREX_GPU_HOST_DEVICE AMREX_INLINE - void CorotatingEfieldSpherical (amrex::Real const r, amrex::Real const theta, + static void CorotatingEfieldSpherical (amrex::Real const r, amrex::Real const theta, amrex::Real const phi, amrex::Real const time, + amrex::Real const omega_star_data, + amrex::Real const ramp_omega_time_data, + amrex::Real const Bstar, + amrex::Real const Rstar, + amrex::Real const dRstar, amrex::Real &Er, amrex::Real &Etheta, amrex::Real &Ephi) { - amrex::Real omega = Omega(time); + amrex::ignore_unused(phi); + amrex::Real omega = Omega(omega_star_data, time, ramp_omega_time_data); amrex::Real c_theta = std::cos(theta); amrex::Real s_theta = std::sin(theta); - amrex::Real c_phi = std::cos(phi); - amrex::Real s_phi = std::sin(phi); amrex::Real r_ratio; if (r > 0) { - r_ratio = R_star/r; + r_ratio = Rstar/r; } else { - r_ratio = R_star/(dR_star*0.5); - } + r_ratio = Rstar/(dRstar*0.5); + } amrex::Real r2 = r_ratio * r_ratio; // Michel and Li -- eq 14 , 15 - Er = B_star * omega * r2 * R_star * s_theta * s_theta; - if (PulsarParm::AddMonopoleInsideRstarOnGrid == 1) { - Er += (2.0/3.0) * omega * B_star * R_star * r_ratio * r_ratio; - } - Etheta = -B_star * omega * r2 * R_star * 2.0 * s_theta * c_theta; + Er = Bstar * omega * r2 * Rstar * s_theta * s_theta; + Etheta = -Bstar * omega * r2 * Rstar * 2.0 * s_theta * c_theta; Ephi = 0.0; // aligned magnetic and rotation axis } AMREX_GPU_HOST_DEVICE AMREX_INLINE - void ExternalEFieldSpherical (amrex::Real const r, amrex::Real const theta, + static void ExternalEFieldSpherical (amrex::Real const r, amrex::Real const theta, amrex::Real const phi, amrex::Real const time, + amrex::Real const omega_star_data, + amrex::Real const ramp_omega_time_data, + amrex::Real const Bstar, amrex::Real const Rstar, + amrex::Real const corotatingE_maxradius, + int const Eexternal_monopole, amrex::Real &Er, amrex::Real &Etheta, amrex::Real &Ephi) { - amrex::Real omega = Omega(time); + amrex::ignore_unused(phi); + amrex::Real omega = Omega(omega_star_data, time, ramp_omega_time_data); amrex::Real c_theta = std::cos(theta); amrex::Real s_theta = std::sin(theta); - amrex::Real c_phi = std::cos(phi); - amrex::Real s_phi = std::sin(phi); - amrex::Real r_ratio = R_star/r; + amrex::Real r_ratio = Rstar/r; // inside pulsar //if ( r <= max_EBcorotating_radius ) { if (r <= corotatingE_maxradius) { amrex::Real r2 = r_ratio * r_ratio; // Michel and Li -- eq 14 , 15 - Er = B_star * omega * r2 * R_star * s_theta * s_theta; - Etheta = -B_star * omega * r2 * R_star * 2.0 * s_theta * c_theta; + Er = Bstar * omega * r2 * Rstar * s_theta * s_theta; + Etheta = -Bstar * omega * r2 * Rstar * 2.0 * s_theta * c_theta; Ephi = 0.0; // aligned magnetic and rotation axis } @@ -135,85 +96,62 @@ namespace PulsarParm if (r > corotatingE_maxradius) { amrex::Real r4 = r_ratio*r_ratio*r_ratio*r_ratio; // Taking derivative of phi given in eq 30 of Michel and Li - Er = B_star * omega * R_star * r4 * (1.0-3.0*c_theta*c_theta); - if (E_external_monopole == 1) { - Er += (2.0/3.0) * omega * B_star * R_star * r_ratio * r_ratio; + Er = Bstar * omega * Rstar * r4 * (1.0-3.0*c_theta*c_theta); + if (Eexternal_monopole == 1) { + Er += (2.0/3.0) * omega * Bstar * Rstar * r_ratio * r_ratio; } - Etheta = (-1.0) * B_star * omega * R_star * r4 * (2.0*s_theta*c_theta); + Etheta = (-1.0) * Bstar * omega * Rstar * r4 * (2.0*s_theta*c_theta); Ephi = 0.0; } } AMREX_GPU_HOST_DEVICE AMREX_INLINE - void ExternalEMonopoleSpherical (amrex::Real const r, amrex::Real const theta, + static void ExternalEMonopoleSpherical (amrex::Real const r, amrex::Real const theta, amrex::Real const phi, amrex::Real const time, + amrex::Real const omega_star_data, + amrex::Real const ramp_omega_time_data, + amrex::Real const Bstar, amrex::Real const Rstar, amrex::Real &Er, amrex::Real &Etheta, amrex::Real &Ephi) { - amrex::Real omega = Omega(time); - amrex::Real c_theta = std::cos(theta); - amrex::Real s_theta = std::sin(theta); - amrex::Real c_phi = std::cos(phi); - amrex::Real s_phi = std::sin(phi); - amrex::Real r_ratio = R_star/r; + amrex::ignore_unused(phi, theta); + amrex::Real omega = Omega(omega_star_data, time, ramp_omega_time_data); + amrex::Real r_ratio = Rstar/r; - Er = (2.0/3.0) * omega * B_star * R_star * r_ratio * r_ratio; + Er = (2.0/3.0) * omega * Bstar * Rstar * r_ratio * r_ratio; Etheta = 0.; Ephi = 0.; } AMREX_GPU_HOST_DEVICE AMREX_INLINE - void ExternalBFieldSpherical (amrex::Real const r, amrex::Real const theta, + static void ExternalBFieldSpherical (amrex::Real const r, amrex::Real const theta, amrex::Real const phi, amrex::Real const time, + amrex::Real const Bstar, amrex::Real const Rstar, + amrex::Real const dRstar, amrex::Real &Br, amrex::Real &Btheta, amrex::Real &Bphi) { - amrex::Real omega = Omega(time); + amrex::ignore_unused(phi, time); amrex::Real c_theta = std::cos(theta); amrex::Real s_theta = std::sin(theta); amrex::Real r_ratio; if (r > 0) { - r_ratio = R_star/r; + r_ratio = Rstar/r; } else { - r_ratio = R_star/(dR_star*0.5); + r_ratio = Rstar/(dRstar*0.5); } amrex::Real r3 = r_ratio*r_ratio*r_ratio; // Michel and Li -- eq 14 and 15 from michel and li // dipole B field inside and outside the pulsar - Br = 2.0*B_star*r3*c_theta; - Btheta = B_star*r3*s_theta; + Br = 2.0*Bstar*r3*c_theta; + Btheta = Bstar*r3*s_theta; Bphi = 0.0; } - namespace Spherical - { - AMREX_GPU_HOST_DEVICE AMREX_INLINE - amrex::Real r(int i, int j, int k, amrex::GeometryData const& geom, amrex::GpuArray const mf_ixType) - { - const auto domain_box = geom.Domain(); - const auto domain_ilo = amrex::lbound(domain_box); - const auto domain_xlo = geom.ProbLo(); - const auto domain_xhi = geom.ProbHi(); - const auto domain_dx = geom.CellSize(); - - const amrex::Real x = domain_xlo[0] + (i ) * domain_dx[0] + (1.0 - mf_ixType[0])*domain_dx[0]*0.5; - const amrex::Real y = domain_xlo[1] + (j ) * domain_dx[1] + (1.0 - mf_ixType[1])*domain_dx[1]*0.5; - const amrex::Real z = domain_xlo[2] + (k ) * domain_dx[2] + (1.0 - mf_ixType[2])*domain_dx[2]*0.5; - - const amrex::Real xc = PulsarParm::center_star[0]; - const amrex::Real yc = PulsarParm::center_star[1]; - const amrex::Real zc = PulsarParm::center_star[2]; - - const amrex::Real r = std::sqrt((x-xc)*(x-xc) + (y-yc)*(y-yc) + (z-zc)*(z-zc)); - - return r; - } - } - /** Compute Cartesian components corresponding to i, j, k based on the staggering, ixType, and the domain_lo and cell size, dx. */ AMREX_GPU_HOST_DEVICE AMREX_INLINE - void ComputeCellCoordinates ( int i, int j, int k, + static void ComputeCellCoordinates ( int i, int j, int k, amrex::GpuArray const mf_type, amrex::GpuArray const domain_lo, amrex::GpuArray const dx, @@ -225,15 +163,14 @@ namespace PulsarParm } AMREX_GPU_HOST_DEVICE AMREX_INLINE - void ConvertCartesianToSphericalCoord ( amrex::Real const x, amrex::Real const y, + static void ConvertCartesianToSphericalCoord ( amrex::Real const x, amrex::Real const y, amrex::Real const z, - amrex::GpuArray const domain_lo, - amrex::GpuArray const domain_hi, + amrex::GpuArray const center_star, amrex::Real &r, amrex::Real &theta, amrex::Real &phi) { - amrex::Real xc = PulsarParm::center_star[0]; - amrex::Real yc = PulsarParm::center_star[1]; - amrex::Real zc = PulsarParm::center_star[2]; + amrex::Real xc = center_star[0]; + amrex::Real yc = center_star[1]; + amrex::Real zc = center_star[2]; r = std::sqrt( (x-xc)*(x-xc) + (y-yc)*(y-yc) + (z-zc)*(z-zc) ); theta = 0.0; @@ -242,38 +179,41 @@ namespace PulsarParm } AMREX_GPU_HOST_DEVICE AMREX_INLINE - void ConvertSphericalToCartesianXComponent (amrex::Real const F_r, + static void ConvertSphericalToCartesianXComponent (amrex::Real const F_r, amrex::Real const F_theta, amrex::Real const F_phi, amrex::Real const r, amrex::Real const theta, amrex::Real const phi, amrex::Real & F_x) { + amrex::ignore_unused(r); F_x = F_r * std::sin(theta) * std::cos(phi) + F_theta * std::cos(theta) * std::cos(phi) - F_phi * std::sin(phi); } AMREX_GPU_HOST_DEVICE AMREX_INLINE - void ConvertSphericalToCartesianYComponent (amrex::Real const F_r, + static void ConvertSphericalToCartesianYComponent (amrex::Real const F_r, amrex::Real const F_theta, amrex::Real const F_phi, amrex::Real const r, amrex::Real const theta, amrex::Real const phi, amrex::Real & F_y) { + amrex::ignore_unused(r); F_y = F_r * std::sin(theta) * std::sin(phi) + F_theta * std::cos(theta) * std::sin(phi) + F_phi * std::cos(phi); } AMREX_GPU_HOST_DEVICE AMREX_INLINE - void ConvertSphericalToCartesianZComponent (amrex::Real const F_r, + static void ConvertSphericalToCartesianZComponent (amrex::Real const F_r, amrex::Real const F_theta, amrex::Real const F_phi, amrex::Real const r, amrex::Real const theta, amrex::Real const phi, amrex::Real & F_z) { + amrex::ignore_unused(r, phi, F_phi); F_z = F_r * std::cos(theta) - F_theta * std::sin(theta); } AMREX_GPU_HOST_DEVICE AMREX_INLINE - void ConvertCartesianToSphericalRComponent (amrex::Real const F_x, + static void ConvertCartesianToSphericalRComponent (amrex::Real const F_x, amrex::Real const F_y, amrex::Real const F_z, amrex::Real const theta, amrex::Real const phi, amrex::Real & F_r) { @@ -283,7 +223,7 @@ namespace PulsarParm } AMREX_GPU_HOST_DEVICE AMREX_INLINE - void ConvertCartesianToSphericalThetaComponent (amrex::Real const F_x, + static void ConvertCartesianToSphericalThetaComponent (amrex::Real const F_x, amrex::Real const F_y, amrex::Real const F_z, amrex::Real const theta, amrex::Real const phi, amrex::Real & F_theta) { @@ -293,26 +233,89 @@ namespace PulsarParm } AMREX_GPU_HOST_DEVICE AMREX_INLINE - void ConvertCartesianToSphericalPhiComponent (amrex::Real const F_x, + static void ConvertCartesianToSphericalPhiComponent (amrex::Real const F_x, amrex::Real const F_y, amrex::Real const F_z, amrex::Real const theta, amrex::Real const phi, amrex::Real &F_phi) { + + amrex::ignore_unused(theta, F_z); F_phi = - F_x * std::sin(phi) + F_y * std::cos(phi); } AMREX_GPU_HOST_DEVICE AMREX_INLINE - void DampField(int i, int j, int k, amrex::GeometryData const& geom, amrex::Array4 const& Efield, amrex::GpuArray const mf_ixtype) + static void DampField(int i, int j, int k, + amrex::Array4 const& Efield, + amrex::GpuArray const mf_ixtype, + amrex::GpuArray const center_star_data, + amrex::GpuArray const domain_xlo, + amrex::GpuArray const domain_dx, + amrex::Real max_EBdamping_radius_data, + amrex::Real damping_scale_data, + amrex::Real Rstar_data) { - amrex::Real r = Spherical::r(i, j, k, geom, mf_ixtype); - if (r <= max_EBdamping_radius) { + amrex::Real x, y, z; + ComputeCellCoordinates(i, j, k, mf_ixtype, domain_xlo, domain_dx, x, y, z); + amrex::Real r, theta, phi; + ConvertCartesianToSphericalCoord (x, y, z, center_star_data, + r, theta, phi); + if (r <= max_EBdamping_radius_data) { // Damping function: Fd = tanh(damping_scale * (r / R_star - 1)) + 1 // for damping_scale >= 10 or so: // Fd(0) ~ 0 // Fd(R_star) ~ 1 - const amrex::Real Fd = std::tanh(damping_scale * (r / R_star - 1.0)) + 1.0; + const amrex::Real Fd = std::tanh(damping_scale_data * (r / Rstar_data - 1.0)) + 1.0; Efield(i, j, k) = Efield(i, j, k) * Fd; } } -} + + + static std::string m_pulsar_type; + static amrex::Real m_omega_star; + static amrex::Real m_R_star; + static amrex::Real m_B_star; + static amrex::Real m_dR_star; + static amrex::Real m_omega_ramp_time; + static amrex::Real m_field_damping_scale; + static int m_do_EB_external; + static int m_do_E_external_monopole; + static amrex::Array m_center_star; + static amrex::Real m_max_ndens; + static amrex::Real m_Ninj_fraction; + static int m_ModifyParticleWtAtInjection; + static amrex::Real m_rhoGJ_scale; + static int m_do_damp_EB_internal; + static amrex::Real m_max_EBcorotating_radius; + static amrex::Real m_max_EBdamping_radius; + static int m_turnoffdeposition; + static amrex::Real m_max_nodepos_radius; + static int m_turnoff_plasmaEB_gather; + static amrex::Real m_max_nogather_radius; + static int m_verbose; + static amrex::Real m_max_particle_absorption_radius; + static amrex::Real m_particle_inject_rmin; + static amrex::Real m_particle_inject_rmax; + static amrex::Real m_corotatingE_maxradius; + static amrex::Real m_enforceDipoleB_maxradius; + static int m_do_InitializeGrid_with_Pulsar_Bfield; + static int m_do_InitializeGrid_with_Pulsar_Efield; + static int m_enforceCorotatingE; + static int m_enforceDipoleB; + static int m_singleParticleTest; + static amrex::Real m_Bdamping_scale; + static int m_do_DampBDipoleInRing; + static amrex::Real m_injection_time; + static int m_continuous_injection; + static amrex::Real m_removeparticle_theta_min; + static amrex::Real m_removeparticle_theta_max; + static int m_use_theoreticalEB; + static amrex::Real m_theory_max_rstar; + static int m_LimitDipoleBInit; + static amrex::Real m_DipoleB_init_maxradius; + static int m_AddExternalMonopoleOnly; + static int m_AddMonopoleInsideRstarOnGrid; + static int m_EnforceTheoreticalEBInGrid; + +private: +}; #endif diff --git a/Source/Particles/PulsarParameters.cpp b/Source/Particles/PulsarParameters.cpp index a50844539ae..9c5a5c21ad7 100644 --- a/Source/Particles/PulsarParameters.cpp +++ b/Source/Particles/PulsarParameters.cpp @@ -6,554 +6,642 @@ #include #include "WarpX.H" -namespace PulsarParm + +std::string Pulsar::m_pulsar_type; +amrex::Real Pulsar::m_omega_star; +amrex::Real Pulsar::m_R_star; +amrex::Real Pulsar::m_B_star; +amrex::Real Pulsar::m_dR_star; +amrex::Real Pulsar::m_omega_ramp_time = 1.0; +amrex::Real Pulsar::m_field_damping_scale; +int Pulsar::m_do_EB_external = 0; +int Pulsar::m_do_E_external_monopole = 0; +amrex::Array Pulsar::m_center_star = {{0.}}; +amrex::Real Pulsar::m_max_ndens; +amrex::Real Pulsar::m_Ninj_fraction = 1.0; +int Pulsar::m_ModifyParticleWtAtInjection = 1.0; +amrex::Real Pulsar::m_rhoGJ_scale; +int Pulsar::m_do_damp_EB_internal = 0; +amrex::Real Pulsar::m_max_EBcorotating_radius; +amrex::Real Pulsar::m_max_EBdamping_radius; +int Pulsar::m_turnoffdeposition = 0; +amrex::Real Pulsar::m_max_nodepos_radius; +int Pulsar::m_turnoff_plasmaEB_gather = 0; +amrex::Real Pulsar::m_max_nogather_radius; +int Pulsar::m_verbose; +amrex::Real Pulsar::m_max_particle_absorption_radius; +amrex::Real Pulsar::m_particle_inject_rmin; +amrex::Real Pulsar::m_particle_inject_rmax; +amrex::Real Pulsar::m_corotatingE_maxradius; +amrex::Real Pulsar::m_enforceDipoleB_maxradius; +int Pulsar::m_do_InitializeGrid_with_Pulsar_Bfield = 1; +int Pulsar::m_do_InitializeGrid_with_Pulsar_Efield = 1; +int Pulsar::m_enforceCorotatingE = 1; +int Pulsar::m_enforceDipoleB = 1; +int Pulsar::m_singleParticleTest = 0; +int Pulsar::m_do_DampBDipoleInRing = 0; +amrex::Real Pulsar::m_Bdamping_scale; +amrex::Real Pulsar::m_injection_time = 0.; +int Pulsar::m_continuous_injection = 1; +amrex::Real Pulsar::m_removeparticle_theta_min = 180.; +amrex::Real Pulsar::m_removeparticle_theta_max = 0.; +int Pulsar::m_use_theoreticalEB = 0; +amrex::Real Pulsar::m_theory_max_rstar; +int Pulsar::m_LimitDipoleBInit = 0; +amrex::Real Pulsar::m_DipoleB_init_maxradius; +int Pulsar::m_AddExternalMonopoleOnly = 1; +int Pulsar::m_AddMonopoleInsideRstarOnGrid = 0; +int Pulsar::m_EnforceTheoreticalEBInGrid = 0; + + +Pulsar::Pulsar () { - std::string pulsar_type; + ReadParameters(); +} + +void +Pulsar::ReadParameters () { + amrex::ParmParse pp("pulsar"); + pp.query("pulsarType",m_pulsar_type); + + pp.get("omega_star",m_omega_star); - AMREX_GPU_DEVICE_MANAGED amrex::Real omega_star; - AMREX_GPU_DEVICE_MANAGED amrex::Real ramp_omega_time = -1.0; - AMREX_GPU_DEVICE_MANAGED amrex::Real B_star; - AMREX_GPU_DEVICE_MANAGED amrex::Real R_star; - AMREX_GPU_DEVICE_MANAGED amrex::Real dR_star; - AMREX_GPU_DEVICE_MANAGED amrex::Real damping_scale = 10.0; - AMREX_GPU_DEVICE_MANAGED int EB_external = 0; - AMREX_GPU_DEVICE_MANAGED int E_external_monopole = 0; - AMREX_GPU_DEVICE_MANAGED - amrex::GpuArray center_star; - AMREX_GPU_DEVICE_MANAGED int damp_EB_internal = 0; - AMREX_GPU_DEVICE_MANAGED int verbose = 0; - AMREX_GPU_DEVICE_MANAGED amrex::Real max_ndens; - AMREX_GPU_DEVICE_MANAGED amrex::Real Ninj_fraction; - AMREX_GPU_DEVICE_MANAGED int ModifyParticleWtAtInjection = 1; - AMREX_GPU_DEVICE_MANAGED amrex::Real rhoGJ_scale; - AMREX_GPU_DEVICE_MANAGED amrex::Real max_EBcorotating_radius; - AMREX_GPU_DEVICE_MANAGED amrex::Real max_EBdamping_radius; - AMREX_GPU_DEVICE_MANAGED int turnoffdeposition = 0; - AMREX_GPU_DEVICE_MANAGED amrex::Real max_nodepos_radius; - AMREX_GPU_DEVICE_MANAGED int turnoff_plasmaEB_gather = 0; - AMREX_GPU_DEVICE_MANAGED amrex::Real max_nogather_radius; - AMREX_GPU_DEVICE_MANAGED amrex::Real max_particle_absorption_radius; - AMREX_GPU_DEVICE_MANAGED amrex::Real particle_inject_rmin; - AMREX_GPU_DEVICE_MANAGED amrex::Real particle_inject_rmax; - AMREX_GPU_DEVICE_MANAGED amrex::Real corotatingE_maxradius; - AMREX_GPU_DEVICE_MANAGED amrex::Real enforceDipoleB_maxradius; - AMREX_GPU_DEVICE_MANAGED amrex::Real InitializeGrid_with_Pulsar_Bfield; - AMREX_GPU_DEVICE_MANAGED amrex::Real InitializeGrid_with_Pulsar_Efield; - AMREX_GPU_DEVICE_MANAGED int enforceCorotatingE; - AMREX_GPU_DEVICE_MANAGED int enforceDipoleB; - AMREX_GPU_DEVICE_MANAGED int singleParticleTest; - AMREX_GPU_DEVICE_MANAGED amrex::Real Bdamping_scale = 10; - AMREX_GPU_DEVICE_MANAGED int DampBDipoleInRing = 0; - AMREX_GPU_DEVICE_MANAGED amrex::Real injection_time = 0; - AMREX_GPU_DEVICE_MANAGED int continuous_injection = 1; - AMREX_GPU_DEVICE_MANAGED amrex::Real removeparticle_theta_min = 90.; - AMREX_GPU_DEVICE_MANAGED amrex::Real removeparticle_theta_max = 0.; - AMREX_GPU_DEVICE_MANAGED int use_theoreticalEB = 0; - AMREX_GPU_DEVICE_MANAGED amrex::Real theory_max_rstar = 0.; - AMREX_GPU_DEVICE_MANAGED int LimitDipoleBInit = 0; - AMREX_GPU_DEVICE_MANAGED amrex::Real DipoleB_init_maxradius; - AMREX_GPU_DEVICE_MANAGED int AddExternalMonopoleOnly = 0; - AMREX_GPU_DEVICE_MANAGED int AddMonopoleInsideRstarOnGrid = 0; - AMREX_GPU_DEVICE_MANAGED int EnforceTheoreticalEBInGrid = 0; + amrex::Vector center_star_v(AMREX_SPACEDIM); + pp.getarr("center_star",center_star_v); + std::copy(center_star_v.begin(),center_star_v.end(),m_center_star.begin()); - void ReadParameters() { - amrex::ParmParse pp("pulsar"); - pp.query("pulsarType",pulsar_type); - pp.get("omega_star",omega_star); - amrex::Vector center_star_v(AMREX_SPACEDIM); - pp.getarr("center_star",center_star_v); - std::copy(center_star_v.begin(),center_star_v.end(),center_star.begin()); - pp.get("R_star",R_star); - pp.get("B_star",B_star); - pp.get("dR",dR_star); - pp.query("verbose",verbose); - pp.query("EB_external",EB_external); - pp.query("E_external_monopole",E_external_monopole); - pp.query("damp_EB_internal",damp_EB_internal); - pp.query("damping_scale",damping_scale); - pp.query("ramp_omega_time",ramp_omega_time); - amrex::Print() << " Pulsar center: " << center_star[0] << " " << center_star[1] << " " << center_star[2] << "\n"; - amrex::Print() << " Pulsar omega: " << omega_star << "\n"; - amrex::Print() << " Pulsar B_star : " << B_star << "\n"; - pp.get("max_ndens", max_ndens); - pp.get("Ninj_fraction",Ninj_fraction); + pp.get("R_star",m_R_star); + pp.get("B_star",m_B_star); + pp.get("dR",m_dR_star); + pp.query("verbose",m_verbose); + pp.query("EB_external",m_do_EB_external); + pp.query("E_external_monopole",m_do_E_external_monopole); + pp.query("damp_EB_internal",m_do_damp_EB_internal); + pp.query("damping_scale",m_field_damping_scale); + pp.query("ramp_omega_time",m_omega_ramp_time); + amrex::Print() << " Pulsar center: " << m_center_star[0] << " " << m_center_star[1] << " " << m_center_star[2] << "\n"; + amrex::Print() << " Pulsar omega: " << m_omega_star << "\n"; + amrex::Print() << " Pulsar B_star : " << m_B_star << "\n"; + pp.get("max_ndens", m_max_ndens); + pp.get("Ninj_fraction",m_Ninj_fraction); - particle_inject_rmin = R_star - dR_star; - particle_inject_rmax = R_star; - pp.query("particle_inj_rmin", particle_inject_rmin); - pp.query("particle_inj_rmax", particle_inject_rmax); - amrex::Print() << " min radius of particle injection : " << particle_inject_rmin << "\n"; - amrex::Print() << " max radius of particle injection : " << particle_inject_rmax << "\n"; + m_particle_inject_rmin = m_R_star - m_dR_star; + m_particle_inject_rmax = m_R_star; + pp.query("particle_inj_rmin", m_particle_inject_rmin); + pp.query("particle_inj_rmax", m_particle_inject_rmax); + amrex::Print() << " min radius of particle injection : " << m_particle_inject_rmin << "\n"; + amrex::Print() << " max radius of particle injection : " << m_particle_inject_rmax << "\n"; - pp.query("ModifyParticleWeight", ModifyParticleWtAtInjection); - // The maximum radius within which particles are absorbed/deleted every timestep. - max_particle_absorption_radius = R_star; - pp.query("max_particle_absorption_radius", max_particle_absorption_radius); - pp.get("rhoGJ_scale",rhoGJ_scale); - amrex::Print() << " pulsar max ndens " << max_ndens << "\n"; - amrex::Print() << " pulsar ninj fraction " << Ninj_fraction << "\n"; - amrex::Print() << " pulsar modify particle wt " << ModifyParticleWtAtInjection << "\n"; - amrex::Print() << " pulsar rhoGJ scaling " << rhoGJ_scale << "\n"; - amrex::Print() << " EB_external : " << EB_external << "\n"; - if (EB_external == 1) { - // Max corotating radius defines the region where the EB field shifts from - // corotating (v X B) to quadrapole. default is R_star. - max_EBcorotating_radius = R_star; - pp.query("EB_corotating_maxradius", max_EBcorotating_radius); - amrex::Print() << " EB coratating maxradius : " << max_EBcorotating_radius << "\n"; - } - if (damp_EB_internal == 1) { - // Radius of the region within which the EB fields are damped - max_EBdamping_radius = R_star; - pp.query("damp_EB_radius", max_EBdamping_radius); - amrex::Print() << " max EB damping radius : " << max_EBdamping_radius << "\n"; - } - // query to see if particle j,rho deposition should be turned off - // within some region interior to the star - // default is 0 - pp.query("turnoffdeposition", turnoffdeposition); - amrex::Print() << " is deposition off ? " << turnoffdeposition << "\n"; - if (turnoffdeposition == 1) { - max_nodepos_radius = R_star; - pp.query("max_nodepos_radius", max_nodepos_radius); - amrex::Print() << " deposition turned off within radius : " << max_nodepos_radius << "\n"; - } - pp.query("turnoff_plasmaEB_gather", turnoff_plasmaEB_gather); - amrex::Print() << " is plasma EB gather off ? " << turnoff_plasmaEB_gather << "\n"; - if (turnoff_plasmaEB_gather == 1) { - max_nogather_radius = R_star; - pp.query("max_nogather_radius", max_nogather_radius); - amrex::Print() << " gather off within radius : " << max_nogather_radius << "\n"; - } - corotatingE_maxradius = R_star; - enforceDipoleB_maxradius = R_star; - pp.query("corotatingE_maxradius", corotatingE_maxradius); - pp.query("enforceDipoleB_maxradius", enforceDipoleB_maxradius); - InitializeGrid_with_Pulsar_Bfield = 1; - InitializeGrid_with_Pulsar_Efield = 1; - pp.query("init_dipoleBfield", InitializeGrid_with_Pulsar_Bfield); - pp.query("init_corotatingEfield", InitializeGrid_with_Pulsar_Efield); - enforceCorotatingE = 1; - enforceDipoleB = 1; - pp.query("enforceCorotatingE", enforceCorotatingE); - pp.query("enforceDipoleB", enforceDipoleB); - singleParticleTest = 0; - pp.query("singleParticleTest", singleParticleTest); - pp.query("DampBDipoleInRing", DampBDipoleInRing); - if (DampBDipoleInRing == 1) pp.query("Bdamping_scale", Bdamping_scale); - pp.query("injection_time", injection_time); - pp.query("continuous_injection", continuous_injection); - pp.query("removeparticle_theta_min",removeparticle_theta_min); - pp.query("removeparticle_theta_max",removeparticle_theta_max); - pp.query("use_theoreticalEB",use_theoreticalEB); - amrex::Print() << "use theory EB " << use_theoreticalEB << "\n"; - pp.query("theory_max_rstar",theory_max_rstar); - amrex::Print() << " theory max rstar : " << theory_max_rstar << "\n"; - pp.query("LimitDipoleBInit", LimitDipoleBInit); - if (LimitDipoleBInit == 1) { - pp.query("DipoleB_init_maxradius", DipoleB_init_maxradius); - } - pp.query("AddExternalMonopoleOnly", AddExternalMonopoleOnly); - pp.query("AddMonopoleInsideRstarOnGrid", AddMonopoleInsideRstarOnGrid); - pp.query("EnforceTheoreticalEBInGrid", EnforceTheoreticalEBInGrid); + pp.query("ModifyParticleWeight", m_ModifyParticleWtAtInjection); + // The maximum radius within which particles are absorbed/deleted every timestep. + m_max_particle_absorption_radius = m_R_star; + pp.query("max_particle_absorption_radius", m_max_particle_absorption_radius); + pp.get("rhoGJ_scale",m_rhoGJ_scale); + amrex::Print() << " pulsar max ndens " << m_max_ndens << "\n"; + amrex::Print() << " pulsar ninj fraction " << m_Ninj_fraction << "\n"; + amrex::Print() << " pulsar modify particle wt " << m_ModifyParticleWtAtInjection << "\n"; + amrex::Print() << " pulsar rhoGJ scaling " << m_rhoGJ_scale << "\n"; + amrex::Print() << " EB_external : " << m_do_EB_external << "\n"; + if (m_do_EB_external == 1) { + // Max corotating radius defines the region where the EB field shifts from + // corotating (v X B) to quadrapole. default is R_star. + m_max_EBcorotating_radius = m_R_star; + pp.query("EB_corotating_maxradius", m_max_EBcorotating_radius); + amrex::Print() << " EB coratating maxradius : " << m_max_EBcorotating_radius << "\n"; } + if (m_do_damp_EB_internal == 1) { + // Radius of the region within which the EB fields are damped + m_max_EBdamping_radius = m_R_star; + pp.query("damp_EB_radius", m_max_EBdamping_radius); + amrex::Print() << " max EB damping radius : " << m_max_EBdamping_radius << "\n"; + } + // query to see if particle j,rho deposition should be turned off + // within some region interior to the star + // default is 0 + pp.query("turnoffdeposition", m_turnoffdeposition); + amrex::Print() << " is deposition off ? " << m_turnoffdeposition << "\n"; + if (m_turnoffdeposition == 1) { + m_max_nodepos_radius = m_R_star; + pp.query("max_nodepos_radius", m_max_nodepos_radius); + amrex::Print() << " deposition turned off within radius : " << m_max_nodepos_radius << "\n"; + } + pp.query("turnoff_plasmaEB_gather", m_turnoff_plasmaEB_gather); + amrex::Print() << " is plasma EB gather off ? " << m_turnoff_plasmaEB_gather << "\n"; + if (m_turnoff_plasmaEB_gather == 1) { + m_max_nogather_radius = m_R_star; + pp.query("max_nogather_radius", m_max_nogather_radius); + amrex::Print() << " gather off within radius : " << m_max_nogather_radius << "\n"; + } + m_corotatingE_maxradius = m_R_star; + m_enforceDipoleB_maxradius = m_R_star; + pp.query("corotatingE_maxradius", m_corotatingE_maxradius); + pp.query("enforceDipoleB_maxradius", m_enforceDipoleB_maxradius); + pp.query("init_dipoleBfield", m_do_InitializeGrid_with_Pulsar_Bfield); + pp.query("init_corotatingEfield", m_do_InitializeGrid_with_Pulsar_Efield); + pp.query("enforceCorotatingE", m_enforceCorotatingE); + pp.query("enforceDipoleB", m_enforceDipoleB); + pp.query("singleParticleTest", m_singleParticleTest); + pp.query("DampBDipoleInRing", m_do_DampBDipoleInRing); + if (m_do_DampBDipoleInRing == 1) pp.query("Bdamping_scale", m_Bdamping_scale); + pp.query("injection_time", m_injection_time); + pp.query("continuous_injection", m_continuous_injection); + pp.query("removeparticle_theta_min",m_removeparticle_theta_min); + pp.query("removeparticle_theta_max",m_removeparticle_theta_max); + pp.query("use_theoreticalEB",m_use_theoreticalEB); + amrex::Print() << "use theory EB " << m_use_theoreticalEB << "\n"; + pp.query("theory_max_rstar",m_theory_max_rstar); + amrex::Print() << " theory max rstar : " << m_theory_max_rstar << "\n"; + pp.query("LimitDipoleBInit", m_LimitDipoleBInit); + if (m_LimitDipoleBInit == 1) { + pp.query("DipoleB_init_maxradius", m_DipoleB_init_maxradius); + } + pp.query("AddExternalMonopoleOnly", m_AddExternalMonopoleOnly); + pp.query("AddMonopoleInsideRstarOnGrid", m_AddMonopoleInsideRstarOnGrid); + pp.query("EnforceTheoreticalEBInGrid", m_EnforceTheoreticalEBInGrid); +} - /** To initialize the grid with dipole magnetic field everywhere and corotating vacuum - * electric field inside the pulsar radius. - */ - void InitializeExternalPulsarFieldsOnGrid ( amrex::MultiFab *mfx, amrex::MultiFab *mfy, - amrex::MultiFab *mfz, const int lev, const bool init_Bfield) - { - auto & warpx = WarpX::GetInstance(); - const auto dx = warpx.Geom(lev).CellSizeArray(); - const auto problo = warpx.Geom(lev).ProbLoArray(); - const auto probhi = warpx.Geom(lev).ProbHiArray(); - const RealBox& real_box = warpx.Geom(lev).ProbDomain(); - amrex::Real cur_time = warpx.gett_new(lev); - amrex::IntVect x_nodal_flag = mfx->ixType().toIntVect(); - amrex::IntVect y_nodal_flag = mfy->ixType().toIntVect(); - amrex::IntVect z_nodal_flag = mfz->ixType().toIntVect(); - GpuArray x_IndexType; - GpuArray y_IndexType; - GpuArray z_IndexType; - for (int idim = 0; idim < 3; ++idim) { - x_IndexType[idim] = x_nodal_flag[idim]; - y_IndexType[idim] = y_nodal_flag[idim]; - z_IndexType[idim] = z_nodal_flag[idim]; - } +/** To initialize the grid with dipole magnetic field everywhere and corotating vacuum + * electric field inside the pulsar radius. + */ +void +Pulsar::InitializeExternalPulsarFieldsOnGrid ( amrex::MultiFab *mfx, amrex::MultiFab *mfy, + amrex::MultiFab *mfz, const int lev, const bool init_Bfield) +{ + auto & warpx = WarpX::GetInstance(); + const auto dx = warpx.Geom(lev).CellSizeArray(); + const auto problo = warpx.Geom(lev).ProbLoArray(); + amrex::Real cur_time = warpx.gett_new(lev); + amrex::IntVect x_nodal_flag = mfx->ixType().toIntVect(); + amrex::IntVect y_nodal_flag = mfy->ixType().toIntVect(); + amrex::IntVect z_nodal_flag = mfz->ixType().toIntVect(); + GpuArray x_IndexType; + GpuArray y_IndexType; + GpuArray z_IndexType; + GpuArray center_star_arr; + for (int idim = 0; idim < 3; ++idim) { + x_IndexType[idim] = x_nodal_flag[idim]; + y_IndexType[idim] = y_nodal_flag[idim]; + z_IndexType[idim] = z_nodal_flag[idim]; + center_star_arr[idim] = m_center_star[idim]; + } + amrex::Real omega_star_data = m_omega_star; + amrex::Real ramp_omega_time_data = m_omega_ramp_time; + amrex::Real Bstar_data = m_B_star; + amrex::Real Rstar_data = m_R_star; + amrex::Real dRstar_data = m_dR_star; + int LimitDipoleBInit_data = m_LimitDipoleBInit; + amrex::Real DipoleB_init_maxradius_data = m_DipoleB_init_maxradius; + amrex::Real corotatingE_maxradius_data = m_corotatingE_maxradius; #ifdef AMREX_USE_OMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif - for (MFIter mfi(*mfx, TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - const amrex::Box& tbx = mfi.tilebox(x_nodal_flag, mfx->nGrowVect() ); - const amrex::Box& tby = mfi.tilebox(y_nodal_flag, mfx->nGrowVect() ); - const amrex::Box& tbz = mfi.tilebox(z_nodal_flag, mfx->nGrowVect() ); - - amrex::Array4 const& mfx_arr = mfx->array(mfi); - amrex::Array4 const& mfy_arr = mfy->array(mfi); - amrex::Array4 const& mfz_arr = mfz->array(mfi); + for (MFIter mfi(*mfx, TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const amrex::Box& tbx = mfi.tilebox(x_nodal_flag, mfx->nGrowVect() ); + const amrex::Box& tby = mfi.tilebox(y_nodal_flag, mfx->nGrowVect() ); + const amrex::Box& tbz = mfi.tilebox(z_nodal_flag, mfx->nGrowVect() ); + + amrex::Array4 const& mfx_arr = mfx->array(mfi); + amrex::Array4 const& mfy_arr = mfy->array(mfi); + amrex::Array4 const& mfz_arr = mfz->array(mfi); - amrex::ParallelFor (tbx, tby, tbz, - [=] AMREX_GPU_DEVICE (int i, int j, int k) { - amrex::Real x, y, z; - amrex::Real r, theta, phi; - amrex::Real Fr, Ftheta, Fphi; - Fr = 0.; Ftheta = 0.; Fphi = 0.; - // compute cell coordinates - ComputeCellCoordinates(i, j, k, x_IndexType, problo, dx, x, y, z); - // convert cartesian to spherical coordinates - ConvertCartesianToSphericalCoord(x, y, z, problo, probhi, r, theta, phi); - // Initialize with Bfield in spherical coordinates - if (init_Bfield == 1) { - if (LimitDipoleBInit == 1) { - if (r < DipoleB_init_maxradius) { - ExternalBFieldSpherical( r, theta, phi, cur_time, - Fr, Ftheta, Fphi); - } - } else { + amrex::ParallelFor (tbx, tby, tbz, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + amrex::Real x, y, z; + amrex::Real r, theta, phi; + amrex::Real Fr, Ftheta, Fphi; + Fr = 0.; Ftheta = 0.; Fphi = 0.; + // compute cell coordinates + ComputeCellCoordinates(i, j, k, x_IndexType, problo, dx, x, y, z); + // convert cartesian to spherical coordinates + ConvertCartesianToSphericalCoord(x, y, z, center_star_arr, + r, theta, phi); + // Initialize with Bfield in spherical coordinates + if (init_Bfield == 1) { + if (LimitDipoleBInit_data == 1) { + if (r < DipoleB_init_maxradius_data) { ExternalBFieldSpherical( r, theta, phi, cur_time, + Bstar_data, Rstar_data, dRstar_data, Fr, Ftheta, Fphi); } + } else { + ExternalBFieldSpherical( r, theta, phi, cur_time, + Bstar_data, Rstar_data, dRstar_data, + Fr, Ftheta, Fphi); } - // Initialize corotating EField in r < corotating - if (init_Bfield == 0) { - if (r <= corotatingE_maxradius) { - CorotatingEfieldSpherical(r, theta, phi, cur_time, - Fr, Ftheta, Fphi); - } + } + // Initialize corotating EField in r < corotating + if (init_Bfield == 0) { + if (r <= corotatingE_maxradius_data) { + CorotatingEfieldSpherical(r, theta, phi, cur_time, + omega_star_data, + ramp_omega_time_data, + Bstar_data, Rstar_data, dRstar_data, + Fr, Ftheta, Fphi); } - // Convert to x component - ConvertSphericalToCartesianXComponent( Fr, Ftheta, Fphi, - r, theta, phi, mfx_arr(i,j,k)); - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) { - amrex::Real x, y, z; - amrex::Real r, theta, phi; - amrex::Real Fr, Ftheta, Fphi; - Fr = 0.; Ftheta = 0.; Fphi = 0.; - // compute cell coordinates - ComputeCellCoordinates(i, j, k, y_IndexType, problo, dx, x, y, z); - // convert cartesian to spherical coordinates - ConvertCartesianToSphericalCoord(x, y, z, problo, probhi, r, theta, phi); - // Initialize with Bfield in spherical coordinates - if (init_Bfield == 1) { - if (LimitDipoleBInit == 1) { - if (r < DipoleB_init_maxradius) { - ExternalBFieldSpherical( r, theta, phi, cur_time, - Fr, Ftheta, Fphi); - } - } else { + } + // Convert to x component + ConvertSphericalToCartesianXComponent( Fr, Ftheta, Fphi, + r, theta, phi, mfx_arr(i,j,k)); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + amrex::Real x, y, z; + amrex::Real r, theta, phi; + amrex::Real Fr, Ftheta, Fphi; + Fr = 0.; Ftheta = 0.; Fphi = 0.; + // compute cell coordinates + ComputeCellCoordinates(i, j, k, y_IndexType, problo, dx, x, y, z); + // convert cartesian to spherical coordinates + ConvertCartesianToSphericalCoord(x, y, z, center_star_arr, + r, theta, phi); + // Initialize with Bfield in spherical coordinates + if (init_Bfield == 1) { + if (LimitDipoleBInit_data == 1) { + if (r < DipoleB_init_maxradius_data) { ExternalBFieldSpherical( r, theta, phi, cur_time, + Bstar_data, Rstar_data, dRstar_data, Fr, Ftheta, Fphi); } + } else { + ExternalBFieldSpherical( r, theta, phi, cur_time, + Bstar_data, Rstar_data, dRstar_data, + Fr, Ftheta, Fphi); } - // Initialize corotating Efield in r < corotating - if (init_Bfield == 0) { - if (r <= corotatingE_maxradius) { - CorotatingEfieldSpherical(r, theta, phi, cur_time, - Fr, Ftheta, Fphi); - } + } + // Initialize corotating Efield in r < corotating + if (init_Bfield == 0) { + if (r <= corotatingE_maxradius_data) { + CorotatingEfieldSpherical(r, theta, phi, cur_time, + omega_star_data, + ramp_omega_time_data, + Bstar_data, Rstar_data, dRstar_data, + Fr, Ftheta, Fphi); } - // convert to y component - ConvertSphericalToCartesianYComponent( Fr, Ftheta, Fphi, - r, theta, phi, mfy_arr(i,j,k)); - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) { - amrex::Real x, y, z; - amrex::Real r, theta, phi; - amrex::Real Fr, Ftheta, Fphi; - Fr = 0.; Ftheta = 0.; Fphi = 0.; - // compute cell coordinates - ComputeCellCoordinates(i, j, k, z_IndexType, problo, dx, x, y, z); - // convert cartesian to spherical coordinates - ConvertCartesianToSphericalCoord(x, y, z, problo, probhi, r, theta, phi); - // Initialize with Bfield in spherical coordinates - if (init_Bfield == 1) { - if (LimitDipoleBInit == 1) { - if (r < DipoleB_init_maxradius) { - ExternalBFieldSpherical( r, theta, phi, cur_time, - Fr, Ftheta, Fphi); - } - } else { + } + // convert to y component + ConvertSphericalToCartesianYComponent( Fr, Ftheta, Fphi, + r, theta, phi, mfy_arr(i,j,k)); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + amrex::Real x, y, z; + amrex::Real r, theta, phi; + amrex::Real Fr, Ftheta, Fphi; + Fr = 0.; Ftheta = 0.; Fphi = 0.; + // compute cell coordinates + ComputeCellCoordinates(i, j, k, z_IndexType, problo, dx, x, y, z); + // convert cartesian to spherical coordinates + ConvertCartesianToSphericalCoord(x, y, z, center_star_arr, + r, theta, phi); + // Initialize with Bfield in spherical coordinates + if (init_Bfield == 1) { + if (LimitDipoleBInit_data == 1) { + if (r < DipoleB_init_maxradius_data) { ExternalBFieldSpherical( r, theta, phi, cur_time, + Bstar_data, Rstar_data, dRstar_data, Fr, Ftheta, Fphi); } + } else { + ExternalBFieldSpherical( r, theta, phi, cur_time, + Bstar_data, Rstar_data, dRstar_data, + Fr, Ftheta, Fphi); } - // Initialize corotating Efield in r < corotating - if (init_Bfield == 0) { - if (r <= corotatingE_maxradius) { - CorotatingEfieldSpherical(r, theta, phi, cur_time, - Fr, Ftheta, Fphi); - } + } + // Initialize corotating Efield in r < corotating + if (init_Bfield == 0) { + if (r <= corotatingE_maxradius_data) { + CorotatingEfieldSpherical(r, theta, phi, cur_time, + omega_star_data, + ramp_omega_time_data, + Bstar_data, Rstar_data, dRstar_data, + Fr, Ftheta, Fphi); } - ConvertSphericalToCartesianZComponent( Fr, Ftheta, Fphi, - r, theta, phi, mfz_arr(i,j,k)); } - ); - } // mfiter loop - } // InitializeExternalField + ConvertSphericalToCartesianZComponent( Fr, Ftheta, Fphi, + r, theta, phi, mfz_arr(i,j,k)); + } + ); + } // mfiter loop +} // InitializeExternalField - void ApplyCorotatingEfield_BC ( std::array< std::unique_ptr, 3> &Efield, - const int lev, const amrex::Real a_dt) - { - amrex::Print() << " applying corotating Efield BC\n"; - auto & warpx = WarpX::GetInstance(); - const auto dx = warpx.Geom(lev).CellSizeArray(); - const auto problo = warpx.Geom(lev).ProbLoArray(); - const auto probhi = warpx.Geom(lev).ProbHiArray(); - const RealBox& real_box = warpx.Geom(lev).ProbDomain(); - amrex::Real cur_time = warpx.gett_new(lev) + a_dt; - amrex::IntVect x_nodal_flag = Efield[0]->ixType().toIntVect(); - amrex::IntVect y_nodal_flag = Efield[1]->ixType().toIntVect(); - amrex::IntVect z_nodal_flag = Efield[2]->ixType().toIntVect(); - GpuArray x_IndexType; - GpuArray y_IndexType; - GpuArray z_IndexType; - for (int idim = 0; idim < 3; ++idim) { - x_IndexType[idim] = x_nodal_flag[idim]; - y_IndexType[idim] = y_nodal_flag[idim]; - z_IndexType[idim] = z_nodal_flag[idim]; - } +void +Pulsar::ApplyCorotatingEfield_BC ( std::array< std::unique_ptr, 3> &Efield, + const int lev, const amrex::Real a_dt) +{ + amrex::Print() << " applying corotating Efield BC\n"; + auto & warpx = WarpX::GetInstance(); + const auto dx = warpx.Geom(lev).CellSizeArray(); + const auto problo = warpx.Geom(lev).ProbLoArray(); + amrex::Real cur_time = warpx.gett_new(lev) + a_dt; + amrex::IntVect x_nodal_flag = Efield[0]->ixType().toIntVect(); + amrex::IntVect y_nodal_flag = Efield[1]->ixType().toIntVect(); + amrex::IntVect z_nodal_flag = Efield[2]->ixType().toIntVect(); + GpuArray x_IndexType; + GpuArray y_IndexType; + GpuArray z_IndexType; + GpuArray center_star_arr; + for (int idim = 0; idim < 3; ++idim) { + x_IndexType[idim] = x_nodal_flag[idim]; + y_IndexType[idim] = y_nodal_flag[idim]; + z_IndexType[idim] = z_nodal_flag[idim]; + center_star_arr[idim] = m_center_star[idim]; + } + amrex::Real omega_star_data = m_omega_star; + amrex::Real ramp_omega_time_data = m_omega_ramp_time; + amrex::Real Bstar_data = m_B_star; + amrex::Real Rstar_data = m_R_star; + amrex::Real dRstar_data = m_dR_star; + amrex::Real corotatingE_maxradius_data = m_corotatingE_maxradius; + int E_external_monopole_data = m_do_E_external_monopole; + int EnforceTheoreticalEBInGrid_data = m_EnforceTheoreticalEBInGrid; #ifdef AMREX_USE_OMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif - for (MFIter mfi(*Efield[0], TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - const amrex::Box& tex = mfi.tilebox(x_nodal_flag); - const amrex::Box& tey = mfi.tilebox(y_nodal_flag); - const amrex::Box& tez = mfi.tilebox(z_nodal_flag); - amrex::Array4 const& Ex_arr = Efield[0]->array(mfi); - amrex::Array4 const& Ey_arr = Efield[1]->array(mfi); - amrex::Array4 const& Ez_arr = Efield[2]->array(mfi); - // loop over cells and set Efield for r < corotating - amrex::ParallelFor(tex, tey, tez, - [=] AMREX_GPU_DEVICE (int i, int j, int k) { - amrex::Real x, y, z; - amrex::Real r, theta, phi; - amrex::Real Fr, Ftheta, Fphi; - Fr = 0.; Ftheta = 0.; Fphi = 0.; - // compute cell coordinates - ComputeCellCoordinates(i, j, k, x_IndexType, problo, dx, x, y, z); - // convert cartesian to spherical coordinates - ConvertCartesianToSphericalCoord(x, y, z, problo, probhi, r, theta, phi); - if (PulsarParm::EnforceTheoreticalEBInGrid == 0) { - if (r <= corotatingE_maxradius) { - CorotatingEfieldSpherical(r, theta, phi, cur_time, - Fr, Ftheta, Fphi); - ConvertSphericalToCartesianXComponent( Fr, Ftheta, Fphi, - r, theta, phi, Ex_arr(i,j,k)); - } - } else { - ExternalEFieldSpherical(r, theta, phi, cur_time, Fr, Ftheta, Fphi); + for (MFIter mfi(*Efield[0], TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const amrex::Box& tex = mfi.tilebox(x_nodal_flag); + const amrex::Box& tey = mfi.tilebox(y_nodal_flag); + const amrex::Box& tez = mfi.tilebox(z_nodal_flag); + amrex::Array4 const& Ex_arr = Efield[0]->array(mfi); + amrex::Array4 const& Ey_arr = Efield[1]->array(mfi); + amrex::Array4 const& Ez_arr = Efield[2]->array(mfi); + // loop over cells and set Efield for r < corotating + amrex::ParallelFor(tex, tey, tez, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + amrex::Real x, y, z; + amrex::Real r, theta, phi; + amrex::Real Fr, Ftheta, Fphi; + Fr = 0.; Ftheta = 0.; Fphi = 0.; + // compute cell coordinates + ComputeCellCoordinates(i, j, k, x_IndexType, problo, dx, x, y, z); + // convert cartesian to spherical coordinates + ConvertCartesianToSphericalCoord(x, y, z, center_star_arr, + r, theta, phi); + if (EnforceTheoreticalEBInGrid_data == 0) { + if (r <= corotatingE_maxradius_data) { + CorotatingEfieldSpherical(r, theta, phi, cur_time, + omega_star_data, + ramp_omega_time_data, + Bstar_data, Rstar_data, dRstar_data, + Fr, Ftheta, Fphi); ConvertSphericalToCartesianXComponent( Fr, Ftheta, Fphi, r, theta, phi, Ex_arr(i,j,k)); } - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) { - amrex::Real x, y, z; - amrex::Real r, theta, phi; - amrex::Real Fr, Ftheta, Fphi; - Fr = 0.; Ftheta = 0.; Fphi = 0.; - // compute cell coordinates - ComputeCellCoordinates(i, j, k, y_IndexType, problo, dx, x, y, z); - // convert cartesian to spherical coordinates - ConvertCartesianToSphericalCoord(x, y, z, problo, probhi, r, theta, phi); - if (PulsarParm::EnforceTheoreticalEBInGrid == 0) { - if (r <= corotatingE_maxradius) { - CorotatingEfieldSpherical(r, theta, phi, cur_time, - Fr, Ftheta, Fphi); - ConvertSphericalToCartesianYComponent( Fr, Ftheta, Fphi, - r, theta, phi, Ey_arr(i,j,k)); - } - } else { - ExternalEFieldSpherical(r, theta, phi, cur_time, Fr, Ftheta, Fphi); + } else { + ExternalEFieldSpherical(r, theta, phi, cur_time, + omega_star_data, + ramp_omega_time_data, + Bstar_data, Rstar_data, + corotatingE_maxradius_data, + E_external_monopole_data, + Fr, Ftheta, Fphi); + ConvertSphericalToCartesianXComponent( Fr, Ftheta, Fphi, + r, theta, phi, Ex_arr(i,j,k)); + } + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + amrex::Real x, y, z; + amrex::Real r, theta, phi; + amrex::Real Fr, Ftheta, Fphi; + Fr = 0.; Ftheta = 0.; Fphi = 0.; + // compute cell coordinates + ComputeCellCoordinates(i, j, k, y_IndexType, problo, dx, x, y, z); + // convert cartesian to spherical coordinates + ConvertCartesianToSphericalCoord(x, y, z, center_star_arr, + r, theta, phi); + if (EnforceTheoreticalEBInGrid_data == 0) { + if (r <= corotatingE_maxradius_data) { + CorotatingEfieldSpherical(r, theta, phi, cur_time, + omega_star_data, + ramp_omega_time_data, + Bstar_data, Rstar_data, dRstar_data, + Fr, Ftheta, Fphi); ConvertSphericalToCartesianYComponent( Fr, Ftheta, Fphi, r, theta, phi, Ey_arr(i,j,k)); } - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) { - amrex::Real x, y, z; - amrex::Real r, theta, phi; - amrex::Real Fr, Ftheta, Fphi; - Fr = 0.; Ftheta = 0.; Fphi = 0.; - // compute cell coordinates - ComputeCellCoordinates(i, j, k, z_IndexType, problo, dx, x, y, z); - // convert cartesian to spherical coordinates - ConvertCartesianToSphericalCoord(x, y, z, problo, probhi, r, theta, phi); - if (PulsarParm::EnforceTheoreticalEBInGrid == 0) { - if (r <= corotatingE_maxradius) { - CorotatingEfieldSpherical(r, theta, phi, cur_time, - Fr, Ftheta, Fphi); - ConvertSphericalToCartesianZComponent( Fr, Ftheta, Fphi, - r, theta, phi, Ez_arr(i,j,k)); - } - } else { - ExternalEFieldSpherical(r, theta, phi, cur_time, Fr, Ftheta, Fphi); + } else { + ExternalEFieldSpherical(r, theta, phi, cur_time, + omega_star_data, + ramp_omega_time_data, + Bstar_data, Rstar_data, + corotatingE_maxradius_data, + E_external_monopole_data, + Fr, Ftheta, Fphi); + ConvertSphericalToCartesianYComponent( Fr, Ftheta, Fphi, + r, theta, phi, Ey_arr(i,j,k)); + } + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + amrex::Real x, y, z; + amrex::Real r, theta, phi; + amrex::Real Fr, Ftheta, Fphi; + Fr = 0.; Ftheta = 0.; Fphi = 0.; + // compute cell coordinates + ComputeCellCoordinates(i, j, k, z_IndexType, problo, dx, x, y, z); + // convert cartesian to spherical coordinates + ConvertCartesianToSphericalCoord(x, y, z, center_star_arr, + r, theta, phi); + if (EnforceTheoreticalEBInGrid_data == 0) { + if (r <= corotatingE_maxradius_data) { + CorotatingEfieldSpherical(r, theta, phi, cur_time, + omega_star_data, + ramp_omega_time_data, + Bstar_data, Rstar_data, dRstar_data, + Fr, Ftheta, Fphi); ConvertSphericalToCartesianZComponent( Fr, Ftheta, Fphi, r, theta, phi, Ez_arr(i,j,k)); } + } else { + ExternalEFieldSpherical(r, theta, phi, cur_time, + omega_star_data, + ramp_omega_time_data, + Bstar_data, Rstar_data, + corotatingE_maxradius_data, + E_external_monopole_data, + Fr, Ftheta, Fphi); + ConvertSphericalToCartesianZComponent( Fr, Ftheta, Fphi, + r, theta, phi, Ez_arr(i,j,k)); } - ); - } - - + } + ); } - void ApplyDipoleBfield_BC ( std::array< std::unique_ptr, 3> &Bfield, - const int lev, const amrex::Real a_dt) - { - amrex::Print() << " applying corotating Bfield BC\n"; - auto & warpx = WarpX::GetInstance(); - const auto dx = warpx.Geom(lev).CellSizeArray(); - const auto problo = warpx.Geom(lev).ProbLoArray(); - const auto probhi = warpx.Geom(lev).ProbHiArray(); - const RealBox& real_box = warpx.Geom(lev).ProbDomain(); - amrex::Real cur_time = warpx.gett_new(lev) + a_dt; - amrex::IntVect x_nodal_flag = Bfield[0]->ixType().toIntVect(); - amrex::IntVect y_nodal_flag = Bfield[1]->ixType().toIntVect(); - amrex::IntVect z_nodal_flag = Bfield[2]->ixType().toIntVect(); - GpuArray x_IndexType; - GpuArray y_IndexType; - GpuArray z_IndexType; - for (int idim = 0; idim < 3; ++idim) { - x_IndexType[idim] = x_nodal_flag[idim]; - y_IndexType[idim] = y_nodal_flag[idim]; - z_IndexType[idim] = z_nodal_flag[idim]; - } + +} + +void +Pulsar::ApplyDipoleBfield_BC ( std::array< std::unique_ptr, 3> &Bfield, + const int lev, const amrex::Real a_dt) +{ + amrex::Print() << " applying corotating Bfield BC\n"; + auto & warpx = WarpX::GetInstance(); + const auto dx = warpx.Geom(lev).CellSizeArray(); + const auto problo = warpx.Geom(lev).ProbLoArray(); + amrex::Real cur_time = warpx.gett_new(lev) + a_dt; + amrex::IntVect x_nodal_flag = Bfield[0]->ixType().toIntVect(); + amrex::IntVect y_nodal_flag = Bfield[1]->ixType().toIntVect(); + amrex::IntVect z_nodal_flag = Bfield[2]->ixType().toIntVect(); + GpuArray x_IndexType; + GpuArray y_IndexType; + GpuArray z_IndexType; + GpuArray center_star_arr; + for (int idim = 0; idim < 3; ++idim) { + x_IndexType[idim] = x_nodal_flag[idim]; + y_IndexType[idim] = y_nodal_flag[idim]; + z_IndexType[idim] = z_nodal_flag[idim]; + center_star_arr[idim] = m_center_star[idim]; + } + amrex::Real Bstar_data = m_B_star; + amrex::Real Rstar_data = m_R_star; + amrex::Real dRstar_data = m_dR_star; + int EnforceTheoreticalEBInGrid_data = m_EnforceTheoreticalEBInGrid; + amrex::Real corotatingE_maxradius_data = m_corotatingE_maxradius; + amrex::Real enforceDipoleB_maxradius_data = m_enforceDipoleB_maxradius; + int DampBDipoleInRing_data = m_do_DampBDipoleInRing; + amrex::Real Bdamping_scale_data = m_Bdamping_scale; #ifdef AMREX_USE_OMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif - for (MFIter mfi(*Bfield[0], TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - const amrex::Box& tbx = mfi.tilebox(x_nodal_flag); - const amrex::Box& tby = mfi.tilebox(y_nodal_flag); - const amrex::Box& tbz = mfi.tilebox(z_nodal_flag); - amrex::Array4 const& Bx_arr = Bfield[0]->array(mfi); - amrex::Array4 const& By_arr = Bfield[1]->array(mfi); - amrex::Array4 const& Bz_arr = Bfield[2]->array(mfi); - // loop over cells and set Efield for r < dipoleB_max_radius - amrex::ParallelFor(tbx, tby, tbz, - [=] AMREX_GPU_DEVICE (int i, int j, int k) { - amrex::Real x, y, z; - amrex::Real r, theta, phi; - amrex::Real Fr, Ftheta, Fphi; - Fr = 0.; Ftheta = 0.; Fphi = 0.; - // compute cell coordinates - ComputeCellCoordinates(i, j, k, x_IndexType, problo, dx, x, y, z); - // convert cartesian to spherical coordinates - ConvertCartesianToSphericalCoord(x, y, z, problo, probhi, r, theta, phi); - if (PulsarParm::EnforceTheoreticalEBInGrid == 0) { - if (r <= enforceDipoleB_maxradius) { - ExternalBFieldSpherical(r, theta, phi, cur_time, - Fr, Ftheta, Fphi); - ConvertSphericalToCartesianXComponent( Fr, Ftheta, Fphi, - r, theta, phi, Bx_arr(i,j,k)); - } - else if ( r > enforceDipoleB_maxradius && r <= corotatingE_maxradius) { - if (DampBDipoleInRing == 1) { - // from inner ring to outer ring - ExternalBFieldSpherical(r, theta, phi, cur_time, - Fr, Ftheta, Fphi); - amrex::Real Bx_dipole; - ConvertSphericalToCartesianXComponent( Fr, Ftheta, Fphi, - r, theta, phi, Bx_dipole); - // Damping Function : Fd = tanh(dampingscale * (1-r/Rinner)) - // where Rinner = enforceDipoleB_maxradius - // is the range where Bdipole is imposed - // Fd(Rinner) ~ 1 - // Fd(R_domainboundary) ~ 0 - amrex::Real Fd = 1._rt - + std::tanh( Bdamping_scale - * (1._rt - r/enforceDipoleB_maxradius) - ); - Bx_arr(i,j,k) += Fd * Bx_dipole; - } - } - } else { + for (MFIter mfi(*Bfield[0], TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const amrex::Box& tbx = mfi.tilebox(x_nodal_flag); + const amrex::Box& tby = mfi.tilebox(y_nodal_flag); + const amrex::Box& tbz = mfi.tilebox(z_nodal_flag); + amrex::Array4 const& Bx_arr = Bfield[0]->array(mfi); + amrex::Array4 const& By_arr = Bfield[1]->array(mfi); + amrex::Array4 const& Bz_arr = Bfield[2]->array(mfi); + // loop over cells and set Efield for r < dipoleB_max_radius + amrex::ParallelFor(tbx, tby, tbz, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + amrex::Real x, y, z; + amrex::Real r, theta, phi; + amrex::Real Fr, Ftheta, Fphi; + Fr = 0.; Ftheta = 0.; Fphi = 0.; + // compute cell coordinates + ComputeCellCoordinates(i, j, k, x_IndexType, problo, dx, x, y, z); + // convert cartesian to spherical coordinates + ConvertCartesianToSphericalCoord(x, y, z, center_star_arr, + r, theta, phi); + if (EnforceTheoreticalEBInGrid_data == 0) { + if (r <= enforceDipoleB_maxradius_data) { ExternalBFieldSpherical(r, theta, phi, cur_time, - Fr, Ftheta, Fphi); + Bstar_data, Rstar_data, dRstar_data, + Fr, Ftheta, Fphi); ConvertSphericalToCartesianXComponent( Fr, Ftheta, Fphi, r, theta, phi, Bx_arr(i,j,k)); } - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) { - amrex::Real x, y, z; - amrex::Real r, theta, phi; - amrex::Real Fr, Ftheta, Fphi; - Fr = 0.; Ftheta = 0.; Fphi = 0.; - // compute cell coordinates - ComputeCellCoordinates(i, j, k, y_IndexType, problo, dx, x, y, z); - // convert cartesian to spherical coordinates - ConvertCartesianToSphericalCoord(x, y, z, problo, probhi, r, theta, phi); - if (PulsarParm::EnforceTheoreticalEBInGrid == 0) { - if (r <= enforceDipoleB_maxradius) { + else if ( r > enforceDipoleB_maxradius_data && r <= corotatingE_maxradius_data) { + if (DampBDipoleInRing_data == 1) { + // from inner ring to outer ring ExternalBFieldSpherical(r, theta, phi, cur_time, + Bstar_data, Rstar_data, dRstar_data, Fr, Ftheta, Fphi); - ConvertSphericalToCartesianYComponent( Fr, Ftheta, Fphi, - r, theta, phi, By_arr(i,j,k)); - } else if ( r > enforceDipoleB_maxradius && r <= corotatingE_maxradius) { - if (DampBDipoleInRing == 1) { - // from inner ring to outer ring - ExternalBFieldSpherical(r, theta, phi, cur_time, - Fr, Ftheta, Fphi); - amrex::Real By_dipole; - ConvertSphericalToCartesianYComponent( Fr, Ftheta, Fphi, - r, theta, phi, By_dipole); - amrex::Real Fd = 1._rt - + std::tanh( Bdamping_scale - * (1._rt - r/enforceDipoleB_maxradius) - ); - By_arr(i,j,k) += Fd * By_dipole; - } + amrex::Real Bx_dipole; + ConvertSphericalToCartesianXComponent( Fr, Ftheta, Fphi, + r, theta, phi, Bx_dipole); + // Damping Function : Fd = tanh(dampingscale * (1-r/Rinner)) + // where Rinner = enforceDipoleB_maxradius + // is the range where Bdipole is imposed + // Fd(Rinner) ~ 1 + // Fd(R_domainboundary) ~ 0 + amrex::Real Fd = 1._rt + + std::tanh( Bdamping_scale_data + * (1._rt - r/enforceDipoleB_maxradius_data) + ); + Bx_arr(i,j,k) += Fd * Bx_dipole; } - } else { + } + } else { + ExternalBFieldSpherical(r, theta, phi, cur_time, + Bstar_data, Rstar_data, dRstar_data, + Fr, Ftheta, Fphi); + ConvertSphericalToCartesianXComponent( Fr, Ftheta, Fphi, + r, theta, phi, Bx_arr(i,j,k)); + } + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + amrex::Real x, y, z; + amrex::Real r, theta, phi; + amrex::Real Fr, Ftheta, Fphi; + Fr = 0.; Ftheta = 0.; Fphi = 0.; + // compute cell coordinates + ComputeCellCoordinates(i, j, k, y_IndexType, problo, dx, x, y, z); + // convert cartesian to spherical coordinates + ConvertCartesianToSphericalCoord(x, y, z, center_star_arr, + r, theta, phi); + if (EnforceTheoreticalEBInGrid_data == 0) { + if (r <= enforceDipoleB_maxradius_data) { ExternalBFieldSpherical(r, theta, phi, cur_time, + Bstar_data, Rstar_data, dRstar_data, Fr, Ftheta, Fphi); ConvertSphericalToCartesianYComponent( Fr, Ftheta, Fphi, r, theta, phi, By_arr(i,j,k)); - } - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) { - amrex::Real x, y, z; - amrex::Real r, theta, phi; - amrex::Real Fr, Ftheta, Fphi; - Fr = 0.; Ftheta = 0.; Fphi = 0.; - // compute cell coordinates - ComputeCellCoordinates(i, j, k, z_IndexType, problo, dx, x, y, z); - // convert cartesian to spherical coordinates - ConvertCartesianToSphericalCoord(x, y, z, problo, probhi, r, theta, phi); - if (PulsarParm::EnforceTheoreticalEBInGrid == 0) { - if (r <= enforceDipoleB_maxradius) { + } else if ( r > enforceDipoleB_maxradius_data && r <= corotatingE_maxradius_data) { + if (DampBDipoleInRing_data == 1) { + // from inner ring to outer ring ExternalBFieldSpherical(r, theta, phi, cur_time, + Bstar_data, Rstar_data, dRstar_data, Fr, Ftheta, Fphi); - ConvertSphericalToCartesianZComponent( Fr, Ftheta, Fphi, - r, theta, phi, Bz_arr(i,j,k)); - } else if ( r > enforceDipoleB_maxradius && r <= corotatingE_maxradius) { - if (DampBDipoleInRing == 1) { - // from inner ring to outer ring - ExternalBFieldSpherical(r, theta, phi, cur_time, - Fr, Ftheta, Fphi); - amrex::Real Bz_dipole; - ConvertSphericalToCartesianZComponent( Fr, Ftheta, Fphi, - r, theta, phi, Bz_dipole); - amrex::Real Fd = 1._rt - + std::tanh( Bdamping_scale - * (1._rt - r/enforceDipoleB_maxradius) - ); - Bz_arr(i,j,k) += Fd * Bz_dipole; - } + amrex::Real By_dipole; + ConvertSphericalToCartesianYComponent( Fr, Ftheta, Fphi, + r, theta, phi, By_dipole); + amrex::Real Fd = 1._rt + + std::tanh( Bdamping_scale_data + * (1._rt - r/enforceDipoleB_maxradius_data) + ); + By_arr(i,j,k) += Fd * By_dipole; } - } else { + } + } else { + ExternalBFieldSpherical(r, theta, phi, cur_time, + Bstar_data, Rstar_data, dRstar_data, + Fr, Ftheta, Fphi); + ConvertSphericalToCartesianYComponent( Fr, Ftheta, Fphi, + r, theta, phi, By_arr(i,j,k)); + } + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + amrex::Real x, y, z; + amrex::Real r, theta, phi; + amrex::Real Fr, Ftheta, Fphi; + Fr = 0.; Ftheta = 0.; Fphi = 0.; + // compute cell coordinates + ComputeCellCoordinates(i, j, k, z_IndexType, problo, dx, x, y, z); + // convert cartesian to spherical coordinates + ConvertCartesianToSphericalCoord(x, y, z, center_star_arr, + r, theta, phi); + if (EnforceTheoreticalEBInGrid_data == 0) { + if (r <= enforceDipoleB_maxradius_data) { ExternalBFieldSpherical(r, theta, phi, cur_time, + Bstar_data, Rstar_data, dRstar_data, Fr, Ftheta, Fphi); ConvertSphericalToCartesianZComponent( Fr, Ftheta, Fphi, r, theta, phi, Bz_arr(i,j,k)); + } else if ( r > enforceDipoleB_maxradius_data && r <= corotatingE_maxradius_data) { + if (DampBDipoleInRing_data == 1) { + // from inner ring to outer ring + ExternalBFieldSpherical(r, theta, phi, cur_time, + Bstar_data, Rstar_data, dRstar_data, + Fr, Ftheta, Fphi); + amrex::Real Bz_dipole; + ConvertSphericalToCartesianZComponent( Fr, Ftheta, Fphi, + r, theta, phi, Bz_dipole); + amrex::Real Fd = 1._rt + + std::tanh( Bdamping_scale_data + * (1._rt - r/enforceDipoleB_maxradius_data) + ); + Bz_arr(i,j,k) += Fd * Bz_dipole; + } } + } else { + ExternalBFieldSpherical(r, theta, phi, cur_time, + Bstar_data, Rstar_data, dRstar_data, + Fr, Ftheta, Fphi); + ConvertSphericalToCartesianZComponent( Fr, Ftheta, Fphi, + r, theta, phi, Bz_arr(i,j,k)); } - ); - } - + } + ); } + } diff --git a/Source/Particles/Pusher/PushSelector.H b/Source/Particles/Pusher/PushSelector.H index 207d51cf90b..c34de0cffd1 100644 --- a/Source/Particles/Pusher/PushSelector.H +++ b/Source/Particles/Pusher/PushSelector.H @@ -116,13 +116,21 @@ void doParticlePush(const GetParticlePosition& GetPosition, if (ion_lev) { qp *= ion_lev; } UpdateMomentumBoris( ux, uy, uz, Ex, Ey, Ez, Bx, - By, Bz, qp, m, dt, PulsarDiag, i, compute_from_theory); + By, Bz, qp, m, dt +#ifdef PULSAR + , PulsarDiag, i, compute_from_theory +#endif + ); +#ifdef PULSAR if (compute_from_theory == 0) { +#endif amrex::ParticleReal x, y, z; GetPosition(i, x, y, z); UpdatePosition(x, y, z, ux, uy, uz, dt ); SetPosition(i, x, y, z); +#ifdef PULSAR } +#endif } else if (pusher_algo == ParticlePusherAlgo::Vay) { amrex::Real qp = q; if (ion_lev){ qp *= ion_lev; } diff --git a/Source/Particles/Pusher/UpdateMomentumBoris.H b/Source/Particles/Pusher/UpdateMomentumBoris.H index 64d9281b8a1..aacd95b5984 100644 --- a/Source/Particles/Pusher/UpdateMomentumBoris.H +++ b/Source/Particles/Pusher/UpdateMomentumBoris.H @@ -17,7 +17,11 @@ void UpdateMomentumBoris( amrex::ParticleReal& ux, amrex::ParticleReal& uy, amrex::ParticleReal& uz, const amrex::ParticleReal Ex, const amrex::ParticleReal Ey, const amrex::ParticleReal Ez, const amrex::ParticleReal Bx, const amrex::ParticleReal By, const amrex::ParticleReal Bz, - const amrex::Real q, const amrex::Real m, const amrex::Real dt, amrex::Real * PulsarDiag = nullptr, const int i=0, const int compute_from_theory = 0 ) + const amrex::Real q, const amrex::Real m, const amrex::Real dt +#ifdef PULSAR + , amrex::Real * PulsarDiag = nullptr, const int i=0, const int compute_from_theory = 0 +#endif + ) { using namespace amrex::literals; @@ -69,10 +73,10 @@ void UpdateMomentumBoris( ux += uy_p*sz - uz_p*sy; uy += uz_p*sx - ux_p*sz; uz += ux_p*sy - uy_p*sx; +#ifdef PULSAR amrex::Real tmp_ux = uy_p*sz - uz_p*sy; amrex::Real tmp_uy = uz_p*sx - ux_p*sz; amrex::Real tmp_uz = ux_p*sy - uy_p*sx; -#ifdef PULSAR // amrex::AllPrintToFile("PulsarParticle") << " Momentum after rotation : " << ux << " " << " uy : " << uy << " uz " << uz << " rot ux : " << tmp_ux << " rot uy " << tmp_uy << " rot uz " << tmp_uz<< "\n"; if (i == 0) { if (compute_from_theory == 0) { diff --git a/Source/Utils/WarpXAlgorithmSelection.H b/Source/Utils/WarpXAlgorithmSelection.H index ba752e012b7..17c7c0d674c 100644 --- a/Source/Utils/WarpXAlgorithmSelection.H +++ b/Source/Utils/WarpXAlgorithmSelection.H @@ -129,6 +129,13 @@ struct ReductionType { }; }; +struct IntegrationType { + enum { + Volume = 0, + Surface = 1 + }; +}; + int GetAlgorithmInteger( amrex::ParmParse& pp, const char* pp_search_key ); diff --git a/Source/Utils/WarpXAlgorithmSelection.cpp b/Source/Utils/WarpXAlgorithmSelection.cpp index 670f195d018..7a4e960680c 100644 --- a/Source/Utils/WarpXAlgorithmSelection.cpp +++ b/Source/Utils/WarpXAlgorithmSelection.cpp @@ -103,7 +103,14 @@ const std::map ParticleBCType_algo_to_enum = const std::map ReductionType_algo_to_int = { {"maximum", ReductionType::Maximum}, {"minimum", ReductionType::Minimum}, - {"integral", ReductionType::Sum} + {"integral", ReductionType::Sum}, + {"surfaceintegral", ReductionType::Sum} +}; + +const std::map IntegrationType_algo_to_int = { + {"volume", IntegrationType::Volume}, + {"surface", IntegrationType::Surface}, + {"default", IntegrationType::Volume} }; int @@ -139,6 +146,8 @@ GetAlgorithmInteger( amrex::ParmParse& pp, const char* pp_search_key ){ algo_to_int = MacroscopicSolver_algo_to_int; } else if (0 == std::strcmp(pp_search_key, "reduction_type")) { algo_to_int = ReductionType_algo_to_int; + } else if (0 == std::strcmp(pp_search_key, "integration_type")) { + algo_to_int = IntegrationType_algo_to_int; } else { std::string pp_search_string = pp_search_key; amrex::Abort("Unknown algorithm type: " + pp_search_string); diff --git a/Source/WarpX.H b/Source/WarpX.H index cd4bce33548..be8b896df9f 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -35,6 +35,9 @@ #include "Particles/WarpXParticleContainer_fwd.H" #include "Utils/IntervalsParser.H" #include "Utils/WarpXAlgorithmSelection.H" +#ifdef PULSAR +# include "Particles/PulsarParameters.H" +#endif #include #include @@ -1050,7 +1053,6 @@ private: // Macroscopic properties std::unique_ptr m_macroscopic_properties; - // Load balancing /** Load balancing intervals that reads the "load_balance_intervals" string int the input file * for getting steps at which load balancing is performed */ diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index fe853a209b6..2840872ff79 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -327,7 +327,6 @@ WarpX::WarpX () m_macroscopic_properties = std::make_unique(); } - // Set default values for particle and cell weights for costs update; // Default values listed here for the case AMREX_USE_GPU are determined // from single-GPU tests on Summit. @@ -1156,8 +1155,9 @@ WarpX::ReadParameters () } } + #ifdef PULSAR - PulsarParm::ReadParameters(); + Pulsar::ReadParameters(); #endif }