[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