[Commits] [svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 618)
rhaas at tapir.caltech.edu
rhaas at tapir.caltech.edu
Tue Apr 15 14:49:07 CDT 2014
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
-
More information about the Commits
mailing list