fix flag, but segfault

IAS-Astrophysics · Sep 9, 2024 · 52fcbf8 · 52fcbf8
1 parent 686be14
commit 52fcbf8
Show file tree

Hide file tree

Showing 3 changed files with 17 additions and 9 deletions.
diff --git a/kokkos b/kokkos
diff --git a/src/z4c/z4c_adm.cpp b/src/z4c/z4c_adm.cpp
@@ -263,6 +263,10 @@ void Z4c::ADMConstraints(MeshBlockPack *pmbp) {
   auto &adm = pmbp->padm->adm;
   auto &tmunu = pmbp->ptmunu->tmunu;
   auto &u_con = pmbp->pz4c->u_con;
+
+  // vacuum or with matter?
+  bool is_vacuum = (pmbp->ptmunu == nullptr) ? true : false;
+
   Kokkos::deep_copy(u_con, 0.);
   auto &con = pmbp->pz4c->con;
   par_for("ADM constraints loop",DevExeSpace(),
@@ -466,14 +470,18 @@ void Z4c::ADMConstraints(MeshBlockPack *pmbp) {
     //
     // Hamiltonian constraint
     //
-    con.H(m,k,j,i) = R + SQR(K) - KK - 16*M_PI * tmunu.E(m,k,j,i);
-
+    con.H(m,k,j,i) = R + SQR(K) - KK;
+    if(!is_vacuum) {
+      con.H(m,k,j,i) -= 16*M_PI * tmunu.E(m,k,j,i);
+    }
     // Momentum constraint (contravariant)
     //
     for(int a = 0; a < 3; ++a) {
       M_u(a) = 0.0;
       for(int b = 0; b < 3; ++b) {
-        M_u(a) -= 8*M_PI * g_uu(a,b) * tmunu.S_d(m,b,k,j,i);
+        if(!is_vacuum) {
+          M_u(a) -= 8*M_PI * g_uu(a,b) * tmunu.S_d(m,b,k,j,i);
+        }
         for(int c = 0; c < 3; ++c) {
           M_u(a) += g_uu(a,b) * DK_udd(c,b,c);
           M_u(a) -= g_uu(b,c) * DK_udd(a,b,c);

diff --git a/src/z4c/z4c_calcrhs.cpp b/src/z4c/z4c_calcrhs.cpp
@@ -372,7 +372,7 @@ TaskStatus Z4c::CalcRHS(Driver *pdriver, int stage) {
     //    S(1) += oopsi4(1) * g_uu(a,b,i) * mat.S_dd(m,a,b,k,j,i);
     //  }
     //}
-    if(is_vacuum) {
+    if(!is_vacuum) {
       for (int a = 0; a < 3; ++a)
       for (int b = 0; b < 3; ++b) {
         S += oopsi4 * g_uu(a,b) * tmunu.S_dd(m,a,b,k,j,i);
@@ -499,15 +499,15 @@ TaskStatus Z4c::CalcRHS(Driver *pdriver, int stage) {
       LKhat + opt.damp_kappa1*(1 - opt.damp_kappa2)
       * z4c.alpha(m,k,j,i) * z4c.vTheta(m,k,j,i);
     // Matter term
-    if(is_vacuum) {
+    if(!is_vacuum) {
       rhs.vKhat(m,k,j,i) += 4.*M_PI * z4c.alpha(m,k,j,i) * (S + tmunu.E(m,k,j,i));
     }
     rhs.chi(m,k,j,i) = Lchi - (1./6.) * opt.chi_psi_power *
       chi_guarded * z4c.alpha(m,k,j,i) * K;
     rhs.vTheta(m,k,j,i) = LTheta + z4c.alpha(m,k,j,i) * (
         0.5*Ht - (2. + opt.damp_kappa2) * opt.damp_kappa1 * z4c.vTheta(m,k,j,i));
     // Matter term
-    if(is_vacuum) {
+    if(!is_vacuum) {
       rhs.vTheta(m,k,j,i) -= 8.*M_PI * z4c.alpha(m,k,j,i) * tmunu.E(m,k,j,i);
     }
     // If BSSN is enabled, theta is disabled.
@@ -520,7 +520,7 @@ TaskStatus Z4c::CalcRHS(Driver *pdriver, int stage) {
       for(int b = 0; b < 3; ++b) {
         rhs.vGam_u(m,a,k,j,i) -= 2. * A_uu(a,b) * dalpha_d(b);
         // Matter term
-        if(is_vacuum) {
+        if(!is_vacuum) {
           rhs.vGam_u(m,a,k,j,i) -= 16.*M_PI * z4c.alpha(m,k,j,i)
                               * g_uu(a,b) * tmunu.S_d(m,b,k,j,i);
         }
@@ -540,7 +540,7 @@ TaskStatus Z4c::CalcRHS(Driver *pdriver, int stage) {
                              - 2.*AA_dd(a,b));
       rhs.vA_dd(m,a,b,k,j,i) += LA_dd(a,b);
       // Matter term
-      if(is_vacuum) {
+      if(!is_vacuum) {
         rhs.vA_dd(m,a,b,k,j,i) -= 8.*M_PI * z4c.alpha(m,k,j,i) *
                 (oopsi4*tmunu.S_dd(m,a,b,k,j,i) - (1./3.)*S*z4c.g_dd(m,a,b,k,j,i));
       }