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
This file provides documentation for running a 'coarse grained' version of the model. It should
be useful for testing code and for learning about the MFNK method. I suggest you start by first
reading the file 'Readme' provided in this directory. It provides detailed instructions for
setting up the various pieces of code you will need. Work through sections (1)-(4) and then come
back and read the rest of this file. One key distinction between the instructions here and those
in Readme are that the coarse grained example, because it is so small, can be run entirely on your
local machine.
Once you have worked through sections (1)-(4) of Readme, do the following:
(5)
As a first step, you need to set the flag 'useCoarseGrainedMatrix' to '1' in the following
Matlab files:
prep_data_for_spinup.m
prep_files_for_petsc_kiel_biogeochem_timestepping_cg.m
That is, change the line 'useCoarseGrainedMatrix=0' to 'useCoarseGrainedMatrix=1'.
(6)
Now, generate grid and geometry information for the coarse grained (CG) model by running the
Matlab file 'do_coarse_graining.m'
This will write out a bunch of .mat files (grid_cg.mat and boxes_cg.mat) containing geometry
information for the CG grid. It will also write out 'coarse_graining_operators_cg.mat' which
contains sparse matrices that can convert variables between the coarse and original grids.
As configured, the CG grid averages in the horizontal over 64 boxes in x and 32 boxes in y.
(The original grid is 128x64.) No averaging is done in the vertical. You might think that this
will result in 4 large CG boxes in the horizontal, but since the function
calc_coarse_graining_operators.m takes into account the various ocean basins (i.e., you don't
want to cut across basins when averaging), in practice you will end up with more than that
(10 for this example). There will still be (a maximum of) 15 boxes in the vertical.
(7) 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.
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).
(8) '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).
(9) 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.
HOWEVER, for the CG example, you can run everything locally on your desktop/laptop.
An example script for running the spinup model is provided:
starttimestepperforspinup_allrelax_10profiles
'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. An example script is provided to run this code:
starttimestepperforjacobian_annualmean_allrelax_10profiles
'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. HOWEVER, the
CG example is small enough to be run locally. See the script: startdirectimplicitpc
'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).
(10) Now, modify the main driver routine 'do_spinup_cg.m'. For the CG example, you really
only need to change the following two variables: locDir and remDir as indicated in the file.
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));
(11) Now launch the various PETSc programs (e.g., using the following provided scripts) and wait
until they have started up and are ready to recieve input.
startdirectimplicitpc
starttimestepperforspinup_allrelax_10profiles
starttimestepperforjacobian_annualmean_allrelax_10profiles
(12) run 'do_spinup_cg.m'