Add oblique angle

CMakeLists.txt

@@ -241,6 +241,9 @@ endif()
 target_compile_definitions(WarpX PRIVATE
     WARPX_PARSER_DEPTH=${WarpX_PARSER_DEPTH})
 
+target_compile_definitions(WarpX PRIVATE PULSAR)
+
+
 
 # Warnings ####################################################################
 #

Examples/Physics_applications/pulsar/analysis_pulsar_vis.py

@@ -0,0 +1,54 @@
+import yt
+import os
+import glob
+import argparse
+
+parser = argparse.ArgumentParser()
+parser.add_argument("-f", "--plotfile", type=str, help="Path to single input plotfile to visualize.")
+parser.add_argument("-d", "--dir", type=str, help="Path to input plotfile directory to visualize.")
+
+args = parser.parse_args()
+
+def visualize(plt, annotate_particles=True):
+    ds = yt.load(plt)
+    sl = yt.SlicePlot(ds, 0, 'Ez', aspect=1) # Create a sliceplot object
+    if annotate_particles:
+        sl.annotate_particles(width=(2.e3, 'm'),ptype="plasma_p",col='b',p_size=5.0)
+        sl.annotate_particles(width=(2.e3, 'm'),ptype="plasma_e",col='r',p_size=5.0)
+    sl.annotate_streamlines("Ey", "Ez", plot_args={"color": "black"})
+    sl.annotate_grids() # Show grids
+    sl.annotate_sphere([90000.0, 90000.0, 90000.0], radius=12000.0,
+                              circle_args={'color':'white', 'linewidth':2})
+    out_name = os.path.join(os.path.dirname(plt), "pulsar_{}.png".format(os.path.basename(plt)))
+    sl.save(out_name)
+
+if __name__ == "__main__":
+    yt.funcs.mylog.setLevel(50)
+    if args.dir:
+        failed = []
+        plotfiles = glob.glob(os.path.join(args.dir, "plt" + "[0-9]"*5))
+        for pf in plotfiles:
+            try:
+                visualize(pf)
+            except:
+                # plotfile 0 may not have particles, so turn off annotations if we first failed
+                try:
+                    visualize(pf, annotate_particles=False)
+                except:
+                    failed.append(pf)
+                    pass
+
+        if len(failed) > 0:
+            print("Visualization failed for the following plotfiles:")
+            for fp in failed:
+                print(fp)
+        else:
+            print("No plots failed, creating a gif ...")
+            input_files = os.path.join(args.dir, "*.png")
+            output_gif  = os.path.join(args.dir, "pulsar.gif")
+            os.system("convert -delay 20 -loop 0 {} {}".format(input_files, output_gif))
+
+    elif args.plotfile:
+        visualize(args.plotfile)
+    else:
+        print("Supply either -f or -d options for visualization. Use the -h option to see help.")

Examples/Physics_applications/pulsar/inputs.corotating.3d.PM

@@ -0,0 +1,120 @@
+# Description:
+#
+# This inputs file sets up a NS with these properties:
+# - E = -(omega*R)[cross]B inside the NS
+# - external E and B outside the star
+#
+# See the "Pulsar Setup" section at the end for the options
+#
+# This initializes the electrons and positrons with a corotating momentum function.
+# Based on the pulsar_type = "active" or "dead", particles are injected continuously or
+#                             until rho_GJ is reached
+
+#################################
+####### GENERAL PARAMETERS ######
+#################################
+max_step = 50000
+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_hi     = 180000 180000 180000
+
+#################################
+############ NUMERICS ###########
+#################################
+#algo.maxwell_fdtd_solver = yee
+warpx.verbose = 1
+warpx.do_dive_cleaning = 0
+warpx.use_filter = 1
+warpx.cfl = .2
+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  = 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.to = 2.e-4
+interpolation.nox = 3
+interpolation.noy = 3
+interpolation.noz = 3
+
+#################################
+############ PLASMA #############
+#################################
+particles.nspecies = 2
+particles.species_names = plasma_e plasma_p
+
+plasma_e.charge = -q_e
+plasma_e.mass = m_e
+plasma_e.injection_style = "NUniformPerCell"
+plasma_e.profile = parse_density_function
+plasma_e.density_function(x,y,z) = "( ((( (z-zc)*(z-zc) + (y-yc)*(y-yc) + (x-xc)*(x-xc) )^(0.5))<=(r_star)) * ((( (z-zc)*(z-zc) + (y-yc)*(y-yc) + (x-xc)*(x-xc) )^(0.5))>=(r_star-dR)) )*dens"
+plasma_e.num_particles_per_cell_each_dim = 4 4 4
+plasma_e.momentum_distribution_type = parse_momentum_function
+
+## Inject stationary pairs
+plasma_e.momentum_function_ux(x,y,z) = "0.0"
+plasma_e.momentum_function_uy(x,y,z) = "0.0"
+## uz is always 0 for the injection methods above
+plasma_e.momentum_function_uz(x,y,z) = "0.0"
+
+plasma_e.do_continuous_injection = 0
+plasma_e.density_min = 1E0
+
+plasma_p.charge = q_e
+plasma_p.mass = m_e
+plasma_p.injection_style = "NUniformPerCell"
+plasma_p.profile = parse_density_function
+plasma_p.density_function(x,y,z) = "( ((( (z-zc)*(z-zc) + (y-yc)*(y-yc) + (x-xc)*(x-xc) )^(0.5))<=(r_star)) * ((( (z-zc)*(z-zc) + (y-yc)*(y-yc) + (x-xc)*(x-xc) )^(0.5))>=(r_star-dR)) )*dens"
+plasma_p.num_particles_per_cell_each_dim = 4 4 4
+plasma_p.momentum_distribution_type = parse_momentum_function
+
+## Inject stationary pairs
+plasma_p.momentum_function_ux(x,y,z) = "0.0"
+plasma_p.momentum_function_uy(x,y,z) = "0.0"
+## uz is always 0 for the injection methods above
+plasma_p.momentum_function_uz(x,y,z) = "0.0"
+
+plasma_p.do_continuous_injection = 0
+plasma_p.density_min = 1E0
+
+diagnostics.diags_names = plt chk
+plt.diag_type = Full
+plt.period = 20
+plt.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz part_per_cell rho divE
+
+chk.format = checkpoint
+chk.diag_type = Full
+chk.period = 1000
+
+#################################
+######### PULSAR SETUP ##########
+#################################
+pulsar.pulsarType = "dead"             # [dead/active]: sets particle injection type
+pulsar.omega_star = 6245.676               # angular velocity of NS (rad/s)
+pulsar.ramp_omega_time = 0.0             # time over which to ramp up NS angular velocity (s)
+                                           # if ramp_omega_time < 0, omega = omega_star for any t
+                                           # consistency requires ramp_omega_time = my_constants.to
+pulsar.center_star = 90000 90000 90000
+pulsar.R_star = 12032                  # radius of NS (m)
+pulsar.B_star = 8.0323e-6                  # magnetic field of NS (T)
+pulsar.dR = 1400                       # thickness on the surface of the pulsar
+pulsar.Chi = 45
+                                       # consistency requires dR = my_constants.dR
+pulsar.verbose = 0                     # [0/1]: turn on verbosity for debugging print statements
+pulsar.EB_external = 1                 # [0/1]: to apply external E and B
+pulsar.E_external_monopole = 1                  # [0/1]
+pulsar.max_ndens = 5.54e6              # max ndens == ndens used in initializing density
+pulsar.Ninj_fraction = 0.2            # fraction of sigma injected
+pulsar.rhoGJ_scale = 1e0              # scaling down of rho_GJ
+pulsar.damp_EB_internal = 1             # Damp internal electric field
+pulsar.damping_scale = 1000.0            # Damping scale factor for internal electric field

Examples/Physics_applications/pulsar/pulsar_vis.py

@@ -0,0 +1,54 @@
+import yt
+import os
+import glob
+import argparse
+
+parser = argparse.ArgumentParser()
+parser.add_argument("-f", "--plotfile", type=str, help="Path to single input plotfile to visualize.")
+parser.add_argument("-d", "--dir", type=str, help="Path to input plotfile directory to visualize.")
+
+args = parser.parse_args()
+
+def visualize(plt, annotate_particles=True):
+    ds = yt.load(plt)
+    sl = yt.SlicePlot(ds, 0, 'Ez', aspect=1) # Create a sliceplot object
+    if annotate_particles:
+        sl.annotate_particles(width=(2.e3, 'm'),ptype="plasma_p",col='b',p_size=5.0)
+        sl.annotate_particles(width=(2.e3, 'm'),ptype="plasma_e",col='r',p_size=5.0)
+    sl.annotate_streamlines("Ey", "Ez", plot_args={"color": "black"})
+    sl.annotate_grids() # Show grids
+    sl.annotate_sphere([90000.0, 90000.0, 90000.0], radius=12000.0,
+                              circle_args={'color':'white', 'linewidth':2})
+    out_name = os.path.join(os.path.dirname(plt), "pulsar_{}.png".format(os.path.basename(plt)))
+    sl.save(out_name)
+
+if __name__ == "__main__":
+    yt.funcs.mylog.setLevel(50)
+    if args.dir:
+        failed = []
+        plotfiles = glob.glob(os.path.join(args.dir, "plt" + "[0-9]"*5))
+        for pf in plotfiles:
+            try:
+                visualize(pf)
+            except:
+                # plotfile 0 may not have particles, so turn off annotations if we first failed
+                try:
+                    visualize(pf, annotate_particles=False)
+                except:
+                    failed.append(pf)
+                    pass
+
+        if len(failed) > 0:
+            print("Visualization failed for the following plotfiles:")
+            for fp in failed:
+                print(fp)
+        else:
+            print("No plots failed, creating a gif ...")
+            input_files = os.path.join(args.dir, "*.png")
+            output_gif  = os.path.join(args.dir, "pulsar.gif")
+            os.system("convert -delay 20 -loop 0 {} {}".format(input_files, output_gif))
+
+    elif args.plotfile:
+        visualize(args.plotfile)
+    else:
+        print("Supply either -f or -d options for visualization. Use the -h option to see help.")

GNUmakefile

@@ -35,6 +35,7 @@ USE_ASCENT_INSITU = FALSE
 USE_OPENPMD = FALSE
 
 WarpxBinDir = Bin
+PULSAR = TRUE
 
 USE_PSATD = FALSE
 USE_PSATD_PICSAR = FALSE

Source/Evolve/WarpXEvolve.cpp

@@ -18,6 +18,9 @@
 #   include "FieldSolver/SpectralSolver/SpectralSolver.H"
 #endif
 
+#ifdef PULSAR
+    #include "Particles/PulsarParameters.H"
+#endif
 #include <cmath>
 #include <limits>
 
@@ -136,7 +139,94 @@ WarpX::Evolve (int numsteps)
             // B : guard cells are NOT up-to-date
         }
 
+#ifdef PULSAR
+        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();
+		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();
+                amrex::GpuArray<int, 3> Ex_stag, Ey_stag, Ez_stag, Bx_stag, By_stag, Bz_stag;
+                for (int idim = 0; idim < 3; ++idim)
+                {
+                    Ex_stag[idim] = ex_type[idim];
+                    Ey_stag[idim] = ey_type[idim];
+                    Ez_stag[idim] = ez_type[idim];
+                    Bx_stag[idim] = bx_type[idim];
+                    By_stag[idim] = by_type[idim];
+                    Bz_stag[idim] = bz_type[idim];
+                }
+                auto geom = Geom(lev).data();
+#ifdef _OPENMP
+#pragma omp parallel if (Gpu::notInLaunchRegion())
+#endif
+                for ( MFIter mfi(*Ex, TilingIfNotGPU()); mfi.isValid(); ++mfi )
+                {
+                    const Box& tex  = mfi.tilebox( Ex->ixType().toIntVect() );
+                    const Box& tey  = mfi.tilebox( Ey->ixType().toIntVect() );
+                    const Box& tez  = mfi.tilebox( Ez->ixType().toIntVect() );
+
+                    auto const& Exfab = Ex->array(mfi);
+                    auto const& Eyfab = Ey->array(mfi);
+                    auto const& Ezfab = Ez->array(mfi);
+
+                    amrex::ParallelFor(tex, tey, tez,
+                    [=] AMREX_GPU_DEVICE (int i, int j, int k)
+                    {
+                        PulsarParm::DampField(i, j, k, geom, Exfab, Ex_stag);
+                    },
+                    [=] AMREX_GPU_DEVICE (int i, int j, int k)
+                    {
+                        PulsarParm::DampField(i, j, k, geom, Eyfab, Ey_stag);
+                    },
+                    [=] AMREX_GPU_DEVICE (int i, int j, int k)
+                    {
+                        PulsarParm::DampField(i, j, k, geom, Ezfab, Ez_stag);
+                    });
+                }
+                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::DampField(i, j, k, geom, Bxfab, Bx_stag);
+                    },
+                    [=] AMREX_GPU_DEVICE (int i, int j, int k)
+                    {
+                        PulsarParm::DampField(i, j, k, geom, Byfab, By_stag);
+                    },
+                    [=] AMREX_GPU_DEVICE (int i, int j, int k)
+                    {
+                        PulsarParm::DampField(i, j, k, geom, Bzfab, Bz_stag);
+                    });
+                }
+            }
+        }
+        FillBoundaryE(guard_cells.ng_FieldSolver, IntVect::TheZeroVector());
+        FillBoundaryB(guard_cells.ng_FieldSolver, IntVect::TheZeroVector());
+#endif
+
+#ifdef WARPX_USE_PY
         if (warpx_py_beforeEsolve) warpx_py_beforeEsolve();
+#endif
 
         if (cur_time + dt[0] >= stop_time - 1.e-3*dt[0] || step == numsteps_max-1) {
             // At the end of last step, push p by 0.5*dt to synchronize
@@ -151,7 +241,9 @@ WarpX::Evolve (int numsteps)
             is_synchronized = true;
         }
 
+#ifdef WARPX_USE_PY
         if (warpx_py_afterEsolve) warpx_py_afterEsolve();
+#endif
 
         for (int lev = 0; lev <= max_level; ++lev) {
             ++istep[lev];
@@ -175,6 +267,16 @@ WarpX::Evolve (int numsteps)
 
         int num_moved = MoveWindow(move_j);
 
+#ifdef PULSAR
+        if (!rho_fp[0]) {
+           amrex::Print() << " no rho -- compute rho! \n";
+        }
+        else {
+           amrex::Print() << " rho is computed \n";
+        }
+        mypc->PulsarParticleRemoval();
+        mypc->PulsarParticleInjection();
+#endif
         mypc->ApplyBoundaryConditions();
 
         // Electrostatic solver: particles can move by an arbitrary number of cells