Apply radiative boundary conditions to Z4c variables

NewRadX/configuration.ccl

@@ -1,3 +1,7 @@
 # Configuration definition for thorn NewRadX
 
 REQUIRES Loop
+
+PROVIDES NewRadX
+{
+}

Z4c/configuration.ccl

@@ -1,3 +1,3 @@
 # Configuration definitions for thorn Z4C
 
-REQUIRES Arith Loop
+REQUIRES Arith Loop NewRadX

Z4c/interface.ccl

@@ -14,6 +14,7 @@ USES INCLUDE HEADER: sum.hxx
 USES INCLUDE HEADER: vec.hxx
 USES INCLUDE HEADER: vect.hxx
 
+USES INCLUDE HEADER: newradx.hxx
 
 
 CCTK_INT FUNCTION GetCallFunctionCount()

Z4c/param.ccl

@@ -13,7 +13,6 @@ BOOLEAN calc_constraints "Calculate constraints" STEERABLE=recover
 } yes
 
 
-
 BOOLEAN set_Theta_zero "set Theta to zero, which converts Z4c to BSSN"
 {
 } no
@@ -57,3 +56,48 @@ CCTK_REAL epsdiss "Dissipation coefficient <arXiv:gr-qc/0610128>" STEERABLE=alwa
 {
   0.0:* :: ""
 } 0.32
+
+KEYWORD boundary_conditions "boundary conditions"{
+  "NewRadX" :: "radiative boundary conditions using NewRadX"
+  "CarpetX" :: "use CarpetX default boundary conditions"
+} "CarpetX"
+
+CCTK_INT n_chi "n power of outgoing boundary r^n fall off rate for chi"
+{
+  0:2           :: "1 is reasonable"
+} 1
+
+CCTK_INT n_gammat "n power of outgoing boundary r^n fall off rate for gammat_ij"
+{
+  0:2           :: "1 is reasonable"
+} 1
+
+CCTK_INT n_Kh "n power of outgoing boundary r^n fall off rate for Kh"
+{
+  0:2           :: "Maybe 1?"
+} 1
+
+CCTK_INT n_At "n power of outgoing boundary r^n fall off rate for At_ij"
+{
+  0:2           :: "Maybe 1?"
+} 1
+
+CCTK_INT n_Gamt "n power of outgoing boundary r^n fall off rate for Gamt^i"
+{
+  0:2           :: "Maybe 1?"
+} 1
+
+CCTK_INT n_Theta "n power of outgoing boundary r^n fall off rate for Theta"
+{
+  0:2           :: "Maybe 1?"
+} 1
+
+CCTK_INT n_alphaG "n power of outgoing boundary r^n fall off rate for alpha"
+{
+  0:2           :: "1 is my guess"
+} 1
+
+CCTK_INT n_betaG "n power of outgoing boundary r^n fall off rate for beta"
+{
+  0:2           :: "1 is my guess"
+} 1

Z4c/schedule.ccl

@@ -247,12 +247,43 @@ SCHEDULE Z4c_RHS IN Z4c_RHSGroup
   WRITES: Theta_rhs(interior)
   WRITES: alphaG_rhs(interior)
   WRITES: betaG_rhs(interior)
-  # SYNC: chi_rhs
-  # SYNC: gamma_tilde_rhs
-  # SYNC: K_hat_rhs
-  # SYNC: A_tilde_rhs
-  # SYNC: Gam_tilde_rhs
-  # SYNC: Theta_rhs
-  # SYNC: alphaG_rhs
-  # SYNC: betaG_rhs
+  SYNC: chi_rhs
+  SYNC: gamma_tilde_rhs
+  SYNC: K_hat_rhs
+  SYNC: A_tilde_rhs
+  SYNC: Gam_tilde_rhs
+  SYNC: Theta_rhs
+  SYNC: alphaG_rhs
+  SYNC: betaG_rhs
 } "Calculate Z4c RHS"
+
+if (CCTK_Equals(boundary_conditions, "NewRadX")) {
+  SCHEDULE Z4c_apply_newradx_boundary_conditions IN Z4c_RHSGroup AFTER Z4c_RHS
+  {
+    LANG: C
+    READS: chi(everywhere)
+    READS: gamma_tilde(everywhere)
+    READS: K_hat(everywhere)
+    READS: A_tilde(everywhere)
+    READS: Gam_tilde(everywhere)
+    READS: Theta(everywhere)
+    READS: alphaG(everywhere)
+    READS: betaG(everywhere)
+    READS: chi_rhs(everywhere)
+    READS: gamma_tilde_rhs(everywhere)
+    READS: K_hat_rhs(everywhere)
+    READS: A_tilde_rhs(everywhere)
+    READS: Gam_tilde_rhs(everywhere)
+    READS: Theta_rhs(everywhere)
+    READS: alphaG_rhs(everywhere)
+    READS: betaG_rhs(everywhere)
+    WRITES: chi_rhs(interior)
+    WRITES: gamma_tilde_rhs(interior)
+    WRITES: K_hat_rhs(interior)
+    WRITES: A_tilde_rhs(interior)
+    WRITES: Gam_tilde_rhs(interior)
+    WRITES: Theta_rhs(interior)
+    WRITES: alphaG_rhs(interior)
+    WRITES: betaG_rhs(interior)
+  } "Apply radiative boundary conditions to Z4c RHS variables using NewRadX"
+}

Z4c/src/apply_newrad.cxx

@@ -0,0 +1,41 @@
+// Apply NewRadX boundary conditions to Z4c variables
+//=============================================================================
+
+#include <cctk.h>
+#include <cctk_Arguments.h>
+#include <cctk_Parameters.h>
+#include <newradx.hxx>
+
+using namespace NewRadX;
+
+namespace Z4c {
+
+extern "C" void Z4c_apply_newradx_boundary_conditions(CCTK_ARGUMENTS) {
+  DECLARE_CCTK_ARGUMENTSX_Z4c_apply_newradx_boundary_conditions;
+  DECLARE_CCTK_PARAMETERS;
+
+  NewRadX_Apply(cctkGH, chi, chi_rhs, 1, 1, n_chi);
+  NewRadX_Apply(cctkGH, gammatxx, gammatxx_rhs, 1, 1, n_gammat);
+  NewRadX_Apply(cctkGH, gammatxy, gammatxy_rhs, 0, 1, n_gammat);
+  NewRadX_Apply(cctkGH, gammatxz, gammatxz_rhs, 0, 1, n_gammat);
+  NewRadX_Apply(cctkGH, gammatyy, gammatyy_rhs, 1, 1, n_gammat);
+  NewRadX_Apply(cctkGH, gammatyz, gammatyz_rhs, 0, 1, n_gammat);
+  NewRadX_Apply(cctkGH, gammatzz, gammatzz_rhs, 1, 1, n_gammat);
+  NewRadX_Apply(cctkGH, Kh, Kh_rhs, 0, 1, n_Kh);
+  NewRadX_Apply(cctkGH, Atxx, Atxx_rhs, 0, 1, n_At);
+  NewRadX_Apply(cctkGH, Atxy, Atxy_rhs, 0, 1, n_At);
+  NewRadX_Apply(cctkGH, Atxz, Atxz_rhs, 0, 1, n_At);
+  NewRadX_Apply(cctkGH, Atyy, Atyy_rhs, 0, 1, n_At);
+  NewRadX_Apply(cctkGH, Atyz, Atyz_rhs, 0, 1, n_At);
+  NewRadX_Apply(cctkGH, Atzz, Atzz_rhs, 0, 1, n_At);
+  NewRadX_Apply(cctkGH, Gamtx, Gamtx_rhs, 0, 1, n_Gamt);
+  NewRadX_Apply(cctkGH, Gamty, Gamty_rhs, 0, 1, n_Gamt);
+  NewRadX_Apply(cctkGH, Gamtz, Gamtz_rhs, 0, 1, n_Gamt);
+  NewRadX_Apply(cctkGH, Theta, Theta_rhs, 0, 1, n_Theta);
+  NewRadX_Apply(cctkGH, alphaG, alphaG_rhs, 1, 1, n_alphaG);
+  NewRadX_Apply(cctkGH, betaGx, betaGx_rhs, 0, 1, n_betaG);
+  NewRadX_Apply(cctkGH, betaGy, betaGy_rhs, 0, 1, n_betaG);
+  NewRadX_Apply(cctkGH, betaGz, betaGz_rhs, 0, 1, n_betaG);
+}
+
+} // namespace Z4c

Z4c/src/make.code.defn

@@ -4,6 +4,7 @@
 SRCS =						\
 	adm.cxx					\
 	adm2.cxx				\
+	apply_newrad.cxx			\
 	constraints.cxx				\
 	enforce.cxx				\
 	initial1.cxx				\

Z4c/test/qc0.par

@@ -0,0 +1,115 @@
+
+# run.me:
+# run.cores: 40
+# run.memory: 1.0e9
+# run.time: 7200.0
+
+ActiveThorns = "
+ADMBaseX
+CarpetX
+CoordinatesX
+ErrorEstimator
+IOUtil
+ODESolvers
+TmunuBaseX
+TwoPuncturesX
+Z4c
+"
+
+Cactus::cctk_show_schedule = no
+
+Cactus::presync_mode = "mixed-error"
+
+Cactus::terminate = "iteration"
+Cactus::cctk_itlast = 3
+
+CarpetX::verbose = no
+CarpetX::poison_undefined_values = no
+
+CarpetX::xmin = -20.0
+CarpetX::ymin = -20.0
+CarpetX::zmin = -20.0
+
+CarpetX::xmax = 20.0
+CarpetX::ymax = 20.0
+CarpetX::zmax = 20.0
+
+CarpetX::ncells_x = 48
+CarpetX::ncells_y = 48
+CarpetX::ncells_z = 24
+
+CarpetX::max_tile_size_x = 1024000
+CarpetX::max_tile_size_y = 4
+CarpetX::max_tile_size_z = 4
+
+CarpetX::reflection_z = yes
+
+CarpetX::boundary_x = "dirichlet"
+CarpetX::boundary_y = "dirichlet"
+CarpetX::boundary_upper_x = "dirichlet"
+CarpetX::boundary_upper_y = "dirichlet"
+CarpetX::boundary_upper_z = "dirichlet"
+
+CarpetX::ghost_size = 3
+
+CarpetX::max_num_levels = 1
+CarpetX::regrid_every = 16
+CarpetX::regrid_error_threshold = 0.03125
+
+ErrorEstimator::region_shape = "cube"
+ErrorEstimator::scale_by_resolution = yes
+
+CarpetX::prolongation_type = "ddf"
+CarpetX::prolongation_order = 5
+
+ODESolvers::verbose = no
+ODESolvers::method = "RK4"
+CarpetX::dtfac = 0.25
+
+ADMBaseX::initial_data = "TwoPunctures"
+ADMBaseX::initial_lapse = "TwoPunctures-averaged"
+
+# QC-0 setup
+TwoPuncturesX::par_b =  1.168642873
+TwoPuncturesX::par_m_plus =  0.453
+TwoPuncturesX::par_m_minus =  0.453
+TwoPuncturesX::par_P_plus [1] = +0.3331917498
+TwoPuncturesX::par_P_minus[1] = -0.3331917498
+
+TwoPuncturesX::npoints_A = 16
+TwoPuncturesX::npoints_B = 16
+TwoPuncturesX::npoints_phi = 8
+
+TwoPuncturesX::grid_setup_method = "evaluation"
+
+TwoPuncturesX::TP_epsilon = 1.0e-2
+TwoPuncturesX::TP_Tiny    = 1.0e-2
+TwoPuncturesX::adm_tol = 1e-6
+TwoPuncturesX::Newton_tol = 1e-6
+
+TwoPuncturesX::verbose = yes
+
+Z4c::calc_ADM_vars = yes             
+Z4c::calc_ADMRHS_vars = yes
+Z4c::calc_constraints = yes
+
+Z4c::chi_floor = 1.0e-6
+Z4c::alphaG_floor = 1.0e-8
+Z4c::epsdiss = 0.32
+Z4c::boundary_conditions = "NewRadX"
+
+IO::out_dir = $parfile
+IO::out_every = 1
+IO::parfile_write = no
+
+CarpetX::out_metadata = no
+CarpetX::out_norm_vars = " "
+CarpetX::out_norm_omit_unstable = yes
+
+CarpetX::out_tsv_vars = "
+ADMBaseX::metric
+ADMBaseX::curv
+ADMBaseX::lapse
+ADMBaseX::shift
+"
+

Z4c/test/qc0/TwoPuncturesX.bbh

@@ -0,0 +1,38 @@
+# ==================================
+# Numerical Relativity Metadata file
+# ==================================
+#
+# This file contains information about the simulation provided by the
+# TwoPuncturesX thorn.  The format is described in the NR Data Format Document
+# http://arxiv.org/abs/0709.0093 [draft SVN r707].
+
+[metadata]
+initial-ADM-energy            = 1.0078763090171796
+initial-ADM-angular-momentumx = 0
+initial-ADM-angular-momentumy = 0
+initial-ADM-angular-momentumz = 0.77876432749233837
+initial-separation            = 2.3372857460000001
+initial-data-type             = Bowen-York
+initial-data-bibtex-keys      = Bowen:1980yu Brandt:1997tf Ansorg:2004ds
+initial-bh-position1x         = 1.168642873
+initial-bh-position1y         = 0
+initial-bh-position1z         = 0
+initial-bh-position2x         = -1.168642873
+initial-bh-position2y         = 0
+initial-bh-position2z         = 0
+initial-bh-momentum1x         = 0
+initial-bh-momentum1y         = 0.33319174979999999
+initial-bh-momentum1z         = 0
+initial-bh-momentum2x         = 0
+initial-bh-momentum2y         = -0.33319174979999999
+initial-bh-momentum2z         = 0
+initial-bh-spin1x             = 0
+initial-bh-spin1y             = 0
+initial-bh-spin1z             = 0
+initial-bh-spin2x             = 0
+initial-bh-spin2y             = 0
+initial-bh-spin2z             = 0
+initial-bh-puncture-adm-mass1 = 0.51681958833218156
+initial-bh-puncture-adm-mass2 = 0.51681958829931418
+initial-bh-puncture-bare-mass1 = 0.45300000000000001
+initial-bh-puncture-bare-mass2 = 0.45300000000000001