[Commits] [svn:einsteintoolkit] IllinoisGRMHD/trunk/src/ (Rev. 26)

[zachetie at gmail.com]

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;