Commits

commits@lists.einsteintoolkit.org

2 participants
6997 discussions

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 574)
by rhaas＠tapir.caltech.edu 13 Aug '13

13 Aug '13

User: rhaas Date: 2013/08/13 09:56 AM Added: /trunk/src/ GRHydro_PPMReconstruct_drv_opt.cc, GRHydro_PPM_opt.h Log: 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 <jfaberrit(a)jfaber3.local> File Changes: Directory: /trunk/src/ ====================== File [added]: GRHydro_PPMReconstruct_drv_opt.cc Delta lines: +521 -0 =================================================================== --- trunk/src/GRHydro_PPMReconstruct_drv_opt.cc (rev 0) +++ trunk/src/GRHydro_PPMReconstruct_drv_opt.cc 2013-08-13 14:56:31 UTC (rev 574) @@ -0,0 +1,521 @@ +#include <cmath> +#include <algorithm> +#include <vector> +#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<CCTK_REAL> vec1d) +{ + int i; + for (i=0; i<npt; i++)vec1d[i] = vec3d[nstart+i*nstride]; + return; +} + +void oned_unslice(int npt,int nstart,int nstride,vector<CCTK_REAL> vec1d,CCTK_REAL* vec3d) +{ + int i; + for (i=0; i<npt; i++)vec3d[nstart+i*nstride]=vec1d[i]; + return; +} + +/* + Cases that must be considered: + * basic hydro + * hydro + temperature + ye + * hydro + ye + * basic mhd + * mhd + temperature + ye + * mhd + ye + * mppm (not supported right now) + * not supporting trivial_rp + * with or without divergence cleaning + */ + +void GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + CCTK_REAL * restrict g11, * restrict g12, *restrict g13; + CCTK_REAL * restrict g22, * restrict g23, *restrict g33; + CCTK_REAL * restrict beta1, * restrict beta2, * restrict beta3; + CCTK_REAL * restrict vup; + + //Multipatch related pointers + if(GRHydro_UseGeneralCoordinates(cctkGH)) { + g11=gaa; + g12=gab; + g13=gac; + g22=gbb; + g23=gbc; + g33=gcc; + beta1=betaa; + beta2=betab; + beta3=betac; + vup=lvel; + } else { + g11=gxx; + g12=gxy; + g13=gxz; + g22=gyy; + g23=gyz; + g33=gzz; + beta1=betax; + beta2=betay; + beta3=betaz; + vup=vel; + } + + int nx=cctk_lsh[0]; + int ny=cctk_lsh[1]; + int nz=cctk_lsh[2]; + int nxy=nx*ny; + 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; + // UNUSED: int trivialx,trivialy,trivialz; + int not_trivialx,not_trivialy,not_trivialz; + + type_bitsx = SpaceMask_GetTypeBits("Hydro_RiemannProblemX"); + // trivialx = SpaceMask_GetStateBits("Hydro_RiemannProblemX","trivial"); + not_trivialx = SpaceMask_GetStateBits("Hydro_RiemannProblemX","not_trivial"); + + type_bitsy = SpaceMask_GetTypeBits("Hydro_RiemannProblemY"); + // trivialy = SpaceMask_GetStateBits("Hydro_RiemannProblemY","trivial"); + not_trivialy = SpaceMask_GetStateBits("Hydro_RiemannProblemY","not_trivial"); + + type_bitsz = SpaceMask_GetTypeBits("Hydro_RiemannProblemZ"); + // trivialz = SpaceMask_GetStateBits("Hydro_RiemannProblemZ","trivial"); + not_trivialz = SpaceMask_GetStateBits("Hydro_RiemannProblemZ","not_trivial"); + + // ENHANCED PPM IS NOT IMPLEMENTED + // logical :: apply_enhanced_ppm + // ! if use_enhanced_ppm, allow old PPM on one level + // if (GRHydro_oppm_reflevel .eq. (-1) .or. & + // GRHydro_reflevel .ne. GRHydro_oppm_reflevel) then + // apply_enhanced_ppm = use_enhanced_ppm .ne. 0 + // else + // apply_enhanced_ppm = .false. + // end if + + + 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 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<do_temp,do_ye,do_mhd,clean_divergence,ppm_detect> + + ppm1d_cxx<true,true,true,true,true> + (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<nx; i++) { + SpaceMask_SetStateBits(space_mask, i+j*nx+k*nx*ny, type_bitsx, not_trivialx); + } + } + } + // !$OMP END PARALLEL DO + + } else if (flux_direction[0]==2) { + + // !$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++) { + + nstart=j+k*nx*ny; + + 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,nx,eps,eps1d); + oned_slice(ny,nstart,nx,press,press1d); + + if(Y_e!=NULL) { + oned_slice(ny,nstart,nx,Y_e,ye1d); + if(temperature!=NULL){ + oned_slice(ny,nstart,nx,temperature,temp1d); + } + } + + 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); + 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"); + } + + //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<true,true,true,true,true> + (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<ny; i++) { + SpaceMask_SetStateBits(space_mask, j+i*nx+k*nx*ny, type_bitsy, not_trivialy); + } + + } + } + // !$OMP END PARALLEL DO + + } else if (flux_direction[0]==3) { + + // !$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++) { + + nstart=j+k*nx; + + 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,nxy,eps,eps1d); + oned_slice(nz,nstart,nxy,press,press1d); + + if(Y_e!=NULL) { + oned_slice(nz,nstart,nxy,Y_e,ye1d); + if(temperature!=NULL){ + oned_slice(nz,nstart,nxy,temperature,temp1d); + } + } + + 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); + 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"); + } + + //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<true,true,true,true,true> + (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<nz; i++) { + SpaceMask_SetStateBits(space_mask, j+k*nx+i*nx*ny, type_bitsz, not_trivialz); + } + + } + } + // !$OMP END PARALLEL DO + + + } +} + + +// 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 + File [added]: GRHydro_PPM_opt.h Delta lines: +38 -0 =================================================================== --- trunk/src/GRHydro_PPM_opt.h (rev 0) +++ trunk/src/GRHydro_PPM_opt.h 2013-08-13 14:56:31 UTC (rev 574) @@ -0,0 +1,38 @@ +template<bool do_temp, bool do_ye, bool do_mhd, + bool dc_flag, bool do_ppm_detect> +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);

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 573)
by rhaas＠tapir.caltech.edu 13 Aug '13

13 Aug '13

User: rhaas Date: 2013/08/13 09:56 AM Added: /trunk/src/ GRHydro_PPM.cc Modified: /trunk/src/ make.code.defn Log: GRHydro: add templated (yay\!) PPM routine; can do old ppm, but for MHD, hot, cold whatever From: Christian Ott <cott(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [added]: GRHydro_PPM.cc Delta lines: +414 -0 =================================================================== --- trunk/src/GRHydro_PPM.cc (rev 0) +++ trunk/src/GRHydro_PPM.cc 2013-08-13 14:56:28 UTC (rev 573) @@ -0,0 +1,414 @@ +#include <cmath> +#include <algorithm> +#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)) + +using namespace std; + +/* + Cases that must be considered: + * basic hydro + * hydro + temperature + ye + * hydro + ye + * basic mhd + * mhd + temperature + ye + * mhd + ye + * mppm (not supported right now) + * not supporting trivial_rp + * with or without divergence cleaning + */ + + +static inline void steep(const double *restrict x, + const double *restrict dx, + double *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]), + MIN(2.0*fabs(x[i]-x[i-1]), + 2.0*fabs(x[i+1]-x[i]))); + } else { + dmx[i] = 0.0; + } +} + +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, + const int i) { + + if ( !(xplus[i]==x[i] && x[i]==xminus[i]) + && ( (xplus[i]-x[i])*(x[i]-xminus[i]) <= 0.0 ) ) + { + xminus[i] = x[i]; + xplus[i] = x[i]; + } else if( 6.0 * (xplus[i]-xminus[i]) * + (x[i]-0.5*(xplus[i]+xminus[i])) > + (xplus[i]-xminus[i])*(xplus[i]-xminus[i]) ) + { + xminus[i] = 3.0*x[i]-2.0*xplus[i]; + } else if( 6.0 * (xplus[i]-xminus[i]) * + (x[i]-0.5*(xplus[i]+xminus[i])) < + -(xplus[i]-xminus[i])*(xplus[i]-xminus[i]) ) + { + xplus[i] = 3.0*x[i]-2.0*xminus[i]; + } + + return; +} + + + +template<bool do_temp, bool do_ye, bool do_mhd, + bool dc_flag, bool do_ppm_detect> +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) +{ + DECLARE_CCTK_PARAMETERS; + + double drho[nx], dpress[nx], deps[nx], d2rho[nx]; + double dmrho[nx], dmeps[nx]; + double dmtemp[nx], dmye[nx]; + double dtemp[nx], dye[nx]; + double dvelx[nx], dvely[nx], dvelz[nx]; + double dmvelx[nx], dmvely[nx], dmvelz[nx]; + double dBvcx[nx], dBvcy[nx], dBvcz[nx], dpsidc[nx]; + double dmBvcx[nx], dmBvcy[nx], dmBvcz[nx], dmpsidc[nx]; + double tilde_flatten[nx]; + + const double onesixth = 1.0/6.0; + + // 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++) { + 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]); + deps[i] = 0.5 * (eps[i+1]-eps[i-1]); + dvelx[i] = 0.5 * (velx[i+1]-velx[i-1]); + dvely[i] = 0.5 * (vely[i+1]-vely[i-1]); + dvelz[i] = 0.5 * (velz[i+1]-velz[i-1]); + if(do_ye) { + if(do_temp) { + dtemp[i] = 0.5 * (temp[i+1]-temp[i-1]); + } + dye[i] = 0.5 * (ye[i+1]-ye[i-1]); + } + if(do_mhd) { + dBvcx[i] = 0.5 * (Bvcx[i+1]-Bvcx[i-1]); + dBvcy[i] = 0.5 * (Bvcy[i+1]-Bvcy[i-1]); + dBvcz[i] = 0.5 * (Bvcz[i+1]-Bvcz[i-1]); + if(dc_flag) { + dpsidc[i] = 0.5 * (psidc[i+1]-psidc[i-1]); + } + } + } + + // Steepened slope. See (1.8) of Colella and Woodward, p.178 + for(int i=1;i<nx-1;i++) { + steep(rho,drho,dmrho,i); + steep(eps,deps,dmeps,i); + steep(velx,dvelx,dmvelx,i); + steep(vely,dvely,dmvely,i); + steep(velz,dvelz,dmvelz,i); + if(do_ye) { + if(do_temp) { + steep(temp,dtemp,dmtemp,i); + } + steep(ye,dye,dmye,i); + } + if(do_mhd) { + steep(Bvcx,dBvcx,dmBvcx,i); + steep(Bvcy,dBvcy,dmBvcy,i); + steep(Bvcz,dBvcz,dmBvcz,i); + if(dc_flag) { + steep(psidc,dpsidc,dmpsidc,i); + } + } + } + + // Initial boundary states. See (1.9) of Colella and Woodward, p.178 + 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]) + + (dmeps[i] - dmeps[i+1]) * onesixth; + velxplus[i] = 0.5 * (velx[i] + velx[i+1]) + + (dmvelx[i] - dmvelx[i+1]) * onesixth; + velyplus[i] = 0.5 * (vely[i] + vely[i+1]) + + (dmvely[i] - dmvely[i+1]) * onesixth; + velzplus[i] = 0.5 * (velz[i] + velz[i+1]) + + (dmvelz[i] - dmvelz[i+1]) * onesixth; + if(do_ye) { + if(do_temp) { + tempplus[i] = 0.5 * (temp[i] + temp[i+1]) + + (dmtemp[i] - dmtemp[i+1]) * onesixth; + } + yeplus[i] = 0.5 * (ye[i] + ye[i+1]) + + (dmye[i] - dmye[i+1]) * onesixth; + } + if(do_mhd) { + Bvcxplus[i] = 0.5 * (Bvcx[i] + Bvcx[i+1]) + + (dmBvcx[i] - dmBvcx[i+1]) * onesixth; + Bvcyplus[i] = 0.5 * (Bvcy[i] + Bvcy[i+1]) + + (dmBvcy[i] - dmBvcy[i+1]) * onesixth; + Bvczplus[i] = 0.5 * (Bvcz[i] + Bvcz[i+1]) + + (dmBvcz[i] - dmBvcz[i+1]) * onesixth; + if(dc_flag) { + psidcplus[i] = 0.5 * (psidc[i] + psidc[i+1]) + + (dmpsidc[i] - dmpsidc[i+1]) * onesixth; + } + } + } + + // fill minus states + for(int i=1;i<nx-2;i++) { + rhominus[i+1] = rhoplus[i]; + epsminus[i+1] = epsplus[i]; + velxminus[i+1] = velxplus[i]; + velyminus[i+1] = velyplus[i]; + velzminus[i+1] = velzplus[i]; + if(do_ye) { + if(do_temp) { + tempminus[i+1] = tempplus[i]; + } + yeminus[i+1] = yeplus[i]; + } + if(do_mhd) { + Bvcxminus[i+1] = Bvcxplus[i]; + Bvcyminus[i+1] = Bvcyplus[i]; + Bvczminus[i+1] = Bvczplus[i]; + if(dc_flag) { + psidcminus[i+1] = psidcplus[i]; + } + } + } + + /* + Discontinuity steepening. See (1.14-17) of C&W. + This is the detect routine which mat be + activated with the ppm_detect parameter + Note that this part really also depends on the grid being even. + Note also that we don''t have access to the gas constant gamma. + So this is just dropped from eq. (3.2) of C&W. + We can get around this by just rescaling the constant k0 (ppm_k0 here). + */ + + if(do_ppm_detect) { + 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 + * MIN(fabs(rho[i+1]), + fabs(rho[i-1])) > 0.0) ) + { + etatilde = (rho[i-2] - rho[i+2] + 4.0 * drho[i]) / (drho[i] * 12.0); + } + double eta = MAX(0.0, MIN(1.0, ppm_eta1 * (etatilde - ppm_eta2))); + if (ppm_k0 * fabs(drho[i]) * MIN(press[i-1],press[i+1]) + < fabs(dpress[i]) * MIN(rho[i-1], rho[i+1])) + { + eta = 0.0; + } + rhominus[i] = rhominus[i] * (1.0 - eta) + + (rho[i-1] + 0.5 * dmrho[i-1]) * eta; + rhoplus[i] = rhoplus[i] * (1.0 - eta) + + (rho[i+1] - 0.5 * dmrho[i+1]) * eta; + } + } + + // 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]; + double w=0.0; + if ( (fabs(dpress[i]) > ppm_epsilon * MIN(press[i-1],press[i+1])) + && (dvel < 0.0) ) + { + w = 1.0; + } + if (fabs(dpress2) < ppm_small) + { + tilde_flatten[i] = 1.0; + } + else + { + tilde_flatten[i] = MAX(0.0, 1.0 - w * MAX(0.0, + ppm_omega2 * (dpress[i] / dpress2 - ppm_omega1))); + } + } + + for(int i=2;i<nx-2;i++) { + 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]; + epsminus[i] = flatten * epsminus[i] + (1.0 - flatten) * eps[i]; + velxplus[i] = flatten * velxplus[i] + (1.0 - flatten) * velx[i]; + velxminus[i] = flatten * velxminus[i] + (1.0 - flatten) * velx[i]; + velyplus[i] = flatten * velyplus[i] + (1.0 - flatten) * vely[i]; + velyminus[i] = flatten * velyminus[i] + (1.0 - flatten) * vely[i]; + velzplus[i] = flatten * velzplus[i] + (1.0 - flatten) * velz[i]; + velzminus[i] = flatten * velzminus[i] + (1.0 - flatten) * velz[i]; + + if(do_ye) { + if(do_temp) { + tempplus[i] = flatten * tempplus[i] + (1.0 - flatten) * temp[i]; + tempminus[i] = flatten * tempminus[i] + (1.0 - flatten) * temp[i]; + } + yeplus[i] = flatten * yeplus[i] + (1.0 - flatten) * ye[i]; + yeminus[i] = flatten * yeminus[i] + (1.0 - flatten) * ye[i]; + } + + if(do_mhd) { + Bvcxplus[i] = flatten * Bvcxplus[i] + (1.0 - flatten) * Bvcx[i]; + Bvcxminus[i] = flatten * Bvcxminus[i] + (1.0 - flatten) * Bvcx[i]; + Bvcyplus[i] = flatten * Bvcyplus[i] + (1.0 - flatten) * Bvcy[i]; + Bvcyminus[i] = flatten * Bvcyminus[i] + (1.0 - flatten) * Bvcy[i]; + Bvczplus[i] = flatten * Bvczplus[i] + (1.0 - flatten) * Bvcz[i]; + Bvczminus[i] = flatten * Bvczminus[i] + (1.0 - flatten) * Bvcz[i]; + if(dc_flag) { + psidcplus[i] = flatten * psidcplus[i] + (1.0 - flatten) * psidc[i]; + psidcminus[i] = flatten * psidcminus[i] + (1.0 - flatten) * psidc[i]; + } + } + } // flattening + + for(int i=2;i<nx-2;i++) { + monotonize(rhominus,rho,rhoplus,i); + monotonize(epsminus,eps,epsplus,i); + monotonize(velxminus,velx,velxplus,i); + monotonize(velyminus,vely,velyplus,i); + monotonize(velzminus,velz,velzplus,i); + if(do_ye) { + if(do_temp) { + monotonize(tempminus,temp,tempplus,i); + } + monotonize(yeminus,ye,yeplus,i); + } + if(do_mhd) { + monotonize(Bvcxminus,Bvcx,Bvcxplus,i); + monotonize(Bvcyminus,Bvcy,Bvcyplus,i); + monotonize(Bvczminus,Bvcz,Bvczplus,i); + if(dc_flag) { + monotonize(psidcminus,psidc,psidcplus,i); + } + } + } + + return; +} + +//template template<bool do_temp, bool do_ye, bool do_mhd, bool dc_flag, bool do_ppm_detect> +#define INSTANTIATE(do_temp,do_ye,do_mhd,dc_flag,do_ppm_detect) \ +template static void ppm1d_cxx<do_temp,do_ye,do_mhd,dc_flag,do_ppm_detect>(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); + +//template template<bool do_temp, bool do_ye, bool do_mhd, bool dc_flag, bool do_ppm_detect> + + +// first stuff without MHD +INSTANTIATE(false,false,false,false,false) +INSTANTIATE(false,false,false,false,true) +// do_temp can only be true of do_ye is also true +INSTANTIATE(true,true,false,false,false) +INSTANTIATE(true,true,false,false,true) +// but do_ye can be true and do_temp can be false +INSTANTIATE(false,true,false,false,false) +INSTANTIATE(false,true,false,false,true) + +// with MHD, but dc_flag false +INSTANTIATE(false,false,true,false,false) +INSTANTIATE(false,false,true,false,true) +// do_temp can only be true of do_ye is also true +INSTANTIATE(true,true,true,false,false) +INSTANTIATE(true,true,true,false,true) +// but do_ye can be true and do_temp can be false +INSTANTIATE(false,true,true,false,false) +INSTANTIATE(false,true,true,false,true) + +// with MHD, dc_flag true +INSTANTIATE(false,false,true,true,false) +INSTANTIATE(false,false,true,true,true) +// do_temp can only be true of do_ye is also true +INSTANTIATE(true,true,true,true,false) +INSTANTIATE(true,true,true,true,true) +// but do_ye can be true and do_temp can be false +INSTANTIATE(false,true,true,true,false) +INSTANTIATE(false,true,true,true,true) File [modified]: make.code.defn Delta lines: +1 -0 =================================================================== --- trunk/src/make.code.defn 2013-08-13 14:56:24 UTC (rev 572) +++ trunk/src/make.code.defn 2013-08-13 14:56:28 UTC (rev 573) @@ -82,6 +82,7 @@ Con2Prim_fortran_interfaces.F90 \ Con2PrimM_fortran_interfaces.F90 \ GRHydro_Prim2Con.c \ + GRHydro_PPM.cc \ GRHydro_Wrappers.F90 # Subdirectories containing source files

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 572)
by rhaas＠tapir.caltech.edu 13 Aug '13

13 Aug '13

User: rhaas Date: 2013/08/13 09:56 AM Modified: /trunk/src/ GRHydro_PPM.c Log: GRHydro: in c PPM routine, remove mistake From: Christian D. Ott <cott(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_PPM.c Delta lines: +0 -8 =================================================================== --- trunk/src/GRHydro_PPM.c 2013-08-13 14:56:21 UTC (rev 571) +++ trunk/src/GRHydro_PPM.c 2013-08-13 14:56:24 UTC (rev 572) @@ -73,14 +73,6 @@ } for(int i=1;i<*nx-2;i++) { - rhoplus[i] = approx_at_cell_interface(rho,i); - epsplus[i] = approx_at_cell_interface(eps,i); - velxplus[i] = approx_at_cell_interface(velx,i); - velyplus[i] = approx_at_cell_interface(vely,i); - velzplus[i] = approx_at_cell_interface(velz,i); - } - - for(int i=1;i<*nx-2;i++) { rhominus[i+1] = rhoplus[i]; epsminus[i+1] = epsplus[i]; velxminus[i+1] = velxplus[i];

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 571)
by rhaas＠tapir.caltech.edu 13 Aug '13

13 Aug '13

User: rhaas Date: 2013/08/13 09:56 AM Modified: /trunk/src/ GRHydro_CalcUpdate.F90 Log: Move Evec calculation in seperate loops this avoids accessing unitinialized values in the ghost zones when computing the hydro varaible updates. This also saves time in running those loops since they are smaller. From: Philipp Moesta <pmoesta(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_CalcUpdate.F90 Delta lines: +40 -24 =================================================================== --- trunk/src/GRHydro_CalcUpdate.F90 2013-08-13 14:56:16 UTC (rev 570) +++ trunk/src/GRHydro_CalcUpdate.F90 2013-08-13 14:56:21 UTC (rev 571) @@ -47,10 +47,46 @@ if (use_weighted_fluxes == 0) then + if(evolve_mhd.ne.0 .and. transport_constraints.ne.0) then + !$OMP PARALLEL DO PRIVATE(i,j,k,alp_l,alp_r,alp_tmp) + do k = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(3) - GRHydro_stencil ! we need to compute Evec on all faces/edges where the fluxes are defined + do j = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(2) - GRHydro_stencil + do i = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(1) - GRHydro_stencil + + alp_l = 0.5d0 * (alp(i,j,k) + & + alp(i-xoffset,j-yoffset,k-zoffset)) + alp_r = 0.5d0 * (alp(i,j,k) + & + alp(i+xoffset,j+yoffset,k+zoffset)) + + ! we have to first compute all components of v\crossB = E and + ! combine them in the last substep into Bconshs + ! Evec lives on edges of cell: Evec(i,j,k,1) is at edge i,j+1/2,k+1/2 ie. the lower-front edge of cell (i,j,k) + if(flux_direction.eq.1) then + alp_tmp = 0.5d0 * (alp(i,j,k+1) + alp(i+xoffset,j+yoffset,k+zoffset+1)) + Evec(i,j,k,2) = Evec(i,j,k,2) + 0.25d0 * (alp_r*Bconszflux(i,j,k) + alp_tmp*Bconszflux(i ,j ,k+1)) + alp_tmp = 0.5d0 * (alp(i,j+1,k) + alp(i+xoffset,j+yoffset+1,k+zoffset)) + Evec(i,j,k,3) = Evec(i,j,k,3) - 0.25d0 * (alp_r*Bconsyflux(i,j,k) + alp_tmp*Bconsyflux(i ,j+1,k )) + elseif(flux_direction.eq.2) then + alp_tmp = 0.5d0 * (alp(i,j,k+1) + alp(i+xoffset,j+yoffset,k+zoffset+1)) + Evec(i,j,k,1) = Evec(i,j,k,1) - 0.25d0 * (alp_r*Bconszflux(i,j,k) + alp_tmp*Bconszflux(i ,j ,k+1)) + alp_tmp = 0.5d0 * (alp(i+1,j,k) + alp(i+xoffset+1,j+yoffset,k+zoffset)) + Evec(i,j,k,3) = Evec(i,j,k,3) + 0.25d0 * (alp_r*Bconsxflux(i,j,k) + alp_tmp*Bconsxflux(i+1,j ,k )) + elseif(flux_direction.eq.3) then + alp_tmp = 0.5d0 * (alp(i,j+1,k) + alp(i+xoffset,j+yoffset+1,k+zoffset)) + Evec(i,j,k,1) = Evec(i,j,k,1) + 0.25d0 * (alp_r*Bconsyflux(i,j,k) + alp_tmp*Bconsyflux(i ,j+1,k )) + alp_tmp = 0.5d0 * (alp(i+1,j,k) + alp(i+xoffset+1,j+yoffset,k+zoffset)) + Evec(i,j,k,2) = Evec(i,j,k,2) - 0.25d0 * (alp_r*Bconsxflux(i,j,k) + alp_tmp*Bconsxflux(i+1,j ,k )) + end if + enddo + enddo + enddo + !$OMP END PARALLEL DO + endif + !$OMP PARALLEL DO PRIVATE(i,j,k,itracer,alp_l,alp_r,alp_tmp,Bcons_l,Bcons_r) - do k = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(3) - GRHydro_stencil ! we need to compute Evec on all faces/edges where the fluxes are defined - do j = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(2) - GRHydro_stencil - do i = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(1) - GRHydro_stencil + do k = GRHydro_stencil + 1, cctk_lsh(3) - GRHydro_stencil ! we need to compute Evec on all faces/edges where the fluxes are defined + do j = GRHydro_stencil + 1, cctk_lsh(2) - GRHydro_stencil + do i = GRHydro_stencil + 1, cctk_lsh(1) - GRHydro_stencil alp_l = 0.5d0 * (alp(i,j,k) + & alp(i-xoffset,j-yoffset,k-zoffset)) @@ -74,27 +110,7 @@ alp_r * tauflux(i,j,k)) * idx if(evolve_mhd.ne.0) then - if(transport_constraints.ne.0) then - ! we have to first compute all components of v\crossB = E and - ! combine them in the last substep into Bconshs - ! Evec lives on edges of cell: Evec(i,j,k,1) is at edge i,j+1/2,k+1/2 ie. the lower-front edge of cell (i,j,k) - if(flux_direction.eq.1) then - alp_tmp = 0.5d0 * (alp(i,j,k+1) + alp(i+xoffset,j+yoffset,k+zoffset+1)) - Evec(i,j,k,2) = Evec(i,j,k,2) + 0.25d0 * (alp_r*Bconszflux(i,j,k) + alp_tmp*Bconszflux(i ,j ,k+1)) - alp_tmp = 0.5d0 * (alp(i,j+1,k) + alp(i+xoffset,j+yoffset+1,k+zoffset)) - Evec(i,j,k,3) = Evec(i,j,k,3) - 0.25d0 * (alp_r*Bconsyflux(i,j,k) + alp_tmp*Bconsyflux(i ,j+1,k )) - elseif(flux_direction.eq.2) then - alp_tmp = 0.5d0 * (alp(i,j,k+1) + alp(i+xoffset,j+yoffset,k+zoffset+1)) - Evec(i,j,k,1) = Evec(i,j,k,1) - 0.25d0 * (alp_r*Bconszflux(i,j,k) + alp_tmp*Bconszflux(i ,j ,k+1)) - alp_tmp = 0.5d0 * (alp(i+1,j,k) + alp(i+xoffset+1,j+yoffset,k+zoffset)) - Evec(i,j,k,3) = Evec(i,j,k,3) + 0.25d0 * (alp_r*Bconsxflux(i,j,k) + alp_tmp*Bconsxflux(i+1,j ,k )) - elseif(flux_direction.eq.3) then - alp_tmp = 0.5d0 * (alp(i,j+1,k) + alp(i+xoffset,j+yoffset+1,k+zoffset)) - Evec(i,j,k,1) = Evec(i,j,k,1) + 0.25d0 * (alp_r*Bconsyflux(i,j,k) + alp_tmp*Bconsyflux(i ,j+1,k )) - alp_tmp = 0.5d0 * (alp(i+1,j,k) + alp(i+xoffset+1,j+yoffset,k+zoffset)) - Evec(i,j,k,2) = Evec(i,j,k,2) - 0.25d0 * (alp_r*Bconsxflux(i,j,k) + alp_tmp*Bconsxflux(i+1,j ,k )) - end if - else + if(transport_constraints.ne.1) then Bconsrhs(i,j,k,1) = Bconsrhs(i,j,k,1) + & (alp_l * Bconsxflux(i-xoffset,j-yoffset,k-zoffset) - & alp_r * Bconsxflux(i,j,k)) * idx

1 0

[svn:einsteintoolkit] GRHydro/trunk/par/ (Rev. 570)
by rhaas＠tapir.caltech.edu 13 Aug '13

13 Aug '13

User: rhaas Date: 2013/08/13 09:56 AM Added: /trunk/par/ TOVstar_simple16d-mag.par Log: GRHydro: Adding parfile for vector potential TOV. From: Philipp Moesta <pmoesta(a)tapir.caltech.edu> File Changes: Directory: /trunk/par/ ====================== File [added]: TOVstar_simple16d-mag.par Delta lines: +305 -0 =================================================================== --- trunk/par/TOVstar_simple16d-mag.par (rev 0) +++ trunk/par/TOVstar_simple16d-mag.par 2013-08-13 14:56:16 UTC (rev 570) @@ -0,0 +1,305 @@ +## mclachlan tov_static +ActiveThorns = "Time MoL" +ActiveThorns = "Coordbase CartGrid3d Boundary StaticConformal" +ActiveThorns = "SymBase ADMBase TmunuBase HydroBase InitBase ADMCoupling ADMMacros" +ActiveThorns = "IOUtil" +ActiveThorns = "Formaline" +ActiveThorns = "SpaceMask CoordGauge Constants LocalReduce aeilocalinterp LoopControl" +ActiveThorns = "Carpet CarpetLib CarpetReduce CarpetRegrid2 CarpetInterp" +ActiveThorns = "CarpetIOASCII CarpetIOScalar CarpetIOHDF5 CarpetIOBasic" + +# grid parameters +CartGrid3D::type = "coordbase" +CartGrid3D::domain = "full" +CartGrid3D::avoid_origin = "no" +CoordBase::domainsize = "minmax" +CoordBase::xmin = 0.0 +CoordBase::ymin = 0.0 +CoordBase::zmin = 0.0 +CoordBase::xmax = 320.0 +CoordBase::ymax = 320.0 +CoordBase::zmax = 320.0 +CoordBase::dx = 8 +CoordBase::dy = 8 +CoordBase::dz = 8 + +CoordBase::boundary_size_x_lower = 3 +CoordBase::boundary_size_y_lower = 3 +CoordBase::boundary_size_z_lower = 3 +CoordBase::boundary_size_x_upper = 3 +CoordBase::boundary_size_y_upper = 3 +CoordBase::boundary_size_z_upper = 3 +CoordBase::boundary_shiftout_x_lower = 1 +CoordBase::boundary_shiftout_y_lower = 1 +CoordBase::boundary_shiftout_z_lower = 1 +CoordBase::boundary_shiftout_x_upper = 0 +CoordBase::boundary_shiftout_y_upper = 0 +CoordBase::boundary_shiftout_z_upper = 0 + +ActiveThorns = "Slab ReflectionSymmetry" +ActiveThorns = "RotatingSymmetry90" + +CarpetRegrid2::ensure_proper_nesting = yes +CarpetRegrid2::symmetry_rotating90 = "yes" + +ReflectionSymmetry::reflection_x = "no" +ReflectionSymmetry::reflection_y = "no" +ReflectionSymmetry::reflection_z = "yes" +ReflectionSymmetry::avoid_origin_x = "no" +ReflectionSymmetry::avoid_origin_y = "no" +ReflectionSymmetry::avoid_origin_z = "no" + +TmunuBase::stress_energy_storage = yes +TmunuBase::stress_energy_at_RHS = yes +TmunuBase::timelevels = 1 +TmunuBase::prolongation_type = none + +HydroBase::timelevels = 3 + +ADMMacros::spatial_order = 4 + +ADMBase::metric_type = "physical" + +SpaceMask::use_mask = "yes" + +Cactus::terminate = "time" +Cactus::cctk_final_time = 1400 + +Carpet::domain_from_coordbase = "yes" +Carpet::enable_all_storage = no +Carpet::use_buffer_zones = "yes" +Carpet::use_overlap_zones = "yes" + +Carpet::poison_new_timelevels = "yes" +Carpet::check_for_poison = "no" + +Carpet::init_3_timelevels = no +Carpet::init_fill_timelevels = "yes" + +CarpetLib::poison_new_memory = "yes" +CarpetLib::poison_value = 114 + +# system specific Carpet paramters +Carpet::max_refinement_levels = 10 +driver::ghost_size = 3 +Carpet::prolongation_order_space = 5 +Carpet::prolongation_order_time = 2 + +CarpetRegrid2::regrid_every = 0 +CarpetRegrid2::num_centres = 1 + +CarpetRegrid2::num_levels_1 = 5 +CarpetRegrid2::radius_1[1] =120.0 +CarpetRegrid2::radius_1[2] = 60.0 +CarpetRegrid2::radius_1[3] = 30.0 +CarpetRegrid2::radius_1[4] = 15.0 + + +time::dtfac = 0.25 +MoL::ODE_Method = "rk4" +MoL::MoL_Intermediate_Steps = 4 +MoL::MoL_Num_Scratch_Levels = 1 + +# check all physical variables for NaNs +ActiveThorns = "NaNChecker" +NaNChecker::check_every = 1 +NaNChecker::action_if_found = "terminate" #"terminate", "just warn", "abort" +NaNChecker::check_vars = "ADMBase::metric ADMBase::lapse ADMBase::shift HydroBase::rho HydroBase::eps HydroBase::press HydroBase::vel HydroBase::Bvec" + +## Lapse Condition: \partial_t alpha = - alpha K +## Shift Condition: \partial_t beta^i = 0 + +# Hydro paramters +ActiveThorns = "EOS_Omni" +ActiveThorns = "GRHydro GRHydro_InitData" + + +HydroBase::evolution_method = "GRHydro" +HydroBase::Bvec_evolution_method = "GRHydro_Avec" + +GRHydro::method_type = "RSA FV" +GRHydro::riemann_solver = "HLLE" +#GRHydro::GRHydro_eos_type = "Polytype" +#GRHydro::GRHydro_eos_table = "2D_Polytrope" +GRHydro::recon_method = "ppm" +GRHydro::recon_vars = "primitive" +GRHydro::GRHydro_stencil = 3 +#GRHydro::bound = "none" +GRHydro::bound = "flat" +GRHydro::rho_abs_min = 1.e-10 +#GRHydro::GRHydro = 18 # Tmunu(10), rho,press,eps,w_lorentz,vel, tau +#GRHydro::GRHydro = 10 # gij(6), alpha, beta(3) + +GRHydro::GRHydro_eos_type = "General" +GRHydro::GRHydro_eos_table = "Ideal_Fluid" + +GRHydro::Grhydro_MaxNumConstrainedVars = 33 +GRHydro::GRHydro_MaxNumEvolvedVars = 10 + +#GRHydro::clean_divergence = "yes" +#GRHydro::transport_constraints = "yes" +#GRHydro::kap_dc = 1.0 +GRHydro::track_divB = "yes" +GRHydro::calculate_bcom = "yes" + +GRHydro::c2p_reset_pressure = "yes" +#GRHydro::GRHydro_atmo_tolerance = 1 + +grhydro::grhydro_c2p_warnlevel = 0 +grhydro::grhydro_c2p_warn_from_reflevel = 2 +grhydro::con2prim_oct_hack = yes + +ActiveThorns = "GenericFD NewRad" +ActiveThorns = "ML_BSSN ML_BSSN_Helper" +ADMBase::evolution_method = "ML_BSSN" +ADMBase::lapse_evolution_method = "ML_BSSN" +ADMBase::shift_evolution_method = "ML_BSSN" +ADMBase::dtlapse_evolution_method= "ML_BSSN" +ADMBase::dtshift_evolution_method= "ML_BSSN" + +TmunuBase::support_old_CalcTmunu_mechanism = "no" + +ML_BSSN::timelevels = 3 + +ML_BSSN::harmonicN = 1 # 1+log +ML_BSSN::harmonicF = 1.0 # 1+log +ML_BSSN::LapseACoeff = 1.0 +ML_BSSN::ShiftBCoeff = 1.0 +ML_BSSN::ShiftGammaCoeff = 0.0 +ML_BSSN::AlphaDriver = 0.0 +ML_BSSN::BetaDriver = 0.0 +ML_BSSN::LapseAdvectionCoeff = 0.0 +ML_BSSN::ShiftAdvectionCoeff = 0.0 + +ML_BSSN::MinimumLapse = 1.0e-8 + +ML_BSSN::my_initial_boundary_condition = "extrapolate-gammas" +ML_BSSN::my_rhs_boundary_condition = "NewRad" + +ML_BSSN::ML_log_confac_bound = "none" +ML_BSSN::ML_metric_bound = "none" +ML_BSSN::ML_Gamma_bound = "none" +ML_BSSN::ML_trace_curv_bound = "none" +ML_BSSN::ML_curv_bound = "none" +ML_BSSN::ML_lapse_bound = "none" +ML_BSSN::ML_dtlapse_bound = "none" +ML_BSSN::ML_shift_bound = "none" +ML_BSSN::ML_dtshift_bound = "none" + +ActiveThorns = "SphericalSurface Dissipation" + +SphericalSurface::nsurfaces = 2 +SphericalSurface::maxntheta = 39 +SphericalSurface::maxnphi = 76 + +SphericalSurface::ntheta [0] = 39 +SphericalSurface::nphi [0] = 76 +SphericalSurface::nghoststheta[0] = 2 +SphericalSurface::nghostsphi [0] = 2 + +SphericalSurface::ntheta [1] = 39 +SphericalSurface::nphi [1] = 76 +SphericalSurface::nghoststheta[1] = 2 +SphericalSurface::nghostsphi [1] = 2 + +Dissipation::verbose = "no" +Dissipation::epsdis = 0.1 +Dissipation::vars = " + ML_BSSN::ML_log_confac + ML_BSSN::ML_metric + ML_BSSN::ML_curv + ML_BSSN::ML_trace_curv + ML_BSSN::ML_Gamma + ML_BSSN::ML_lapse + ML_BSSN::ML_shift + HydroBase::Avec +" + + +# init parameters +InitBase::initial_data_setup_method = "init_some_levels" + +ActiveThorns = "TOVSolver" + +ADMBase::initial_data = "tov" +ADMBase::initial_lapse = "tov" +ADMBase::initial_shift = "tov" +ADMBase::initial_dtlapse = "zero" +ADMBase::initial_dtshift = "zero" + +TOVSolver::TOV_Rho_Central[0] = 1.28e-3 +#TOVSolver::TOV_Mass[0] = 1.4 +TOVSolver::TOV_Gamma[0] = 2.0 +TOVSolver::TOV_K[0] = 100 + +HydroBase::initial_hydro = "tov" +HydroBase::initial_Bvec = "poloidalmagfield" +HydroBase::initial_Avec = "poloidalmagfield" +GRHydro::Avec_gauge = "algebraic" +GRHydro_InitData::poloidal_A_b = 0.1 +GRHydro_InitData::poloidal_n_p = 1 +GRHydro_InitData::poloidal_P_cut = 1.0e-6 +GRHydro_InitData::poloidal_rho_max = 1.221209344680609e-03 + +EOS_Omni::poly_K = 1.019575500488281e+02 +EOS_Omni::poly_gamma = 2.0 +EOS_Omni::gl_K = 1.019575500488281e+02 +EOS_Omni::gl_gamma = 2.0 +GRHydro::initial_Gamma = 2.0 + +# I/O + +IO::out_dir = $parfile + +IOBasic::outInfo_every = 8 +IOBasic::outInfo_vars = "Carpet::physical_time_per_hour HydroBase::rho ML_BSSN::ML_Ham HydroBase::bvec GRHydro::divb" + +IOScalar::outScalar_every = 8 +IOScalar::one_file_per_group = yes +IOScalar::outScalar_vars = " + HydroBase::rho + HydroBase::press + HydroBase::eps + HydroBase::vel + ADMBase::lapse + ADMBase::metric + ADMBase::curv + ML_BSSN::ML_Ham + ML_BSSN::ML_mom + GRHydro::dens +" + +IOASCII::out1D_every = 64 +IOASCII::one_file_per_group = yes +IOASCII::output_symmetry_points = no +IOASCII::out3D_ghosts = no +IOASCII::out3D_outer_ghosts = no +IOASCII::out1D_vars = " + HydroBase::rho + HydroBase::press + HydroBase::eps + HydroBase::vel + HydroBase::Bvec + ADMBase::lapse + ADMBase::metric + ADMBase::curv + ML_BSSN::ML_Ham + ML_BSSN::ML_mom + GRHydro::divb + GRHydro::bcons + GRHydro::bcom + GRHydro::bcom0 + GRHydro::bcom_sq +" + +IOHDF5::checkpoint = "yes" +IO::checkpoint_dir = $parfile +IO::checkpoint_ID = "yes" +IO::checkpoint_every = 2048 +#IO::checkpoint_every_walltime_hours = 6.0 +IO::checkpoint_keep = 3 +IO::checkpoint_on_terminate = "yes" + +IO::recover = "no" +IO::recover_dir = $parfile +

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 569)
by rhaas＠tapir.caltech.edu 13 Aug '13

13 Aug '13

User: rhaas Date: 2013/08/13 09:56 AM Modified: /trunk/ interface.ccl Log: add interface use OMP version of ESO_Omni_press function From: Christian Ott <cott(a)tapir.caltech.edu> File Changes: Directory: /trunk/ ================== File [modified]: interface.ccl Delta lines: +14 -0 =================================================================== --- trunk/interface.ccl 2013-08-13 14:56:11 UTC (rev 568) +++ trunk/interface.ccl 2013-08-13 14:56:13 UTC (rev 569) @@ -258,6 +258,20 @@ USES FUNCTION EOS_Omni_GetHandle +void FUNCTION EOS_Omni_pressOMP(CCTK_INT IN eoskey, \ + CCTK_INT IN havetemp, \ + CCTK_REAL IN rf_precision, \ + CCTK_INT IN npoints, \ + CCTK_REAL IN ARRAY rho, \ + CCTK_REAL INOUT ARRAY eps, \ + CCTK_REAL INOUT ARRAY temp, \ + CCTK_REAL IN ARRAY ye, \ + CCTK_REAL OUT ARRAY press, \ + CCTK_INT OUT ARRAY keyerr, \ + CCTK_INT OUT anyerr) + +USES FUNCTION EOS_Omni_pressOMP + void FUNCTION EOS_Omni_press(CCTK_INT IN eoskey, \ CCTK_INT IN havetemp, \ CCTK_REAL IN rf_precision, \

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 568)
by rhaas＠tapir.caltech.edu 13 Aug '13

13 Aug '13

User: rhaas Date: 2013/08/13 09:56 AM Modified: /trunk/src/ GRHydro_Prim2Con.c Log: GRHydro: avoid allocating large arrays on the stack From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Prim2Con.c Delta lines: +12 -8 =================================================================== --- trunk/src/GRHydro_Prim2Con.c 2013-08-13 14:56:08 UTC (rev 567) +++ trunk/src/GRHydro_Prim2Con.c 2013-08-13 14:56:11 UTC (rev 568) @@ -69,8 +69,9 @@ // EOS calls (now GF-wide) if(!*evolve_temper) { - int n = cctk_lsh[0]*cctk_lsh[1]*cctk_lsh[2]; - int keyerr[n]; int anyerr = 0; + int n = cctk_ash[0]*cctk_ash[1]*cctk_ash[2]; + int *keyerr = malloc(sizeof(*keyerr)*n); + int anyerr = 0; int keytemp = 0; // don't need special error handling for analytic EOS @@ -80,12 +81,12 @@ EOS_Omni_pressOMP(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n, rhoplus,epsplus,NULL,NULL,pressplus,keyerr,&anyerr); + free(keyerr); } else { if(reconstruct_temper) { - int n = cctk_lsh[0]*cctk_lsh[1]*cctk_lsh[2]; - // Linux usually offers ~8MB of stack which is good for ~1M points. - // We should consider malloc() if we ever smash the stack. - int keyerr[n]; int anyerr = 0; + int n = cctk_ash[0]*cctk_ash[1]*cctk_ash[2]; + int *keyerr = malloc(sizeof(*keyerr)*n); + int anyerr = 0; int keytemp = 1; // ensure Y_e and temperature within bounds @@ -140,10 +141,12 @@ CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, "Aborting!"); } + free(keyerr); } else { // ******************** EPS RECONSTRUCTION BRANCH ****************** - int n = cctk_lsh[0]*cctk_lsh[1]*cctk_lsh[2]; - int keyerr[n]; int anyerr = 0; + int n = cctk_ash[0]*cctk_ash[1]*cctk_ash[2]; + int *keyerr = malloc(sizeof(*keyerr)*n); + int anyerr = 0; int keytemp = 0; // ensure Y_e and temperature within bounds @@ -236,6 +239,7 @@ } } // for i, i<n error checking } + free(keyerr); } // end branch for no temp reconsturction } // end of evolve temper branch

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 567)
by rhaas＠tapir.caltech.edu 13 Aug '13

13 Aug '13

User: rhaas Date: 2013/08/13 09:56 AM Modified: /trunk/src/ GRHydro_Con2PrimHot.F90 Log: expand warnings in GRHydro_Con2PrimHot * compute reflevel in GRHydro_Con2PrimHot to be sure it is current * only output warnings if we want them From: Christian Ott <cott(a)tapir.caltech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Con2PrimHot.F90 Delta lines: +20 -16 =================================================================== --- trunk/src/GRHydro_Con2PrimHot.F90 2013-08-13 14:56:05 UTC (rev 566) +++ trunk/src/GRHydro_Con2PrimHot.F90 2013-08-13 14:56:08 UTC (rev 567) @@ -28,6 +28,7 @@ CCTK_REAL :: local_min_tracer CCTK_REAL :: local_perc_ptol + integer :: reflevel ! save memory when MP is not used CCTK_INT :: GRHydro_UseGeneralCoordinates @@ -42,6 +43,7 @@ ! end EOS Omni vars myproc = CCTK_MyProc(cctkGH) + reflevel = int(log10(dble(cctk_levfac(1)))/log10(2.0d0)) ! save memory when MP is not used if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then @@ -155,7 +157,7 @@ press(i,j,k),w_lorentz(i,j,k), & uxx,uxy,uxz,uyy,uyz,uzz,sdet,x(i,j,k),y(i,j,k), & z(i,j,k),r(i,j,k),epsnegative,GRHydro_rho_min,pmin, epsmin, & - GRHydro_reflevel, GRHydro_C2P_failed(i,j,k), GRHydro_perc_ptol) + reflevel, GRHydro_C2P_failed(i,j,k), GRHydro_perc_ptol) if(temperature(i,j,k).gt.GRHydro_max_temp) then !$OMP CRITICAL @@ -172,7 +174,7 @@ call CCTK_WARN(1,warnline) write(warnline,"(A7,i8)") "code: ",keyerr(1) call CCTK_WARN(1,warnline) - write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel + write(warnline,"(A10,i5)") "reflevel: ", reflevel call CCTK_WARN(1,warnline) call CCTK_WARN(0,"Aborting!!!") !$OMP END CRITICAL @@ -200,12 +202,12 @@ temperature(i,j,k),y_e(i,j,k),press(i,j,k),w_lorentz(i,j,k), & uxx,uxy,uxz,uyy,uyz,uzz,sdet,x(i,j,k),y(i,j,k), & z(i,j,k),r(i,j,k),epsnegative,GRHydro_rho_min,pmin, epsmin, & - GRHydro_reflevel, GRHydro_C2P_failed(i,j,k), local_perc_ptol) + reflevel, GRHydro_C2P_failed(i,j,k), local_perc_ptol) if( abs(GRHydro_C2P_failed(i,j,k)-2.0d0) .lt. 1.0d-10) then !$OMP CRITICAL - if (GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then + if (reflevel.ge.GRHydro_c2p_warn_from_reflevel) then call CCTK_WARN(1,"Convergence problem in c2p") - write(warnline,"(A10,i5)") "reflevel: ",GRHydro_reflevel + write(warnline,"(A10,i5)") "reflevel: ",reflevel call CCTK_WARN(1,warnline) write(warnline,"(A10,i10)") "iteration: ",cctk_iteration call CCTK_WARN(1,warnline) @@ -228,17 +230,19 @@ ! the temperature of the atmosphere. ! Now reset this point to minimum temperature. GRHydro_C2P_failed(i,j,k) = 0 - !$OMP CRITICAL - write(warnline,"(A10,i7,4i5,1P10E15.6)") "reset T:",& - cctk_iteration,GRHydro_Reflevel,& - i,j,k,rho(i,j,k),& - eps(i,j,k),temperature(i,j,k),y_e(i,j,k) - call CCTK_WARN(1,warnline) - write(warnline,"(A10,i7,4i5,1P10E15.6)") "reset T:",& - cctk_iteration,GRHydro_Reflevel,& - i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),r(i,j,k) - call CCTK_WARN(1,warnline) - !$OMP END CRITICAL + if (reflevel.ge.GRHydro_c2p_warn_from_reflevel) then + !$OMP CRITICAL + write(warnline,"(A10,i7,4i5,1P10E15.6)") "reset T:",& + cctk_iteration,GRHydro_Reflevel,& + i,j,k,rho(i,j,k),& + eps(i,j,k),temperature(i,j,k),y_e(i,j,k) + call CCTK_WARN(1,warnline) + write(warnline,"(A10,i7,4i5,1P10E15.6)") "reset T:",& + cctk_iteration,GRHydro_Reflevel,& + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),r(i,j,k) + call CCTK_WARN(1,warnline) + !$OMP END CRITICAL + endif temperature(i,j,k) = GRHydro_hot_atmo_temp keytemp = 1 call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 566)
by rhaas＠tapir.caltech.edu 13 Aug '13

13 Aug '13

User: rhaas Date: 2013/08/13 09:56 AM Modified: /trunk/src/ GRHydro_Boundaries.F90 Log: don't specify boundaries twice From: Frank Loeffler <knarf(a)cct.lsu.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Boundaries.F90 Delta lines: +0 -7 =================================================================== --- trunk/src/GRHydro_Boundaries.F90 2013-08-13 14:56:03 UTC (rev 565) +++ trunk/src/GRHydro_Boundaries.F90 2013-08-13 14:56:05 UTC (rev 566) @@ -235,13 +235,6 @@ if (CCTK_EQUALS(Bvec_evolution_method, "GRHydro_Avec")) then ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "HydroBase::Avec", "Flat") - if (general_coordinates .ne. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::lBvec", "Flat") - else - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::Bvec", "Flat") - endif endif if(evolve_tracer.ne.0) then

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 565)
by rhaas＠tapir.caltech.edu 13 Aug '13

13 Aug '13

User: rhaas Date: 2013/08/13 09:56 AM Modified: /trunk/ schedule.ccl Log: GRHydro: explicitly recompute SqrtDetg after RECOVERY this is usually not required, unless MoL_run_MoL_PostStep_in_PostRecoverVariables is set to false. From: Roland Haas <rhaas(a)tapir.caltech.edu> File Changes: Directory: /trunk/ ================== File [modified]: schedule.ccl Delta lines: +5 -0 =================================================================== --- trunk/schedule.ccl 2013-08-13 14:56:00 UTC (rev 564) +++ trunk/schedule.ccl 2013-08-13 14:56:03 UTC (rev 565) @@ -1084,7 +1084,12 @@ LANG: Fortran } "Calculate sdetg" +schedule GRHydro_SqrtSpatialDeterminant AT CCTK_POST_RECOVER_VARIABLES +{ + LANG: Fortran +} "Calculate sdetg" + #debug ###schedule GRHydro_Debug IN HydroBase_PostStep ###{

1 0

Jump to page:

2026

2025

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

Commits