wiki:PDAF_debugging

PDAF Debugging Information

Implementation Guide

  1. Main page
  2. Adaptation of the parallelization
  3. Initialization of PDAF
  4. Modifications for ensemble integration
  5. Implementation of the analysis step
  6. PDAF-OMI, the Observation Module Infrastructure
  7. Memory and timing information
  8. PDAF debugging output
  9. Ensemble Generation
  10. Ensemble Diagnostics
  11. Auxiliary routines
  12. Generate synthetic observations
  13. Filter-specific options

Contents of this page

  1. Overview
  2. Activating Debugging Output
  3. Understanding the Debugging Output
    1. Debug output from PDAF_init
    2. Debug output for global ensemble filters
    3. Debug output for local ensemble filters
    4. Debug output for 3D-Var

Overview

The PDAF debugging output functionality was added with PDAF V2.1.

When implementing PDAF with your model or adding functionality like extending the state vector or adding support for new observations, it cab useful to check whether the inputs to the PDAF routines are correctly used inside PDAF. For this purpose, PDAF provides a debugging functionality. It allows you to activate debug output e.g. for a single local analysis domain on a single process of a complex application of a local filter like LEKSTF.

Usually, issues arise inside the user-supplied call-back routines or inside PDAF caused by the values that are provided by the user-supplied call-back routines. Thus, a particular aim of the PDAF debug output is to show values of arrays inside PDAF that depend on operations performed inside the user-supplied call-back routines.

Note on PDAF-OMI: When you use PDAF-OMI for the observation handling, there is a separate debug functionality for the observation part. For details see the documentation on PDAF-OMI debuggging. One can activate the debugging for PDAF and PDAF-OMI separately or can used both in combination.

Activating Debugging Output

Debugging output is activated by a call 

   INTEGER :: dbg_id  ! Debugging flag: >0 for debug output; =0 for no debug output

   CALL PDAF_set_debug_flag(dbg_id)

Setting the single argument of PDAF_set_debug_flag to a value larger 0 will active the output, while =0 will deactivate it.

The call to PDAF_set_debug_flag call can be inserted in different routines. Some recommended places are

  • init_pdaf: If the call is inserted before PDAF_init is called, one obtains debug output from PDAF_init. This allows one to check if the control parameters are handed over correctly and whether the values of the initial ensemble are correct inside PDAF
  • assimilate_pdaf: If the call is inserted before PDAF_assimilate_* or PDAF_put_state_* are called, one obtains debug output for the full analysis step. This variant is particularly useful for the global filters or the 3D-Var schemes.
  • init_dim_l_pdaf: If the call is used in this routine one can activate the debug output for a single local analysis domain. For example, to activate debug output only for the local analysis domain domain_p=4 and filter process 0 (mype_filter=0), one uses 
    SUBROUTINE init_dim_l_pdaf(step, domain_p, dim_l))

  USE mod_parallel_pdaf, ONLY: mype_filter

  ...

  IF (domain_p==4 .AND. mype_filter==0) THEN
    CALL PDAF_set_debug_flag(domain_p)
  ELSE
    CALL PDAF_set_debug_flag(0)
  ENDIF

  ...

For the debugging it is useful to keep the amount of output low. In particular for a large number of observations, a large ensemble or a large local state vector, the output will be very lengthy. To this end, it can be useful to intentionally reduce the size of the local state vector, number of observatios, or the ensemble size for the debugging. For the localized filters it is recommended to only activate the debugging for a single local analysis domain as shown above. The particular domain index can be chosen e.g. based on the coordinates of the domain, which are usually determined in init_dim_l_pdaf. Thus, if one encounters an issue in the analysis step for some coordinates one can find the corresponding coordinates in init_dim_l_pdaf and subsequently use PDAF_set_debug_flag to activate the debug output for this location.

For the domain-local filters, vectors like the local state vector are shown with their full length. Likewise vector or matrices whose dimension is the ensemble size are shown with their full length.

For the global filters, the state vector or the innovation can have a high dimension. To keep the debug output limited, the outputs of the global filters and the 3D-Var methods are generally limited to the first 6 elements of vectors that have the size of the state vector (like single ensemble states) or the size of the observation vector (like the innovation). Only vectors or matrices whose dimension is the ensemble size are shown with their full length.

Understanding the Debugging Output

The debug output prints information on configuration parameter values inside PDAF, dimensions and the values of actual arrays inside PDAF. Here, it depends for which routine or part of a data assimilation method the debugging is activated. We discuss a few options here for ensemble filters based on output obtain when activating debug output in the case /tutorial/offline_2D_serial. In addition, we discuss the debug output obtained for 3D-Var based on output obtain when activating debug output in the case /tutorial/3dvar/offline_2D_serial. WE discuss the debug output for an implementation using PDAF-OMI. However, its use is analogous for the case of the implementation not using PDAF-OMI.

Generally all debug output lines start with PDAF-debug, thus one can use grep PDAF-debug when running the data assimilation or on a text file containing the output to see only the debug output. Each line also shows the name of the PDAF routine from which the debug output is written.

Here we discuss some aspects of the debug output from different typical cases:

Debug output from PDAF_init

When we activate the PDAF debug output just before the call the PDAF_init and deactivate it directly afterwards (both in init_pdaf_offline.F90) and then run ./PDAF_offline -dim_ens 3|grep debug we obtain the following output (with line numbers added to support the description): 

 1: ++ PDAF-debug set_debug_flag: mype_world           0 activate           1
 2: ++ PDAF-debug:            1 PDAF_init -- START
 3: ++ PDAF-debug PDAF_init:           1 Use MPI_COMM_WORLD for COMM_PDAF
 4: ++ PDAF-debug PDAF_init:           1 param_int of size           7 values:         648           3           0           0           0           0           0
 5: ++ PDAF-debug PDAF_init:           1 param_real of size           1 values:   1.0000000000000000     
 6: ++ PDAF-debug PDAF_init:           1   Note: If REAL values appear incorrect, please check if you provide them with the correct precision
 7: ++ PDAF-debug PDAF_init:           1 ensemble member           1  values (1:min(dim_p,6)):  0.96592582999999999       0.99619469999999999       0.99619469999999999       0.96592582999999999       0.90630778999999995       0.81915203999999997     
 8: ++ PDAF-debug PDAF_init:           1 ensemble member           2  values (1:min(dim_p,6)):  0.99619469999999999       0.99619469999999999       0.96592582999999999       0.90630778999999995       0.81915203999999997       0.70710678000000005     
 9: ++ PDAF-debug PDAF_init:           1 ensemble member           3  values (1:min(dim_p,6)):  0.99619469999999999       0.96592582999999999       0.90630778999999995       0.81915203999999997       0.70710678000000005       0.57357643999999997     
10: ++ PDAF-debug:            1 PDAF_init -- END
11: ++ PDAF-debug set_debug_flag: mype_filter           0 deactivate

In the debug output we see the following:

  • Lines 1 and 11 show that the debug output is activated, respectively deactivated
  • Lines 2 and 10 mark the beginning and end of the debug output from the routine PDAF_init
  • Line 3 shows that MPI_COMM_WORLD is used for PDAF. This is the used in most cases (it can be change using a call to PDAF_set_comm_pdaf)
  • Lines 4 and 5 show the configuration parameter values (param_int and param_real. Here, in particular it can happen tha the value of the forgetting factor (param_real(1)) is displayed incorrectly. In this case it is likely that this REAL value is not specified with the correct precision or that the user routines and the PDAF library are not compiled for the same precision of REALs. (see the reminder in line 6)
  • Lines 7 - 9 show the first 6 elements of the 3 ensemble states. As for param_real it might be the values would be incorrect if the precision of reals in the user code and the PDAF library are inconsistent.

Debug output for global ensemble filters

To demonstrate the typical output from a global ensemble Kalman filter we use the debug output from the ESTKF. When we activate the PDAF debug output just before the call the PDAF_put_state_global in assimilation_pdaf_offline.F90 and then run ./PDAF_offline -dim_ens 3|grep debug we obtain a longer output. The first lines are 

 1: ++ PDAF-debug set_debug_flag: mype_filter           0 activate           1
 2: ++ PDAF-debug:            1 PDAF_estkf_update -- START
 3: ++ PDAF-debug PDAF_estkf_update:           1 ensemble member           1  forecast values (1:min(dim_p,6)):  0.96592582999999999       0.99619469999999999       0.99619469999999999       0.96592582999999999       0.90630778999999995       0.81915203999999997     
 4: ++ PDAF-debug PDAF_estkf_update:           1 ensemble member           2  forecast values (1:min(dim_p,6)):  0.99619469999999999       0.99619469999999999       0.96592582999999999       0.90630778999999995       0.81915203999999997       0.70710678000000005     
 5: ++ PDAF-debug PDAF_estkf_update:           1 ensemble member           3  forecast values (1:min(dim_p,6)):  0.99619469999999999       0.96592582999999999       0.90630778999999995       0.81915203999999997       0.70710678000000005       0.57357643999999997     
 6: ++ PDAF-debug PDAF_estkf_update           1 Configuration: param_int(3) dim_lag                0
 7: ++ PDAF-debug PDAF_estkf_update           1 Configuration: param_int(4) -not used-  
 8: ++ PDAF-debug PDAF_estkf_update           1 Configuration: param_int(5) type_forget            0
 9: ++ PDAF-debug PDAF_estkf_update           1 Configuration: param_int(6) type_trans             0
10: ++ PDAF-debug PDAF_estkf_update           1 Configuration: param_int(7) type_sqrt              0
11: ++ PDAF-debug PDAF_estkf_update           1 Configuration: param_int(8) observe_ens            F
12: ++ PDAF-debug PDAF_estkf_update           1 Configuration: param_real(1) forget        1.0000000000000000

In these lines we see that the output is written by PDAF_estkf_update which is an interal routine of PDAF controlling the analysis step of the ESTKF filter method. In the following we see

  • Lines 3-5 shows the first 6 elements of each forecast ensemble state. If these values are not realistic something failed suring the forecast phase and in collect_state_pdaf.
  • Lines 6-12 show the configuration parameters used in the filter. They should be consistent with what you provided to PDAF_init.

In the following lines we see the output from PDAF_estkf_analysis, which is the routine computing the actual analysis update: 

13: ++ PDAF-debug:            1 PDAF_estkf_analysis -- START
14: ++ PDAF-debug PDAF_estkf_analysis:           1   dim_p         648
15: ++ PDAF-debug:            1 PDAF_estkf_analysis -- call init_dim_obs
16: ++ PDAF-debug PDAF_estkf_analysis:           1   dim_obs_p          28
17: ++ PDAF-debug:            1 PDAF_estkf_analysis -- observe_ens F
18: ++ PDAF-debug:            1 PDAF_estkf_analysis -- call obs_op for ensemble mean
19: ++ PDAF-debug:            1 PDAF_estkf_analysis -- call init_obs
20: ++ PDAF-debug PDAF_estkf_analysis:           1 innovation d(1:min(dim_obs_p,10)) -0.49477208666666656       0.78162598333333333        1.9351931600000001       0.74712977333333319       0.81521189333333322        1.5003081066666666        1.1976581533333333       0.56667845666666650       0.87435301666666665        1.2913450866666667     
21: ++ PDAF-debug PDAF_estkf_analysis:           1 MIN/MAX of innovation  -1.2727408933333333        1.9351931600000001     
22: ++ PDAF-debug:            1 PDAF_estkf_analysis -- call obs_op           3 times
23: ++ PDAF-debug:            1 PDAF_estkf_analysis -- call prodRinvA_l
24: ++ PDAF-debug PDAF_estkf_analysis:           1   A^-1   5.1392136433237630       0.83883187576573037       0.83883187576573037        25.2334147876546644     
26: ++ PDAF-debug PDAF_estkf_analysis:           1   (HXT R^-1)^T d   11.810299761345098        2.6840904064519067     
27: ++ PDAF-debug:            1 PDAF_estkf_analysis -- type_sqrt           0
28: ++ PDAF-debug PDAF_estkf_analysis:           1   Compute eigenvalue decomposition of A^-1
29: ++ PDAF-debug PDAF_estkf_analysis:           1   eigenvalues   2.0086504639577534        5.3639779670206744     
30: ++ PDAF-debug PDAF_estkf_analysis:           1   A(HXT R^-1)^T d   2.2391867950151436       0.36078795190823887     
31: ++ PDAF-debug PDAF_estkf_analysis:           1   Omega^T  0.18261819498697551      -0.68154038207791268      -0.41706167726071641      -0.11175133961597850        0.0000000000000000        0.0000000000000000     
32: ++ PDAF-debug PDAF_estkf_analysis:           1   transform   2.1613652149063913      -0.52126615293618861       -1.6400990619702032        1.3571326901707161       0.47657822124568616       -1.8337109114164025        1.5507445396169155      -0.52126615293618861       -1.0294783866807269     
33: ++ PDAF-debug:            1 PDAF_estkf_analysis -- END

The lines 13-33 show different information:

  • Lines 14 and 16 show the dimension of the state vector and observation vector, respectively.
  • Line 18 (and some others) show an information like call obs_op for ensemble mean. This shows that an observation-related routine is called. When you use PDAF-OMI and in addition activate the debug output of PDAF-OMI you will see additional information about the observation-related values.
  • Line 20 shows the first 10 elements of the innovation vector, thus the difference between the vector of observations and the mean of the observed ensemble. Here one can check for extreme values which might result from incorrect observation or ensemble values.
  • Line 21 shows the minimum and maximum values of the innovation vector. Here, extreme values might also point to errors
  • Line 24 shows the inverse of matrix A of the ESTKF algorithm (see [wiki:PublicationsandPresentations Nerger et al., 2012 for the description of the ESTKF and ETKF algorithm in the notation used here). This is the full matrix of size dim_ens*dim_ens. If there are NaNs or extreme values, the analysis update will or can fail
  • Lines 26-31 show intermediate steps of the computation of the transform matrix. In particular, line 30 shows the vector that is used to transform the ensemble mean value from the forecast to the analysis. The shown steps allow you to trace back possible errors.
  • Lines 27 and 28 show that the eigenvalue decomposition of the inverse of matrix A is computed and the eigenvalues obtained. If the eigenvalue decomposition fails, there will be the error message PDAF-ERROR(1): Problem in computation of analysis weights!!!. In this case you should check the values of the inverse of matrix A which are shown in line 24.
  • Line 31 shows the actual transform matrix. Thus, the forecast ensemble matrix is multiplied with this matrix to obtain the analysis ensemble.

After the analysis ensemble is computed, the first 6 elements of each ensemble state are shown as: 

34: ++ PDAF-debug PDAF_estkf_update:           1 ensemble member           1  analysis values (1:min(dim_p,6)):  0.92068299395414255        1.0357490219605643        1.1193443475263585        1.1689289749905400        1.1829962688199513        1.1611187995043710     
35: ++ PDAF-debug PDAF_estkf_update:           1 ensemble member           2  analysis values (1:min(dim_p,6)):  0.94502620369513857        1.0416094338619108        1.1065438967701202        1.1378565941089185        1.1345960699625004        1.0968614132544858     
35: ++ PDAF-debug PDAF_estkf_update:           1 ensemble member           3  analysis values (1:min(dim_p,6)):  0.93916579179379245        1.0172662241209154        1.0644575418424744        1.0793058642259448        1.0613600135754058        1.0111652607977530     
36: ++ PDAF-debug:            1 PDAF_estkf_update -- END

This completes the analysis step of the ESTKF.

Debug output for local ensemble filters

To demonstrate the typical output from a local ensemble Kalman filter we use the debug output from the LESTKF. In this case we do not activate the PDAF debug output in assimilation_pdaf_offline.F90. Instead we activate the debug output in init_dim_l_pdaf for domain_p=4 as shown in the example on activating debug output above (we omit the check for mype_filter since the example is run using a single process). Then we run ./PDAF_offline -dim_ens 3 -filtertype 7 -cradius 4|grep debug. In this case we obtain 40 lines of output.

The first lines are: 

 1: ++ PDAF-debug set_debug_flag: mype_filter           0 activate           4
 2: ++ PDAF-debug PDAF_lestkf_update:           4   dim_l           1
 3: ++ PDAF-debug:            4 PDAF_lestkf_update -- call init_dim_obs_l
 4: ++ PDAF-debug PDAF_lestkf_update:           4   dim_obs_l           1
 5: ++ PDAF-debug:            4 PDAF_lestkf_update -- call g2l_state for ensemble member           1
 6: ++ PDAF-debug:            4 PDAF_lestkf_update --    Note: if ens_l is incorrect check user-defined indices in g2l_state!
 7: ++ PDAF-debug PDAF_lestkf_update:           4   ens_l  0.96592582999999999     
 8: ++ PDAF-debug:            4 PDAF_lestkf_update -- call g2l_state for ensemble member           2
 9: ++ PDAF-debug PDAF_lestkf_update:           4   ens_l  0.90630778999999995     
10: ++ PDAF-debug:            4 PDAF_lestkf_update -- call g2l_state for ensemble member           3
11: ++ PDAF-debug PDAF_lestkf_update:           4   ens_l  0.81915203999999997     
12: ++ PDAF-debug:            4 PDAF_lestkf_update -- call g2l_state for ensemble mean
13: ++ PDAF-debug PDAF_lestkf_update:           4   meanens_l  0.89712855333333319

These lines are analogous to the output for the global ESTKF. However, here the local dimensions are ensemble values are shown:

  • Lines 1 and 4 show the local state dimension (dim_l) and local observation dimension (dim_obs_l)
  • Lines 7-11 show the values of the local ensemble states
  • Line 13 shows the values of the local ensemble mean state

Note that the values of the configuration parameters are not shown if one activates the debug output in init_dim_l_pdaf, but only for the case that one activates the debug output in assimilate_pdaf (or assimilation_pdaf_offline). To get this information for the local filter while avoiding all the debug output form the local analysis, one can activate the debug ouptut in assimilate_pdaf and directly deactivate if in init_dim_l_pdaf.

In the following lines we see the output from PDAF_lestkf_analysis, which is the routine computing the actual local analysis update. 

 ++ PDAF-debug:            4 PDAF_lestkf_analysis -- START
 ++ PDAF-debug:            4 PDAF_lestkf_analysis -- call g2l_obs for mean
 ++ PDAF-debug:            4 PDAF_lestkf_analysis -- call init_obs_l
 ++ PDAF-debug PDAF_lestkf_analysis:           4   innovation d_l -0.49477208666666656     
 ++ PDAF-debug:            4 PDAF_lestkf_analysis -- call g2l_obs           3 times
 ++ PDAF-debug PDAF_lestkf_analysis:           4   HXT_l  0.16546108825519312        5.3415828255193198E-002
 ++ PDAF-debug:            4 PDAF_lestkf_analysis -- call prodRinvA_l
 ++ PDAF-debug PDAF_lestkf_analysis:           4   R^-1(HXT_l)  0.66184435302077249       0.21366331302077279     
 ++ PDAF-debug PDAF_lestkf_analysis:           4   A^-1_l   2.1095094869063713        3.5352964292627041E-002   3.5352964292627041E-002   2.0114130028327533     
 ++ PDAF-debug PDAF_lestkf_analysis:           4   (HXT_l R^-1)^T d_l -0.32746211159263749      -0.10571464322740090     
 ++ PDAF-debug:            4 PDAF_lestkf_analysis -- type_sqrt           0
 ++ PDAF-debug PDAF_lestkf_analysis:           4   Compute eigenvalue decomposition of A^-1_l
 ++ PDAF-debug PDAF_lestkf_analysis:           4   eigenvalues   2.0000000000000004        2.1209224897391241     
 ++ PDAF-debug PDAF_lestkf_analysis:           4   wbar_l -0.15439607679058354       -4.9843708923282673E-002
 ++ PDAF-debug PDAF_lestkf_analysis:           4   transform  0.54183477899890109      -0.34083296043308642      -0.20100181856581475      -0.44538532830038724       0.65993482760105171      -0.21454949930066455      -0.43015484539116355      -0.33915015825871242       0.76930500364987597     
 ++ PDAF-debug:            4 PDAF_lestkf_analysis -- END

These outputs are analogous to those for the global analysis, but for the vectors are matrices corresponding to the local analysis domain. A difference is that for the local analysis the local arrays are shown in their full length, while for the global analysis only the first six elements are shown.

After the analysis ensemble is computed, the full local analysis ensemble states are shown as: 

 ++ PDAF-debug:            4 PDAF_lestkf_update -- call l2g_state for ensemble member           1
 ++ PDAF-debug PDAF_lestkf_update:           4   ens_l  0.94695014510954834     
 ++ PDAF-debug:            4 PDAF_lestkf_update -- call l2g_state for ensemble member           2
 ++ PDAF-debug PDAF_lestkf_update:           4   ens_l  0.88927477553898138     
 ++ PDAF-debug:            4 PDAF_lestkf_update -- call l2g_state for ensemble member           3
 ++ PDAF-debug PDAF_lestkf_update:           4   ens_l  0.80443420998275139     
 ++ PDAF-debug:            4 PDAF_lestkf_update -- call l2g_state for ensemble mean
 ++ PDAF-debug PDAF_lestkf_update:           4   meanens_l  0.89712855333333319     
 ++ PDAF-debug:            4 PDAF_lestkf_update -- local analysis for domain_p           5
 ++ PDAF-debug:            4 PDAF_lestkf_update -- call init_dim_l
 ++ PDAF-debug set_debug_flag: mype_filter           0 deactivate

Debug output for 3D-Var

To demonstrate the typical output from 3D-Var methods we use the debug output from the 3D-Var with parameterized covariances. The outputs for ensemble 3D-Var and hybrid 3D-Var are similar but, in addition, contain output for the ensemble Kalman filter step used to update the ensemble.

When we activate the PDAF debug output just before the call the PDAF_put_state_3dvar in assimilation_pdaf_offline.F90 and then run ./PDAF_offline -dim_cvec 3|grep debug we obtain 277 lines of output. This long output is due to the iterations performed by the solver algorithm used to minimize the cost function.

The first lines come from before the itervative solver is executed: 

 1: ++ PDAF-debug set_debug_flag: mype_filter           0 activate           1
 2: ++ PDAF-debug:            1 PDAF_3dvar_update -- START
 3: ++ PDAF-debug PDAF_3dvar_update:           1  forecast state (1:min(dim_p,6)):  0.98610507666666669       0.98610507666666669       0.95614277333333331       0.89712855333333330       0.81085553666666654       0.69994508666666655     
 4: ++ PDAF-debug PDAF_3dvar_update           1 Configuration: param_int(3) solver                 1
 5: ++ PDAF-debug PDAF_3dvar_update           1 Configuration: param_int(4) dim_cvec               3
 6: ++ PDAF-debug PDAF_3dvar_update           1 Configuration: param_int(5) -not used-
 7: ++ PDAF-debug:            1 PDAF_3dvar_analysis -- START
 8: ++ PDAF-debug PDAF_3dvar_analysis:           1   dim_p         648
 9: ++ PDAF-debug:            1 PDAF_3dvar_analysis -- call init_dim_obs
10: ++ PDAF-debug PDAF_3dvar_analysis:           1   dim_obs_p          28
11: ++ PDAF-debug:            1 PDAF_3dvar_analysis -- call obs_op
12: ++ PDAF-debug:            1 PDAF_3dvar_analysis -- call init_obs
13: ++ PDAF-debug PDAF_3dvar_analysis:           1 innovation d(1:min(dim_obs_p,6)) -0.49477208666666656       0.78162598333333333        1.9351931599999999       0.74712977333333319       0.81521189333333322        1.5003081066666666     
14: ++ PDAF-debug PDAF_3dvar_analysis:           1 MIN/MAX of innovation  -1.2727408933333333        1.9351931599999999

This output provide the following information:

  • Line 3 shows the first 6 elements for the forecast (background) state vector
  • Lines 4-6 show the configuration parameters for the 3D-Var scheme.
  • Lines 8 and 10 show the dimension of the state vector and observation vector, respectively
  • Line 13 shows the first 6 elements of the innovation vector (observation - observed forecast state)
  • Line 14 shows the minimum and maximum values for the innovation vector.

For all these value one can check whether they are as pextect or show extreme values or NaNs

The next lines initialize the solver (in this case the LBFGS solver) and show the configuration variables for the chosen LBFGS solver (m, factr, pgtol). 

15: ++ PDAF-debug:            1 PDAF_3dvar_optim_LBFGS -- START
16: ++ PDAF-debug PDAF_3dvar_optim_LBFGS           1 Solver config: m               5
17: ++ PDAF-debug PDAF_3dvar_optim_LBFGS           1 Solver config: factr   10000000.000000000     
18: ++ PDAF-debug PDAF_3dvar_optim_LBFGS           1 Solver config: pgtol   1.0000000000000001E-005

Subsequently the iterations are performed. The debug output that is now shown shows steps in the calculation of the cost function.

For the first iteration the output is 

19: ++ PDAF-debug:            1 PDAF_3dvar_costf_cvt -- START: iteration           1
20: ++ PDAF-debug:            1 PDAF_3dvar_costf_cvt -- call cvt
21: ++ PDAF-debug PDAF_3dvar_costf_cvt:           1 state increment after CVT dX(1:min(dim_p,6))   0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000     
22: ++ PDAF-debug PDAF_3dvar_costf_cvt:           1 MIN/MAX dX   0.0000000000000000        0.0000000000000000     
23: ++ PDAF-debug:            1 PDAF_3dvar_costf_cvt -- call obs_op_lin
24: ++ PDAF-debug PDAF_3dvar_costf_cvt:           1 observed state after CVT (1:min(dim_obs_p,6))   0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000     
25: ++ PDAF-debug PDAF_3dvar_costf_cvt:           1 MIN/MAX HdX   0.0000000000000000        0.0000000000000000     
26: ++ PDAF-debug PDAF_3dvar_costf_cvt:           1 process local residual d (1:min(dim_obs_p,6))  0.49477208666666656      -0.78162598333333333       -1.9351931599999999      -0.74712977333333319      -0.81521189333333322       -1.5003081066666666     
27: ++ PDAF-debug PDAF_3dvar_costf_cvt:           1 MIN/MAX d  -1.9351931599999999        1.2727408933333333     
28: ++ PDAF-debug:            1 PDAF_3dvar_costf_cvt -- call prodRinvA
29: ++ PDAF-debug PDAF_3dvar_costf_cvt:           1 R^-1 d (1:min(dim_obs_p,6))   1.9790883466666662       -3.1265039333333333       -7.7407726399999994       -2.9885190933333328       -3.2608475733333329       -6.0012324266666663     
30: ++ PDAF-debug PDAF_3dvar_costf_cvt:           1 MIN/MAX R^-1 d  -7.7407726399999994        5.0909635733333332     
31: ++ PDAF-debug PDAF_3dvar_costf_cvt:           1 process local observation cost J_obs   46.955000904374110     
32: ++ PDAF-debug PDAF_3dvar_costf_cvt:           1 global observation cost J_obs   46.955000904374110     
33: ++ PDAF-debug PDAF_3dvar_costf_cvt:           1 process local CV cost J_B   0.0000000000000000     
34: ++ PDAF-debug PDAF_3dvar_costf_cvt:           1 global CV cost J_B   0.0000000000000000     
35: ++ PDAF-debug:            1 PDAF_3dvar_costf_cvt -- call obs_op_adj
36: ++ PDAF-debug PDAF_3dvar_costf_cvt:           1 H^TR^-1 d (1:min(dim_p,6))   0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000     
37: ++ PDAF-debug PDAF_3dvar_costf_cvt:           1 MIN/MAX H^TR^-1 d  -7.7407726399999994        5.0909635733333332     
38: ++ PDAF-debug:            1 PDAF_3dvar_costf_cvt -- call cvt_adj
39: ++ PDAF-debug PDAF_3dvar_costf_cvt:           1 CVT(H^TR^-1 d) (1:min(dim_cv_p,6))  -6.1852572643462160       0.26794725702686673        5.9173100073193554     
40: ++ PDAF-debug PDAF_3dvar_costf_cvt:           1 MIN/MAX CVT(H^TR^-1 d)  -6.1852572643462160        5.9173100073193554     
41: ++ PDAF-debug PDAF_3dvar_costf_cvt:           1 process local gradient gradJ (1:min(dim_cv_p,6))  -6.1852572643462160       0.26794725702686673        5.9173100073193554     
42: ++ PDAF-debug PDAF_3dvar_costf_cvt:           1 MIN/MAX gradJ  -6.1852572643462160        5.9173100073193554     
43: ++ PDAF-debug:            1 PDAF_3dvar_costf_cvt -- END

Noteable lines are here:

  • Lines 20, 23, 28, 35, 38 indicate that call-back routines are called. When PDAF-OMI is used one can activate the OMI debug output which will show information for the observation-related routines. cvt and cvt_adj are the control vector transform routine and its adjoint, which are both user-supplied. Thus, PDAF itself can only show the result returned from these routines.
  • Line 21 shows the first 6 elements of the state increment vector after appling the CVT. The values are all zero at iteration 1 since the initial increment is zero.
  • Lines 24 and 25 shows the observed increment vector. Again, all values are 0 at iteration 1
  • Lines 26 and 27 show the initial residual (observed increment minus innovation) and it minimum and maximum values
  • Lines 31 and 32 show the observation part of the cost function. In case of parallelization, line 31 shows the process-local part and line 32 the global sum
  • Lines 33 and 34 show the background part of the cost function, again for the process-local part and the global sum (both are zero at the initial iteration because no increment is computed so far)
  • Line 36 shows the residual scaled by the inverse observation error covariance matrix and projected back to the state space by the adjoint observation operator. Here the values are 0 because the first six elements of the state vector are not observed. (One can see that these values remain 0 during all iterations of the minimization.
  • Line 37 shows the minimum and maximum values of the vector for which line 36 shows the first 6 elements.
  • Line 39 shows the three elements of the vector after the adjoin CVT is applied (at maximum 6 elements are shown, but we run here with dim_cvec=3)
  • Lines 41 and 42 show the first 6 elements of the gradient of the cost function at its minimum and maximum values.

Particular aspects for debugging are here on the one hand the outcome of the observation-related operations (obs_op_lin, prodRinvA, obs_op_adj) and the outcome of the CVT and adjoint CVT. All these have to result in meaningful values.

The further lines of the debug output are for the different iterations of the minimization. Overall 10 iterations are computed durig which the value of J_obs decreases and the value of J_B increases. At the last iteration the elements in the gradient are reduced to an amplitude of avbout 10E-8: 

 ++ PDAF-debug PDAF_3dvar_costf_cvt:           1 MIN/MAX gradJ  -1.1445232317441878E-008   1.1542620637072787E-008

and the convergence is reached. At this point, the iterative solver is ended and the final control vector is projected back by the CVT to obtain the state vector increment as shown in the final output lines: 

 1: ++ PDAF-debug:            1 PDAF_3dvar_optim_LBFGS -- END
 2: ++ PDAF-debug PDAF_3dvar_analysis:           1 control vector (1:min(dim_p,6))   2.3896638011774280      -0.26679331415014473       -2.1228704870272916
 3: ++ PDAF-debug PDAF_3dvar_analysis:           1 MIN/MAX of control vector  -2.1228704870272916        2.3896638011774280     
 4: ++ PDAF-debug:            1 PDAF_3dvar_analysis -- call cvt for final state increment
 5: ++ PDAF-debug PDAF_3dvar_analysis:           1 state vector increment (1:min(dim_p,6))  -5.1146746761620225E-002   4.5436483221699978E-002  0.14063915510625419       0.23156859066274435       0.31546191351477898       0.38977007044061401     
 6: ++ PDAF-debug PDAF_3dvar_analysis:           1 MIN/MAX of state vector increment -0.55408418302780982       0.55408418302780982     
 7: ++ PDAF-debug:            1 PDAF_3dvar_analysis -- END
 8: ++ PDAF-debug PDAF_3dvar_update:           1  analysis state (1:min(dim_p,6)):  0.93495832990504646        1.0315415598883666        1.0967819284395874        1.1286971439960776        1.1263174501814455        1.0897151571072805     
 9: ++ PDAF-debug:            1 PDAF_3dvar_update -- END

The infromation shows is the following:

  • Lines 2 and 3 show the three elements of the final incremental control vector as well as its minmum and maximum values
  • Lines 5 and 6 show the first six elements of the state vector increment and its minimumand maximum values
  • Line 8 shows the first six elements of the analysis state vector

Other solvers
The debug outputs produced when using the other solvers are similar since the overall steps, like applying the CVT and adjoint CVT or obs_op_lin and obs_op_adj are similar.

Ensemble 3D-Var and Hybrid 3D-Var
For these ensemble-based variational methods the debug outputs are also similar to the outputs from the parameterized 3D-Var. However, there are now outputs relating to the ensemble-based CVT (cvt_ens and cvt_adj_ens). In addition there is will be debug output from the execution of the ESTKF or LESTKF filters that update the ensemble perturbations after the ensemble mean state is updated by the minimization of the cost function.

Last modified 10 months ago Last modified on Feb 5, 2024, 4:13:35 PM