-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathrun_cycle_filter.py
executable file
·119 lines (102 loc) · 3.77 KB
/
run_cycle_filter.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
#!/usr/bin/env python
import numpy as np
import L96_model as L96
import data_assimilation as DA
import misc
import sys
###CONFIG
###random seed fixed so that results are repeatable
np.random.seed(0)
### Model setup
nx = 40
dt = 0.2
F = np.round(float(sys.argv[5]), 2) #8.0
### Cycling experiment setup
nt = 100000
cycle_period = 1
time_window = 0 ##smoother analysis window (+-cycles)
time_space_ratio = 1.0 ##ratio of dt/dx
### Observation network setup
obs_err = np.round(float(sys.argv[3]), 2) #1.0
L = np.round(float(sys.argv[2]), 2) #spatial corr in R
obs_thin = 1
obs_ind = np.tile(np.arange(0, nx, obs_thin), (nt, 1)).T
### Ensemble filter tuning parameters
filter_kind = sys.argv[1]
nens = int(sys.argv[4])
ROI = int(sys.argv[6]) #localization in space (grid points)
alpha = 0.0 #np.round(float(sys.argv[7]), 2) #0.0 ##relaxation coef
inflation = np.round(float(sys.argv[7]), 2) ##multiplicative inflation
#for MS algorithms
dk = int(sys.argv[8])
if (dk < int(nx/2)) :
krange = np.arange(dk, int(nx/2), dk)
R = DA.R_matrix(nx, obs_ind, np.array([0]), 1, 5, 0)
Lo = misc.matrix_spec(R)
obs_err_inf = np.ones(krange.size+1)
obs_err_inf[0] = np.sqrt(np.mean(Lo[0:krange[0]+1]**2))
obs_err_inf[krange.size] = np.sqrt(np.mean(Lo[krange[-1]+1:]**2))
for i in range(1, krange.size):
obs_err_inf[i] = np.sqrt(np.mean(Lo[krange[i-1]+1:krange[i]+1]**2))
ROI_adjust = np.ones(krange.size+1)
casename = filter_kind+"/dk{}".format(dk)+"/L{}_s{}".format(L, obs_err)+"/N{}_F{}".format(nens, F)+"/ROI{}".format(ROI)+"_inf{:4.2f}".format(inflation)
print(casename)
# read in initial data
xt = np.load("output/truth.npy")
obs = np.load("output/obs.npy")
xens0 = np.load("output/initial_ensemble.npy")
xens = np.zeros([nx, nens, nt])
xens[:, :, 0] = xens0[:, 0:nens]
xens1 = np.copy(xens)
##### assimilation cycle
for tt in range(nt-1):
# analysis step
if(np.mod(tt, cycle_period) == 0) and tt>0:
xb = xens1[:, :, tt].copy()
#####prior inflation
xb_mean = np.mean(xb, axis=1)
for k in np.arange(nens):
xb[:, k] = inflation*(xb[:, k]-xb_mean) + xb_mean
######Define obs network
t1 = max(0, tt-time_window)
t2 = min(nt, tt+time_window)
t_ind = np.arange(t1, t2+1)
# print(t_ind)
x = xb.copy()
dx = np.zeros(x.shape)
##assimilate obs from within window
for t in t_ind:
# print(t)
yo = obs[:, t]
H = DA.H_matrix(nx, obs_ind, np.array([t]), 0)
R = DA.R_matrix(nx, obs_ind, np.array([t]), obs_err, L, 0)
rho = DA.local_matrix(nx, np.array([tt]), ROI, 0)
###run filter
if filter_kind == "EnKF":
x1 = DA.EnKF_perturbed_obs(x, yo, H, R, rho, tt)
if filter_kind == "EnSRF":
x1 = DA.EnSRF(x, yo, H, R, rho)
if filter_kind == "EnSRF_serial":
x1 = DA.EnSRF_serial(x, x, yo, obs_err, obs_ind[:, t], ROI)
if filter_kind == "MultiscaleObs":
x1 = x.copy()
x_s = np.zeros((nx, nens))
for s in range(krange.size+1):
for k in range(nens):
x_s[:, k] = misc.spec_bandpass(x1[:, k], krange, s)
yo_s = misc.spec_bandpass(yo, krange, s)
x1 = DA.EnSRF_serial(x1, x_s, yo_s, obs_err*obs_err_inf[s], obs_ind[:, t], ROI*ROI_adjust[s])
dx += x1 - x
xa = xb + dx
#####posterior inflation/relaxation
xb_mean = np.mean(xb, axis=1)
xa_mean = np.mean(xa, axis=1)
for k in np.arange(nens):
xens1[:, k, tt] = (1-alpha)*(xa[:, k]-xa_mean) + alpha*(xb[:, k]-xb_mean) + xa_mean
# xens1[:, k, tt] = inflation*(xa[:, k]-xa_mean) + xa_mean
# forecast step
xens1[:, :, tt+1] = L96.M_nl(xens1[:, :, tt], nx, F, dt)
xens[:, :, tt+1] = L96.M_nl(xens1[:, :, tt], nx, F, dt)
outdir = "/glade/scratch/mying/L96_DA"
np.save(outdir+"/"+casename+"/ensemble_prior", xens)
np.save(outdir+"/"+casename+"/ensemble_post", xens1)