Running KPP with an example stratospheric mechanism

Here we consider as an example a very simple Chapman-like mechanism for stratospheric chemistry:

\[\begin{split}\begin{aligned} O_2 + h\nu & \longrightarrow & 2 O & ~~~~~~~~~~ (1)\\ O + O_2 & \longrightarrow & O_3 & ~~~~~~~~~~ (2)\\ O_3 + h\nu & \longrightarrow & O + O_2 & ~~~~~~~~~~ (3)\\ O + O_3 & \longrightarrow & 2 O_2 & ~~~~~~~~~~ (4)\\ O_3 + h\nu & \longrightarrow & O(^1D) + O_2 & ~~~~~~~~~~ (5)\\ O(^1D) + M & \longrightarrow & O + M & ~~~~~~~~~~ (6)\\ O(^1D) + O_3 & \longrightarrow & 2 O_2 & ~~~~~~~~~~ (7)\\ NO + O_3 & \longrightarrow & NO_2 + O_2 & ~~~~~~~~~~ (8)\\ NO_2 + O & \longrightarrow & NO + O_2 & ~~~~~~~~~~ (9)\\ NO_2 + h\nu & \longrightarrow & NO + O & ~~~~~~~~~~ (10) \end{aligned}\end{split}\]

We use the mechanism with the purpose of illustrating the KPP capabilities. However, the software tools are general and can be applied to virtually any kinetic mechanism.

We focus on Fortran90. Particularities of the C and Matlab languages are discussed in the #LANGUAGE section.

Important

Most of the recent KPP developments described in this manual have been added into the Fortran90 language. We look to members of the KPP user community to spearhead development in C, Matlab, and other languages.

The KPP input files (with suffix .kpp) specify the target model, the target language, the integrator the driver program. etc. The file name (without the .kpp) serves as the root name for the simulation. Here we will refer to this name as ROOT. Since the root name will be incorporated into Fortran90 module names, it can only contain valid Fortran90 characters, i.e. letters, numbers, and the underscore.

The sections below outline the steps necessary to build an run a “box-model” simulation with an example mechanism.

1. Create a folder for the example

Create a folder in which to build and run the example mechanism:

$ cd $HOME
$ mkdir small_strato_example
$ cd small_strato_example

In the following sections we will refer to $HOME/small_strato_example as “the example folder”.

2. Create a KPP Definition File

Create a KPP definition file in the example folder. The name of this file will always be ROOT.kpp, where ROOT is the name of the chemical mechanism.

For this example, write the following lines into a file named small_strato.kpp in the example folder:

#MODEL      small_strato
#LANGUAGE   Fortran90
#DOUBLE     ON
#INTEGRATOR rosenbrock
#DRIVER     general
#JACOBIAN   SPARSE_LU_ROW
#HESSIAN    ON
#STOICMAT   ON

Important

KPP will look for the relevant files (e.g. mechanism definition, driver, etc.) in the proper subfolders of KPP_HOME. Therefore you won’t need to copy these manually to the example folder.

Also note, KPP command options can be either uppercase or lowercase (i.e. INTEGRATOR ON or INTEGRATOR on are identical).

We will now look at the following KPP commands in small_strato.kpp.

#MODEL small_strato

The #MODEL command selects a specific kinetic mechanism (in this example, small_strato). KPP will look in the path $KPP_HOME/models/ for the model definition file small_strato.def. This file contains the following code in the KPP language.

#include small_strato.spc       { Includes file w/ species definitons     }
#include small_strato.eqn       { Includes file w/ chemical equations     }

#LOOKATALL                      { Output all species to small_strato.dat}
#MONITOR O3;N;O2;O;NO;O1D;NO2;  { Print selected species to screen        }

#CHECK O; N;                    { Check Mass Balance of oxygen & nitrogen }

#INITVALUES                     { Set initial values of species           }
  CFACTOR = 1.    ;             { and set units conversion factor to 1    }
  O1D = 9.906E+01 ;
  O   = 6.624E+08 ;
  O3  = 5.326E+11 ;
  O2  = 1.697E+16 ;
  NO  = 8.725E+08 ;
  NO2 = 2.240E+08 ;
  M   = 8.120E+16 ;

{ Fortran code to be inlined into ROOT_Global }
#INLINE F90_INIT
  TSTART = (12*3600)
  TEND = TSTART + (3*24*3600)
  DT = 0.25*3600
  TEMP = 270
#ENDINLINE

{ Matlab code to be inlined into ROOT_Global }
#INLINE MATLAB_INIT
  global TSTART TEND DT TEMP
  TSTART = (12*3600);
  TEND = TSTART + (3*24*3600);
  DT = 0.25*3600;
  TEMP = 270;
#ENDINLINE

{ C code to be inlined into ROOT_GLOBAL }
#INLINE C_INIT
  TSTART = (12*3600);
  TEND = TSTART + (3*24*3600);
  DT = 0.25*3600;
  TEMP = 270;
#ENDINLINE

File (small_strato.def) #INCLUDE-s the species file and the equation file. It also specifies parameters for running a “box-model” simualation, such as species initial values, start time, stop, time, and timestep (cf. Inlined Code).

The species file (small_strato.spc) file lists all the species in the model. Some of them are variable, meaning that their concentrations change according to the law of mass action kinetics. Others are fixed, with the concentrations determined by physical and not chemical factors (cf. #DEFVAR and #DEFFIX). For each species its atomic composition is given (unless the user chooses to ignore it).

#INCLUDE atoms.kpp
#DEFVAR
  O   = O;
  O1D = O;
  O3  = O + O + O;
  NO  = N + O;
  NO2 = N + O + O;
#DEFFIX
  M   = IGNORE;
  O2  = O + O;

The species file also includes the atoms file (atoms.kpp), which lists the periodic table of elements in an ATOM section (cf. #ATOMS).

The equation file (small_strato.eqn) contains the description of the equations in an #EQUATIONS section. The chemical kinetic mechanism is specified in the KPP language. Each reaction is described as “the sum of reactants equals the sum of products” and is followed by its rate coefficient. SUN is the normalized sunlight intensity, equal to one at noon and zero at night. Equation tags, e.g. <R1>, are optional.

#EQUATIONS { Small Stratospheric Mechanism }


<R1>  O2   + hv = 2O            : (2.643E-10) * SUN*SUN*SUN;
<R2>  O    + O2 = O3            : (8.018E-17);
<R3>  O3   + hv = O   + O2      : (6.120E-04) * SUN;
<R4>  O    + O3 = 2O2           : (1.576E-15);
<R5>  O3   + hv = O1D + O2      : (1.070E-03) * SUN*SUN;
<R6>  O1D  + M  = O   + M       : (7.110E-11);
<R7>  O1D  + O3 = 2O2           : (1.200E-10);
<R8>  NO   + O3 = NO2 + O2      : (6.062E-15);
<R9>  NO2  + O  = NO  + O2      : (1.069E-11);
<R10> NO2  + hv = NO  + O       : (1.289E-02) * SUN;

#LANGUAGE Fortran90

The #LANGUAGE command selects the language for the KPP-generated solver code. In this example we are using Fortran90.

#DOUBLE ON

The data type of the generated model can be switched between single/double precision with the #DOUBLE command. We recommend using double-precision in order to avoid integrator errors caused by roundoff or underflow/overflow.

#INTEGRATOR rosenbrock

The #INTEGRATOR command selects a numerical integration routine from the templates provided in the $KPP_HOME/int folder, or implemented by the user.

In this example, the Rosenbrock integrator and the Fortran90 language have been been specified. Therefore it will use the file $KPP_HOME/int/rosenbrock.f90.

#DRIVER general

The #DRIVER command selects a specific main program (located in the $KPP_HOME/drv folder):

  1. general_adj.f90 : Used with integrators that use the discrete adjoint method

  2. general_tlm.f90 : Used with integrators that use the tangent-linear method

  3. general.f90 : Used with all other integrators.

In this example, the rosenbrock.f90 integrator does not use either adjoint or tangent-linear methods, so the $KPP_HOME/drv/general.f90 will be used.

Other options

The other options listed control internal aspects of the integration (cf. ROOT_Jacobian and ROOT_JacobianSP), as well as activating optional outputs (cf. ROOT_Hessian and ROOT_HessianSP and ROOT_Stoichiom and ROOT_StoichiomSP).

3. Build the mechanism with KPP

Now that all the necessary files have been copied to the example folder, the small_strato mechanism can be built.

Type:

$ kpp small_strato.kpp

You should see output similar to:

This is KPP-2.6.0.

KPP is parsing the equation file.
KPP is computing Jacobian sparsity structure.
KPP is starting the code generation.
KPP is initializing the code generation.
KPP is generating the monitor data:
    - small_strato_Monitor
KPP is generating the utility data:
    - small_strato_Util
KPP is generating the global declarations:
    - small_strato_Main
KPP is generating the ODE function:
    - small_strato_Function
KPP is generating the ODE Jacobian:
    - small_strato_Jacobian
    - small_strato_JacobianSP
KPP is generating the linear algebra routines:
    - small_strato_LinearAlgebra
KPP is generating the Hessian:
    - small_strato_Hessian
    - small_strato_HessianSP
KPP is generating the utility functions:
    - small_strato_Util
KPP is generating the rate laws:
    - small_strato_Rates
KPP is generating the parameters:
    - small_strato_Parameters
KPP is generating the global data:
    - small_strato_Global
KPP is generating the stoichiometric description files:
    - small_strato_Stoichiom
    - small_strato_StoichiomSP
KPP is generating the driver from general.f90:
    - small_strato_Main
KPP is starting the code post-processing.

KPP has succesfully created the model "small_strato".

This will generate the Fortran90 code needed to solve the small_strato mechanism. Get a file listing:

ls

and you should see output similar to:

atoms.kpp                     small_strato.kpp
general.f90                   small_strato_LinearAlgebra.f90
Makefile_small_strato         small_strato_Main.f90
rosenbrock.def                small_strato_mex_Fun.f90
rosenbrock.f90                small_strato_mex_Hessian.f90
small_strato.def              small_strato_mex_Jac_SP.f90
small_strato.eqn              small_strato_Model.f90
small_strato_Function.f90     small_strato_Monitor.f90
small_strato_Global.f90       small_strato_Parameters.f90
small_strato_Hessian.f90      small_strato_Precision.f90
small_strato_HessianSP.f90    small_strato_Rates.f90
small_strato_Initialize.f90   small_strato.spc@
small_strato_Integrator.f90   small_strato_Stoichiom.f90
small_strato_Jacobian.f90     small_strato_StoichiomSP.f90
small_strato_JacobianSP.f90   small_strato_Util.f90

KPP creates Fortran90 beginning with the mechanism name (which is small_strato_ in this example). KPP also generates a human-readable summary of the mechanism (small_strato.map) as well as the Makefile_small_strato) that can be used to build the executable.

4. Build and run the small_strato example

To compile the Fortran90 code generated by KPP into an executable, type:

$ make -f Makefile_small_strato

You will see output similar to this:

gfortran -cpp -O -g  -c small_strato_Precision
gfortran -cpp -O -g  -c small_strato_Precision.f90
gfortran -cpp -O -g  -c small_strato_Parameters.f90
gfortran -cpp -O -g  -c small_strato_Global.f90
gfortran -cpp -O -g  -c small_strato_Function.f90
gfortran -cpp -O -g  -c small_strato_JacobianSP.f90
gfortran -cpp -O -g  -c small_strato_Jacobian.f90
gfortran -cpp -O -g  -c small_strato_HessianSP.f90
gfortran -cpp -O -g  -c small_strato_Hessian.f90
gfortran -cpp -O -g  -c small_strato_StoichiomSP.f90
gfortran -cpp -O -g  -c small_strato_Stoichiom.f90
gfortran -cpp -O -g  -c small_strato_Rates.f90
gfortran -cpp -O -g  -c small_strato_Monitor.f90
gfortran -cpp -O -g  -c small_strato_Util.f90
gfortran -cpp -O -g  -c small_strato_LinearAlgebra.f90
gfortran -cpp -O -g  -c small_strato_Initialize.f90
gfortran -cpp -O -g  -c small_strato_Integrator.f90
gfortran -cpp -O -g  -c small_strato_Model.f90
gfortran -cpp -O -g  -c small_strato_Main.f90
gfortran -cpp -O -g  small_strato_Precision.o    small_strato_Parameters.o    small_strato_Global.o small_strato_Function.o small_strato_JacobianSP.o small_strato_Jacobian.o small_strato_HessianSP.o small_strato_Hessian.o small_strato_Stoichiom.o small_strato_StoichiomSP.o small_strato_Rates.o   small_strato_Util.o   small_strato_Monitor.o small_strato_LinearAlgebra.o small_strato_Main.o          small_strato_Initialize.o small_strato_Integrator.o    small_strato_Model.o  -o small_strato.exe

Once compilation has finished, you can run the small_strato example by typing:

$ ./small_strato.exe | tee small_strato.log

This will run a “box-model” simulation forward several steps in time. You will see the concentrations of selected species at several timesteps displayed to the screen (aka the Unix stdout stream) as well as to a log file (small_strato.log).

If your simulation results exits abruptly with the Killed error, you will need to increase your stack memory limit. On most Linux systems the default stacksize limit is 8 kIb = or 8192 kB. You can max this out with the following commands:

If you are using bash, type:

$ ulimit -s unlimited

If you are using csh, type:

$ limit stacksize unlimited

5. Cleanup

If you wish to remove the executable (small_strato.exe), as well as the object (*.o) and module (*.mod) files generated by the Fortran compiler, type:

$ make clean

If you also wish to remove all the files that were generated by KPP (i.e. small_strato.map and small_strato_*.f90), type:

$ make distclean