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'