Skip to content

Commit

Permalink
damp b field and include staggering. fix compilation errors
Browse files Browse the repository at this point in the history
  • Loading branch information
@RevathiJambunathan
RevathiJambunathan committed Jan 6, 2021
1 parent 31bf403 commit 21a552a
Show file tree
Hide file tree
Showing 5 changed files with 67 additions and 15 deletions.
4 changes: 2 additions & 2 deletions Examples/Physics_applications/pulsar/inputs.corotating.3d.PM
View file
Original file line number Diff line number Diff line change
Expand Up @@ -146,9 +146,9 @@ pulsar.B_star = 8.0323e-6 # magnetic field of NS (T)
pulsar.dR = 1075
pulsar.verbose = 0 # [0/1]: turn on verbosity for debugging print statements
pulsar.EB_external = 1 # [0/1]: to apply external E and B
pulsar.E_external_monopole = 1 # [0/1]
pulsar.E_external_monopole = 0 # [0/1]
pulsar.max_ndens = 5.54e6 # max ndens == ndens used in initializing density
pulsar.Ninj_fraction = 0.1 # fraction of sigma injected
pulsar.rhoGJ_scale = 1e0 # scaling down of rho_GJ
pulsar.damp_E_internal = 0 # Damp internal electric field
pulsar.damp_EB_internal = 1 # Damp internal electric field
pulsar.damping_scale = 10.0 # Damping scale factor for internal electric field
64 changes: 59 additions & 5 deletions Source/Evolve/WarpXEvolve.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -140,13 +140,43 @@ WarpX::Evolve (int numsteps)
}

#ifdef PULSAR
if (PulsarParm::damp_E_internal) {
if (PulsarParm::damp_EB_internal) {
MultiFab *Ex, *Ey, *Ez;
MultiFab *Bx, *By, *Bz;
for (int lev = 0; lev <= finest_level; ++lev) {
Ex = Efield_fp[lev][0].get();
Ey = Efield_fp[lev][1].get();
Ez = Efield_fp[lev][2].get();

Bx = Bfield_fp[lev][0].get();
By = Bfield_fp[lev][1].get();
Bz = Bfield_fp[lev][2].get();
Gpu::ManagedVector<int> Ex_stag, Ey_stag, Ez_stag, Bx_stag, By_stag, Bz_stag;
Ex_stag.resize(3);
Ey_stag.resize(3);
Ez_stag.resize(3);
Bx_stag.resize(3);
By_stag.resize(3);
Bz_stag.resize(3);
amrex::IntVect ex_type = Ex->ixType().toIntVect();
amrex::IntVect ey_type = Ey->ixType().toIntVect();
amrex::IntVect ez_type = Ez->ixType().toIntVect();
amrex::IntVect bx_type = Bx->ixType().toIntVect();
amrex::IntVect by_type = By->ixType().toIntVect();
amrex::IntVect bz_type = Bz->ixType().toIntVect();
for (int idim = 0; idim < AMREX_SPACEDIM-1; ++idim) {
Ex_stag[idim] = ex_type[idim];
Ey_stag[idim] = ey_type[idim];
Ez_stag[idim] = ez_type[idim];
Bx_stag[idim] = bx_type[idim];
By_stag[idim] = by_type[idim];
Bz_stag[idim] = bz_type[idim];
}
int const* const AMREX_RESTRICT Ex_stag_ptr = Ex_stag.data();
int const* const AMREX_RESTRICT Ey_stag_ptr = Ey_stag.data();
int const* const AMREX_RESTRICT Ez_stag_ptr = Ez_stag.data();
int const* const AMREX_RESTRICT Bx_stag_ptr = Bx_stag.data();
int const* const AMREX_RESTRICT By_stag_ptr = By_stag.data();
int const* const AMREX_RESTRICT Bz_stag_ptr = Bz_stag.data();
auto geom = Geom(lev).data();
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
Expand All @@ -164,15 +194,39 @@ WarpX::Evolve (int numsteps)
amrex::ParallelFor(tex, tey, tez,
[=] AMREX_GPU_DEVICE (int i, int j, int k)
{
PulsarParm::DampEField(i, j, k, geom, Exfab);
PulsarParm::DampEField(i, j, k, geom, Exfab, Ex_stag_ptr);
},
[=] AMREX_GPU_DEVICE (int i, int j, int k)
{
PulsarParm::DampEField(i, j, k, geom, Eyfab, Ey_stag_ptr);
},
[=] AMREX_GPU_DEVICE (int i, int j, int k)
{
PulsarParm::DampEField(i, j, k, geom, Ezfab, Ez_stag_ptr);
});
}
for ( MFIter mfi(*Bx, TilingIfNotGPU()); mfi.isValid(); ++mfi )
{
const Box& tex = mfi.tilebox( Bx->ixType().toIntVect() );
const Box& tey = mfi.tilebox( By->ixType().toIntVect() );
const Box& tez = mfi.tilebox( Bz->ixType().toIntVect() );

auto const& Bxfab = Bx->array(mfi);
auto const& Byfab = By->array(mfi);
auto const& Bzfab = Bz->array(mfi);

amrex::ParallelFor(tex, tey, tez,
[=] AMREX_GPU_DEVICE (int i, int j, int k)
{
PulsarParm::DampEField(i, j, k, geom, Bxfab, Bx_stag_ptr);
},
[=] AMREX_GPU_DEVICE (int i, int j, int k)
{
PulsarParm::DampEField(i, j, k, geom, Eyfab);
PulsarParm::DampEField(i, j, k, geom, Byfab, By_stag_ptr);
},
[=] AMREX_GPU_DEVICE (int i, int j, int k)
{
PulsarParm::DampEField(i, j, k, geom, Ezfab);
PulsarParm::DampEField(i, j, k, geom, Bzfab, Bz_stag_ptr);
});
}
}
Expand Down
2 changes: 0 additions & 2 deletions Source/Particles/PhysicalParticleContainer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1110,7 +1110,6 @@ PhysicalParticleContainer::Evolve (int lev,
tmp_particle_data[t_lev][index][i].resize(np);
}
}
amrex::Print() << " we ppc evolve \n";
#ifdef _OPENMP
#pragma omp parallel
#endif
Expand All @@ -1123,7 +1122,6 @@ PhysicalParticleContainer::Evolve (int lev,

FArrayBox filtered_Ex, filtered_Ey, filtered_Ez;
FArrayBox filtered_Bx, filtered_By, filtered_Bz;
amrex::Print() << " par iter loop \n" ;
for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti)
{
if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers)
Expand Down
8 changes: 4 additions & 4 deletions Source/Particles/PulsarParameters.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ namespace PulsarParm
extern AMREX_GPU_DEVICE_MANAGED amrex::Real max_ndens;
extern AMREX_GPU_DEVICE_MANAGED amrex::Real Ninj_fraction;
extern AMREX_GPU_DEVICE_MANAGED amrex::Real rhoGJ_scale;
extern AMREX_GPU_DEVICE_MANAGED int damp_E_internal;
extern AMREX_GPU_DEVICE_MANAGED int damp_EB_internal;

extern AMREX_GPU_DEVICE_MANAGED int verbose;

Expand Down Expand Up @@ -219,7 +219,7 @@ namespace PulsarParm
namespace Spherical
{
AMREX_GPU_HOST_DEVICE AMREX_INLINE
amrex::Real r(int i, int j, int k, amrex::GeometryData const& geom, const amrex::IntVect mf_ixType)
amrex::Real r(int i, int j, int k, amrex::GeometryData const& geom, int const* const AMREX_RESTRICT mf_ixType)
{
const auto domain_box = geom.Domain();
const auto domain_ilo = amrex::lbound(domain_box);
Expand All @@ -246,9 +246,9 @@ namespace PulsarParm


AMREX_GPU_HOST_DEVICE AMREX_INLINE
void DampEField(int i, int j, int k, amrex::GeometryData const& geom, amrex::Array4<amrex::Real> const& Efield, const amrex::IntVect mf_ixtype)
void DampEField(int i, int j, int k, amrex::GeometryData const& geom, amrex::Array4<amrex::Real> const& Efield, int const* const AMREX_RESTRICT mf_ixtype)
{
const amrex::Real r = Spherical::r(i, j, k, geom, mf_ixType);
const amrex::Real r = Spherical::r(i, j, k, geom, mf_ixtype);

if (r < R_star) {
// Damping function: Fd = tanh(damping_scale * (r / R_star - 1)) + 1
Expand Down
4 changes: 2 additions & 2 deletions Source/Particles/PulsarParameters.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@ namespace PulsarParm
AMREX_GPU_DEVICE_MANAGED int E_external_monopole = 0;
AMREX_GPU_DEVICE_MANAGED
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> center_star;
AMREX_GPU_DEVICE_MANAGED int damp_E_internal = 0;
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;
Expand All @@ -38,7 +38,7 @@ namespace PulsarParm
pp.query("verbose",verbose);
pp.query("EB_external",EB_external);
pp.query("E_external_monopole",E_external_monopole);
pp.query("damp_E_internal",damp_E_internal);
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";
Expand Down

0 comments on commit 21a552a

Please sign in to comment.