Development

AsterAnalysis/interface.ccl

@@ -8,14 +8,10 @@ USES INCLUDE HEADER: aster_fd.hxx aster_interp.hxx aster_utils.hxx
 
 CCTK_REAL mhd_analysis TYPE=gf CENTERING={ccc} tags='checkpoint="no"'
 {
-divA, A_norm, divB, B_norm, b2small, mean_B_numerator, mean_B_denominator, mhd_press_ratio, mhd_energy_ratio, magnetic_energy, kinetic_energy, poynting_flux_x, poynting_flux_y, poynting_flux_z, poynting_flux_norm, poynting_vector_x, poynting_vector_y, poynting_vector_z, poynting_scalar, lambda_MRI_x, lambda_MRI_y, lambda_MRI_z
+divA, A_norm, divB, mean_B_numerator, mean_B_denominator, mhd_press_ratio, mhd_energy_ratio, magnetic_energy, kinetic_energy, poynting_flux_x, poynting_flux_y, poynting_flux_z, poynting_flux_norm, poynting_vector_x, poynting_vector_y, poynting_vector_z, poynting_scalar, lambda_MRI_x, lambda_MRI_y, lambda_MRI_z
 } "Variables used for analysis purposes when magnetic fields are active"
 
 CCTK_REAL mhd_spherical TYPE=gf CENTERING={ccc} tags='checkpoint="no"'
 {
 Brho, Btheta, Bphi, B_Poloidal_norm, B_Toroidal_norm, mean_B_Tor_numerator, mean_B_Tor_denominator, mean_B_Pol_numerator
 } "For decomposition of B field into toroidal (phi) and poloidal parts"
-
-public:
-
-CCTK_REAL w_lorentz TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' "Lorentz factor"

AsterAnalysis/schedule.ccl

@@ -13,15 +13,16 @@ if (compute_MHD) {
     READS: HydroBaseX::vel(everywhere)
     READS: HydroBaseX::Bvec(everywhere)
     READS: AsterX::zvec(everywhere)
+    READS: AsterX::B_norm(everywhere)
+    READS: AsterX::b2small(everywhere)
+    READS: AsterX::w_lorentz(everywhere)
     READS: AsterX::Avec_x(everywhere)
     READS: AsterX::Avec_y(everywhere)
     READS: AsterX::Avec_z(everywhere)
     READS: AsterX::dBx_stag(everywhere)
     READS: AsterX::dBy_stag(everywhere)
     READS: AsterX::dBz_stag(everywhere)
-    WRITES: w_lorentz(interior)
     WRITES: mhd_analysis(interior)
-    SYNC: w_lorentz
   }
   "Compute MHD analysis quantities"
 }
@@ -30,7 +31,7 @@ if (compute_MHD_spherical) {
   SCHEDULE AsterAnalysis_MHD_Spherical AT analysis AFTER AsterAnalysis_MHD
   {
     LANG: C
-    READS: w_lorentz(interior)
+    READS: AsterX::w_lorentz(interior)
     READS: ADMBaseX::metric(interior)
     READS: HydroBaseX::Bvec(interior)
     READS: HydroBaseX::rho(interior)

AsterAnalysis/src/analysis_mhd.cxx

@@ -56,13 +56,11 @@ extern "C" void AsterAnalysis_MHD(CCTK_ARGUMENTS) {
           v_up(1) = vely(p.I);
           v_up(2) = velz(p.I);
           v_low = calc_contraction(g_avg, v_up);
-          w_lorentz(p.I) = calc_wlorentz(v_low, v_up);
 
         } else {
 
           const vec<CCTK_REAL, 3> z_up{zvec_x(p.I), zvec_y(p.I), zvec_z(p.I)};
           const vec<CCTK_REAL, 3> z_low = calc_contraction(g_avg, z_up);
-          w_lorentz(p.I) = calc_wlorentz_zvec(z_low, z_up);
           v_up = z_up / w_lorentz(p.I);
           v_low = z_low / w_lorentz(p.I);
         }
@@ -76,11 +74,12 @@ extern "C" void AsterAnalysis_MHD(CCTK_ARGUMENTS) {
         const vec<CCTK_REAL, 3> B_up{Bvecx(p.I), Bvecy(p.I), Bvecz(p.I)};
         const vec<CCTK_REAL, 3> B_low = calc_contraction(g_avg, B_up);
         const CCTK_REAL B2big = calc_contraction(B_up, B_low);
-        B_norm(p.I) = sqrt(B2big);
 
         /* cell-centered A^i and A_i */
-        const vec<CCTK_REAL, 3> A_low(
-            [&](int i) ARITH_INLINE { return calc_avg_e2c(gf_Avec(i), p, i); });
+        const vec<CCTK_REAL, 3> A_low{calc_avg_e2c<0>(gf_Avec(0), p),
+                                      calc_avg_e2c<1>(gf_Avec(1), p),
+                                      calc_avg_e2c<2>(gf_Avec(2), p)};
+
         const vec<CCTK_REAL, 3> A_up = calc_contraction(ug_avg, A_low);
 
         A_norm(p.I) = sqrt(calc_contraction(A_low, A_up));
@@ -91,10 +90,6 @@ extern "C" void AsterAnalysis_MHD(CCTK_ARGUMENTS) {
         const vec<CCTK_REAL, 3> bs =
             (B_up / w_lorentz(p.I) + alp_avg * bs0 * velshift);
 
-        /* b^2 */
-        b2small(p.I) = (calc_contraction(B_up, B_low) + pow2(alp_avg * bs0)) /
-                       pow2(w_lorentz(p.I));
-
         /* inverse Beta and magnetization */
         mhd_press_ratio(p.I) = b2small(p.I) / (2.0 * press(p.I));
         mhd_energy_ratio(p.I) = b2small(p.I) / (2.0 * rho(p.I));

AsterMasks/interface.ccl

@@ -8,3 +8,10 @@ PUBLIC:
 
 CCTK_REAL aster_mask_vc CENTERING={vvv} TYPE=gf TAGS='checkpoint="no"' "Vertex-centered mask"
 CCTK_REAL aster_mask_cc CENTERING={ccc} TYPE=gf TAGS='checkpoint="no"' "Cell-centered mask"
+
+# This scalar can be used by external thorns to indicate 
+# whether we reset the gridfunctions above as 
+# AsterMasks_SetToOne/PostStep is always scheduled. 
+# We set the gfs to one if: MaskBool = 1 
+# MaskBool = 0 by default most of the time
+CCTK_REAL MaskBool TYPE=scalar "Whether we set the mask gridfunctions to one"

AsterMasks/schedule.ccl

@@ -1,21 +1,58 @@
 #Schedule definitions for thorn AsterMasks
 
 STORAGE: aster_mask_vc
+STORAGE: aster_mask_cc
 
-SCHEDULE AsterMasks_Initialize IN CCTK_INITIAL BEFORE ODESolvers_Initial {
+# ----- CCTK_INITIAL
+
+SCHEDULE AsterMasks_SetToOne IN CCTK_INITIAL BEFORE ODESolvers_Initial {
+LANG: C
+WRITES:  aster_mask_vc(everywhere) aster_mask_cc(everywhere)
+} "Set masks = 1"
+
+SCHEDULE AsterMasks_SetMaskBoolZero IN CCTK_INITIAL BEFORE ODESolvers_Initial AFTER AsterMasks_Init {
+LANG: C
+OPTIONS: global
+WRITES:  MaskBool
+} "Set MaskBool = 0"
+
+# ----- ODESolvers_PostStep
+
+SCHEDULE AsterMasks_SetToOnePostStep IN ODESolvers_PostStep BEFORE ADMBaseX_SetADMVars {
+LANG: C
+READS:  MaskBool
+READS:  aster_mask_vc(everywhere) aster_mask_cc(everywhere)
+WRITES: aster_mask_vc(everywhere) aster_mask_cc(everywhere)
+} "If (MaskBool = 1) masks = 1"
+
+SCHEDULE AsterMasks_SetMaskBoolZero IN ODESolvers_PostStep BEFORE ADMBaseX_SetADMVars AFTER AsterMasks_SetToOnePostStep {
+LANG: C
+OPTIONS: global
+WRITES:  MaskBool
+} "Set MaskBool = 0"
+
+# ----- CCTK_POSTREGRID
+
+SCHEDULE AsterMasks_SetToOne IN CCTK_POSTREGRID BEFORE ODESolvers_PostStep {
 LANG: C
-WRITES: aster_mask_vc(everywhere)
-WRITES: aster_mask_cc(everywhere)
-} "Initialize VC and CC mask to true"
+WRITES: aster_mask_vc(everywhere) aster_mask_cc(everywhere)
+} "Set masks = 1"
+
+SCHEDULE AsterMasks_SetMaskBoolZero IN CCTK_POSTREGRID BEFORE ODESolvers_PostStep AFTER AsterMasks_SetToOne {
+LANG: C
+OPTIONS: global
+WRITES:  MaskBool
+} "Set MaskBool = 0"
+
+# ----- CCTK_POST_RECOVER_VARIABLES
 
-SCHEDULE AsterMasks_Initialize IN CCTK_POSTREGRID BEFORE ODESolvers_PostStep {
+SCHEDULE AsterMasks_SetToOne IN CCTK_POST_RECOVER_VARIABLES BEFORE ODESolvers_PostStep {
 LANG: C
-WRITES: aster_mask_vc(everywhere)
-WRITES: aster_mask_cc(everywhere)
-} "Initialize VC and CC mask to true"
+WRITES: aster_mask_vc(everywhere) aster_mask_cc(everywhere)
+} "Set masks = 1"
 
-SCHEDULE AsterMasks_Initialize IN CCTK_POST_RECOVER_VARIABLES {
+SCHEDULE AsterMasks_SetMaskBoolZero IN CCTK_POST_RECOVER_VARIABLES BEFORE ODESolvers_PostStep AFTER AsterMasks_SetToOne {
 LANG: C
-WRITES: aster_mask_vc(everywhere)
-WRITES: aster_mask_cc(everywhere)
-} "Initialize VC and CC mask to true"
+OPTIONS: global
+WRITES:  MaskBool
+} "Set MaskBool = 0"

AsterMasks/src/Masks.cxx

@@ -0,0 +1,62 @@
+#include <loop.hxx>
+#include <loop_device.hxx>
+
+#include <cctk.h>
+#include <cctk_Arguments.h>
+#include <cctk_Parameters.h>
+
+extern "C" void AsterMasks_SetToOne(CCTK_ARGUMENTS) {
+  DECLARE_CCTK_ARGUMENTSX_AsterMasks_SetToOne;
+  DECLARE_CCTK_PARAMETERS;
+
+  grid.loop_all_device<0, 0, 0>(
+        grid.nghostzones,
+        [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
+
+          aster_mask_vc(p.I) = 1.0;
+
+  });
+
+  grid.loop_all_device<1, 1, 1>(
+        grid.nghostzones,
+        [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
+
+          aster_mask_cc(p.I) = 1.0;
+
+  });
+}
+
+extern "C" void AsterMasks_SetToOnePostStep(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTSX_AsterMasks_SetToOnePostStep;
+  DECLARE_CCTK_PARAMETERS;
+
+  if (*MaskBool == 1.0) {
+
+    grid.loop_all_device<0, 0, 0>(
+          grid.nghostzones,
+          [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
+
+            aster_mask_vc(p.I) = 1.0;
+
+    });
+
+    grid.loop_all_device<1, 1, 1>(
+          grid.nghostzones,
+          [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
+
+            aster_mask_cc(p.I) = 1.0;
+
+    });
+
+  }
+}
+
+extern "C" void AsterMasks_SetMaskBoolZero(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTSX_AsterMasks_SetMaskBoolZero;
+  DECLARE_CCTK_PARAMETERS;
+
+  *MaskBool = 0.0;
+
+}

AsterMasks/src/make.code.defn

@@ -2,7 +2,7 @@
 
 # Source files in this directory
 SRCS =	\
-  InitializeMasks.cxx
+  Masks.cxx
 
 # Subdirectories containing source files
 SUBDIRS = 

AsterSeeds/param.ccl

@@ -26,6 +26,8 @@ KEYWORD test_case "Name of the testcase" STEERABLE=never
   "cylindrical blast"        :: ""
   "magnetic rotor"           :: ""
   "magnetic loop advection"  :: ""
+  "KHI"                      :: ""
+  "magKHI"                   :: ""
 
   "spherical shock"          :: ""
   "magTOV"                   :: ""
@@ -75,7 +77,7 @@ CCTK_REAL press_atmosphere "floor pressure in the atmosphere"
 
 # parameters for 1D tests
 
-CCTK_REAL amplitude "Wave amplitude" STEERABLE=never
+CCTK_REAL amplitude "Amplitude" STEERABLE=never
 {
   0.0:* :: ""
 } 1.0e-3
@@ -111,6 +113,46 @@ CCTK_KEYWORD mag_loop_adv_axial_vel "Zero or non-zero fluid velocity along the d
   "non-zero" :: "Non-zero fluid velocity along the direction of the magnetic loop's axis"
 } "zero"
 
+# parameters for 2D KHI
+
+CCTK_REAL p_val "initial pressure" STEERABLE=always
+{ 
+  0.0:* :: ""
+} 2.5
+
+CCTK_REAL w0 "amplitude of pertubation applied to v^y" STEERABLE=always
+{ 
+  0.0:* :: ""
+} 0.1
+
+CCTK_REAL sigma "parameter for the exponential in v^y" STEERABLE=always
+{ 
+  0.0:* :: ""
+} 0.0353553390593274
+
+# Initial data for region above the slip surface
+
+CCTK_REAL rhoUp "initial density" STEERABLE=always
+{ 
+  0.0:* :: ""
+} 1.0
+
+CCTK_REAL vxUp "initial x component of velocity" STEERABLE=always
+{ 
+  *:* :: ""
+} -0.5
+
+# Initial data for region below the slip surface
+
+CCTK_REAL rhoLow "initial density" STEERABLE=always
+{ 
+  0.0:* :: ""
+} 2.0
+
+CCTK_REAL vxLow "initial x component of velocity" STEERABLE=always
+{ 
+  *:* :: ""
+} 0.5
 
 # ---------------------------------------
 

AsterSeeds/src/1D_tests.cxx

@@ -44,7 +44,7 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) {
           vely(p.I) = 0.0;
           velz(p.I) = 0.0;
           press(p.I) = 1.0;
-          eps(p.I) = eos_3p_ig->eps_from_valid_rho_press_ye(
+          eps(p.I) = eos_3p_ig->eps_from_rho_press_ye(
               rho(p.I), press(p.I), dummy_ye);
         });
 
@@ -72,7 +72,7 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) {
           vely(p.I) = 0.0;
           velz(p.I) = 0.0;
           press(p.I) = 1.0; // should add kinetic energy here
-          eps(p.I) = eos_3p_ig->eps_from_valid_rho_press_ye(
+          eps(p.I) = eos_3p_ig->eps_from_rho_press_ye(
               rho(p.I), press(p.I), dummy_ye);
         });
 
@@ -104,7 +104,7 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) {
           vely(p.I) = -va * A0 * cos(k * p.x);
           velz(p.I) = -va * A0 * sin(k * p.x);
           press(p.I) = 0.5; // should add kinetic energy here
-          eps(p.I) = eos_3p_ig->eps_from_valid_rho_press_ye(
+          eps(p.I) = eos_3p_ig->eps_from_rho_press_ye(
               rho(p.I), press(p.I), dummy_ye);
         });
 
@@ -143,7 +143,7 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) {
             velz(p.I) = 0.0;
             press(p.I) = 1.0;
           }
-          eps(p.I) = eos_3p_ig->eps_from_valid_rho_press_ye(
+          eps(p.I) = eos_3p_ig->eps_from_rho_press_ye(
               rho(p.I), press(p.I), dummy_ye);
         });
 
@@ -377,7 +377,7 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) {
             velz(p.I) = oldvrs(2);
             press(p.I) = pressr;
           }
-          eps(p.I) = eos_3p_ig->eps_from_valid_rho_press_ye(
+          eps(p.I) = eos_3p_ig->eps_from_rho_press_ye(
               rho(p.I), press(p.I), dummy_ye);
         });