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'