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

8 files changed, 934 insertions, 5 deletions

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 <iostream>
+#include <cassert>
+#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 <class RECONSTRUCT>
+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 <false, false, false, true, false, RECONSTRUCT > (cctkGH);
+       else if (do_mhd && !do_Ye && !do_temp && reconstruct_Wv && !clean_divergence)
+          GRHydro_Reconstruct_drv_cxx <true, false, false, true, false, RECONSTRUCT > (cctkGH);
+       else if (do_mhd && do_Ye && do_temp && reconstruct_Wv && !clean_divergence)
+          GRHydro_Reconstruct_drv_cxx <true, true, true, true, false, RECONSTRUCT  > (cctkGH);
+       else if (do_mhd && !do_Ye && !do_temp && reconstruct_Wv && clean_divergence)
+          GRHydro_Reconstruct_drv_cxx <true, false, false, true, true, RECONSTRUCT  > (cctkGH);
+       else if (do_mhd && do_Ye && do_temp && reconstruct_Wv && clean_divergence)
+          GRHydro_Reconstruct_drv_cxx <true, true, true, true, true, RECONSTRUCT  > (cctkGH);
+       else if (!do_mhd && do_Ye && do_temp && reconstruct_Wv)
+          GRHydro_Reconstruct_drv_cxx <false, true, true, true, false, RECONSTRUCT  > (cctkGH);
+       else if (!do_mhd && !do_Ye && !do_temp && !reconstruct_Wv)
+          GRHydro_Reconstruct_drv_cxx <false, false, false, false, false, RECONSTRUCT  > (cctkGH);
+       else if (do_mhd && !do_Ye && !do_temp && !reconstruct_Wv && !clean_divergence)
+          GRHydro_Reconstruct_drv_cxx <true, false, false, false, false, RECONSTRUCT  > (cctkGH);
+       else if (do_mhd && do_Ye && do_temp && !reconstruct_Wv && !clean_divergence)
+          GRHydro_Reconstruct_drv_cxx <true, true, true, false, false, RECONSTRUCT  > (cctkGH);
+       else if (do_mhd && !do_Ye && !do_temp && !reconstruct_Wv && clean_divergence)
+          GRHydro_Reconstruct_drv_cxx <true, false, false, false, true, RECONSTRUCT  > (cctkGH);
+       else if (do_mhd && do_Ye && do_temp && !reconstruct_Wv && clean_divergence)
+          GRHydro_Reconstruct_drv_cxx <true, true, true, false, true, RECONSTRUCT  > (cctkGH);
+       else if (!do_mhd && do_Ye && do_temp && !reconstruct_Wv)
+          GRHydro_Reconstruct_drv_cxx <false, true, true, false, false, RECONSTRUCT  > (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 <false, false, false, false, false, GRHydro_TrivialReconstruct1d > (cctkGH);
+  else if (do_mhd && !do_Ye && !do_temp && !clean_divergence)
+      GRHydro_Reconstruct_drv_cxx <true, false, false, false, false, GRHydro_TrivialReconstruct1d > (cctkGH);
+  else if (do_mhd && !do_Ye && !do_temp && clean_divergence)
+      GRHydro_Reconstruct_drv_cxx <true, false, false, false, true, GRHydro_TrivialReconstruct1d > (cctkGH);
+  else if (do_mhd && do_Ye && do_temp && !clean_divergence)
+      GRHydro_Reconstruct_drv_cxx <true, true, true, false, false, GRHydro_TrivialReconstruct1d > (cctkGH);
+  else if (do_mhd && do_Ye && do_temp && clean_divergence)
+      GRHydro_Reconstruct_drv_cxx <true, true, true, false, true, GRHydro_TrivialReconstruct1d > (cctkGH);
+  else if (!do_mhd && do_Ye && do_temp)
+      GRHydro_Reconstruct_drv_cxx <false, true, true, false, false, GRHydro_TrivialReconstruct1d > (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<GRHydro_TVDReconstruct1d<0> >::select(do_mhd, do_Ye, do_temp, reconstruct_Wv, clean_divergence, cctkGH);
+       
+    } else if (CCTK_EQUALS(tvd_limiter, "vanleerMC2")) {
+       
+       reconstruct<GRHydro_TVDReconstruct1d<1> >::select(do_mhd, do_Ye, do_temp, reconstruct_Wv, clean_divergence, cctkGH);
+       
+    } else if (CCTK_EQUALS(tvd_limiter, "superbee")) {
+       
+       reconstruct<GRHydro_TVDReconstruct1d<2> >::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<GRHydro_WENOReconstruct1d_cxx<false, true> >::select(do_mhd, do_Ye, do_temp, reconstruct_Wv, clean_divergence, cctkGH);
+          
+    } else {
+       
+       reconstruct<GRHydro_WENOReconstruct1d_cxx<false, false> >::select(do_mhd, do_Ye, do_temp, reconstruct_Wv, clean_divergence, cctkGH);
+    
+    }
+
+  } else if (CCTK_EQUALS(recon_method, "weno-z")) {
+
+       reconstruct<GRHydro_WENOReconstruct1d_cxx<true, false> >::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<GRHydro_MP5Reconstruct1d_cxx<true> >::select(do_mhd, do_Ye, do_temp, reconstruct_Wv, clean_divergence, cctkGH);
+    else
+       reconstruct<GRHydro_MP5Reconstruct1d_cxx<false> >::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 <bool do_mhd,
+          bool do_Ye,
+          bool do_temp,
+          bool do_reconstruct_Wv,
+          bool do_clean_divergence,
+          class RECONSTRUCT>                   // 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 <false, false, false, true, false, RECONSTRUCT> (  \
+  cGH const * const restrict & restrict cctkGH);                               \
+template void                                                                  \
+GRHydro_Reconstruct_drv_cxx <true, false, false, true, false, RECONSTRUCT> (   \
+  cGH const * const restrict & restrict cctkGH);                               \
+template void                                                                  \
+GRHydro_Reconstruct_drv_cxx <true, true, true, true, false, RECONSTRUCT> (     \
+  cGH const * const restrict & restrict cctkGH);                               \
+template void                                                                  \
+GRHydro_Reconstruct_drv_cxx <true, false, false, true, true, RECONSTRUCT> (    \
+  cGH const * const restrict & restrict cctkGH);                               \
+template void                                                                  \
+GRHydro_Reconstruct_drv_cxx <true, true, true, true, true, RECONSTRUCT> (      \
+  cGH const * const restrict & restrict cctkGH);                               \
+template void                                                                  \
+GRHydro_Reconstruct_drv_cxx <false, true, true, true, false, RECONSTRUCT> (    \
+  cGH const * const restrict & restrict cctkGH);                               \
+template void                                                                  \
+GRHydro_Reconstruct_drv_cxx <false, false, false, false, false, RECONSTRUCT> ( \
+  cGH const * const restrict & restrict cctkGH);                               \
+template void                                                                  \
+GRHydro_Reconstruct_drv_cxx <true, false, false, false, false, RECONSTRUCT> (  \
+  cGH const * const restrict & restrict cctkGH);                               \
+template void                                                                  \
+GRHydro_Reconstruct_drv_cxx <true, true, true, false, false, RECONSTRUCT> (    \
+  cGH const * const restrict & restrict cctkGH);                               \
+template void                                                                  \
+GRHydro_Reconstruct_drv_cxx <true, false, false, false, true, RECONSTRUCT> (   \
+  cGH const * const restrict & restrict cctkGH);                               \
+template void                                                                  \
+GRHydro_Reconstruct_drv_cxx <true, true, true, false, true, RECONSTRUCT> (     \
+  cGH const * const restrict & restrict cctkGH);                               \
+template void                                                                  \
+GRHydro_Reconstruct_drv_cxx <false, true, true, false, false, RECONSTRUCT> (   \
+  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 <cmath>
+#include <vector>
+#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 <int dir>
+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 <bool do_mhd,
+          bool do_Ye,
+          bool do_temp,
+          bool do_reconstruct_Wv,
+          bool do_clean_divergence,
+          class RECONSTRUCT>                   // 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<CCTK_REAL> Wvx(0);
+  vector<CCTK_REAL> Wvy(0);
+  vector<CCTK_REAL> 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 T> 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<nx; ++i) {
+          SpaceMask_SetStateBits(space_mask, i+j*nx+k*nx*ny, type_bitsx, not_trivialx);
+        }
+      }
+    }  
+
+  } else if (flux_direction[0]==2) {
+    //Make sure to doublecheck the bounds here!
+    #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 < nx - GRHydro_stencil+1+transport_constraints; ++j) {
+        
+        RECONSTRUCT::template apply<1>(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<ny; ++i) {
+          SpaceMask_SetStateBits(space_mask, j+i*nx+k*nx*ny, type_bitsy, not_trivialy);
+        }
+        
+      }
+    }
+    
+  } else if (flux_direction[0]==3) {
+    //Make sure to doublecheck the bounds here!
+    #pragma omp parallel for
+    for(int k = GRHydro_stencil-1; k < ny - GRHydro_stencil+1+transport_constraints; ++k) {
+      for(int j = GRHydro_stencil-1; j < nx - GRHydro_stencil+1+transport_constraints; ++j) {
+        
+        RECONSTRUCT::template apply<2>(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<nz; ++i) {
+          SpaceMask_SetStateBits(space_mask, j+k*nx+i*nx*ny, type_bitsz, not_trivialz);
+        }
+        
+      }
+    }
+    
+  }
+}
+
+#endif // _GRHYDRO_RECONSTRUCT_DRV_IMPL_H
diff --git a/src/make.code.defn b/src/make.code.defn
index df16d8c..f39f91c 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -28,6 +28,7 @@ SRCS = 	Utils.F90 \
 	GRHydro_PreLoop.F90 \
 	GRHydro_Prim2Con.F90 \
 	GRHydro_Reconstruct.F90 \
+	GRHydro_Reconstruct.cc \
 	GRHydro_ReconstructPoly.F90 \
 	GRHydro_RegisterMask.c \
 	GRHydro_RegisterVars.cc \