8.3 Framework for Joint Process Sets for Model and Filter
8.3.3 The Functionality of the Framework Routines
To gain further insight in the functionality of the data assimilation framework, we discuss here the operations which are performed in its main routines. The filter algo-rithms are hidden behind the three subroutines Filter Init, Get State, and Put State.
Due to this, the filter main routine, which was discussed in section 3.3.1, is split into two parts. These parts reside in Get State and Put State. Some additional operations are contained in these routines which are required for the parallel execution of the data assimilation framework.
The interface to the routine Filter Init has been shown in algorithm 8.1. Al-gorithm 8.4 sketches the operations which are performed in this routine when the
SEIK filter is used with a mode-decomposed ensemble matrix. The routine is called by all processes. Here several parameters are initialized, like the chosen filter algorithm or the ensemble size. These parameters are shared between the filter routines using Fortran modules. All subsequent operations in Filter Init are only performed by the filter processes. First, the sizes of sub-ensembles are computed. Subsequently, the arrays required for the filter are allocated. These are the state vector xand the local ensemble matrices Xp. In addition, the full ensemble matrix X is allocated on the fil-ter process with rank 0. Affil-ter the allocation of the fields, the user-supplied subroutine Init Ensemble is called. For the SEIK filter, this routine initializes the ensemble ma-trix. If a parallelization with mode-decomposition is used, Init Ensemble is only called by the process with rank 0. Here the full ensemble matrix is initialized. Subsequently, it is necessary to distribute sub-ensembles to all filter processes. This is performed by MPI communication operations.
In the case of domain-decomposed states, the routineInit Ensemble is called by all filter processes. The routine has to provide the full state ensemble for the local domain of each process. Since the state ensembles are readily initialized by all filter processes no further distribution of the ensembles is performed inFilter Init. A similar technique could be used for a mode-decomposed ensemble matrix. That is,Init Ensemble is called by each filter process with the local sub-ensemble as argument. Then Init Ensemble initializes only this local sub-ensemble. Since the sub-ensembles are readily initialized on the filter processes, no distribution of sub-ensembles would be required inFilter Init.
Using this variant would avoid the storage of the full ensemble matrix on a single pro-cess. On the other hand the user would be obliged to implement Init Ensemble such that all sub-ensembles are initialized correctly. From this point of view, the first vari-ant, which initializes the full ensemble matrix on a single process, is simpler to use.
If memory limitations render the allocation of the full ensemble matrix on a single process impossible, the initialization should directly operate on the sub-ensembles. To allow for this flexibility,Filter Init contains both variants.
The subroutineGet State is called prior to each model state evolution. Its structure is sketched in algorithm 8.5 for the SEIK and EnKF filters. If the routine is called for the very first time, it calls the user analysis routine User Analysis. This permits to analyze the initial ensemble consistently with the calls to User Analysis which are performed during the assimilation. Also the ensemble countermember is set to one at the very first call to Get State. For the remainder of the routine, this signals that a new forecast phase has to be performed.
IfGet State is called in the beginning of a forecast phase (i.e., with member = 1), the routine N ext Observation is called by the process of rank 0 in COMM FILTER.
N ext Observation initializes the number of time steps nsteps for the next forecast phase and the current model timetime. Subsequently, the value ofnstepsis distributed to all processes. If nsteps > 0, also the variable time is distributed to all processes by a broadcast operation. If the number of filter processes is smaller than the number of model tasks, as was the case in figure 8.3, the sub-ensemble of each filter process is further distributed such that each model task holds several ensemble members. This concludes the initialization of a forecast phase.
Subroutine Filter Init(. . .) ...
int mype f ilter {Rank of process in COMM FILTER} int npes f ilter {Number of processes in COMM FILTER} ...
1: initialize parameters
2: if f ilterpe== 1 then
3: initialize local ensemble sizesNp
4: allocate fields: Xp(n, Np),x(n)
5: if mype f ilter == 0 then
6: allocate ensemble matrix X(n, N)
7: call Init Ensemble(X) {Initialize full ensemble matrix}
8: for i= 1, npes f ilter do
9: send sub-ensemble X(jp :jp+rp−1) to filter process i {With MPI Send}
10: end for
11: deallocate field X
12: else if mype f ilter >0 then
13: receive sub-ensemble Xp {With MPI operation MPI Recv}
14: end if
15: end if
Algorithm 8.4: Sketch of the operations which are executed in the routineFilter Init for the case of mode-decomposition. The interface to this routine is shown as algorithm 8.1 When Get State is called during a forecast phase, it calls the user-supplied routine Distribute State. Here the model fields are initialized from the state vector which is provided to Distribute State as a subroutine argument. Since the state vector is only initialized on a single process of a model task, it might also be necessary to distribute the state information to the other processes of the model task.
Distribute Stateis not called directly by the model routines. Accordingly, the model fields or information on the model grid cannot be supplied as subroutine arguments.
Thus,Distribute State requires that the model fields are available via Fortran modules or ’common’ blocks. We will discuss this issue in section 8.5.
The routine Put State is called after a model state has been evolved by the model time stepper. Algorithm 8.6 sketches the operations which are performed in this routine for the SEIK filter. During the forecast, the user-supplied routineCollect State is called with the current ensemble state vector as argument. Also the ensemble countermember is incremented. Collect State initializes the forecasted state vector from the evolved model fields. This is the inverse operation to that performed by Distribute State. We will discuss Collect State in section 8.5.
If the forecast of all ensemble members is not yet finished, the program exits Put State and loops back to Get State in order to evolve the next ensemble member.
If the ensemble forecast is completed, the filter processes proceed in routine Put State
Subroutine Get State(. . .) ...
int f irsttime {Flag whether routine is called the very first time} int member {ensemble counter; shared using Fortran module} ...
1: if f irsttime == 1 then
2: call User Analysis(. . .)
3: f irsttime ←0
4: member ←1
5: end if
6: if member == 1 then
7: if mype f ilter == 0 then
8: call Next Observation(step,nsteps,time) {User supplied routine}
9: end if
10: broadcast nsteps to all processes {With operation MPI Bcast}
11: if nsteps >0then
12: broadcast time to all processes {With operation MPI Bcast}
13: distribute sub-ensembles {With operations MPI Send and MPI Recv}
14: end if
15: end if
16: if nsteps >0 then
17: call Distribute State(n,Xp(:, member)) {User supplied routine}
18: end if
Algorithm 8.5: Sketch of the operations which are executed in the routine Get State.
The interface to this routine is shown as algorithm 8.2
to perform the analysis and resampling phases of the filter algorithm. If there are less filter processes than model tasks, all ensemble members are gathered by the filter pro-cesses. Consecutively, the filter update phases are performed by callingSEIK Analysis and SEIK Resample and the user supplied analysis routine User Analysis. After the update, the ensemble counter member is reset to one and the filter processes exit Put State. Only the filter processes perform the update. The remaining processes re-set the ensemble counter and proceed directly to the routine Get State. Here, they wait to receive the variable nsteps which is send from the filter process with rank 0 in COMM FILTER to all processes by a broadcast operation (line 10 of algorithm 8.5).