Commits April 2014

commits@lists.einsteintoolkit.org

8 participants
175 discussions

[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

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

15 Apr '14

User: rhaas Date: 2014/04/15 02:48 PM Modified: /trunk/src/ GRHydro_PPM.cc, GRHydro_PPMReconstruct_drv_opt.cc Log: GRHydro: adjust loop bounds in C++ PPM to match F90 code some are just affecting data in ghost zones, some are required File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_PPM.cc Delta lines: +1 -1 =================================================================== --- trunk/src/GRHydro_PPM.cc 2014-04-15 19:48:44 UTC (rev 609) +++ trunk/src/GRHydro_PPM.cc 2014-04-15 19:48:47 UTC (rev 610) @@ -316,7 +316,7 @@ } } // flattening - for(int i=2;i<nx-2;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: +8 -6 =================================================================== --- trunk/src/GRHydro_PPMReconstruct_drv_opt.cc 2014-04-15 19:48:44 UTC (rev 609) +++ trunk/src/GRHydro_PPMReconstruct_drv_opt.cc 2014-04-15 19:48:47 UTC (rev 610) @@ -269,9 +269,11 @@ if (flux_direction[0] == 1){ // !$OMP PARALLEL DO PRIVATE(i, j, k) - for(k = GRHydro_stencil-1; k < nz - GRHydro_stencil; k++) { - for(j = GRHydro_stencil-1; j < ny - GRHydro_stencil; j++) { + for(k = GRHydro_stencil-1; k < nz - GRHydro_stencil + 1; k++) { + for(j = GRHydro_stencil-1; j < ny - GRHydro_stencil + 1; 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); @@ -361,8 +363,8 @@ // !$OMP PARALLEL DO PRIVATE(i, j, k) //Make sure to doublecheck the bounds here! - for(k = GRHydro_stencil-1; k < nz - GRHydro_stencil; k++) { - for(j = GRHydro_stencil-1; j < nx - GRHydro_stencil; j++) { + for(k = GRHydro_stencil-1; k < nz - GRHydro_stencil + 1; k++) { + for(j = GRHydro_stencil-1; j < nx - GRHydro_stencil + 1; j++) { nstart=j+k*nx*ny; @@ -453,8 +455,8 @@ // !$OMP PARALLEL DO PRIVATE(i, j, k) //Make sure to doublecheck the bounds here! - for(k = GRHydro_stencil-1; k < ny - GRHydro_stencil; k++) { - for(j = GRHydro_stencil-1; j < nx - GRHydro_stencil; j++) { + for(k = GRHydro_stencil-1; k < ny - GRHydro_stencil + 1; k++) { + for(j = GRHydro_stencil-1; j < nx - GRHydro_stencil + 1; j++) { nstart=j+k*nx;

1 0

Jump to page:

2026

2025

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

Commits April 2014