Example illustrating application of MFNK to a simple biogeochemical model with two 
tracers (PO4 and DOP). The biogeochemical model was devised and implemented by Iris Kriest 
(ikriest@ifm-geomar.de) who has very generously allowed me to post it here. Also, some 
code fragments (notably the insolation routine and the O2 saturation calculation) are 
taken from the MITgcm biogeochemical model written by Stephanie Dutkiewicz and 
Mick Follows at MIT. My thanks to them for allowing me to include it in this distribution.
Email me if you want more details or look in the routine KielBiogeochem_NDOP_Matrix5_4/BGC_MODEL.F

References:

[1] Khatiwala et al, (2005). Accelerated simulation of passive tracers in ocean circulation 
    models. Ocean Modell., 9:51Ð69, 2005.
[2] Khatiwala, S. (2007). A computational framework for simulation of biogeochemical tracers 
    in the ocean. Glob. Biogeochem. Cy., 21 (doi: 10.1029/2007GB002923), 2007.
[3] Khatiwala, S. (2008). Fast spin up of ocean biogeochemical models using matrix-free 
    Newton-Krylov. Ocean Modell., 23:121Ð129, 2008.

You can download these from: http://www.ldeo.columbia.edu/~spk/

General notes:

NOTE-1: There are two main components to this example. The first is a biogeochemical 
model implemented using the PETSc framework (see Ref. [2]). The second is a set of driver 
routines written in Matlab to call the KINSOL Newton-Krylov solver and to communicate between 
that solver and the biogeochemical model. The typical situation is one in which the PETSc 
code runs on a parallel computer/cluster, and the Matlab code runs on the users desktop machine 
or (better) the 'frontend' of the cluster. 

NOTE-2: If you only want to peform a forward integration, build the model using:
make tracercalckielbiondop
and then modify the example job script 'runpbs_forwardrun_rzc'.
(see step 7 below).

Specific steps:

(1) Download and install PETSc (http://www.mcs.anl.gov/petsc/petsc-as/). Make sure to install 
the optional MUMPS and SuperLU/SuperLU_DIST packages. These are sparse direct solvers. We will 
use these for the implicit preconditioning. To enable these packages, you can use the following 
*additional* flags when configuring PETSc:

--download-superlu=1 --download-superlu_dist=1 --download-blacs=1  \
--download-scalapack=1 --download-mumps=1 

(2) Download and install the SUNDIALS package (https://computation.llnl.gov/casc/sundials/main.html). 
Also install the sundialsTB Matlab interface. 

NOTE-1: Install SUNDIALS and the sundialsTB on the machine on which you will be running 
Matlab. 

NOTE-2: I have made several modifications to the SUNDIALS code that include several  
bug fixes and the facility to write out a lot of useful diagnostic information. If you 
want my private version, email me.

(3) Download the following Matlab toolboxes from http://www.ldeo.columbia.edu/~spk/Code/Matlab/ 
and add them to your Matlab path.

Misc
PETSC
ClusterTools
MITgcm
Matrix

Compile the mex file 'writePetscBinaryMat_mex.c' in PETSC/ with the command:
mex -largeArrayDims writePetscBinaryMat_mex.c

(4) Follow the download instructions given here: http://www.ldeo.columbia.edu/~spk/Research/TMM/tmm.html
Once you uncompress and untar the file MIT_Matrix_Global_2.8deg.tar.gz, you should have the directory 
MIT_Matrix_Global_2.8deg, which we will call the "TOPLEVEL". Within it you will find the following 
(plus some others which are not necessary for this example):

KielBiogeochem_NDOP_Matrix5_4/
Matrix5/
Matrix10/
MIT/
BiogeochemData/
grid.mat

The biogeochemical model and associated files are in KielBiogeochem_NDOP_Matrix5_4/

In Matlab, create a string variable 'base_path' that contains the full 
path to the TOPLEVEL directory. Save this variable in 'basepath.mat' within 
the TOPLEVEL directory.

Next, in KielBiogeochem_NDOP_Matrix5_4, make a soft link to basepath.mat, e.g.:
cd KielBiogeochem_NDOP_Matrix5_4
ln -s ../basepath.mat .

(5) In KielBiogeochem_NDOP_Matrix5_4, generate input data by running 
'prep_files_for_petsc_kiel_biogeochem_timestepping_cg.m'. You need to run this 
file twice, with different sets of parameters as follows:

(a) set the following parameters in prep_files_for_petsc_kiel_biogeochem_timestepping_cg.m 
and run:
writePCFiles=0
periodicForcing=1
periodicMatrix=1

This will write out the monthly mean transport matrices (Ae_00, Ae_01, ..., Ai_00, Ai_01, ...) 
and other monthly data (in PETSc and generic binary format) required to run the time-stepping 
biogeochemical model. (You will have to copy these files to the the remote machine 
on which the computation is actually being performed. See (7) below.)

To run the PETSc spinup and Jacobian code you will need to know the number of vertical 
profiles for your model grid. This information is computed when you run 
prep_files_for_petsc_kiel_biogeochem_timestepping_cg.m, provided the flag 
rearrangeProfiles is set to 1 (as it should be). What this does is 'rearrange' the 
transport matrices and all other variables by vertical profiles. Once you run this 
file, this number is written out to the screen. (You can also get the number of profiles via the 
variable 'numProfiles'.)

(b) 'clear all' and repeat with the following paramaeters: 
writePCFiles=1
periodicForcing=0
periodicMatrix=0

This will write out the annual mean transport matrices (Ae.petsc and Ai.petsc) and 
annual mean data required to run the time-stepping biogeochemical model in 
'Jacobian mode' (see 5b below).

(6) 'clear all' and run 'prep_data_for_spinup.m'. This will write out various files containing 
data required for the MFNK method (including information about the preconditioner).

If you are running the PETSc code (see (7) below) on a remote machine, then you need to copy 
the data generated by steps 4-6 to that machine. After completing step 6, I suggest you 
simply copy or rsync the entire directory KielBiogeochem_NDOP_Matrix5_4 over to the remote 
machine.

(7) Compile the various bits of PETSc code on the machine on which you will be running it. 
In general, you need 3 executable programs:
a) The time-stepping model (hereafter "spinup model")
b) A version of the above model used to compute the biogeochemical model Jacobian matrix
   (hereafter, "jacobian model"), and
c) A PETSc program to apply the implicit preconditioner

Since we need to call these programs multiple times to perform whatever 'action' 
they are written to perform, each time with a different input, all three programs 
are designed to run continuously. Each program waits for an external signal telling 
it to read the next input. It then reads the input, performs its action, writes 
the output, sends a signal that it is done, and then waits for the next input. 
There is a program-specific Matlab 'wrapper' function that does the actual signaling and 
exchange of input/output between Matlab and the running programs. 

Some specifics of the above:

(a)
The "spinup model" reads an initial condition as input, time steps the tracer 
equations for 1 period, and writes out the final state, and sends a signal 
to the user that it is done. Then, it waits for an external trigger from the 
user signalling it to read in the next initial condition, and so forth. 
Note that the difference between the final and initial condition is what we call 
a "function evaluation". We want to minimize the norm of this difference. 
In the current example, you build the executable with the command: 
'make tracercalckielbiondopspinup'

VERY IMPORTANT: do a 'make clean' before compiling each program!!!

Typically, you run the spinup model on a remote cluster (or at least something 
with a decent number of CPUs). Since this is the most expensive part of the MFNK 
calculation, the more resources you can throw at running the model, the better.
An example script for running the spinup model on a cluster is provided: 
runpbs_spinup0_allrelax_rzc

'calc_solution_with_running_time_stepper.m' is a wrapper function that takes a 
state vector (initial condition) as input, and returns the solution from the 
spinup model. 

(b)
The "jacobian model" is derived from the exact same time stepping code, but is 
compiled with a different CPP flag (see the Makefile). It works just like the 
spinup model, but instead of time stepping for a period, it takes a single time step, 
and instead of writing out the final model state, it writes the biogeochemical 
model tendency term, "q". We use repeated evaluations of q to compute the 
Jacobian dq/du (where u is the model state about which q is evaluated). Graph 
coloring is used to do this with a minimal number of "q-evaluations".
In the current example, you build the executable with the command: 
'make tracercalckielbiondopjacobian'

VERY IMPORTANT: do a 'make clean' before compiling each program!!!

Since it only takes a small number of time steps to compute the Jacobian, and as the 
Jacobian is only infrequently computed, you typically run the jacobian model on your 
local machine. Two example scripts are provided for running this code. To run locally, 
see: starttimestepperforjacobian_annualmean_allrelax. To run on a cluster, see:
runpbs_jacobiancalc_rzc

'calc_q_with_running_time_stepper.m' is a wrapper function that takes one (or more)  
state vector(s) (initial condition) as input, and returns the solution(s) from the 
Jacobian model. 

(c)
The 3d executable required is a PETSc program to apply the implicit preconditioner. 
This is essentially an implementation of eq. 10 of Ref [3], i.e., it computes 
w = P*inv(Q)*v, for a given v (see eqs. 10 and 11). This code too works like the 
above two in that it waits for an external trigger from the user signalling it to 
read in a new vector v, etc. (The matrices P and Q are computed using the dq/du 
computed by (b) above. See 'calc_preconditioner.m' for details.)
You build this executable with the command:
make directimplicitpc

VERY IMPORTANT: do a 'make clean' before compiling each program!!!

Typically, you run this program on a remote cluster. This program is very memory intensive 
(when using a sparse direct solver) so you will need substantial resources. 
An example script for running this code on a cluster is provided:
runpbs_directimplicitpc_rzc

'apply_direct_inverse_with_implicit_model_pc.m' is a wrapper function that takes a 
vector v as input, and returns the solution w computed by the preconditioner 
(i.e., the solution w to J^{\hat}*w = v, where J^{\hat} is the approximate Jacobian).

(8) Now, modify the main driver routine 'do_spinup.m' by defining the struct variables 
'spinupoptions', 'pcoptions', and 'jacobianoptions' that specify where/how the above 
executables (step 7) are being run. This is required so that the MFNK solver knows how to 
exchange information between Matlab and the three running programs. I use the ClusterTools 
toolbox to perform this communication (specifically to check whether the running job is ready 
to receive new data and if it is done writing data to file (so that it can be read into 
Matlab)). I use a script called 'make_cluster_options.m' to simply the task of constructing 
these variables. This file has hardwired information for the various machines on which I 
routinely run my code. I suggest you modify this file to match your own hardware. 
For now, however, you can use its 'default' mode to construct the appropriate 
options structures. Specifically, modify do_spinup.m to change the following three variables:
locDir, remDir, and remoteMachine, as indicated in the file. You also have a choice 
between selecting whether the PETSc programs are running on a different (remote) machine 
from the one on which Matlab is running (the default), or on the same machine.

IMPORTANT NOTE: 

The variable 'jacobianType' in do_spinup.m controls how the Jacobian is computed. There are 
several options for this, but for now you should set jacobianType to either 5 or 6. The 
difference is as follows:

jacobianType=5: 

This option computes dq/du by calling the 'q function' one input vector at a time. Thus, 
if the number of structurally independent columns of dq/du is M, the 'q function' is 
called M+1 times. Each time, the 'q function' sends this information to the running PETSc 
Jacobian model, reads the results and passes it back. If you select this option, you MUST 
set itout in the Jacobian model run script to '1'. However, this can be rather slow because of 
the overhead of communicating between the Matlab 'q function' and the running PETSc program. 
The preferred alternative is to use the following option:

jacobianType=6: 

This option will call the 'q function' ONCE with multiple vectors which are all passed on to the 
PETSc Jacobian model. Communication overhead is therefore minimized. If you select this option, 
you MUST set itout in the Jacobian model run script to M+1. The number M is available when you 
execute do_spinup.m in the variable 'nqvecs'. The number M+1 is also displayed to the screen 
when you execute do_spinup.m. Alternatively, to find out this number do the following in Matlab:

  load sparsity_data g
  nqvecs=length(unique(g));
  
(9) Now launch the various PETSc programs and wait until they have started up and 
are ready to recieve input.

(10) run 'do_spinup.m'