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'