[Commits] [svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 633)
rhaas at tapir.caltech.edu
rhaas at tapir.caltech.edu
Tue Apr 15 14:49:50 CDT 2014
User: rhaas
Date: 2014/04/15 02:49 PM
Modified:
/trunk/src/
GRHydro_Prim2Con.c, GRHydro_Prim2ConM.F90
Log:
GRHydro: make routine to compute speed of sound callable from F90 MHD p2c
the code is long and complicated enough that it should better not be
copied and pasted.
File Changes:
Directory: /trunk/src/
======================
File [modified]: GRHydro_Prim2Con.c
Delta lines: +77 -62
===================================================================
--- trunk/src/GRHydro_Prim2Con.c 2014-04-15 19:49:47 UTC (rev 632)
+++ trunk/src/GRHydro_Prim2Con.c 2014-04-15 19:49:50 UTC (rev 633)
@@ -25,17 +25,21 @@
const double rho, const double vx,
const double vy, const double vz,
const double eps, const double press);
+static inline void GRHydro_SpeedOfSound(CCTK_ARGUMENTS);
-
void Primitive2ConservativeC(CCTK_ARGUMENTS);
CCTK_FCALL void CCTK_FNAME(Primitive2ConservativeCforF)(cGH ** p_cctkGH) {
Primitive2ConservativeC(*p_cctkGH);
}
+CCTK_FCALL void CCTK_FNAME(GRHydro_SpeedOfSound)(cGH ** p_cctkGH) {
+ GRHydro_SpeedOfSound(*p_cctkGH);
+}
+
void Primitive2ConservativeC(CCTK_ARGUMENTS)
{
DECLARE_CCTK_ARGUMENTS;
@@ -61,12 +65,73 @@
g33 = gzz;
}
+ GRHydro_SpeedOfSound(CCTK_PASS_CTOC);
+
// if padding is used the the "extra" vector elements could contain junk
// which we do not know how to handle
assert(cctk_lsh[0] == cctk_ash[0]);
assert(cctk_lsh[1] == cctk_ash[1]);
assert(cctk_lsh[2] == cctk_ash[2]);
+#pragma omp parallel for
+ for(int k = GRHydro_stencil-1; k < cctk_lsh[2]-GRHydro_stencil+1; k++)
+ for(int j = GRHydro_stencil-1; j < cctk_lsh[1]-GRHydro_stencil+1; j++)
+#pragma ivdep // force compiler to vectorize the goddamn loop
+ for(int i = GRHydro_stencil-1; i < cctk_lsh[0]-GRHydro_stencil+1; i++) {
+
+ const int idx = CCTK_GFINDEX3D(cctkGH,i,j,k);
+
+ const int idxl = CCTK_GFINDEX3D(cctkGH,i-*xoffset,j-*yoffset,k-*zoffset);
+ const int idxr = CCTK_GFINDEX3D(cctkGH,i+*xoffset,j+*yoffset,k+*zoffset);
+
+ const double g11l = 0.5 * (g11[idx] + g11[idxl]);
+ const double g12l = 0.5 * (g12[idx] + g12[idxl]);
+ const double g13l = 0.5 * (g13[idx] + g13[idxl]);
+ const double g22l = 0.5 * (g22[idx] + g22[idxl]);
+ const double g23l = 0.5 * (g23[idx] + g23[idxl]);
+ const double g33l = 0.5 * (g33[idx] + g33[idxl]);
+
+ const double g11r = 0.5 * (g11[idx] + g11[idxr]);
+ const double g12r = 0.5 * (g12[idx] + g12[idxr]);
+ const double g13r = 0.5 * (g13[idx] + g13[idxr]);
+ const double g22r = 0.5 * (g22[idx] + g22[idxr]);
+ const double g23r = 0.5 * (g23[idx] + g23[idxr]);
+ const double g33r = 0.5 * (g33[idx] + g33[idxr]);
+
+ const double savg_detl =
+ sqrt(SpatialDeterminantC(g11l,g12l,g13l,g22l,g23l,g33l));
+ const double savg_detr =
+ sqrt(SpatialDeterminantC(g11r,g12r,g13r,g22r,g23r,g33r));
+
+ // minus call to p2c
+ prim2conC(&w_lorentzminus[idx], &densminus[idx], &sxminus[idx],
+ &syminus[idx], &szminus[idx], &tauminus[idx],
+ g11l,g12l,g13l,g22l,g23l,g33l,
+ savg_detl,rhominus[idx], velxminus[idx], velyminus[idx],
+ velzminus[idx], epsminus[idx], pressminus[idx]);
+
+
+ // plus call to p2c
+ prim2conC(&w_lorentzplus[idx], &densplus[idx], &sxplus[idx],
+ &syplus[idx], &szplus[idx], &tauplus[idx],
+ g11r,g12r,g13r,g22r,g23r,g33r,
+ savg_detr,rhoplus[idx], velxplus[idx], velyplus[idx],
+ velzplus[idx], epsplus[idx], pressplus[idx]);
+ }
+
+} // end function Conservative2PrimitiveC
+
+static inline void GRHydro_SpeedOfSound(CCTK_ARGUMENTS)
+{
+ DECLARE_CCTK_PARAMETERS;
+ DECLARE_CCTK_ARGUMENTS;
+
+ // if padding is used the the "extra" vector elements could contain junk
+ // which we do not know how to handle
+ assert(cctk_lsh[0] == cctk_ash[0]);
+ assert(cctk_lsh[1] == cctk_ash[1]);
+ assert(cctk_lsh[2] == cctk_ash[2]);
+
// EOS calls (now GF-wide)
if(!*evolve_temper) {
// n needs to be computed using ash since ash is used when computing the
@@ -92,7 +157,7 @@
} else {
if(reconstruct_temper) {
int n = cctk_ash[0]*cctk_ash[1]*cctk_ash[2];
- int nx = cctk_lsh[0] - (GRHydro_stencil - 1) - (GRHydro_stencil) + 1;
+ int nx = cctk_lsh[0] - (GRHydro_stencil - 1) - (GRHydro_stencil) + 1 + (transport_constraints && !*xoffset);
int *keyerr = malloc(sizeof(*keyerr)*n);
int keytemp = 1;
// ensure Y_e and temperature within bounds
@@ -107,8 +172,8 @@
// call the EOS with slices
#pragma omp parallel for
- for(int k=GRHydro_stencil-1;k<cctk_lsh[2]-GRHydro_stencil+1;k++)
- for(int j=GRHydro_stencil-1;j<cctk_lsh[1]-GRHydro_stencil+1;j++) {
+ for(int k=GRHydro_stencil-1;k<cctk_lsh[2]-GRHydro_stencil+1 + (transport_constraints && !*zoffset);k++)
+ for(int j=GRHydro_stencil-1;j<cctk_lsh[1]-GRHydro_stencil+1 + (transport_constraints && !*yoffset);j++) {
int anyerr = 0;
int i = CCTK_GFINDEX3D(cctkGH,GRHydro_stencil-1,j,k);
EOS_Omni_press_cs2(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,nx,
@@ -136,8 +201,8 @@
#pragma omp parallel for
- for(int k=GRHydro_stencil-1;k<cctk_lsh[2]-GRHydro_stencil+1;k++)
- for(int j=GRHydro_stencil-1;j<cctk_lsh[1]-GRHydro_stencil+1;j++) {
+ for(int k=GRHydro_stencil-1;k<cctk_lsh[2]-GRHydro_stencil+1 + (transport_constraints && !*zoffset);k++)
+ for(int j=GRHydro_stencil-1;j<cctk_lsh[1]-GRHydro_stencil+1 + (transport_constraints && !*yoffset);j++) {
int anyerr = 0;
int i = CCTK_GFINDEX3D(cctkGH,GRHydro_stencil-1,j,k);
EOS_Omni_press_cs2(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,nx,
@@ -166,7 +231,7 @@
} else {
// ******************** EPS RECONSTRUCTION BRANCH ******************
int n = cctk_ash[0]*cctk_ash[1]*cctk_ash[2];
- int nx = cctk_lsh[0] - (GRHydro_stencil - 1) - (GRHydro_stencil) + 1;
+ int nx = cctk_lsh[0] - (GRHydro_stencil - 1) - (GRHydro_stencil) + 1 + (transport_constraints && !*xoffset);
int *keyerr = malloc(sizeof(*keyerr)*n);
int keytemp = 0;
@@ -183,8 +248,8 @@
// call the EOS with slices for minus states
#pragma omp parallel for
- for(int k=GRHydro_stencil-1;k<cctk_lsh[2]-GRHydro_stencil+1;k++)
- for(int j=GRHydro_stencil-1;j<cctk_lsh[1]-GRHydro_stencil+1;j++) {
+ for(int k=GRHydro_stencil-1;k<cctk_lsh[2]-GRHydro_stencil+1 + (transport_constraints && !*zoffset);k++)
+ for(int j=GRHydro_stencil-1;j<cctk_lsh[1]-GRHydro_stencil+1 + (transport_constraints && !*yoffset);j++) {
int i = CCTK_GFINDEX3D(cctkGH,GRHydro_stencil-1,j,k);
int anyerr = 0;
EOS_Omni_press_cs2(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,nx,
@@ -233,8 +298,8 @@
// call the EOS with slices for plus states
#pragma omp parallel for
- for(int k=GRHydro_stencil-1;k<cctk_lsh[2]-GRHydro_stencil+1;k++)
- for(int j=GRHydro_stencil-1;j<cctk_lsh[1]-GRHydro_stencil+1;j++) {
+ for(int k=GRHydro_stencil-1;k<cctk_lsh[2]-GRHydro_stencil+1 + (transport_constraints && !*zoffset);k++)
+ for(int j=GRHydro_stencil-1;j<cctk_lsh[1]-GRHydro_stencil+1 + (transport_constraints && !*yoffset);j++) {
int i = CCTK_GFINDEX3D(cctkGH,GRHydro_stencil-1,j,k);
int anyerr = 0;
EOS_Omni_press_cs2(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,nx,
@@ -278,64 +343,14 @@
} // omp critical pragma
} // if keyerr
} // loop ii
- CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
- "Aborting!");
} // if (anyerr)
} // big loop over plus states
free(keyerr);
} // end branch for no temp reconsturction
} // end of evolve temper branch
+}
-
-#pragma omp parallel for
- for(int k = GRHydro_stencil-1; k < cctk_lsh[2]-GRHydro_stencil+1; k++)
- for(int j = GRHydro_stencil-1; j < cctk_lsh[1]-GRHydro_stencil+1; j++)
-#pragma ivdep // force compiler to vectorize the goddamn loop
- for(int i = GRHydro_stencil-1; i < cctk_lsh[0]-GRHydro_stencil+1; i++) {
-
- const int idx = CCTK_GFINDEX3D(cctkGH,i,j,k);
-
- const int idxl = CCTK_GFINDEX3D(cctkGH,i-*xoffset,j-*yoffset,k-*zoffset);
- const int idxr = CCTK_GFINDEX3D(cctkGH,i+*xoffset,j+*yoffset,k+*zoffset);
-
- const double g11l = 0.5 * (g11[idx] + g11[idxl]);
- const double g12l = 0.5 * (g12[idx] + g12[idxl]);
- const double g13l = 0.5 * (g13[idx] + g13[idxl]);
- const double g22l = 0.5 * (g22[idx] + g22[idxl]);
- const double g23l = 0.5 * (g23[idx] + g23[idxl]);
- const double g33l = 0.5 * (g33[idx] + g33[idxl]);
-
- const double g11r = 0.5 * (g11[idx] + g11[idxr]);
- const double g12r = 0.5 * (g12[idx] + g12[idxr]);
- const double g13r = 0.5 * (g13[idx] + g13[idxr]);
- const double g22r = 0.5 * (g22[idx] + g22[idxr]);
- const double g23r = 0.5 * (g23[idx] + g23[idxr]);
- const double g33r = 0.5 * (g33[idx] + g33[idxr]);
-
- const double savg_detl =
- sqrt(SpatialDeterminantC(g11l,g12l,g13l,g22l,g23l,g33l));
- const double savg_detr =
- sqrt(SpatialDeterminantC(g11r,g12r,g13r,g22r,g23r,g33r));
-
- // minus call to p2c
- prim2conC(&w_lorentzminus[idx], &densminus[idx], &sxminus[idx],
- &syminus[idx], &szminus[idx], &tauminus[idx],
- g11l,g12l,g13l,g22l,g23l,g33l,
- savg_detl,rhominus[idx], velxminus[idx], velyminus[idx],
- velzminus[idx], epsminus[idx], pressminus[idx]);
-
-
- // plus call to p2c
- prim2conC(&w_lorentzplus[idx], &densplus[idx], &sxplus[idx],
- &syplus[idx], &szplus[idx], &tauplus[idx],
- g11r,g12r,g13r,g22r,g23r,g33r,
- savg_detr,rhoplus[idx], velxplus[idx], velyplus[idx],
- velzplus[idx], epsplus[idx], pressplus[idx]);
- }
-
-} // end function Conservative2PrimitiveC
-
static inline double SpatialDeterminantC(double gxx, double gxy,
double gxz, double gyy,
double gyz, double gzz)
File [modified]: GRHydro_Prim2ConM.F90
Delta lines: +11 -1
===================================================================
--- trunk/src/GRHydro_Prim2ConM.F90 2014-04-15 19:49:47 UTC (rev 632)
+++ trunk/src/GRHydro_Prim2ConM.F90 2014-04-15 19:49:50 UTC (rev 633)
@@ -85,6 +85,17 @@
pBprim = loc(Bvec)
end if
+ ! currently only the C++ code uses the speed-of-sound grid function, however
+ ! in the MHD code, both C++ and F90 code share the same prim2con routine (this
+ ! one). GRHydro_SpeedOfSound is a misnomer since it computes both cs2 and
+ ! pressure so we call it first so that the actual prim2con below overwrites
+ ! the pressur again. This is certainly not the "right" way of doing this but
+ ! it lets us have a working version for now.
+ ! TODO: extend prim2con in GRHydro_Prim2Con.c to handle MHD.
+ if (use_cxx_code .ne. 0) then
+ call GRHydro_SpeedOfSound(cctkGH)
+ end if
+
! constraint transport needs to be able to average fluxes in the directions
! other that flux_direction, which in turn need the primitives on interfaces
@@ -211,7 +222,6 @@
!$OMP END PARALLEL DO
endif
-
end subroutine primitive2conservativeM
/*@@
More information about the Commits
mailing list