VolumeIntegrals_GRMHDX Updates

MHD_Diagnostics/param.ccl

@@ -13,6 +13,9 @@
 # USES CCTK_INT metric_timelevels
 #############################################################################
 
+shares: AsterX
+USES CCTK_INT mag_correction_order
+
 private:
 
 CCTK_INT diagnostics_compute_every "How often to compute diagnostics?" STEERABLE=ALWAYS

MHD_Diagnostics/src/master.cxx

@@ -387,29 +387,11 @@ extern "C" void compute_mhd_derivs(CCTK_ARGUMENTS) {
     double invdet = 1.0 / det;
 
     // Calculate also divB from AsterX staggered field
-    // Use fourth-order ECHO expressions
-    // TODO: check if face index is correct, how to handle densitized B?
-
-    double Bx_stagg = 13. * dBx_stag(p.I + p.DI[0]) / 12. -
-                      (dBx_stag(p.I) + dBx_stag(p.I + 2 * p.DI[0])) / 24.;
-    double By_stagg = 13. * dBy_stag(p.I + p.DI[1]) / 12. -
-                      (dBy_stag(p.I) + dBy_stag(p.I + 2 * p.DI[1])) / 24.;
-    double Bz_stagg = 13. * dBz_stag(p.I + p.DI[2]) / 12. -
-                      (dBz_stag(p.I) + dBz_stag(p.I + 2 * p.DI[2])) / 24.;
-
-    double Bx_stagg_m1 =
-        13. * dBx_stag(p.I) / 12. -
-        (dBx_stag(p.I + p.DI[0]) + dBx_stag(p.I - p.DI[0])) / 24.;
-    double By_stagg_m1 =
-        13. * dBy_stag(p.I) / 12. -
-        (dBy_stag(p.I + p.DI[1]) + dBy_stag(p.I - p.DI[1])) / 24.;
-    double Bz_stagg_m1 =
-        13. * dBz_stag(p.I) / 12. -
-        (dBz_stag(p.I + p.DI[2]) + dBz_stag(p.I - p.DI[2])) / 24.;
-
-    divB(p.I) = dxi[1] * (Bx_stagg - Bx_stagg_m1) +
-                dxi[2] * (By_stagg - By_stagg_m1) +
-                dxi[3] * (Bz_stagg - Bz_stagg_m1);
+    // using fourth-order ECHO expressions
+
+    divB(p.I) = calc_fd_forward_midpoint<0>(dBx_stag, p, mag_correction_order) +
+      calc_fd_forward_midpoint<1>(dBy_stag, p, mag_correction_order) +
+      calc_fd_forward_midpoint<2>(dBz_stag, p, mag_correction_order);
 
     // ===========================================
 

MeudonBNSX/src/Bin_NS.cxx

@@ -28,7 +28,7 @@ using namespace Arith;
 using namespace Loop;
 
 // define namespace here for old versions of Lorene that don't do so
-namespace Lorene {}
+// namespace Lorene {}
 using namespace Lorene;
 
 namespace {
@@ -75,9 +75,6 @@ extern "C" void MeudonBNSX_initialise(CCTK_ARGUMENTS) {
   int const npoints = nx * ny * nz;
   vector<double> xx(npoints), yy(npoints), zz(npoints);
 
-  // TODO: currently works only with polytropic EOS
-  CCTK_INFO("MeudonBNSX will use the polytropic equation of state.");
-
   grid.loop_all<1, 1, 1>(grid.nghostzones,
                          [&](const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
                            CCTK_INT idx = CCC_layout.linear(p.I);
@@ -210,7 +207,9 @@ extern "C" void MeudonBNSX_initialise(CCTK_ARGUMENTS) {
           rho(p.I) = bin_ns.nbar[idx] / rho_unit;
           if (!recalculate_eps)
             eps(p.I) = bin_ns.ener_spec[idx];
-          // TODO: currently works only with polytropic EOS
+          // TODO: this only works for simple polytropic, or tabulated EOS.
+          // The expressions are valid for polytropes. Eps and press can and must be
+          // recalculated in a separate thorn for tabulated EOS.
           else
             eps(p.I) = K * pow(rho(p.I), bin_ns.gamma_poly1 - 1.) /
                        (bin_ns.gamma_poly1 - 1.);

VolumeIntegrals_GRMHDX/interface.ccl

@@ -6,6 +6,22 @@ inherits: ADMBaseX HydroBaseX AsterX BoxInBox
 
 USES INCLUDE HEADER: loop_device.hxx
 
+public:
+CCTK_REAL GlobalCoM type=SCALAR tags='checkpoint="no"'
+{
+  comx, comy, comz
+} "Record CoM to set origin for toroidal/poloidal decompositions"
+
+CCTK_REAL CoM1 type=SCALAR tags='checkpoint="no"'
+{
+  comx1, comy1, comz1
+} "Record CoM of NS 1"
+
+CCTK_REAL CoM2 type=SCALAR tags='checkpoint="no"'
+{
+  comx2, comy2, comz2
+} "Record CoM of NS 2"
+
 private:
 CCTK_REAL VolIntegrands TYPE=GF CENTERING={ccc} tags='checkpoint="no"'
 {
@@ -33,7 +49,3 @@ CCTK_REAL VolIntegrals_vacuum_time type=SCALAR tags='checkpoint="no"'
   physical_time
 } "keeps track of the physical time, in case time coordinate is reparameterized, a la http://arxiv.org/abs/1404.6523"
 
-CCTK_REAL GlobalCoM type=SCALAR tags='checkpoint="no"'
-{
-  comx, comy, comz
-} "Record CoM to set origin for toroidal/poloidal decompositions"

VolumeIntegrals_GRMHDX/param.ccl

@@ -86,8 +86,16 @@ keyword Integration_quantity_keyword[101] "Which quantity to integrate" STEERABL
 "thermal_energy"                :: "Thermal energy"
 "volume_average_norm_B_ab"      :: "Volume average of norm B in density interval (dens_a,dens_b]"       
 "volume_average_norm_B_cd"      :: "Volume average of norm B in density interval (dens_c,dens_d]"
+"volume_average_norm_B_posz"      :: "Volume average of norm B where z > CoM z"
+"volume_average_norm_B_negz"      :: "Volume average of norm B where z < CoM z"
+"volume_average_norm_B_m1_cd"      :: "Volume average of m=1 mode of tor., pol. B in (dens_c,dens_d]"
+"volume_average_norm_B_m2_cd"      :: "Volume average of m=2 mode of tor., pol. B in (dens_c,dens_d]"
+"volume_average_norm_B_m3_cd"      :: "Volume average of m=3 mode of tor., pol. B in (dens_c,dens_d]"
+"volume_average_norm_B_m4_cd"      :: "Volume average of m=4 mode of tor., pol. B in (dens_c,dens_d]"
 "em_energy_ab"                  :: "Total energy in electro-magnetic sector in density interval (dens_a,dens_b]"       
 "em_energy_cd"                  :: "Total energy in electro-magnetic sector in density interval (dens_c,dens_d]"
+"em_energy_posz"                  :: "Total energy in electro-magnetic sector where z > CoM z"
+"em_energy_negz"                  :: "Total energy in electro-magnetic sector where z < CoM z"
 "coordvolume"                  		:: "Coordinate volume of region with rho > dens_atmo"
 } "nothing"
 

VolumeIntegrals_GRMHDX/schedule.ccl

@@ -21,11 +21,11 @@ SCHEDULE VI_GRMHDX_file_output_routine_Startup AT CCTK_POST_RECOVER_VARIABLES
 }
 
 #############################
-SCHEDULE VI_GRMHDX_InitializeIntegralCounterToZero AT CCTK_INITIAL
+SCHEDULE VI_GRMHDX_InitializeIntegralCounterToZero AT CCTK_INITIAL BEFORE AsterX_InitialGroup
 {
   LANG: C
   OPTIONS: GLOBAL
-  WRITES: IntegralCounter(everywhere) GlobalCoM(everywhere) 
+  WRITES: IntegralCounter(everywhere) GlobalCoM(everywhere) CoM1(everywhere) CoM2(everywhere) 
   WRITES: VolIntegral(everywhere)
   WRITES: volintegral_inside_sphere__center_x(everywhere)
   WRITES: volintegral_inside_sphere__center_y(everywhere)
@@ -39,7 +39,7 @@ SCHEDULE VI_GRMHDX_InitializeIntegralCounterToZero AT CCTK_POST_RECOVER_VARIABLE
 {
   LANG: C
   OPTIONS: GLOBAL
-  WRITES: IntegralCounter(everywhere) GlobalCoM(everywhere)
+  WRITES: IntegralCounter(everywhere) GlobalCoM(everywhere) CoM1(everywhere) CoM2(everywhere) 
   WRITES: VolIntegral(everywhere)
   WRITES: volintegral_inside_sphere__center_x(everywhere)
   WRITES: volintegral_inside_sphere__center_y(everywhere)
@@ -57,7 +57,7 @@ SCHEDULE VI_GRMHDX_InitializeIntegralCounter before VI_GRMHDX_VolumeIntegralGrou
   WRITES: IntegralCounter(everywhere)
 } "Initialize IntegralCounter variable"
 ##################
-SCHEDULE GROUP VI_GRMHDX_VolumeIntegralGroup AT poststep BEFORE EstimateError WHILE VolumeIntegrals_GRMHDX::IntegralCounter
+SCHEDULE GROUP VI_GRMHDX_VolumeIntegralGroup AT poststep BEFORE (EstimateError Seed_During_Evolution) WHILE VolumeIntegrals_GRMHDX::IntegralCounter
 {
 } "Evaluate all volume integrals"
 
@@ -105,7 +105,7 @@ SCHEDULE VI_GRMHDX_DoSum in VI_GRMHDX_VolumeIntegralGroup AFTER VI_GRMHDX_Comput
   READS: BoxInBox::position_x(everywhere)
   READS: BoxInBox::position_y(everywhere)
   READS: BoxInBox::position_z(everywhere)
-  WRITES: VolIntegral(everywhere) GlobalCoM(everywhere) 
+  WRITES: VolIntegral(everywhere) GlobalCoM(everywhere) CoM1(everywhere) CoM2(everywhere) 
   WRITES: BoxInBox::position_x(everywhere)
   WRITES: BoxInBox::position_y(everywhere)
   WRITES: BoxInBox::position_z(everywhere)

VolumeIntegrals_GRMHDX/src/Initialize_or_Decrement_IntegralCounter.cxx

@@ -17,6 +17,12 @@ extern "C" void VI_GRMHDX_InitializeIntegralCounterToZero(CCTK_ARGUMENTS) {
   *comx = 0;
   *comy = 0;
   *comz = 0;
+  *comx1 = 0;
+  *comy1 = 0;
+  *comz1 = 0;
+  *comx2 = 0;
+  *comy2 = 0;
+  *comz2 = 0;
 
   // Init Array that Holds Results
   for (int ii = 0; ii < 101; ii++) {

VolumeIntegrals_GRMHDX/src/evaluate_integrands_local.cxx

@@ -213,6 +213,50 @@ extern "C" void VI_GRMHDX_ComputeIntegrand(CCTK_ARGUMENTS) {
                           velx, vely, velz, Bvecx, Bvecy, Bvecz, gxx, gxy, gxz,
                           gyy, gyz, gzz, cms_x, cms_y);
         });
+  } else if (CCTK_EQUALS(Integration_quantity_keyword[which_integral],
+                         "em_energy_posz")) {
+    grid.loop_all_device<1, 1, 1>(
+        grid.nghostzones,
+        [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
+          // TODO: Get BNS COM from --BNSTrackerGen-- somewhere
+          // Make sure we have updated values
+          // before this integral is called!
+
+          double cms_x = 0.0;
+          double cms_y = 0.0;
+          double cms_z = 0.0;
+          if (set_origin_with_VIX) {
+            cms_x = *comx;
+            cms_y = *comy;
+            cms_z = *comz;
+          }
+
+          magnetic_tot_zsplit(VolIntegrand1, VolIntegrand2, p,
+                          velx, vely, velz, Bvecx, Bvecy, Bvecz, gxx, gxy, gxz,
+                          gyy, gyz, gzz, cms_x, cms_y, cms_z, 1);
+        });
+  } else if (CCTK_EQUALS(Integration_quantity_keyword[which_integral],
+                         "em_energy_negz")) {
+    grid.loop_all_device<1, 1, 1>(
+        grid.nghostzones,
+        [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
+          // TODO: Get BNS COM from --BNSTrackerGen-- somewhere
+          // Make sure we have updated values
+          // before this integral is called!
+
+          double cms_x = 0.0;
+          double cms_y = 0.0;
+          double cms_z = 0.0;
+          if (set_origin_with_VIX) {
+            cms_x = *comx;
+            cms_y = *comy;
+            cms_z = *comz;
+          }
+
+          magnetic_tot_zsplit(VolIntegrand1, VolIntegrand2, p, 
+                          velx, vely, velz, Bvecx, Bvecy, Bvecz, gxx, gxy, gxz,
+                          gyy, gyz, gzz, cms_x, cms_y, cms_z, -1);
+        });
   } else if (CCTK_EQUALS(Integration_quantity_keyword[which_integral],
                          "volume_average_norm_B_ab")) {
     grid.loop_all_device<1, 1, 1>(
@@ -254,12 +298,128 @@ extern "C" void VI_GRMHDX_ComputeIntegrand(CCTK_ARGUMENTS) {
                            Bvecz, gxx, gxy, gxz, gyy, gyz, gzz, cms_x, cms_y);
         });
   } else if (CCTK_EQUALS(Integration_quantity_keyword[which_integral],
-                         "magnetic_energy_comov")) {
+                         "volume_average_norm_B_posz")) {
     grid.loop_all_device<1, 1, 1>(
         grid.nghostzones,
         [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
-          magnetic_co(VolIntegrand1, p, b2small, w_lorentz, gxx, gxy, gxz, gyy,
-                      gyz, gzz, gamma_lim);
+          // TODO: Get BNS COM from --BNSTrackerGen-- somewhere
+          // Make sure we have updated values
+          // before this integral is called!
+
+          double cms_x = 0.0;
+          double cms_y = 0.0;
+          double cms_z = 0.0;
+          if (set_origin_with_VIX) {
+            cms_x = *comx;
+            cms_y = *comy;
+            cms_z = *comz;
+          }
+
+          volume_norm_B_zsplit(VolIntegrand1, VolIntegrand2, VolIntegrand3,
+                           VolIntegrand4, p, Bvecx, Bvecy,
+                           Bvecz, gxx, gxy, gxz, gyy, gyz, gzz, cms_x, cms_y, cms_z, 1);
+        });
+  } else if (CCTK_EQUALS(Integration_quantity_keyword[which_integral],
+                         "volume_average_norm_B_negz")) {
+    grid.loop_all_device<1, 1, 1>(
+        grid.nghostzones,
+        [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
+          // TODO: Get BNS COM from --BNSTrackerGen-- somewhere
+          // Make sure we have updated values
+          // before this integral is called!
+
+          double cms_x = 0.0;
+          double cms_y = 0.0;
+          double cms_z = 0.0;
+          if (set_origin_with_VIX) {
+            cms_x = *comx;
+            cms_y = *comy;
+            cms_z = *comz;
+          }
+
+          volume_norm_B_zsplit(VolIntegrand1, VolIntegrand2, VolIntegrand3,
+                           VolIntegrand4, p, Bvecx, Bvecy,
+                           Bvecz, gxx, gxy, gxz, gyy, gyz, gzz, cms_x, cms_y, cms_z, -1);
+        });
+  } else if (CCTK_EQUALS(Integration_quantity_keyword[which_integral],
+                         "volume_average_norm_B_m1_cd")) {
+    grid.loop_all_device<1, 1, 1>(
+        grid.nghostzones,
+        [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
+          // TODO: Get BNS COM from --BNSTrackerGen-- somewhere
+          // Make sure we have updated values
+          // before this integral is called!
+
+          double cms_x = 0.0;
+          double cms_y = 0.0;
+          if (set_origin_with_VIX) {
+            cms_x = *comx;
+            cms_y = *comy;
+          }
+
+          volume_norm_B_m_12(VolIntegrand1, VolIntegrand2, VolIntegrand3,
+                           VolIntegrand4, p, rho, dens_c, dens_d, Bvecx, Bvecy,
+                           Bvecz, gxx, gxy, gxz, gyy, gyz, gzz, cms_x, cms_y, 1);
+        });
+  } else if (CCTK_EQUALS(Integration_quantity_keyword[which_integral],
+                         "volume_average_norm_B_m2_cd")) {
+    grid.loop_all_device<1, 1, 1>(
+        grid.nghostzones,
+        [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
+          // TODO: Get BNS COM from --BNSTrackerGen-- somewhere
+          // Make sure we have updated values
+          // before this integral is called!
+
+          double cms_x = 0.0;
+          double cms_y = 0.0;
+          if (set_origin_with_VIX) {
+            cms_x = *comx;
+            cms_y = *comy;
+          }
+
+          volume_norm_B_m_12(VolIntegrand1, VolIntegrand2, VolIntegrand3,
+                           VolIntegrand4, p, rho, dens_c, dens_d, Bvecx, Bvecy,
+                           Bvecz, gxx, gxy, gxz, gyy, gyz, gzz, cms_x, cms_y, 2);
+        });
+  } else if (CCTK_EQUALS(Integration_quantity_keyword[which_integral],
+                         "volume_average_norm_B_m3_cd")) {
+    grid.loop_all_device<1, 1, 1>(
+        grid.nghostzones,
+        [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
+          // TODO: Get BNS COM from --BNSTrackerGen-- somewhere
+          // Make sure we have updated values
+          // before this integral is called!
+
+          double cms_x = 0.0;
+          double cms_y = 0.0;
+          if (set_origin_with_VIX) {
+            cms_x = *comx;
+            cms_y = *comy;
+          }
+
+          volume_norm_B_m_12(VolIntegrand1, VolIntegrand2, VolIntegrand3,
+                           VolIntegrand4, p, rho, dens_c, dens_d, Bvecx, Bvecy,
+                           Bvecz, gxx, gxy, gxz, gyy, gyz, gzz, cms_x, cms_y, 3);
+        });
+  } else if (CCTK_EQUALS(Integration_quantity_keyword[which_integral],
+                         "volume_average_norm_B_m4_cd")) {
+    grid.loop_all_device<1, 1, 1>(
+        grid.nghostzones,
+        [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
+          // TODO: Get BNS COM from --BNSTrackerGen-- somewhere
+          // Make sure we have updated values
+          // before this integral is called!
+
+          double cms_x = 0.0;
+          double cms_y = 0.0;
+          if (set_origin_with_VIX) {
+            cms_x = *comx;
+            cms_y = *comy;
+          }
+
+          volume_norm_B_m_12(VolIntegrand1, VolIntegrand2, VolIntegrand3,
+                           VolIntegrand4, p, rho, dens_c, dens_d, Bvecx, Bvecy,
+                           Bvecz, gxx, gxy, gxz, gyy, gyz, gzz, cms_x, cms_y, 4);
         });
   } else if (CCTK_EQUALS(Integration_quantity_keyword[which_integral],
                          "magnetic_energy_comov")) {