horizon find and cartesian grid dump

src/CMakeLists.txt

@@ -147,8 +147,10 @@ add_executable(
         utils/show_config.cpp
         utils/lagrange_interpolator.cpp
         utils/tr_table.cpp
+        utils/cart_grid.cpp
 
         z4c/compact_object_tracker.cpp
+        z4c/horizon_dump.cpp
         z4c/tmunu.cpp
         z4c/z4c.cpp
         z4c/z4c_adm.cpp

src/athena.hpp

@@ -130,6 +130,10 @@ template <typename T>
 using DualArray2D = Kokkos::DualView<T **, LayoutWrapper, DevMemSpace>;
 template <typename T>
 using DualArray3D = Kokkos::DualView<T ***, LayoutWrapper, DevMemSpace>;
+template <typename T>
+using DualArray4D = Kokkos::DualView<T ****, LayoutWrapper, DevMemSpace>;
+template <typename T>
+using DualArray5D = Kokkos::DualView<T *****, LayoutWrapper, DevMemSpace>;
 
 // template declarations for construction of Kokkos::View in scratch memory
 template <typename T>
@@ -477,4 +481,3 @@ struct reduction_identity< array_sum::GlobalSum > {
 }
 
 #endif // ATHENA_HPP_
-

src/driver/driver.cpp

@@ -258,7 +258,7 @@ Driver::Driver(ParameterInput *pin, Mesh *pmesh, Real wtlim, Kokkos::Timer* ptim
     } else {
       std::cout << "### FATAL ERROR in " << __FILE__ << " at line " << __LINE__
          << std::endl << "integrator=" << integrator << " not implemented. "
-         << "Valid choices are [rk1,rk2,rk3,imex2,imex3]." << std::endl;
+         << "Valid choices are [rk1,rk2,rk3,rk4,imex2,imex3]." << std::endl;
       exit(EXIT_FAILURE);
     }
   }

src/outputs/restart.cpp

@@ -225,7 +225,7 @@ void RestartOutput::WriteOutputFile(Mesh *pm, ParameterInput *pin) {
     // output puncture tracker data
     if (nco > 0) {
       for (auto & pt : pz4c->ptracker) {
-        resfile.Write_any_type(pt.GetPos(), 3*sizeof(Real), "byte");
+        resfile.Write_any_type(pt->GetPos(), 3*sizeof(Real), "byte");
       }
     }
     // turbulence driver internal RNG

src/pgen/gr_analytic/kerr_schild.cpp

@@ -0,0 +1,125 @@
+//========================================================================================
+// AthenaXXX astrophysical plasma code
+// Copyright(C) 2020 James M. Stone <[email protected]> and the Athena code team
+// Licensed under the 3-clause BSD License (the "LICENSE")
+//========================================================================================
+//! \file z4c_one_puncture.cpp
+//  \brief Problem generator for a single puncture placed at the origin of the domain
+
+#include <algorithm>
+#include <cmath>
+#include <sstream>
+#include <iomanip>
+#include <iostream>   // endl
+#include <limits>     // numeric_limits::max()
+#include <memory>
+#include <string>     // c_str(), string
+#include <vector>
+
+#include "athena.hpp"
+#include "parameter_input.hpp"
+#include "globals.hpp"
+#include "mesh/mesh.hpp"
+#include "z4c/z4c.hpp"
+#include "z4c/z4c_amr.hpp"
+#include "coordinates/adm.hpp"
+#include "coordinates/cell_locations.hpp"
+#include "coordinates/cartesian_ks.hpp"
+
+
+void KerrSchild(MeshBlockPack *pmbp, ParameterInput *pin);
+void RefinementCondition(MeshBlockPack* pmbp);
+KOKKOS_INLINE_FUNCTION
+bool invertSymmetricMatrix4x4(const double g_dd[4][4], double g_dd_inv[4][4]);
+
+//----------------------------------------------------------------------------------------
+//! \fn ProblemGenerator::UserProblem_()
+//! \brief Problem Generator for single puncture
+void ProblemGenerator::UserProblem(ParameterInput *pin, const bool restart) {
+  user_ref_func  = RefinementCondition;
+  MeshBlockPack *pmbp = pmy_mesh_->pmb_pack;
+  auto &indcs = pmy_mesh_->mb_indcs;
+
+  if (pmbp->pz4c == nullptr) {
+    std::cout << "### FATAL ERROR in " << __FILE__ << " at line " << __LINE__ << std::endl
+              << "One Puncture test can only be run in Z4c, but no <z4c> block "
+              << "in input file" << std::endl;
+    exit(EXIT_FAILURE);
+  }
+
+  KerrSchild(pmbp, pin);
+  switch (indcs.ng) {
+    case 2: pmbp->pz4c->ADMToZ4c<2>(pmbp, pin);
+            break;
+    case 3: pmbp->pz4c->ADMToZ4c<3>(pmbp, pin);
+            break;
+    case 4: pmbp->pz4c->ADMToZ4c<4>(pmbp, pin);
+            break;
+  }
+  pmbp->pz4c->Z4cToADM(pmbp);
+  pmbp->pz4c->GaugePreCollapsedLapse(pmbp, pin);
+  switch (indcs.ng) {
+    case 2: pmbp->pz4c->ADMConstraints<2>(pmbp);
+            break;
+    case 3: pmbp->pz4c->ADMConstraints<3>(pmbp);
+            break;
+    case 4: pmbp->pz4c->ADMConstraints<4>(pmbp);
+            break;
+  }
+  std::cout<<"Kerr Schild initialized."<<std::endl;
+  return;
+}
+
+//----------------------------------------------------------------------------------------
+//! \fn void ADMOnePuncture(MeshBlockPack *pmbp, ParameterInput *pin)
+//! \brief Initialize ADM vars to single puncture (no spin)
+
+void KerrSchild(MeshBlockPack *pmbp, ParameterInput *pin) {
+  Real ADM_mass = pin->GetOrAddReal("problem", "punc_ADM_mass", 1.);
+  Real center_x1 = pin->GetOrAddReal("problem", "punc_center_x1", 0.);
+  Real center_x2 = pin->GetOrAddReal("problem", "punc_center_x2", 0.);
+  Real center_x3 = pin->GetOrAddReal("problem", "punc_center_x3", 0.);
+  Real a = 0;
+  bool minkowski = false;
+
+  // capture variables for the kernel
+  auto &indcs = pmbp->pmesh->mb_indcs;
+  auto &size = pmbp->pmb->mb_size;
+  int &is = indcs.is; int &ie = indcs.ie;
+  int &js = indcs.js; int &je = indcs.je;
+  int &ks = indcs.ks; int &ke = indcs.ke;
+  int nmb = pmbp->nmb_thispack;
+  auto &adm = pmbp->padm->adm;
+  int &ng = indcs.ng;
+  int n1 = indcs.nx1 + 2*ng;
+  int n2 = indcs.nx2 + 2*ng;
+  int n3 = indcs.nx3 + 2*ng;
+  par_for("pgen_adm_vars", DevExeSpace(), 0,nmb-1,0,(n3-1),0,(n2-1),0,(n1-1),
+  KOKKOS_LAMBDA(int m, int k, int j, int i) {
+    Real &x1min = size.d_view(m).x1min;
+    Real &x1max = size.d_view(m).x1max;
+    Real x1v = CellCenterX(i-is, indcs.nx1, x1min, x1max);
+
+    Real &x2min = size.d_view(m).x2min;
+    Real &x2max = size.d_view(m).x2max;
+    Real x2v = CellCenterX(j-js, indcs.nx2, x2min, x2max);
+
+    Real &x3min = size.d_view(m).x3min;
+    Real &x3max = size.d_view(m).x3max;
+    Real x3v = CellCenterX(k-ks, indcs.nx3, x3min, x3max);
+
+    ComputeADMDecomposition(x1v, x2v, x3v, minkowski, a,
+      &adm.alpha(m,k,j,i),
+      &adm.beta_u(m,0,k,j,i), &adm.beta_u(m,1,k,j,i), &adm.beta_u(m,2,k,j,i),
+      &adm.psi4(m,k,j,i),
+      &adm.g_dd(m,0,0,k,j,i), &adm.g_dd(m,0,1,k,j,i), &adm.g_dd(m,0,2,k,j,i),
+      &adm.g_dd(m,1,1,k,j,i), &adm.g_dd(m,1,2,k,j,i), &adm.g_dd(m,2,2,k,j,i),
+      &adm.vK_dd(m,0,0,k,j,i), &adm.vK_dd(m,0,1,k,j,i), &adm.vK_dd(m,0,2,k,j,i),
+      &adm.vK_dd(m,1,1,k,j,i), &adm.vK_dd(m,1,2,k,j,i), &adm.vK_dd(m,2,2,k,j,i));
+  });
+}
+
+// how decide the refinement
+void RefinementCondition(MeshBlockPack* pmbp) {
+  pmbp->pz4c->pamr->Refine(pmbp);
+}

src/pgen/pgen.cpp

@@ -202,7 +202,7 @@ ProblemGenerator::ProblemGenerator(ParameterInput *pin, Mesh *pm, IOWrapper resf
 #if MPI_PARALLEL_ENABLED
       MPI_Bcast(&pos[0], 3*sizeof(Real), MPI_CHAR, 0, MPI_COMM_WORLD);
 #endif
-      pt.SetPos(&pos[0]);
+      pt->SetPos(&pos[0]);
     }
   }
 

src/tasklist/numerical_relativity.hpp

@@ -82,6 +82,7 @@ enum TaskName {
   Z4c_ClearRW,
   Z4c_Wave,
   Z4c_PT,
+  Z4c_DumpHorizon,
   Z4c_NTASKS
 };