Skip to content

Commit 2b4ac3e

Browse files
StevePnyStevePny
authored
Merge pull request #1 from UMD-AOSC/day2-dev
corrections from day1 sessions
2 parents 8af51c0 + 6b5ba6b commit 2b4ac3e

7 files changed

Lines changed: 23 additions & 217 deletions

generate_analysis_init.py renamed to analysis_init.py

Lines changed: 10 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -70,16 +70,16 @@
7070
# Set the linear observation operator matrix as the identity by default
7171
H = I
7272

73-
#print('B = ')
74-
#print(B)
75-
#print('R = ')
76-
#print(R)
77-
#print('H = ')
78-
#print(H)
73+
print('B = ')
74+
print(B)
75+
print('R = ')
76+
print(R)
77+
print('H = ')
78+
print(H)
7979

80-
das.setB(sigma_b**2*I)
81-
das.setR(sigma_r**2*I)
82-
das.setH(I)
80+
das.setB(B)
81+
das.setR(R)
82+
das.setH(H)
8383

8484
#-----------------------------------------------------------------------
8585
# Initialize the timesteps
@@ -123,7 +123,7 @@
123123
# Ensemble methods
124124
#-----------
125125
# EnKF
126-
method='ETKF'
126+
#method='ETKF'
127127
# Particle filter
128128
#method='PF'
129129

‎generate_nature_Jhist.py‎

Lines changed: 0 additions & 49 deletions
This file was deleted.

‎plot_analysis_plus_obs.py‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -12,12 +12,12 @@
1212
#----------------
1313
# Choose method:
1414
#----------------
15-
#method = 'skip'
15+
method = 'skip'
1616
#method = 'nudging'
1717
#method='OI'
1818
#method = '3DVar'
1919
#method = 'ETKF'
20-
method = 'PF'
20+
#method = 'PF'
2121
#method = 'Hybrid'
2222
infile=name+'_'+method+'.pkl'
2323
das = da_system()

‎plot_analysis_vs_analysis.py‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -9,14 +9,14 @@
99
#------------------------------------------------------------------
1010
name = 'x_analysis'
1111

12-
method1 = '3DVar'
12+
method1 = 'OI'
1313
infile1 = name+'_'+method1+'.pkl'
1414
das = da_system()
1515
das = das.load(infile1)
1616
sv1 = das.getStateVector()
1717
print(sv1)
1818

19-
method2 = 'Hybrid'
19+
method2 = '3DVar'
2020
infile2 = name+'_'+method2+'.pkl'
2121
das = da_system()
2222
das = das.load(infile2)

‎plot_analysis_vs_nature.py‎

Lines changed: 9 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -7,17 +7,17 @@
77
#------------------------------------------------------------------
88
# Read state vector objects
99
#------------------------------------------------------------------
10-
name = 'x_analysis'
11-
method = 'ETKF'
12-
infile1 = name+'_'+method+'.pkl'
13-
das = da_system()
14-
das = das.load(infile1)
15-
sv1 = das.getStateVector()
10+
infile1='x_nature.pkl'
11+
sv1 = state_vector()
12+
sv1 = sv1.load(infile1)
1613
print(sv1)
1714

18-
infile2='x_nature.pkl'
19-
sv2 = state_vector()
20-
sv2 = sv2.load(infile2)
15+
name = 'x_analysis'
16+
method = 'skip'
17+
infile2 = name+'_'+method+'.pkl'
18+
das = da_system()
19+
das = das.load(infile2)
20+
sv2 = das.getStateVector()
2121
print(sv2)
2222

2323
#------------------------------------------------------------------

‎tutorial_4.py‎

Lines changed: 0 additions & 144 deletions
Original file line numberDiff line numberDiff line change
@@ -21,147 +21,3 @@
2121
# (c) Change the forecast length
2222
#-----------------------------------------------------------------------
2323

24-
#-----------------------------------------------------------------------
25-
# Read the L63 nature run
26-
#-----------------------------------------------------------------------
27-
infile = 'x_nature.pkl'
28-
sv = state_vector()
29-
sv = sv.load(infile)
30-
x_nature = sv.getTrajectory()
31-
32-
#-----------------------------------------------------------------------
33-
# Initialize the timesteps
34-
#-----------------------------------------------------------------------
35-
t_nature = sv.getTimes()
36-
ainc_step = 1 # (how frequently to perform an analysis)
37-
dtau = (t_nature[ainc_step] - t_nature[0])
38-
tsteps=10 * ainc_step
39-
dt = dtau/tsteps
40-
maxit,xdim = np.shape(x_nature)
41-
42-
#-----------------------------------------------------------------------
43-
# Read the L63 observations
44-
#-----------------------------------------------------------------------
45-
infile = 'y_obs.pkl'
46-
obs = obs_data()
47-
obs = obs.load(infile)
48-
y_obs = obs.getVal()
49-
y_pts = obs.getPos()
50-
y_err = obs.getErr()
51-
print('y_obs = ')
52-
print(y_obs[0,:])
53-
print('y_pts = ')
54-
print(y_pts[0,:])
55-
56-
#-----------------------------------------------------------------------
57-
# Initialize the model
58-
#-----------------------------------------------------------------------
59-
l63 = lorenz63()
60-
61-
#-----------------------------------------------------------------------
62-
# Initialize the da system
63-
#-----------------------------------------------------------------------
64-
das = da_system()
65-
I = np.identity(xdim)
66-
67-
sigma_b = 0.9
68-
das.setB(sigma_b**2*I)
69-
print('B = ')
70-
print(das.B)
71-
72-
sigma_r = 1.0
73-
das.setR(sigma_r**2*I)
74-
print('R = ')
75-
print(das.R)
76-
77-
das.setH(I)
78-
print('H = ')
79-
print(das.H)
80-
81-
#-----------------------------------------------------------------------
82-
# Choose DA method:
83-
#-----------------------------------------------------------------------
84-
85-
#-----------
86-
# Session 0:
87-
# Test basic functionality
88-
#-----------
89-
#method='skip'
90-
91-
#-----------
92-
# Session 4:
93-
# Hybrid methods
94-
#-----------
95-
# Hybrid gain
96-
method='Hybrid'
97-
98-
das.setMethod(method)
99-
100-
#-----------------------------------------------------------------------
101-
# Initialize the ensemble
102-
#-----------------------------------------------------------------------
103-
xa = sv.x0
104-
bias_init = 0
105-
sigma_init = 0.1
106-
edim = 20
107-
Xa = das.initEns(xa,mu=bias_init,sigma=sigma_init,edim=edim)
108-
109-
#-----------------------------------------------------------------------
110-
# Test assimilation methods:
111-
#-----------------------------------------------------------------------
112-
xa_history = np.zeros_like(x_nature)
113-
KH_history = []
114-
for i in range(0,maxit,ainc_step):
115-
116-
#----------------------------------------------
117-
# Run forecast model for this analysis cycle:
118-
#----------------------------------------------
119-
t = np.arange(t_nature[i],t_nature[i]+dtau+dt,dt)
120-
# print('t = ', t)
121-
# Run the model ensemble forecast
122-
Xf = np.zeros_like(Xa)
123-
for k in range(edim):
124-
xf_4D = l63.run(Xa[:,k].A1,t)
125-
# Get last timestep of the forecast
126-
Xf[:,k] = np.transpose(np.matrix(xf_4D[-1,:]))
127-
128-
#----------------------------------------------
129-
# Get the observations for this analysis cycle
130-
#----------------------------------------------
131-
yo = y_obs[i,:]
132-
133-
#----------------------------------------------
134-
# Update the error covariances
135-
#----------------------------------------------
136-
# das.setB(sigma_b**2*I)
137-
# das.setR(sigma_r**2*I)
138-
# das.setH(I)
139-
140-
#----------------------------------------------
141-
# Compute analysis
142-
#----------------------------------------------
143-
Xa, KH = das.compute_analysis(Xf,yo)
144-
xa = np.mean(Xa,axis=1)
145-
146-
# print('Xa = ')
147-
# print(Xa)
148-
# print('xa = ')
149-
# print(xa)
150-
# print('xn = ')
151-
# print(x_nature[i,:].T)
152-
# print('KH = ')
153-
# print(KH)
154-
155-
# Archive the analysis
156-
xa_history[i,:] = xa.T
157-
158-
# Archive the KH matrix
159-
KH_history.append(KH)
160-
161-
# if ((i+1)%5==0):
162-
# exit()
163-
164-
sv.setTrajectory(xa_history)
165-
outfile='x_analysis_'+method+'.pkl'
166-
sv.save(outfile)
167-

‎tutorial_5.py‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,4 @@
11

2-
32
# Compute the Lyapunov exponents of the DA system
43

54

0 commit comments

Comments
 (0)