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