add two apple-with-apple tests (#625)

* fixing ricci scalar

* spin for kerr-schild

* strong stability test

* input script for stability test

* fix bug and cleanup

* delete unused variables

* gauge wave test

* add input for gauge wave, test still failing

* add initial lapse, still does not fix the error

* update input files

* moving the awa tests to the pgen folder

* fix default build

---------

Co-authored-by: Hengrui Zhu <[email protected]>

inputs/z4c/awa/z4c_gauge_wave.athinput

@@ -0,0 +1,76 @@
+# AthenaXXX input file for Z4c linear wave tests
+
+<comment>
+problem   = z4c linear waves
+reference = e.g. Daverio et al. arxiv:1810.12346 (2018)
+
+<job>
+basename  = gauge    # problem ID: basename of output filenames
+
+<mesh>
+nghost    = 2          # Number of ghost cells
+nx1       = 50         # Number of zones in X1-direction
+x1min     = 0.0        # minimum value of X1
+x1max     = 1.0        # maximum value of X1
+ix1_bc    = periodic   # inner-X1 boundary flag
+ox1_bc    = periodic   # outer-X1 boundary flag
+
+nx2       = 4         # Number of zones in X2-direction
+x2min     = 0.0        # minimum value of X2
+x2max     = 1.0        # maximum value of X2
+ix2_bc    = periodic   # inner-X2 boundary flag
+ox2_bc    = periodic   # outer-X2 boundary flag
+
+nx3       = 4         # Number of zones in X3-direction
+x3min     = 0.0        # minimum value of X3
+x3max     = 1.0        # maximum value of X3
+ix3_bc    = periodic   # inner-X3 boundary flag
+ox3_bc    = periodic   # outer-X3 boundary flag
+
+<meshblock>
+nx1       = 50         # Number of cells in each MeshBlock, X1-dir
+nx2       = 4         # Number of cells in each MeshBlock, X2-dir
+nx3       = 4         # Number of cells in each MeshBlock, X3-dir
+
+<time>
+evolution  = dynamic   # dynamic/kinematic/static
+integrator = rk4       # time integration algorithm
+cfl_number = 0.5       # The Courant, Friedrichs, & Lewy (CFL) Number
+nlim       = -1        # cycle limit (no limit if <0)
+tlim       = 1000.0       # time limit
+ndiag      = 1         # cycles between diagostic output
+
+<z4c>
+lapse_harmonic  = 0.0       # Harmonic lapse parameter mu_L
+lapse_oplog     = 2.0       # 1+log lapse parameter
+shift_eta       = 2.0       # Shift damping term
+diss            = 0.02       # Kreiss-Oliger dissipation parameter
+chi_div_floor   = 0.00001   # Floor on the conformal factor
+damp_kappa1     = 0.02       # Constraint damping factor 1
+damp_kappa2     = 0.0       # Constraint damping factor 2
+
+<problem>
+pgen_name = z4c_gauge_wave # problem generator name
+amp = 0.01
+
+<output1>
+file_type   = tab       # legacy VTK output
+variable    = z4c   # variables to be output
+dt          = 1      # time increment between outputs
+slice_x2    = 0.5
+slice_x3    = 0.5
+ghost_zones = true     # switch to output ghost cells
+
+<output2>
+file_type   = hst       # history data dump
+user_hist   = false
+data_format = %12.5e    # Optional data format string
+dt          = 1       # time increment between outputs
+
+#<output3>
+#file_type   = tab       # legacy VTK output
+#variable    = z4c   # variables to be output
+#dt          = 0.025      # time increment between outputs
+#slice_x2    = 0.5
+#slice_x3    = 0.5
+#ghost_zones = true     # switch to output ghost cells

inputs/z4c/awa/z4c_stability.athinput

@@ -0,0 +1,76 @@
+# AthenaXXX input file for Z4c linear wave tests
+
+<comment>
+problem   = z4c linear waves
+reference = e.g. Daverio et al. arxiv:1810.12346 (2018)
+
+<job>
+basename  = stability    # problem ID: basename of output filenames
+
+<mesh>
+nghost    = 2          # Number of ghost cells
+nx1       = 4         # Number of zones in X1-direction
+x1min     = 0.0        # minimum value of X1
+x1max     = 1.0        # maximum value of X1
+ix1_bc    = periodic   # inner-X1 boundary flag
+ox1_bc    = periodic   # outer-X1 boundary flag
+
+nx2       = 4         # Number of zones in X2-direction
+x2min     = 0.0        # minimum value of X2
+x2max     = 1.0        # maximum value of X2
+ix2_bc    = periodic   # inner-X2 boundary flag
+ox2_bc    = periodic   # outer-X2 boundary flag
+
+nx3       = 50         # Number of zones in X3-direction
+x3min     = 0.0        # minimum value of X3
+x3max     = 1.0        # maximum value of X3
+ix3_bc    = periodic   # inner-X3 boundary flag
+ox3_bc    = periodic   # outer-X3 boundary flag
+
+<meshblock>
+nx1       = 4         # Number of cells in each MeshBlock, X1-dir
+nx2       = 4         # Number of cells in each MeshBlock, X2-dir
+nx3       = 50         # Number of cells in each MeshBlock, X3-dir
+
+<time>
+evolution  = dynamic   # dynamic/kinematic/static
+integrator = rk4       # time integration algorithm
+cfl_number = 0.5       # The Courant, Friedrichs, & Lewy (CFL) Number
+nlim       = -1        # cycle limit (no limit if <0)
+tlim       = 1000.0       # time limit
+ndiag      = 1         # cycles between diagostic output
+
+<z4c>
+diss       = 0.02
+lapse_harmonic  = 0.0                 # Harmonic lapse parameter mu_L
+lapse_oplog     = 1.0                 # 1+log lapse parameter
+shift_eta       = 2.0                 # Shift damping term
+chi_div_floor   = 0.00001
+damp_kappa1     = 0.00      # Constraint damping factor 1
+damp_kappa2     = 0.02
+
+
+<problem>
+pgen_name = z4c_stability # problem generator name
+user_hist = true      # enroll user-defined history function
+rho = 1
+
+<output1>
+file_type   = tab       # legacy VTK output
+variable    = z4c   # variables to be output
+dt          = 100      # time increment between outputs
+slice_x1    = 0.5
+slice_x2    = 0.5
+ghost_zones = true     # switch to output ghost cells
+
+<output2>
+file_type   = hst       # history data dump
+user_hist   = true
+data_format = %12.5e    # Optional data format string
+dt          = 1       # time increment between outputs
+
+<output2>
+file_type   = hst       # history data dump
+user_hist   = false
+data_format = %12.5e    # Optional data format string
+dt          = 1       # time increment between outputs

src/outputs/history.cpp

@@ -161,14 +161,14 @@ void HistoryOutput::LoadHydroHistoryData(HistoryData *pdata, Mesh *pm) {
 void HistoryOutput::LoadZ4cHistoryData(HistoryData *pdata, Mesh *pm) {
   // set number of and names of history variables for z4c
   pdata->nhist = 8;
-  pdata->label[0] = "H-norm2";
-  pdata->label[1] = "M-norm2";
-  pdata->label[2] = "Mx-norm2";
-  pdata->label[3] = "My-norm2";
-  pdata->label[4] = "Mz-norm2";
-  pdata->label[5] = "Z-norm2";
-  pdata->label[6] = "Theta-norm2";
-  pdata->label[7] = "C-norm2";
+  pdata->label[0] = "C-norm2";
+  pdata->label[1] = "H-norm2";
+  pdata->label[2] = "M-norm2";
+  pdata->label[3] = "Z-norm2";
+  pdata->label[4] = "Mx-norm2";
+  pdata->label[5] = "My-norm2";
+  pdata->label[6] = "Mz-norm2";
+  pdata->label[7] = "Theta-norm2";
 
   // capture class variabels for kernel
   auto &u0_ = pm->pmb_pack->pz4c->u0;
@@ -201,14 +201,14 @@ void HistoryOutput::LoadZ4cHistoryData(HistoryData *pdata, Mesh *pm) {
 
     // Hydro conserved variables:
     array_sum::GlobalSum hvars;
-    hvars.the_array[0] = vol*SQR(u_con_(m,0,k,j,i)); // ||H||^2
-    hvars.the_array[1] = vol*u_con_(m,1,k,j,i);      // ||M||^2 (comes already squared)
-    hvars.the_array[2] = vol*SQR(u_con_(m,2,k,j,i)); // ||Mx||^2
-    hvars.the_array[3] = vol*SQR(u_con_(m,3,k,j,i)); // ||My||^2
-    hvars.the_array[4] = vol*u_con_(m,4,k,j,i);      // ||Mz||^2
-    hvars.the_array[5] = vol*u_con_(m,5,k,j,i);      // ||Z||^2 (comes already squared)
-    hvars.the_array[6] = vol*SQR(u0_(m,I_Z4c_Theta_,k,j,i)); // ||Theta||^2
-    hvars.the_array[7] = vol*u_con_(m,6,k,j,i);      // ||C||^2 (comes already squared)
+    hvars.the_array[0] = vol*u_con_(m,0,k,j,i); // ||C||^2 (comes already squared)
+    hvars.the_array[1] = vol*SQR(u_con_(m,1,k,j,i)); //||H||^2
+    hvars.the_array[2] = vol*u_con_(m,2,k,j,i); // ||M||^2 (comes already squared)
+    hvars.the_array[3] = vol*u_con_(m,3,k,j,i); // ||Z||^2 (comes already squared)
+    hvars.the_array[4] = vol*SQR(u_con_(m,4,k,j,i));      // ||Mx||^2
+    hvars.the_array[5] = vol*SQR(u_con_(m,5,k,j,i));      // ||My||^2
+    hvars.the_array[6] = vol*SQR(u_con_(m,6,k,j,i));      // ||Mz||^2
+    hvars.the_array[7] = vol*SQR(u0_(m,I_Z4c_Theta_,k,j,i)); // ||Theta||^2
 
     // fill rest of the_array with zeros, if nhist < NHISTORY_VARIABLES
     for (int n=nhist_; n<NHISTORY_VARIABLES; ++n) {

src/pgen/z4c_gauge_wave.cpp

@@ -0,0 +1,119 @@
+//========================================================================================
+// 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_gauge_wave.cpp
+//! \brief z4c gauge wave test
+
+
+// C/C++ headers
+#include <algorithm> // min, max
+#include <cmath>     // sqrt()
+#include <cstdio>    // fopen(), fprintf(), freopen()
+#include <iostream>  // endl
+#include <sstream>   // stringstream
+#include <string>    // c_str()
+#include <limits>
+
+// Athena++ headers
+#include "athena.hpp"
+#include "globals.hpp"
+#include "parameter_input.hpp"
+#include "coordinates/cell_locations.hpp"
+#include "coordinates/adm.hpp"
+#include "mesh/mesh.hpp"
+#include "z4c/z4c.hpp"
+#include "driver/driver.hpp"
+#include "pgen/pgen.hpp"
+
+// function to compute errors in solution at end of run
+void Z4cGaugeWaveErrors(ParameterInput *pin, Mesh *pm);
+
+//----------------------------------------------------------------------------------------
+//! \fn void ProblemGenerator::UserProblem()
+//! \brief Sets initial conditions for gw linear wave tests
+
+void ProblemGenerator::UserProblem(ParameterInput *pin, const bool restart) {
+  if (restart)
+    return;
+
+  MeshBlockPack *pmbp = pmy_mesh_->pmb_pack;
+  if (pmbp->pz4c == nullptr) {
+    std::cout << "### FATAL ERROR in " << __FILE__ << " at line " << __LINE__ << std::endl
+              << "Z4c Wave test can only be run in Z4c, but no <z4c> block "
+              << "in input file" << std::endl;
+    exit(EXIT_FAILURE);
+  }
+
+  // Prepare Initial Data
+
+  // 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 isg = is-indcs.ng; int ieg = ie+indcs.ng;
+  int jsg = js-indcs.ng; int jeg = je+indcs.ng;
+  int ksg = ks-indcs.ng; int keg = ke+indcs.ng;
+
+  auto &pz4c = pmbp->pz4c;
+  auto &adm = pmbp->padm->adm;
+  Real x1size = pmy_mesh_->mesh_size.x1max - pmy_mesh_->mesh_size.x1min;
+  std::cout << x1size << std::endl;
+
+  // Wave amplitude
+  Real amp = pin->GetOrAddReal("problem", "amp", 0.001);
+
+  par_for("pgen_gauge_wave", DevExeSpace(), 0, (pmbp->nmb_thispack - 1),
+      ksg, keg, jsg, jeg, isg, ieg, KOKKOS_LAMBDA(int m, int k, int j, int i) {
+    Real &x1min = size.d_view(m).x1min;
+    Real &x1max = size.d_view(m).x1max;
+
+    int nx1 = indcs.nx1;
+
+    Real x1v = CellCenterX(i - is, nx1, x1min, x1max);
+    Real H = amp * sin(2 * M_PI * x1v / x1size);
+    Real dH_dt = - amp * 2 * M_PI / x1size *
+                    cos(2 * M_PI * x1v / x1size);
+
+    adm.g_dd(m,0,0,k,j,i) = 1 - H;
+    adm.g_dd(m,0,1,k,j,i) = 0;
+    adm.g_dd(m,0,2,k,j,i) = 0;
+    adm.g_dd(m,1,1,k,j,i) = 1;
+    adm.g_dd(m,1,2,k,j,i) = 0;
+    adm.g_dd(m,2,2,k,j,i) = 1;
+
+    adm.vK_dd(m,0,0,k,j,i) = 0.5 * dH_dt / sqrt(adm.g_dd(m,0,0,k,j,i));
+    adm.vK_dd(m,0,1,k,j,i) = 0;
+    adm.vK_dd(m,0,2,k,j,i) = 0;
+    adm.vK_dd(m,1,1,k,j,i) = 0;
+    adm.vK_dd(m,1,2,k,j,i) = 0;
+    adm.vK_dd(m,2,2,k,j,i) = 0;
+
+    adm.alpha(m,k,j,i) = sqrt(1 - H);
+  });
+  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);
+  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;
+  }
+  return;
+}