From 08cbb62fe2fedf3b01c1fe5acedf3cb23b0e5034 Mon Sep 17 00:00:00 2001 From: rhaas Date: Tue, 13 Aug 2013 14:56:31 +0000 Subject: GRHydro: First draft of templated version of PPM Reconstruction call, Mostly involves slicing and unslicing 3-d grids into 1-d segments and calling the new templated routines in GRHydro_PPM.cc From: Josh Faber git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@574 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45 --- src/GRHydro_PPMReconstruct_drv_opt.cc | 521 ++++++++++++++++++++++++++++++++++ src/GRHydro_PPM_opt.h | 38 +++ 2 files changed, 559 insertions(+) create mode 100644 src/GRHydro_PPMReconstruct_drv_opt.cc create mode 100644 src/GRHydro_PPM_opt.h diff --git a/src/GRHydro_PPMReconstruct_drv_opt.cc b/src/GRHydro_PPMReconstruct_drv_opt.cc new file mode 100644 index 0000000..43df1b6 --- /dev/null +++ b/src/GRHydro_PPMReconstruct_drv_opt.cc @@ -0,0 +1,521 @@ +#include +#include +#include +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#define MIN(a,b) (((a)<(b))?(a):(b)) +#define MAX(a,b) (((a)>(b))?(a):(b)) + +#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); + +void oned_slice(int npt,int nstart,int nstride,CCTK_REAL* vec3d,vector vec1d) +{ + int i; + for (i=0; i vec1d,CCTK_REAL* vec3d) +{ + int i; + for (i=0; i rho1d(nmax_xyz); + + vector velx1d(nmax_xyz); + vector vely1d(nmax_xyz); + vector velz1d(nmax_xyz); + vector eps1d(nmax_xyz); + vector temp1d(nmax_xyz); + vector ye1d(nmax_xyz); + vector press1d(nmax_xyz); + vector Bvcx1d(nmax_xyz); + vector Bvcy1d(nmax_xyz); + vector Bvcz1d(nmax_xyz); + vector psidc1d(nmax_xyz); + + vector rho1d_plus(nmax_xyz); + vector velx1d_plus(nmax_xyz); + vector vely1d_plus(nmax_xyz); + vector velz1d_plus(nmax_xyz); + vector eps1d_plus(nmax_xyz); + vector temp1d_plus(nmax_xyz); + vector ye1d_plus(nmax_xyz); + vector Bvcx1d_plus(nmax_xyz); + vector Bvcy1d_plus(nmax_xyz); + vector Bvcz1d_plus(nmax_xyz); + vector psidc1d_plus(nmax_xyz); + + vector rho1d_minus(nmax_xyz); + vector velx1d_minus(nmax_xyz); + vector vely1d_minus(nmax_xyz); + vector velz1d_minus(nmax_xyz); + vector eps1d_minus(nmax_xyz); + vector temp1d_minus(nmax_xyz); + vector ye1d_minus(nmax_xyz); + vector Bvcx1d_minus(nmax_xyz); + vector Bvcy1d_minus(nmax_xyz); + vector Bvcz1d_minus(nmax_xyz); + vector psidc1d_minus(nmax_xyz); + + int npt,nstart,nstride; + + //index = i + j*nx + k*nx*ny; + +// !!$ PPM starts: + 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 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); + + 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,1,eps,eps1d); + oned_slice(nx,nstart,1,press,press1d); + + if(Y_e!=NULL) { + oned_slice(nx,nstart,1,Y_e,ye1d); + if(temperature!=NULL){ + oned_slice(nx,nstart,1,temperature,temp1d); + } + } + + 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); + 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"); + } + + //This needs to be fixed into a loop for the the template version! + //WE'll need booleans arguments corresponding to the following + //ppm1d_cxx + + ppm1d_cxx + (nx,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], + &ye1d_minus[0],&Bvcx1d_minus[0],&Bvcy1d_minus[0],&Bvcz1d_minus[0],&psidc1d_minus[0], + &rho1d_plus[0],&velx1d_plus[0],&vely1d_plus[0],&velz1d_plus[0],&eps1d_plus[0],&temp1d_plus[0], + &ye1d_plus[0],&Bvcx1d_plus[0],&Bvcy1d_plus[0],&Bvcz1d_plus[0],&psidc1d_plus[0]); + + + oned_unslice(nx,nstart,1,rho1d_minus,rhominus); + oned_unslice(nx,nstart,1,rho1d_plus,rhoplus); + + oned_unslice(nx,nstart,1,velx1d_minus,velxminus); + oned_unslice(nx,nstart,1,vely1d_minus,velyminus); + oned_unslice(nx,nstart,1,velz1d_minus,velzminus); + oned_unslice(nx,nstart,1,velx1d_plus,velxplus); + oned_unslice(nx,nstart,1,vely1d_plus,velyplus); + oned_unslice(nx,nstart,1,velz1d_plus,velzplus); + + oned_unslice(nx,nstart,1,eps1d_minus,epsminus); + oned_unslice(nx,nstart,1,eps1d_plus,epsplus); + + if(CCTK_Equals(Y_e_evolution_method,"GRHydro")) { + 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")) { + oned_unslice(nx,nstart,1,temp1d_minus,tempminus); + oned_unslice(nx,nstart,1,temp1d_plus,tempplus); + } + } + + if(CCTK_Equals(Bvec_evolution_method,"GRHydro")) { + oned_unslice(nx,nstart,1,Bvcx1d_minus,Bvecxminus); + oned_unslice(nx,nstart,1,Bvcy1d_minus,Bvecyminus); + oned_unslice(nx,nstart,1,Bvcz1d_minus,Bveczminus); + oned_unslice(nx,nstart,1,Bvcx1d_plus,Bvecxplus); + oned_unslice(nx,nstart,1,Bvcy1d_plus,Bvecyplus); + oned_unslice(nx,nstart,1,Bvcz1d_plus,Bveczplus); + if(clean_divergence) { + oned_unslice(nx,nstart,1,psidc1d_minus,psidcminus); + oned_unslice(nx,nstart,1,psidc1d_plus,psidcplus); + } + } + for (i=0; i + + ppm1d_cxx + (ny,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], + &ye1d_minus[0],&Bvcy1d_minus[0],&Bvcz1d_minus[0],&Bvcx1d_minus[0],&psidc1d_minus[0], + &rho1d_plus[0],&vely1d_plus[0],&velz1d_plus[0],&velx1d_plus[0],&eps1d_plus[0],&temp1d_plus[0], + &ye1d_plus[0],&Bvcy1d_plus[0],&Bvcz1d_plus[0],&Bvcx1d_plus[0],&psidc1d_plus[0]); + + + oned_unslice(ny,nstart,nx,rho1d_minus,rhominus); + oned_unslice(ny,nstart,nx,rho1d_plus,rhoplus); + + oned_unslice(ny,nstart,nx,velx1d_minus,velxminus); + oned_unslice(ny,nstart,nx,vely1d_minus,velyminus); + oned_unslice(ny,nstart,nx,velz1d_minus,velzminus); + oned_unslice(ny,nstart,nx,velx1d_plus,velxplus); + oned_unslice(ny,nstart,nx,vely1d_plus,velyplus); + oned_unslice(ny,nstart,nx,velz1d_plus,velzplus); + + oned_unslice(ny,nstart,nx,eps1d_minus,epsminus); + oned_unslice(ny,nstart,nx,eps1d_plus,epsplus); + + if(CCTK_Equals(Y_e_evolution_method,"GRHydro")) { + 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")) { + oned_unslice(ny,nstart,nx,temp1d_minus,tempminus); + oned_unslice(ny,nstart,nx,temp1d_plus,tempplus); + } + } + + if(CCTK_Equals(Bvec_evolution_method,"GRHydro")) { + oned_unslice(ny,nstart,nx,Bvcx1d_minus,Bvecxminus); + oned_unslice(ny,nstart,nx,Bvcy1d_minus,Bvecyminus); + oned_unslice(ny,nstart,nx,Bvcz1d_minus,Bveczminus); + oned_unslice(ny,nstart,nx,Bvcx1d_plus,Bvecxplus); + oned_unslice(ny,nstart,nx,Bvcy1d_plus,Bvecyplus); + oned_unslice(ny,nstart,nx,Bvcz1d_plus,Bveczplus); + if(clean_divergence) { + oned_unslice(ny,nstart,nx,psidc1d_minus,psidcminus); + oned_unslice(ny,nstart,nx,psidc1d_plus,psidcplus); + } + } + for (i=0; i + + ppm1d_cxx + (ny,CCTK_DELTA_SPACE(1), + &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], + &ye1d_minus[0],&Bvcz1d_minus[0],&Bvcx1d_minus[0],&Bvcy1d_minus[0],&psidc1d_minus[0], + &rho1d_plus[0],&velz1d_plus[0],&velx1d_plus[0],&vely1d_plus[0],&eps1d_plus[0],&temp1d_plus[0], + &ye1d_plus[0],&Bvcz1d_plus[0],&Bvcx1d_plus[0],&Bvcy1d_plus[0],&psidc1d_plus[0]); + + + oned_unslice(nz,nstart,nxy,rho1d_minus,rhominus); + oned_unslice(nz,nstart,nxy,rho1d_plus,rhoplus); + + oned_unslice(nz,nstart,nxy,velx1d_minus,velxminus); + oned_unslice(nz,nstart,nxy,vely1d_minus,velyminus); + oned_unslice(nz,nstart,nxy,velz1d_minus,velzminus); + oned_unslice(nz,nstart,nxy,velx1d_plus,velxplus); + oned_unslice(nz,nstart,nxy,vely1d_plus,velyplus); + oned_unslice(nz,nstart,nxy,velz1d_plus,velzplus); + + oned_unslice(nz,nstart,nxy,eps1d_minus,epsminus); + oned_unslice(nz,nstart,nxy,eps1d_plus,epsplus); + + if(CCTK_Equals(Y_e_evolution_method,"GRHydro")) { + 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")) { + oned_unslice(nz,nstart,nxy,temp1d_minus,tempminus); + oned_unslice(nz,nstart,nxy,temp1d_plus,tempplus); + } + } + + if(CCTK_Equals(Bvec_evolution_method,"GRHydro")) { + oned_unslice(nz,nstart,nxy,Bvcx1d_minus,Bvecxminus); + oned_unslice(nz,nstart,nxy,Bvcy1d_minus,Bvecyminus); + oned_unslice(nz,nstart,nxy,Bvcz1d_minus,Bveczminus); + oned_unslice(nz,nstart,nxy,Bvcx1d_plus,Bvecxplus); + oned_unslice(nz,nstart,nxy,Bvcy1d_plus,Bvecyplus); + oned_unslice(nz,nstart,nxy,Bvcz1d_plus,Bveczplus); + if(clean_divergence) { + oned_unslice(nz,nstart,nxy,psidc1d_minus,psidcminus); + oned_unslice(nz,nstart,nxy,psidc1d_plus,psidcplus); + } + } + for (i=0; i +static void 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); -- cgit v1.2.3