Merge pull request ECP-WarpX#4 from dwillcox/damp_internal_Efield

Damp the electric field inside the star for the pulsar setup

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.")

Source/Particles/PulsarParameters.H

@@ -124,6 +124,31 @@ namespace PulsarParm
         }
    }
 
+   namespace Spherical
+   {
+        AMREX_GPU_HOST_DEVICE AMREX_INLINE
+        amrex::Real r(int i, int j, int k, amrex::GeometryData const& geom)
+        {
+            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_ilo.x + 0.5) * domain_dx[0];
+            const amrex::Real y = domain_xlo[1] + (j - domain_ilo.y + 0.5) * domain_dx[1];
+            const amrex::Real z = domain_xlo[2] + (k - domain_ilo.z + 0.5) * domain_dx[2];
+
+            const amrex::Real xc = 0.5 * (domain_xlo[0] + domain_xhi[0]);
+            const amrex::Real yc = 0.5 * (domain_xlo[1] + domain_xhi[1]);
+            const amrex::Real zc = 0.5 * (domain_xlo[2] + domain_xhi[2]);
+
+            const amrex::Real r = std::sqrt((x-xc)*(x-xc) + (y-yc)*(y-yc) + (z-zc)*(z-zc));
+
+            return r;
+        }
+   }
+
    AMREX_GPU_HOST_DEVICE AMREX_INLINE
    amrex::Real Omega(const amrex::Real time)
    {
@@ -146,7 +171,7 @@ namespace PulsarParm
            // Fd(0) ~ 0
            // Fd(R_star) ~ 1
 
-           const amrex::Real Fd = std::tanh(damping_scale * (r / R_star - 1)) + 1;
+           const amrex::Real Fd = std::tanh(damping_scale * (r / R_star - 1.0)) + 1.0;
            Efield(i, j, k) = Efield(i, j, k) * Fd;
        }
    }