User: zetienne Date: 2014/04/17 12:15 PM
Modified: /trunk/src/ compute_B_and_Bstagger_from_A.C, driver_conserv_to_prims.C, reconstruct_set_of_prims_PPM.C
Log: Adjust sensitivity of shock detection mechanism. Reduces explosion in roundoff-error differences between different codes.
File Changes:
Directory: /trunk/src/ ======================
File [modified]: compute_B_and_Bstagger_from_A.C Delta lines: +1 -1 =================================================================== --- trunk/src/compute_B_and_Bstagger_from_A.C 2014-04-16 16:00:47 UTC (rev 25) +++ trunk/src/compute_B_and_Bstagger_from_A.C 2014-04-17 17:15:36 UTC (rev 26) @@ -180,7 +180,7 @@ IllinoisGRMHD_set_symmetry_gzs_staggered(cctkGH,cctk_lsh,x,y,z,Symmetry, Bx , gridfunc_syms_Bx,0,0,0); IllinoisGRMHD_set_symmetry_gzs_staggered(cctkGH,cctk_lsh,x,y,z,Symmetry, Bx_stagger, gridfunc_syms_Bx,1,0,0); double gridfunc_syms_By[3] = { 1,-1,-Sym_Bz}; - IllinoisGRMHD_set_symmetry_gzs_staggered(cctkGH,cctk_lsh,x,y,z,Symmetry, By , gridfunc_syms_Bx,0,0,0); + IllinoisGRMHD_set_symmetry_gzs_staggered(cctkGH,cctk_lsh,x,y,z,Symmetry, By , gridfunc_syms_By,0,0,0); IllinoisGRMHD_set_symmetry_gzs_staggered(cctkGH,cctk_lsh,x,y,z,Symmetry, By_stagger, gridfunc_syms_By,0,1,0); double gridfunc_syms_Bz[3] = { 1, 1, Sym_Bz}; IllinoisGRMHD_set_symmetry_gzs_staggered(cctkGH,cctk_lsh,x,y,z,Symmetry, Bz , gridfunc_syms_Bz,0,0,0);
File [modified]: driver_conserv_to_prims.C Delta lines: +3 -3 =================================================================== --- trunk/src/driver_conserv_to_prims.C 2014-04-16 16:00:47 UTC (rev 25) +++ trunk/src/driver_conserv_to_prims.C 2014-04-17 17:15:36 UTC (rev 26) @@ -186,7 +186,7 @@ PRIMS[ww] = By[index]; ww++; PRIMS[ww] = Bz[index]; ww++; double u0L = u0[index]; - + double CONSERVS[NUM_CONSERVS] = {rho_star[index], mhd_st_x[index],mhd_st_y[index],mhd_st_z[index],tau[index]};
double METRIC_LAP_PSI4[NUMVARS_METRIC_AUX]; @@ -284,7 +284,6 @@ for(int ii=0;ii<3;ii++) { check = harm_primitives_gammalaw_lowlevel(index,i,j,k,x,y,z,METRIC,METRIC_PHYS,METRIC_LAP_PSI4, CONSERVS,PRIMS,u0L, g4dn,g4up, stats,eos); - if(check==0) ii=4; } } else { @@ -299,6 +298,7 @@ rho_star_fix_applied++; }
+ // Enforce limits on primitive variables and recompute conservatives. static const int already_computed_physical_metric_and_inverse=1; enforce_limits_on_primitives_and_recompute_conservs(already_computed_physical_metric_and_inverse,PRIMS,u0L,eos,METRIC,g4dn,g4up, TUPMUNU,TDNMUNU,CONSERVS); @@ -352,7 +352,7 @@ useconds = end.tv_usec - start.tv_usec;
mtime = ((seconds) * 1000 + useconds/1000.0) + 0.999; // We add 0.999 since mtime is a long int; this rounds up the result before setting the value. Here, rounding down is incorrect. - printf("Pointcount: %d Font fixes: %d VL: %d Failures: %d, InHoriz: %d / %d = %.3e\t%.2f solutions/second, Error: %e, ErrDenom: %.15e, rho*fixes: %d Lev: %d\n", + printf("Pointcount: %d Font fixes: %d VL: %d Failures: %d, InHoriz: %d / %d = %.3e\t%.0f solutions/second, Error: %e, ErrDenom: %.15e, rho*fixes: %d Lev: %d\n", pointcount,font_fixes,vel_limited_ptcount, failures, failures_inhoriz,pointcount_inhoriz,failures_inhoriz/((double)pointcount_inhoriz+1e-10),
File [modified]: reconstruct_set_of_prims_PPM.C Delta lines: +3 -3 =================================================================== --- trunk/src/reconstruct_set_of_prims_PPM.C 2014-04-16 16:00:47 UTC (rev 25) +++ trunk/src/reconstruct_set_of_prims_PPM.C 2014-04-17 17:15:36 UTC (rev 26) @@ -339,11 +339,11 @@ double dP2 = U[PRESSURE][PLUS2] - U[PRESSURE][MINUS2];
// MODIFICATION TO STANDARD PPM: - // Cure roundoff error issues when dP1==0 or dP2==0 to 7 or more significant digits. + // Cure roundoff error issues when dP1==0 or dP2==0 to 6 or more significant digits. double avg1=0.5*(U[PRESSURE][PLUS1] + U[PRESSURE][MINUS1]); double avg2=0.5*(U[PRESSURE][PLUS2] + U[PRESSURE][MINUS2]); - if(fabs(dP1)/avg1<1e-7) dP1=0.0; /* If this is triggered, there is NO shock */ - if(fabs(dP2)/avg2<1e-7) dP2=0.0; /* If this is triggered alone, there may be a shock. Otherwise if triggered with above, NO shock. */ + if(fabs(dP1)/avg1<1e-6) dP1=0.0; /* If this is triggered, there is NO shock */ + if(fabs(dP2)/avg2<1e-6) dP2=0.0; /* If this is triggered alone, there may be a shock. Otherwise if triggered with above, NO shock. */
double dP1_over_dP2=1.0; if (dP2 != 0.0) dP1_over_dP2 = dP1/dP2;