From e43f9b8b52738f683eec0e5a04b793e3c777f938 Mon Sep 17 00:00:00 2001
From: RevathiJambunathan <revanathan@gmail.com>
Date: Fri, 15 Sep 2023 17:48:16 -0700
Subject: [PATCH] pml cubic sigma

---
 Source/BoundaryConditions/PML.cpp     | 32 +++++++++++++++++++++++++++
 Source/Particles/PulsarParameters.H   |  1 +
 Source/Particles/PulsarParameters.cpp | 10 +++++----
 3 files changed, 39 insertions(+), 4 deletions(-)

diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp
index 0739175cbc2..fb2354b38ed 100644
--- a/Source/BoundaryConditions/PML.cpp
+++ b/Source/BoundaryConditions/PML.cpp
@@ -19,6 +19,7 @@
 #include "Utils/WarpXProfilerWrapper.H"
 #include "Utils/Parser/ParserUtils.H"
 #include "WarpX.H"
+#include "Particles/PulsarParameters.H"
 
 #include <ablastr/utils/Communication.H>
 
@@ -71,6 +72,9 @@ namespace
         Real* p_sigma_cumsum = sigma_cumsum.data();
         Real* p_sigma_star = sigma_star.data();
         Real* p_sigma_star_cumsum = sigma_star_cumsum.data();
+#ifdef PULSAR
+        int pml_cubic_sigma = Pulsar::m_pml_cubic_sigma;
+#endif
         amrex::ParallelFor(N, [=] AMREX_GPU_DEVICE (int i) noexcept
         {
             i += olo;
@@ -78,11 +82,24 @@ namespace
             p_sigma[i-slo] = fac*(offset*offset);
             // sigma_cumsum is the analytical integral of sigma function at same points than sigma
             p_sigma_cumsum[i-slo] = (fac*(offset*offset*offset)/3._rt)/v_sigma;
+#ifdef PULSAR
+            if (pml_cubic_sigma == 1) {
+                p_sigma[i-slo] = fac*(offset*offset*offset);
+                // sigma_cumsum is the analytical integral of sigma function at same points than sigma
+                p_sigma_cumsum[i-slo] = (fac*(offset*offset*offset*offset)/4._rt)/v_sigma;
+            }
+#endif
             if (i <= ohi+1) {
                 offset = static_cast<Real>(glo-i) - 0.5_rt;
                 p_sigma_star[i-sslo] = fac*(offset*offset);
                 // sigma_star_cumsum is the analytical integral of sigma function at same points than sigma_star
                 p_sigma_star_cumsum[i-sslo] = (fac*(offset*offset*offset)/3._rt)/v_sigma;
+#ifdef PULSAR
+                if (pml_cubic_sigma == 1) {
+                    p_sigma_star[i-sslo] = fac*(offset*offset*offset);
+                    p_sigma_star_cumsum[i-sslo] = (fac*(offset*offset*offset*offset)/4._rt)/v_sigma;
+                }
+#endif
             }
         });
     }
@@ -100,16 +117,31 @@ namespace
         Real* p_sigma_cumsum = sigma_cumsum.data();
         Real* p_sigma_star = sigma_star.data();
         Real* p_sigma_star_cumsum = sigma_star_cumsum.data();
+#ifdef PULSAR
+        int pml_cubic_sigma = Pulsar::m_pml_cubic_sigma;
+#endif
         amrex::ParallelFor(N, [=] AMREX_GPU_DEVICE (int i) noexcept
         {
             i += olo;
             Real offset = static_cast<Real>(i-ghi-1);
             p_sigma[i-slo] = fac*(offset*offset);
             p_sigma_cumsum[i-slo] = (fac*(offset*offset*offset)/3._rt)/v_sigma;
+#ifdef PULSAR
+            if (pml_cubic_sigma == 1) {
+                p_sigma[i-slo] = fac*(offset*offset*offset);
+                p_sigma_cumsum[i-slo] = (fac*(offset*offset*offset*offset)/4._rt)/v_sigma;
+            }
+#endif
             if (i <= ohi+1) {
                 offset = static_cast<Real>(i-ghi) - 0.5_rt;
                 p_sigma_star[i-sslo] = fac*(offset*offset);
                 p_sigma_star_cumsum[i-sslo] = (fac*(offset*offset*offset)/3._rt)/v_sigma;
+#ifdef PULSAR
+            if (pml_cubic_sigma == 1) {
+                p_sigma_star[i-sslo] = fac*(offset*offset*offset);
+                p_sigma_star_cumsum[i-sslo] = (fac*(offset*offset*offset*offset)/4._rt)/v_sigma;
+            }
+#endif
             }
         });
     }
diff --git a/Source/Particles/PulsarParameters.H b/Source/Particles/PulsarParameters.H
index c4fc6d1d839..a208541f787 100644
--- a/Source/Particles/PulsarParameters.H
+++ b/Source/Particles/PulsarParameters.H
@@ -701,6 +701,7 @@ public:
     static amrex::Real m_gatherbuffer_max;
     static amrex::Real m_depositbuffer_min;
     static amrex::Real m_depositbuffer_max;
+    static int m_pml_cubic_sigma;
 private:
 };
 
diff --git a/Source/Particles/PulsarParameters.cpp b/Source/Particles/PulsarParameters.cpp
index f37e394a0c9..7408e75cf4f 100644
--- a/Source/Particles/PulsarParameters.cpp
+++ b/Source/Particles/PulsarParameters.cpp
@@ -138,10 +138,11 @@ amrex::Real Pulsar::m_limit_GJfactor = 1.;
 int Pulsar::m_use_BC_smoothening;
 amrex::Real Pulsar::m_min_BC_radius;
 amrex::Real Pulsar::m_BC_width;
-amrex::Real Pulsar::m_gatherbuffer_min;
-amrex::Real Pulsar::m_gatherbuffer_max;
-amrex::Real Pulsar::m_depositbuffer_min;
-amrex::Real Pulsar::m_depositbuffer_max;
+amrex::Real Pulsar::m_gatherbuffer_min = 0.;
+amrex::Real Pulsar::m_gatherbuffer_max = 0.;
+amrex::Real Pulsar::m_depositbuffer_min = 0.;
+amrex::Real Pulsar::m_depositbuffer_max = 0.;
+int Pulsar::m_pml_cubic_sigma;
 
 Pulsar::Pulsar ()
 {
@@ -372,6 +373,7 @@ Pulsar::ReadParameters () {
     pp.query("gatherbuffer_max",m_gatherbuffer_max);
     pp.query("depositbuffer_min",m_depositbuffer_min);
     pp.query("depositbuffer_max",m_depositbuffer_max);
+    pp.get("pml_cubic_sigma", m_pml_cubic_sigma);
 }