Changes between Initial Version and Version 1 of Lorenz_63_model


Ignore:
Timestamp:
Jul 4, 2019, 3:19:56 PM (15 months ago)
Author:
lnerger
Comment:

--

Legend:

Unmodified
Added
Removed
Modified
  • Lorenz_63_model

    v1 v1  
     1= Lorenz-63 model with PDAF =
     2
     3[[PageOutline(2-3,Contents of this page)]]
     4
     5The implementation of the Lorenz-63 model coupled to PDAF is in the directory `testsuite/src/lorenz63/` of the PDAF package.
     6Provided is a full implementation of PDAF with the nonlinear Lorenz63 model (E. N. Lorenz (1963) Deterministic non-periodic flows. J. Atmos. Sci. 20, 130-141) providing various filter and smoother methods. The model has state dimenion 3. So it's not usable for
     7localized filters, but it's a good test case for the particle filter.
     8
     9Next to the implementation of the Lorenz-63 model with PDAF, the test case provides tool programs and scripts that
     10allow to run a test case and to display the outputs. Note, that this implementation runs without parallelization.
     11
     12The data assimilation with the Lorenz-63 model has been added in version 1.14 of PDAF.
     13
     14
     15== Running the test case  ==
     16
     17Runnning a data assimilation experiment with the Lorenz-63 model is a two step process: First one runs the model without PDAF to generate a file holding the trajectory of a forward run. Then one generates files with observations and a covariance matrix for the initialization of the initial ensemble. In the second step, one compiles the Lorenz-63 model with activated coupling to PDAF and runs the data assimilation experiments.
     18
     19=== 1. Compile and run the forward model without assimilation ===
     20
     21First change in the file `make.arch/linux_gfortran.h` the line
     22{{{
     23CPP_DEFS = -DUSE_PDAF
     24}}}
     25to
     26{{{
     27CPP_DEFS = #-DUSE_PDAF
     28}}}
     29to deactivate the coupling to PDAF in the model code. Now build the forward model program with
     30{{{
     31  cd testsuite/src
     32  make lorenz_63 PDAF_ARCH=linux_gfortran
     33}}}
     34in the directory `testsuite/src/` of the PDAF package. You have to ensure
     35that in the machine-specific make include file `linux_gfortran.h` `-DUSE_PDAF` is not defined
     36for CPP_DEFS (such that calls to PDAF are not active). You can replace `linux_gfortran` by any other make include file from `make.arch/`, e.g. specify `macos_gfortran` for compiling on MacOS. The executable
     37is generated in `testsuite/bin/`.
     38
     39'''Note''': The implementation uses the NetCDF library for file outputs. If the compilation above fails, please ensure the netcdf-library ist installed. On computers running Linux, it is usually available as a package of the operating system. On MacOS one can install the netcdf library e.g. using Fink or !MacPorts. NetCDF is a self-describing binary output format, but here it is not required that you know details about it. Anyway, if you like to look 'into' a NetCDF file, please try to use `ncdump FILENAME | less`. 
     40
     41To run the forward model use
     42{{{
     43  cd ../bin
     44  ./lorenz_63 -total_steps 10000
     45}}}
     46This runs the Lorenz-63 model for 10000 time steps and the trajectory is written into a file `state_l63.nc`.
     47
     48=== 2. Generate observations ===
     49
     50To build the executable to generate observations from the state trajectore use
     51{{{
     52  cd ../src/lorenz63/tools
     53  make generate_obs PDAF_ARCH=linux_gfortran
     54}}}
     55
     56Now run
     57{{{
     58  ./generate_obs
     59}}}
     60to generate a file holding observations (`obs_l63.nc` in testsuite/bin/).
     61
     62
     63=== 3. Build and run the assimilation program ===
     64
     65Now do
     66{{{
     67  cd ../../../../make.arch
     68}}}
     69and change in the file `linux_gfortran.h` the line
     70{{{
     71CPP_DEFS = #-DUSE_PDAF
     72}}}
     73back to
     74{{{
     75CPP_DEFS = -DUSE_PDAF
     76}}}
     77to activate the coupling to PDAF in the model code.
     78
     79Now compile the Lorenz-63 model with activated
     80coupling to PDAF. First clean the directories for the main driver and the Lorenz-63 model using
     81{{{
     82  cd ../testsuite/src
     83  make cleandriver PDAF_ARCH=linux_gfortran
     84  make cleanlorenz_63 PDAF_ARCH=linux_gfortran
     85}}}
     86(This removes object files that were compiled without support for PDAF)
     87Then build the executable using
     88{{{
     89  make pdaf_lorenz_63 PDAF_ARCH=linux_gfortran
     90}}}
     91The program `pdaf_lorenz_63` is generated in testsuite/bin.
     92
     93
     94To run the assimilation program, do
     95{{{
     96  cd ../bin
     97  ../src/lorenz63/tools/run_ESTKF.sh
     98}}}
     99The script run_ESTKF.sh runs an experiment with the ESTKF filter method for ensmeble size 20 in which only the variable X os the model is observed at each 10th time step. The execution should only take seconds.
     100Further you can run
     101{{{
     102  cd ../bin
     103  ../src/lorenz63/tools/run_PF.sh
     104}}}
     105to run a similar experiment with the particle filter.
     106
     107
     108
     109=== 4. Plot output from the assimilation experiments ===
     110
     111To display the output of the assimilation experiments we provide plotting scripts for Matlab, Octave and Python. To use them do
     112{{{
     113cd ../src/lorenz63/plotting/
     114}}}
     115
     116The following plot functions are available:
     117||= Script =||= Description =||
     118|| plot_eofs || Plot the mean state and the eigenvectors (EOFs) from the file holdig the covariance matrix that is generated with tools/generate_covar.F90. ||
     119|| plot_obs || Plot the observations for a specified time step. The observation file is generated using tools/generate_obs.F90. ||
     120|| plot_obs_series || Plot a time series of the observations for a specified variable of the model. The observation file is generated using tools/generate_obs.F90. ||
     121|| plot_rms || Plot the true and estimated RMS errors over time from an assimilation experiment. ||
     122|| plot_sigma || Plot the singular values from the file holdig the covariance matrix that is generated with tools/generate_covar.F90. ||
     123|| plot_state || Plot the state from the file holding a true state trajectory (state.nc) or the estimated forecast or analysis states from an assimilation experiment. ||
     124|| plot_state_series || Plot a time series of the state from the file holding a true state trajectory (state._l63nc) or the estimated forecast or analysis states from an assimilation experiment for a selected state variable. ||
     125
     126The scripts need the specification of the input file and the time step or variable to plot. The syntax of the functions is identical for Matlab and Python (except for `plot_rms`):
     127
     128||= Matlab =||= Python =||= Comment =||
     129|| plot_eofs(filename, index) || plot_eofs(filename, index) || index=0 for mean state, index>1 for EOFs ||
     130|| plot_obs(filename, timestep) || plot_obs(filename, timestep) || ||
     131|| plot_obs_series(filename, variable) || plot_obs_series(filename, variable) || ||
     132|| plot_rms(filename [, plot_forecast=1, plot_analysis=1]) || plot_rms(filename [, plot_forecast=True, plot_analysis=True]) || plot_forecast/plot_analysis are optional and switch on/off a plot ||
     133|| plot_sigma(filename) || plot_sigma(filename) || ||
     134|| plot_state(filename, timestep [, choice='t']) || plot_state(filename, timestep [, choice='t']) || `choice` is optional; options for `choice`: 't' true (use with truth file `state.nc`); for use with assimilation experiment output files: 'f' forecast, 'a' analysis, 'i' initial ||
     135|| plot_state_series(filename, variable [, choice='t']) || plot_state(filename, variable [, choice='t']) || `choice` is optional; options for `choice`: 't' true (use with truth file `state_l63.nc`); for use with assimilation experiment output files: 'f' forecast, 'a' analysis, 'i' initial ||
     136
     137Here `filename` is the name of the file including its path.
     138In Matlab use 'help' to display the information about required input.
     139
     140'''!Matlab/Octave plotting examples:'''[[BR]]
     141`plot_obs('../../../bin/obs.nc',100)` plots the observation at time step 100[[BR]]
     142`plot_obs_series('../../../bin/obs.nc',1)` plots the time series of observations for variable X (likely set 2=Y, 3=Z) [[BR]]
     143`plot_state('../../../bin/ESTKF_N20.nc',100,'f')` plots the forecast state estimate from the ESTKF run example at the 100th analysis step[[BR]]
     144`plot_rms('../../../bin/t1_N30_f0.97.nc')` plots the true and estimated RMS errors over time for the chosen experiment[[BR]]
     145`plot_state('../../../bin/state_l63.nc',1101)` plots the true state at model time step 1101 (= analysis step 100)[[BR]]
     146`plot_state_series('../../../bin/ESTKF_N20.nc',1,'f')` plots the time series for the ESTKF run example of the forecast state for variable X[[BR]]
     147
     148'''Python plotting'''[[BR]]
     149For Python the scripts are provided by `plot_l63.py`.
     150
     151The module can either be imported, e.g. to use its functions interactively:[[BR]]
     152`>>> import plot_l63`[[BR]]
     153`>>> plot_l63.plot_obs('../../../bin/obs_l63.nc', 4)`[[BR]]
     154
     155Alternatively the script can be run from the command line, providing the function
     156name and its argument as command line parameters:[[BR]]
     157`./plot_l63.py plot_obs ../../../bin/obs_l63.nc 4`[[BR]]
     158(If this fails you can also try to run the script as
     159`python plot_l63.py plot_obs ../../../bin/obs_l63.nc 4`). To plot only the Analysis RMS error one can use[[BR]]
     160`plot_rms ../../../bin/ESTKF_N20.nc False True`
     161
     162
     163== Run options ==
     164
     165The implementation for the Lorenz-63 model has a wide range of options. For the full set, we recommend to check the list in the file init_pdaf.F90 in the directory lorenz63. Essentially all of the options can be specified on the command line in the format '-VARIABLE VALUE', where VARIABLE is the name of the variable in the program and VALUE is its value. The experiment in `run_ESTKF.sh` is executed like
     166{{{
     167./pdaf_lorenz_63 -total_steps 5000 -step_null 1000 -dim_ens 20 -forget 0.8 -filtertype 6 -file_asml ESTKF_N20.nc
     168}}}
     169and the experiment in `run_PF.sh` is executed like
     170{{{
     171./pdaf_lorenz_63 -total_steps 5000 -step_null 1000 -dim_ens 20 -pf_res_type 2 -pf_noise_type 2 -pf_noise_amp 0.2 -filtertype 12 -file_asml PF_N20.nc
     172}}}
     173The meaning of the different options is the following:
     174
     175||= Variable =||= Description =||= Default value =||
     176|| dim_ens || ensemble size || 20
     177|| file_asml || the name of the output file || ESTKF_N20.nc ||
     178|| filtertype || Specifies the filter to use; 6 sets the ESTKF filter, 12 the particle filter || 6 ||
     179|| forget  || value of the forgetting factor for multiplicative variance inflation || 1.0 ||
     180|| step_null || initial time step of experiment. This specifies the offset of the observations which are generated from a model forward run. A value 1000 avoids the spin-up phase of the model || 0 ||
     181|| total_steps || number of time step in the experiment || 5000 ||
     182|| pf_res_type || type of resmpling scheme in the particle filter (2 set stochastic unversal resampling) || 1 ||
     183|| pf_noise_type || type of noise to perturb resampled particles (2 sets noise with amplitude relative to ensemble standard deviation) || 0 ||
     184|| pf_noise_amp || amplitude of noise to perturb resampled particles || 0.0 ||
     185
     186Some further options:
     187
     188||= Variable =||= Description =||= Default value =||
     189|| delt_obs  || time interval between observations, i.e. the forecast length || 10 ||
     190|| dim_lag   || set the time lag for the smoother (the smoother is active for dim_lag>0) || 0 ||
     191|| dx_obs    || distance between observations on the grid; let's you specify an incompletely observated state (only used when `use_obs_mask=.true.`) || 3 ||
     192|| numobs || number of observed grid points; the points 1 to numobs are observed (can be combined with dx_obs; only used when `use_obs_mask=.true.`) || 1 ||
     193|| use_obs_mask || whether to use incomplete observations (need to be .true. for any other settings on the observation density to be used) || .true. ||
     194
     195For the full set of options, please see init_pdaf.F90 in the lorenz63 directory. There, also the possible settings e.g. for `filtertype` are described.
     196
     197