Implementation of a stable solver for the full (3+1D) vacuum Einstein equations

core/zero/gkyl_eqn_type.h

@@ -20,6 +20,8 @@ enum gkyl_eqn_type {
   GKYL_EQN_GR_MAXWELL_TETRAD, // General relativistic Maxwell equations in the tetrad basis.
   GKYL_EQN_GR_MEDIUM, // Coupled fluid-Einstein equations in plane-symmetric spacetimes.
   GKYL_EQN_GR_TWOFLUID, // General relativistic two-fluid equations.
+  GKYL_EQN_VACUUM_EINSTEIN, // Vacuum Einstein equations using the Bona-Masso formalism.
+  GKYL_EQN_VACUUM_EINSTEIN_CONFORMAL, // Vacuum Einstein equations using the conformal Bona-Masso formalism.
   GKYL_EQN_GR_TWOFLUID_TETRAD, // General relativistic two-fluid equations in the tetrad basis.
   GKYL_EQN_GR_MHD, // General relativistic magnetohydrodynamics equations.
   GKYL_EQN_GR_MHD_TETRAD, // General relativistic magnetohydrodynamics equations in the tetrad basis.

moments/apps/gkyl_moment.h

@@ -8,6 +8,7 @@
 #include <gkyl_util.h>
 #include <gkyl_wave_prop.h>
 #include <gkyl_wv_eqn.h>
+#include <gkyl_gr_spacetime.h>
 #include <gkyl_wv_embed_geo.h>
 
 #include <time.h>
@@ -62,6 +63,16 @@ struct gkyl_moment_species {
   double gr_twofluid_gas_gamma_ion; // Adiabatic index for ions in general relativistic two-fluid equations.
   double gr_twofluid_e_fact; // Electric field divergence error propagation speed for general relativistic two-fluid equations.
 
+  bool has_vacuum_einstein; // Run with vacuum Einstein sources using the Bona-Masso formalism.
+  double vacuum_einstein_excision_threshold; // Excision threshold (lapse) for vacuum Einstein equations using the Bona-Masso formalism.
+  enum gkyl_spacetime_slicing vacuum_einstein_spacetime_slicing; // Spacetime slicing condition for vacuum Einstein equations using the Bona-Masso formalism.
+  enum gkyl_spacetime_evolution vacuum_einstein_spacetime_evolution; // Spacetime evolution system for vacuum Einstein equations using the Bona-Masso formalism.
+
+  bool has_vacuum_einstein_conformal; // Run with vacuum Einstein sources using the conformal Bona-Masso formalism.
+  double vacuum_einstein_conformal_excision_threshold; // Excision threshold (lapse) for vacuum Einstein equations using the conformal Bona-Masso formalism.
+  enum gkyl_spacetime_slicing vacuum_einstein_conformal_spacetime_slicing; // Spacetime slicing condition for vacuum Einstein equations using the conformal Bona-Masso formalism.
+  enum gkyl_spacetime_evolution vacuum_einstein_conformal_spacetime_evolution; // Spacetime evolution system for vacuum Einstein equations using the conformal Bona-Masso formalism.
+
   bool has_gr_mhd; // Run with general relativistic source terms (general relativistic magnetohydrodynamics equations).
   double gr_mhd_gas_gamma; // Adiabatic index for general relativistic magnetohydrodynamics equations.
 

moments/apps/gkyl_moment_lw.h

@@ -67,4 +67,20 @@ gkyl_register_braginskii_types(lua_State *L);
 void
 gkyl_register_spacetime_gauge_types(lua_State *L);
 
+/**
+* Add spacetime slicing type flags for Einstein equations into Lua interpreter.
+*
+* @param L Lua state to use.
+ */
+void
+gkyl_register_spacetime_slicing_types(lua_State *L);
+
+/**
+* Add spacetime evolution type flags for Einstein equations into Lua interpreter.
+*
+* @param L Lua state to use.
+ */
+void
+gkyl_register_spacetime_evolution_types(lua_State *L);
+
 #endif

moments/apps/gkyl_moment_priv.h

@@ -44,6 +44,7 @@
 #include <gkyl_wv_maxwell.h>
 #include <gkyl_wv_mhd.h>
 #include <gkyl_wv_ten_moment.h>
+#include <gkyl_gr_spacetime.h>
 
 // number of components that various applied functions should return
 enum {
@@ -112,6 +113,16 @@ struct moment_species {
   double gr_twofluid_gas_gamma_ion; // Adiabatic index for ions in general relativistic two-fluid equations.
   double gr_twofluid_e_fact; // Electric field divergence error propagation speed for general relativistic two-fluid equations.
 
+  bool has_vacuum_einstein; // Run with vacuum Einstein sources using the Bona-Masso formalism.
+  double vacuum_einstein_excision_threshold; // Excision threshold (lapse) for vacuum Einstein equations using the Bona-Masso formalism.
+  enum gkyl_spacetime_slicing vacuum_einstein_spacetime_slicing; // Spacetime slicing condition for vacuum Einstein equations using the Bona-Masso formalism.
+  enum gkyl_spacetime_evolution vacuum_einstein_spacetime_evolution; // Spacetime evolution system for vacuum Einstein equations using the Bona-Masso formalism.
+
+  bool has_vacuum_einstein_conformal; // Run with vacuum Einstein sources using the conformal Bona-Masso formalism.
+  double vacuum_einstein_conformal_excision_threshold; // Excision threshold (lapse) for vacuum Einstein equations using the conformal Bona-Masso formalism.
+  enum gkyl_spacetime_slicing vacuum_einstein_conformal_spacetime_slicing; // Spacetime slicing condition for vacuum Einstein equations using the conformal Bona-Masso formalism.
+  enum gkyl_spacetime_evolution vacuum_einstein_conformal_spacetime_evolution; // Spacetime evolution system for vacuum Einstein equations using the conformal Bona-Masso formalism.
+
   bool has_gr_mhd; // Run with general relativistic source terms (general relativistic magnetohydrodynamics equations).
   double gr_mhd_gas_gamma; // Adiabatic index for general relativistic magnetohydrodynamics equations.
 

moments/apps/mom_coupling.c

@@ -164,6 +164,28 @@ moment_coupling_init(const struct gkyl_moment_app *app, struct moment_coupling *
     }
   }
 
+  src_inp.has_vacuum_einstein_sources = false;
+  for (int i = 0; i < app->num_species; i++) {
+    if (app->species[i].has_vacuum_einstein) {
+      src_inp.has_vacuum_einstein_sources = true;
+
+      src_inp.vacuum_einstein_excision_threshold = app->species[i].vacuum_einstein_excision_threshold;
+      src_inp.vacuum_einstein_spacetime_slicing = app->species[i].vacuum_einstein_spacetime_slicing;
+      src_inp.vacuum_einstein_spacetime_evolution = app->species[i].vacuum_einstein_spacetime_evolution;
+    }
+  }
+
+  src_inp.has_vacuum_einstein_conformal_sources = false;
+  for (int i = 0; i < app->num_species; i++) {
+    if (app->species[i].has_vacuum_einstein_conformal) {
+      src_inp.has_vacuum_einstein_conformal_sources = true;
+
+      src_inp.vacuum_einstein_conformal_excision_threshold = app->species[i].vacuum_einstein_conformal_excision_threshold;
+      src_inp.vacuum_einstein_conformal_spacetime_slicing = app->species[i].vacuum_einstein_conformal_spacetime_slicing;
+      src_inp.vacuum_einstein_conformal_spacetime_evolution = app->species[i].vacuum_einstein_conformal_spacetime_evolution;
+    }
+  }
+
   // save the use-rel bool
   src_inp.use_rel = use_rel;
 

moments/apps/mom_species.c

@@ -120,6 +120,26 @@ moment_species_init(const struct gkyl_moment *mom, const struct gkyl_moment_spec
     sp->gr_twofluid_e_fact = mom_sp->gr_twofluid_e_fact;
   }
 
+  sp->has_vacuum_einstein = false;
+  if (mom_sp->has_vacuum_einstein) {
+    sp->update_sources = true;
+    sp->has_vacuum_einstein = true;
+
+    sp->vacuum_einstein_excision_threshold = mom_sp->vacuum_einstein_excision_threshold;
+    sp->vacuum_einstein_spacetime_slicing = mom_sp->vacuum_einstein_spacetime_slicing;
+    sp->vacuum_einstein_spacetime_evolution = mom_sp->vacuum_einstein_spacetime_evolution;
+  }
+
+  sp->has_vacuum_einstein_conformal = false;
+  if (mom_sp->has_vacuum_einstein_conformal) {
+    sp->update_sources = true;
+    sp->has_vacuum_einstein_conformal = true;
+
+    sp->vacuum_einstein_conformal_excision_threshold = mom_sp->vacuum_einstein_conformal_excision_threshold;
+    sp->vacuum_einstein_conformal_spacetime_slicing = mom_sp->vacuum_einstein_conformal_spacetime_slicing;
+    sp->vacuum_einstein_conformal_spacetime_evolution = mom_sp->vacuum_einstein_conformal_spacetime_evolution;
+  }
+
   sp->has_gr_mhd = false;
   if (mom_sp->has_gr_mhd) {
     sp->update_sources = true;