From dba1b065acbee3070ef1158335a0f717a9019d5e Mon Sep 17 00:00:00 2001
From: Philipp Griewank <philipp.griewank@uni-koeln.de>
Date: Mon, 17 May 2021 18:21:28 +0200
Subject: [PATCH] Added Gaspari Cohn localization
Keep in mind that in comparison to the gaussian the gaspari cohn is slimmer than the gaussian with the same loc_length, so they shouldn't be compared directly against each other.
---
da_functions.py | 34 ++++++++++++++++++++++++++++------
1 file changed, 28 insertions(+), 6 deletions(-)
diff --git a/da_functions.py b/da_functions.py
index 97299fd..06e0e50 100644
--- a/da_functions.py
+++ b/da_functions.py
@@ -262,8 +262,7 @@ def run_linear_advection_EnKF(m_const,da_const):
#Construct localization matrix C if loc!=None
if da_const["loc"]:
- if da_const["loc_type"]=='gaussian':
- C = gaussian_loc_matrix(da_const,m_const)
+ C = loc_matrix(da_const,m_const)
else:
C=np.ones([m_const["nx"],m_const["nx"]])
@@ -712,12 +711,35 @@ def var_reduction_estimate_iterative(states,m_const,da_const,j=0,dJ_update_flag=
-def gaussian_loc_matrix(da_const,m_const):
+def loc_matrix(da_const,m_const):
+ """
+ Creates the localization matrix, using either a gaussian or gaspari cohn function
+ """
C = np.zeros([m_const['nx'],m_const['nx']])
- #The gaussian is slightly less trivial because of the periodic boundary domain, so what I do is I use the same function that i use to generate the initial conditions and then step it once for all other locations
- #And of course i am not sure if i have my dimensions switched
- C[:,0] = gaussian_initial_condition(m_const["x_grid"],da_const["loc_length"])
+ # I cheat a bit by mirroring the functions to accomodate for the repeating boundary conditions, but this should only lead to a maximum dx/2 error.
+ if da_const['loc_type']=='gaussian':
+ C[:,0] = gaussian_initial_condition(m_const["x_grid"],da_const["loc_length"])
+ if da_const['loc_type']=='gaspari_cohn':
+ C[:,0] = gaspari_cohn(m_const["x_grid"],da_const["loc_length"])
+
for i in range(1,m_const['nx']):
C[:,i] = np.hstack([C[-1,i-1],C[:-1,i-1]])
return C
+def gaspari_cohn(x,loc_length):
+ """Gaspari-Cohn function."""
+
+ ra = x/loc_length
+ gp = np.zeros_like(ra)
+ i=np.where(ra<=1.)[0]
+ gp[i]=-0.25*ra[i]**5+0.5*ra[i]**4+0.625*ra[i]**3-5./3.*ra[i]**2+1.
+ i=np.where((ra>1.)*(ra<=2.))[0]
+ gp[i]=1./12.*ra[i]**5-0.5*ra[i]**4+0.625*ra[i]**3+5./3.*ra[i]**2-5.*ra[i]+4.-2./3./ra[i]
+
+ #Now we mirror things to zero for the periodic boundary domain
+ half_idx = int(ra.shape[0]/2)
+ gp[-half_idx:] = gp[half_idx-1::-1]
+ return gp
+
+
+
--
GitLab