AsterSeeds: add code to perturb neutron star

AsterSeeds/param.ccl

@@ -42,6 +42,28 @@ KEYWORD test_case "Name of the testcase" STEERABLE=never
 
 private:
 
+# parameters for perturbing neutron-star initial data
+
+KEYWORD ns_perturbation_mode "Neutron-star perturbation applied after hydro initial data are read" STEERABLE=never
+{
+  "none" :: "Do not perturb the neutron-star initial data"
+  "m0"   :: "Deplete pressure by ns_perturbation_strength inside the star"
+  "m1"   :: "rho -> rho * (1 + B * (x + y) / Req)"
+  "m2"   :: "rho -> rho * (1 + B * (x^2 - y^2) / Req^2)"
+} "none"
+
+CCTK_REAL ns_perturbation_strength "Fractional perturbation strength B"
+{
+  *:* :: ""
+} 0.0
+
+CCTK_REAL ns_perturbation_Req "Equatorial radius Req used for m=1 and m=2 perturbations"
+{
+  (0.0:* :: "Positive"
+} 1.0
+
+# ---------------------------------------
+
 # parameters for rotating intial data
 REAL rotate_angle_x "Rotation angle about x-axis (in units of pi)" STEERABLE=never
 {

AsterSeeds/schedule.ccl

@@ -252,7 +252,7 @@ if (CCTK_Equals(test_type, "TabEOSTest")) {
 
 #Initial conditions for temperature, entropy and Ye
 
-if (set_Ye_postinitial) {
+if (set_Ye_postinitial && !CCTK_Equals(evolution_eos, "Tabulated3d")) {
   SCHEDULE SetYe IN HydroBaseX_PostInitial
   {
     LANG: C
@@ -261,22 +261,89 @@ if (set_Ye_postinitial) {
   } "Set initial Ye"
 }
 
-if (set_temperature_postinitial) {
-  SCHEDULE SetTemp IN HydroBaseX_PostInitial AFTER SetYe
-  {
-    LANG: C 
-    READS: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::Ye(everywhere)
-    WRITES: HydroBaseX::temperature(everywhere) HydroBaseX::eps(everywhere)
-
-  } "Set initial temperature"
+if (set_temperature_postinitial && !CCTK_Equals(evolution_eos, "Tabulated3d")) {
+  if (set_Ye_postinitial) {
+    SCHEDULE SetTemp IN HydroBaseX_PostInitial AFTER SetYe
+    {
+      LANG: C
+      READS: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::Ye(everywhere)
+      WRITES: HydroBaseX::temperature(everywhere) HydroBaseX::eps(everywhere)
+
+    } "Set initial temperature"
+  } else {
+    SCHEDULE SetTemp IN HydroBaseX_PostInitial
+    {
+      LANG: C
+      READS: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::Ye(everywhere)
+      WRITES: HydroBaseX::temperature(everywhere) HydroBaseX::eps(everywhere)
+
+    } "Set initial temperature"
+  }
 }
 
-if (set_entropy_postinitial) {
-  SCHEDULE SetEntropy IN HydroBaseX_PostInitial AFTER SetTemp
-  {
-    LANG: C
-    READS: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::Ye(everywhere)
-    WRITES: HydroBaseX::entropy(everywhere) HydroBaseX::eps(everywhere)
-
-  } "Set initial entropy"
+if (set_entropy_postinitial && !CCTK_Equals(evolution_eos, "Tabulated3d")) {
+  if (set_temperature_postinitial) {
+    SCHEDULE SetEntropy IN HydroBaseX_PostInitial AFTER SetTemp
+    {
+      LANG: C
+      READS: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::Ye(everywhere)
+      WRITES: HydroBaseX::entropy(everywhere) HydroBaseX::eps(everywhere)
+
+    } "Set initial entropy"
+  } else if (set_Ye_postinitial) {
+    SCHEDULE SetEntropy IN HydroBaseX_PostInitial AFTER SetYe
+    {
+      LANG: C
+      READS: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::Ye(everywhere)
+      WRITES: HydroBaseX::entropy(everywhere) HydroBaseX::eps(everywhere)
+
+    } "Set initial entropy"
+  } else {
+    SCHEDULE SetEntropy IN HydroBaseX_PostInitial
+    {
+      LANG: C
+      READS: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::Ye(everywhere)
+      WRITES: HydroBaseX::entropy(everywhere) HydroBaseX::eps(everywhere)
+
+    } "Set initial entropy"
+  }
+}
+
+if (!CCTK_Equals(ns_perturbation_mode, "none")) {
+  if (CCTK_Equals(evolution_eos, "Tabulated3d")) {
+    SCHEDULE AsterSeeds_PerturbNS IN HydroBaseX_PostInitial AFTER ID_TabEOS_HydroQuantities_recompute_HydroBase_variables
+    {
+      LANG: C
+      READS: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::press(everywhere) HydroBaseX::temperature(everywhere) HydroBaseX::Ye(everywhere)
+      WRITES: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::press(everywhere) HydroBaseX::temperature(everywhere) HydroBaseX::entropy(everywhere)
+    } "Apply neutron-star pressure or density perturbation after tabulated-EOS hydro quantities"
+  } else if (set_entropy_postinitial) {
+    SCHEDULE AsterSeeds_PerturbNS IN HydroBaseX_PostInitial AFTER SetEntropy
+    {
+      LANG: C
+      READS: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::press(everywhere) HydroBaseX::temperature(everywhere) HydroBaseX::Ye(everywhere)
+      WRITES: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::press(everywhere) HydroBaseX::temperature(everywhere) HydroBaseX::entropy(everywhere)
+    } "Apply neutron-star pressure or density perturbation after hydro initial data"
+  } else if (set_temperature_postinitial) {
+    SCHEDULE AsterSeeds_PerturbNS IN HydroBaseX_PostInitial AFTER SetTemp
+    {
+      LANG: C
+      READS: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::press(everywhere) HydroBaseX::temperature(everywhere) HydroBaseX::Ye(everywhere)
+      WRITES: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::press(everywhere) HydroBaseX::temperature(everywhere) HydroBaseX::entropy(everywhere)
+    } "Apply neutron-star pressure or density perturbation after hydro initial data"
+  } else if (set_Ye_postinitial) {
+    SCHEDULE AsterSeeds_PerturbNS IN HydroBaseX_PostInitial AFTER SetYe
+    {
+      LANG: C
+      READS: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::press(everywhere) HydroBaseX::temperature(everywhere) HydroBaseX::Ye(everywhere)
+      WRITES: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::press(everywhere) HydroBaseX::temperature(everywhere) HydroBaseX::entropy(everywhere)
+    } "Apply neutron-star pressure or density perturbation after hydro initial data"
+  } else {
+    SCHEDULE AsterSeeds_PerturbNS IN HydroBaseX_PostInitial
+    {
+      LANG: C
+      READS: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::press(everywhere) HydroBaseX::temperature(everywhere) HydroBaseX::Ye(everywhere)
+      WRITES: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::press(everywhere) HydroBaseX::temperature(everywhere) HydroBaseX::entropy(everywhere)
+    } "Apply neutron-star pressure or density perturbation after hydro initial data"
+  }
 }

AsterSeeds/src/SetTempEntropyYe.cxx

@@ -1,44 +1,48 @@
-#include <loop.hxx>
 #include <loop_device.hxx>
 
 #include <cctk.h>
 #include <cctk_Arguments.h>
 #include <cctk_Parameters.h>
 
-#include <cstdio>
-#include <cstdbool>
-#include <cmath>
+#include <cassert>
 
 #include <setup_eos.hxx>
 
+namespace AsterSeeds {
+using namespace Loop;
 using namespace EOSX;
+
 enum class eos_3param { IdealGas, Hybrid, Tabulated };
 
 template <typename EOSType>
-void SetTemp_typeEoS(CCTK_ARGUMENTS, EOSType *eos_3p) {
+void AsterSeeds_SetTemp_typeEoS(CCTK_ARGUMENTS, EOSType *eos_3p) {
 
   DECLARE_CCTK_ARGUMENTSX_SetTemp;
   DECLARE_CCTK_PARAMETERS;
 
   grid.loop_all_device<1, 1, 1>(
       grid.nghostzones,
-      [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
+      [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
+        CCTK_REAL epsL = eps(p.I);
         temperature(p.I) =
-            eos_3p->temp_from_rho_eps_ye(rho(p.I), eps(p.I), Ye(p.I));
+            eos_3p->temp_from_rho_eps_ye(rho(p.I), epsL, Ye(p.I));
+        eps(p.I) = epsL;
       });
 }
 
 template <typename EOSType>
-void SetEntropy_typeEoS(CCTK_ARGUMENTS, EOSType *eos_3p) {
+void AsterSeeds_SetEntropy_typeEoS(CCTK_ARGUMENTS, EOSType *eos_3p) {
 
   DECLARE_CCTK_ARGUMENTSX_SetEntropy;
   DECLARE_CCTK_PARAMETERS;
 
   grid.loop_all_device<1, 1, 1>(
       grid.nghostzones,
-      [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
+      [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
+        CCTK_REAL epsL = eps(p.I);
         entropy(p.I) =
-            eos_3p->kappa_from_rho_eps_ye(rho(p.I), eps(p.I), Ye(p.I));
+            eos_3p->kappa_from_rho_eps_ye(rho(p.I), epsL, Ye(p.I));
+        eps(p.I) = epsL;
       });
 }
 
@@ -74,16 +78,16 @@ extern "C" void SetTemp(CCTK_ARGUMENTS) {
   switch (eos_3p_type) {
   case eos_3param::IdealGas: {
     auto eos_3p_ig = global_eos_3p_ig;
-    SetTemp_typeEoS(CCTK_PASS_CTOC, eos_3p_ig);
+    AsterSeeds_SetTemp_typeEoS(CCTK_PASS_CTOC, eos_3p_ig);
     break;
   }
   case eos_3param::Hybrid: {
-    if (global_eos_3p_hyb_poly) {
-      auto eos_3p_hyb = global_eos_3p_hyb_poly;
-      SetTemp_typeEoS(CCTK_PASS_CTOC, eos_3p_hyb);
-    } else if (global_eos_3p_hyb_pwpoly) {
+    if (global_eos_3p_hyb_pwpoly) {
       auto eos_3p_hyb = global_eos_3p_hyb_pwpoly;
-      SetTemp_typeEoS(CCTK_PASS_CTOC, eos_3p_hyb);
+      AsterSeeds_SetTemp_typeEoS(CCTK_PASS_CTOC, eos_3p_hyb);
+    } else if (global_eos_3p_hyb_poly) {
+      auto eos_3p_hyb = global_eos_3p_hyb_poly;
+      AsterSeeds_SetTemp_typeEoS(CCTK_PASS_CTOC, eos_3p_hyb);
     } else {
       CCTK_ERROR(
           "Hybrid EOS selected but no hybrid EOS object was initialized");
@@ -92,7 +96,7 @@ extern "C" void SetTemp(CCTK_ARGUMENTS) {
   }
   case eos_3param::Tabulated: {
     auto eos_3p_tab3d = global_eos_3p_tab3d;
-    SetTemp_typeEoS(CCTK_PASS_CTOC, eos_3p_tab3d);
+    AsterSeeds_SetTemp_typeEoS(CCTK_PASS_CTOC, eos_3p_tab3d);
     break;
   }
   default:
@@ -121,16 +125,16 @@ extern "C" void SetEntropy(CCTK_ARGUMENTS) {
   switch (eos_3p_type) {
   case eos_3param::IdealGas: {
     auto eos_3p_ig = global_eos_3p_ig;
-    SetEntropy_typeEoS(CCTK_PASS_CTOC, eos_3p_ig);
+    AsterSeeds_SetEntropy_typeEoS(CCTK_PASS_CTOC, eos_3p_ig);
     break;
   }
   case eos_3param::Hybrid: {
-    if (global_eos_3p_hyb_poly) {
-      auto eos_3p_hyb = global_eos_3p_hyb_poly;
-      SetEntropy_typeEoS(CCTK_PASS_CTOC, eos_3p_hyb);
-    } else if (global_eos_3p_hyb_pwpoly) {
+    if (global_eos_3p_hyb_pwpoly) {
       auto eos_3p_hyb = global_eos_3p_hyb_pwpoly;
-      SetEntropy_typeEoS(CCTK_PASS_CTOC, eos_3p_hyb);
+      AsterSeeds_SetEntropy_typeEoS(CCTK_PASS_CTOC, eos_3p_hyb);
+    } else if (global_eos_3p_hyb_poly) {
+      auto eos_3p_hyb = global_eos_3p_hyb_poly;
+      AsterSeeds_SetEntropy_typeEoS(CCTK_PASS_CTOC, eos_3p_hyb);
     } else {
       CCTK_ERROR(
           "Hybrid EOS selected but no hybrid EOS object was initialized");
@@ -139,10 +143,12 @@ extern "C" void SetEntropy(CCTK_ARGUMENTS) {
   }
   case eos_3param::Tabulated: {
     auto eos_3p_tab3d = global_eos_3p_tab3d;
-    SetEntropy_typeEoS(CCTK_PASS_CTOC, eos_3p_tab3d);
+    AsterSeeds_SetEntropy_typeEoS(CCTK_PASS_CTOC, eos_3p_tab3d);
     break;
   }
   default:
     assert(0);
   }
 }
+
+} // namespace AsterSeeds

AsterSeeds/src/magTOV.cxx

@@ -52,6 +52,34 @@ extern "C" void AsterSeeds_InitializeCenteredAvec_TOV(CCTK_ARGUMENTS) {
           Avec_z_cent(p.I) = 0.0;
         });
 
+  } else if (CCTK_EQUALS(Afield_config, "external dipole UIUC")) {
+
+    /* computing cell centered vector potential components */
+    grid.loop_all<1, 1, 1>(
+        grid.nghostzones,
+        [=] CCTK_HOST(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
+          const CCTK_REAL pi = 2 * acos(0.0);
+
+          CCTK_REAL x_local = p.x - dipole_x[0];
+          CCTK_REAL y_local = p.y - dipole_y[0];
+          CCTK_REAL z_local = p.z - dipole_z[0];
+          CCTK_REAL cylrad2 = x_local * x_local + y_local * y_local;
+          CCTK_REAL sphrad2 =
+              x_local * x_local + y_local * y_local + z_local * z_local;
+          CCTK_REAL r02 = r0 * r0;
+
+          // See e.g. Ruiz+ 2020, Eq. 1. Here, B0 is I0 from Eq. 1, and we
+          // have preemptively canceled out the factor of cylrad2.
+          CCTK_REAL Aphi_local =
+              pi * B0 * r02 / pow(r02 + sphrad2, 1.5) *
+              (1.0 + (15.0 * r02 * (r02 + cylrad2) /
+                      (8.0 * pow(r02 + sphrad2, 2.0))));
+
+          Avec_x_cent(p.I) = -y_local * Aphi_local;
+          Avec_y_cent(p.I) = x_local * Aphi_local;
+          Avec_z_cent(p.I) = 0.0;
+        });
+
   } else {
     CCTK_ERROR("Vector potential configuration not defined");
   }

AsterSeeds/src/make.code.defn

@@ -9,6 +9,7 @@ SRCS =	\
   atmosphere.cxx \
   bbhcloud.cxx \
   initial_betafloor.cxx \
+  perturbNS.cxx \
   magTOV.cxx \
   magBNS.cxx \
   SetTempEntropyYe.cxx