villanul
commented
2 months ago Hello guys,
I've worked a bit on that matter. Actually, maybe we don't have to do anything, non-linear updates seems good.
Consider a cell in an overlapping region. On the first sub-EnKF (noted 1), the value of the velocity at the cell is updated just like a normal EnKF. Now what happens for the next sub-EnKF (noted 2) for the two methods ?
The update for the second sub_EnKF looks like that : $x^{a2} = x^{a1} + K_2(y - s) = x^f + K_1(y - s) + K_2(y - s)$
with $s = Hx$.
The innovation term $(y - s)$ stays the same, only the Kalman gain changes. The value of the Kalman gain is given by : $K = SX/((S^2+R)$ in the case of an EnKF with 1 observation only. The only thing that is affected by a state update here is the anomaly matrix of the state $X$. Let's see the value of $X$ in $K_2$ for both cases :
-
Non-linear method (the one we use) : $X = x_i^{a1} - \overline{X^{a1}}/\sqrt{Ne-1}$
-
Addition of corrections (like intended at the beginning): $X = x_i^{f} - \overline{X^{f}}/\sqrt{Ne-1}$
In the Non-linear case, the term $x_i - \overline{x}$ is supposed to tend to 0 as the $x_i$ are supposed to get closer to the mean with each update. In the case of Addition of corrections, it never tends to 0.
So even if the covariance localization plays a great deal in weighting the values, in the addition of corrections the state may tend to a high number for a high number of observations, but for the Non-linear case the state tends to the ensemble mean for an infinite number of observations. On the other hand, addition of corrections is supposed to be how the EnKF works.
So i wanted to test the results for the channel flow RMSE. I was expecting slightly better results in the Non-linear case than for the Additions of corrections method, but not this weird and bad one, so you may consider i have a bug . If the results are real, Additions of corrections is incredibly bad and I did a genius mistake at the beginning but i don't think so :D
Giving those results i also tested doing non-Linear calculations for the parameter, but it's worse since without localization, the first sub-EnKF update gives a huge importance to the final value, individual sub-EnKF gives bad updates in the case of the parameter, doing the average of all the sub-EnKF updates like we do is much better.
If one of you have the time to verify this Addition of corrections behavior his own way it would be great !
Code snippet with the changes i made (comments are the previous piece of code, note that inflation don't work with this modification i did quickly, and i may have an error giving the results) :
Declare a correction matrix next to an existing declaration :
MatrixXf stateMatrixUpt = stateMatrix; // Save orginial state matrix to reconstruct it after hyperlocalization process
MatrixXf correctionsSum(nb_cells*3 + cwipiParams,cwipiMembers); // Save orginial state matrix to reconstruct it after hyperlocalization processModification of the loops and final addition :
//** Update temporary state
// temp_state_upt = temp_state + temp_K*(temp_obsMatrix - temp_sampMatrix);
temp_state_upt = temp_K*(temp_obsMatrix - temp_sampMatrix);
// if (cwipiVerbose) std::cout << "temp_state_upt = \n " << temp_state_upt << std::endl;
//** Perform inflation on temp_sate withtout parameter inflation **
temp_state_upt = inflation(temp_state_upt, cwipiMembers, nb_clipCells, cwipiParams, stateInfl, 0, typeInfl, paramEstSwitch, cwipiVerbose, stringRootPath);
// if (cwipiVerbose) std::cout << "temp_state_upt_infl = \n " << temp_state_upt << std::endl;
//** Save updated params to perform an average on the values updated by each clipping **
for (int j = 0; j < cwipiMembers; ++j){
for (int k = 0; k < cwipiParams; ++k){
temp_params(k, j) = temp_params(k, j) + temp_state_upt(3*nb_clipCells + k, j);
}
}
// if (cwipiVerbose) std::cout << "Params saved = " << temp_params << std::endl;
//** Update the values giving the updated temp state **
for (int j = 0; j < cwipiMembers; j++){
for (int k = 0; k < nb_clipCells; k++){
// if (cwipiVerbose) std::cout << "k when updating the values = " << k << std::endl;
cellIndex = clippingCells(k + count_cellsClipsDone + 1);
// if (cwipiVerbose) std::cout << "cellIndex when updating the values = " << cellIndex << std::endl;
// stateMatrixUpt(cellIndex, j) = temp_state_upt(k, j);
// stateMatrixUpt(nb_cells + cellIndex, j) = temp_state_upt(nb_clipCells + k, j);
// stateMatrixUpt(2*nb_cells + cellIndex, j) = temp_state_upt(2*nb_clipCells + k, j);
correctionsSum(cellIndex, j) = correctionsSum(cellIndex, j) + temp_state_upt(k, j);
correctionsSum(nb_cells + cellIndex, j) = correctionsSum(nb_cells + cellIndex, j) + temp_state_upt(nb_clipCells + k, j);
correctionsSum(2*nb_cells + cellIndex, j) = correctionsSum(2*nb_cells + cellIndex, j) + temp_state_upt(2*nb_clipCells + k, j);
}
}
if (i<(count_clippingCells-1)){ // If clippings are not done
count_cellsClipsDone = count_cellsClipsDone + nb_clipCells +1; // Update the count of the clipping already performed
nb_clipCells = clippingCells(i+1); // Update the number of cells in next clipping
count_obs = count_obs + 1; // New observation index
temp_state.resize(nb_clipCells*3 + cwipiParams,cwipiMembers); // Resize temp_state over the new number of cells in the clipping
}
}
else{
IDclippingCell = clippingCells(i+1);
//** Retrieve the correct values until reaching the next clipping **
if (cwipiVerbose) std::cout << "IDclippingCell = " << IDclippingCell << std::endl;
for (int j = 0; j < cwipiMembers; j++){
// temp_state(i - count_cellsClipsDone, j) = stateMatrixUpt(IDclippingCell, j);
// temp_state(nb_clipCells + i - count_cellsClipsDone, j) = stateMatrixUpt(nb_cells + IDclippingCell, j);
// temp_state(2*nb_clipCells + i - count_cellsClipsDone, j) = stateMatrixUpt(2*nb_cells + IDclippingCell, j);
temp_state(i - count_cellsClipsDone, j) = stateMatrix(IDclippingCell, j);
temp_state(nb_clipCells + i - count_cellsClipsDone, j) = stateMatrix(nb_cells + IDclippingCell, j);
temp_state(2*nb_clipCells + i - count_cellsClipsDone, j) = stateMatrix(2*nb_cells + IDclippingCell, j);
}
}
}
if (cwipiVerbose) std::cout << "Hyperlocalization loop done, beginning to finalize the parameters " << std::endl;
temp_params = temp_params/count_nbClips;
// if (cwipiVerbose) std::cout << "Params saved avg = " << temp_params << std::endl;
for (int j = 0; j < cwipiMembers; j++){
for (int k = 0; k < cwipiParams; k++){
stateMatrixUpt(3*nb_cells + k, j) = temp_params(k, j);
}
}
stateMatrixUpt = stateMatrixUpt + correctionsSum;
//** Perform inflation on stateMatrixUpt withtout state inflation **
stateMatrixUpt = inflation(stateMatrixUpt, cwipiMembers, nb_cells, cwipiParams, 0, paramsInfl, typeInfl, paramEstSwitch, cwipiVerbose, stringRootPath);
return stateMatrixUpt;
In the case of overlapping hyper-localization regions, DA is carried out in the second hyper-localization region by taking into account, the updated state matrix from the first hyper-localized region. This does not correspond to a linear addition of the contributions from different hyper-localization regions.