Merge pull request #124 from brittonsmith/issue112

Convert iteration limit to a runtime parameter

doc/source/Parameters.rst

@@ -500,6 +500,33 @@ For all on/off integer flags, 0 is off and 1 is on.
    On/off flag to toggle calculation of rate coefficients corresponding to bremsstrahlung cooling 
    (``brem``). Default: 1
 
+.. c:var:: int max_iterations
+
+   The maximum subcycle iterations allowed when evolving the chemistry
+   network and internal energy in :c:func:`solve_chemistry`. The
+   default of 10000 should be sufficient in most situations. However,
+   certain physical conditions, such as dense gas that has been
+   photo-ionized, can lead to extremely short timesteps. In cases like
+   this, increasing the iteration limit by a factor of 10 or more is
+   enough to work through difficult conditions. For other situations
+   which are currently not well understood, this will not help. If you
+   encounter a situation where the iteration limit continues to be
+   exceeded for extremely high values, please report it. The behavior
+   of the code when the iteration limit is exceeded is controled by
+   the :c:data:`exit_after_iterations_exceeded` parameter. Default:
+   10000.
+
+.. c:var:: int exit_after_iterations_exceeded
+
+   Flag controlling the behavior when the maximum subcycle iterations
+   (set by :c:data:`max_iterations`) is exceeded in
+   :c:func:`solve_chemistry`. If set to 1, an error message will be
+   printed and the function will immediately exit with a return value
+   of 0, indicating failure. If set to 0, the message will be
+   produced and there will be no further integration, but the function
+   will proceed to a clean exit such that the simulation can be
+   continued. Default: 0.
+
 .. c:var:: int omp_nthreads
 
    Sets the number of OpenMP threads.  If not set, this will be set to

src/clib/grackle_chemistry_data.h

@@ -164,6 +164,12 @@ typedef struct
   int recombination_cooling_rates; //Recombination cooling
   int bremsstrahlung_cooling_rates; //Bremsstrahlung cooling
 
+  /* maximum number of subcycle iterations for solve_chemistry */
+  int max_iterations;
+
+  /* flag to exit if max iterations exceeded */
+  int exit_after_iterations_exceeded;
+
   /* number of OpenMP threads, if supported */
 # ifdef _OPENMP
   int omp_nthreads;

src/clib/grackle_fortran_interface.def

@@ -136,6 +136,8 @@ c     This is the fortran definition of grackle_chemistry_data
          INTEGER(C_INT) :: collisional_ionisation_rates
          INTEGER(C_INT) :: recombination_cooling_rates
          INTEGER(C_INT) :: bremsstrahlung_cooling_rates
+         INTEGER(C_INT) :: max_iterations
+         INTEGER(C_INT) :: exit_after_iterations_exceeded
 cc       INTEGER(C_INT) :: omp_nthreads // not supported in fortran
       END TYPE
 

src/clib/initialize_chemistry_data.c

@@ -357,6 +357,10 @@ void show_parameters(FILE *fp, chemistry_data *my_chemistry)
           my_chemistry->H2_custom_shielding);
   fprintf(fp, "H2_self_shielding                 = %d\n",
           my_chemistry->H2_self_shielding);
+  fprintf(fp, "max_iterations                    = %d\n",
+          my_chemistry->max_iterations);
+  fprintf(fp, "exit_after_iterations_exceeded    = %d\n",
+          my_chemistry->exit_after_iterations_exceeded);
 # ifdef _OPENMP
   fprintf(fp, "omp_nthreads                      = %d\n",
           my_chemistry->omp_nthreads);

src/clib/set_default_chemistry_parameters.c

@@ -134,6 +134,10 @@ chemistry_data _set_default_chemistry_parameters(void)
   my_chemistry.H2_self_shielding                      = 0;
   my_chemistry.H2_custom_shielding                    = 0;
 
+  /* max subcycle iterations for solve_rate_cool and whether to exit if exceeded */
+  my_chemistry.max_iterations                 = 10000;
+  my_chemistry.exit_after_iterations_exceeded = FALSE;
+
 //number of OpenMP threads
 # ifdef _OPENMP
   my_chemistry.omp_nthreads = omp_get_max_threads(); // maximum allowed number

src/clib/solve_chemistry.c

@@ -85,7 +85,8 @@ extern void FORTRAN_NAME(solve_rate_cool_g)(
         double *metHeating, int *clnew,
         int *iVheat, int *iMheat, gr_float *Vheat, gr_float *Mheat,
         int *iisrffield, gr_float* isrf_habing, 
-        int *iH2shieldcustom, gr_float* f_shield_custom);
+        int *iH2shieldcustom, gr_float* f_shield_custom,
+        int *itmax, int *exititmax);
 
 int local_solve_chemistry(chemistry_data *my_chemistry,
                           chemistry_data_storage *my_rates,
@@ -362,11 +363,17 @@ int local_solve_chemistry(chemistry_data *my_chemistry,
     my_fields->volumetric_heating_rate,
     my_fields->specific_heating_rate,
     &my_chemistry->use_isrf_field,
-    my_fields->isrf_habing, 
+    my_fields->isrf_habing,
     &my_chemistry->H2_custom_shielding,
-    my_fields->H2_custom_shielding_factor);
+    my_fields->H2_custom_shielding_factor,
+    &my_chemistry->max_iterations,
+    &my_chemistry->exit_after_iterations_exceeded);
 
-  return SUCCESS;
+  if (ierr == FAIL) {
+    fprintf(stderr, "Error in solve_rate_cool_g.\n");
+  }
+
+  return ierr;
 
 }
 

src/clib/solve_rate_cool_g.F

@@ -46,7 +46,8 @@ subroutine solve_rate_cool_g(icool, d, e, u, v, w, de,
      &                metDataSize, metCooling, metHeating, clnew,
      &                iVheat, iMheat, Vheat, Mheat,
      &                iisrffield, isrf_habing, 
-     &                iH2shieldcustom, f_shield_custom)
+     &                iH2shieldcustom, f_shield_custom,
+     &                itmax, exititmax)
 
 !
 !  SOLVE MULTI-SPECIES RATE EQUATIONS AND RADIATIVE COOLING
@@ -143,12 +144,13 @@ subroutine solve_rate_cool_g(icool, d, e, u, v, w, de,
 !    iMheat      - flag for using specific heating rate
 !    Vheat       - array of volumetric heating rates
 !    Mheat       - array of specific heating rates
+!    itmax       - maximum allowed sub-cycle iterations
+!    exititmax   - flag to exit if max iterations exceeded
 !
 !  OUTPUTS:
 !    update chemical rate densities (HI, HII, etc)
 !
 !  PARAMETERS:
-!    itmax   - maximum allowed sub-cycle iterations
 !    mh      - H mass in cgs units
 !
 !-----------------------------------------------------------------------
@@ -167,7 +169,7 @@ subroutine solve_rate_cool_g(icool, d, e, u, v, w, de,
      &        igammah, ih2optical, iciecool, ithreebody,
      &        ndratec, clnew, iVheat, iMheat, iH2shield, iradshield,
      &        iradtrans, iradcoupled, iradstep, irt_honly,
-     &        iisrffield, iH2shieldcustom
+     &        iisrffield, iH2shieldcustom, itmax, exititmax
 
       real*8  dx, dt, aye, temstart, temend, gamma,
      &        utim, uxyz, uaye, urho, utem, fh, dtoh, z_solar, 
@@ -250,9 +252,6 @@ subroutine solve_rate_cool_g(icool, d, e, u, v, w, de,
 
 !  Parameters
 
-      integer itmax
-      parameter (itmax = 10000)
-
 #ifdef GRACKLE_FLOAT_4
       R_PREC tolerance
       parameter (tolerance = 1.0e-05_RKIND)
@@ -321,7 +320,7 @@ subroutine solve_rate_cool_g(icool, d, e, u, v, w, de,
 
 !     Set error indicator
 
-      ierr = 0
+      ierr = 1
 
 !     Set flag for dust-related options
 
@@ -833,6 +832,10 @@ subroutine solve_rate_cool_g(icool, d, e, u, v, w, de,
             write(0,'((16(1pe8.1)))') (ttot(i),i=is+1,ie+1)
             write(0,'((16(1pe8.1)))') (edot(i),i=is+1,ie+1)
             write(0,'((16(l3)))') (itmask(i),i=is+1,ie+1)
+
+            if (exititmax .eq. 1) then
+               ierr = 0
+            endif
 #ifdef _OPENMP
 !$omp end critical
 #endif
@@ -856,6 +859,12 @@ subroutine solve_rate_cool_g(icool, d, e, u, v, w, de,
 !$omp end parallel do
 #endif
 
+!     If an error has been produced, return now.
+
+      if (ierr .eq. 0) then
+         return
+      endif
+
 !     Convert densities back to comoving from proper
 
       if (iexpand .eq. 1) then