Seed External Dipole Field in Inspiral

AsterSeeds/interface.ccl

@@ -1,7 +1,7 @@
 # Interface definition for thorn AsterSeeds
 
 IMPLEMENTS: AsterSeeds
-INHERITS: HydroBaseX AsterX
+INHERITS: HydroBaseX AsterX VolumeIntegrals_GRMHDX
 
 USES INCLUDE HEADER: loop_device.hxx
 USES INCLUDE HEADER: setup_eos.hxx
@@ -14,3 +14,31 @@ CCTK_REAL Avec_cent TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"'
   Avec_z_cent
 } "Cell-centered vector potential components"
 
+CCTK_REAL vel_NS1 TYPE=ARRAY DISTRIB=CONSTANT DIM=1 SIZE=3
+{
+  vel_NS1
+} "Helper for dynamic seeding of external dipole field. Records velocity of NS CoM"
+
+CCTK_REAL vel_NS2 TYPE=ARRAY DISTRIB=CONSTANT DIM=1 SIZE=3
+{
+  vel_NS2
+} "Helper for dynamic seeding of external dipole field. Records velocity of NS CoM"
+
+CCTK_INT SeedNow type=SCALAR tags='checkpoint="yes"'
+CCTK_INT DoneSeeding type=SCALAR tags='checkpoint="yes"'
+
+CCTK_INT FUNCTION DriverInterpolate(
+  CCTK_POINTER_TO_CONST IN cctkGH,
+  CCTK_INT IN N_dims,
+  CCTK_INT IN local_interp_handle,
+  CCTK_INT IN param_table_handle,
+  CCTK_INT IN coord_system_handle,
+  CCTK_INT IN N_interp_points,
+  CCTK_INT IN interp_coords_type,
+  CCTK_POINTER_TO_CONST ARRAY IN interp_coords,
+  CCTK_INT IN N_input_arrays,
+  CCTK_INT ARRAY IN input_array_indices,
+  CCTK_INT IN N_output_arrays,
+  CCTK_INT ARRAY IN output_array_types,
+  CCTK_POINTER ARRAY IN output_arrays)
+REQUIRES FUNCTION DriverInterpolate

AsterSeeds/param.ccl

@@ -165,6 +165,7 @@ KEYWORD Afield_config "Definition of the initial vector potential"
   "none"             :: "Nothing set here"
   "internal dipole"  :: "Dipole field according to Ciolfi+2017"
   "external dipole"  :: "Dipole field according to Moesta+2020"
+  "external dipole UIUC"  :: "Dipole field according to Paschalidis+2013"
 } "none"
 
 # parameters for internal dipolar magnetic field
@@ -218,6 +219,65 @@ REAL dipole_z[2] "z-coordinate of the dipole center" STEERABLE=ALWAYS
   *:*     :: "Anything"
 } 0.0
 
+# Initial low density magnetosphere parameters
+
+BOOLEAN initial_beta "Pad pressure to set initial plasma beta floor?" STEERABLE=ALWAYS
+{
+} "no"
+
+REAL initial_beta_min "Initial plasma beta floor" STEERABLE=ALWAYS
+{
+  0:*     :: "Positive"
+} 0.0
+
+REAL initial_beta_rhocut "Pad pressure below this density" STEERABLE=ALWAYS
+{
+  0:*     :: "Positive"
+} 0.0
+
+BOOLEAN set_coorbiting_vel "After padding pressure, set coorbiting velocity? This feature is not well tested, use with caution." STEERABLE=ALWAYS
+{
+} "no"
+
+REAL radius_NS1 "Approximate radius of initial left NS, used only for tapering coorbiting vel." STEERABLE=ALWAYS
+{
+  0:*     :: "Positive"
+} 9.0
+
+REAL radius_NS2 "Approximate radius of initial right NS, used only for tapering coorbiting vel." STEERABLE=ALWAYS
+{
+  0:*     :: "Positive"
+} 9.0
+
+REAL vtol "Overwrite velocity with coorbiting velocity if all |v^i| are below tolerance" STEERABLE=ALWAYS
+{
+  0:*     :: "Positive"
+} 0.0
+
+# BNS inspiral initialization parameters
+
+BOOLEAN use_separation "Seeding time determined by separation of stars from VolumeIntegrals_GRMHDX" STEERABLE=ALWAYS
+{
+} "yes"
+
+BOOLEAN use_time "Seeding time determined by evolution time" STEERABLE=ALWAYS
+{ 
+} "no"
+
+REAL seeding_separation "Coordinate separation of stars below which magnetic field is seeded" STEERABLE=ALWAYS
+{
+  (0.0:*) :: ""
+} 5.0
+
+REAL seeding_time "Time [Cactus units] after which magnetic field is seeded" STEERABLE=ALWAYS
+{
+  [0.0:*) :: ""
+} 0.0
+
+INT seed_every "Check if conditions for seeding are met every <seed_every> iterations" STEERABLE=ALWAYS
+{
+  0:* :: "0 means dont seed in CCTK_POSTSTEP"
+} 0
 
 # ---------------------------------------
 
@@ -311,3 +371,7 @@ USES KEYWORD evolution_eos
 
 SHARES: Con2PrimFactory
 USES CCTK_REAL Ye_atmo
+USES CCTK_REAL rho_abs_min
+USES CCTK_REAL atmo_tol
+USES CCTK_REAL r_atmo
+USES CCTK_REAL n_rho_atmo

AsterSeeds/schedule.ccl

@@ -87,22 +87,148 @@ if (CCTK_Equals(test_type, "3DTest")) {
   if (CCTK_Equals(test_case, "magBNS")) {
 
     STORAGE: Avec_cent
-    SCHEDULE AsterSeeds_InitializeCenteredAvec_BNS AT initial IN ODESolvers_Initial AFTER HydroBaseX_PostInitial BEFORE AsterX_InitialGroup {
-    LANG: C
-    READS: HydroBaseX::press(everywhere)
-    WRITES: Avec_x_cent(everywhere) Avec_y_cent(everywhere) Avec_z_cent(everywhere)
+    SCHEDULE AsterSeeds_Initialize_Seeding_Flags AT CCTK_INITIAL
+    {
+       LANG: C
+       OPTIONS: GLOBAL
+       WRITES:  SeedNow(everywhere), DoneSeeding(everywhere)
+       WRITES: vel_NS1(everywhere) vel_NS2(everywhere)
+    } "Initialize seeding variables."
+
+    if (seed_every != 0) {
+      SCHEDULE AsterSeeds_InitializeAvectoZero AT initial IN ODESolvers_Initial BEFORE AsterX_InitialGroup 
+      {
+        LANG: C
+        WRITES: AsterX::Avec_x(everywhere) AsterX::Avec_y(everywhere) AsterX::Avec_z(everywhere)
+      } "Set up zero vector potential"
+
+      SCHEDULE GROUP Seed_During_Evolution at CCTK_POSTSTEP
+      {
+      } "Routines for seeding magnetic field during evolution."
+
+      SCHEDULE AsterSeeds_Set_Seeding_Flags in Seed_During_Evolution
+      {
+         LANG: C
+         OPTIONS: GLOBAL
+         READS: VolumeIntegrals_GRMHDX::CoM1(everywhere) VolumeIntegrals_GRMHDX::CoM2(everywhere)
+         READS: DoneSeeding(everywhere)
+         WRITES: SeedNow(everywhere)
+      } "Evaluate whether to seed now or not."
+
+      SCHEDULE GROUP Seeding_Tasks IN Seed_During_Evolution IF AsterSeeds::SeedNow AFTER AsterSeeds_Set_Seeding_Flags 
+      {
+      } "If conditions are met, seed magnetic field."
+
+      SCHEDULE AsterSeeds_InitializeCenteredAvec_BNS IN Seeding_Tasks 
+      {
+        LANG: C
+        READS: HydroBaseX::press(everywhere)
+        READS: VolumeIntegrals_GRMHDX::CoM1(everywhere) VolumeIntegrals_GRMHDX::CoM2(everywhere)
+        WRITES: Avec_x_cent(everywhere) Avec_y_cent(everywhere) Avec_z_cent(everywhere)
+        WRITES: SeedNow(everywhere), DoneSeeding(everywhere)
+      } "Set up initial conditions for the cell-centered vector potential"
+
+      SCHEDULE AsterSeeds_InitializeStagAvec_BNS IN Seeding_Tasks AFTER AsterSeeds_InitializeCenteredAvec_BNS BEFORE AsterX_InitialGroup 
+      {
+        LANG: C
+        READS: Avec_x_cent(everywhere) Avec_y_cent(everywhere) Avec_z_cent(everywhere)
+        WRITES: AsterX::Avec_x(interior) AsterX::Avec_y(interior) AsterX::Avec_z(interior)
+        SYNC: AsterX::Avec_x AsterX::Avec_y AsterX::Avec_z
+      } "Set up initial conditions for the vector potential"
+
+      SCHEDULE GROUP AsterX_Con2PrimGroup IN Seeding_Tasks AFTER AsterSeeds_InitializeStagAvec_BNS 
+      {
+      } "Calc B from A, run Con2Prim."
+
+      if (initial_beta) {
+        if (set_coorbiting_vel) {
+          SCHEDULE AsterSeeds_InterpolateNSVelocity IN Seeding_Tasks BEFORE AsterSeeds_SetInitialBetaFloor 
+          {
+            LANG: C
+            OPTIONS: GLOBAL
+            READS: HydroBaseX::vel(everywhere)
+            READS: VolumeIntegrals_GRMHDX::CoM1(everywhere) VolumeIntegrals_GRMHDX::CoM2(everywhere)
+            WRITES: vel_NS1(everywhere) vel_NS2(everywhere)
+          } "Get velocity of NS CoMs"
+        }
+
+        SCHEDULE AsterSeeds_SetInitialBetaFloor IN Seeding_Tasks AFTER AsterX_Con2PrimGroup 
+        {
+          LANG: C
+          READS: HydroBaseX::vel(everywhere) HydroBaseX::Bvec(everywhere) HydroBaseX::press(everywhere) HydroBaseX::temperature(everywhere) HydroBaseX::Ye(everywhere)
+          READS: VolumeIntegrals_GRMHDX::CoM1(everywhere) VolumeIntegrals_GRMHDX::CoM2(everywhere)
+          READS: ADMBaseX::metric(everywhere)
+          READS: vel_NS1(everywhere) vel_NS2(everywhere)
+          WRITES: HydrobaseX::rho(everywhere) HydroBaseX::press(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::entropy(everywhere)  HydrobaseX::vel(everywhere)
+        } "Enforce initial minimum plasma beta"
+
+        SCHEDULE AsterX_Prim2Con_Initial IN Seeding_Tasks AFTER AsterSeeds_SetInitialBetaFloor
+        {
+          LANG: C
+          READS: ADMBaseX::metric(interior)
+          READS: HydroBaseX::rho(interior) HydroBaseX::vel(interior) HydroBaseX::eps(interior) HydroBaseX::press(interior) HydroBaseX::Bvec(interior)
+          READS: HydroBaseX::entropy(interior) HydroBaseX::Ye(interior)
+          WRITES: AsterX::cons_vector(interior) AsterX::dB(interior)
+          SYNC: AsterX::cons_vector AsterX::dB
+        } "Update cons after setting initial beta"
+      }
+
+      SCHEDULE AsterX_Tmunu IN Seeding_Tasks AFTER AsterX_Con2PrimGroup AFTER AsterX_ApplyOuterBCOnPrim
+      {
+        LANG: C
+        READS: ADMBaseX::metric(everywhere) ADMBaseX::lapse(everywhere) ADMBaseX::shift(everywhere)
+        READS: HydroBaseX::rho(everywhere) HydroBaseX::press(everywhere) HydroBaseX::eps(everywhere)
+        READS: HydroBaseX::vel(everywhere)
+        READS: AsterX::zvec_x(everywhere), AsterX::zvec_y(everywhere), AsterX::zvec_z(everywhere)
+        READS: AsterX::svec_x(everywhere), AsterX::svec_y(everywhere), AsterX::svec_z(everywhere)
+        READS: HydroBaseX::Bvec(everywhere)
+        READS: TmunuBaseX::eTtt(interior) TmunuBaseX::eTti(interior) TmunuBaseX::eTij(interior)
+        WRITES: TmunuBaseX::eTtt(interior) TmunuBaseX::eTti(interior) TmunuBaseX::eTij(interior)
+        SYNC: TmunuBaseX::eTtt TmunuBaseX::eTti TmunuBaseX::eTij
+      } "Compute the energy-momentum tensor"
+
     }
-    "Set up initial conditions for the cell-centered vector potential"
+    else {
+      SCHEDULE AsterSeeds_InitializeCenteredAvec_BNS AT initial IN ODESolvers_Initial AFTER HydroBaseX_PostInitial BEFORE AsterX_InitialGroup {
+        LANG: C
+        READS: HydroBaseX::press(everywhere)
+        READS: VolumeIntegrals_GRMHDX::CoM1(everywhere) VolumeIntegrals_GRMHDX::CoM2(everywhere)
+        WRITES: Avec_x_cent(everywhere) Avec_y_cent(everywhere) Avec_z_cent(everywhere)
+        WRITES:  SeedNow(everywhere), DoneSeeding(everywhere)
+      } "Set up initial conditions for the cell-centered vector potential"
 
-    SCHEDULE AsterSeeds_InitializeStagAvec_BNS AT initial IN ODESolvers_Initial AFTER AsterSeeds_InitializeCenteredAvec_BNS BEFORE AsterX_InitialGroup {
-    LANG: C
-    READS: Avec_x_cent(everywhere) Avec_y_cent(everywhere) Avec_z_cent(everywhere)
-    WRITES: AsterX::Avec_x(interior) AsterX::Avec_y(interior) AsterX::Avec_z(interior)
-    SYNC: AsterX::Avec_x AsterX::Avec_y AsterX::Avec_z
+      SCHEDULE AsterSeeds_InitializeStagAvec_BNS AT initial IN ODESolvers_Initial AFTER AsterSeeds_InitializeCenteredAvec_BNS BEFORE AsterX_InitialGroup 
+      {
+        LANG: C
+        READS: Avec_x_cent(everywhere) Avec_y_cent(everywhere) Avec_z_cent(everywhere)
+        WRITES: AsterX::Avec_x(interior) AsterX::Avec_y(interior) AsterX::Avec_z(interior)
+        SYNC: AsterX::Avec_x AsterX::Avec_y AsterX::Avec_z
+      } "Set up initial conditions for the vector potential"
+
+      if (initial_beta) {
+        if (set_coorbiting_vel) {
+          SCHEDULE AsterSeeds_InterpolateNSVelocity IN Seeding_Tasks BEFORE AsterSeeds_SetInitialBetaFloor 
+          {
+            LANG: C
+            OPTIONS: GLOBAL
+            READS: HydroBaseX::vel(everywhere)
+            READS: VolumeIntegrals_GRMHDX::CoM1(everywhere) VolumeIntegrals_GRMHDX::CoM2(everywhere)
+            WRITES: vel_NS1(everywhere) vel_NS2(everywhere)
+          } "Get velocity of NS CoMs"
+        }
+
+        SCHEDULE AsterSeeds_SetInitialBetaFloor AT initial IN AsterX_InitialGroup AFTER AsterX_ComputeBFromdB BEFORE AsterX_CheckPrims 
+        {
+          LANG: C
+          READS: HydroBaseX::vel(everywhere) HydroBaseX::Bvec(everywhere) HydroBaseX::press(everywhere) HydroBaseX::temperature(everywhere) HydroBaseX::Ye(everywhere)
+          READS: ADMBaseX::metric(everywhere)
+          READS: VolumeIntegrals_GRMHDX::CoM1(everywhere) VolumeIntegrals_GRMHDX::CoM2(everywhere)
+          READS: vel_NS1(everywhere) vel_NS2(everywhere)
+          WRITES: HydrobaseX::rho(everywhere) HydroBaseX::press(everywhere) HydroBaseX::eps(everywhere) HydroBaseX::entropy(everywhere) HydrobaseX::vel(everywhere)
+        } "Enforce initial minimum plasma beta"
+      }
     }
-    "Set up initial conditions for the vector potential"
   }
-
 }
 
 #Initial conditions for TabEOS tests