From 995467a062526c666e69d6c53d23686da30086cb Mon Sep 17 00:00:00 2001 From: rhaas Date: Tue, 15 Apr 2014 19:49:50 +0000 Subject: GRHydro: make routine to compute speed of sound callable from F90 MHD p2c the code is long and complicated enough that it should better not be copied and pasted. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@633 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45 --- src/GRHydro_Prim2Con.c | 139 +++++++++++++++++++++++++--------------------- src/GRHydro_Prim2ConM.F90 | 12 +++- 2 files changed, 88 insertions(+), 63 deletions(-) diff --git a/src/GRHydro_Prim2Con.c b/src/GRHydro_Prim2Con.c index 7bf8ec4..805584c 100644 --- a/src/GRHydro_Prim2Con.c +++ b/src/GRHydro_Prim2Con.c @@ -25,7 +25,7 @@ static inline __attribute__((always_inline)) void prim2conC(double *w, double *d const double rho, const double vx, const double vy, const double vz, const double eps, const double press); - +static inline void GRHydro_SpeedOfSound(CCTK_ARGUMENTS); @@ -35,6 +35,10 @@ CCTK_FCALL void CCTK_FNAME(Primitive2ConservativeCforF)(cGH ** p_cctkGH) { Primitive2ConservativeC(*p_cctkGH); } +CCTK_FCALL void CCTK_FNAME(GRHydro_SpeedOfSound)(cGH ** p_cctkGH) { + GRHydro_SpeedOfSound(*p_cctkGH); +} + void Primitive2ConservativeC(CCTK_ARGUMENTS) { @@ -61,6 +65,67 @@ void Primitive2ConservativeC(CCTK_ARGUMENTS) g33 = gzz; } + GRHydro_SpeedOfSound(CCTK_PASS_CTOC); + + // if padding is used the the "extra" vector elements could contain junk + // which we do not know how to handle + assert(cctk_lsh[0] == cctk_ash[0]); + assert(cctk_lsh[1] == cctk_ash[1]); + assert(cctk_lsh[2] == cctk_ash[2]); + +#pragma omp parallel for + for(int k = GRHydro_stencil-1; k < cctk_lsh[2]-GRHydro_stencil+1; k++) + for(int j = GRHydro_stencil-1; j < cctk_lsh[1]-GRHydro_stencil+1; j++) +#pragma ivdep // force compiler to vectorize the goddamn loop + for(int i = GRHydro_stencil-1; i < cctk_lsh[0]-GRHydro_stencil+1; i++) { + + const int idx = CCTK_GFINDEX3D(cctkGH,i,j,k); + + const int idxl = CCTK_GFINDEX3D(cctkGH,i-*xoffset,j-*yoffset,k-*zoffset); + const int idxr = CCTK_GFINDEX3D(cctkGH,i+*xoffset,j+*yoffset,k+*zoffset); + + const double g11l = 0.5 * (g11[idx] + g11[idxl]); + const double g12l = 0.5 * (g12[idx] + g12[idxl]); + const double g13l = 0.5 * (g13[idx] + g13[idxl]); + const double g22l = 0.5 * (g22[idx] + g22[idxl]); + const double g23l = 0.5 * (g23[idx] + g23[idxl]); + const double g33l = 0.5 * (g33[idx] + g33[idxl]); + + const double g11r = 0.5 * (g11[idx] + g11[idxr]); + const double g12r = 0.5 * (g12[idx] + g12[idxr]); + const double g13r = 0.5 * (g13[idx] + g13[idxr]); + const double g22r = 0.5 * (g22[idx] + g22[idxr]); + const double g23r = 0.5 * (g23[idx] + g23[idxr]); + const double g33r = 0.5 * (g33[idx] + g33[idxr]); + + const double savg_detl = + sqrt(SpatialDeterminantC(g11l,g12l,g13l,g22l,g23l,g33l)); + const double savg_detr = + sqrt(SpatialDeterminantC(g11r,g12r,g13r,g22r,g23r,g33r)); + + // minus call to p2c + prim2conC(&w_lorentzminus[idx], &densminus[idx], &sxminus[idx], + &syminus[idx], &szminus[idx], &tauminus[idx], + g11l,g12l,g13l,g22l,g23l,g33l, + savg_detl,rhominus[idx], velxminus[idx], velyminus[idx], + velzminus[idx], epsminus[idx], pressminus[idx]); + + + // plus call to p2c + prim2conC(&w_lorentzplus[idx], &densplus[idx], &sxplus[idx], + &syplus[idx], &szplus[idx], &tauplus[idx], + g11r,g12r,g13r,g22r,g23r,g33r, + savg_detr,rhoplus[idx], velxplus[idx], velyplus[idx], + velzplus[idx], epsplus[idx], pressplus[idx]); + } + +} // end function Conservative2PrimitiveC + +static inline void GRHydro_SpeedOfSound(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_PARAMETERS; + DECLARE_CCTK_ARGUMENTS; + // if padding is used the the "extra" vector elements could contain junk // which we do not know how to handle assert(cctk_lsh[0] == cctk_ash[0]); @@ -92,7 +157,7 @@ void Primitive2ConservativeC(CCTK_ARGUMENTS) } else { if(reconstruct_temper) { int n = cctk_ash[0]*cctk_ash[1]*cctk_ash[2]; - int nx = cctk_lsh[0] - (GRHydro_stencil - 1) - (GRHydro_stencil) + 1; + int nx = cctk_lsh[0] - (GRHydro_stencil - 1) - (GRHydro_stencil) + 1 + (transport_constraints && !*xoffset); int *keyerr = malloc(sizeof(*keyerr)*n); int keytemp = 1; // ensure Y_e and temperature within bounds @@ -107,8 +172,8 @@ void Primitive2ConservativeC(CCTK_ARGUMENTS) // call the EOS with slices #pragma omp parallel for - for(int k=GRHydro_stencil-1;k