Commits

commits@lists.einsteintoolkit.org

2 participants
6997 discussions

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 620)
by rhaas＠tapir.caltech.edu 15 Apr '14

15 Apr '14

User: rhaas Date: 2014/04/15 02:49 PM Added: /trunk/src/ GRHydro_HLLE.cc Modified: /trunk/ interface.ccl, schedule.ccl /trunk/src/ GRHydro_RiemannSolve.F90, make.code.defn Log: GRHydro: C++ version of HLLE combined result of work by Christian Ott, Christian Reisswig, Roland Haas File Changes: Directory: /trunk/src/ ====================== File [added]: GRHydro_HLLE.cc Delta lines: +717 -0 =================================================================== --- trunk/src/GRHydro_HLLE.cc (rev 0) +++ trunk/src/GRHydro_HLLE.cc 2014-04-15 19:49:13 UTC (rev 620) @@ -0,0 +1,717 @@ +#include <cmath> +#include <vector> +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#include "SpaceMask.h" + +#define MIN(a,b) (((a)<(b))?(a):(b)) +#define MAX(a,b) (((a)>(b))?(a):(b)) + +using namespace std; + +// some prototypes +extern "C" +CCTK_INT GRHydro_UseGeneralCoordinates(const cGH * cctkGH); + +static inline CCTK_REAL SpatialDeterminantC(CCTK_REAL gxx, CCTK_REAL gxy, + CCTK_REAL gxz, CCTK_REAL gyy, + CCTK_REAL gyz, CCTK_REAL gzz); +static inline void UpperMetricClim(CCTK_REAL * restrict uxx, + CCTK_REAL * restrict uyy, + CCTK_REAL * restrict uzz, + const CCTK_REAL sdet, + const CCTK_REAL gxx, + const CCTK_REAL gxy, + const CCTK_REAL gxz, + const CCTK_REAL gyy, + const CCTK_REAL gyz, + const CCTK_REAL gzz); + +static inline CCTK_REAL max10(const CCTK_REAL * restrict a) +{ + CCTK_REAL maxval=0; + for(int i=0;i<10;++i) maxval=MAX(maxval,a[i]); + return maxval; +} + +static inline CCTK_REAL max9(const CCTK_REAL a, + const CCTK_REAL b, + const CCTK_REAL c, + const CCTK_REAL d, + const CCTK_REAL e, + const CCTK_REAL f, + const CCTK_REAL g, + const CCTK_REAL h, + const CCTK_REAL i) + +{ + return MAX(a,MAX(b,MAX(c,MAX(d,MAX(e,MAX(f,MAX(g,MAX(h,i)))))))); +} + +static inline CCTK_REAL min10(const CCTK_REAL * restrict a) +{ + CCTK_REAL minval=0; + for(int i=0;i<10;++i) minval=MIN(minval,a[i]); + return minval; +} + +static inline void num_x_flux_cc(const CCTK_REAL dens, + const CCTK_REAL sx, + const CCTK_REAL sy, + const CCTK_REAL sz, + const CCTK_REAL tau, + CCTK_REAL *restrict densflux, + CCTK_REAL *restrict sxflux, + CCTK_REAL *restrict syflux, + CCTK_REAL *restrict szflux, + CCTK_REAL *restrict tauflux, + const CCTK_REAL vel, + const CCTK_REAL press, + const CCTK_REAL sdet, + const CCTK_REAL alp, + const CCTK_REAL beta) +{ + + const CCTK_REAL velmbetainvalp = vel - beta / alp; + const CCTK_REAL sqrtdetpress = sdet*press; + + *densflux = dens * velmbetainvalp; + *sxflux = sx * velmbetainvalp + sqrtdetpress; + *syflux = sy * velmbetainvalp; + *szflux = sz * velmbetainvalp; + *tauflux = tau * velmbetainvalp + sqrtdetpress*vel; + + return; +} + + +static inline void eigenvalues_cc(const CCTK_REAL rho, + const CCTK_REAL velx, + const CCTK_REAL vely, + const CCTK_REAL velz, + const CCTK_REAL eps, + const CCTK_REAL press, + const CCTK_REAL cs2, + const CCTK_REAL w, + const CCTK_REAL gxx, + const CCTK_REAL gxy, + const CCTK_REAL gxz, + const CCTK_REAL gyy, + const CCTK_REAL gyz, + const CCTK_REAL gzz, + const CCTK_REAL alp, + const CCTK_REAL beta, + const CCTK_REAL u, + CCTK_REAL *restrict lam) +{ + + const CCTK_REAL vlowx = gxx*velx + gxy*vely + gxz*velz; + const CCTK_REAL vlowy = gxy*velx + gyy*vely + gyz*velz; + const CCTK_REAL vlowz = gxz*velx + gyz*vely + gzz*velz; + const CCTK_REAL v2 = vlowx*velx + vlowy*vely + vlowz*velz; + + const CCTK_REAL boa = beta/alp; + lam[1] = velx - boa; + lam[2] = lam[1]; + lam[3] = lam[1]; + const CCTK_REAL lam_tmp1 = 1.0e0/(1.0-v2*cs2); + const CCTK_REAL lam_tmp2 = sqrt(cs2*(1.0-v2)* + (u*(1.0-v2*cs2) - velx*velx*(1.0-cs2))); + const CCTK_REAL lam_tmp3 = velx*(1.0-cs2); + lam[0] = (lam_tmp3 - lam_tmp2)*lam_tmp1 - boa; + lam[4] = (lam_tmp3 + lam_tmp2)*lam_tmp1 - boa; + + return; +} + +// helpers for H viscosity +static inline CCTK_REAL etaX(const cGH *cctkGH, + const int i, const int j, const int k, + const CCTK_REAL *restrict vel, + const CCTK_REAL *restrict cs) { + + const int vidxr = CCTK_VECTGFINDEX3D(cctkGH,i+1,j,k,0); + const int vidx = CCTK_VECTGFINDEX3D(cctkGH,i,j,k,0); + const int idxr = CCTK_GFINDEX3D(cctkGH,i+1,j,k); + const int idx = CCTK_GFINDEX3D(cctkGH,i,j,k); + + return 0.5 * (fabs(vel[idxr]-vel[idx]) + fabs(cs[idxr]-cs[idx])); +} + +static inline CCTK_REAL etaY(const cGH *cctkGH, + const int i, const int j, const int k, + const CCTK_REAL *restrict vel, + const CCTK_REAL *restrict cs) { + + const int vidxr = CCTK_VECTGFINDEX3D(cctkGH,i,j+1,k,1); + const int vidx = CCTK_VECTGFINDEX3D(cctkGH,i,j,k,1); + const int idxr = CCTK_GFINDEX3D(cctkGH,i,j+1,k); + const int idx = CCTK_GFINDEX3D(cctkGH,i,j,k); + + return 0.5 * (fabs(vel[idxr]-vel[idx]) + fabs(cs[idxr]-cs[idx])); +} + +static inline CCTK_REAL etaZ(const cGH *cctkGH, + const int i, const int j, const int k, + const CCTK_REAL *restrict vel, + const CCTK_REAL *restrict cs) { + + int vidxr = CCTK_VECTGFINDEX3D(cctkGH,i,j,k+1,2); + int vidx = CCTK_VECTGFINDEX3D(cctkGH,i,j,k,2); + int idxr = CCTK_GFINDEX3D(cctkGH,i,j,k+1); + int idx = CCTK_GFINDEX3D(cctkGH,i,j,k); + + return 0.5 * (fabs(vel[vidxr]-vel[vidx]) + fabs(cs[idxr]-cs[idx])); +} + + +template<const int fdir, const bool do_ye, const bool do_hvisc, const bool do_tracers> +void GRHydro_HLLE_CC_LL(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + // save memory when multipatch is not used + CCTK_REAL * restrict g11, * restrict g12, * restrict g13; + CCTK_REAL * restrict g22, * restrict g23, * restrict g33; + CCTK_REAL * restrict beta1, * restrict beta2, * restrict beta3; + CCTK_REAL * restrict vup; + const int offsetx = fdir==1, offsety = fdir==2, offsetz = fdir==3; + assert(offsetx+offsety+offsetz == 1); + + if(GRHydro_UseGeneralCoordinates(cctkGH)) { + g11 = gaa; + g12 = gab; + g13 = gac; + g22 = gbb; + g23 = gbc; + g33 = gcc; + beta1 = betaa; + beta2 = betab; + beta3 = betac; + vup = lvel; + } else { + g11 = gxx; + g12 = gxy; + g13 = gxz; + g22 = gyy; + g23 = gyz; + g33 = gzz; + beta1 = betax; + beta2 = betay; + beta3 = betaz; + vup = vel; + } + + int type_bits, trivial; + if(fdir==1) { + type_bits = SpaceMask_GetTypeBits("Hydro_RiemannProblemX"); + trivial = SpaceMask_GetStateBits("Hydro_RiemannProblemX", + "trivial"); + } else if (fdir==2) { + type_bits = SpaceMask_GetTypeBits("Hydro_RiemannProblemY"); + trivial = SpaceMask_GetStateBits("Hydro_RiemannProblemY", + "trivial"); + } else if (fdir==3) { + type_bits = SpaceMask_GetTypeBits("Hydro_RiemannProblemZ"); + trivial = SpaceMask_GetStateBits("Hydro_RiemannProblemZ", + "trivial"); + } else { + CCTK_VError(__LINE__, __FILE__, CCTK_THORNSTRING, + "Flux direction %d not 1,2,3.", fdir); + } + +#pragma omp parallel for + for(int k = GRHydro_stencil-1; k < cctk_lsh[2]-GRHydro_stencil; k++) + for(int j = GRHydro_stencil-1; j < cctk_lsh[1]-GRHydro_stencil; j++) + for(int i = GRHydro_stencil-1; i < cctk_lsh[0]-GRHydro_stencil; i++) { + + const int idx = CCTK_GFINDEX3D(cctkGH,i,j,k); + const int idxr = + CCTK_GFINDEX3D(cctkGH,i+offsetx,j+offsety,k+offsetz); + + CCTK_REAL consp[5],consm_i[5]; + CCTK_REAL fplus[5],fminus[5], f1[5]; + CCTK_REAL lamplusminus[10]; + CCTK_REAL * const lamminus = lamplusminus; + CCTK_REAL * const lamplus = lamplusminus+5; + CCTK_REAL qdiff[5]; + CCTK_REAL avg_beta; + + consp[0] = densplus[idx]; + consp[1] = sxplus[idx]; + consp[2] = syplus[idx]; + consp[3] = szplus[idx]; + consp[4] = tauplus[idx]; + + consm_i[0] = densminus[idxr]; + consm_i[1] = sxminus[idxr]; + consm_i[2] = syminus[idxr]; + consm_i[3] = szminus[idxr]; + consm_i[4] = tauminus[idxr]; + + if(fdir==1) + avg_beta = 0.5 * (beta1[idxr]+beta1[idx]); + else if(fdir==2) + avg_beta = 0.5 * (beta2[idxr]+beta2[idx]); + else + avg_beta = 0.5 * (beta3[idxr]+beta3[idx]); + + const CCTK_REAL avg_alp = 0.5 * (alp[idxr]+alp[idx]); + + const CCTK_REAL gxxh = 0.5 * (g11[idx] + g11[idxr]); + const CCTK_REAL gxyh = 0.5 * (g12[idx] + g12[idxr]); + const CCTK_REAL gxzh = 0.5 * (g13[idx] + g13[idxr]); + const CCTK_REAL gyyh = 0.5 * (g22[idx] + g22[idxr]); + const CCTK_REAL gyzh = 0.5 * (g23[idx] + g23[idxr]); + const CCTK_REAL gzzh = 0.5 * (g33[idx] + g33[idxr]); + + const CCTK_REAL savg_detr = + sqrt(SpatialDeterminantC(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh)); + + + // If the Riemann problem is trivial, just calculate the fluxes + // from the left state and skip to the next cell + if ( SpaceMask_CheckStateBits(space_mask, idx, type_bits, trivial) ) { + + if(fdir==1) + num_x_flux_cc(consp[0],consp[1],consp[2],consp[3],consp[4], + &f1[0],&f1[1],&f1[2],&f1[3],&f1[4], + velxplus[idx],pressplus[idx],savg_detr, + avg_alp,avg_beta); + else if (fdir==2) + num_x_flux_cc(consp[0],consp[2],consp[3],consp[1],consp[4], + &f1[0],&f1[2],&f1[3],&f1[1],&f1[4], + velyplus[idx],pressplus[idx],savg_detr, + avg_alp,avg_beta); + else + num_x_flux_cc(consp[0],consp[3],consp[1],consp[2],consp[4], + &f1[0],&f1[3],&f1[1],&f1[2],&f1[4], + velzplus[idx],pressplus[idx],savg_detr, + avg_alp,avg_beta); + + } else { + + // normal case -- full HLLE Riemann solution + CCTK_REAL uxxh,uyyh,uzzh; + UpperMetricClim(&uxxh,&uyyh,&uzzh,savg_detr,gxxh,gxyh,gxzh, + gyyh, gyzh, gzzh); + + CCTK_REAL usendh; + if(fdir==1) + usendh = uxxh; + else if(fdir==2) + usendh = uyyh; + else + usendh = uzzh; + + // calculate the jumps in the conserved variables + qdiff[0] = consm_i[0] - consp[0]; + qdiff[1] = consm_i[1] - consp[1]; + qdiff[2] = consm_i[2] - consp[2]; + qdiff[3] = consm_i[3] - consp[3]; + qdiff[4] = consm_i[4] - consp[4]; + + if (fdir==1) { + eigenvalues_cc(rhominus[idxr], + velxminus[idxr], + velyminus[idxr], + velzminus[idxr], + epsminus[idxr], + pressminus[idxr], + cs2minus[idxr], + w_lorentzminus[idxr], + gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, + avg_alp,avg_beta,usendh, + lamminus); + eigenvalues_cc(rhoplus[idx], + velxplus[idx], + velyplus[idx], + velzplus[idx], + epsplus[idx], + pressplus[idx], + cs2plus[idx], + w_lorentzplus[idx], + gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, + avg_alp,avg_beta,usendh, + lamplus); + num_x_flux_cc(consp[0],consp[1],consp[2],consp[3],consp[4], + &fplus[0],&fplus[1],&fplus[2],&fplus[3],&fplus[4], + velxplus[idx],pressplus[idx],savg_detr, + avg_alp,avg_beta); + num_x_flux_cc(consm_i[0],consm_i[1],consm_i[2],consm_i[3],consm_i[4], + &fminus[0],&fminus[1],&fminus[2],&fminus[3],&fminus[4], + velxminus[idxr],pressminus[idxr],savg_detr, + avg_alp,avg_beta); + } + else if(fdir==2) { + eigenvalues_cc(rhominus[idxr], + velyminus[idxr], + velzminus[idxr], + velxminus[idxr], + epsminus[idxr], + pressminus[idxr], + cs2minus[idxr], + w_lorentzminus[idxr], + gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, + avg_alp,avg_beta,usendh, + lamminus); + eigenvalues_cc(rhoplus[idx], + velyplus[idx], + velzplus[idx], + velxplus[idx], + epsplus[idx], + pressplus[idx], + cs2plus[idx], + w_lorentzplus[idx], + gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, + avg_alp,avg_beta,usendh, + lamplus); + num_x_flux_cc(consp[0],consp[2],consp[3],consp[1],consp[4], + &fplus[0],&fplus[2],&fplus[3],&fplus[1],&fplus[4], + velyplus[idx],pressplus[idx],savg_detr, + avg_alp,avg_beta); + num_x_flux_cc(consm_i[0],consm_i[2],consm_i[3],consm_i[1],consm_i[4], + &fminus[0],&fminus[2],&fminus[3],&fminus[1],&fminus[4], + velyminus[idxr],pressminus[idxr],savg_detr, + avg_alp,avg_beta); + } else { + eigenvalues_cc(rhominus[idxr], + velzminus[idxr], + velxminus[idxr], + velyminus[idxr], + epsminus[idxr], + pressminus[idxr], + cs2minus[idxr], + w_lorentzminus[idxr], + gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, + avg_alp,avg_beta,usendh, + lamminus); + eigenvalues_cc(rhoplus[idx], + velzplus[idx], + velxplus[idx], + velyplus[idx], + epsplus[idx], + pressplus[idx], + cs2plus[idx], + w_lorentzplus[idx], + gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, + avg_alp,avg_beta,usendh, + lamplus); + num_x_flux_cc(consp[0],consp[3],consp[1],consp[2],consp[4], + &fplus[0],&fplus[3],&fplus[1],&fplus[2],&fplus[4], + velzplus[idx],pressplus[idx],savg_detr, + avg_alp,avg_beta); + num_x_flux_cc(consm_i[0],consm_i[3],consm_i[1],consm_i[2],consm_i[4], + &fminus[0],&fminus[3],&fminus[1],&fminus[2],&fminus[4], + velzminus[idxr],pressminus[idxr],savg_detr, + avg_alp,avg_beta); + } + + // H-viscosity + if(do_hvisc) { + CCTK_REAL etabar = 0.0; + + if(fdir==1) { + etabar = max9(etaX(cctkGH,i,j,k,vup,eos_c), + etaY(cctkGH,i,j,k,vup,eos_c), + etaY(cctkGH,i+1,j,k,vup,eos_c), + etaY(cctkGH,i,j-1,k,vup,eos_c), + etaY(cctkGH,i+1,j-1,k,vup,eos_c), + etaZ(cctkGH,i,j,k,vup,eos_c), + etaZ(cctkGH,i+1,j,k,vup,eos_c), + etaZ(cctkGH,i,j,k-1,vup,eos_c), + etaZ(cctkGH,i+1,j,k-1,vup,eos_c)); + } else if(fdir==2) { + etabar = max9(etaY(cctkGH,i,j,k,vup,eos_c), + etaX(cctkGH,i,j,k,vup,eos_c), + etaX(cctkGH,i,j+1,k,vup,eos_c), + etaX(cctkGH,i-1,j,k,vup,eos_c), + etaX(cctkGH,i-1,j+1,k,vup,eos_c), + etaZ(cctkGH,i,j,k,vup,eos_c), + etaZ(cctkGH,i,j+1,k,vup,eos_c), + etaZ(cctkGH,i,j,k-1,vup,eos_c), + etaZ(cctkGH,i,j+1,k-1,vup,eos_c)); + } else if(fdir==3) { + etabar = max9(etaZ(cctkGH,i,j,k,vup,eos_c), + etaX(cctkGH,i,j,k,vup,eos_c), + etaX(cctkGH,i,j,k+1,vup,eos_c), + etaX(cctkGH,i-1,j,k,vup,eos_c), + etaX(cctkGH,i-1,j,k+1,vup,eos_c), + etaY(cctkGH,i,j,k,vup,eos_c), + etaY(cctkGH,i,j,k+1,vup,eos_c), + etaY(cctkGH,i,j-1,k,vup,eos_c), + etaY(cctkGH,i,j-1,k+1,vup,eos_c)); + } else { + CCTK_ERROR("Flux direction not x,y,z"); + } + // modify eigenvalues of Roe's matrix by computed H viscosity + for(int m=0;m<10;m++) { + lamplusminus[m] = copysign(1.0,lamplusminus[m])*MAX(fabs(lamplusminus[m]),etabar); + // copysign is C99 and the equivalent of FORTRAN sign(a,b) + } + + } + + // Find minimum and maximum wavespeeds + CCTK_REAL charmax = max10(lamplusminus); + CCTK_REAL charmin = min10(lamplusminus); + + CCTK_REAL icharpm = 1.0 / (charmax-charmin); + + for(int m=0;m<5;m++) { + f1[m] = (charmax * fplus[m] - charmin * fminus[m] + + charmax*charmin * qdiff[m]) * icharpm; + } + + } // else trivial + + densflux[idx] = f1[0]; + sxflux[idx] = f1[1]; + syflux[idx] = f1[2]; + szflux[idx] = f1[3]; + tauflux[idx] = f1[4]; + + if(do_ye) { + if(densflux[idx] > 0.0) + Y_e_con_flux[idx] = Y_e_plus[idx] * densflux[idx]; + else + Y_e_con_flux[idx] = Y_e_minus[idxr] * densflux[idx]; + } + + if(do_tracers) { + if(densflux[idx] > 0.0) + for(int m=0;m<number_of_tracers;m++) { + const int idx4 = CCTK_VECTGFINDEX3D(cctkGH,i,j,k,m); + cons_tracerflux[idx4] = cons_tracerplus[idx4] + * densflux[idx]; + } + else + for(int m=0;m<number_of_tracers;m++) { + const int idx4 = CCTK_VECTGFINDEX3D(cctkGH,i,j,k,m); + const int idx4r = + CCTK_VECTGFINDEX3D(cctkGH,i+offsetx,j+offsety,k+offsetz,m); + cons_tracerflux[idx4] = cons_tracerminus[idx4r] + * densflux[idx]; + } + } // do tracers + + + } // end big loop +} + +// instantiate the templates explicitely +template void GRHydro_HLLE_CC_LL<1,false,false,false>(CCTK_ARGUMENTS); +template void GRHydro_HLLE_CC_LL<2,false,false,false>(CCTK_ARGUMENTS); +template void GRHydro_HLLE_CC_LL<3,false,false,false>(CCTK_ARGUMENTS); +template void GRHydro_HLLE_CC_LL<1,true, false,false>(CCTK_ARGUMENTS); +template void GRHydro_HLLE_CC_LL<2,true, false,false>(CCTK_ARGUMENTS); +template void GRHydro_HLLE_CC_LL<3,true, false,false>(CCTK_ARGUMENTS); +template void GRHydro_HLLE_CC_LL<1,false,true,false>(CCTK_ARGUMENTS); +template void GRHydro_HLLE_CC_LL<2,false,true,false>(CCTK_ARGUMENTS); +template void GRHydro_HLLE_CC_LL<3,false,true,false>(CCTK_ARGUMENTS); +template void GRHydro_HLLE_CC_LL<1,true, true,false>(CCTK_ARGUMENTS); +template void GRHydro_HLLE_CC_LL<2,true, true,false>(CCTK_ARGUMENTS); +template void GRHydro_HLLE_CC_LL<3,true, true,false>(CCTK_ARGUMENTS); + +template void GRHydro_HLLE_CC_LL<1,false,false,true>(CCTK_ARGUMENTS); +template void GRHydro_HLLE_CC_LL<2,false,false,true>(CCTK_ARGUMENTS); +template void GRHydro_HLLE_CC_LL<3,false,false,true>(CCTK_ARGUMENTS); +template void GRHydro_HLLE_CC_LL<1,true, false,true>(CCTK_ARGUMENTS); +template void GRHydro_HLLE_CC_LL<2,true, false,true>(CCTK_ARGUMENTS); +template void GRHydro_HLLE_CC_LL<3,true, false,true>(CCTK_ARGUMENTS); +template void GRHydro_HLLE_CC_LL<1,false,true,true>(CCTK_ARGUMENTS); +template void GRHydro_HLLE_CC_LL<2,false,true,true>(CCTK_ARGUMENTS); +template void GRHydro_HLLE_CC_LL<3,false,true,true>(CCTK_ARGUMENTS); +template void GRHydro_HLLE_CC_LL<1,true, true,true>(CCTK_ARGUMENTS); +template void GRHydro_HLLE_CC_LL<2,true, true,true>(CCTK_ARGUMENTS); +template void GRHydro_HLLE_CC_LL<3,true, true,true>(CCTK_ARGUMENTS); + + + +void GRHydro_HLLE_CC(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + // flux direction 1 + if(*flux_direction==1) + if(apply_H_viscosity) + if(*evolve_Y_e) + if(evolve_tracer) + GRHydro_HLLE_CC_LL<1,true,true,true>(CCTK_PASS_CTOC); + else + GRHydro_HLLE_CC_LL<1,true,true,false>(CCTK_PASS_CTOC); + else + if(evolve_tracer) + GRHydro_HLLE_CC_LL<1,false,true,true>(CCTK_PASS_CTOC); + else + GRHydro_HLLE_CC_LL<1,false,true,false>(CCTK_PASS_CTOC); + else + if(*evolve_Y_e) + if(evolve_tracer) + GRHydro_HLLE_CC_LL<1,true,false,true>(CCTK_PASS_CTOC); + else + GRHydro_HLLE_CC_LL<1,true,false,false>(CCTK_PASS_CTOC); + else + if(evolve_tracer) + GRHydro_HLLE_CC_LL<1,false,false,true>(CCTK_PASS_CTOC); + else + GRHydro_HLLE_CC_LL<1,false,false,false>(CCTK_PASS_CTOC); + + // flux direction 2 + else if(*flux_direction==2) + if(apply_H_viscosity) + if(*evolve_Y_e) + if(evolve_tracer) + GRHydro_HLLE_CC_LL<2,true,true,true>(CCTK_PASS_CTOC); + else + GRHydro_HLLE_CC_LL<2,true,true,false>(CCTK_PASS_CTOC); + else + if(evolve_tracer) + GRHydro_HLLE_CC_LL<2,false,true,true>(CCTK_PASS_CTOC); + else + GRHydro_HLLE_CC_LL<2,false,true,false>(CCTK_PASS_CTOC); + else + if(*evolve_Y_e) + if(evolve_tracer) + GRHydro_HLLE_CC_LL<2,true,false,true>(CCTK_PASS_CTOC); + else + GRHydro_HLLE_CC_LL<2,true,false,false>(CCTK_PASS_CTOC); + else + if(evolve_tracer) + GRHydro_HLLE_CC_LL<2,false,false,true>(CCTK_PASS_CTOC); + else + GRHydro_HLLE_CC_LL<2,false,false,false>(CCTK_PASS_CTOC); + + // flux direction 3 + else if(*flux_direction==3) + if(apply_H_viscosity) + if(*evolve_Y_e) + if(evolve_tracer) + GRHydro_HLLE_CC_LL<3,true,true,true>(CCTK_PASS_CTOC); + else + GRHydro_HLLE_CC_LL<3,true,true,false>(CCTK_PASS_CTOC); + else + if(evolve_tracer) + GRHydro_HLLE_CC_LL<3,false,true,true>(CCTK_PASS_CTOC); + else + GRHydro_HLLE_CC_LL<3,false,true,false>(CCTK_PASS_CTOC); + else + if(*evolve_Y_e) + if(evolve_tracer) + GRHydro_HLLE_CC_LL<3,true,false,true>(CCTK_PASS_CTOC); + else + GRHydro_HLLE_CC_LL<3,true,false,false>(CCTK_PASS_CTOC); + else + if(evolve_tracer) + GRHydro_HLLE_CC_LL<3,false,false,true>(CCTK_PASS_CTOC); + else + GRHydro_HLLE_CC_LL<3,false,false,false>(CCTK_PASS_CTOC); + + else + CCTK_ERROR("Illegal flux direction!"); + + return; +} + +extern "C" +CCTK_FCALL void CCTK_FNAME(GRHydro_HLLE_CC_F2C)(cGH ** p_cctkGH) { + GRHydro_HLLE_CC(*p_cctkGH); +} + + +extern "C" +void H_viscosity_calc_cs_cc(CCTK_ARGUMENTS) { + + // This is a preperatory routine for H viscosity. + // All it does is fill a speed of sound helper array + + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + size_t n = size_t(cctk_ash[0]*cctk_ash[1]*cctk_ash[2]); + std::vector<int> keyerr(n); + int anyerr = 0; + int keytemp = 0; + + if(!*evolve_temper) { +#pragma omp parallel for + for(int k=0;k<cctk_lsh[2];k++) + for(int j=0;j<cctk_lsh[1];j++) { + int i = CCTK_GFINDEX3D(cctkGH,0,j,k); + CCTK_REAL cs2[cctk_ash[0]]; + EOS_Omni_cs2(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,cctk_lsh[0], + &(rho[i]),&(eps[i]),NULL,NULL,&cs2[i], + &(keyerr[i]),&anyerr); + + for(int ii=0;ii<cctk_lsh[0];i++) { + eos_c[ii+i] = sqrt(cs2[ii]); + } + } + } else { +#pragma omp parallel for + for(int k=0;k<cctk_lsh[2];k++) + for(int j=0;j<cctk_lsh[1];j++) { + int i = CCTK_GFINDEX3D(cctkGH,0,j,k); + CCTK_REAL cs2[cctk_ash[0]]; + EOS_Omni_cs2(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,cctk_lsh[0], + &rho[i],&eps[i],&temperature[i],&Y_e[i], + &cs2[i],&keyerr[i],&anyerr); + for(int ii=0;ii<cctk_lsh[0];i++) { + eos_c[ii+i] = sqrt(cs2[ii]); + } + } // for + if(anyerr) { + for(int i=0;i<n;i++) { + if(keyerr[i] != 0) { +#pragma omp critical + { + CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, + "rl: %d i,x,y,z: %d %15.6E %15.6E %15.6E, keyerr: %d", + *GRHydro_reflevel, i, x[i], y[i], z[i], keyerr[i]); + CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, + "rl: %d r,t,ye: %15.6E %15.6E %15.6E, keyerr: %d", + *GRHydro_reflevel, rho[i], temperature[i], Y_e[i], keyerr[i]); + } + } + } // for i, i<n + CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, + "Aborting!"); + } // anyerr + } // if evolve_temper + +} + +static inline void UpperMetricClim(CCTK_REAL * restrict uxx, + CCTK_REAL * restrict uyy, + CCTK_REAL * restrict uzz, + const CCTK_REAL sdet, + const CCTK_REAL gxx, + const CCTK_REAL gxy, + const CCTK_REAL gxz, + const CCTK_REAL gyy, + const CCTK_REAL gyz, + const CCTK_REAL gzz) +{ + const CCTK_REAL invdet = 1.0 / (sdet*sdet); + *uxx = (-gyz*gyz + gyy*gzz)*invdet; + // *uxy = (gxz*gyz - gxy*gzz)*invdet; + // *uxz = (-gxz*gyy + gxy*gyz)*invdet; + *uyy = (-gxz*gxz + gxx*gzz)*invdet; + // *uyz = (gxy*gxz - gxx*gyz)*invdet; + *uzz = (-gxy*gxy + gxx*gyy)*invdet; + + return; +} + +static inline CCTK_REAL SpatialDeterminantC(CCTK_REAL gxx, CCTK_REAL gxy, + CCTK_REAL gxz, CCTK_REAL gyy, + CCTK_REAL gyz, CCTK_REAL gzz) +{ + return -gxz*gxz*gyy + 2.0*gxy*gxz*gyz - gxx*gyz*gyz + - gxy*gxy*gzz + gxx*gyy*gzz; +} + File [modified]: GRHydro_RiemannSolve.F90 Delta lines: +5 -1 =================================================================== --- trunk/src/GRHydro_RiemannSolve.F90 2014-04-15 19:49:10 UTC (rev 619) +++ trunk/src/GRHydro_RiemannSolve.F90 2014-04-15 19:49:13 UTC (rev 620) @@ -38,7 +38,11 @@ if (CCTK_EQUALS(riemann_solver,"HLLE")) then - call GRHydro_HLLE(CCTK_PASS_FTOF) + if (use_cxx_code.eq.0) then + call GRHydro_HLLE(CCTK_PASS_FTOF) + else + call GRHydro_HLLE_CC_F2C(cctkGH) + endif if (evolve_tracer .ne. 0) then File [modified]: make.code.defn Delta lines: +2 -1 =================================================================== --- trunk/src/make.code.defn 2014-04-15 19:49:10 UTC (rev 619) +++ trunk/src/make.code.defn 2014-04-15 19:49:13 UTC (rev 620) @@ -86,7 +86,8 @@ GRHydro_PPM.cc \ GRHydro_ePPM.cc \ GRHydro_PPMReconstruct_drv_opt.cc \ - GRHydro_Wrappers.F90 + GRHydro_Wrappers.F90 \ + GRHydro_HLLE.cc # Subdirectories containing source files SUBDIRS = Directory: /trunk/ ================== File [modified]: interface.ccl Delta lines: +1 -1 =================================================================== --- trunk/interface.ccl 2014-04-15 19:49:10 UTC (rev 619) +++ trunk/interface.ccl 2014-04-15 19:49:13 UTC (rev 620) @@ -295,7 +295,7 @@ CCTK_REAL INOUT ARRAY temp, \ CCTK_REAL IN ARRAY ye, \ CCTK_REAL OUT ARRAY press, \ - CCTK_REAL OUT ARRAY cs2, \ + CCTK_REAL OUT ARRAY cs2, \ CCTK_INT OUT ARRAY keyerr, \ CCTK_INT OUT anyerr) File [modified]: schedule.ccl Delta lines: +13 -5 =================================================================== --- trunk/schedule.ccl 2014-04-15 19:49:10 UTC (rev 619) +++ trunk/schedule.ccl 2014-04-15 19:49:13 UTC (rev 620) @@ -1616,10 +1616,18 @@ if (apply_H_viscosity) { - SCHEDULE H_viscosity in GRHydroRHS BEFORE FluxTerms - { - LANG: FORTRAN - OPTION: LOCAL - } "Compute local temporaries for H viscosity" + if (use_cxx_code) { + SCHEDULE H_viscosity_calc_cs_cc in GRHydroRHS BEFORE FluxTerms + { + LANG: C + OPTION: LOCAL + } "Compute local temporaries for H viscosity - C++ version" + } else { + SCHEDULE H_viscosity in GRHydroRHS BEFORE FluxTerms + { + LANG: FORTRAN + OPTION: LOCAL + } "Compute local temporaries for H viscosity" + } }

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 619)
by rhaas＠tapir.caltech.edu 15 Apr '14

15 Apr '14

User: rhaas Date: 2014/04/15 02:49 PM Modified: /trunk/src/ GRHydro_PPMReconstruct_drv_opt.cc Log: GRHydro: reconstruct_Wv for ePPM/oPPM. File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_PPMReconstruct_drv_opt.cc Delta lines: +73 -11 =================================================================== --- trunk/src/GRHydro_PPMReconstruct_drv_opt.cc 2014-04-15 19:49:07 UTC (rev 618) +++ trunk/src/GRHydro_PPMReconstruct_drv_opt.cc 2014-04-15 19:49:10 UTC (rev 619) @@ -10,7 +10,12 @@ #include "SpaceMask.h" #include "GRHydro_PPM_opt.h" +#include "GRHydro_Reconstruct_drv_cxx.hh" +#undef velx +#undef vely +#undef velz + using namespace std; extern "C" CCTK_INT GRHydro_UseGeneralCoordinates(const cGH * cctkGH); @@ -49,12 +54,25 @@ //Multipatch related pointers 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; + if(GRHydro_UseGeneralCoordinates(cctkGH)) { vup=lvel; Bprim=lBvec; + g11=gaa; g12=gab; g13=gac; + g22=gbb; g23=gbc; + g33=gcc; } else { vup=vel; Bprim=Bvec; + g11=gxx; g12=gxy; g13=gxz; + g22=gyy; g23=gyz; + g33=gzz; } const int nx=cctk_ash[0]; @@ -217,6 +235,18 @@ //index = i + j*nx + k*nx*ny; + vector<CCTK_REAL> Wvup(0); + // allocate temp memory for Wv reconstruction + if (reconstruct_Wv) { + Wvup.resize(3*nxyz, 0); + #pragma omp parallel for + for (int i=0; i < nxyz; ++i) { + Wvup[i] = w_lorentz[i]*vup[i]; + Wvup[i+nxyz] = w_lorentz[i]*vup[i+nxyz]; + Wvup[i+2*nxyz] = w_lorentz[i]*vup[i+2*nxyz]; + } + } + // PPM starts: if (*flux_direction == 1) { #pragma omp parallel @@ -257,10 +287,16 @@ oned_slice(nx,nstart,1,rho,rho1d); - oned_slice(nx,nstart, 1,vup,velx1d); - oned_slice(nx,nstart+nxyz, 1,vup,vely1d); - oned_slice(nx,nstart+2*nxyz,1,vup,velz1d); - + // TODO: reconstruct_Wv is not yet templated! + if (!reconstruct_Wv) { + oned_slice(nx,nstart, 1,vup,velx1d); + oned_slice(nx,nstart+nxyz, 1,vup,vely1d); + oned_slice(nx,nstart+2*nxyz,1,vup,velz1d); + } else { + oned_slice(nx,nstart, 1,&Wvup.front(),velx1d); + oned_slice(nx,nstart+nxyz, 1,&Wvup.front(),vely1d); + oned_slice(nx,nstart+2*nxyz,1,&Wvup.front(),velz1d); + } oned_slice(nx,nstart,1,eps,eps1d); oned_slice(nx,nstart,1,press,press1d); @@ -300,6 +336,11 @@ oned_unslice(nx,nstart,1,velx1d_plus,velxplus); oned_unslice(nx,nstart,1,vely1d_plus,velyplus); oned_unslice(nx,nstart,1,velz1d_plus,velzplus); + if (reconstruct_Wv) + undo_Wv<0>(nx, velxminus, velyminus, velzminus, + velxplus, velyplus, velzplus, + g11, g12, g13, g22, g23, g33, + cctkGH, j, k); oned_unslice(nx,nstart,1,eps1d_minus,epsminus); oned_unslice(nx,nstart,1,eps1d_plus,epsplus); @@ -368,10 +409,15 @@ oned_slice(ny,nstart,nx,rho,rho1d); - oned_slice(ny,nstart, nx,vup,velx1d); - oned_slice(ny,nstart+nxyz, nx,vup,vely1d); - oned_slice(ny,nstart+2*nxyz,nx,vup,velz1d); - + if (!reconstruct_Wv) { + oned_slice(ny,nstart, nx,vup,velx1d); + oned_slice(ny,nstart+nxyz, nx,vup,vely1d); + oned_slice(ny,nstart+2*nxyz,nx,vup,velz1d); + } else { + oned_slice(ny,nstart, nx,&Wvup.front(),velx1d); + oned_slice(ny,nstart+nxyz, nx,&Wvup.front(),vely1d); + oned_slice(ny,nstart+2*nxyz,nx,&Wvup.front(),velz1d); + } oned_slice(ny,nstart,nx,eps,eps1d); oned_slice(ny,nstart,nx,press,press1d); @@ -412,6 +458,11 @@ oned_unslice(ny,nstart,nx,velx1d_plus,velxplus); oned_unslice(ny,nstart,nx,vely1d_plus,velyplus); oned_unslice(ny,nstart,nx,velz1d_plus,velzplus); + if (reconstruct_Wv) + undo_Wv<1>(ny, velyminus, velzminus, velxminus, + velyplus, velzplus, velxplus, + g22, g23, g12, g33, g13, g11, + cctkGH, j, k); oned_unslice(ny,nstart,nx,eps1d_minus,epsminus); oned_unslice(ny,nstart,nx,eps1d_plus,epsplus); @@ -481,9 +532,15 @@ oned_slice(nz,nstart,nxy,rho,rho1d); - oned_slice(nz,nstart, nxy,vup,velx1d); - oned_slice(nz,nstart+nxyz, nxy,vup,vely1d); - oned_slice(nz,nstart+2*nxyz,nxy,vup,velz1d); + if (!reconstruct_Wv) { + oned_slice(nz,nstart, nxy,vup,velx1d); + oned_slice(nz,nstart+nxyz, nxy,vup,vely1d); + oned_slice(nz,nstart+2*nxyz,nxy,vup,velz1d); + } else { + oned_slice(nz,nstart, nxy,&Wvup.front(),velx1d); + oned_slice(nz,nstart+nxyz, nxy,&Wvup.front(),vely1d); + oned_slice(nz,nstart+2*nxyz,nxy,&Wvup.front(),velz1d); + } oned_slice(nz,nstart,nxy,eps,eps1d); oned_slice(nz,nstart,nxy,press,press1d); @@ -525,6 +582,11 @@ oned_unslice(nz,nstart,nxy,velx1d_plus,velxplus); oned_unslice(nz,nstart,nxy,vely1d_plus,velyplus); oned_unslice(nz,nstart,nxy,velz1d_plus,velzplus); + if (reconstruct_Wv) + undo_Wv<2>(nz, velzminus, velxminus, velyminus, + velzplus, velxplus, velyplus, + g33, g13, g23, g11, g12, g22, + cctkGH, j, k); oned_unslice(nz,nstart,nxy,eps1d_minus,epsminus); oned_unslice(nz,nstart,nxy,eps1d_plus,epsplus);

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 618)
by rhaas＠tapir.caltech.edu 15 Apr '14

15 Apr '14

User: rhaas Date: 2014/04/15 02:49 PM Modified: /trunk/src/ GRHydro_PPM.cc, GRHydro_PPMReconstruct_drv_opt.cc Log: GRHydro: C++ PPM: enable OMP, bugfixes * Cleanup of code. * Const-correctness. * correct upper bound for PPM. File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_PPM.cc Delta lines: +52 -55 =================================================================== --- trunk/src/GRHydro_PPM.cc 2014-04-15 19:49:04 UTC (rev 617) +++ trunk/src/GRHydro_PPM.cc 2014-04-15 19:49:07 UTC (rev 618) @@ -23,9 +23,9 @@ */ -static inline void steep(const double *restrict x, - const double *restrict dx, - double *restrict dmx, +static inline void steep(const double * const restrict x, + const double * const restrict dx, + double * const restrict dmx, const int i) { if ( (x[i+1] - x[i]) * (x[i]-x[i-1]) > 0.0 ) { dmx[i] = copysign(1.0,dx[i]) * MIN(fabs(dx[i]), @@ -36,13 +36,10 @@ } } -static inline double approx_at_cell_interface(double* a, const int i) { - return 7.0/12.0*(a[i]+a[i+1]) - 1.0/12.0*(a[i-1]+a[i+2]); -} -static inline void monotonize(double* restrict xminus, - const double* restrict x, - double* restrict xplus, +static inline void monotonize(double* const restrict xminus, + const double* const restrict x, + double* const restrict xplus, const int i) { if ( !(xplus[i]==x[i] && x[i]==xminus[i]) @@ -63,7 +60,7 @@ } return; -} +} @@ -71,40 +68,40 @@ bool dc_flag, bool do_ppm_detect> void GRHydro_ppm1d_cxx(const int nx, const double dx, - const double* restrict rho, - const double* restrict velx, - const double* restrict vely, - const double* restrict velz, - const double* restrict eps, - const double* restrict press, - const double* restrict temp, - const double* restrict ye, - const double* restrict Bvcx, - const double* restrict Bvcy, - const double* restrict Bvcz, - const double* restrict psidc, - double* restrict rhominus, - double* restrict velxminus, - double* restrict velyminus, - double* restrict velzminus, - double* restrict epsminus, - double* restrict tempminus, - double* restrict yeminus, - double* restrict Bvcxminus, - double* restrict Bvcyminus, - double* restrict Bvczminus, - double* restrict psidcminus, - double* restrict rhoplus, - double* restrict velxplus, - double* restrict velyplus, - double* restrict velzplus, - double* restrict epsplus, - double* restrict tempplus, - double* restrict yeplus, - double* restrict Bvcxplus, - double* restrict Bvcyplus, - double* restrict Bvczplus, - double* restrict psidcplus) + const double* const restrict rho, + const double* const restrict velx, + const double* const restrict vely, + const double* const restrict velz, + const double* const restrict eps, + const double* const restrict press, + const double* const restrict temp, + const double* const restrict ye, + const double* const restrict Bvcx, + const double* const restrict Bvcy, + const double* const restrict Bvcz, + const double* const restrict psidc, + double* const restrict rhominus, + double* const restrict velxminus, + double* const restrict velyminus, + double* const restrict velzminus, + double* const restrict epsminus, + double* const restrict tempminus, + double* const restrict yeminus, + double* const restrict Bvcxminus, + double* const restrict Bvcyminus, + double* const restrict Bvczminus, + double* const restrict psidcminus, + double* const restrict rhoplus, + double* const restrict velxplus, + double* const restrict velyplus, + double* const restrict velzplus, + double* const restrict epsplus, + double* const restrict tempplus, + double* const restrict yeplus, + double* const restrict Bvcxplus, + double* const restrict Bvcyplus, + double* const restrict Bvczplus, + double* const restrict psidcplus) { DECLARE_CCTK_PARAMETERS; @@ -122,7 +119,7 @@ // Average slopes delta_m(a). See (1.7) of Colella and Woodward, p.178 // This is the expression for an even grid. - for(int i=1;i < nx-1;i++) { + for(int i=1; i < nx-1; ++i) { drho[i] = 0.5 * (rho[i+1]-rho[i-1]); d2rho[i] = rho[i+1] - 2.0 * rho[i] + rho[i-1]; dpress[i] = (press[i+1]-press[i-1]); @@ -147,7 +144,7 @@ } // Steepened slope. See (1.8) of Colella and Woodward, p.178 - for(int i=1;i<nx-1;i++) { + for(int i=1; i<nx-1; ++i) { steep(rho,drho,dmrho,i); steep(eps,deps,dmeps,i); steep(velx,dvelx,dmvelx,i); @@ -170,7 +167,7 @@ } // Initial boundary states. See (1.9) of Colella and Woodward, p.178 - for(int i=1;i<nx-2;i++) { + for(int i=1; i<nx-2; ++i) { rhoplus[i] = 0.5 * (rho[i] + rho[i+1]) + (dmrho[i] - dmrho[i+1]) * onesixth; epsplus[i] = 0.5 * (eps[i] + eps[i+1]) + @@ -204,7 +201,7 @@ } // fill minus states - for(int i=1;i<nx-2;i++) { + for(int i=1; i<nx-2; ++i) { rhominus[i+1] = rhoplus[i]; epsminus[i+1] = epsplus[i]; velxminus[i+1] = velxplus[i]; @@ -237,7 +234,7 @@ */ if(do_ppm_detect) { - for(int i=2;i<nx-2;i++) { + for(int i=2; i<nx-2; ++i) { double etatilde = 0.0; if ( d2rho[i+1]*d2rho[i-1] < 0.0 && ( fabs(rho[i+1]-rho[i-1]) - ppm_epsilon_shock @@ -260,9 +257,9 @@ } // flattening - for(int i=2;i<nx-2;i++) { - double dpress2 = press[i+2] - press[i-2]; - double dvel = velx[i+1] - velx[i-1]; + for(int i=2; i<nx-2; ++i) { + const double dpress2 = press[i+2] - press[i-2]; + const double dvel = velx[i+1] - velx[i-1]; double w=0.0; if ( (fabs(dpress[i]) > ppm_epsilon * MIN(press[i-1],press[i+1])) && (dvel < 0.0) ) @@ -280,8 +277,8 @@ } } - for(int i=2;i<nx-2;i++) { - double flatten = tilde_flatten[i]; + for(int i=2; i<nx-2; ++i) { + const double flatten = tilde_flatten[i]; rhoplus[i] = flatten * rhoplus[i] + (1.0 - flatten) * rho[i]; rhominus[i] = flatten * rhominus[i] + (1.0 - flatten) * rho[i]; epsplus[i] = flatten * epsplus[i] + (1.0 - flatten) * eps[i]; @@ -316,7 +313,7 @@ } } // flattening - for(int i=GRHydro_stencil-1;i<nx-GRHydro_stencil+1;i++) { + for(int i=GRHydro_stencil-1; i<nx-GRHydro_stencil+1; ++i) { monotonize(rhominus,rho,rhoplus,i); monotonize(epsminus,eps,epsplus,i); monotonize(velxminus,velx,velxplus,i); File [modified]: GRHydro_PPMReconstruct_drv_opt.cc Delta lines: +117 -168 =================================================================== --- trunk/src/GRHydro_PPMReconstruct_drv_opt.cc 2014-04-15 19:49:04 UTC (rev 617) +++ trunk/src/GRHydro_PPMReconstruct_drv_opt.cc 2014-04-15 19:49:07 UTC (rev 618) @@ -11,27 +11,20 @@ #include "SpaceMask.h" #include "GRHydro_PPM_opt.h" -//#define velx(i,j,k) vup(i,j,k,1) -//#define vely(i,j,k) vup(i,j,k,2) -//#define velz(i,j,k) vup(i,j,k,3) -//#define sx(i,j,k) scon(i,j,k,1) -//#define sy(i,j,k) scon(i,j,k,2) -//#define sz(i,j,k) scon(i,j,k,3) - using namespace std; extern "C" CCTK_INT GRHydro_UseGeneralCoordinates(const cGH * cctkGH); -static void oned_slice(int npt,int nstart,int nstride,CCTK_REAL const * restrict vec3d,vector<CCTK_REAL>& vec1d) +static inline void oned_slice(const int npt, const int nstart, const int nstride, CCTK_REAL const * const restrict vec3d, vector<CCTK_REAL>& vec1d) { - int i; - for (i=0; i<npt; i++)vec1d[i] = vec3d[nstart+i*nstride]; + for (int i=0; i<npt; ++i) + vec1d[i] = vec3d[nstart+i*nstride]; } -static void oned_unslice(int npt,int nstart,int nstride,vector<CCTK_REAL>& vec1d,CCTK_REAL * restrict vec3d) +static inline void oned_unslice(const int npt, const int nstart, const int nstride, const vector<CCTK_REAL>& vec1d, CCTK_REAL * const restrict vec3d) { - int i; - for (i=0; i<npt; i++)vec3d[nstart+i*nstride]=vec1d[i]; + for (int i=0; i<npt; ++i) + vec3d[nstart+i*nstride]=vec1d[i]; } /* @@ -53,10 +46,9 @@ DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; - CCTK_REAL * restrict vup; - CCTK_REAL * restrict Bprim; - //Multipatch related pointers + const CCTK_REAL * restrict vup; + const CCTK_REAL * restrict Bprim; if(GRHydro_UseGeneralCoordinates(cctkGH)) { vup=lvel; Bprim=lBvec; @@ -71,23 +63,23 @@ const int nxy=nx*ny; const int nxyz=nxy*nz; - int nmax_xyz=MAX(nx,MAX(ny,nz)); + const int nmax_xyz=MAX(nx,MAX(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"); + // bit masks and patterns for trivial Riemann Problem spacemask + const int type_bitsx = SpaceMask_GetTypeBits("Hydro_RiemannProblemX"); // trivialx = SpaceMask_GetStateBits("Hydro_RiemannProblemX","trivial"); - not_trivialx = SpaceMask_GetStateBits("Hydro_RiemannProblemX","not_trivial"); + const int not_trivialx = SpaceMask_GetStateBits("Hydro_RiemannProblemX", + "not_trivial"); - type_bitsy = SpaceMask_GetTypeBits("Hydro_RiemannProblemY"); + const int type_bitsy = SpaceMask_GetTypeBits("Hydro_RiemannProblemY"); // trivialy = SpaceMask_GetStateBits("Hydro_RiemannProblemY","trivial"); - not_trivialy = SpaceMask_GetStateBits("Hydro_RiemannProblemY","not_trivial"); + const int not_trivialy = SpaceMask_GetStateBits("Hydro_RiemannProblemY", + "not_trivial"); - type_bitsz = SpaceMask_GetTypeBits("Hydro_RiemannProblemZ"); + const int type_bitsz = SpaceMask_GetTypeBits("Hydro_RiemannProblemZ"); // trivialz = SpaceMask_GetStateBits("Hydro_RiemannProblemZ","trivial"); - not_trivialz = SpaceMask_GetStateBits("Hydro_RiemannProblemZ","not_trivial"); + const int not_trivialz = SpaceMask_GetStateBits("Hydro_RiemannProblemZ", + "not_trivial"); // if use_enhanced_ppm, allow old PPM on one level bool apply_enhanced_ppm = false; @@ -98,47 +90,6 @@ apply_enhanced_ppm = false; - int i,j,k; - vector<CCTK_REAL> rho1d(nmax_xyz); - - vector<CCTK_REAL> velx1d(nmax_xyz); - vector<CCTK_REAL> vely1d(nmax_xyz); - vector<CCTK_REAL> velz1d(nmax_xyz); - vector<CCTK_REAL> eps1d(nmax_xyz); - vector<CCTK_REAL> temp1d(nmax_xyz); - vector<CCTK_REAL> ye1d(nmax_xyz); - vector<CCTK_REAL> press1d(nmax_xyz); - vector<CCTK_REAL> Bvcx1d(nmax_xyz); - vector<CCTK_REAL> Bvcy1d(nmax_xyz); - vector<CCTK_REAL> Bvcz1d(nmax_xyz); - vector<CCTK_REAL> psidc1d(nmax_xyz); - - vector<CCTK_REAL> rho1d_plus(nmax_xyz); - vector<CCTK_REAL> velx1d_plus(nmax_xyz); - vector<CCTK_REAL> vely1d_plus(nmax_xyz); - vector<CCTK_REAL> velz1d_plus(nmax_xyz); - vector<CCTK_REAL> eps1d_plus(nmax_xyz); - vector<CCTK_REAL> temp1d_plus(nmax_xyz); - vector<CCTK_REAL> ye1d_plus(nmax_xyz); - vector<CCTK_REAL> Bvcx1d_plus(nmax_xyz); - vector<CCTK_REAL> Bvcy1d_plus(nmax_xyz); - vector<CCTK_REAL> Bvcz1d_plus(nmax_xyz); - vector<CCTK_REAL> psidc1d_plus(nmax_xyz); - - vector<CCTK_REAL> rho1d_minus(nmax_xyz); - vector<CCTK_REAL> velx1d_minus(nmax_xyz); - vector<CCTK_REAL> vely1d_minus(nmax_xyz); - vector<CCTK_REAL> velz1d_minus(nmax_xyz); - vector<CCTK_REAL> eps1d_minus(nmax_xyz); - vector<CCTK_REAL> temp1d_minus(nmax_xyz); - vector<CCTK_REAL> ye1d_minus(nmax_xyz); - vector<CCTK_REAL> Bvcx1d_minus(nmax_xyz); - vector<CCTK_REAL> Bvcy1d_minus(nmax_xyz); - vector<CCTK_REAL> Bvcz1d_minus(nmax_xyz); - vector<CCTK_REAL> psidc1d_minus(nmax_xyz); - - int nstart; - typedef void (*ppm1d_cxx_ptr_t)(const int nx, \ const double dx, \ @@ -266,18 +217,43 @@ //index = i + j*nx + k*nx*ny; -// !!$ PPM starts: - if (*flux_direction == 1){ - // !$OMP PARALLEL DO PRIVATE(i, j, k) - - for(k = GRHydro_stencil-1; k < nz - GRHydro_stencil + 1; k++) { - for(j = GRHydro_stencil-1; j < ny - GRHydro_stencil + 1; j++) { +// PPM starts: + if (*flux_direction == 1) { + #pragma omp parallel + { + vector<CCTK_REAL> rho1d(nmax_xyz), eps1d(nmax_xyz), press1d(nmax_xyz); + vector<CCTK_REAL> velx1d(nmax_xyz), vely1d(nmax_xyz), velz1d(nmax_xyz); + vector<CCTK_REAL> temp1d(nmax_xyz), ye1d(nmax_xyz); + vector<CCTK_REAL> Bvcx1d(nmax_xyz), Bvcy1d(nmax_xyz), Bvcz1d(nmax_xyz); + vector<CCTK_REAL> psidc1d(nmax_xyz); + + vector<CCTK_REAL> rho1d_plus(nmax_xyz), eps1d_plus(nmax_xyz), + press1d_plus(nmax_xyz); + vector<CCTK_REAL> velx1d_plus(nmax_xyz), vely1d_plus(nmax_xyz), + velz1d_plus(nmax_xyz); + vector<CCTK_REAL> temp1d_plus(nmax_xyz), ye1d_plus(nmax_xyz); + vector<CCTK_REAL> Bvcx1d_plus(nmax_xyz), Bvcy1d_plus(nmax_xyz), + Bvcz1d_plus(nmax_xyz); + vector<CCTK_REAL> psidc1d_plus(nmax_xyz); + + vector<CCTK_REAL> rho1d_minus(nmax_xyz), eps1d_minus(nmax_xyz), + press1d_minus(nmax_xyz); + vector<CCTK_REAL> velx1d_minus(nmax_xyz), vely1d_minus(nmax_xyz), + velz1d_minus(nmax_xyz); + vector<CCTK_REAL> temp1d_minus(nmax_xyz), ye1d_minus(nmax_xyz); + vector<CCTK_REAL> Bvcx1d_minus(nmax_xyz), Bvcy1d_minus(nmax_xyz), + Bvcz1d_minus(nmax_xyz); + vector<CCTK_REAL> psidc1d_minus(nmax_xyz); + + #pragma omp 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) { // TODO: there is no need to slice in the x direction, we can just point into the array //For reference, the slicing call is: //void oned_slice(int npt,int nstart,int nstride,CCTK_REAL* vec3d,CCTK_REAL* vec1d) - nstart=nx*(j+k*ny); + const int nstart=nx*(j+k*ny); oned_slice(nx,nstart,1,rho,rho1d); @@ -349,21 +325,46 @@ oned_unslice(nx,nstart,1,psidc1d_plus,psidcplus); } } - for (i=0; i<nx; i++) { + for (int i=0; i<nx; ++i) { SpaceMask_SetStateBits(space_mask, CCTK_GFINDEX3D(cctkGH, i,j,k), type_bitsx, not_trivialx); } } } - // !$OMP END PARALLEL DO + } // pragma omp parallel } else if (*flux_direction==2) { - // !$OMP PARALLEL DO PRIVATE(i, j, k) - //Make sure to doublecheck the bounds here! - for(k = GRHydro_stencil-1; k < nz - GRHydro_stencil + 1; k++) { - for(j = GRHydro_stencil-1; j < nx - GRHydro_stencil + 1; j++) { + #pragma omp parallel + { + vector<CCTK_REAL> rho1d(nmax_xyz), eps1d(nmax_xyz), press1d(nmax_xyz); + vector<CCTK_REAL> velx1d(nmax_xyz), vely1d(nmax_xyz), velz1d(nmax_xyz); + vector<CCTK_REAL> temp1d(nmax_xyz), ye1d(nmax_xyz); + vector<CCTK_REAL> Bvcx1d(nmax_xyz), Bvcy1d(nmax_xyz), Bvcz1d(nmax_xyz); + vector<CCTK_REAL> psidc1d(nmax_xyz); + + vector<CCTK_REAL> rho1d_plus(nmax_xyz), eps1d_plus(nmax_xyz), + press1d_plus(nmax_xyz); + vector<CCTK_REAL> velx1d_plus(nmax_xyz), vely1d_plus(nmax_xyz), + velz1d_plus(nmax_xyz); + vector<CCTK_REAL> temp1d_plus(nmax_xyz), ye1d_plus(nmax_xyz); + vector<CCTK_REAL> Bvcx1d_plus(nmax_xyz), Bvcy1d_plus(nmax_xyz), + Bvcz1d_plus(nmax_xyz); + vector<CCTK_REAL> psidc1d_plus(nmax_xyz); + + vector<CCTK_REAL> rho1d_minus(nmax_xyz), eps1d_minus(nmax_xyz), + press1d_minus(nmax_xyz); + vector<CCTK_REAL> velx1d_minus(nmax_xyz), vely1d_minus(nmax_xyz), + velz1d_minus(nmax_xyz); + vector<CCTK_REAL> temp1d_minus(nmax_xyz), ye1d_minus(nmax_xyz); + vector<CCTK_REAL> Bvcx1d_minus(nmax_xyz), Bvcy1d_minus(nmax_xyz), + Bvcz1d_minus(nmax_xyz); + vector<CCTK_REAL> psidc1d_minus(nmax_xyz); + + #pragma omp 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) { - nstart=j+k*nx*ny; + const int nstart=j+k*nx*ny; oned_slice(ny,nstart,nx,rho,rho1d); @@ -436,22 +437,47 @@ oned_unslice(ny,nstart,nx,psidc1d_plus,psidcplus); } } - for (i=0; i<ny; i++) { + for (int i=0; i<ny; ++i) { SpaceMask_SetStateBits(space_mask, CCTK_GFINDEX3D(cctkGH, j,i,k), type_bitsy, not_trivialy); } } } - // !$OMP END PARALLEL DO + } // pragma omp parallel } else if (*flux_direction==3) { - // !$OMP PARALLEL DO PRIVATE(i, j, k) - //Make sure to doublecheck the bounds here! - for(k = GRHydro_stencil-1; k < ny - GRHydro_stencil + 1; k++) { - for(j = GRHydro_stencil-1; j < nx - GRHydro_stencil + 1; j++) { + #pragma omp parallel + { + vector<CCTK_REAL> rho1d(nmax_xyz), eps1d(nmax_xyz), press1d(nmax_xyz); + vector<CCTK_REAL> velx1d(nmax_xyz), vely1d(nmax_xyz), velz1d(nmax_xyz); + vector<CCTK_REAL> temp1d(nmax_xyz), ye1d(nmax_xyz); + vector<CCTK_REAL> Bvcx1d(nmax_xyz), Bvcy1d(nmax_xyz), Bvcz1d(nmax_xyz); + vector<CCTK_REAL> psidc1d(nmax_xyz); + + vector<CCTK_REAL> rho1d_plus(nmax_xyz), eps1d_plus(nmax_xyz), + press1d_plus(nmax_xyz); + vector<CCTK_REAL> velx1d_plus(nmax_xyz), vely1d_plus(nmax_xyz), + velz1d_plus(nmax_xyz); + vector<CCTK_REAL> temp1d_plus(nmax_xyz), ye1d_plus(nmax_xyz); + vector<CCTK_REAL> Bvcx1d_plus(nmax_xyz), Bvcy1d_plus(nmax_xyz), + Bvcz1d_plus(nmax_xyz); + vector<CCTK_REAL> psidc1d_plus(nmax_xyz); + + vector<CCTK_REAL> rho1d_minus(nmax_xyz), eps1d_minus(nmax_xyz), + press1d_minus(nmax_xyz); + vector<CCTK_REAL> velx1d_minus(nmax_xyz), vely1d_minus(nmax_xyz), + velz1d_minus(nmax_xyz); + vector<CCTK_REAL> temp1d_minus(nmax_xyz), ye1d_minus(nmax_xyz); + vector<CCTK_REAL> Bvcx1d_minus(nmax_xyz), Bvcy1d_minus(nmax_xyz), + Bvcz1d_minus(nmax_xyz); + vector<CCTK_REAL> psidc1d_minus(nmax_xyz); + + #pragma omp 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) { - nstart=j+k*nx; + const int nstart=j+k*nx; oned_slice(nz,nstart,nxy,rho,rho1d); @@ -524,90 +550,13 @@ oned_unslice(nz,nstart,nxy,psidc1d_plus,psidcplus); } } - for (i=0; i<nz; i++) { + for (int i=0; i<nz; ++i) { SpaceMask_SetStateBits(space_mask, CCTK_GFINDEX3D(cctkGH, j,k,i), type_bitsz, not_trivialz); } } } - // !$OMP END PARALLEL DO - + } // pragma omp parallel } } - - -// trivial_rp(:,j,k) = .false. -// call PPMtyc(nx,CCTK_DELTA_SPACE(1),rho(:,j,k),velx(:,j,k),& -// vely(:,j,k),velz(:,j,k),eps(:,j,k),temperature(:,j,k),Y_e(:,j,k),& -// press(:,j,k),rhominus(:,j,k),& -// velxminus(:,j,k),velyminus(:,j,k),velzminus(:,j,k),& -// epsminus(:,j,k),tempminus(:,j,k), y_e_minus(:,j,k),& -// rhoplus(:,j,k),& -// velxplus(:,j,k),velyplus(:,j,k),velzplus(:,j,k),& -// epsplus(:,j,k),tempplus(:,j,k),y_e_plus(:,j,k)) -// do i = 1, nx -// if (trivial_rp(i,j,k)) then -// SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, trivialx) -// else -// SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, not_trivialx) -// end if -// end do -// end do -// end do -// !$OMP END PARALLEL DO -// else if (flux_direction == 2) then -// !$OMP PARALLEL DO PRIVATE(i, j, k) -// do k = GRHydro_stencil, nz - GRHydro_stencil + 1 -// do j = GRHydro_stencil, nx - GRHydro_stencil + 1 -// trivial_rp(j,:,k) = .false. -// call PPMtyc(ny,CCTK_DELTA_SPACE(2),rho(j,:,k),& -// vely(j,:,k),velz(j,:,k),velx(j,:,k),& -// eps(j,:,k),temperature(j,:,k),Y_e(j,:,k), & -// press(j,:,k),rhominus(j,:,k),& -// velyminus(j,:,k),velzminus(j,:,k),velxminus(j,:,k),& -// epsminus(j,:,k),tempminus(j,:,k), y_e_minus(j,:,k),rhoplus(j,:,k),& -// velyplus(j,:,k),velzplus(j,:,k),velxplus(j,:,k),& -// epsplus(j,:,k),tempplus(j,:,k),y_e_plus(j,:,k)) -// do i = 1, ny -// if (trivial_rp(j,i,k)) then -// SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, trivialy) -// else -// SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, not_trivialy) -// end if -// end do -// end do -// end do -// !$OMP END PARALLEL DO - -// else if (flux_direction == 3) then -// !$OMP PARALLEL DO PRIVATE(i, j, k) -// do k = GRHydro_stencil, ny - GRHydro_stencil + 1 -// do j = GRHydro_stencil, nx - GRHydro_stencil + 1 -// trivial_rp(j,k,:) = .false. -// call PPMtyc(nz,CCTK_DELTA_SPACE(3),rho(j,k,:),& -// velz(j,k,:),velx(j,k,:),vely(j,k,:),& -// eps(j,k,:),temperature(j,k,:),Y_e(j,k,:), press(j,k,:),rhominus(j,k,:),& -// velzminus(j,k,:),velxminus(j,k,:),velyminus(j,k,:),& -// epsminus(j,k,:),tempminus(j,k,:),y_e_minus(j,k,:),rhoplus(j,k,:),& -// velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),& -// epsplus(j,k,:),tempplus(j,k,:),y_e_plus(j,k,:)) -// do i = 1, nz -// if (trivial_rp(j,k,i)) then -// SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, trivialz) -// else -// SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, not_trivialz) -// end if -// end do -// end do -// end do -// !$OMP END PARALLEL DO -// else -// call CCTK_WARN(0, "Flux direction not x,y,z") -// end if -// !!$ PPM ends. - -// deallocate(trivial_rp) - -// end subroutine GRHydro_PPMReconstruct_drv -

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 617)
by rhaas＠tapir.caltech.edu 15 Apr '14

15 Apr '14

User: rhaas Date: 2014/04/15 02:49 PM Modified: /trunk/src/ GRHydro_PPMReconstruct_drv_opt.cc Log: GRHydro: unify usage of a[0] vs *a File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_PPMReconstruct_drv_opt.cc Delta lines: +3 -3 =================================================================== --- trunk/src/GRHydro_PPMReconstruct_drv_opt.cc 2014-04-15 19:49:02 UTC (rev 616) +++ trunk/src/GRHydro_PPMReconstruct_drv_opt.cc 2014-04-15 19:49:04 UTC (rev 617) @@ -267,7 +267,7 @@ //index = i + j*nx + k*nx*ny; // !!$ PPM starts: - if (flux_direction[0] == 1){ + if (*flux_direction == 1){ // !$OMP PARALLEL DO PRIVATE(i, j, k) for(k = GRHydro_stencil-1; k < nz - GRHydro_stencil + 1; k++) { @@ -356,7 +356,7 @@ } // !$OMP END PARALLEL DO - } else if (flux_direction[0]==2) { + } else if (*flux_direction==2) { // !$OMP PARALLEL DO PRIVATE(i, j, k) //Make sure to doublecheck the bounds here! @@ -444,7 +444,7 @@ } // !$OMP END PARALLEL DO - } else if (flux_direction[0]==3) { + } else if (*flux_direction==3) { // !$OMP PARALLEL DO PRIVATE(i, j, k) //Make sure to doublecheck the bounds here!

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 616)
by rhaas＠tapir.caltech.edu 15 Apr '14

15 Apr '14

User: rhaas Date: 2014/04/15 02:49 PM Modified: /trunk/src/ GRHydro_PPMReconstruct_drv_opt.cc Log: GRHydro: remove string comparisons from C++ PPM driver File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_PPMReconstruct_drv_opt.cc Delta lines: +25 -24 =================================================================== --- trunk/src/GRHydro_PPMReconstruct_drv_opt.cc 2014-04-15 19:48:59 UTC (rev 615) +++ trunk/src/GRHydro_PPMReconstruct_drv_opt.cc 2014-04-15 19:49:02 UTC (rev 616) @@ -256,6 +256,7 @@ const bool do_mhd = *evolve_MHD; const bool dc_flag = clean_divergence; const bool do_ppm_detect = ppm_detect; + const bool do_Avec = CCTK_Equals(Bvec_evolution_method,"GRHydro_Avec"); ppm1d_cxx_func = ppm1d_cxx_funcs[MKINDEX(do_temp,do_ye,do_mhd,dc_flag,do_ppm_detect,apply_enhanced_ppm)]; if(ppm1d_cxx_func == NULL) { CCTK_VError(__LINE__, __FILE__, CCTK_THORNSTRING, @@ -287,22 +288,22 @@ oned_slice(nx,nstart,1,eps,eps1d); oned_slice(nx,nstart,1,press,press1d); - if(Y_e!=NULL) { + if(do_ye) { oned_slice(nx,nstart,1,Y_e,ye1d); - if(temperature!=NULL){ + if(do_temp){ oned_slice(nx,nstart,1,temperature,temp1d); } } - if(Bvec!=NULL) { + if(do_mhd) { oned_slice(nx,nstart, 1,Bprim,Bvcx1d); oned_slice(nx,nstart+nxyz, 1,Bprim,Bvcy1d); oned_slice(nx,nstart+2*nxyz,1,Bprim,Bvcz1d); if(clean_divergence) { oned_slice(nx,nstart,1,psidc,psidc1d); } - } else if (CCTK_Equals(Bvec_evolution_method,"GRHydro_Avec")) { - CCTK_WARN(0, "Someone needs to figure out which variables to pass for Avec"); + if (do_Avec) + CCTK_WARN(0, "Someone needs to figure out which variables to pass for Avec"); } ppm1d_cxx_func(cctk_lsh[0],CCTK_DELTA_SPACE(0), @@ -327,16 +328,16 @@ oned_unslice(nx,nstart,1,eps1d_minus,epsminus); oned_unslice(nx,nstart,1,eps1d_plus,epsplus); - if(CCTK_Equals(Y_e_evolution_method,"GRHydro")) { + if(do_ye) { oned_unslice(nx,nstart,1,ye1d_minus,Y_e_minus); oned_unslice(nx,nstart,1,ye1d_plus,Y_e_plus); - if(CCTK_Equals(temperature_evolution_method,"GRHydro")) { + if(do_temp) { oned_unslice(nx,nstart,1,temp1d_minus,tempminus); oned_unslice(nx,nstart,1,temp1d_plus,tempplus); } } - if(CCTK_Equals(Bvec_evolution_method,"GRHydro")) { + if(do_mhd) { oned_unslice(nx,nstart,1,Bvcx1d_minus,Bvecxminus); oned_unslice(nx,nstart,1,Bvcy1d_minus,Bvecyminus); oned_unslice(nx,nstart,1,Bvcz1d_minus,Bveczminus); @@ -373,22 +374,22 @@ oned_slice(ny,nstart,nx,eps,eps1d); oned_slice(ny,nstart,nx,press,press1d); - if(Y_e!=NULL) { + if(do_ye) { oned_slice(ny,nstart,nx,Y_e,ye1d); - if(temperature!=NULL){ + if(do_temp){ oned_slice(ny,nstart,nx,temperature,temp1d); } } - if(Bvec!=NULL) { + if(do_mhd) { oned_slice(ny,nstart, nx,Bprim,Bvcx1d); oned_slice(ny,nstart+nxyz, nx,Bprim,Bvcy1d); oned_slice(ny,nstart+2*nxyz,nx,Bprim,Bvcz1d); if(clean_divergence) { oned_slice(ny,nstart,nx,psidc,psidc1d); } - } else if (CCTK_Equals(Bvec_evolution_method,"GRHydro_Avec")) { - CCTK_WARN(0, "Someone needs to figure out which variables to pass for Avec"); + if (do_Avec) + CCTK_WARN(0, "Someone needs to figure out which variables to pass for Avec"); } ppm1d_cxx_func @@ -414,16 +415,16 @@ oned_unslice(ny,nstart,nx,eps1d_minus,epsminus); oned_unslice(ny,nstart,nx,eps1d_plus,epsplus); - if(CCTK_Equals(Y_e_evolution_method,"GRHydro")) { + if(do_ye) { oned_unslice(ny,nstart,nx,ye1d_minus,Y_e_minus); oned_unslice(ny,nstart,nx,ye1d_plus,Y_e_plus); - if(CCTK_Equals(temperature_evolution_method,"GRHydro")) { + if(do_temp) { oned_unslice(ny,nstart,nx,temp1d_minus,tempminus); oned_unslice(ny,nstart,nx,temp1d_plus,tempplus); } } - if(CCTK_Equals(Bvec_evolution_method,"GRHydro")) { + if(do_mhd) { oned_unslice(ny,nstart,nx,Bvcx1d_minus,Bvecxminus); oned_unslice(ny,nstart,nx,Bvcy1d_minus,Bvecyminus); oned_unslice(ny,nstart,nx,Bvcz1d_minus,Bveczminus); @@ -461,22 +462,22 @@ oned_slice(nz,nstart,nxy,eps,eps1d); oned_slice(nz,nstart,nxy,press,press1d); - if(Y_e!=NULL) { + if(do_ye) { oned_slice(nz,nstart,nxy,Y_e,ye1d); - if(temperature!=NULL){ + if(do_temp){ oned_slice(nz,nstart,nxy,temperature,temp1d); } } - if(Bvec!=NULL) { + if(do_mhd) { oned_slice(nz,nstart, nxy,Bprim,Bvcx1d); oned_slice(nz,nstart+nxyz, nxy,Bprim,Bvcy1d); oned_slice(nz,nstart+2*nxyz,nxy,Bprim,Bvcz1d); if(clean_divergence) { oned_slice(nz,nstart,nxy,psidc,psidc1d); } - } else if (CCTK_Equals(Bvec_evolution_method,"GRHydro_Avec")) { - CCTK_WARN(0, "Someone needs to figure out which variables to pass for Avec"); + if (do_Avec) + CCTK_WARN(0, "Someone needs to figure out which variables to pass for Avec"); } ppm1d_cxx_func @@ -502,16 +503,16 @@ oned_unslice(nz,nstart,nxy,eps1d_minus,epsminus); oned_unslice(nz,nstart,nxy,eps1d_plus,epsplus); - if(CCTK_Equals(Y_e_evolution_method,"GRHydro")) { + if(do_ye) { oned_unslice(nz,nstart,nxy,ye1d_minus,Y_e_minus); oned_unslice(nz,nstart,nxy,ye1d_plus,Y_e_plus); - if(CCTK_Equals(temperature_evolution_method,"GRHydro")) { + if(do_temp) { oned_unslice(nz,nstart,nxy,temp1d_minus,tempminus); oned_unslice(nz,nstart,nxy,temp1d_plus,tempplus); } } - if(CCTK_Equals(Bvec_evolution_method,"GRHydro")) { + if(do_mhd) { oned_unslice(nz,nstart,nxy,Bvcx1d_minus,Bvecxminus); oned_unslice(nz,nstart,nxy,Bvcy1d_minus,Bvecyminus); oned_unslice(nz,nstart,nxy,Bvcz1d_minus,Bveczminus);

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 615)
by rhaas＠tapir.caltech.edu 15 Apr '14

15 Apr '14

User: rhaas Date: 2014/04/15 02:48 PM Modified: /trunk/src/ GRHydro_PPMReconstruct_drv_opt.cc Log: GRHydro: support array padding in optimized C++ PPM File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_PPMReconstruct_drv_opt.cc Delta lines: +11 -24 =================================================================== --- trunk/src/GRHydro_PPMReconstruct_drv_opt.cc 2014-04-15 19:48:56 UTC (rev 614) +++ trunk/src/GRHydro_PPMReconstruct_drv_opt.cc 2014-04-15 19:48:59 UTC (rev 615) @@ -65,13 +65,12 @@ Bprim=Bvec; } - int nx=cctk_lsh[0]; - int ny=cctk_lsh[1]; - int nz=cctk_lsh[2]; - int nxy=nx*ny; - int nxyz=nxy*nz; + const int nx=cctk_ash[0]; + const int ny=cctk_ash[1]; + const int nz=cctk_ash[2]; + const int nxy=nx*ny; + const int nxyz=nxy*nz; - //Is there a more efficient way? int nmax_xyz=MAX(nx,MAX(ny,nz)); int type_bitsx,type_bitsy,type_bitsz; @@ -306,11 +305,7 @@ CCTK_WARN(0, "Someone needs to figure out which variables to pass for Avec"); } - //This needs to be fixed into a loop for the the template version! - //WE'll need booleans arguments corresponding to the following - //ppm1d_cxx<do_temp,do_ye,do_mhd,clean_divergence,ppm_detect> - - ppm1d_cxx_func(nx,CCTK_DELTA_SPACE(0), + ppm1d_cxx_func(cctk_lsh[0],CCTK_DELTA_SPACE(0), &rho1d[0],&velx1d[0],&vely1d[0],&velz1d[0],&eps1d[0],&press1d[0],&temp1d[0], &ye1d[0],&Bvcx1d[0],&Bvcy1d[0],&Bvcz1d[0],&psidc1d[0], &rho1d_minus[0],&velx1d_minus[0],&vely1d_minus[0],&velz1d_minus[0],&eps1d_minus[0],&temp1d_minus[0], @@ -354,7 +349,7 @@ } } for (i=0; i<nx; i++) { - SpaceMask_SetStateBits(space_mask, i+j*nx+k*nx*ny, type_bitsx, not_trivialx); + SpaceMask_SetStateBits(space_mask, CCTK_GFINDEX3D(cctkGH, i,j,k), type_bitsx, not_trivialx); } } } @@ -396,12 +391,8 @@ CCTK_WARN(0, "Someone needs to figure out which variables to pass for Avec"); } - //This needs to be fixed into a loop for the the template version! - //WE'll need booleans arguments corresponding to the following - //ppm1d_cxx<do_temp,do_ye,do_mhd,clean_divergence,ppm_detect> - ppm1d_cxx_func - (ny,CCTK_DELTA_SPACE(1), + (cctk_lsh[1],CCTK_DELTA_SPACE(1), &rho1d[0],&vely1d[0],&velz1d[0],&velx1d[0],&eps1d[0],&press1d[0],&temp1d[0], &ye1d[0],&Bvcy1d[0],&Bvcz1d[0],&Bvcx1d[0],&psidc1d[0], &rho1d_minus[0],&vely1d_minus[0],&velz1d_minus[0],&velx1d_minus[0],&eps1d_minus[0],&temp1d_minus[0], @@ -445,7 +436,7 @@ } } for (i=0; i<ny; i++) { - SpaceMask_SetStateBits(space_mask, j+i*nx+k*nx*ny, type_bitsy, not_trivialy); + SpaceMask_SetStateBits(space_mask, CCTK_GFINDEX3D(cctkGH, j,i,k), type_bitsy, not_trivialy); } } @@ -488,12 +479,8 @@ CCTK_WARN(0, "Someone needs to figure out which variables to pass for Avec"); } - //This needs to be fixed into a loop for the the template version! - //WE'll need booleans arguments corresponding to the following - //ppm1d_cxx<do_temp,do_ye,do_mhd,clean_divergence,ppm_detect> - ppm1d_cxx_func - (nz,CCTK_DELTA_SPACE(2), + (cctk_lsh[2],CCTK_DELTA_SPACE(2), &rho1d[0],&velz1d[0],&velx1d[0],&vely1d[0],&eps1d[0],&press1d[0],&temp1d[0], &ye1d[0],&Bvcz1d[0],&Bvcx1d[0],&Bvcy1d[0],&psidc1d[0], &rho1d_minus[0],&velz1d_minus[0],&velx1d_minus[0],&vely1d_minus[0],&eps1d_minus[0],&temp1d_minus[0], @@ -537,7 +524,7 @@ } } for (i=0; i<nz; i++) { - SpaceMask_SetStateBits(space_mask, j+k*nx+i*nx*ny, type_bitsz, not_trivialz); + SpaceMask_SetStateBits(space_mask, CCTK_GFINDEX3D(cctkGH, j,k,i), type_bitsz, not_trivialz); } }

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 614)
by rhaas＠tapir.caltech.edu 15 Apr '14

15 Apr '14

User: rhaas Date: 2014/04/15 02:48 PM Modified: /trunk/src/ GRHydro_PPMReconstruct_drv_opt.cc Log: GRHydro: add MHD multipatch support to optimized C++ PPM File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_PPMReconstruct_drv_opt.cc Delta lines: +12 -9 =================================================================== --- trunk/src/GRHydro_PPMReconstruct_drv_opt.cc 2014-04-15 19:48:54 UTC (rev 613) +++ trunk/src/GRHydro_PPMReconstruct_drv_opt.cc 2014-04-15 19:48:56 UTC (rev 614) @@ -54,12 +54,15 @@ DECLARE_CCTK_PARAMETERS; CCTK_REAL * restrict vup; + CCTK_REAL * restrict Bprim; //Multipatch related pointers if(GRHydro_UseGeneralCoordinates(cctkGH)) { vup=lvel; + Bprim=lBvec; } else { vup=vel; + Bprim=Bvec; } int nx=cctk_lsh[0]; @@ -293,9 +296,9 @@ } if(Bvec!=NULL) { - oned_slice(nx,nstart, 1,Bvec,Bvcx1d); - oned_slice(nx,nstart+nxyz, 1,Bvec,Bvcy1d); - oned_slice(nx,nstart+2*nxyz,1,Bvec,Bvcz1d); + oned_slice(nx,nstart, 1,Bprim,Bvcx1d); + oned_slice(nx,nstart+nxyz, 1,Bprim,Bvcy1d); + oned_slice(nx,nstart+2*nxyz,1,Bprim,Bvcz1d); if(clean_divergence) { oned_slice(nx,nstart,1,psidc,psidc1d); } @@ -383,9 +386,9 @@ } if(Bvec!=NULL) { - oned_slice(ny,nstart, nx,Bvec,Bvcx1d); - oned_slice(ny,nstart+nxyz, nx,Bvec,Bvcy1d); - oned_slice(ny,nstart+2*nxyz,nx,Bvec,Bvcz1d); + oned_slice(ny,nstart, nx,Bprim,Bvcx1d); + oned_slice(ny,nstart+nxyz, nx,Bprim,Bvcy1d); + oned_slice(ny,nstart+2*nxyz,nx,Bprim,Bvcz1d); if(clean_divergence) { oned_slice(ny,nstart,nx,psidc,psidc1d); } @@ -475,9 +478,9 @@ } if(Bvec!=NULL) { - oned_slice(nz,nstart, nxy,Bvec,Bvcx1d); - oned_slice(nz,nstart+nxyz, nxy,Bvec,Bvcy1d); - oned_slice(nz,nstart+2*nxyz,nxy,Bvec,Bvcz1d); + oned_slice(nz,nstart, nxy,Bprim,Bvcx1d); + oned_slice(nz,nstart+nxyz, nxy,Bprim,Bvcy1d); + oned_slice(nz,nstart+2*nxyz,nxy,Bprim,Bvcz1d); if(clean_divergence) { oned_slice(nz,nstart,nxy,psidc,psidc1d); }

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 613)
by rhaas＠tapir.caltech.edu 15 Apr '14

15 Apr '14

User: rhaas Date: 2014/04/15 02:48 PM Modified: /trunk/src/ GRHydro_PPMReconstruct_drv_opt.cc Log: GRHydro: fix z direction in optimized C++ PPM File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_PPMReconstruct_drv_opt.cc Delta lines: +1 -1 =================================================================== --- trunk/src/GRHydro_PPMReconstruct_drv_opt.cc 2014-04-15 19:48:52 UTC (rev 612) +++ trunk/src/GRHydro_PPMReconstruct_drv_opt.cc 2014-04-15 19:48:54 UTC (rev 613) @@ -490,7 +490,7 @@ //ppm1d_cxx<do_temp,do_ye,do_mhd,clean_divergence,ppm_detect> ppm1d_cxx_func - (ny,CCTK_DELTA_SPACE(1), + (nz,CCTK_DELTA_SPACE(2), &rho1d[0],&velz1d[0],&velx1d[0],&vely1d[0],&eps1d[0],&press1d[0],&temp1d[0], &ye1d[0],&Bvcz1d[0],&Bvcx1d[0],&Bvcy1d[0],&psidc1d[0], &rho1d_minus[0],&velz1d_minus[0],&velx1d_minus[0],&vely1d_minus[0],&eps1d_minus[0],&temp1d_minus[0],

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 612)
by rhaas＠tapir.caltech.edu 15 Apr '14

15 Apr '14

User: rhaas Date: 2014/04/15 02:48 PM Modified: /trunk/src/ GRHydro_PPMReconstruct_drv_opt.cc Log: GRhydro: correct argument types for oned_slice and unslice File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_PPMReconstruct_drv_opt.cc Delta lines: +2 -4 =================================================================== --- trunk/src/GRHydro_PPMReconstruct_drv_opt.cc 2014-04-15 19:48:49 UTC (rev 611) +++ trunk/src/GRHydro_PPMReconstruct_drv_opt.cc 2014-04-15 19:48:52 UTC (rev 612) @@ -22,18 +22,16 @@ extern "C" CCTK_INT GRHydro_UseGeneralCoordinates(const cGH * cctkGH); -static void oned_slice(int npt,int nstart,int nstride,CCTK_REAL* vec3d,vector<CCTK_REAL> vec1d) +static void oned_slice(int npt,int nstart,int nstride,CCTK_REAL const * restrict vec3d,vector<CCTK_REAL>& vec1d) { int i; for (i=0; i<npt; i++)vec1d[i] = vec3d[nstart+i*nstride]; - return; } -static void oned_unslice(int npt,int nstart,int nstride,vector<CCTK_REAL> vec1d,CCTK_REAL* vec3d) +static void oned_unslice(int npt,int nstart,int nstride,vector<CCTK_REAL>& vec1d,CCTK_REAL * restrict vec3d) { int i; for (i=0; i<npt; i++)vec3d[nstart+i*nstride]=vec1d[i]; - return; } /*

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 611)
by rhaas＠tapir.caltech.edu 15 Apr '14

15 Apr '14

User: rhaas Date: 2014/04/15 02:48 PM Modified: /trunk/src/ GRHydro_PPMReconstruct_drv_opt.cc Log: GRHydro: correct vel slicing in optimized C++ PPM reconstruction File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_PPMReconstruct_drv_opt.cc Delta lines: +6 -6 =================================================================== --- trunk/src/GRHydro_PPMReconstruct_drv_opt.cc 2014-04-15 19:48:47 UTC (rev 610) +++ trunk/src/GRHydro_PPMReconstruct_drv_opt.cc 2014-04-15 19:48:49 UTC (rev 611) @@ -281,8 +281,8 @@ oned_slice(nx,nstart,1,rho,rho1d); oned_slice(nx,nstart, 1,vup,velx1d); - oned_slice(nx,nstart+nxyz, 1,vup,velx1d); - oned_slice(nx,nstart+2*nxyz,1,vup,velx1d); + oned_slice(nx,nstart+nxyz, 1,vup,vely1d); + oned_slice(nx,nstart+2*nxyz,1,vup,velz1d); oned_slice(nx,nstart,1,eps,eps1d); oned_slice(nx,nstart,1,press,press1d); @@ -371,8 +371,8 @@ oned_slice(ny,nstart,nx,rho,rho1d); oned_slice(ny,nstart, nx,vup,velx1d); - oned_slice(ny,nstart+nxyz, nx,vup,velx1d); - oned_slice(ny,nstart+2*nxyz,nx,vup,velx1d); + oned_slice(ny,nstart+nxyz, nx,vup,vely1d); + oned_slice(ny,nstart+2*nxyz,nx,vup,velz1d); oned_slice(ny,nstart,nx,eps,eps1d); oned_slice(ny,nstart,nx,press,press1d); @@ -463,8 +463,8 @@ oned_slice(nz,nstart,nxy,rho,rho1d); oned_slice(nz,nstart, nxy,vup,velx1d); - oned_slice(nz,nstart+nxyz, nxy,vup,velx1d); - oned_slice(nz,nstart+2*nxyz,nxy,vup,velx1d); + oned_slice(nz,nstart+nxyz, nxy,vup,vely1d); + oned_slice(nz,nstart+2*nxyz,nxy,vup,velz1d); oned_slice(nz,nstart,nxy,eps,eps1d); oned_slice(nz,nstart,nxy,press,press1d);

1 0

Jump to page:

2026

2025

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

Commits