Fortran

Workflow

Testing for reproducibility of the inflation algorithms has always been a bit problematic. Here is what I have done to see if the algorithms reproduce.

  1. Start with the lorenz_96 `workshop_setup.csh` script to generate the necessary input states. We are not interested in the output, we just need the `filter_input.nc`.
  2. Replace the obs_seq.out with the one attached here (`obs_seq.out.inf_test.txt`).
  3. Replace the input.nml with the one attached here (`input.nml.inf_test`)
  4. Edit the `inf_flavor` (and only the inf_flavor) for the test configuration - the RTPS test has one additional change, explained later.
  5. Run filter
  6. Make a directory for the output - I name them `inf24` for example - to reflect the fact the prior inflation was 2, the posterior inflation was 4.
  7. Move the filter output to the new directory.
  8. Repeat for other inflation configurations.
  9. Do the same thing in the 'gold standard' branch (Manhattan?)
  10. Edit the `CheckInflation.m` script to reference both sets of inflation directories and run it in matlab.


Comments

The `obs_seq.out.inf_test` observes state variable 10 at three times - separated by 600 seconds. I ran the obs_seq.in through `perfect_model_obs` and then manually biased the observations to create the necessary separation between the filter input and the observation. The observation error variance is set to be pretty small.  The `filter_output.nc`  has the values for the state-space inflation. `ncview` shows the inflation values (option 3 - spatially-constant - is kind of boring). The `CheckInflation.m` script compares all numeric variables in both files and writes out a file `Compare_netCDF_files_results.txt` with a complete summary. The Matlab command window will only show variables that are different, and will report the min/max differences.


RTPS

There is special consideration for testing the RTPS (Relaxation to prior spread, inf_flavor = 4). First, it can only be applied as posterior inflation. Second, the meaning (and value) of the `inf_initial` is distinctly different than for other inflation algorithms. The value should be between 0.0 and 1.0, I have tested with 0.3. Third, there is a bug in Manhattan (9.9.0) that causes the inflation values written to `filter_output.nc` to be incorrect. The algorithm is applied correctly, so the `obs_seq.final` should be considered the definitive answer. Since RTPS does not depend on the inflation values from previous cycles, the bug is believed to be a reporting error only. This reporting bug has been fixed, but comparison of 9.9.0 and current versions are expected to report something similar to that reported in the table below.



>> CheckInflation      
Comparing inflation case inf00 
                  MemberMetadata is a char variable.
                        inputnml is a char variable.

Comparing inflation case inf20 
                  MemberMetadata is a char variable.
                        inputnml is a char variable.

Comparing inflation case inf22 
                  MemberMetadata is a char variable.
                        inputnml is a char variable.

Comparing inflation case inf33 
                  MemberMetadata is a char variable.
                        inputnml is a char variable.

Comparing inflation case inf24 
                  MemberMetadata is a char variable.
                        inputnml is a char variable.
              state_postinf_mean has min/max differences of -3.22252 0.199984

Comparing inflation case inf55 
                  MemberMetadata is a char variable.
                        inputnml is a char variable.


The files used for testing and comparison are attached.


CheckInflation.m


Matlab

Discussion

The reproducibility between the Matlab implementation and the Fortran implementation is achieved by using identical input for both the Fortran and Matlab routines. The (adaptive_inflate_mod.f90) has a section of code in the update_inflation() routine that must be enabled and will then write out a tiny Matlab file with the variables needed to test the core matlab routines. There is a comment immediately prior to the block of code that must be enabled:


! Write minimum test data to a matlab file that can be explored with 'inflation_test.m'
! This test should only be run with the Lorenz_63 model ... and even then it will
! keep (over)writing the same file for each state variable. Since L63 only has 3 ...
! Also - there is no need to test for more than 1 cycle.
! Intentionally putting this before the constraint test.


The intent is to run the Lorenz_63 filter with the settings in input.inflationtest.nml using the input data from the workflow_setup.csh.  Run workflow_setup.csh as normal and then (essentially) modify the input.nml so that only the first assimilation cycle is used. The docs/DART_LAB/matlab/private/inflation_test.m script can then be used to explore the Matlab implementation.  The Fortran nugget writes out the solutions before the namelist bounds are enforced. (This is because the Matlab routine update_inflate() does not enforce the bounds, they are enforced by the Matlab calling routines.) There is no need to advance the model, and there is no need to run it for more than one location. Since it is a bit of a hack to write out the matlab nugget (Fortran creates a new file every entry into the update_inflation() routine) it is  efficient to do this as few times as possible - hence the L63 model with no advances.

Workflow

The testing process is to :

  1. enable adaptive_inflate_mod.f90 to write out inflation_input.m,
  2. compile the L63 filter
  3. replace the input.nml with input.inflationtest.nml
  4. run filter
  5.  run docs/DART_LAB/matlab/private/inflation_test.m
  6. choose another inflation algorithm in input.nml and repeat 4,5


Special note that Moha can expound upon ... there is one value needed for the Matlab implementation (ss_inflate_base) that is not directly attainable from the Fortran code that writes out the nugget. Consequently, the test is constructed such that the unavailable value (the prior inflation value before inflation) is known ... to be 1.0 ... by ensuring that the initial value of the inflation is 1.0 and that the model does not advance ... only a single (the first) timestep is used, filter does not advance the model.


  • No labels