Commits April 2014

commits@lists.einsteintoolkit.org

8 participants
175 discussions

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 609)
by rhaas＠tapir.caltech.edu 15 Apr '14

15 Apr '14

User: rhaas Date: 2014/04/15 02:48 PM Added: /trunk/src/ GRHydro_ePPM.cc Modified: /trunk/src/ GRHydro_PPMReconstruct_drv_opt.cc, GRHydro_PPM_opt.h, make.code.defn Log: GRHydro: New C++ ePPM algorithm Tested with TOV star. Improved comments and pointed out a difference to F90 code which makes the C++ code more accurate. File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_PPMReconstruct_drv_opt.cc Delta lines: +71 -36 =================================================================== --- trunk/src/GRHydro_PPMReconstruct_drv_opt.cc 2014-04-15 19:48:40 UTC (rev 608) +++ trunk/src/GRHydro_PPMReconstruct_drv_opt.cc 2014-04-15 19:48:44 UTC (rev 609) @@ -89,15 +89,13 @@ // 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 + // if use_enhanced_ppm, allow old PPM on one level + bool apply_enhanced_ppm = false; + if (GRHydro_oppm_reflevel == (-1) || + *GRHydro_reflevel != GRHydro_oppm_reflevel) + apply_enhanced_ppm = use_enhanced_ppm != 0; + else + apply_enhanced_ppm = false; int i,j,k; @@ -178,54 +176,91 @@ double* restrict Bvcyplus, \ double* restrict Bvczplus, \ double* restrict psidcplus); - static ppm1d_cxx_ptr_t ppm1d_cxx_funcs[32]; + static ppm1d_cxx_ptr_t ppm1d_cxx_funcs[64]; ppm1d_cxx_ptr_t ppm1d_cxx_func; - #define MKINDEX(do_temp,do_ye,do_mhd,dc_flag,do_ppm_detect) ((do_temp<<4)|(do_ye<<3)|(do_mhd<<2)|(dc_flag<<1)|(do_ppm_detect)) - #define SETPOINTER(do_temp,do_ye,do_mhd,dc_flag,do_ppm_detect) \ - ppm1d_cxx_funcs[MKINDEX(do_temp,do_ye,do_mhd,dc_flag,do_ppm_detect)] = \ - GRHydro_ppm1d_cxx<do_temp,do_ye,do_mhd,dc_flag,do_ppm_detect> + #define MKINDEX(do_temp,do_ye,do_mhd,dc_flag,do_ppm_detect,do_eppm) ((do_eppm<<5)|(do_temp<<4)|(do_ye<<3)|(do_mhd<<2)|(dc_flag<<1)|(do_ppm_detect)) + #define SETPOINTER(do_temp,do_ye,do_mhd,dc_flag,do_ppm_detect,do_eppm) \ + if (!do_eppm) \ + ppm1d_cxx_funcs[MKINDEX(do_temp,do_ye,do_mhd,dc_flag,do_ppm_detect, do_eppm)] = \ + GRHydro_ppm1d_cxx<do_temp,do_ye,do_mhd,dc_flag,do_ppm_detect>; \ + else \ + ppm1d_cxx_funcs[MKINDEX(do_temp,do_ye,do_mhd,dc_flag,do_ppm_detect, do_eppm)] = \ + GRHydro_eppm1d_cxx<do_temp,do_ye,do_mhd,dc_flag,do_ppm_detect>; + + // --- original PPM + // first stuff without MHD + SETPOINTER(false,false,false,false,false,false); + SETPOINTER(false,false,false,false,true,false); + // do_temp can only be true of do_ye is also true + SETPOINTER(true,true,false,false,false,false); + SETPOINTER(true,true,false,false,true,false); + // but do_ye can be true and do_temp can be false + SETPOINTER(false,true,false,false,false,false); + SETPOINTER(false,true,false,false,true,false); + // with MHD, but dc_flag false + SETPOINTER(false,false,true,false,false,false); + SETPOINTER(false,false,true,false,true,false); + // do_temp can only be true of do_ye is also true + SETPOINTER(true,true,true,false,false,false); + SETPOINTER(true,true,true,false,true,false); + // but do_ye can be true and do_temp can be false + SETPOINTER(false,true,true,false,false,false); + SETPOINTER(false,true,true,false,true,false); + + // with MHD, dc_flag true + SETPOINTER(false,false,true,true,false,false); + SETPOINTER(false,false,true,true,true,false); + // do_temp can only be true of do_ye is also true + SETPOINTER(true,true,true,true,false,false); + SETPOINTER(true,true,true,true,true,false); + // but do_ye can be true and do_temp can be false + SETPOINTER(false,true,true,true,false,false); + SETPOINTER(false,true,true,true,true,false); + + // ---- enhanced PPM // first stuff without MHD - SETPOINTER(false,false,false,false,false); - SETPOINTER(false,false,false,false,true); + SETPOINTER(false,false,false,false,false,true); + SETPOINTER(false,false,false,false,true,true); // do_temp can only be true of do_ye is also true - SETPOINTER(true,true,false,false,false); - SETPOINTER(true,true,false,false,true); + SETPOINTER(true,true,false,false,false,true); + SETPOINTER(true,true,false,false,true,true); // but do_ye can be true and do_temp can be false - SETPOINTER(false,true,false,false,false); - SETPOINTER(false,true,false,false,true); + SETPOINTER(false,true,false,false,false,true); + SETPOINTER(false,true,false,false,true,true); // with MHD, but dc_flag false - SETPOINTER(false,false,true,false,false); - SETPOINTER(false,false,true,false,true); + SETPOINTER(false,false,true,false,false,true); + SETPOINTER(false,false,true,false,true,true); // do_temp can only be true of do_ye is also true - SETPOINTER(true,true,true,false,false); - SETPOINTER(true,true,true,false,true); + SETPOINTER(true,true,true,false,false,true); + SETPOINTER(true,true,true,false,true,true); // but do_ye can be true and do_temp can be false - SETPOINTER(false,true,true,false,false); - SETPOINTER(false,true,true,false,true); + SETPOINTER(false,true,true,false,false,true); + SETPOINTER(false,true,true,false,true,true); // with MHD, dc_flag true - SETPOINTER(false,false,true,true,false); - SETPOINTER(false,false,true,true,true); + SETPOINTER(false,false,true,true,false,true); + SETPOINTER(false,false,true,true,true,true); // do_temp can only be true of do_ye is also true - SETPOINTER(true,true,true,true,false); - SETPOINTER(true,true,true,true,true); + SETPOINTER(true,true,true,true,false,true); + SETPOINTER(true,true,true,true,true,true); // but do_ye can be true and do_temp can be false - SETPOINTER(false,true,true,true,false); - SETPOINTER(false,true,true,true,true); - + SETPOINTER(false,true,true,true,false,true); + SETPOINTER(false,true,true,true,true,true); + + const bool do_temp = *evolve_temper; const bool do_ye = *evolve_Y_e; const bool do_mhd = *evolve_MHD; const bool dc_flag = clean_divergence; const bool do_ppm_detect = ppm_detect; - ppm1d_cxx_func = ppm1d_cxx_funcs[MKINDEX(do_temp,do_ye,do_mhd,dc_flag,do_ppm_detect)]; + ppm1d_cxx_func = ppm1d_cxx_funcs[MKINDEX(do_temp,do_ye,do_mhd,dc_flag,do_ppm_detect,apply_enhanced_ppm)]; if(ppm1d_cxx_func == NULL) { CCTK_VError(__LINE__, __FILE__, CCTK_THORNSTRING, - "Combination do_temp = %d do_ye = %d do_mhd = %d dc_flag = %d do_ppm_detect = %d is not supported for PPM. Please instatate the respective template", - do_temp, do_ye, do_mhd, dc_flag, do_ppm_detect); + "Combination do_temp = %d do_ye = %d do_mhd = %d dc_flag = %d do_ppm_detect = %d do_eppm = %d is not supported for PPM. Please instantiate the respective template", + do_temp, do_ye, do_mhd, dc_flag, do_ppm_detect, apply_enhanced_ppm); } //index = i + j*nx + k*nx*ny; File [modified]: GRHydro_PPM_opt.h Delta lines: +41 -0 =================================================================== --- trunk/src/GRHydro_PPM_opt.h 2014-04-15 19:48:40 UTC (rev 608) +++ trunk/src/GRHydro_PPM_opt.h 2014-04-15 19:48:44 UTC (rev 609) @@ -36,3 +36,44 @@ double* restrict Bvcyplus, double* restrict Bvczplus, double* restrict psidcplus); + +template<bool do_temp, bool do_ye, bool do_mhd, + bool dc_flag, bool do_ppm_detect> +void GRHydro_eppm1d_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); + + File [added]: GRHydro_ePPM.cc Delta lines: +367 -0 =================================================================== --- trunk/src/GRHydro_ePPM.cc (rev 0) +++ trunk/src/GRHydro_ePPM.cc 2014-04-15 19:48:44 UTC (rev 609) @@ -0,0 +1,367 @@ +#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)) + +#define MIN3(a,b,c) (MIN(a, MIN(b, c))) +#define MIN4(a,b,c,d) (MIN(a, MIN(b, MIN(c, d)))) + +#define MAX5(a,b,c,d,e) (MAX(a, MAX(b, MAX(c, MAX(d, e))))) + +using namespace std; + + +/** + ePPM Reconstruction (see Reisswig et al. 2013) based on McCorquodale & Colella (2011). + Note: This routine differs from the Fortran routine (and the algorithm described in Reisswig et al. 2013) + in the following way: We do not treat the specific internal energy "epsilon" in any special way. + Rather, to avoid reconstruction to eps<0 for gamma-laws, we rely on + 'eps<0'-error checks after reconstruction (a similar treatment is also necessary for other reconstruction methods such as WENO5 and MP5). +*/ + + +/* + 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 double approx_at_cell_interface(const double* const restrict 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 double limit(const double* const restrict a, const double ah, const double C, const int i) +{ + if ((MIN(a[i],a[i+1]) <= ah) && (ah <= MAX(a[i],a[i+1]))) { + return ah; + } else { + const double D2a = 3.0 * ((a[i] + a[i+1]) - 2.0*ah ); + const double D2aL = ((a[i-1] + a[i+1]) - 2.0*a[i] ); + const double D2aR = ((a[i] + a[i+2]) - 2.0*a[i+1]); + const double D2aLim = copysign(1.0, D2a)*MIN3(C*fabs(D2aL), C*fabs(D2aR), fabs(D2a)) * 1.0/3.0; + if (D2a*D2aR >= 0 && D2a*D2aL >= 0) + return 0.5*(a[i]+a[i+1]) - D2aLim; + else + return 0.5*(a[i]+a[i+1]); + } + return 0; +} + + +/// Monotonicity of PPM 2011 (McCorquodale & Colella 2011), Sec. 2.4.1, Eq. 23-34. +/// This does not use the check for deviations from a cubic. Thus, it gets away with only 3 stencil points! +static inline void monotonize(double* const restrict aminus, + const double* const restrict a, + double* const restrict aplus, + const double C, + const int i) +{ + double D2aLim = 0; + double rhi = 0; + const double daplus = aplus[i]-a[i]; + const double daminus = a[i]-aminus[i]; + if (daplus*daminus <= 0 || (a[i-2]-a[i])*(a[i]-a[i+2]) <= 0) { + const double D2a = - (12.0*a[i] - 6.0*(aminus[i]+aplus[i])); + const double D2aC = (a[i-1] + a[i+1]) - 2.0*a[i]; + const double D2aL = (a[i-2] + a[i] ) - 2.0*a[i-1]; + const double D2aR = (a[i] + a[i+2]) - 2.0*a[i+1]; + if (copysign(1.0, D2a) == copysign(1.0, D2aC) && copysign(1.0, D2a) == copysign(1.0, D2aL) && copysign(1.0, D2a) == copysign(1.0, D2aR)) + D2aLim = copysign(1.0, D2a) * MIN4(C*fabs(D2aL), C*fabs(D2aR), C*fabs(D2aC), fabs(D2a)); + if (!(fabs(D2a) <= 1e-12*MAX5(fabs(a[i-2]), fabs(a[i-1]), fabs(a[i]), fabs(a[i+1]), fabs(a[i+2])))) + rhi = D2aLim / D2a; + if (! (rhi >= 1.0 - 1e-12)) { + if (daplus*daminus < 0) { + aplus[i] = a[i] + daplus * rhi; + aminus[i] = a[i] - daminus * rhi; + } else if (fabs(daminus) >= 2.0*fabs(daplus)) { + aminus[i] = a[i] - (2.0*(1.0-rhi)*daplus + rhi*daminus); + } else if (fabs(daplus) >= 2.0*fabs(daminus)) { + aplus[i] = a[i] + (2.0*(1.0-rhi)*daminus + rhi*daplus); + } + } + //trivial_rp(i-1) = .false. + //trivial_rp(i) = .false. + } else { + if (fabs(daplus) >= 2.0*fabs(daminus)) + aplus[i] = a[i] + 2.0*daminus; + else if (fabs(daminus) >= 2.0*fabs(daplus)) + aminus[i] = a[i] - 2.0*daplus; + //trivial_rp(i-1) = .false. + //trivial_rp(i) = .false. + } + + return; +} + + +template<bool do_temp, bool do_ye, bool do_mhd, + bool dc_flag, bool do_ppm_detect> +void GRHydro_eppm1d_cxx(const int nx, + const double dx, + const double* const restrict rho, + const double* const restrict velx, + const double* const restrict vely, + const double* const restrict velz, + const double* const restrict eps, + const double* const restrict press, + const double* const restrict temp, + const double* const restrict ye, + const double* const restrict Bvcx, + const double* const restrict Bvcy, + const double* const restrict Bvcz, + const double* const restrict psidc, + double* const restrict rhominus, + double* const restrict velxminus, + double* const restrict velyminus, + double* const restrict velzminus, + double* const restrict epsminus, + double* const restrict tempminus, + double* const restrict yeminus, + double* const restrict Bvcxminus, + double* const restrict Bvcyminus, + double* const restrict Bvczminus, + double* const restrict psidcminus, + double* const restrict rhoplus, + double* const restrict velxplus, + double* const restrict velyplus, + double* const restrict velzplus, + double* const restrict epsplus, + double* const restrict tempplus, + double* const restrict yeplus, + double* const restrict Bvcxplus, + double* const restrict Bvcyplus, + double* const restrict Bvczplus, + double* const restrict psidcplus) +{ + DECLARE_CCTK_PARAMETERS; + + + // We initialize "plus" \equiv a_j+1/2 with (16) via APPROX_AT_CELL_INTERFACE, + // then checking for (13) of Colella & Sekora 2008 and applying + // (18) and (19) if (13) is not satisfied. This is done with LIMIT. + for (int i=1; i < nx-2; ++i) { + rhoplus[i] = approx_at_cell_interface(rho, i); + rhoplus[i] = limit(rho, rhoplus[i], enhanced_ppm_C2, i); + rhominus[i+1] = rhoplus[i]; + + epsplus[i] = approx_at_cell_interface(eps, i); + epsplus[i] = limit(eps, epsplus[i], enhanced_ppm_C2, i); + epsminus[i+1] = epsplus[i]; + + velxplus[i] = approx_at_cell_interface(velx, i); + velxplus[i] = limit(velx, velxplus[i], enhanced_ppm_C2, i); + velxminus[i+1] = velxplus[i]; + + velyplus[i] = approx_at_cell_interface(vely, i); + velyplus[i] = limit(vely, velyplus[i], enhanced_ppm_C2, i); + velyminus[i+1] = velyplus[i]; + + velzplus[i] = approx_at_cell_interface(velz, i); + velzplus[i] = limit(velz, velzplus[i], enhanced_ppm_C2, i); + velzminus[i+1] = velzplus[i]; + + if(do_ye) { + if(do_temp) { + tempplus[i] = approx_at_cell_interface(temp, i); + tempplus[i] = limit(temp, tempplus[i], enhanced_ppm_C2, i); + tempminus[i+1] = tempplus[i]; + } + yeplus[i] = approx_at_cell_interface(ye, i); + yeplus[i] = limit(ye, yeplus[i], enhanced_ppm_C2, i); + yeminus[i+1] = yeplus[i]; + } + if(do_mhd) { + Bvcxplus[i] = approx_at_cell_interface(Bvcx, i); + Bvcxplus[i] = limit(Bvcx, Bvcxplus[i], enhanced_ppm_C2, i); + Bvcxminus[i+1] = Bvcxplus[i]; + + Bvcyplus[i] = approx_at_cell_interface(Bvcy, i); + Bvcyplus[i] = limit(Bvcy, Bvcyplus[i], enhanced_ppm_C2, i); + Bvcyminus[i+1] = Bvcyplus[i]; + + Bvczplus[i] = approx_at_cell_interface(Bvcz, i); + Bvczplus[i] = limit(Bvcz, Bvczplus[i], enhanced_ppm_C2, i); + Bvczminus[i+1] = Bvczplus[i]; + + if(dc_flag) { + psidcplus[i] = approx_at_cell_interface(psidc, i); + psidcplus[i] = limit(psidc, psidcplus[i], enhanced_ppm_C2, i); + psidcminus[i+1] = psidcplus[i]; + } + } + + } + + // Monotonize + for(int i=GRHydro_stencil-1; i<nx-GRHydro_stencil+1; ++i) { + monotonize(rhominus,rho,rhoplus,enhanced_ppm_C2,i); + monotonize(epsminus,eps,epsplus,enhanced_ppm_C2,i); + monotonize(velxminus,velx,velxplus,enhanced_ppm_C2,i); + monotonize(velyminus,vely,velyplus,enhanced_ppm_C2,i); + monotonize(velzminus,velz,velzplus,enhanced_ppm_C2,i); + if(do_ye) { + if(do_temp) { + monotonize(tempminus,temp,tempplus,enhanced_ppm_C2,i); + } + monotonize(yeminus,ye,yeplus,enhanced_ppm_C2,i); + } + if(do_mhd) { + monotonize(Bvcxminus,Bvcx,Bvcxplus,enhanced_ppm_C2,i); + monotonize(Bvcyminus,Bvcy,Bvcyplus,enhanced_ppm_C2,i); + monotonize(Bvczminus,Bvcz,Bvczplus,enhanced_ppm_C2,i); + if(dc_flag) { + monotonize(psidcminus,psidc,psidcplus,enhanced_ppm_C2,i); + } + } + } + + + // Finally, apply flattening + vector<double> dpress(nx,0); + for(int i=1; i<nx-1; ++i) + dpress[i] = (press[i+1]-press[i-1]); + + vector<double> tilde_flatten(nx,0); + for(int i=2; i<nx-2; ++i) { + const double dpress2 = press[i+2] - press[i-2]; + const 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) { + const 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 + + 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 void GRHydro_eppm1d_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 2014-04-15 19:48:40 UTC (rev 608) +++ trunk/src/make.code.defn 2014-04-15 19:48:44 UTC (rev 609) @@ -84,6 +84,7 @@ Con2PrimM_fortran_interfaces.F90 \ GRHydro_Prim2Con.c \ GRHydro_PPM.cc \ + GRHydro_ePPM.cc \ GRHydro_PPMReconstruct_drv_opt.cc \ GRHydro_Wrappers.F90

1 0

[svn:einsteintoolkit] GRHydro/trunk/test/ (Rev. 608)
by rhaas＠tapir.caltech.edu 15 Apr '14

15 Apr '14

User: rhaas Date: 2014/04/15 02:48 PM Added: /trunk/test/ GRHydro_test_shock_hllc.par /trunk/test/GRHydro_test_shock_hllc/ dens_3D_diagonal.xg, dens_x_[7][7].xg, dens_y_[7][7].xg, dens_z_[7][7].xg, eps_3D_diagonal.xg, eps_x_[7][7].xg, eps_y_[7][7].xg, eps_z_[7][7].xg, press_3D_diagonal.xg, press_x_[7][7].xg, press_y_[7][7].xg, press_z_[7][7].xg, rho_3D_diagonal.xg, rho_x_[7][7].xg, rho_y_[7][7].xg, rho_z_[7][7].xg, scon[0]_3D_diagonal.xg, scon[0]_x_[7][7].xg, scon[0]_y_[7][7].xg, scon[0]_z_[7][7].xg, scon[1]_3D_diagonal.xg, scon[1]_x_[7][7].xg, scon[1]_y_[7][7].xg, scon[1]_z_[7][7].xg, scon[2]_3D_diagonal.xg, scon[2]_x_[7][7].xg, scon[2]_y_[7][7].xg, scon[2]_z_[7][7].xg, tau_3D_diagonal.xg, tau_x_[7][7].xg, tau_y_[7][7].xg, tau_z_[7][7].xg, vel[0]_3D_diagonal.xg, vel[0]_x_[7][7].xg, vel[0]_y_[7][7].xg, vel[0]_z_[7][7].xg, vel[1]_3D_diagonal.xg, vel[1]_x_[7][7].xg, vel[1]_y_[7][7].xg, vel[1]_z_[7][7].xg, vel[2]_3D_diagonal.xg, vel[2]_x_[7][7].xg, vel[2]_y_[7][7].xg, vel[2]_z_[7][7].xg, w_lorentz_3D_diagonal.xg, w_lorentz_x_[7][7].xg, w_lorentz_y_[7][7].xg, w_lorentz_z_[7][7].xg Log: GRHydro: add hllc shock tube test File Changes: Directory: /trunk/test/GRHydro_test_shock_hllc/ =============================================== File [added]: dens_3D_diagonal.xg Delta lines: <Unavailable> File [added]: dens_x_[7][7].xg Delta lines: <Unavailable> File [added]: dens_y_[7][7].xg Delta lines: <Unavailable> File [added]: dens_z_[7][7].xg Delta lines: <Unavailable> File [added]: eps_3D_diagonal.xg Delta lines: <Unavailable> File [added]: eps_x_[7][7].xg Delta lines: <Unavailable> File [added]: eps_y_[7][7].xg Delta lines: <Unavailable> File [added]: eps_z_[7][7].xg Delta lines: <Unavailable> File [added]: press_3D_diagonal.xg Delta lines: <Unavailable> File [added]: press_x_[7][7].xg Delta lines: <Unavailable> File [added]: press_y_[7][7].xg Delta lines: <Unavailable> File [added]: press_z_[7][7].xg Delta lines: <Unavailable> File [added]: rho_3D_diagonal.xg Delta lines: <Unavailable> File [added]: rho_x_[7][7].xg Delta lines: <Unavailable> File [added]: rho_y_[7][7].xg Delta lines: <Unavailable> File [added]: rho_z_[7][7].xg Delta lines: <Unavailable> File [added]: scon[0]_3D_diagonal.xg Delta lines: <Unavailable> File [added]: scon[0]_x_[7][7].xg Delta lines: <Unavailable> File [added]: scon[0]_y_[7][7].xg Delta lines: <Unavailable> File [added]: scon[0]_z_[7][7].xg Delta lines: <Unavailable> File [added]: scon[1]_3D_diagonal.xg Delta lines: <Unavailable> File [added]: scon[1]_x_[7][7].xg Delta lines: <Unavailable> File [added]: scon[1]_y_[7][7].xg Delta lines: <Unavailable> File [added]: scon[1]_z_[7][7].xg Delta lines: <Unavailable> File [added]: scon[2]_3D_diagonal.xg Delta lines: <Unavailable> File [added]: scon[2]_x_[7][7].xg Delta lines: <Unavailable> File [added]: scon[2]_y_[7][7].xg Delta lines: <Unavailable> File [added]: scon[2]_z_[7][7].xg Delta lines: <Unavailable> File [added]: tau_3D_diagonal.xg Delta lines: <Unavailable> File [added]: tau_x_[7][7].xg Delta lines: <Unavailable> File [added]: tau_y_[7][7].xg Delta lines: <Unavailable> File [added]: tau_z_[7][7].xg Delta lines: <Unavailable> File [added]: vel[0]_3D_diagonal.xg Delta lines: <Unavailable> File [added]: vel[0]_x_[7][7].xg Delta lines: <Unavailable> File [added]: vel[0]_y_[7][7].xg Delta lines: <Unavailable> File [added]: vel[0]_z_[7][7].xg Delta lines: <Unavailable> File [added]: vel[1]_3D_diagonal.xg Delta lines: <Unavailable> File [added]: vel[1]_x_[7][7].xg Delta lines: <Unavailable> File [added]: vel[1]_y_[7][7].xg Delta lines: <Unavailable> File [added]: vel[1]_z_[7][7].xg Delta lines: <Unavailable> File [added]: vel[2]_3D_diagonal.xg Delta lines: <Unavailable> File [added]: vel[2]_x_[7][7].xg Delta lines: <Unavailable> File [added]: vel[2]_y_[7][7].xg Delta lines: <Unavailable> File [added]: vel[2]_z_[7][7].xg Delta lines: <Unavailable> File [added]: w_lorentz_3D_diagonal.xg Delta lines: <Unavailable> File [added]: w_lorentz_x_[7][7].xg Delta lines: <Unavailable> File [added]: w_lorentz_y_[7][7].xg Delta lines: <Unavailable> File [added]: w_lorentz_z_[7][7].xg Delta lines: <Unavailable> Directory: /trunk/test/ ======================= File [added]: GRHydro_test_shock_hllc.par Delta lines: <Unavailable>

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 607)
by rhaas＠tapir.caltech.edu 15 Apr '14

15 Apr '14

User: rhaas Date: 2014/04/15 02:48 PM Added: /trunk/src/ GRHydro_HLLC.F90 Modified: /trunk/ param.ccl /trunk/src/ GRHydro_RiemannSolve.F90, make.code.defn Log: GRHydro: add hllc solver * fix sign of eigenvalue of E coefficient in HLLC, probably still incorrect for beta!=0 File Changes: Directory: /trunk/src/ ====================== File [added]: GRHydro_HLLC.F90 Delta lines: +666 -0 =================================================================== --- trunk/src/GRHydro_HLLC.F90 (rev 0) +++ trunk/src/GRHydro_HLLC.F90 2014-04-15 19:48:34 UTC (rev 607) @@ -0,0 +1,666 @@ + /*@@ + @file GRHydro_HLLC.F90 + @date Sat Jan 26 01:40:14 2002 + @author Ian Hawke, Pedro Montero, Toni Font, Roland Haas + @desc + The HLLC solver. Called from the wrapper function, so works in + all directions. + @enddesc + @@*/ + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" + +#include "GRHydro_Macros.h" +#include "SpaceMask.h" + + /*@@ + @routine GRHydro_HLLC + @date Sat Jan 26 01:41:02 2002 + @author Ian Hawke, Pedro Montero, Toni Font, Roland Haas + @desc + The HLLC solver. Sufficiently simple that its just one big routine. + @enddesc + @calls + @calledby + @history + Based entirely on Matt Duez notes on HLL and HLLC + @endhistory + +@@*/ + +module GRHydro_hllc_helpers + + implicit none + + contains + ! using nested subroutines seems to not work (either because gfortran has a + ! bug or because they are not allowed with OpenMP). So we have to pass each + ! argument explictly. ok for now since I plan to re-implement this in C + ! anyway. + subroutine compute_stardata(charstar,pressstar,consm_i,consp,sqrt_avg_det,pressminus,pressplus,& + uxxh,uxyh,uxzh,uyyh,uyzh,uzzh,usendh,avg_beta,avg_alp,vxR,vxL,charmin,charmax,i,j,k,& + cctk_ash1,cctk_ash2,cctk_ash3,xoffset,yoffset,zoffset,flux_direction) + + implicit none + + CCTK_REAL, intent(out) :: charstar, pressstar + + integer, intent(in) :: cctk_ash1, cctk_ash2, cctk_ash3 + CCTK_REAL, intent(in), DIMENSION(5) :: consm_i, consp + CCTK_REAL, intent(in) :: sqrt_avg_det, avg_beta, avg_alp + CCTK_REAL, intent(in), DIMENSION(cctk_ash1,cctk_ash2,cctk_ash3) :: pressminus, pressplus + CCTK_REAL, intent(in) :: uxxh,uxyh,uxzh,uyyh,uyzh,uzzh,usendh + CCTK_REAL, intent(in) :: charmin, charmax + CCTK_REAL, intent(in) :: vxR, vxL + integer, intent(in) :: i,j,k + CCTK_INT, intent(in) :: xoffset, yoffset, zoffset, flux_direction + + ! naming convention: in this routine we use the convention of Duez (and by + ! extenstion Mignone et al.MNRAS 364 126 (2005) which means that L<->plus, + ! R<->minus + ! for convenience, use s^k = D h u^k = D h g^{k \mu} u_\mu + ! ie. this is NOT s^k = gamma^{kj} s_j + ! doing the math, I find: + ! s^k = gamma^{kj} s_j - beta^k/alpha sqrt(detgamma) rho h W**2 + ! where the last part is rho h W**2 = tau+D+press + CCTK_REAL :: sXL, sXR + ! the coefficients of the quadratic equaion determining lambdastar + CCTK_REAL :: AL,AR,BL,BR,CL,CR,EL,ER + + ! the combination D*h*W which appear more than once + CCTK_REAL :: DhWL,DhWR + CCTK_REAL :: lamstar_tilde ! eq (22) + + + ! compute sX + DhWR = consm_i(1)+consm_i(5)+ & + sqrt_avg_det*pressminus(i+xoffset,j+yoffset,k+zoffset) + DhWL = consp (1)+consp (5)+ & + sqrt_avg_det*pressplus (i ,j ,k ) + if (flux_direction == 1) then + sXR = uxxh*consm_i(2)+uxyh*consm_i(3)+uxzh*consm_i(4) - & + avg_beta/avg_alp*DhWR + sXL = uxxh*consp (2)+uxyh*consp (3)+uxzh*consp (4) - & + avg_beta/avg_alp*DhWL + else if (flux_direction == 2) then + sXR = uxyh*consm_i(2)+uyyh*consm_i(3)+uyzh*consm_i(4) - & + avg_beta/avg_alp*DhWR + sXL = uxyh*consp (2)+uyyh*consp (3)+uyzh*consp (4) - & + avg_beta/avg_alp*DhWL + else if (flux_direction == 3) then + sXR = uxzh*consm_i(2)+uyzh*consm_i(3)+uzzh*consm_i(4) - & + avg_beta/avg_alp*DhWR + sXL = uxzh*consp (2)+uyzh*consp (3)+uzzh*consp (4) - & + avg_beta/avg_alp*DhWL + else + call CCTK_ERROR("Flux direction not x,y,z") + end if + + ! compute the coefficients A-D of Duez (23)-(26) + AR = (consm_i(5)+consm_i(1))*(charmax+avg_beta) - (avg_alp*sXR+DhWR*avg_beta) + AL = (consp (5)+consp (1))*(charmin+avg_beta) - (avg_alp*sXL+DhWL*avg_beta) + CR = avg_alp*usendh + CL = CR + BR = CR*sqrt_avg_det*pressminus(i+xoffset,j+yoffset,k+zoffset) - & + (sXR+DhWR*avg_beta/avg_alp)*(charmax-vXR) + BL = CL*sqrt_avg_det*pressplus (i ,j ,k ) - & + (sXL+DhWL*avg_beta/avg_alp)*(charmin-vXL) + ER = -(avg_beta + charmax) ! sing ERROR in Duez (26) + EL = -(avg_beta + charmin) + + ! we require the negative root only (see Mignone et al), this is (27) of + ! Duez + call solve_quadratic(EL*AR-AL*ER, -(AL*CR+BL*ER-CL*AR-EL*BR), BR*CL-BL*CR, & + lamstar_tilde) + + charstar = avg_alp * lamstar_tilde - avg_beta ! eq (22) + ! pressstar is equal left and right of contact discontinuity + pressstar = (AL*lamstar_tilde + BL)/((CL+EL*lamstar_tilde)*sqrt_avg_det) + + end subroutine + + ! this function is a straightforward transcription of the GSL solve_quadratic + ! functions. Both codes are GPL. + ! Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough + subroutine solve_quadratic(a,b,c, res) + + implicit none + + CCTK_REAL, intent(in) :: a,b,c + CCTK_REAL, intent(out) :: res + + CCTK_REAL, parameter :: one = 1 + CCTK_REAL :: disc, r, sgnb, temp, r1, r2 + character(len=256) :: warnline + + disc = b * b - 4 * a * c + + if (a == 0) then! Handle linear case + if (b == 0) then + write (warnline, '("Could not solve quadratic equation with coefficients ",3g22.15)') a,b,c + call CCTK_ERROR(warnline) + else + res = -c / b + end if + return + end if + + if (disc > 0) then + if (b == 0) then + r = abs (0.5d0 * sqrt (disc) / a) + res = -r + !*x1 = r ! only want negative root + else + sgnb = sign(one, b) + temp = -0.5d0 * (b + sgnb * sqrt (disc)) + r1 = temp / a + r2 = c / temp + + res = r2 ! only want negative (second) root + + !if (r1 < r2) then + ! res = r1 + ! !*x1 = r2 ! only want negative root + !else + ! res = r2 + ! !*x1 = r1 ! only want negative root + !end if + end if + else if (disc == 0) then + res = -0.5d0 * b / a + !*x1 = -0.5 * b / a ! only want one root + else + write (warnline, '("Could not solve quadratic equation with coefficients ",3g22.15)') a,b,c + call CCTK_ERROR(warnline) + end if + + end subroutine +end module GRHydro_hllc_helpers + + +subroutine GRHydro_HLLC(CCTK_ARGUMENTS) + USE GRHydro_Eigenproblem + USE GRHydro_hllc_helpers + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + integer :: i, j, k + integer :: keytemp + CCTK_REAL, dimension(5) :: consp,consm_i,fplus,fminus,lamplus + CCTK_REAL, dimension(5) :: f1,lamminus, consstar + CCTK_REAL :: charmin, charmax, charstar, avg_alp,avg_det + CCTK_REAL :: gxxh, gxyh, gxzh, gyyh, gyzh, gzzh, uxxh, uxyh, & + uxzh, uyyh, uyzh, uzzh, avg_beta, usendh + CCTK_REAL :: pressstar ! NOTE: this name is unfortuante for MHD where it + ! conflicts + CCTK_REAL :: sqrt_avg_det + ! Duez und GRHydro define different quantities v^i, we have + ! v^i_Duez = alp*v^i_GRHydro - beta^i + CCTK_REAL :: vXL, vXR + CCTK_REAL :: velstar ! GRhydro's velocity in the star region + + CCTK_INT :: type_bits, trivial + + ! sign requires its arguments to be of identical KIND + CCTK_REAL, parameter :: one = 1d0 + + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + CCTK_REAL, DIMENSION(cctk_ash1,cctk_ash2,cctk_ash3) :: g11, g12, g13, g22, g23, g33 + pointer (pg11,g11), (pg12,g12), (pg13,g13), (pg22,g22), (pg23,g23), (pg33,g33) + CCTK_REAL, DIMENSION(cctk_ash1,cctk_ash2,cctk_ash3) :: beta1, beta2, beta3 + pointer (pbeta1,beta1), (pbeta2,beta2), (pbeta3,beta3) + CCTK_REAL, DIMENSION(cctk_ash1,cctk_ash2,cctk_ash3,3) :: vup + pointer (pvup,vup) + + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + pg11 = loc(gaa) + pg12 = loc(gab) + pg13 = loc(gac) + pg22 = loc(gbb) + pg23 = loc(gbc) + pg33 = loc(gcc) + pbeta1 = loc(betaa) + pbeta2 = loc(betab) + pbeta3 = loc(betac) + pvup = loc(lvel) + else + pg11 = loc(gxx) + pg12 = loc(gxy) + pg13 = loc(gxz) + pg22 = loc(gyy) + pg23 = loc(gyz) + pg33 = loc(gzz) + pbeta1 = loc(betax) + pbeta2 = loc(betay) + pbeta3 = loc(betaz) + pvup = loc(vel) + end if + + if (flux_direction == 1) then + call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemX") + call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemX", & + &"trivial") + else if (flux_direction == 2) then + call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemY") + call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemY", & + &"trivial") + else if (flux_direction == 3) then + call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemZ") + call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemZ", & + &"trivial") + else + call CCTK_ERROR("Flux direction not x,y,z") + end if + + if(evolve_temper.eq.1.and.reconstruct_temper.eq.1) then + keytemp = 1 + else + keytemp = 0 + endif + + !$OMP PARALLEL DO PRIVATE(k,j,i,f1,lamminus,lamplus,consp,consm_i,consstar, & + !$OMP fplus,fminus,avg_beta,avg_alp, & + !$OMP avg_det,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, & + !$OMP uxxh,uxyh,uxzh,uyyh,uyzh,uzzh, & + !$OMP usendh, charmin, charmax, charstar, & + !$OMP pressstar, sqrt_avg_det, vXL, vXR, velstar) + do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil + do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil + do i = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + + f1 = 0.d0 + lamminus = 0.d0 + lamplus = 0.d0 + consp = 0.d0 + consm_i = 0.d0 + fplus = 0.d0 + fminus = 0.d0 + +!!$ Set the left (p for plus) and right (m_i for minus, i+1) states + + consp(1) = densplus(i,j,k) + consp(2) = sxplus(i,j,k) + consp(3) = syplus(i,j,k) + consp(4) = szplus(i,j,k) + consp(5) = tauplus(i,j,k) + + consm_i(1) = densminus(i+xoffset,j+yoffset,k+zoffset) + consm_i(2) = sxminus(i+xoffset,j+yoffset,k+zoffset) + consm_i(3) = syminus(i+xoffset,j+yoffset,k+zoffset) + consm_i(4) = szminus(i+xoffset,j+yoffset,k+zoffset) + consm_i(5) = tauminus(i+xoffset,j+yoffset,k+zoffset) + +!!$ Calculate various metric terms here. +!!$ Note also need the average of the lapse at the +!!$ left and right points. + + if (flux_direction == 1) then + avg_beta = 0.5d0 * (beta1(i+xoffset,j+yoffset,k+zoffset) + & + beta1(i,j,k)) + else if (flux_direction == 2) then + avg_beta = 0.5d0 * (beta2(i+xoffset,j+yoffset,k+zoffset) + & + beta2(i,j,k)) + else if (flux_direction == 3) then + avg_beta = 0.5d0 * (beta3(i+xoffset,j+yoffset,k+zoffset) + & + beta3(i,j,k)) + else + call CCTK_ERROR("Flux direction not x,y,z") + end if + + avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset)) + + gxxh = 0.5d0 * (g11(i+xoffset,j+yoffset,k+zoffset) + g11(i,j,k)) + gxyh = 0.5d0 * (g12(i+xoffset,j+yoffset,k+zoffset) + g12(i,j,k)) + gxzh = 0.5d0 * (g13(i+xoffset,j+yoffset,k+zoffset) + g13(i,j,k)) + gyyh = 0.5d0 * (g22(i+xoffset,j+yoffset,k+zoffset) + g22(i,j,k)) + gyzh = 0.5d0 * (g23(i+xoffset,j+yoffset,k+zoffset) + g23(i,j,k)) + gzzh = 0.5d0 * (g33(i+xoffset,j+yoffset,k+zoffset) + g33(i,j,k)) + + avg_det = SPATIAL_DETERMINANT(gxxh,gxyh,gxzh,\ + gyyh,gyzh,gzzh) + sqrt_avg_det = sqrt(avg_det) + +!!$ If the Riemann problem is trivial, just calculate the fluxes from the +!!$ left state and skip to the next cell + + if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, trivial)) then + + if (flux_direction == 1) then + call num_x_flux(consp(1),consp(2),consp(3),consp(4),consp(5),& + f1(1),f1(2),f1(3),& + f1(4),f1(5),& + velxplus(i,j,k),pressplus(i,j,k),& + avg_det,avg_alp,avg_beta) + else if (flux_direction == 2) then + call num_x_flux(consp(1),consp(3),consp(4),consp(2),consp(5),& + f1(1),f1(3),f1(4),& + f1(2),f1(5),& + velyplus(i,j,k),pressplus(i,j,k),& + avg_det,avg_alp,avg_beta) + else if (flux_direction == 3) then + call num_x_flux(consp(1),consp(4),consp(2),consp(3),consp(5),& + f1(1),f1(4),f1(2),& + f1(3),f1(5),& + velzplus(i,j,k),pressplus(i,j,k),& + avg_det,avg_alp,avg_beta) + else + call CCTK_ERROR("Flux direction not x,y,z") + end if + + else !!! The end of this branch is right at the bottom of the routine + + call UpperMetric(uxxh, uxyh, uxzh, uyyh, uyzh, uzzh, & + avg_det,gxxh, gxyh, gxzh, & + gyyh, gyzh, gzzh) + + if (flux_direction == 1) then + usendh = uxxh + vXR = avg_alp*velxminus(i+xoffset,j+yoffset,k+zoffset) - avg_beta + vXL = avg_alp*velxplus (i ,j ,k ) - avg_beta + else if (flux_direction == 2) then + usendh = uyyh + vXR = avg_alp*velyminus(i+xoffset,j+yoffset,k+zoffset) - avg_beta + vXL = avg_alp*velyplus (i ,j ,k ) - avg_beta + else if (flux_direction == 3) then + usendh = uzzh + vXR = avg_alp*velzminus(i+xoffset,j+yoffset,k+zoffset) - avg_beta + vXL = avg_alp*velzplus (i ,j ,k ) - avg_beta + else + call CCTK_ERROR("Flux direction not x,y,z") + end if + +!!$ Eigenvalues and fluxes either side of the cell interface + + if (flux_direction == 1) then + if(evolve_temper.ne.1) then + call eigenvalues(GRHydro_eos_handle,& + rhominus(i+xoffset,j+yoffset,k+zoffset),& + velxminus(i+xoffset,j+yoffset,k+zoffset),& + velyminus(i+xoffset,j+yoffset,k+zoffset),& + velzminus(i+xoffset,j+yoffset,k+zoffset),& + epsminus(i+xoffset,j+yoffset,k+zoffset),& + w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),& + lamminus,gxxh,gxyh,gxzh,gyyh,& + gyzh,gzzh,& + usendh,avg_alp,avg_beta) + call eigenvalues(GRHydro_eos_handle,rhoplus(i,j,k),& + velxplus(i,j,k),velyplus(i,j,k),& + velzplus(i,j,k),epsplus(i,j,k),w_lorentzplus(i,j,k),& + lamplus,gxxh,gxyh,gxzh,gyyh,& + gyzh,gzzh,& + usendh,avg_alp,avg_beta) + else + call eigenvalues_hot(GRHydro_eos_handle,keytemp,& + i,j,k,& + rhominus(i+xoffset,j+yoffset,k+zoffset),& + velxminus(i+xoffset,j+yoffset,k+zoffset),& + velyminus(i+xoffset,j+yoffset,k+zoffset),& + velzminus(i+xoffset,j+yoffset,k+zoffset),& + epsminus(i+xoffset,j+yoffset,k+zoffset),& + tempminus(i+xoffset,j+yoffset,k+zoffset),& + y_e_minus(i+xoffset,j+yoffset,k+zoffset),& + w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),& + lamminus,gxxh,gxyh,gxzh,gyyh,& + gyzh,gzzh,& + usendh,avg_alp,avg_beta) + call eigenvalues_hot(GRHydro_eos_handle,keytemp,& + i,j,k,& + rhoplus(i,j,k),& + velxplus(i,j,k),velyplus(i,j,k),& + velzplus(i,j,k),epsplus(i,j,k), & + tempplus(i,j,k),& + y_e_plus(i,j,k),& + w_lorentzplus(i,j,k),& + lamplus,gxxh,gxyh,gxzh,gyyh,& + gyzh,gzzh,& + usendh,avg_alp,avg_beta) + endif + call num_x_flux(consp(1),consp(2),consp(3),consp(4),consp(5),& + fplus(1),fplus(2),fplus(3),fplus(4),& + fplus(5),velxplus(i,j,k),pressplus(i,j,k),& + avg_det,avg_alp,avg_beta) + call num_x_flux(consm_i(1),consm_i(2),consm_i(3),& + consm_i(4),consm_i(5),fminus(1),fminus(2),fminus(3),& + fminus(4), fminus(5),& + velxminus(i+xoffset,j+yoffset,k+zoffset),& + pressminus(i+xoffset,j+yoffset,k+zoffset),& + avg_det,avg_alp,avg_beta) + else if (flux_direction == 2) then + if(evolve_temper.ne.1) then + call eigenvalues(GRHydro_eos_handle,& + rhominus(i+xoffset,j+yoffset,k+zoffset),& + velyminus(i+xoffset,j+yoffset,k+zoffset),& + velzminus(i+xoffset,j+yoffset,k+zoffset),& + velxminus(i+xoffset,j+yoffset,k+zoffset),& + epsminus(i+xoffset,j+yoffset,k+zoffset),& + w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),& + lamminus,gyyh,gyzh,gxyh,gzzh,& + gxzh,gxxh,& + usendh,avg_alp,avg_beta) + call eigenvalues(GRHydro_eos_handle,rhoplus(i,j,k),& + velyplus(i,j,k),velzplus(i,j,k),& + velxplus(i,j,k),epsplus(i,j,k),w_lorentzplus(i,j,k),& + lamplus,gyyh,gyzh,gxyh,gzzh,& + gxzh,gxxh,& + usendh,avg_alp,avg_beta) + else + call eigenvalues_hot(GRHydro_eos_handle,keytemp,i,j,k,& + rhominus(i+xoffset,j+yoffset,k+zoffset),& + velyminus(i+xoffset,j+yoffset,k+zoffset),& + velzminus(i+xoffset,j+yoffset,k+zoffset),& + velxminus(i+xoffset,j+yoffset,k+zoffset),& + epsminus(i+xoffset,j+yoffset,k+zoffset),& + tempminus(i+xoffset,j+yoffset,k+zoffset),& + y_e_minus(i+xoffset,j+yoffset,k+zoffset),& + w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),& + lamminus,gyyh,gyzh,gxyh,gzzh,& + gxzh,gxxh,& + usendh,avg_alp,avg_beta) + call eigenvalues_hot(GRHydro_eos_handle,keytemp,i,j,k,& + rhoplus(i,j,k),& + velyplus(i,j,k),velzplus(i,j,k),& + velxplus(i,j,k),epsplus(i,j,k),& + tempplus(i,j,k),y_e_plus(i,j,k),& + w_lorentzplus(i,j,k),& + lamplus,gyyh,gyzh,gxyh,gzzh,& + gxzh,gxxh,& + usendh,avg_alp,avg_beta) + endif + call num_x_flux(consp(1),consp(3),consp(4),consp(2),consp(5),& + fplus(1),fplus(3),fplus(4),fplus(2),& + fplus(5),velyplus(i,j,k),pressplus(i,j,k),& + avg_det,avg_alp,avg_beta) + call num_x_flux(consm_i(1),consm_i(3),consm_i(4),& + consm_i(2),consm_i(5),fminus(1),fminus(3),fminus(4),& + fminus(2), fminus(5),& + velyminus(i+xoffset,j+yoffset,k+zoffset),& + pressminus(i+xoffset,j+yoffset,k+zoffset),& + avg_det,avg_alp,avg_beta) + else if (flux_direction == 3) then + if(evolve_temper.ne.1) then + call eigenvalues(GRHydro_eos_handle,& + rhominus(i+xoffset,j+yoffset,k+zoffset),& + velzminus(i+xoffset,j+yoffset,k+zoffset),& + velxminus(i+xoffset,j+yoffset,k+zoffset),& + velyminus(i+xoffset,j+yoffset,k+zoffset),& + epsminus(i+xoffset,j+yoffset,k+zoffset),& + w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),& + lamminus,gzzh,gxzh,gyzh,& + gxxh,gxyh,gyyh,& + usendh,avg_alp,avg_beta) + call eigenvalues(GRHydro_eos_handle,rhoplus(i,j,k),& + velzplus(i,j,k),velxplus(i,j,k),& + velyplus(i,j,k),epsplus(i,j,k),w_lorentzplus(i,j,k),& + lamplus,gzzh,gxzh,gyzh,& + gxxh,gxyh,gyyh,& + usendh,avg_alp,avg_beta) + else + call eigenvalues_hot(GRHydro_eos_handle,keytemp,i,j,k,& + rhominus(i+xoffset,j+yoffset,k+zoffset),& + velzminus(i+xoffset,j+yoffset,k+zoffset),& + velxminus(i+xoffset,j+yoffset,k+zoffset),& + velyminus(i+xoffset,j+yoffset,k+zoffset),& + epsminus(i+xoffset,j+yoffset,k+zoffset),& + tempminus(i+xoffset,j+yoffset,k+zoffset),& + y_e_minus(i+xoffset,j+yoffset,k+zoffset),& + w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),& + lamminus,gzzh,gxzh,gyzh,& + gxxh,gxyh,gyyh,& + usendh,avg_alp,avg_beta) + call eigenvalues_hot(GRHydro_eos_handle,keytemp,& + i,j,k,& + rhoplus(i,j,k),& + velzplus(i,j,k),velxplus(i,j,k),& + velyplus(i,j,k),epsplus(i,j,k),& + tempplus(i,j,k),y_e_plus(i,j,k),& + w_lorentzplus(i,j,k),& + lamplus,gzzh,gxzh,gyzh,& + gxxh,gxyh,gyyh,& + usendh,avg_alp,avg_beta) + endif + call num_x_flux(consp(1),consp(4),consp(2),consp(3),consp(5),& + fplus(1),fplus(4),fplus(2),fplus(3),& + fplus(5),velzplus(i,j,k),pressplus(i,j,k),avg_det,& + avg_alp,avg_beta) + call num_x_flux(consm_i(1),consm_i(4),consm_i(2),& + consm_i(3),consm_i(5),fminus(1),fminus(4),fminus(2),& + fminus(3), fminus(5),& + velzminus(i+xoffset,j+yoffset,k+zoffset),& + pressminus(i+xoffset,j+yoffset,k+zoffset),& + avg_det,avg_alp,avg_beta) + else + call CCTK_ERROR("Flux direction not x,y,z") + end if + + +!!$ we need characteristic speeds for leftmost and rightmost state only. + + ! Duez and GRHydro differ by a factor of lapse that we pull out of the + ! fluxes. To aid the overall confusion, we are going to adhere to + ! Duez' convetion for a while until we compute the fluxes. Note that + ! the contained function compute_lamstar has access to all of our + ! local variables and does indeed look at charmin and charmax. + charmin = min(lamminus(1),lamplus(1))*avg_alp + charmax = max(lamminus(5),lamplus(5))*avg_alp + call compute_stardata(charstar, pressstar,consm_i,consp,sqrt_avg_det,pressminus,pressplus,& + uxxh,uxyh,uxzh,uyyh,uyzh,uzzh,usendh,avg_beta,avg_alp,vxR,vxL,charmin,charmax,i,j,k,& + cctk_ash1,cctk_ash2,cctk_ash3,xoffset,yoffset,zoffset,flux_direction) + +!!$ Calculate flux by standard formula + + ! TODO: possibly one can safe a bit of time but oly computing charstar + ! when needed? No need to compute it when 0.lt.charmin or 0.gt.charmax + ! TODO: only compute fplus and fminus when we need them + ! TODO: clean this up. This is horrible. Consider reviving the GT C + ! HLL code for this. + if (0.lt.charmin) then + f1 = fplus + else if (0.lt.charstar) then + ! solve for state vector ULstar using (12)-(16) + ! note the change in sign in both numerator and denom + consstar = (vXL - charmin)/(charstar-charmin)*consp + ! sx, sy, sz + consstar(1+flux_direction) = consstar(1+flux_direction) + & + avg_alp*sqrt_avg_det* & + (pressplus (i ,j ,k )-pressstar) / & + (charstar-charmin) + ! tau + consstar(5) = consstar(5) + sqrt_avg_det * & + ((vXL+avg_beta)*pressplus (i ,j ,k ) - & + (charstar+avg_beta)*pressstar) / & + (charstar-charmin) + ! need to pass in GRHydro velocity computed from Duez velocity + ! which is charstar in the star region + velstar = (charstar+avg_beta)/avg_alp + if (flux_direction == 1) then + call num_x_flux(consstar(1),consstar(2),consstar(3),& + consstar(4),consstar(5),& + f1(1),f1(2),f1(3),f1(4),f1(5),& + velstar,pressstar,avg_det,avg_alp,avg_beta) + else if (flux_direction == 2) then + call num_x_flux(consstar(1),consstar(3),consstar(4),& + consstar(2),consstar(5),& + f1(1),f1(3),f1(4),f1(2),f1(5),& + velstar,pressstar,avg_det,avg_alp,avg_beta) + else if (flux_direction == 3) then + call num_x_flux(consstar(1),consstar(4),consstar(2),& + consstar(3),consstar(5),& + f1(1),f1(4),f1(2),f1(3),f1(5),& + velstar,pressstar,avg_det,avg_alp,avg_beta) + else + call CCTK_ERROR("Flux direction not x,y,z") + end if + else if (0.lt.charmax) then + ! solve for state vector URstar using (12)-(16) + ! note the change in sign in both numerator and denom + consstar = (vXR - charmax)/(charstar-charmax)*consm_i + ! sx, sy, sz + consstar(1+flux_direction) = consstar(1+flux_direction) + & + avg_alp*sqrt_avg_det* & + (pressminus(i+xoffset,j+yoffset,k+zoffset)-pressstar) / & + (charstar-charmax) + ! tau + consstar(5) = consstar(5) + sqrt_avg_det * & + ((vXR+avg_beta)*pressminus(i+xoffset,j+yoffset,k+zoffset) - & + (charstar+avg_beta)*pressstar) / & + (charstar-charmax) + ! need to pass in GRHydro velocity computed from Duez velocity + ! which is charstar in the star region + velstar = (charstar+avg_beta)/avg_alp + if (flux_direction == 1) then + call num_x_flux(consstar(1),consstar(2),consstar(3),& + consstar(4),consstar(5),& + f1(1),f1(2),f1(3),f1(4),f1(5),& + velstar,pressstar,avg_det,avg_alp,avg_beta) + else if (flux_direction == 2) then + call num_x_flux(consstar(1),consstar(3),consstar(4),& + consstar(2),consstar(5),& + f1(1),f1(3),f1(4),f1(2),f1(5),& + velstar,pressstar,avg_det,avg_alp,avg_beta) + else if (flux_direction == 3) then + call num_x_flux(consstar(1),consstar(4),consstar(2),& + consstar(3),consstar(5),& + f1(1),f1(4),f1(2),f1(3),f1(5),& + velstar,pressstar,avg_det,avg_alp,avg_beta) + else + call CCTK_ERROR("Flux direction not x,y,z") + end if + else + f1 = fminus + end if + + end if !!! The end of the SpaceMask check for a trivial RP. + + densflux(i, j, k) = f1(1) + sxflux(i, j, k) = f1(2) + syflux(i, j, k) = f1(3) + szflux(i, j, k) = f1(4) + tauflux(i, j, k) = f1(5) + + if(evolve_Y_e.ne.0) then + if (densflux(i, j, k) > 0.d0) then + Y_e_con_flux(i, j, k) = & + Y_e_plus(i, j, k) * & + densflux(i, j, k) + else + Y_e_con_flux(i, j, k) = & + Y_e_minus(i + xoffset, j + yoffset, k + zoffset) * & + densflux(i, j, k) + endif + endif + + end do + end do + end do + !$OMP END PARALLEL DO + + +end subroutine GRHydro_HLLC + File [modified]: GRHydro_RiemannSolve.F90 Delta lines: +10 -0 =================================================================== --- trunk/src/GRHydro_RiemannSolve.F90 2014-04-15 19:48:30 UTC (rev 606) +++ trunk/src/GRHydro_RiemannSolve.F90 2014-04-15 19:48:34 UTC (rev 607) @@ -46,6 +46,16 @@ end if + else if (CCTK_EQUALS(riemann_solver,"HLLC")) then + + call GRHydro_HLLC(CCTK_PASS_FTOF) + + if (evolve_tracer .ne. 0) then + + call GRHydro_HLLE_Tracer(CCTK_PASS_FTOF) + + end if + else if (CCTK_EQUALS(riemann_solver,"Roe")) then call GRHydro_RoeSolve(CCTK_PASS_FTOF) File [modified]: make.code.defn Delta lines: +1 -0 =================================================================== --- trunk/src/make.code.defn 2014-04-15 19:48:30 UTC (rev 606) +++ trunk/src/make.code.defn 2014-04-15 19:48:34 UTC (rev 607) @@ -18,6 +18,7 @@ GRHydro_Flux.F90 \ GRHydro_FluxSplit.F90 \ GRHydro_HLLE.F90 \ + GRHydro_HLLC.F90 \ GRHydro_Loop.F90 \ GRHydro_Minima.F90 \ GRHydro_Minima.cc \ Directory: /trunk/ ================== File [modified]: param.ccl Delta lines: +1 -0 =================================================================== --- trunk/param.ccl 2014-04-15 19:48:30 UTC (rev 606) +++ trunk/param.ccl 2014-04-15 19:48:34 UTC (rev 607) @@ -300,6 +300,7 @@ "Roe" :: "Standard Roe solver" "Marquina" :: "Marquina flux" "HLLE" :: "HLLE" + "HLLC" :: "HLLC" "LLF" :: "Local Lax-Friedrichs (MHD only at the moment)" } "HLLE"

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 606)
by rhaas＠tapir.caltech.edu 15 Apr '14

15 Apr '14

User: rhaas Date: 2014/04/15 02:48 PM Modified: /trunk/src/ GRHydro_Prim2Con.c Log: GRHydro: compute pressure in prim2con only where needed * change Prim2Con handling for realistic EOS: now compute pressure and cs2 only where really needed (saves some EOS calls and fixes issue with initialization with "save" values in Reconstruction") * may be useful to still set save values *everywhere* via memcopy in reconstruct rather than call TrivalReconstruct. File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Prim2Con.c Delta lines: +149 -141 =================================================================== --- trunk/src/GRHydro_Prim2Con.c 2014-04-15 19:48:28 UTC (rev 605) +++ trunk/src/GRHydro_Prim2Con.c 2014-04-15 19:48:30 UTC (rev 606) @@ -69,6 +69,8 @@ // EOS calls (now GF-wide) if(!*evolve_temper) { + // n needs to be computed using ash since ash is used when computing the + // linear index in CCTK_GFINDEX3D int n = cctk_ash[0]*cctk_ash[1]*cctk_ash[2]; int *keyerr = malloc(sizeof(*keyerr)*n); int anyerr = 0; @@ -76,13 +78,13 @@ // don't need special error handling for analytic EOS #pragma omp parallel for - for(int k=0;k<cctk_ash[2];k++) - for(int j=0;j<cctk_ash[1];j++) { + for(int k=0;k<cctk_lsh[2];k++) + for(int j=0;j<cctk_lsh[1];j++) { int i = CCTK_GFINDEX3D(cctkGH,0,j,k); - EOS_Omni_press_cs2(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,cctk_ash[0], + EOS_Omni_press_cs2(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,cctk_lsh[0], &(rhominus[i]),&(epsminus[i]),NULL,NULL,&(pressminus[i]),&(cs2minus[i]), &(keyerr[i]),&anyerr); - EOS_Omni_press_cs2(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,cctk_ash[0], + EOS_Omni_press_cs2(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,cctk_lsh[0], &(rhoplus[i]),&(epsplus[i]),NULL,NULL,&(pressplus[i]),&(cs2plus[i]), &(keyerr[i]),&anyerr); } @@ -90,87 +92,85 @@ } 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 *keyerr = malloc(sizeof(*keyerr)*n); - int anyerr = 0; int keytemp = 1; - // ensure Y_e and temperature within bounds #pragma omp parallel for - for(int i=0;i<n;i++) { + for(int i=0;i<n;i++) { // walks over slightly too many elements but cannot fail Y_e_minus[i] = MAX(MIN(Y_e_minus[i],GRHydro_Y_e_max),GRHydro_Y_e_min); Y_e_plus[i] = MAX(MIN(Y_e_plus[i],GRHydro_Y_e_max),GRHydro_Y_e_min); tempminus[i] = MIN(MAX(tempminus[i],GRHydro_hot_atmo_temp),GRHydro_max_temp); tempplus[i] = MIN(MAX(tempplus[i],GRHydro_hot_atmo_temp),GRHydro_max_temp); + keyerr[i] = 0; } // call the EOS with slices -#pragma omp parallel for - for(int k=0;k<cctk_ash[2];k++) - for(int j=0;j<cctk_ash[1];j++) { - int i = CCTK_GFINDEX3D(cctkGH,0,j,k); - EOS_Omni_press_cs2(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,cctk_ash[0], +#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++) { + 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, &rhominus[i],&epsminus[i],&tempminus[i],&Y_e_minus[i], &pressminus[i],&cs2minus[i], &keyerr[i],&anyerr); - } - - if(anyerr!=0) { -#pragma omp parallel for - for(int i=0;i<n;i++) { - if(keyerr[i] != 0) { + if(anyerr) { + for(int ii=GRHydro_stencil-1;ii<cctk_lsh[0]-GRHydro_stencil+1;ii++) { + int idx = CCTK_GFINDEX3D(cctkGH,ii,j,k); + if(keyerr[idx]!=0) { #pragma omp critical - { - CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, - "rl: %d i,x,y,z: %d %15.6E %15.6E %15.6E, keyerr: %d", - *GRHydro_reflevel, i, x[i], y[i], z[i], keyerr[i]); - CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, - "rl: %d r,t,ye: %15.6E %15.6E %15.6E, keyerr: %d", - *GRHydro_reflevel, rhominus[i], tempminus[i], Y_e_minus[i], keyerr[i]); + { + CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, + "rl: %d i,j,k,x,y,z: %d,%d,%d %15.6E %15.6E %15.6E, keyerr: %d", + *GRHydro_reflevel, ii, j, k, x[idx], y[idx], z[idx], keyerr[idx]); + CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, + "rl: %d r,t,ye: %15.6E %15.6E %15.6E, keyerr: %d", + *GRHydro_reflevel, rhominus[idx], tempminus[idx], Y_e_minus[idx], keyerr[idx]); + } + } } - } - } // for i, i<n - CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, - "Aborting!"); - } + CCTK_ERROR("Aborting!"); + } // if (anyerr) + } // loop + -#pragma omp parallel for - for(int k=0;k<cctk_ash[2];k++) - for(int j=0;j<cctk_ash[1];j++) { - int i = CCTK_GFINDEX3D(cctkGH,0,j,k); - EOS_Omni_press_cs2(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,cctk_ash[0], +#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++) { + 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, &rhoplus[i],&epsplus[i],&tempplus[i],&Y_e_plus[i], &pressplus[i],&cs2plus[i], &keyerr[i],&anyerr); - } - - - if(anyerr!=0) { -#pragma omp parallel for - for(int i=0;i<n;i++) { - if(keyerr[i] != 0) { + if(anyerr) { + for(int ii=GRHydro_stencil-1;ii<cctk_lsh[0]-GRHydro_stencil+1;ii++) { + int idx = CCTK_GFINDEX3D(cctkGH,ii,j,k); + if(keyerr[idx]!=0) { #pragma omp critical - { - CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, - "rl: %d i,x,y,z: %d %15.6E %15.6E %15.6E, keyerr: %d", - *GRHydro_reflevel, i, x[i], y[i], z[i], keyerr[i]); - CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, - "rl: %d r,t,ye: %15.6E %15.6E %15.6E, keyerr: %d", - *GRHydro_reflevel, rhoplus[i], tempplus[i], Y_e_plus[i], keyerr[i]); + { + CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, + "rl: %d i,j,k,x,y,z: %d,%d,%d %15.6E %15.6E %15.6E, keyerr: %d", + *GRHydro_reflevel, ii, j, k, x[idx], y[idx], z[idx], keyerr[idx]); + CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, + "rl: %d r,t,ye: %15.6E %15.6E %15.6E, keyerr: %d", + *GRHydro_reflevel, rhoplus[idx], tempplus[idx], Y_e_plus[idx], keyerr[idx]); + } + } } - } - } // for i, i<n - CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, - "Aborting!"); - } + CCTK_ERROR("Aborting!"); + } // if (anyerr) + } // loop free(keyerr); } 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 *keyerr = malloc(sizeof(*keyerr)*n); - int anyerr = 0; int keytemp = 0; - // ensure Y_e and temperature within bounds + // ensure Y_e and temperature within bounds #pragma omp parallel for for(int i=0;i<n;i++) { Y_e_minus[i] = MAX(MIN(Y_e_minus[i],GRHydro_Y_e_max),GRHydro_Y_e_min); @@ -178,103 +178,111 @@ tempminus[i] = MIN(MAX(tempminus[i],GRHydro_hot_atmo_temp),GRHydro_max_temp); tempplus[i] = MIN(MAX(tempplus[i],GRHydro_hot_atmo_temp),GRHydro_max_temp); temperature[i] = MIN(MAX(temperature[i],GRHydro_hot_atmo_temp),GRHydro_max_temp); + keyerr[i] = 0; } - // call the EOS with slices -#pragma omp parallel for - for(int k=0;k<cctk_ash[2];k++) - for(int j=0;j<cctk_ash[1];j++) { - int i = CCTK_GFINDEX3D(cctkGH,0,j,k); - EOS_Omni_press_cs2(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,cctk_ash[0], + // 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++) { + 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, &rhominus[i],&epsminus[i],&tempminus[i],&Y_e_minus[i], &pressminus[i],&cs2minus[i], &keyerr[i],&anyerr); - } - - if(anyerr!=0) { -#pragma omp parallel for - for(int i=0;i<n;i++) { - if(keyerr[i] != 0) { + if(anyerr) { + for(int ii=GRHydro_stencil-1;ii<cctk_lsh[0]-GRHydro_stencil+1;ii++) { + int idx = CCTK_GFINDEX3D(cctkGH,ii,j,k); + if(keyerr[idx]!=0) { #pragma omp critical - { - CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, - "rl: %d x,y,z: %15.6E %15.6E %15.6E, keyerr: %d", - *GRHydro_reflevel, x[i], y[i], z[i], keyerr[i]); - CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, - "rl: %d r,t,ye: %15.6E %15.6E %15.6E, keyerr: %d", - *GRHydro_reflevel, rhominus[i], tempminus[i], - Y_e_minus[i], keyerr[i]); - if(keyerr[i] == 668) { - // This means the temperature came back negative. - // We'll try using piecewise constant for the temperature - tempminus[i] = temperature[i]; - const int ln=1; - int lkeyerr[1]; - int lanyerr = 0; - int lkeytemp = 1; - EOS_Omni_press_cs2(*GRHydro_eos_handle,lkeytemp,GRHydro_eos_rf_prec,ln, - &rhominus[i],&epsminus[i],&tempminus[i], - &Y_e_minus[i],&pressminus[i],&cs2minus[i],lkeyerr,&lanyerr); - if(lanyerr !=0) { - CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, - "Aborting! keyerr=%d, r=%15.6E, t=%15.6E, ye=%15.6E", - lkeyerr[0],rhominus[i],tempminus[i],Y_e_minus[i]); - } - } else { - CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, - "Aborting!"); - } - } - } - } // for i, i<n - } + { + CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, + "rl: %d i,j,k,x,y,z: %d,%d,%d %15.6E %15.6E %15.6E, keyerr: %d", + *GRHydro_reflevel, ii, j, k, x[idx], y[idx], z[idx], keyerr[idx]); + CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, + "rl: %d r,t,ye: %15.6E %15.6E %15.6E, keyerr: %d", + *GRHydro_reflevel, rhominus[idx], tempminus[idx], Y_e_minus[idx], keyerr[idx]); -#pragma omp parallel for - for(int k=0;k<cctk_ash[2];k++) - for(int j=0;j<cctk_ash[1];j++) { - int i = CCTK_GFINDEX3D(cctkGH,0,j,k); - EOS_Omni_press_cs2(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,cctk_ash[0], + if(keyerr[idx] == 668) { + // This means the temperature came back negative. + // We'll try using piecewise constant for the temperature + tempminus[idx] = temperature[idx]; + const int ln=1; + int lkeyerr[1]; + int lanyerr = 0; + int lkeytemp = 1; + EOS_Omni_press_cs2(*GRHydro_eos_handle,lkeytemp,GRHydro_eos_rf_prec,ln, + &rhominus[i],&epsminus[i],&tempminus[i], + &Y_e_minus[i],&pressminus[i],&cs2minus[i],lkeyerr,&lanyerr); + if(lanyerr !=0) { + CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, + "Aborting! keyerr=%d, r=%15.6E, t=%15.6E, ye=%15.6E", + lkeyerr[0],rhominus[i],tempminus[i],Y_e_minus[i]); + } + } else { + CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, + "Aborting! keyerr=%d, r=%15.6E, t=%15.6E, ye=%15.6E", + keyerr[i],rhominus[i],tempminus[i],Y_e_minus[i]); + } + } // omp critical pragma + } // if keyerr + } // loop ii + } // if (anyerr) + } // big loop + + // 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++) { + 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, &rhoplus[i],&epsplus[i],&tempplus[i],&Y_e_plus[i], &pressplus[i],&cs2plus[i], &keyerr[i],&anyerr); - } - - if(anyerr!=0) { -#pragma omp parallel for - for(int i=0;i<n;i++) { - if(keyerr[i] != 0) { + if(anyerr) { + for(int ii=GRHydro_stencil-1;ii<cctk_lsh[0]-GRHydro_stencil+1;ii++) { + int idx = CCTK_GFINDEX3D(cctkGH,ii,j,k); + if(keyerr[idx]!=0) { #pragma omp critical - { - CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, - "rl: %d x,y,z: %15.6E %15.6E %15.6E, keyerr: %d", - *GRHydro_reflevel, x[i], y[i], z[i], keyerr[i]); - CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, - "rl: %d r,t,ye: %15.6E %15.6E %15.6E, keyerr: %d", - *GRHydro_reflevel, rhoplus[i], tempplus[i], Y_e_plus[i], keyerr[i]); - if(keyerr[i] == 668) { - // This means the temperature came back negative. - // We'll try using piecewise constant for the temperature - tempplus[i] = temperature[i]; - const int ln=1; - int lkeyerr[1]; - int lanyerr = 0; - int lkeytemp = 1; - EOS_Omni_press_cs2(*GRHydro_eos_handle,lkeytemp,GRHydro_eos_rf_prec,ln, - &rhoplus[i],&epsplus[i],&tempplus[i], - &Y_e_plus[i],&pressplus[i],&cs2plus[i],lkeyerr,&lanyerr); - if(lanyerr !=0) { + { + CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, + "rl: %d i,j,k,x,y,z: %d,%d,%d %15.6E %15.6E %15.6E, keyerr: %d", + *GRHydro_reflevel, ii, j, k, x[idx], y[idx], z[idx], keyerr[idx]); CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, - "Aborting! keyerr=%d, r=%15.6E, t=%15.6E, ye=%15.6E", - lkeyerr[0],rhoplus[i],tempplus[i],Y_e_plus[i]); - } - } else { - CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, - "Aborting!"); - } - } // end critical - } - } // for i, i<n error checking - } + "rl: %d r,t,ye: %15.6E %15.6E %15.6E, keyerr: %d", + *GRHydro_reflevel, rhoplus[idx], tempplus[idx], Y_e_plus[idx], keyerr[idx]); + + if(keyerr[idx] == 668) { + // This means the temperature came back negative. + // We'll try using piecewise constant for the temperature + tempplus[idx] = temperature[idx]; + const int ln=1; + int lkeyerr[1]; + int lanyerr = 0; + int lkeytemp = 1; + EOS_Omni_press_cs2(*GRHydro_eos_handle,lkeytemp,GRHydro_eos_rf_prec,ln, + &rhoplus[i],&epsplus[i],&tempplus[i], + &Y_e_plus[i],&pressplus[i],&cs2plus[i],lkeyerr,&lanyerr); + if(lanyerr !=0) { + CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, + "Aborting! keyerr=%d, r=%15.6E, t=%15.6E, ye=%15.6E", + lkeyerr[0],rhoplus[i],tempplus[i],Y_e_plus[i]); + } + } else { + CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, + "Aborting! keyerr=%d, r=%15.6E, t=%15.6E, ye=%15.6E", + keyerr[i],rhoplus[i],tempplus[i],Y_e_plus[i]); + } + } // 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

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 605)
by rhaas＠tapir.caltech.edu 15 Apr '14

15 Apr '14

User: rhaas Date: 2014/04/15 02:48 PM Modified: /trunk/ interface.ccl /trunk/src/ GRHydro_Prim2Con.c Log: GRHydro: start changes for improved EOS and HLLE handling File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Prim2Con.c Delta lines: +16 -16 =================================================================== --- trunk/src/GRHydro_Prim2Con.c 2014-04-15 19:48:25 UTC (rev 604) +++ trunk/src/GRHydro_Prim2Con.c 2014-04-15 19:48:28 UTC (rev 605) @@ -79,11 +79,11 @@ for(int k=0;k<cctk_ash[2];k++) for(int j=0;j<cctk_ash[1];j++) { int i = CCTK_GFINDEX3D(cctkGH,0,j,k); - EOS_Omni_press(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,cctk_ash[0], - &(rhominus[i]),&(epsminus[i]),NULL,NULL,&(pressminus[i]), + EOS_Omni_press_cs2(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,cctk_ash[0], + &(rhominus[i]),&(epsminus[i]),NULL,NULL,&(pressminus[i]),&(cs2minus[i]), &(keyerr[i]),&anyerr); - EOS_Omni_press(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,cctk_ash[0], - &(rhoplus[i]),&(epsplus[i]),NULL,NULL,&(pressplus[i]), + EOS_Omni_press_cs2(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,cctk_ash[0], + &(rhoplus[i]),&(epsplus[i]),NULL,NULL,&(pressplus[i]),&(cs2plus[i]), &(keyerr[i]),&anyerr); } free(keyerr); @@ -108,9 +108,9 @@ for(int k=0;k<cctk_ash[2];k++) for(int j=0;j<cctk_ash[1];j++) { int i = CCTK_GFINDEX3D(cctkGH,0,j,k); - EOS_Omni_press(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,cctk_ash[0], + EOS_Omni_press_cs2(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,cctk_ash[0], &rhominus[i],&epsminus[i],&tempminus[i],&Y_e_minus[i], - &pressminus[i], + &pressminus[i],&cs2minus[i], &keyerr[i],&anyerr); } @@ -137,9 +137,9 @@ for(int k=0;k<cctk_ash[2];k++) for(int j=0;j<cctk_ash[1];j++) { int i = CCTK_GFINDEX3D(cctkGH,0,j,k); - EOS_Omni_press(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,cctk_ash[0], + EOS_Omni_press_cs2(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,cctk_ash[0], &rhoplus[i],&epsplus[i],&tempplus[i],&Y_e_plus[i], - &pressplus[i], + &pressplus[i],&cs2plus[i], &keyerr[i],&anyerr); } @@ -185,9 +185,9 @@ for(int k=0;k<cctk_ash[2];k++) for(int j=0;j<cctk_ash[1];j++) { int i = CCTK_GFINDEX3D(cctkGH,0,j,k); - EOS_Omni_press(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,cctk_ash[0], + EOS_Omni_press_cs2(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,cctk_ash[0], &rhominus[i],&epsminus[i],&tempminus[i],&Y_e_minus[i], - &pressminus[i], + &pressminus[i],&cs2minus[i], &keyerr[i],&anyerr); } @@ -212,9 +212,9 @@ int lkeyerr[1]; int lanyerr = 0; int lkeytemp = 1; - EOS_Omni_press(*GRHydro_eos_handle,lkeytemp,GRHydro_eos_rf_prec,ln, + EOS_Omni_press_cs2(*GRHydro_eos_handle,lkeytemp,GRHydro_eos_rf_prec,ln, &rhominus[i],&epsminus[i],&tempminus[i], - &Y_e_minus[i],&pressminus[i],lkeyerr,&lanyerr); + &Y_e_minus[i],&pressminus[i],&cs2minus[i],lkeyerr,&lanyerr); if(lanyerr !=0) { CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, "Aborting! keyerr=%d, r=%15.6E, t=%15.6E, ye=%15.6E", @@ -233,9 +233,9 @@ for(int k=0;k<cctk_ash[2];k++) for(int j=0;j<cctk_ash[1];j++) { int i = CCTK_GFINDEX3D(cctkGH,0,j,k); - EOS_Omni_press(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,cctk_ash[0], + EOS_Omni_press_cs2(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,cctk_ash[0], &rhoplus[i],&epsplus[i],&tempplus[i],&Y_e_plus[i], - &pressplus[i], + &pressplus[i],&cs2plus[i], &keyerr[i],&anyerr); } @@ -259,9 +259,9 @@ int lkeyerr[1]; int lanyerr = 0; int lkeytemp = 1; - EOS_Omni_press(*GRHydro_eos_handle,lkeytemp,GRHydro_eos_rf_prec,ln, + EOS_Omni_press_cs2(*GRHydro_eos_handle,lkeytemp,GRHydro_eos_rf_prec,ln, &rhoplus[i],&epsplus[i],&tempplus[i], - &Y_e_plus[i],&pressplus[i],lkeyerr,&lanyerr); + &Y_e_plus[i],&pressplus[i],&cs2plus[i],lkeyerr,&lanyerr); if(lanyerr !=0) { CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, "Aborting! keyerr=%d, r=%15.6E, t=%15.6E, ye=%15.6E", Directory: /trunk/ ================== File [modified]: interface.ccl Delta lines: +17 -2 =================================================================== --- trunk/interface.ccl 2014-04-15 19:48:25 UTC (rev 604) +++ trunk/interface.ccl 2014-04-15 19:48:28 UTC (rev 605) @@ -286,7 +286,22 @@ USES FUNCTION EOS_Omni_press +void FUNCTION EOS_Omni_press_cs2(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_REAL OUT ARRAY cs2, \ + CCTK_INT OUT ARRAY keyerr, \ + CCTK_INT OUT anyerr) +USES FUNCTION EOS_Omni_press_cs2 + + void FUNCTION EOS_Omni_DPressByDEps(CCTK_INT IN eoskey, \ CCTK_INT IN havetemp, \ CCTK_REAL IN rf_precision, \ @@ -455,8 +470,8 @@ real GRHydro_prim_bext type = GF Timelevels = 1 tags='Prolongation="None" checkpoint="no"' { rhoplus, velxplus, velyplus, velzplus, pressplus, epsplus, w_lorentzplus, - rhominus, velxminus, velyminus, velzminus, pressminus, epsminus, - w_lorentzminus + cs2plus, rhominus, velxminus, velyminus, velzminus, pressminus, epsminus, + w_lorentzminus, cs2minus } "Primitive variables extended to the cell boundaries" int GRHydro_scalars type = SCALAR tags='checkpoint="no"'

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 604)
by rhaas＠tapir.caltech.edu 15 Apr '14

15 Apr '14

User: rhaas Date: 2014/04/15 02:48 PM Modified: /trunk/src/ GRHydro_Prim2Con.c Log: GRHydro: change prim2con EOS calls to calls for x-slices File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Prim2Con.c Delta lines: +51 -15 =================================================================== --- trunk/src/GRHydro_Prim2Con.c 2014-04-15 19:48:22 UTC (rev 603) +++ trunk/src/GRHydro_Prim2Con.c 2014-04-15 19:48:25 UTC (rev 604) @@ -75,12 +75,17 @@ int keytemp = 0; // don't need special error handling for analytic EOS - EOS_Omni_pressOMP(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n, - rhominus,epsminus,NULL,NULL,pressminus,keyerr,&anyerr); - - EOS_Omni_pressOMP(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n, - rhoplus,epsplus,NULL,NULL,pressplus,keyerr,&anyerr); - +#pragma omp parallel for + for(int k=0;k<cctk_ash[2];k++) + for(int j=0;j<cctk_ash[1];j++) { + int i = CCTK_GFINDEX3D(cctkGH,0,j,k); + EOS_Omni_press(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,cctk_ash[0], + &(rhominus[i]),&(epsminus[i]),NULL,NULL,&(pressminus[i]), + &(keyerr[i]),&anyerr); + EOS_Omni_press(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,cctk_ash[0], + &(rhoplus[i]),&(epsplus[i]),NULL,NULL,&(pressplus[i]), + &(keyerr[i]),&anyerr); + } free(keyerr); } else { if(reconstruct_temper) { @@ -98,9 +103,17 @@ tempplus[i] = MIN(MAX(tempplus[i],GRHydro_hot_atmo_temp),GRHydro_max_temp); } - EOS_Omni_pressOMP(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n, - rhominus,epsminus,tempminus,Y_e_minus,pressminus,keyerr,&anyerr); - + // call the EOS with slices +#pragma omp parallel for + for(int k=0;k<cctk_ash[2];k++) + for(int j=0;j<cctk_ash[1];j++) { + int i = CCTK_GFINDEX3D(cctkGH,0,j,k); + EOS_Omni_press(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,cctk_ash[0], + &rhominus[i],&epsminus[i],&tempminus[i],&Y_e_minus[i], + &pressminus[i], + &keyerr[i],&anyerr); + } + if(anyerr!=0) { #pragma omp parallel for for(int i=0;i<n;i++) { @@ -120,9 +133,17 @@ "Aborting!"); } - EOS_Omni_pressOMP(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n, - rhoplus,epsplus,tempplus,Y_e_plus,pressplus,keyerr,&anyerr); +#pragma omp parallel for + for(int k=0;k<cctk_ash[2];k++) + for(int j=0;j<cctk_ash[1];j++) { + int i = CCTK_GFINDEX3D(cctkGH,0,j,k); + EOS_Omni_press(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,cctk_ash[0], + &rhoplus[i],&epsplus[i],&tempplus[i],&Y_e_plus[i], + &pressplus[i], + &keyerr[i],&anyerr); + } + if(anyerr!=0) { #pragma omp parallel for for(int i=0;i<n;i++) { @@ -159,8 +180,16 @@ temperature[i] = MIN(MAX(temperature[i],GRHydro_hot_atmo_temp),GRHydro_max_temp); } - EOS_Omni_pressOMP(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n, - rhominus,epsminus,tempminus,Y_e_minus,pressminus,keyerr,&anyerr); + // call the EOS with slices +#pragma omp parallel for + for(int k=0;k<cctk_ash[2];k++) + for(int j=0;j<cctk_ash[1];j++) { + int i = CCTK_GFINDEX3D(cctkGH,0,j,k); + EOS_Omni_press(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,cctk_ash[0], + &rhominus[i],&epsminus[i],&tempminus[i],&Y_e_minus[i], + &pressminus[i], + &keyerr[i],&anyerr); + } if(anyerr!=0) { #pragma omp parallel for @@ -200,8 +229,15 @@ } // for i, i<n } - EOS_Omni_pressOMP(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n, - rhoplus,epsplus,tempplus,Y_e_plus,pressplus,keyerr,&anyerr); +#pragma omp parallel for + for(int k=0;k<cctk_ash[2];k++) + for(int j=0;j<cctk_ash[1];j++) { + int i = CCTK_GFINDEX3D(cctkGH,0,j,k); + EOS_Omni_press(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,cctk_ash[0], + &rhoplus[i],&epsplus[i],&tempplus[i],&Y_e_plus[i], + &pressplus[i], + &keyerr[i],&anyerr); + } if(anyerr!=0) { #pragma omp parallel for

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 603)
by rhaas＠tapir.caltech.edu 15 Apr '14

15 Apr '14

User: rhaas Date: 2014/04/15 02:48 PM Removed: /trunk/src/ GRHydro_PPM.c Log: GRHydro: remove unused source file GRHydro_PPM.c File Changes: Directory: /trunk/src/ ====================== File [removed]: GRHydro_PPM.c Delta lines: +0 -429 =================================================================== --- trunk/src/GRHydro_PPM.c 2014-04-15 19:48:21 UTC (rev 602) +++ trunk/src/GRHydro_PPM.c 2014-04-15 19:48:22 UTC (rev 603) @@ -1,429 +0,0 @@ -#include <assert.h> -#include <math.h> -#include <stdio.h> -#include <stdlib.h> - -#if 1 -#include "cctk.h" -#include "cctk_Parameters.h" -#include "cctk_Arguments.h" -#include "cctk_Functions.h" -#endif - -#include "GRHydro_PPM.h" - - -void SimplePPM_1d_c(int* restrict nx, - double* restrict dx, - double* restrict rho, - double* restrict velx, - double* restrict vely, - double* restrict velz, - double* restrict eps, - double* restrict press, - double* restrict rhominus, - double* restrict velxminus, - double* restrict velyminus, - double* restrict velzminus, - double* restrict epsminus, - double* restrict rhoplus, - double* restrict velxplus, - double* restrict velyplus, - double* restrict velzplus, - double* restrict epsplus) - -{ - double drho[*nx], dpress[*nx], deps[*nx], d2rho[*nx]; - double dmrho[*nx], dmpress[*nx], dmeps[*nx]; - double dvelx[*nx], dvely[*nx], dvelz[*nx]; - double dmvelx[*nx], dmvely[*nx], dmvelz[*nx]; - double tilde_flatten[*nx]; - - const double onesixth = 1.0/6.0; - - 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]); - } - - 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); - } - - for(int i=1;i<*nx-1;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; - } - - 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]; - } - - const double ppm_epsilon_shock = 0.01; - const double ppm_eta1 = 20.0; - const double ppm_eta2 = 0.05; - const double ppm_k0 = 0.2; - const double ppm_epsilon = 0.33; - const double ppm_omega1 = 0.75; - const double ppm_omega2 = 10.0; - const double ppm_small = 1.0e-7; - - // steepening - 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]; - } // 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); - } - - return; -} - -CCTK_FCALL void CCTK_FNAME(PPMc)(int* restrict nx, - double* restrict dx, - double* restrict rho, - double* restrict velx, - double* restrict vely, - double* restrict velz, - double* restrict eps, - double* restrict press, - double* restrict rhominus, - double* restrict velxminus, - double* restrict velyminus, - double* restrict velzminus, - double* restrict epsminus, - double* restrict rhoplus, - double* restrict velxplus, - double* restrict velyplus, - double* restrict velzplus, - double* restrict epsplus) { - -void SimplePPM_1d_c(nx, dx, - rho, - velx, - vely, - velz, - eps, - press, - rhominus, - velxminus, - velyminus, - velzminus, - epsminus, - rhoplus, - velxplus, - velyplus, - velzplus, - epsplus); - -} - - -void SimplePPM_1d_tyc(int* restrict nx, - double* restrict dx, - double* restrict rho, - double* restrict velx, - double* restrict vely, - double* restrict velz, - double* restrict eps, - double* restrict temp, - double* restrict ye, - double* restrict press, - double* restrict rhominus, - double* restrict velxminus, - double* restrict velyminus, - double* restrict velzminus, - double* restrict epsminus, - double * restrict tempminus, - double * restrict yeminus, - double* restrict rhoplus, - double* restrict velxplus, - double* restrict velyplus, - double* restrict velzplus, - double* restrict epsplus, - double* restrict tempplus, - double* restrict yeplus) - -{ - double drho[*nx], dpress[*nx], deps[*nx], d2rho[*nx]; - double dtemp[*nx], dye[*nx]; - double dmrho[*nx], dmpress[*nx], dmeps[*nx]; - double dmtemp[*nx], dmye[*nx]; - double dvelx[*nx], dvely[*nx], dvelz[*nx]; - double dmvelx[*nx], dmvely[*nx], dmvelz[*nx]; - double tilde_flatten[*nx]; - - const double onesixth = 1.0/6.0; - - 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]); - dtemp[i] = 0.5 * (temp[i+1]-temp[i-1]); - dye[i] = 0.5 * (ye[i+1]-ye[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]); - } - - for(int i=1;i<*nx-1;i++) { - steep(rho,drho,dmrho,i); - steep(eps,deps,dmeps,i); - steep(temp,dtemp,dmtemp,i); - steep(ye,dye,dmye,i); - steep(velx,dvelx,dmvelx,i); - steep(vely,dvely,dmvely,i); - steep(velz,dvelz,dmvelz,i); - } - - for(int i=1;i<*nx-1;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; - 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; - 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; - } - - for(int i=1;i<*nx-2;i++) { - rhoplus[i] = approx_at_cell_interface(rho,i); - epsplus[i] = approx_at_cell_interface(eps,i); - tempplus[i] = approx_at_cell_interface(temp,i); - yeplus[i] = approx_at_cell_interface(ye,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]; - tempminus[i+1] = tempplus[i]; - yeminus[i+1] = yeplus[i]; - velxminus[i+1] = velxplus[i]; - velyminus[i+1] = velyplus[i]; - velzminus[i+1] = velzplus[i]; - } - - const double ppm_epsilon_shock = 0.01; - const double ppm_eta1 = 20.0; - const double ppm_eta2 = 0.05; - const double ppm_k0 = 0.2; - const double ppm_epsilon = 0.33; - const double ppm_omega1 = 0.75; - const double ppm_omega2 = 10.0; - const double ppm_small = 1.0e-7; - - // steepening - 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]; - 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]; - 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]; - } // flattening - - for(int i=2;i<*nx-2;i++) { - monotonize(rhominus,rho,rhoplus,i); - monotonize(epsminus,eps,epsplus,i); - monotonize(tempminus,temp,tempplus,i); - monotonize(yeminus,ye,yeplus,i); - monotonize(velxminus,velx,velxplus,i); - monotonize(velyminus,vely,velyplus,i); - monotonize(velzminus,velz,velzplus,i); - } - - return; -} - - - -CCTK_FCALL void CCTK_FNAME(PPMtyc)(int* restrict nx, - double* restrict dx, - double* restrict rho, - double* restrict velx, - double* restrict vely, - double* restrict velz, - double* restrict eps, - double* restrict temp, - double* restrict ye, - double* restrict press, - double* restrict rhominus, - double* restrict velxminus, - double* restrict velyminus, - double* restrict velzminus, - double* restrict epsminus, - double* restrict tempminus, - double* restrict yeminus, - double* restrict rhoplus, - double* restrict velxplus, - double* restrict velyplus, - double* restrict velzplus, - double* restrict epsplus, - double* restrict tempplus, - double* restrict yeplus) { - -void SimplePPM_1d_tyc(nx, dx, - rho, - velx, - vely, - velz, - eps, - temp, - ye, - press, - rhominus, - velxminus, - velyminus, - velzminus, - epsminus, - tempminus, - yeminus, - rhoplus, - velxplus, - velyplus, - velzplus, - epsplus, - tempplus, - yeplus); -} -

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 602)
by rhaas＠tapir.caltech.edu 15 Apr '14

15 Apr '14

User: rhaas Date: 2014/04/15 02:48 PM Modified: /trunk/src/ GRHydro_LastMoLPostStep.c Log: GRHydro: Simplify logical expression File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_LastMoLPostStep.c Delta lines: +1 -4 =================================================================== --- trunk/src/GRHydro_LastMoLPostStep.c 2014-04-15 19:48:18 UTC (rev 601) +++ trunk/src/GRHydro_LastMoLPostStep.c 2014-04-15 19:48:21 UTC (rev 602) @@ -28,10 +28,7 @@ // If counter becomes zero, the only thing left to do is to call PostStep! // NOTE: MoL's counter is uninitialised during initial data setup! - if (cctk_iteration == 0 || *MoL_Intermediate_Step == 0) - *InLastMoLPostStep = 1; - else - *InLastMoLPostStep = 0; + *InLastMoLPostStep = cctk_iteration == 0 || *MoL_Intermediate_Step == 0; } void GRHydro_ClearLastMoLPostStep(CCTK_ARGUMENTS)

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 601)
by rhaas＠tapir.caltech.edu 15 Apr '14

15 Apr '14

User: rhaas Date: 2014/04/15 02:48 PM Modified: /trunk/src/ GRHydro_HLLE.F90, GRHydro_HLLE_AM.F90 Log: GRHydro: add STOP statement after CCTK_ERROR this is to let the compiler know that CCTK_ERROR never returns. File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_HLLE.F90 Delta lines: +1 -0 =================================================================== --- trunk/src/GRHydro_HLLE.F90 2014-04-15 19:48:15 UTC (rev 600) +++ trunk/src/GRHydro_HLLE.F90 2014-04-15 19:48:18 UTC (rev 601) @@ -529,6 +529,7 @@ etaY(i,j-1,k+1, vup, eos_c)) else call CCTK_ERROR("Flux direction not x,y,z") + STOP end if ! modify eigenvalues of Roe's matrix by computed H viscosity File [modified]: GRHydro_HLLE_AM.F90 Delta lines: +1 -0 =================================================================== --- trunk/src/GRHydro_HLLE_AM.F90 2014-04-15 19:48:15 UTC (rev 600) +++ trunk/src/GRHydro_HLLE_AM.F90 2014-04-15 19:48:18 UTC (rev 601) @@ -433,6 +433,7 @@ else call CCTK_ERROR("Flux direction not x,y,z") + STOP end if

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 600)
by rhaas＠tapir.caltech.edu 15 Apr '14

15 Apr '14

User: rhaas Date: 2014/04/15 02:48 PM Modified: /trunk/src/ GRHydro_HLLE_AM.F90 Log: GRHydro: remove divergence cleaning code from Avec evolution File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_HLLE_AM.F90 Delta lines: +2 -123 =================================================================== --- trunk/src/GRHydro_HLLE_AM.F90 2014-04-15 19:48:05 UTC (rev 599) +++ trunk/src/GRHydro_HLLE_AM.F90 2014-04-15 19:48:15 UTC (rev 600) @@ -53,9 +53,6 @@ CCTK_REAL :: Bvecxlowp,Bvecxlowm,Bvecylowp,Bvecylowm,Bveczlowp,Bveczlowm CCTK_REAL :: pressstarp,pressstarm,velxlowp,velxlowm,velylowp,velylowm,velzlowp,velzlowm - CCTK_REAL :: psidcp, psidcm, psidcf, psidcdiff, psidcfp, psidcfm - CCTK_REAL :: charmax_dc, charmin_dc, charpm_dc - CCTK_INT :: type_bits, trivial CCTK_REAL :: xtemp @@ -118,7 +115,7 @@ ! constraint transport needs to be able to average fluxes in the directions ! other that flux_direction - !$OMP PARALLEL DO PRIVATE(k,j,i,f1,lamminus,lamplus,cons_p,cons_m,fplus,fminus,qdiff,psidcf,psidcp,psidcm,prim_p,prim_m,& + !$OMP PARALLEL DO PRIVATE(k,j,i,f1,lamminus,lamplus,cons_p,cons_m,fplus,fminus,qdiff,prim_p,prim_m,& !$OMP avg_betax,avg_betay,avg_betaz,avg_beta,avg_alp,& !$OMP gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,avg_det,sdet,uxxh,uxyh,uxzh,uyyh,uyzh,uzzh,& !$OMP vxtp,vxtm,vytp,vytm,vztp,vztm,& @@ -127,7 +124,7 @@ !$OMP velzlowp,velzlowm,Bveczlowp,Bveczlowm,& !$OMP bdotvp,bdotvm,b2p,b2m,v2p,v2m,wp,wm,& !$OMP bxlowp,bxlowm,bylowp,bylowm,bzlowp,bzlowm,& - !$OMP rhoenth_p,rhoenth_m,ab0p,ab0m,vA2p,vA2m,pressstarp,pressstarm,usendh,psidcdiff,psidcfp,psidcfm,charmin,charmax,chartop,charpm,m,xtemp) + !$OMP rhoenth_p,rhoenth_m,ab0p,ab0m,vA2p,vA2m,pressstarp,pressstarm,usendh,charmin,charmax,chartop,charpm,m,xtemp) do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil + transport_constraints*(1-zoffset) do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil + transport_constraints*(1-yoffset) do i = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + transport_constraints*(1-xoffset) @@ -141,11 +138,6 @@ fminus = 0.d0 qdiff = 0.d0 - if(clean_divergence.ne.0) then - psidcp = 0.d0 - psidcm = 0.d0 - endif - !!$ Set the left (p for plus) and right (m_i for minus, i+1) states cons_p(1) = densplus(i,j,k) @@ -188,11 +180,6 @@ prim_m(9) = Bvecyminus(i+xoffset,j+yoffset,k+zoffset) prim_m(10)= Bveczminus(i+xoffset,j+yoffset,k+zoffset) - if(clean_divergence.ne.0) then - psidcp = psidcplus(i,j,k) - psidcm = psidcminus(i+xoffset,j+yoffset,k+zoffset) - endif - !!$ Calculate various metric terms here. !!$ Note also need the average of the lapse at the !!$ left and right points. @@ -277,12 +264,6 @@ f1(1),f1(2),f1(3),f1(4),f1(5),f1(6),f1(7),f1(8),& vxtp,vytp,vztp,pressstarp,bxlowp,bylowp,bzlowp,ab0p,wp, & avg_det,avg_alp,avg_beta) - if(clean_divergence.ne.0) then - f1(6)=f1(6) + 1.0d0*sdet*uxxh*psidcp - cons_p(6)*avg_betax/avg_alp - f1(7)=f1(7) + 1.0d0*sdet*uxyh*psidcp - cons_p(6)*avg_betay/avg_alp - f1(8)=f1(8) + 1.0d0*sdet*uxzh*psidcp - cons_p(6)*avg_betaz/avg_alp - psidcf = cons_p(6)/sdet-psidcp*avg_betax/avg_alp - endif else if (flux_direction == 2) then call num_x_fluxM(cons_p(1),cons_p(3),cons_p(4),cons_p(2),cons_p(5),& @@ -290,12 +271,6 @@ f1(1),f1(3),f1(4),f1(2),f1(5),f1(7),f1(8),f1(6),& vytp,vztp,vxtp,pressstarp,bylowp,bzlowp,bxlowp,ab0p,wp, & avg_det,avg_alp,avg_beta) - if(clean_divergence.ne.0) then - f1(6)=f1(6) + 1.0d0*sdet*uxyh*psidcp - cons_p(7)*avg_betax/avg_alp - f1(7)=f1(7) + 1.0d0*sdet*uyyh*psidcp - cons_p(7)*avg_betay/avg_alp - f1(8)=f1(8) + 1.0d0*sdet*uyzh*psidcp - cons_p(7)*avg_betaz/avg_alp - psidcf = cons_p(7)/sdet-psidcp*avg_betay/avg_alp - endif else if (flux_direction == 3) then call num_x_fluxM(cons_p(1),cons_p(4),cons_p(2),cons_p(3),cons_p(5),& @@ -303,12 +278,6 @@ f1(1),f1(4),f1(2),f1(3),f1(5),f1(8),f1(6),f1(7), & vztp,vxtp,vytp,pressstarp,bzlowp,bxlowp,bylowp,ab0p,wp, & avg_det,avg_alp,avg_beta) - if(clean_divergence.ne.0) then - f1(6)=f1(6) + 1.0d0*sdet*uxzh*psidcp - cons_p(8)*avg_betax/avg_alp - f1(7)=f1(7) + 1.0d0*sdet*uyzh*psidcp - cons_p(8)*avg_betay/avg_alp - f1(8)=f1(8) + 1.0d0*sdet*uzzh*psidcp - cons_p(8)*avg_betaz/avg_alp - psidcf = cons_p(8)/sdet-psidcp*avg_betaz/avg_alp - endif else call CCTK_ERROR("Flux direction not x,y,z") @@ -337,10 +306,6 @@ qdiff(7) = cons_m(7) - cons_p(7) qdiff(8) = cons_m(8) - cons_p(8) - if (clean_divergence.ne.0) then - psidcdiff = psidcm - psidcp - endif - !!$ Eigenvalues and fluxes either side of the cell interface if (flux_direction == 1) then @@ -384,17 +349,6 @@ vxtm,vytm,vztm,pressstarm,bxlowm,bylowm,bzlowm,ab0m,wm, & avg_det,avg_alp,avg_beta) - if(clean_divergence.ne.0) then - fminus(6)=fminus(6) + 1.0d0*sdet*uxxh*psidcm - cons_m(6)*avg_betax/avg_alp - fminus(7)=fminus(7) + 1.0d0*sdet*uxyh*psidcm - cons_m(6)*avg_betay/avg_alp - fminus(8)=fminus(8) + 1.0d0*sdet*uxzh*psidcm - cons_m(6)*avg_betaz/avg_alp - fplus(6)=fplus(6) + 1.0d0*sdet*uxxh*psidcp - cons_p(6)*avg_betax/avg_alp - fplus(7)=fplus(7) + 1.0d0*sdet*uxyh*psidcp - cons_p(6)*avg_betay/avg_alp - fplus(8)=fplus(8) + 1.0d0*sdet*uxzh*psidcp - cons_p(6)*avg_betaz/avg_alp - psidcfp = cons_p(6)/sdet-avg_betax*psidcp/avg_alp - psidcfm = cons_m(6)/sdet-avg_betax*psidcm/avg_alp - endif - else if (flux_direction == 2) then if(evolve_temper.ne.1) then call eigenvaluesM(GRHydro_eos_handle,& @@ -436,17 +390,6 @@ vytm,vztm,vxtm,pressstarm,bylowm,bzlowm,bxlowm,ab0m,wm, & avg_det,avg_alp,avg_beta) - if(clean_divergence.ne.0) then - fminus(6)=fminus(6) + 1.0d0*sdet*uxyh*psidcm - cons_m(7)*avg_betax/avg_alp - fminus(7)=fminus(7) + 1.0d0*sdet*uyyh*psidcm - cons_m(7)*avg_betay/avg_alp - fminus(8)=fminus(8) + 1.0d0*sdet*uyzh*psidcm - cons_m(7)*avg_betaz/avg_alp - fplus(6)=fplus(6) + 1.0d0*sdet*uxyh*psidcp - cons_p(7)*avg_betax/avg_alp - fplus(7)=fplus(7) + 1.0d0*sdet*uyyh*psidcp - cons_p(7)*avg_betay/avg_alp - fplus(8)=fplus(8) + 1.0d0*sdet*uyzh*psidcp - cons_p(7)*avg_betaz/avg_alp - psidcfp = cons_p(7)/sdet-avg_betay*psidcp/avg_alp - psidcfm = cons_m(7)/sdet-avg_betay*psidcm/avg_alp - endif - else if (flux_direction == 3) then if(evolve_temper.ne.1) then call eigenvaluesM(GRHydro_eos_handle,& @@ -488,17 +431,6 @@ vztm,vxtm,vytm,pressstarm,bzlowm,bxlowm,bylowm,ab0m,wm, & avg_det,avg_alp,avg_beta) - if(clean_divergence.ne.0) then - fminus(6)=fminus(6) + 1.0d0*sdet*uxzh*psidcm - cons_m(8)*avg_betax/avg_alp - fminus(7)=fminus(7) + 1.0d0*sdet*uyzh*psidcm - cons_m(8)*avg_betay/avg_alp - fminus(8)=fminus(8) + 1.0d0*sdet*uzzh*psidcm - cons_m(8)*avg_betaz/avg_alp - fplus(6)=fplus(6) + 1.0d0*sdet*uxzh*psidcp - cons_p(8)*avg_betax/avg_alp - fplus(7)=fplus(7) + 1.0d0*sdet*uyzh*psidcp - cons_p(8)*avg_betay/avg_alp - fplus(8)=fplus(8) + 1.0d0*sdet*uzzh*psidcp - cons_p(8)*avg_betaz/avg_alp - psidcfp = cons_p(8)/sdet-avg_betaz*psidcp/avg_alp - psidcfm = cons_m(8)/sdet-avg_betaz*psidcm/avg_alp - endif - else call CCTK_ERROR("Flux direction not x,y,z") end if @@ -533,55 +465,6 @@ end do - if(clean_divergence.ne.0) then - - psidcdiff = psidcm - psidcp - select case(whichpsidcspeed) - case(0) - if (HLLE) then - psidcf = (charmax * psidcfp - charmin * psidcfm + & - charmax * charmin * psidcdiff) / charpm - else if (LLF) then - psidcf = 0.5d0 * (psidcfp + psidcfm - chartop * psidcdiff) - end if - case(1) - !!$ Wavespeeds for psidc are +/-c, not Fast Magnetosonic? - !!$ psidcf = 0.5d0 * (1.d0 * psidcfp - (-1.d0) * psidcfm + & - !!$ 1.d0 * (-1.d0) * psidcdiff) - - !!$ The fastest speed for psidc comes from the condition - !!$ that the normal vector to the characteristic hypersurface - !!$ be spacelike (Eq. 60 of Anton et al.) - - charmax_dc = sqrt(usendh) - avg_beta/avg_alp - charmin_dc = -1.d0*sqrt(usendh) - avg_beta/avg_alp - charpm_dc = charmax_dc - charmin_dc - - psidcf = (charmax_dc * psidcfp - charmin_dc * psidcfm + & - charmax_dc * charmin_dc * psidcdiff) / charpm_dc - - if(decouple_normal_Bfield .ne. 0) then ! same expression for HLLE and LLF - !!$ B^i field decouples from the others and has same propagation - !!$ speed as divergence -null direction, - !!$ \pm sqrt(g^{xx}} - beta^x/alpha - f1(5+flux_direction) = (charmax_dc * fplus(5+flux_direction) & - - charmin_dc * fminus(5+flux_direction) + & - charmax_dc * charmin_dc * qdiff(5+flux_direction)) / charpm_dc - end if - case(2) - charmax = setcharmax - charmin = setcharmin - if (HLLE) then - psidcf = (charmax * psidcfp - charmin * psidcfm + & - charmax * charmin * psidcdiff) / charpm - else if (LLF) then - chartop = max(-charmin,charmax) - psidcf = 0.5d0 * (psidcfp + psidcfm - chartop * psidcdiff) - end if - end select - endif - - end if !!! The end of the SpaceMask check for a trivial RP. densflux(i, j, k) = f1(1) @@ -602,10 +485,6 @@ Aveczflux(i,j,k) = 0.d0 end if - if(clean_divergence.ne.0) then - psidcflux(i,j,k) = psidcf - endif - if(evolve_Y_e.ne.0) then if (densflux(i, j, k) > 0.d0) then Y_e_con_flux(i, j, k) = &

1 0

Jump to page:

2026

2025

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

Commits April 2014