Merge pull request #19 from jaykalinani/frontier

Bug fixes and enhancements from Frontier hackathon

AsterX/par/magTOV_AMR_frontier.par

@@ -0,0 +1,150 @@
+###############################
+# Simple test of magnetised TOV neutron star
+# Same neutron star as the gallery example
+# K=100, rho_c = 1.28e-3 => M = 1.4, M_b = 1.506
+# evolve for iterations = 10
+##############################
+ActiveThorns = "
+    ADMBase
+    CarpetX
+    ErrorEstimator
+    HydroBase
+    IOUtil
+    ODESolvers
+    TimerReport
+    TmunuBase
+    AsterX
+    AsterSeeds
+    Z4c
+"
+
+$nlevels	= 4
+$ncells 	= 256
+
+CarpetX::verbose = no
+
+Cactus::presync_mode = "mixed-error"
+CarpetX::poison_undefined_values = no
+
+CarpetX::xmin = -20
+CarpetX::ymin = -20
+CarpetX::zmin = -20
+
+CarpetX::xmax = 20
+CarpetX::ymax = 20
+CarpetX::zmax = 20
+
+CarpetX::ncells_x = $ncells
+CarpetX::ncells_y = $ncells
+CarpetX::ncells_z = $ncells
+
+CarpetX::boundary_x =  "neumann"
+CarpetX::boundary_y =  "neumann"
+CarpetX::boundary_z =  "neumann"
+CarpetX::boundary_upper_x =  "neumann"
+CarpetX::boundary_upper_y =  "neumann"
+CarpetX::boundary_upper_z =  "neumann"
+
+CarpetX::max_num_levels = $nlevels
+CarpetX::regrid_every = 16
+CarpetX::blocking_factor_x = 8
+CarpetX::blocking_factor_y = 8
+CarpetX::blocking_factor_z = 8
+
+
+# the regrid_error_threshold should be in the same units of the 
+# parameter in comparison, see AsterX/src/estimate_error.cxx 
+# for which parameter to use
+CarpetX::regrid_error_threshold = 0.1 
+
+ErrorEstimator::region_shape = "cube"
+ErrorEstimator::scale_by_resolution = yes
+
+CarpetX::prolongation_type = "ddf"
+CarpetX::ghost_size = 3
+CarpetX::dtfac = 0.25
+
+#ADMBase::set_adm_variables_during_evolution = "yes"
+ADMBase::initial_data       = "tov"
+ADMBase::initial_lapse      = "tov"
+ADMBase::initial_shift      = "tov"
+ADMBase::initial_dtlapse    = "zero"
+ADMBase::initial_dtshift    = "zero"
+
+Z4c::calc_ADM_vars = yes                  # for Weyl
+Z4c::calc_ADMRHS_vars = no  #yes               # for Weyl
+Z4c::calc_constraints = yes
+Z4c::chi_floor = 1.0e-6
+Z4c::alphaG_floor = 1.0e-8
+Z4c::epsdiss = 0.32
+
+
+ActiveThorns = "TOVSolver"
+TOVSolver::TOV_Rho_Central[0] = 1.28e-3
+TOVSolver::TOV_Gamma          = 2.0
+TOVSolver::TOV_K              = 100.0
+TOVSolver::TOV_Cowling = no
+
+
+AsterSeeds::test_type = "3DTest"
+AsterSeeds::test_case = "magTOV"
+AsterSeeds::Afield_config = "internal dipole"
+AsterSeeds::Ab = 10000.0
+AsterSeeds::press_cut = 0.04
+AsterSeeds::press_max = 1.638e-4
+AsterSeeds::Avec_kappa = 2.0
+
+#AsterSeeds::Afield_config = "external dipole"
+#AsterSeeds::B0 = 1e-7
+#AsterSeeds::r0 = 5.0
+
+AsterX::rho_abs_min = 1e-13
+AsterX::atmo_tol = 1e-3
+ReconX::reconstruction_method = "ppm"
+AsterX::flux_type = "HLLE"
+AsterX::max_iter = 100
+AsterX::c2p_tol = 1e-8
+AsterX::vector_potential_gauge = "algebraic"
+AsterX::local_estimate_error = no
+
+Con2PrimFactory::unit_test = "no"
+
+EOSX::evol_eos_name = "IdealGas"
+EOSX::gl_gamma = 2
+EOSX::poly_gamma = 2.0
+EOSX::poly_k = 100
+
+Cactus::terminate = "iteration"
+Cactus::cctk_itlast = 10
+ODESolvers::method = "RK4"
+
+IO::out_dir = $parfile
+IO::out_every = 10
+CarpetX::out_silo_vars = "
+    HydroBase::Bvec
+    HydroBase::rho
+    HydroBase::vel
+    HydroBase::eps
+    HydroBase::press
+    CarpetX::regrid_error
+    ADMBase::lapse
+    ADMBase::shift
+    ADMBase::metric
+"
+CarpetX::out_tsv_vars = "
+    AsterX::dB
+    HydroBase::Bvec
+    HydroBase::rho
+    HydroBase::vel
+    HydroBase::press
+    CarpetX::regrid_error
+    ADMBase::lapse
+    ADMBase::shift
+    ADMBase::metric
+"
+TimerReport::out_every = 10
+TimerReport::out_filename = "TimerReport"
+TimerReport::output_all_timers_together = yes
+TimerReport::output_all_timers_readable = yes
+TimerReport::n_top_timers = 50
+

AsterX/par/magTOV_unigrid_frontier.par

@@ -0,0 +1,145 @@
+###############################
+# Simple test of magnetised TOV neutron star
+# Same neutron star as the gallery example
+# K=100, rho_c = 1.28e-3 => M = 1.4, M_b = 1.506
+# evolve for iterations = 10
+##############################
+ActiveThorns = "
+    ADMBase
+    CarpetX
+    HydroBase
+    IOUtil
+    ODESolvers
+    TimerReport
+    TmunuBase
+    AsterX
+    AsterSeeds
+    Z4c
+"
+
+$nlevels	= 1
+$ncells 	=  512
+
+CarpetX::verbose = no
+
+Cactus::presync_mode = "mixed-error"
+CarpetX::poison_undefined_values = no
+
+CarpetX::xmin = -20
+CarpetX::ymin = -20
+CarpetX::zmin = -20
+
+CarpetX::xmax = 20
+CarpetX::ymax = 20
+CarpetX::zmax = 20
+
+CarpetX::ncells_x = $ncells
+CarpetX::ncells_y = $ncells
+CarpetX::ncells_z = $ncells
+
+CarpetX::boundary_x =  "neumann"
+CarpetX::boundary_y =  "neumann"
+CarpetX::boundary_z =  "neumann"
+CarpetX::boundary_upper_x =  "neumann"
+CarpetX::boundary_upper_y =  "neumann"
+CarpetX::boundary_upper_z =  "neumann"
+
+CarpetX::max_num_levels = $nlevels
+CarpetX::regrid_every = 100000
+CarpetX::blocking_factor_x = 1
+CarpetX::blocking_factor_y = 1
+CarpetX::blocking_factor_z = 1
+
+
+# the regrid_error_threshold should be in the same units of the 
+# parameter in comparison, see AsterX/src/estimate_error.cxx 
+# for which parameter to use
+CarpetX::regrid_error_threshold = 5.0e-5
+
+CarpetX::prolongation_type = "ddf"
+CarpetX::ghost_size = 3
+CarpetX::dtfac = 0.25
+
+#ADMBase::set_adm_variables_during_evolution = "yes"
+ADMBase::initial_data       = "tov"
+ADMBase::initial_lapse      = "tov"
+ADMBase::initial_shift      = "tov"
+ADMBase::initial_dtlapse    = "zero"
+ADMBase::initial_dtshift    = "zero"
+
+Z4c::calc_ADM_vars = yes                  # for Weyl
+Z4c::calc_ADMRHS_vars = no  #yes               # for Weyl
+Z4c::calc_constraints = yes
+Z4c::chi_floor = 1.0e-6
+Z4c::alphaG_floor = 1.0e-8
+Z4c::epsdiss = 0.32
+
+
+ActiveThorns = "TOVSolver"
+TOVSolver::TOV_Rho_Central[0] = 1.28e-3
+TOVSolver::TOV_Gamma          = 2.0
+TOVSolver::TOV_K              = 100.0
+TOVSolver::TOV_Cowling = no
+
+
+AsterSeeds::test_type = "3DTest"
+AsterSeeds::test_case = "magTOV"
+AsterSeeds::Afield_config = "internal dipole"
+AsterSeeds::Ab = 10000.0
+AsterSeeds::press_cut = 0.04
+AsterSeeds::press_max = 1.638e-4
+AsterSeeds::Avec_kappa = 2.0
+
+#AsterSeeds::Afield_config = "external dipole"
+#AsterSeeds::B0 = 1e-7
+#AsterSeeds::r0 = 5.0
+
+AsterX::rho_abs_min = 1e-13
+AsterX::atmo_tol = 1e-3
+ReconX::reconstruction_method = "ppm"
+AsterX::flux_type = "HLLE"
+AsterX::max_iter = 100
+AsterX::c2p_tol = 1e-8
+AsterX::vector_potential_gauge = "algebraic"
+
+Con2PrimFactory::unit_test = "no"
+
+EOSX::evol_eos_name = "IdealGas"
+EOSX::gl_gamma = 2
+EOSX::poly_gamma = 2.0
+EOSX::poly_k = 100
+
+Cactus::terminate = "iteration"
+Cactus::cctk_itlast = 10
+ODESolvers::method = "RK4"
+
+IO::out_dir = $parfile
+IO::out_every = 10
+CarpetX::out_silo_vars = "
+    HydroBase::Bvec
+    HydroBase::rho
+    HydroBase::vel
+    HydroBase::eps
+    HydroBase::press
+    CarpetX::regrid_error
+    ADMBase::lapse
+    ADMBase::shift
+    ADMBase::metric
+"
+CarpetX::out_tsv_vars = "
+    AsterX::dB
+    HydroBase::Bvec
+    HydroBase::rho
+    HydroBase::vel
+    HydroBase::press
+    CarpetX::regrid_error
+    ADMBase::lapse
+    ADMBase::shift
+    ADMBase::metric
+"
+TimerReport::out_every = 10
+TimerReport::out_filename = "TimerReport"
+TimerReport::output_all_timers_together = yes
+TimerReport::output_all_timers_readable = yes
+TimerReport::n_top_timers = 50
+

AsterX/param.ccl

@@ -24,6 +24,9 @@ KEYWORD flux_type "Flux solver" STEERABLE=never
   "HLLE" :: ""
 } "LxF"
 
+BOOLEAN local_estimate_error "Use error estimation criteria of this thorn"
+{
+} yes
 
 KEYWORD regrid_method "Regrid method" STEERABLE=never
 {

AsterX/schedule.ccl

@@ -22,12 +22,14 @@ if(unit_test)
   } "Self-test"
 }
 
-SCHEDULE AsterX_EstimateError_Setup AT wragh
-{
-  LANG: C
-  OPTIONS: GLOBAL
-} "set up for calculating regrid error"
-
+if(local_estimate_error)
+{  
+  SCHEDULE AsterX_EstimateError_Setup AT wragh
+  {
+    LANG: C
+    OPTIONS: GLOBAL
+  } "set up for calculating regrid error"
+}
 # Initial conditions
 
 SCHEDULE AsterX_ComputedBstagFromA AT initial AFTER HydroBase_PostInitial BEFORE AsterX_ComputedBFromdBstag
@@ -66,13 +68,15 @@ SCHEDULE AsterX_Prim2Con_Initial AT initial AFTER HydroBase_PostInitial
   SYNC: saved_prims
 } "Compute conserved variables from primitive variables at initial"
 
-SCHEDULE AsterX_EstimateError IN ODESolvers_EstimateError
+if(local_estimate_error)
 {
-  LANG: C
-  READS: HydroBase::rho(everywhere) HydroBase::vel(everywhere) HydroBase::eps(everywhere) HydroBase::press(everywhere)
-  WRITES: CarpetX::regrid_error(interior)
-} "Estimate local error for regridding initial conditions"
-
+  SCHEDULE AsterX_EstimateError IN ODESolvers_EstimateError
+  {
+    LANG: C
+    READS: HydroBase::rho(everywhere) HydroBase::vel(everywhere) HydroBase::eps(everywhere) HydroBase::press(everywhere)
+    WRITES: CarpetX::regrid_error(interior)
+  } "Estimate local error for regridding initial conditions"
+}
 # Regridding
 
 SCHEDULE AsterX_Sync AT postregrid

AsterX/src/fd.hxx

@@ -13,6 +13,7 @@
 #include <algorithm>
 #include <array>
 #include <cmath>
+#include <type_traits>
 
 namespace AsterX {
 using namespace std;
@@ -45,9 +46,10 @@ calc_fd2_v2e(const GF3D2<const T> &gf, const PointDesc &p, const int dir) {
 }
 
 // FD2: vertex centered input, cell centered output
-template <typename T>
-CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline T
-calc_fd2_v2c(const GF3D2<const T> &gf, const PointDesc &p, int dir) {
+template <int FDORDER, typename T>
+CCTK_DEVICE CCTK_HOST
+    CCTK_ATTRIBUTE_ALWAYS_INLINE inline std::enable_if_t<FDORDER == 2, T>
+    calc_fd_v2c(const GF3D2<const T> &gf, const PointDesc &p, int dir) {
   T dgf1, dgf2, dgf3, dgf4;
   const int dir1 = (dir == 0) ? 1 : ((dir == 1) ? 2 : 0);
   const int dir2 = (dir == 0) ? 2 : ((dir == 1) ? 0 : 1);
@@ -72,9 +74,10 @@ calc_fd4_c2c(const GF3D2<const T> &gf, const PointDesc &p, const int dir) {
 
 // FD4: vertex centered input, cell centered output
 // Interpolation from edges to cell centers is 2nd order
-template <typename T>
-CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline T
-calc_fd4_v2c(const GF3D2<const T> &gf, const PointDesc &p, int dir) {
+template <int FDORDER, typename T>
+CCTK_DEVICE CCTK_HOST
+    CCTK_ATTRIBUTE_ALWAYS_INLINE inline std::enable_if_t<FDORDER == 4, T>
+    calc_fd_v2c(const GF3D2<const T> &gf, const PointDesc &p, int dir) {
   T dgf1, dgf2, dgf3, dgf4;
 
   const int dir1 = (dir == 0) ? 1 : ((dir == 1) ? 2 : 0);