From 049af2cc59378aa11c3277dd6b62958ac247702e Mon Sep 17 00:00:00 2001 From: rhaas Date: Tue, 15 Apr 2014 19:49:20 +0000 Subject: GRHydro: New reconstruction interface. All in template C++. Tested with TOV and compared to old implementation. Seems to work. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@622 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45 --- param.ccl | 4 + schedule.ccl | 18 +- src/GRHydro_Functions.h | 16 ++ src/GRHydro_PPMReconstruct_drv_opt.cc | 2 +- src/GRHydro_Reconstruct.cc | 504 ++++++++++++++++++++++++++++++++++ src/GRHydro_Reconstruct_drv_cxx.hh | 64 +++++ src/GRHydro_Reconstruct_drv_impl.hh | 330 ++++++++++++++++++++++ src/make.code.defn | 1 + 8 files changed, 934 insertions(+), 5 deletions(-) create mode 100644 src/GRHydro_Functions.h create mode 100644 src/GRHydro_Reconstruct.cc create mode 100644 src/GRHydro_Reconstruct_drv_cxx.hh create mode 100644 src/GRHydro_Reconstruct_drv_impl.hh diff --git a/param.ccl b/param.ccl index cb899c0..b0cd1d8 100644 --- a/param.ccl +++ b/param.ccl @@ -803,4 +803,8 @@ BOOLEAN constrain_to_1D "Set fluid velocities to zero for non-radial motion" +# use_cxx_code is used in schedule.ccl so cannot be STEERABLE=ALWAYS +BOOLEAN use_cxx_code "Use experimental C++ code?" STEERABLE = RECOVER +{ +} no diff --git a/schedule.ccl b/schedule.ccl index f1cfe7a..5e08cfb 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -893,11 +893,21 @@ if (CCTK_Equals(method_type, "RSA FV")) if (CCTK_Equals(GRHydro_eos_type,"General")) { - schedule Reconstruction IN FluxTerms AS Reconstruct - { - LANG: Fortran - } "Reconstruct the functions at the cell boundaries" + if (use_cxx_code) { + + schedule Reconstruction_cxx IN FluxTerms AS Reconstruct + { + LANG: C + } "Reconstruct the functions at the cell boundaries" + } else { + + schedule Reconstruction IN FluxTerms AS Reconstruct + { + LANG: Fortran + } "Reconstruct the functions at the cell boundaries" + + } } else if (CCTK_Equals(GRHydro_eos_type,"Polytype")) { diff --git a/src/GRHydro_Functions.h b/src/GRHydro_Functions.h new file mode 100644 index 0000000..be78afa --- /dev/null +++ b/src/GRHydro_Functions.h @@ -0,0 +1,16 @@ +#ifndef _GRHYDRO_FUNCTIONS_H_ +#define _GRHYDRO_FUNCTIONS_H_ + +#ifdef __cplusplus +extern "C" { +#endif + +#include "cctk.h" + +void Primitive2ConservativeC(CCTK_ARGUMENTS); + +#ifdef __cplusplus +} +#endif + +#endif // _GRHYDRO_FUNCTIONS_H_ diff --git a/src/GRHydro_PPMReconstruct_drv_opt.cc b/src/GRHydro_PPMReconstruct_drv_opt.cc index 90206c7..6cf2e3e 100644 --- a/src/GRHydro_PPMReconstruct_drv_opt.cc +++ b/src/GRHydro_PPMReconstruct_drv_opt.cc @@ -10,7 +10,7 @@ #include "SpaceMask.h" #include "GRHydro_PPM_opt.h" -#include "GRHydro_Reconstruct_drv_cxx.hh" +#include "GRHydro_Reconstruct_drv_impl.hh" #undef velx #undef vely diff --git a/src/GRHydro_Reconstruct.cc b/src/GRHydro_Reconstruct.cc new file mode 100644 index 0000000..4d2974d --- /dev/null +++ b/src/GRHydro_Reconstruct.cc @@ -0,0 +1,504 @@ + /*@@ + @file GRHydro_Reconstruct.cc + @date Thu Dec 2013 + @author Christian Reisswig + @desc + Wrapper routine to perform the reconstruction.using tmeplate metaprogramming + @enddesc + @@*/ + +#include +#include +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" +#include "cctk_Functions.h" + +// for calling Fortran scheduled routines from C +#include "cctki_FortranWrappers.h" + +#include "SpaceMask.h" + +#include "GRHydro_Functions.h" + +#include "GRHydro_Reconstruct_drv_cxx.hh" +#include "GRHydro_WENOReconstruct.hh" +#include "GRHydro_TrivialReconstruct.hh" +#include "GRHydro_TVDReconstruct.hh" +#include "GRHydro_MP5Reconstruct.hh" + +#define velx (&vup[0]) +#define vely (&vup[N]) +#define velz (&vup[2*N]) + +#define Bvecx (&Bprim[0]) +#define Bvecy (&Bprim[N]) +#define Bvecz (&Bprim[2*N]) + +/* + External routines +*/ +extern "C" +CCTK_INT GRHydro_UseGeneralCoordinates(const cGH * cctkGH); + +/* + This is needed right now because we want to also call this from Fotran +*/ +extern "C" +void CCTK_FCALL CCTK_FNAME(GRHydro_PPMReconstruct_drv_opt)(cGH const * const restrict & restrict cctkGH); + +extern "C" +void CCTK_FCALL CCTK_FNAME(primitive2conservativeM)(cGH const * const restrict & restrict cctkGH); + + +/** + This class selects a reconstruction variant. +*/ +template +struct reconstruct +{ + static inline void select(const bool do_mhd, + const bool do_Ye, + const bool do_temp, + const bool reconstruct_Wv, + const bool clean_divergence, + cGH const * const restrict & restrict cctkGH) + { + if (!do_mhd && !do_Ye && !do_temp && reconstruct_Wv) + GRHydro_Reconstruct_drv_cxx (cctkGH); + else if (do_mhd && !do_Ye && !do_temp && reconstruct_Wv && !clean_divergence) + GRHydro_Reconstruct_drv_cxx (cctkGH); + else if (do_mhd && do_Ye && do_temp && reconstruct_Wv && !clean_divergence) + GRHydro_Reconstruct_drv_cxx (cctkGH); + else if (do_mhd && !do_Ye && !do_temp && reconstruct_Wv && clean_divergence) + GRHydro_Reconstruct_drv_cxx (cctkGH); + else if (do_mhd && do_Ye && do_temp && reconstruct_Wv && clean_divergence) + GRHydro_Reconstruct_drv_cxx (cctkGH); + else if (!do_mhd && do_Ye && do_temp && reconstruct_Wv) + GRHydro_Reconstruct_drv_cxx (cctkGH); + else if (!do_mhd && !do_Ye && !do_temp && !reconstruct_Wv) + GRHydro_Reconstruct_drv_cxx (cctkGH); + else if (do_mhd && !do_Ye && !do_temp && !reconstruct_Wv && !clean_divergence) + GRHydro_Reconstruct_drv_cxx (cctkGH); + else if (do_mhd && do_Ye && do_temp && !reconstruct_Wv && !clean_divergence) + GRHydro_Reconstruct_drv_cxx (cctkGH); + else if (do_mhd && !do_Ye && !do_temp && !reconstruct_Wv && clean_divergence) + GRHydro_Reconstruct_drv_cxx (cctkGH); + else if (do_mhd && do_Ye && do_temp && !reconstruct_Wv && clean_divergence) + GRHydro_Reconstruct_drv_cxx (cctkGH); + else if (!do_mhd && do_Ye && do_temp && !reconstruct_Wv) + GRHydro_Reconstruct_drv_cxx (cctkGH); + else + assert(0&&"Unsupported template arguments"); + } +}; + + +/** + This is a new template metaprogramming way of calling various reconstruction methods. + It tries to avoid duplicated code as much as possible +*/ + +extern "C" void Reconstruction_cxx(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + const bool do_mhd = CCTK_EQUALS(Bvec_evolution_method, "GRHydro"); + const bool do_Ye = CCTK_EQUALS(Y_e_evolution_method, "GRHydro"); + const bool do_temp = CCTK_EQUALS(temperature_evolution_method, "GRHydro"); + const bool do_clean_divergence = clean_divergence; + + /* + Set up multipatch stuff + */ + + const CCTK_REAL * restrict vup; + const CCTK_REAL * restrict Bprim; + + //Multipatch related pointers + if(GRHydro_UseGeneralCoordinates(cctkGH)) { + vup=lvel; + Bprim=lBvec; + } else { + vup=vel; + Bprim=Bvec; + } + + const int nx=cctk_lsh[0]; + const int ny=cctk_lsh[1]; + const int nz=cctk_lsh[2]; + + const int N = nx*ny*nz; + + int type_bits; + int not_trivial; + int trivial; + + if (*flux_direction == 1) { + type_bits = SpaceMask_GetTypeBits("Hydro_RiemannProblemX"); + not_trivial = SpaceMask_GetStateBits("Hydro_RiemannProblemX","not_trivial"); + trivial = SpaceMask_GetStateBits("Hydro_RiemannProblemX","trivial"); + } else if (*flux_direction == 2) { + type_bits = SpaceMask_GetTypeBits("Hydro_RiemannProblemY"); + not_trivial = SpaceMask_GetStateBits("Hydro_RiemannProblemY","not_trivial"); + trivial = SpaceMask_GetStateBits("Hydro_RiemannProblemY","trivial"); + } else if (*flux_direction == 3) { + type_bits = SpaceMask_GetTypeBits("Hydro_RiemannProblemZ"); + not_trivial = SpaceMask_GetStateBits("Hydro_RiemannProblemZ","not_trivial"); + trivial = SpaceMask_GetStateBits("Hydro_RiemannProblemZ","trivial"); + } + + CCTK_REAL local_min_tracer = 1e42; + if (evolve_tracer) { + if (use_min_tracer) + local_min_tracer = min_tracer; + else + local_min_tracer = 0.0; + } + + /* + First initialize all variables by first order trivial reconstruction! + */ + + + if (!do_mhd && !do_Ye && !do_temp) + GRHydro_Reconstruct_drv_cxx (cctkGH); + else if (do_mhd && !do_Ye && !do_temp && !clean_divergence) + GRHydro_Reconstruct_drv_cxx (cctkGH); + else if (do_mhd && !do_Ye && !do_temp && clean_divergence) + GRHydro_Reconstruct_drv_cxx (cctkGH); + else if (do_mhd && do_Ye && do_temp && !clean_divergence) + GRHydro_Reconstruct_drv_cxx (cctkGH); + else if (do_mhd && do_Ye && do_temp && clean_divergence) + GRHydro_Reconstruct_drv_cxx (cctkGH); + else if (!do_mhd && do_Ye && do_temp) + GRHydro_Reconstruct_drv_cxx (cctkGH); + else + CCTK_ERROR("Don't know how to initialize reconstructed variables"); + + + /* + Now do the real reconstruction! + */ + + if (CCTK_EQUALS(recon_method,"tvd")) { + // this handles MHD and non-MHD + + if (CCTK_EQUALS(tvd_limiter, "minmod")) { + + reconstruct >::select(do_mhd, do_Ye, do_temp, reconstruct_Wv, clean_divergence, cctkGH); + + } else if (CCTK_EQUALS(tvd_limiter, "vanleerMC2")) { + + reconstruct >::select(do_mhd, do_Ye, do_temp, reconstruct_Wv, clean_divergence, cctkGH); + + } else if (CCTK_EQUALS(tvd_limiter, "superbee")) { + + reconstruct >::select(do_mhd, do_Ye, do_temp, reconstruct_Wv, clean_divergence, cctkGH); + + } else { + CCTK_ERROR("TVD limiter not recognized!"); + } + + } else if (CCTK_EQUALS(recon_method,"ppm")) { + // PPM is special. It cannot reconstruct individual functions. + // Rather, all function are reconstructed at once! + // Hence we have a specialized calling function. + + CCTK_FNAME(GRHydro_PPMReconstruct_drv_opt)(cctkGH); + + } else if (CCTK_EQUALS(recon_method,"eno")) { + // this handles MHD and non-MHD + + // TODO!! + CCTK_ERROR("Sorry. Eno reconstruction not yet implemented!"); + + } else if (CCTK_EQUALS(recon_method,"weno")) { + // this handles MHD and non-MHD + + if (weno_adaptive_epsilon) { + + reconstruct >::select(do_mhd, do_Ye, do_temp, reconstruct_Wv, clean_divergence, cctkGH); + + } else { + + reconstruct >::select(do_mhd, do_Ye, do_temp, reconstruct_Wv, clean_divergence, cctkGH); + + } + + } else if (CCTK_EQUALS(recon_method, "weno-z")) { + + reconstruct >::select(do_mhd, do_Ye, do_temp, reconstruct_Wv, clean_divergence, cctkGH); + + } else if (CCTK_EQUALS(recon_method,"mp5")) { + // this handles MHD and non-MHD + + if (mp5_adaptive_eps) + reconstruct >::select(do_mhd, do_Ye, do_temp, reconstruct_Wv, clean_divergence, cctkGH); + else + reconstruct >::select(do_mhd, do_Ye, do_temp, reconstruct_Wv, clean_divergence, cctkGH); + + } else + CCTK_ERROR("Reconstruction method not recognized!"); + + + /* + Here goes the excision treatment. + If any of the points in the reconstruction stencil is excised, + we switch to trivial reconstruction. + If the current point itself is excised, we mark the Riemann problem as "trivial". + If the current point is just one point away from an excision boundary, we use the current cell values + to provide an "extrapolated" guess for the reconstructed state at the excision boundary + */ + if (GRHydro_enable_internal_excision) { + + #pragma omp parallel for + for (int k=GRHydro_stencil-1; k < cctk_lsh[2] - GRHydro_stencil+1+transport_constraints*(1-*zoffset); ++k) + for (int j=GRHydro_stencil-1; j < cctk_lsh[1] - GRHydro_stencil+1+transport_constraints*(1-*yoffset); ++j) + for (int i=GRHydro_stencil-1; i < cctk_lsh[0] - GRHydro_stencil+1+transport_constraints*(1-*xoffset); ++i) { + + const int ijk = CCTK_GFINDEX3D(cctkGH, i, j, k); + const int ijk_m = CCTK_GFINDEX3D(cctkGH, i-*xoffset, j-*yoffset, k-*zoffset); + const int ijk_p = CCTK_GFINDEX3D(cctkGH, i+*xoffset, j+*yoffset, k+*zoffset); + + bool excise = false; + bool first_order = false; + bool extrapolate = false; + int extrap_ijk = 0; + + // Whenever we modify the previously reconstructed value, + // we need to change the trivial flag! + bool change_riemann_trivial_flag = false; + + // for an excised point in the reconstruction stencil radius of the current point, we need to switch down to first order + for (int l=-(GRHydro_stencil-1); l < GRHydro_stencil-1; ++l) { + const int ijk_l = CCTK_GFINDEX3D(cctkGH, i+*xoffset*l, j+*yoffset*l, k+*zoffset*l); + if (hydro_excision_mask[ijk_l]) { + first_order = true; + break; + } + } + + // if the current point is excised, then we will mark the Riemann problem as trivial. + if (hydro_excision_mask[ijk]) + excise = true; + + if (!excise && hydro_excision_mask[ijk_m]) { + extrapolate = true; + extrap_ijk = ijk_m; + } + + if (!excise && hydro_excision_mask[ijk_p]) { + extrapolate = true; + extrap_ijk = ijk_p; + } + + if (first_order) + { + rhominus[ijk] = rho[ijk]; + rhoplus[ijk] = rho[ijk]; + + epsminus[ijk] = eps[ijk]; + epsplus[ijk] = eps[ijk]; + + velxminus[ijk] = velx[ijk]; + velxplus[ijk] = velx[ijk]; + velyminus[ijk] = vely[ijk]; + velyplus[ijk] = vely[ijk]; + velzminus[ijk] = velz[ijk]; + velzplus[ijk] = velz[ijk]; + + if (do_Ye) { + Y_e_plus[ijk] = Y_e[ijk]; + Y_e_minus[ijk] = Y_e[ijk]; + } + if (do_temp) { + tempminus[ijk] = temperature[ijk]; + tempplus[ijk] = temperature[ijk]; + } + if (do_mhd) { + Bvecxminus[ijk] = Bvecx[ijk]; + Bvecxplus[ijk] = Bvecx[ijk]; + Bvecyminus[ijk] = Bvecy[ijk]; + Bvecyplus[ijk] = Bvecy[ijk]; + Bveczminus[ijk] = Bvecz[ijk]; + Bveczplus[ijk] = Bvecz[ijk]; + if (do_clean_divergence) { + psidcminus[ijk] = psidc[ijk]; + psidcplus[ijk] = psidc[ijk]; + } + } + if (evolve_tracer) { + for (int l=0; l < number_of_tracers; ++l) { + tracerminus[ijk + l*N] = tracer[ijk + l*N]; + tracerplus[ijk + l*N] = tracer[ijk + l*N]; + } + } + + // Riemann problem might not be trivial anymore! + change_riemann_trivial_flag = true; + } + + if (extrapolate) + { + rhominus[extrap_ijk] = rho[ijk]; + rhoplus[extrap_ijk] = rho[ijk]; + + epsminus[extrap_ijk] = eps[ijk]; + epsplus[extrap_ijk] = eps[ijk]; + + velxminus[extrap_ijk] = velx[ijk]; + velxplus[extrap_ijk] = velx[ijk]; + velyminus[extrap_ijk] = vely[ijk]; + velyplus[extrap_ijk] = vely[ijk]; + velzminus[extrap_ijk] = velz[ijk]; + velzplus[extrap_ijk] = velz[ijk]; + + if (do_Ye) { + Y_e_plus[extrap_ijk] = Y_e[ijk]; + Y_e_minus[extrap_ijk] = Y_e[ijk]; + } + if (do_temp) { + tempminus[extrap_ijk] = temperature[ijk]; + tempplus[extrap_ijk] = temperature[ijk]; + } + if (do_mhd) { + Bvecxminus[extrap_ijk] = Bvecx[ijk]; + Bvecxplus[extrap_ijk] = Bvecx[ijk]; + Bvecyminus[extrap_ijk] = Bvecy[ijk]; + Bvecyplus[extrap_ijk] = Bvecy[ijk]; + Bveczminus[extrap_ijk] = Bvecz[ijk]; + Bveczplus[extrap_ijk] = Bvecz[ijk]; + if (do_clean_divergence) { + psidcminus[extrap_ijk] = psidc[ijk]; + psidcplus[extrap_ijk] = psidc[ijk]; + } + } + if (evolve_tracer) { + for (int l=0; l < number_of_tracers; ++l) { + tracerminus[extrap_ijk + l*N] = tracer[ijk + l*N]; + tracerplus[extrap_ijk + l*N] = tracer[ijk + l*N]; + } + } + + // Riemann problem might not be trivial anymore! + change_riemann_trivial_flag = true; + } + + if (change_riemann_trivial_flag && !excise) { + SpaceMask_SetStateBits(space_mask, ijk_m, type_bits, not_trivial); + SpaceMask_SetStateBits(space_mask, ijk, type_bits, not_trivial); + } else if (excise) { + SpaceMask_SetStateBits(space_mask, ijk_m, type_bits, trivial); + SpaceMask_SetStateBits(space_mask, ijk, type_bits, trivial); + } + } + + + } + + + /* + Here goes the atmosphere treatment! + */ + #pragma omp parallel for + for (int k=GRHydro_stencil-1; k < cctk_lsh[2] - GRHydro_stencil+1+transport_constraints*(1-*zoffset); ++k) + for (int j=GRHydro_stencil-1; j < cctk_lsh[1] - GRHydro_stencil+1+transport_constraints*(1-*yoffset); ++j) + for (int i=GRHydro_stencil-1; i < cctk_lsh[0] - GRHydro_stencil+1+transport_constraints*(1-*xoffset); ++i) { + const int ijk = CCTK_GFINDEX3D(cctkGH, i, j, k); + const int ijk_m = CCTK_GFINDEX3D(cctkGH, i-*xoffset, j-*yoffset, k-*zoffset); + + // Whenever we modify a previously reconstructed value, + // we need to change the trivial flag! + bool change_riemann_trivial_flag = false; + + if (rhominus[ijk] < *GRHydro_rho_min || rhoplus[ijk] < *GRHydro_rho_min) + { + rhominus[ijk] = rho[ijk]; + rhoplus[ijk] = rho[ijk]; + + epsminus[ijk] = eps[ijk]; + epsplus[ijk] = eps[ijk]; + + velxminus[ijk] = velx[ijk]; + velxplus[ijk] = velx[ijk]; + velyminus[ijk] = vely[ijk]; + velyplus[ijk] = vely[ijk]; + velzminus[ijk] = velz[ijk]; + velzplus[ijk] = velz[ijk]; + // TODO: add reconstruct_Wv treatment + + if (do_Ye) { + Y_e_plus[ijk] = Y_e[ijk]; + Y_e_minus[ijk] = Y_e[ijk]; + } + if (do_temp) { + tempminus[ijk] = temperature[ijk]; + tempplus[ijk] = temperature[ijk]; + } + if (do_mhd) { + Bvecxminus[ijk] = Bvecx[ijk]; + Bvecxplus[ijk] = Bvecx[ijk]; + Bvecyminus[ijk] = Bvecy[ijk]; + Bvecyplus[ijk] = Bvecy[ijk]; + Bveczminus[ijk] = Bvecz[ijk]; + Bveczplus[ijk] = Bvecz[ijk]; + if (do_clean_divergence) { + psidcminus[ijk] = psidc[ijk]; + psidcplus[ijk] = psidc[ijk]; + } + } + if (evolve_tracer) { + for (int l=0; l < number_of_tracers; ++l) { + if (tracerminus[ijk + l*N] < local_min_tracer || tracerplus[ijk + l*N] < local_min_tracer) { + tracerminus[ijk + l*N] = tracer[ijk + l*N]; + tracerplus[ijk + l*N] = tracer[ijk + l*N]; + } + } + } + + change_riemann_trivial_flag = true; + } + + if (!do_temp) { + if (epsplus[ijk] < 0.0) { + epsplus[ijk] = eps[ijk]; + change_riemann_trivial_flag = true; + } + if (epsminus[ijk] < 0.0) { + epsminus[ijk] = eps[ijk]; + change_riemann_trivial_flag = true; + } + } + + if (change_riemann_trivial_flag) { + // Riemann problem might not be trivial anymore! + SpaceMask_SetStateBits(space_mask, ijk_m, type_bits, not_trivial); + SpaceMask_SetStateBits(space_mask, ijk, type_bits, not_trivial); + } + } + + + /* + Here goes the prim2con call for reconstructed variables! + */ + + if (CCTK_EQUALS(recon_vars, "primitive") || CCTK_EQUALS(recon_method, "ppm")) { + if (do_mhd) { + // the alternatives to avoid using internals of Cactus are currently + // to move this whole block into schedule.ccl or to construct a struct + // cFunction and call CallScheduledFunction + int (*wrapper)(void*,void*) = CCTKi_FortranWrapper(CCTK_THORNSTRING); + wrapper(cctkGH, (void*)CCTK_FNAME(primitive2conservativeM)); + } else { + Primitive2ConservativeC(CCTK_PASS_CTOC); + } + } else if (CCTK_EQUALS(recon_vars,"conservative")) { + CCTK_ERROR("Conservative reconstruction currently not supported!"); + } else + CCTK_ERROR("Variable type to reconstruct not recognized."); + + +} + diff --git a/src/GRHydro_Reconstruct_drv_cxx.hh b/src/GRHydro_Reconstruct_drv_cxx.hh new file mode 100644 index 0000000..57fb0bf --- /dev/null +++ b/src/GRHydro_Reconstruct_drv_cxx.hh @@ -0,0 +1,64 @@ +#ifndef _GRHYDRO_RECONSTRUCT_DRV_CXX_H +#define _GRHYDRO_RECONSTRUCT_DRV_CXX_H + +#include "cctk.h" + +/* + Cases that must be considered: + * basic hydro + * hydro + temperature + ye + * hydro + ye + * basic mhd + * mhd + temperature + ye + * mhd + ye + */ + +template // the reconstruction operator +void GRHydro_Reconstruct_drv_cxx(cGH const * const restrict & restrict cctkGH); + +// helper macro to instantiate all required permuations of the template options +// this must match GRHydro_Reconstruct's reconstruct::select routine +#define INSTANTIATE_RECONSTRUCTION_OPERATOR(RECONSTRUCT) \ +template void \ +GRHydro_Reconstruct_drv_cxx ( \ + cGH const * const restrict & restrict cctkGH); \ +template void \ +GRHydro_Reconstruct_drv_cxx ( \ + cGH const * const restrict & restrict cctkGH); \ +template void \ +GRHydro_Reconstruct_drv_cxx ( \ + cGH const * const restrict & restrict cctkGH); \ +template void \ +GRHydro_Reconstruct_drv_cxx ( \ + cGH const * const restrict & restrict cctkGH); \ +template void \ +GRHydro_Reconstruct_drv_cxx ( \ + cGH const * const restrict & restrict cctkGH); \ +template void \ +GRHydro_Reconstruct_drv_cxx ( \ + cGH const * const restrict & restrict cctkGH); \ +template void \ +GRHydro_Reconstruct_drv_cxx ( \ + cGH const * const restrict & restrict cctkGH); \ +template void \ +GRHydro_Reconstruct_drv_cxx ( \ + cGH const * const restrict & restrict cctkGH); \ +template void \ +GRHydro_Reconstruct_drv_cxx ( \ + cGH const * const restrict & restrict cctkGH); \ +template void \ +GRHydro_Reconstruct_drv_cxx ( \ + cGH const * const restrict & restrict cctkGH); \ +template void \ +GRHydro_Reconstruct_drv_cxx ( \ + cGH const * const restrict & restrict cctkGH); \ +template void \ +GRHydro_Reconstruct_drv_cxx ( \ + cGH const * const restrict & restrict cctkGH); + +#endif // _GRHYDRO_RECONSTRUCT_DRV_CXX_H diff --git a/src/GRHydro_Reconstruct_drv_impl.hh b/src/GRHydro_Reconstruct_drv_impl.hh new file mode 100644 index 0000000..07402db --- /dev/null +++ b/src/GRHydro_Reconstruct_drv_impl.hh @@ -0,0 +1,330 @@ +#ifndef _GRHYDRO_RECONSTRUCT_DRV_IMPL_H +#define _GRHYDRO_RECONSTRUCT_DRV_IMPL_H + +#include +#include +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#define MIN(a,b) (((a)<(b))?(a):(b)) +#define MAX(a,b) (((a)>(b))?(a):(b)) + +#include "SpaceMask.h" + +#include "GRHydro_Reconstruct_drv_cxx.hh" + +#define velx (&vup[0]) +#define vely (&vup[N]) +#define velz (&vup[2*N]) + +#define Bvecx (&Bprim[0]) +#define Bvecy (&Bprim[N]) +#define Bvecz (&Bprim[2*N]) + + +using namespace std; + +extern "C" CCTK_INT GRHydro_UseGeneralCoordinates(const cGH * cctkGH); + +/// transform back to vel if Wv was used +template +static inline void undo_Wv(const int nx, + CCTK_REAL* const velxminus, + CCTK_REAL* const velyminus, + CCTK_REAL* const velzminus, + CCTK_REAL* const velxplus, + CCTK_REAL* const velyplus, + CCTK_REAL* const velzplus, + const CCTK_REAL* const gxx, + const CCTK_REAL* const gxy, + const CCTK_REAL* const gxz, + const CCTK_REAL* const gyy, + const CCTK_REAL* const gyz, + const CCTK_REAL* const gzz, + const cGH* const cctkGH, + const int j, const int k) +{ + DECLARE_CCTK_PARAMETERS; + + // TODO: check index bounds!! + for (int i = GRHydro_stencil-1; i < nx - GRHydro_stencil+1; ++i) { + + const int ijk = dir == 0 ? CCTK_GFINDEX3D(cctkGH, i, j, k) : dir == 1 ? CCTK_GFINDEX3D(cctkGH, j, i, k) : CCTK_GFINDEX3D(cctkGH, j, k, i); + + // divide out the Loretnz factor obtained from w_lorentz = + // sqrt(1+g_{ij} w^i w^j) for both the + // plus and minus quantities this should by construction ensure + // that any Lorentz factor calculated from them later on is + // physical (ie. > 1.d0) + { + const int ijk_m = dir == 0 ? CCTK_GFINDEX3D(cctkGH, i-1, j, k) : dir == 1 ? CCTK_GFINDEX3D(cctkGH, j, i-1, k) : CCTK_GFINDEX3D(cctkGH, j, k, i-1); + + const CCTK_REAL agxx = 0.5*( gxx[ijk] + gxx[ijk_m] ); + const CCTK_REAL agxy = 0.5*( gxy[ijk] + gxy[ijk_m] ); + const CCTK_REAL agxz = 0.5*( gxz[ijk] + gxz[ijk_m] ); + const CCTK_REAL agyy = 0.5*( gyy[ijk] + gyy[ijk_m] ); + const CCTK_REAL agyz = 0.5*( gyz[ijk] + gyz[ijk_m] ); + const CCTK_REAL agzz = 0.5*( gzz[ijk] + gzz[ijk_m] ); + const CCTK_REAL w = sqrt( 1.0 + agxx*velxminus[ijk]*velxminus[ijk] + + agyy*velyminus[ijk]*velyminus[ijk] + + agzz*velzminus[ijk]*velzminus[ijk] + + 2.0*agxy*velxminus[ijk]*velyminus[ijk] + + 2.0*agxz*velxminus[ijk]*velzminus[ijk] + + 2.0*agyz*velyminus[ijk]*velzminus[ijk] ); + velxminus[ijk] = velxminus[ijk]/w; + velyminus[ijk] = velyminus[ijk]/w; + velzminus[ijk] = velzminus[ijk]/w; + } + + { + const int ijk_p = dir == 0 ? CCTK_GFINDEX3D(cctkGH, i+1, j, k) : dir == 1 ? CCTK_GFINDEX3D(cctkGH, j, i+1, k) : CCTK_GFINDEX3D(cctkGH, j, k, i+1); + + const CCTK_REAL agxx = 0.5*( gxx[ijk] + gxx[ijk_p] ); + const CCTK_REAL agxy = 0.5*( gxy[ijk] + gxy[ijk_p] ); + const CCTK_REAL agxz = 0.5*( gxz[ijk] + gxz[ijk_p] ); + const CCTK_REAL agyy = 0.5*( gyy[ijk] + gyy[ijk_p] ); + const CCTK_REAL agyz = 0.5*( gyz[ijk] + gyz[ijk_p] ); + const CCTK_REAL agzz = 0.5*( gzz[ijk] + gzz[ijk_p] ); + const CCTK_REAL w = sqrt( 1.0 + agxx*velxplus[ijk]*velxplus[ijk] + + agyy*velyplus[ijk]*velyplus[ijk] + + agzz*velzplus[ijk]*velzplus[ijk] + + 2.0*agxy*velxplus[ijk]*velyplus[ijk] + + 2.0*agxz*velxplus[ijk]*velzplus[ijk] + + 2.0*agyz*velyplus[ijk]*velzplus[ijk] ); + velxplus[ijk] = velxplus[ijk]/w; + velyplus[ijk] = velyplus[ijk]/w; + velzplus[ijk] = velzplus[ijk]/w; + } + } +} + + +/* + Cases that must be considered: + * basic hydro + * hydro + temperature + ye + * hydro + ye + * basic mhd + * mhd + temperature + ye + * mhd + ye + */ + +// TODO: remove all templatization from this reoutine since they are only used +// in the 2d loops and bloat the source code by a factor of 12 +// TODO: measure if 1d slicing slows things down since this would buy a factor +// of 3 in source code size again +template // the reconstruction operator +void GRHydro_Reconstruct_drv_cxx(cGH const * const restrict & restrict cctkGH) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + const CCTK_REAL * restrict vup; + const CCTK_REAL * restrict Bprim; + const CCTK_REAL * restrict g11; + const CCTK_REAL * restrict g12; + const CCTK_REAL * restrict g13; + const CCTK_REAL * restrict g22; + const CCTK_REAL * restrict g23; + const CCTK_REAL * restrict g33; + + //Multipatch related pointers + if(GRHydro_UseGeneralCoordinates(cctkGH)) { + vup=lvel; + g11=gaa; g12=gab; g13=gac; + g22=gbb; g23=gbc; + g33=gcc; + Bprim=lBvec; + } else { + vup=vel; + g11=gxx; g12=gxy; g13=gxz; + g22=gyy; g23=gyz; + g33=gzz; + Bprim=Bvec; + } + + const int nx=cctk_lsh[0]; + const int ny=cctk_lsh[1]; + const int nz=cctk_lsh[2]; + + const int N = nx*ny*nz; + + int type_bitsx,type_bitsy,type_bitsz; + // UNUSED: int trivialx,trivialy,trivialz; + int not_trivialx,not_trivialy,not_trivialz; + + type_bitsx = SpaceMask_GetTypeBits("Hydro_RiemannProblemX"); + // trivialx = SpaceMask_GetStateBits("Hydro_RiemannProblemX","trivial"); + not_trivialx = SpaceMask_GetStateBits("Hydro_RiemannProblemX","not_trivial"); + + type_bitsy = SpaceMask_GetTypeBits("Hydro_RiemannProblemY"); + // trivialy = SpaceMask_GetStateBits("Hydro_RiemannProblemY","trivial"); + not_trivialy = SpaceMask_GetStateBits("Hydro_RiemannProblemY","not_trivial"); + + type_bitsz = SpaceMask_GetTypeBits("Hydro_RiemannProblemZ"); + // trivialz = SpaceMask_GetStateBits("Hydro_RiemannProblemZ","trivial"); + not_trivialz = SpaceMask_GetStateBits("Hydro_RiemannProblemZ","not_trivial"); + + + vector Wvx(0); + vector Wvy(0); + vector Wvz(0); + // allocate temp memory for Wv reconstruction + // TODO: measure impact of C++ intializing the vector to zero and write + // custom wrapper around new/delete if required + // template class pod_vector { + // pod_vector() {data = NULL;}; + // ~pod_vecotr() {delete[] data;}; + // T& operator[](size_t i) {return data[i];}; + // void allocate(size_t sz) {data = new T[sz];}; + // private: + // pod_vector(const pod_vector&); + // pod_vector& operator=(const pod_vector&); + // }; + if (do_reconstruct_Wv) { + Wvx.resize(N, 0); + Wvy.resize(N, 0); + Wvz.resize(N, 0); + #pragma omp parallel for + for (int i=0; i < N; ++i) { + Wvx[i] = w_lorentz[i]*velx[i]; + Wvy[i] = w_lorentz[i]*vely[i]; + Wvz[i] = w_lorentz[i]*velz[i]; + } + } + +// Reconstruction starts: + if (flux_direction[0] == 1) { + #pragma omp parallel for + for(int k = GRHydro_stencil-1; k < nz - GRHydro_stencil+1+transport_constraints; ++k) { + for(int j = GRHydro_stencil-1; j < ny - GRHydro_stencil+1+transport_constraints; ++j) { + + RECONSTRUCT::template apply<0>(nx, rho, rhominus, rhoplus, cctkGH, j, k); + RECONSTRUCT::template apply<0>(nx, eps, epsminus, epsplus, cctkGH, j, k); + + if (do_reconstruct_Wv) { + RECONSTRUCT::template apply<0>(nx, &Wvx[0], velxminus, velxplus, cctkGH, j, k); + RECONSTRUCT::template apply<0>(nx, &Wvy[0], velyminus, velyplus, cctkGH, j, k); + RECONSTRUCT::template apply<0>(nx, &Wvz[0], velzminus, velzplus, cctkGH, j, k); + undo_Wv<0>(nx, velxminus, velyminus, velzminus, + velxplus, velyplus, velzplus, + g11, g12, g13, g22, g23, g33, + cctkGH, j, k); + } else { + RECONSTRUCT::template apply<0>(nx, velx, velxminus, velxplus, cctkGH, j, k); + RECONSTRUCT::template apply<0>(nx, vely, velyminus, velyplus, cctkGH, j, k); + RECONSTRUCT::template apply<0>(nx, velz, velzminus, velzplus, cctkGH, j, k); + } + + if (do_Ye) + RECONSTRUCT::template apply<0>(nx, Y_e, Y_e_minus, Y_e_plus, cctkGH, j, k); + if (do_temp) + RECONSTRUCT::template apply<0>(nx, temperature, tempminus, tempplus, cctkGH, j, k); + if (do_mhd) { + RECONSTRUCT::template apply<0>(nx, Bvecx, Bvecxminus, Bvecxplus, cctkGH, j, k); + RECONSTRUCT::template apply<0>(nx, Bvecy, Bvecyminus, Bvecyplus, cctkGH, j, k); + RECONSTRUCT::template apply<0>(nx, Bvecz, Bveczminus, Bveczplus, cctkGH, j, k); + if (do_clean_divergence) + RECONSTRUCT::template apply<0>(nx, psidc, psidcminus, psidcplus, cctkGH, j, k); + } + + for (int i=0; i(ny, rho, rhominus, rhoplus, cctkGH, j, k); + RECONSTRUCT::template apply<1>(ny, eps, epsminus, epsplus, cctkGH, j, k); + + if (do_reconstruct_Wv) { + RECONSTRUCT::template apply<1>(ny, &Wvx[0], velxminus, velxplus, cctkGH, j, k); + RECONSTRUCT::template apply<1>(ny, &Wvy[0], velyminus, velyplus, cctkGH, j, k); + RECONSTRUCT::template apply<1>(ny, &Wvz[0], velzminus, velzplus, cctkGH, j, k); + undo_Wv<1>(ny, velyminus, velzminus, velxminus, + velyplus, velzplus, velxplus, + g22, g23, g12, g33, g13, g11, + cctkGH, j, k); + } else { + RECONSTRUCT::template apply<1>(ny, velx, velxminus, velxplus, cctkGH, j, k); + RECONSTRUCT::template apply<1>(ny, vely, velyminus, velyplus, cctkGH, j, k); + RECONSTRUCT::template apply<1>(ny, velz, velzminus, velzplus, cctkGH, j, k); + } + + if (do_Ye) + RECONSTRUCT::template apply<1>(ny, Y_e, Y_e_minus, Y_e_plus, cctkGH, j, k); + if (do_temp) + RECONSTRUCT::template apply<1>(ny, temperature, tempminus, tempplus, cctkGH, j, k); + if (do_mhd) { + RECONSTRUCT::template apply<1>(ny, Bvecx, Bvecxminus, Bvecxplus, cctkGH, j, k); + RECONSTRUCT::template apply<1>(ny, Bvecy, Bvecyminus, Bvecyplus, cctkGH, j, k); + RECONSTRUCT::template apply<1>(ny, Bvecz, Bveczminus, Bveczplus, cctkGH, j, k); + if (do_clean_divergence) + RECONSTRUCT::template apply<1>(ny, psidc, psidcminus, psidcplus, cctkGH, j, k); + } + + for (int i=0; i(nz, rho, rhominus, rhoplus, cctkGH, j, k); + RECONSTRUCT::template apply<2>(nz, eps, epsminus, epsplus, cctkGH, j, k); + + if (do_reconstruct_Wv) { + RECONSTRUCT::template apply<2>(nz, &Wvx[0], velxminus, velxplus, cctkGH, j, k); + RECONSTRUCT::template apply<2>(nz, &Wvy[0], velyminus, velyplus, cctkGH, j, k); + RECONSTRUCT::template apply<2>(nz, &Wvz[0], velzminus, velzplus, cctkGH, j, k); + undo_Wv<2>(nz, velzminus, velxminus, velyminus, + velzplus, velxplus, velyplus, + g33, g13, g23, g11, g12, g22, + cctkGH, j, k); + } else { + RECONSTRUCT::template apply<2>(nz, velx, velxminus, velxplus, cctkGH, j, k); + RECONSTRUCT::template apply<2>(nz, vely, velyminus, velyplus, cctkGH, j, k); + RECONSTRUCT::template apply<2>(nz, velz, velzminus, velzplus, cctkGH, j, k); + } + + if (do_Ye) + RECONSTRUCT::template apply<2>(nz, Y_e, Y_e_minus, Y_e_plus, cctkGH, j, k); + if (do_temp) + RECONSTRUCT::template apply<2>(nz, temperature, tempminus, tempplus, cctkGH, j, k); + if (do_mhd) { + RECONSTRUCT::template apply<2>(nz, Bvecx, Bvecxminus, Bvecxplus, cctkGH, j, k); + RECONSTRUCT::template apply<2>(nz, Bvecy, Bvecyminus, Bvecyplus, cctkGH, j, k); + RECONSTRUCT::template apply<2>(nz, Bvecz, Bveczminus, Bveczplus, cctkGH, j, k); + if (do_clean_divergence) + RECONSTRUCT::template apply<2>(nz, psidc, psidcminus, psidcplus, cctkGH, j, k); + } + + for (int i=0; i