Difference between revisions of "KPP solvers FAQ"

From Geos-chem
Jump to: navigation, search
(My simulation with KPP died. Can I get some debug output?)
(Example: Tightening tolerances for SOA simulations)
(14 intermediate revisions by 2 users not shown)
Line 1: Line 1:
Welcome to the "KPP and GEOS-Chem - Frequently Asked Questions" page.
+
== Current KPP documentation ==
  
== What is KPP? ==
+
We have migrated our KPP-for-GEOS-Chem documentation to '''https://kpp.readthedocs.io'''.
  
The Kinetic PreProcessor is a software that automatically generates F90 code that solves the chemistry (defined by input files) in a box model.
+
--[[User:Bmy|Bob Yantosca]] ([[User talk:Bmy|talk]]) 16:45, 2 March 2021 (UTC)
  
The KPP input files have been generated from the globchem.dat GEOS-Chem file with a Perl script, and KPP was run. The Kpp-generated code has been interfaced with GEOS-Chem.
+
== How do I choose the absolute and relative tolerance? ==
  
KPP has been run twice for the standard NOx-Ox-VOC GEOS-Chem simulation: once for the 43 transported tracers model, and once for the 54 tracers version (ie with in-line secondary aerosols), since the two models have a slightly different globchem.dat.
+
The tolerance are set in the <tt>DO_FLEXCHEM</tt> routine (in <tt>flexchem_mod.F90</tt>). GEOS-Chem 13.0.0 uses these default values:
  
== What are the pros of using KPP? ==
+
    !%%%%% CONVERGENCE CRITERIA %%%%%
 
+
In short, '''faster and better'''. In many points:
+
 
+
# It gives you a choice of different solvers that can have relative advantages depending on the level of desired accuracy [Sandu and Sander, 2006].
+
# One of the solver (radau5) is also recommended for "accurate reference solutions" [Sandu and Sander, 2006].
+
# By default, GEOS-Chem is compiled with the Rosenbrock Rodas3 solver:
+
## Tests reported in Henze et al. [2007] and Eller et al. [2009] show that the Rosenbrock solver is twice faster than SMVGEAR for the same accuracy. However this accuracy is with respect to a reference output obtained with the same solver.
+
## "At the very least, [using the Rosenbrock solver] results in an improvement in the numerical solution [...] for slightly less computational cost." [Henze et al., 2007]
+
## At low NOx regimes, like paleo simulations, the Rosenbrock solver is found more stable than SMVGEAR (L. Murray, personal communication, GEOS-Chem v8-01-04, 2009).
+
 
   
 
   
To get a measure of (3.2), a detailed comparison between SMVGEAR, Rosenbrock (Rodas3), and a reference solution was performed [http://acmg.seas.harvard.edu/geos/word_pdf_docs/rodas3_smvgear_comp.pdf (pdf here)]. It shows that Rodas3 can be almost 10% faster for a more accurate solution. SMVGEAR is
+
    ! Absolute tolerance
particularly bad at reproducing isoprene and the inorganic sulfur nitrates.
+
     ATOL     = 1e-2_dp
 
+
Here is the full list of transported species that show an RMS error larger than 10% for SMVGEAR and/or Rodas3. The numbers are RMS of relative error in %, from a comparison with a reference run, obtained with a Runge-Kutta order 5 integrator and very tight tolerances. All runs are for 13 days without transport on 4x5, the comparison is done with concentrations averaged over the 13th day.
+
 
+
          SMVGEAR    Rodas3      
+
        ====================
+
  ISOP 222.22 
+
  NIT   87.36     11.82
+
  MVK   33.26 
+
  MACR   30.40 
+
  PMN   23.91 
+
  NH3   19.10 
+
  DMS   11.43 
+
  ALK4   10.72 
+
 
+
See also "How does the default Rodas-3 Rosenbrock integrator compare with the SMVGEAR ?"
+
 
+
== What are the cons of using KPP? ==
+
 
+
If you make ANY modification to the <tt>globchem.dat</tt> chemistry mechanism file, including:
+
 
+
* Adding or deleting species
+
* Adding or deleting reactions
+
* Modifying rate constants
+
* Switching species from active to inactive or dead (or vice-versa)
+
 
+
Then you must do the following:
+
 
+
# Run KPP to generate the <tt>gckpp*.f90</tt> files
+
# Adapt the new <tt>gckpp*.f90</tt> files for GEOS-Chem
+
# Replace the existing <tt>gckpp*.f90</tt> files in the GEOS-Chem <tt>KPP</tt> directory structure with the new files you just created.
+
 
+
The whole process is described in [[Interfacing GEOS-Chem with KPP|"Interfacing GEOS-Chem with KPP"]].
+
 
+
Generally, if you are developing GEOS-Chem, be aware that the interface with KPP variables is done through both CSPEC and RRATE arrays (see <tt>physproc.f</tt>, <tt>calcrate.f</tt>, <tt>chemdr.f</tt>, and <tt>chemistry_mod.f</tt>).
+
 
+
Also, those of you who are using the Intel Fortran Compiler ("IFORT") version 9.1 should note that the [[GEOS-Chem_v8-02-03#KPP_is_not_compatible_with_IFORT_9.1|KPP solver package has problems with this compiler version]].  You should upgrade to IFORT 10 or higher.
+
 
+
== Can I still use the KPP-generated solver if I modify GEOS-Chem? ==
+
 
+
Yes.  See the section: [[#What are the cons of using KPP.3F|What are the cons of using KPP?]]
+
 
+
== Which GC simulation can use the KPP solvers? ==
+
 
+
For now, only the standard NOx-Ox-VOC simulations, with in- or off-line secondary aerosols. As of v8-02-03, the difference is in the number of transported tracers: 43 and 54.  ('''''NOTE:''''' We will add a third option into KPP, the [[New isoprene scheme|new Caltech isoprene chemistry scheme]] in versions [[GEOS-Chem v8-02-05|v8-02-05]] and higher.)
+
 
+
By default, the code is compiled for 43 tracers.
+
 
+
When switching between 43 and 54 tracers, you need to call:
+
 
+
make clean
+
 
+
For the 54 tracers simulation, you can use the <tt>CHEM</tt> flag when calling make (this is recommended):
+
 
+
make CHEM=SOA
+
 
+
You may also use the deprecated <tt>NTRAC</tt> option:
+
 
+
make NTRAC=54
+
 
+
== Can I still use SMVGEAR? ==
+
 
+
Yes. There is a switch in the chemistry menu of the input.geos to choose between SMVGEAR and the KPP-derived solver at runtime.
+
 
+
== Which KPP solvers can I use? ==
+
 
+
As of GEOS-Chem v8-02-03, you can choose between several stiff numerical solvers that depend on either the backward differentiation formula (BDF), Runge Kutta, or Rosenbrock methods.
+
 
+
Runge-Kutta (Fully Implicit 3-stage Runge-Kutta methods), based on several quadratures (*):
+
 
   
 
   
  Radau-2A  quadrature (order 5) [1] (default, stiffly accurate, robust, expensive)
+
    ! Relative tolerance
  Radau-1A  quadrature (order 5) [4]
+
    IF ( Input_Opt%LUCX  ) THEN
  Lobatto-3C quadrature (order 4) [2]
+
      ! UCX-based mechanisms
  Gauss     quadrature (order 6) [3]
+
      !RTOL      = 2e-2_dp
+
      !RTOL     = 1e-2_dp
LSODES (BDF, Livermore ODE solver, sparse array version)
+
      RTOL      = 0.5e-2_dp
+
    ELSE
Rosenbrock methods (*):
+
      ! Non-UCX mechanisms
+
      RTOL      = 1e-2_dp
  Ros-2    [1]
+
    ENDIF
  Ros-3    [2]
+
  Ros-4    [3]
+
  Rodas-3  [4] (default)
+
  Rodas-4  [5]
+
  
(*) For the Rosenbrock and Runge-Kutta methods, the default method can be overwritten by setting (before compilation, in chemistry_mod.f) the variable ictrl(3) to the number in square brackets.
+
For the Rodas-3 Rosenbrock solver, runtime and accuracy of the solution for various set of tolerances are [http://acmg.seas.harvard.edu/geos/word_pdf_docs/rodas3_smvgear_comp.pdf described in this pdf].  
  
For details on those methods and their implementation, see Sandu and Sander [2006] and its electronic supplement, and references therein.
+
If you switch to Livermore solver (LSODE), you do need to set ATOL between 1d4 and 1d6, to have a performance similar to SMVGEAR.
  
== How do I switch solvers? ==
+
For an accurate solution with Radau5, you need small tolerances: RTOL=1d-8 for example.
  
Solvers are chosen at three steps: ''before compilation'', at ''compilation'', and at ''runtime''.
+
=== Solver convergence issues posted on Github  ===
  
'''At runtime''', you can switch between SMVGEAR and the KPP-derived solver with a single flag in the input.geos file. See the GEOS-Chem manual.
+
See these issues:
  
The KKP solver is chosen '''at compilation'''. When calling make, you can overwrite the default rosenbrock solver, with the KPPSOLVER flag:
+
#[https://github.com/geoschem/geos-chem/issues/66 geoschem/geos-chem #66]
 
+
make KPPSOLVER=<solver>
+
+
where <solver> must be one of: '''radau5''', '''lsodes''', '''runge_kutta''', '''rosenbrock'''
+
+
(radau5 is a standalone Runge-Kutta solver of order 5 with Radau-2A quadrature)
+
 
+
You can further fine tune your choice of integrator by setting ictrl(3) in chemistry_mod.f '''before compilation'''. See "Which KPP solvers can I use?" above.
+
 
+
== How do I choose the absolute and relative tolerance? ==
+
 
+
The tolerances used by solvers to check convergence are set differently for SMVGEAR and the KPP solvers.
+
 
+
For SMVGEAR, the range of absolute tolerance (ATOL) is defined in the mglob.dat file by YLOWU and YHIU. The relative tolerance (RTOL) is also set in mglob.dat with ERRMAXU, but it is overwritten in the code if you set it above 1d-3 (this is its max value, in other words).
+
 
+
For KPP-derived solvers, the tolerance are set in the <tt>GCKPP_DRIVER</tt> routine (in <tt>chemistry_mod.f</tt>). The code is distributed with the suggested values of:
+
 
+
*RTOL = 2d-1
+
*ATOL = 1d-2
+
 
+
which give a solution slightly more accurate but ~8% faster than SMVGEAR.
+
 
+
For the Rodas-3 Rosenbrock solver, runtime and accuracy of the solution for various set of tolerances are [http://acmg.seas.harvard.edu/geos/word_pdf_docs/rodas3_smvgear_comp.pdf described in this pdf]. Some results are reproduced in the section [[#How_does_the_default_Rodas-3_Rosenbrock_integrator_compare_with_the_SMVGEAR_.3F|How does the default Rodas-3 Rosenbrock integrator compare with the SMVGEAR?]] below.
+
 
+
If you switch to Livermore solver (LSODE), you do need to set ATOL between 1d4 and 1d6, to have a performance similar to SMVGEAR.
+
 
+
For an accurate solution with Radau5, you need small tolerances: RTOL=1d-8 for example.
+
  
 
=== Example: Tightening tolerances for SOA simulations ===
 
=== Example: Tightening tolerances for SOA simulations ===
Line 149: Line 41:
 
Here is an instance in which the relative error tolerance <tt>RTOL</tt> had to be tightened when performing a GEOS-Chem simulation with the secondary organic aerosol tracers:
 
Here is an instance in which the relative error tolerance <tt>RTOL</tt> had to be tightened when performing a GEOS-Chem simulation with the secondary organic aerosol tracers:
  
'''''[mailto:tai@seas.harvard.edu Amos P.K. Tai] wrote:'''''
+
'''''Amos P.K. Tai] wrote:'''''
  
:I've performed a few 3-month 4x5 GEOS-5 runs with the new v8-02-03 KPP solver with two different relative tolerance levels (<tt>RTOL</tt>) in the <tt>GCKPP_DRIVER</tt> subroutine in the <tt>chemistry_mod.f</tt> module.  Then I counted the number of times I got the following error:
+
<blockquote>I've performed a few 3-month 4x5 GEOS-5 runs with the v8-02-03 KPP solver with two different relative tolerance levels (<tt>RTOL</tt>) in the <tt>GCKPP_DRIVER</tt> subroutine in the <tt>chemistry_mod.f</tt> module.  Then I counted the number of times I got the following error:</blockquote>
  
 
     Forced exit from Rosenbrock due to the following error:
 
     Forced exit from Rosenbrock due to the following error:
Line 157: Line 49:
 
     T=  2690.36541862269      and H=  6.025073049527794E-013
 
     T=  2690.36541862269      and H=  6.025073049527794E-013
  
:When this error occurs once in a row, the run would continue; but if it occurs twice in a row, the run crashes.  I counted the number of this error in each of my 3-month runs that went without crashing.
+
<blockquote>When this error occurs once in a row, the run would continue; but if it occurs twice in a row, the run crashes.  I counted the number of this error in each of my 3-month runs that went without crashing.
  
:*For RTOL = 2d-1 (recommended), I've done four 3-month runs; the number of error ranges from 92 to 114.
+
*For RTOL = 2d-1 (recommended), I've done four 3-month runs; the number of error ranges from 92 to 114.
:*For RTOL = 5d-2 (reduced), I've done two 3-month runs (otherwise all the same as before); no errors were detected in both runs.
+
*For RTOL = 5d-2 (reduced), I've done two 3-month runs (otherwise all the same as before); no errors were detected in both runs.
  
:So, reducing the tolerance level, as you have suggested before, appears to solve the problem of KPP solver leading to a crash.
+
So, reducing the tolerance level, as you have suggested before, appears to solve the problem of KPP solver leading to a crash.</blockquote>
  
 
--[[User:Bmy|Bob Y.]] 10:22, 10 November 2009 (EST)
 
--[[User:Bmy|Bob Y.]] 10:22, 10 November 2009 (EST)
Line 170: Line 62:
 
'''''[mailto:Tao.Zeng@dnr.state.ga.us Tao Zeng] wrote:'''''
 
'''''[mailto:Tao.Zeng@dnr.state.ga.us Tao Zeng] wrote:'''''
  
:By reducing RTOL from 0.2 to 0.05, you were able to run the model for longer time without errors.
+
<blockquote>
 +
By reducing RTOL from 0.2 to 0.05, you were able to run the model for longer time without errors.
  
:I tried your method and found it works. But I am quite confused with the setting. To my understanding, we enlarge the tolerance to make it not so strict, it's like to open the door a little wider. But why does it work in the opposite way in the KPP case?
+
I tried your method and found it works. But I am quite confused with the setting. To my understanding, we enlarge the tolerance to make it not so strict, it's like to open the door a little wider. But why does it work in the opposite way in the KPP case?
 +
</blockquote>
  
'''''[mailto:ccarouge@seas.harvard.edu Claire Carouge] wrote:'''''
+
'''''Claire Carouge wrote:'''''
 +
<blockquote>I am not sure of what is happening but I think I have an idea.
  
:I am not sure of what is happening but I think I have an idea.  
+
The errors occur because the system arrives to a state too far from the truth to be able to converge. By tightening the tolerances, you make sure the system stays closer to the truth at every time step. Then, the problematic time steps start the chemistry with a system closer to the true state, enabling the chemistry to converge even with a tighter tolerance.
  
:The errors occur because the system arrives to a state too far from the truth to be able to converge. By tightening the tolerances, you make sure the system stays closer to the truth at every time step. Then, the problematic time steps start the chemistry with a system closer to the true state, enabling the chemistry to converge even with a tighter tolerance.
+
Typically, if the first time step of chemistry couldn't converge, tightening the tolerances wouldn't work but loosening the tolerance would. Since here it is not the first time step, tightening the tolerances can work by keeping the system errors smaller before.</blockquote>
 
+
:Typically, if the first time step of chemistry couldn't converge, tightening the tolerances wouldn't work but loosening the tolerance would. Since here it is not the first time step, tightening the tolerances can work by keeping the system errors smaller before.
+
  
 
--[[User:Bmy|Bob Y.]] 10:35, 11 March 2010 (EST)
 
--[[User:Bmy|Bob Y.]] 10:35, 11 March 2010 (EST)
 
== How does the default Rodas-3 Rosenbrock integrator compare with the SMVGEAR? ==
 
 
Faster by up to 10% for a slightly more accurate solution:
 
*at one end, there is a RTOL between 0.2-0.3 that gives a solution as accurate but ~10% faster
 
*at the other end, there is a RTOL between 0.03-0.04 that gives a solution as fast but an order of magnitude more accurate (~one significant digit better)
 
See [http://acmg.seas.harvard.edu/geos/word_pdf_docs/rodas3_smvgear_comp.pdf this pdf for details]. Here we summarize some of the findings:
 
 
Absolute tolerance ATOL has little to no effect in the [1e-2,1e+2] examined range. Using
 
ATOL=1d-2, here is a comparison between SMVGEAR and Rodas-3 with several relative tolerances
 
RTOL. The RMS of relative error in % of species with an error
 
larger than 10%, and the runtime with respect to SMVGEAR are:
 
 
                         
 
        Rodas3          SMVGEAR          Rodas3          Rodas3
 
        RTol=0.3                          RTol=0.2        RTol=0.1
 
        =======================================================================
 
  ISOP         222.22 
 
  NIT 14.75            87.36            16.33          11.82
 
  MVK   33.26 
 
  MACR   30.40
 
  PMN   23.91 
 
  NH3   19.10
 
  DMS   11.43 
 
  ALK4   10.72
 
  H2O2  2.23e9 
 
  HNO4 26.12
 
  PAN 23.43
 
  PPN 11.47
 
 
runtime  87%              100%            92%            93%
 
  notes                                      recommended    used in the adjoint-code
 
 
== Can I run KPP and generate a box model myself ? ==
 
Yes. Installing, and running KPP, as well as interfacing the generated code with GEOS-Chem is described in [[Interfacing GEOS-Chem with KPP|"Interfacing GEOS-Chem with KPP"]].
 
 
== My simulation with KPP died.  Can I get some debug output? ==
 
 
Yes.  The trick is to add some debug print statements into the Integrator module file.
 
 
Let's assume we are using the Rosenbrock solver option (i.e. the most popular option).  The Rosenbrock integrator file is stored in the GEOS-Chem subdirectory:
 
 
KPP/int/gckpp_Integrator_rosenbrock.f90
 
 
Open this file with a text editor and look for a subroutine named <tt>ros_ErrorMsg</tt>.  Then add the following WRITE statements into this subroutine, such th
 
 
  SUBROUTINE ros_ErrorMsg(Code,T,H,IERR)
 
 
!    Handles all error messages
 
 
 
 
 
    REAL(kind=dp), INTENT(IN) :: T, H
 
    INTEGER, INTENT(IN)  :: Code
 
    INTEGER, INTENT(OUT) :: IERR 
 
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
!%%% BMY DEBUG FOR KPP ERROR PRINTOUT (bmy, 2/28/11)
 
!%%%
 
    INTEGER :: N, IUNIT
 
    IUNIT = 65
 
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 
    IERR = Code
 
    PRINT * , &
 
      'Forced exit from Rosenbrock due to the following error:'
 
     
 
    SELECT CASE (Code)
 
    CASE (-1)   
 
      PRINT * , '--> Improper value for maximal no of steps'
 
    CASE (-2)   
 
      PRINT * , '--> Selected Rosenbrock method not implemented'
 
    CASE (-3)   
 
      PRINT * , '--> Hmin/Hmax/Hstart must be positive'
 
    CASE (-4)   
 
      PRINT * , '--> FacMin/FacMax/FacRej must be positive'
 
    CASE (-5)
 
      PRINT * , '--> Improper tolerance values'
 
    CASE (-6)
 
      PRINT * , '--> No of steps exceeds maximum bound'
 
    CASE (-7)
 
      PRINT * , '--> Step size too small: T + 10*H = T', &
 
            ' or H < Roundoff'
 
    CASE (-8)   
 
      PRINT * , '--> Matrix is repeatedly singular'
 
    CASE DEFAULT
 
      PRINT *, 'Unknown Error code: ', Code
 
    END SELECT
 
   
 
    PRINT *, "T=", T, "and H=", H
 
     
 
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
!%%% BMY DEBUG FOR KPP ERROR PRINTOUT (bmy, 2/28/11)
 
!%%%
 
 
!$OMP CRITICAL
 
 
    ! Print KPP settings
 
    OPEN( IUNIT, FILE='kpp_debug.txt', FORM='formatted' )
 
    WRITE( IUNIT, '(a8,1x,es15.8)' ) 'TIME    ', TIME
 
    WRITE( IUNIT, '(a8,1x,es15.8)' ) 'SUN    ', SUN
 
    WRITE( IUNIT, '(a8,1x,es15.8)' ) 'TEMP    ', TEMP
 
    WRITE( IUNIT, '(a8,1x,es15.8)' ) 'RTOLS  ', RTOLS
 
    WRITE( IUNIT, '(a8,1x,es15.8)' ) 'TSTART  ', TSTART
 
    WRITE( IUNIT, '(a8,1x,es15.8)' ) 'TEND    ', TEND
 
    WRITE( IUNIT, '(a8,1x,es15.8)' ) 'DT      ', DT
 
    WRITE( IUNIT, '(a8,1x,es15.8)' ) 'STEPMIN ', STEPMIN
 
    WRITE( IUNIT, '(a8,1x,es15.8)' ) 'STEPMAX ', STEPMAX
 
    WRITE( IUNIT, '(a8,1x,es15.8)' ) 'CFACTOR ', CFACTOR
 
    WRITE( IUNIT, '(a8,1x,es15.8)' ) 'DDMTYPE ', DDMTYPE
 
    WRITE( IUNIT, '(a8,1x,es15.8)' ) 'NTT    ', NTT
 
    WRITE( IUNIT, '(a8,1x,es15.8)' ) 'JLOOP  ', JLOOP
 
    call flush(IUNIT)
 
 
    WRITE( IUNIT, '(/,"Absolute tolerance ATOL")' )
 
    DO N = 1, NVAR
 
      WRITE( IUNIT, '(i4,1x,es15.8)' ) N, ATOL(N)
 
    ENDDO
 
    call flush(IUNIT)
 
 
    WRITE( IUNIT, '(/,"Relative tolernaRTOL")' )
 
    DO N = 1, NVAR
 
      WRITE( IUNIT, '(i4,1x,es15.8)' ) N, RTOL(N)
 
    ENDDO
 
    call flush(IUNIT)
 
 
    WRITE( IUNIT, '(/,"Varying concentrations VAR")' )
 
    DO N = 1, NVAR
 
      WRITE( IUNIT, '(i4,1x,es15.8)' ) N, VAR(N)
 
    ENDDO
 
    call flush(IUNIT)
 
 
    WRITE( IUNIT, '(/,"Fixed concentrations FIX")' )
 
    DO N = 1, NFIX
 
      WRITE( IUNIT, '(i4,1x,es15.8)' ) N, FIX(N)
 
    ENDDO
 
    call flush(IUNIT)
 
 
    WRITE( IUNIT, '(/,"Rate constants RCONST")' )
 
    DO N = 1, NREACT
 
      WRITE( IUNIT, '(i4,1x,es15.8)' ) N, RCONST(N)
 
    ENDDO
 
    call flush(IUNIT)
 
 
    WRITE( IUNIT, '(/,"VAR_ADJ")' )
 
    DO N = 1, NVAR
 
      WRITE( IUNIT, '(i4,1x,es15.8)' ) N, V_CSPEC(N)
 
    ENDDO
 
    call flush(IUNIT)
 
 
    WRITE( IUNIT, '(/,"VAR_R_ADJ")' )
 
    DO N = 1, NCOEFF
 
      WRITE( IUNIT, '(i4,1x,es15.8)' ) N, VAR_R_ADJ(N)
 
    ENDDO
 
    call flush(IUNIT)
 
 
    WRITE( IUNIT, '(/,"V_CSPEC")' )
 
    DO N = 1, NVAR
 
      WRITE( IUNIT, '(i4,1x,es15.8)' ) N, V_CSPEC(N)
 
    ENDDO
 
    call flush(IUNIT)
 
 
    WRITE( IUNIT, '(/,"V_CSPEC_ADJ")' )
 
    DO N = 1, NVAR
 
      WRITE( IUNIT, '(i4,1x,es15.8)' ) N, V_CSPEC_ADJ(N)
 
    ENDDO
 
 
    CLOSE( IUNIT )
 
 
!$OMP END CRITICAL
 
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 
  END SUBROUTINE ros_ErrorMsg
 
 
This will write the values of several important KPP variables to a file named <tt>kpp_debug.txt</tt>.  This may give you some indication of what is causing your simulation to die (i.e. a concentration is going negative, NaN, div-by-zero, etc).
 
 
''Thanks to Adrian Sandu and Claire Carouge for their suggestions.''
 
 
--[[User:Bmy|Bob Y.]] 17:03, 28 February 2011 (EST)
 
 
== Known issues ==
 
 
=== KPP is not compatible with IFORT 9.1 ===
 
 
Please see [[GEOS-Chem v8-02-03#KPP is not compatible with IFORT 9.1|this wiki post]] about compilation problems KPP with IFORT 9.1.
 
 
=== Prod/loss diagnostic does not work with KPP ===
 
 
'''''[mailto:michael.barkley@ed.ac.uk Mike Barkley] wrote:'''''
 
 
:I want to use the prod/loss diagnostic to look at HCHO. However, I did a quick basic run to look at 4 families: POX & LOX, and PCH2O & LCH2O, with the KPP solver switched on & off. When I use smvgear enabled the ND65 output looks normal, but when I use the KPP  solver only output for POX is produced, all the other families (LOX, PCH2O & LCH2O) are zero arrays. Is the ND65 diagnostic only compatible with SMVGEAR?
 
 
'''''[mailto:yantosca@seas.harvard.edu Bob Yantosca] wrote:'''''
 
 
:Yes, at present ND65 is only compatible with SMVGEAR.  KPP is 3rd-party code (which is generated automatically), and there are no facilities at present to compute prod/loss from it.
 
 
:We are working on a fix for this, possibly for GEOS-Chem v8-03-02 or later.
 
 
--[[User:Bmy|Bob Y.]] 16:15, 28 April 2010 (EDT)
 
 
== References ==
 
#Eller et al., Geosci. Model Dev., 2, 89-96, 2009
 
#Henze et al., ACP, 7 ,2413-2433, 2007
 
#Sandu and Sander, ACP, 6, 187-195, 2006
 
 
--[[User:Plesager|phs]] 16:54, 17 September 2009 (EDT)
 

Revision as of 16:54, 2 March 2021

Current KPP documentation

We have migrated our KPP-for-GEOS-Chem documentation to https://kpp.readthedocs.io.

--Bob Yantosca (talk) 16:45, 2 March 2021 (UTC)

How do I choose the absolute and relative tolerance?

The tolerance are set in the DO_FLEXCHEM routine (in flexchem_mod.F90). GEOS-Chem 13.0.0 uses these default values:

   !%%%%% CONVERGENCE CRITERIA %%%%%

   ! Absolute tolerance
   ATOL      = 1e-2_dp

   ! Relative tolerance
   IF ( Input_Opt%LUCX  ) THEN
      ! UCX-based mechanisms
      !RTOL      = 2e-2_dp
      !RTOL      = 1e-2_dp
      RTOL      = 0.5e-2_dp
   ELSE
      ! Non-UCX mechanisms
      RTOL      = 1e-2_dp
   ENDIF

For the Rodas-3 Rosenbrock solver, runtime and accuracy of the solution for various set of tolerances are described in this pdf.

If you switch to Livermore solver (LSODE), you do need to set ATOL between 1d4 and 1d6, to have a performance similar to SMVGEAR.

For an accurate solution with Radau5, you need small tolerances: RTOL=1d-8 for example.

Solver convergence issues posted on Github

See these issues:

  1. geoschem/geos-chem #66

Example: Tightening tolerances for SOA simulations

Here is an instance in which the relative error tolerance RTOL had to be tightened when performing a GEOS-Chem simulation with the secondary organic aerosol tracers:

Amos P.K. Tai] wrote:

I've performed a few 3-month 4x5 GEOS-5 runs with the v8-02-03 KPP solver with two different relative tolerance levels (RTOL) in the GCKPP_DRIVER subroutine in the chemistry_mod.f module. Then I counted the number of times I got the following error:
   Forced exit from Rosenbrock due to the following error:
   --> Step size too small: T + 10*H = T or H < Roundoff
   T=   2690.36541862269      and H=  6.025073049527794E-013
When this error occurs once in a row, the run would continue; but if it occurs twice in a row, the run crashes. I counted the number of this error in each of my 3-month runs that went without crashing.
  • For RTOL = 2d-1 (recommended), I've done four 3-month runs; the number of error ranges from 92 to 114.
  • For RTOL = 5d-2 (reduced), I've done two 3-month runs (otherwise all the same as before); no errors were detected in both runs.
So, reducing the tolerance level, as you have suggested before, appears to solve the problem of KPP solver leading to a crash.

--Bob Y. 10:22, 10 November 2009 (EST)

More about tolerances

Tao Zeng wrote:

By reducing RTOL from 0.2 to 0.05, you were able to run the model for longer time without errors.

I tried your method and found it works. But I am quite confused with the setting. To my understanding, we enlarge the tolerance to make it not so strict, it's like to open the door a little wider. But why does it work in the opposite way in the KPP case?

Claire Carouge wrote:

I am not sure of what is happening but I think I have an idea.

The errors occur because the system arrives to a state too far from the truth to be able to converge. By tightening the tolerances, you make sure the system stays closer to the truth at every time step. Then, the problematic time steps start the chemistry with a system closer to the true state, enabling the chemistry to converge even with a tighter tolerance.

Typically, if the first time step of chemistry couldn't converge, tightening the tolerances wouldn't work but loosening the tolerance would. Since here it is not the first time step, tightening the tolerances can work by keeping the system errors smaller before.

--Bob Y. 10:35, 11 March 2010 (EST)