pushing IC files for longer simulations

Author: ray-chew

RKLM_Python/data_assimilation/params.py

@@ -13,10 +13,11 @@ def __init__(self,N,da_type='rloc'):
 
         self._da_times = np.arange(0.0,3.25,0.25)[1:]
         # self._da_times = np.arange(5.0,10.5,0.5)/10.0
+        # self._da_times = [0.1]
         self._da_times = np.around(self.da_times,3)
 
         self.obs_attributes = ['rho','rhou', 'rhov', 'rhoY', 'p2_nodes']
-        # self.obs_attributes = ['rho']
+        # self.obs_attributes = ['rhou','rhov']
 
         # which attributes to inflate in ensemble inflation?
         self.attributes = ['rho', 'rhou', 'rhov']
@@ -45,7 +46,7 @@ def __init__(self,N,da_type='rloc'):
         if self.sparse_obs_by_attr:
             assert(0, "Not yet implemented.")
 
-        self.obs_frac = 0.10 # fraction of the observations to pick.
+        self.obs_frac = 0.90 # fraction of the observations to pick.
         self.gen_obs_sparse()
 
         ############################################

RKLM_Python/inputs/travelling_vortex_2D.py

@@ -134,7 +134,8 @@ def __init__(self):
         if self.continuous_blending == True:
             self.output_suffix = "_%i_%i_%.1f" %(self.inx-1,self.iny-1,self.tout[-1])
 
-        aux = 'pos_debug'
+        aux = 'analytical_truth'
+        aux = 'debug_cov'
         self.aux = aux
 
         self.output_suffix = "_%i_%i_%.1f_%s" %(self.inx-1,self.iny-1,self.tout[-1],aux)
@@ -166,11 +167,11 @@ def sol_init(Sol, mpv, elem, node, th, ud, seed=None):
 
     # gradient wind eqn:
     # p = Int[( alpha_const * 1/2 * rho0 + alpha * 1/2 * rho(r)) * u^2 / r, r]
-    # factors of 1/2 * rho(r)
+    # factors of 1/2 * rho(r) * u^2 / r
     alpha = -1.0
     # factors of 1/2 rho0 * u^2 / r
-    alpha_const = 3.0
-    rho0 = 1.5
+    alpha_const = 2.0
+    rho0 = 1.0
     del_rho = -0.5
     R0 = 0.4
     fac = 1. * 1024.0
@@ -181,13 +182,13 @@ def sol_init(Sol, mpv, elem, node, th, ud, seed=None):
         np.random.seed(seed)
         xc += (np.random.random() - 0.5) / 5.0
         yc += (np.random.random() - 0.5) / 5.0
-        print(seed, xc, yc)
+
 
     if 'truth' in ud.aux or 'obs' in ud.aux:
         np.random.seed(2233)
         xc += (np.random.random() - 0.5) / 5.0
         yc += (np.random.random() - 0.5) / 5.0
-        print(seed, xc, yc)
+
         ud.xc = xc
         ud.yc = yc
 

RKLM_Python/results_gen.py

@@ -42,7 +42,7 @@
 ##########################################
 # generate initial blending for section 5a
 gen_5a1_euler = False
-gen_5a2_rb = False
+gen_5a2_rb = True
 # generate obs and truth for section 5b
 gen_5b_obs_truth_euler = False
 # generate Euler ensemble simulations for
@@ -60,7 +60,7 @@
 gen_5b_5c2_rb = False
 
 gen_loc_errors_tv = False
-gen_loc_errors_rb = True
+gen_loc_errors_rb = False
 
 # Otherwise, if gen_all = True, generate all
 # results
@@ -90,6 +90,8 @@
     # Do not output in between, but output
     # after each timestep.
     tout = [1.0]
+    # tout = np.arange(0.0,1.005,0.005)[1:]
+    # tout = tout.tolist()
 
     # JSON dumping does not accept ndarrays.
     # tout = tout.tolist()
@@ -289,6 +291,8 @@
     rp.N = 1
     rp.tc = 'rb'
     tout = [1.0]
+    # tout = np.arange(0.05,1.05,0.05)
+    # tout = list(tout)
     # No data assimilation.
     dap = {
         'None' : None,
@@ -447,7 +451,7 @@
     # run simulation
     rp.ud = json.dumps(ud)
     rp.dap = json.dumps(dap)
-    # rp.queue_run()
+    rp.queue_run()
 
     ##########################################
 
@@ -467,7 +471,7 @@
     # run simulation
     rp.ud = json.dumps(ud)
     rp.dap = json.dumps(dap)
-    # rp.queue_run()
+    rp.queue_run()
 
 
 if gen_5c1_euler_full or gen_all:
@@ -480,8 +484,8 @@
 
     # da_times = np.arange(5.0,10.5,0.5)/10.0
     # da_times = np.array([0.25,10.0])
-    da_times = np.arange(0.0,10.25,0.25)[1:]
-    # da_times = np.arange(0.0,3.25,0.25)[1:]
+    # da_times = np.arange(0.0,10.25,0.25)[1:]
+    da_times = np.arange(0.0,3.25,0.25)[1:]
     da_times = np.around(da_times,3)
     da_times = da_times.tolist()
 
@@ -501,7 +505,7 @@
         # Assimilate all fields
         'obs_attrs' : ['rho', 'rhou', 'rhov', 'rhoY', 'p2_nodes'],
         # Path to the generated observation
-        'obs_path' : path_to_obs + 'output_travelling_vortex/output_travelling_vortex_ensemble=1_64_64_10.0_obs.h5',
+        'obs_path' : path_to_obs + 'output_travelling_vortex/output_travelling_vortex_ensemble=1_64_64_3.0_obs.h5',
         'loc_setter' : (11,11)
     }
 

RKLM_Python/tv_enda_tests.py

@@ -46,9 +46,10 @@
 # section 5c1 with assimilation of all
 # fields
 generate_noda = False
-gen_5c1_euler_full = False
+generate_oneda = False
+gen_5c1_euler_full = True
 gen_5c1_euler_momenta = False
-gen_loc_errors_tv = True
+gen_loc_errors_tv = False
 
 # Otherwise, if gen_all = True, generate all
 # results
@@ -58,7 +59,9 @@
 # path_to_obs = '/srv/public/ray/'
 path_to_obs = './'
 
-enda_sfx = 'wda_varcov'
+enda_sfx = 'wda_obsconv'
+noda_sfx = ''
+ref_aux = '_truthgen'
 
 if gen_5b_obs_truth_euler or gen_all:
     ##########################################
@@ -76,7 +79,8 @@
 
     # simulation parameters for the observation
     ud = {
-        'aux' : 'obs',
+        'aux' : 'obs%s' %ref_aux,
+        'initial_blending' : False
     }
 
     # data assimilation parameters for the
@@ -94,7 +98,7 @@
 
     # simulation parameters for the truth
     ud = {
-        'aux' : 'truth',
+        'aux' : 'truth%s' %ref_aux,
         # Do blending for the initial time-step
         'initial_blending' : True
     }
@@ -113,8 +117,8 @@
     #
     ##########################################
 
-    da_times = np.arange(0.0,10.25,0.25)[1:]
-    da_times = np.arange(0.0,3.25,0.25)[1:]
+    da_times = np.arange(0.0,10.05,0.05)[1:]
+    # da_times = np.arange(0.0,3.25,0.25)[1:]
     da_times = np.around(da_times,3)
     da_times = da_times.tolist()
 
@@ -125,7 +129,7 @@
 
     # Generate ensemble with no DA
     ud = {
-        'aux' : 'noda',
+        'aux' : 'noda%s' %noda_sfx,
         # Do blending for initial time-step
         'initial_blending' : True
     }
@@ -137,7 +141,7 @@
         'obs_attrs' : ['rho', 'rhou', 'rhov', 'rhoY', 'p2_nodes'],
         'loc_setter' : (11, 11),
         # Path to the generated observation
-        'obs_path' : path_to_obs + 'output_travelling_vortex/output_travelling_vortex_ensemble=1_64_64_3.0_obs.h5'
+        'obs_path' : path_to_obs + 'output_travelling_vortex/output_travelling_vortex_ensemble=1_64_64_10.0_obs%s.h5' %ref_aux
     }
 
     # run simulation
@@ -146,6 +150,32 @@
     if generate_noda:
         rp.queue_run()
 
+
+    # Generate ensemble with one DA
+    ud = {
+        'aux' : 'oneda%s' %noda_sfx,
+        # Do blending for initial time-step
+        'initial_blending' : True
+    }
+
+    # For the data assimilation parameters,
+    # do not do DA at any time-point.
+    dap = {
+        'da_times' : [0.25],
+        'obs_attrs' : ['rho', 'rhou', 'rhov', 'rhoY', 'p2_nodes'],
+        'loc_setter' : (11, 11),
+        # Path to the generated observation
+        'obs_path' : path_to_obs + 'output_travelling_vortex/output_travelling_vortex_ensemble=1_64_64_10.0_obs%s.h5' %ref_aux
+    }
+
+    # run simulation
+    rp.ud = json.dumps(ud)
+    rp.dap = json.dumps(dap)
+    if generate_oneda:
+        rp.queue_run()
+
+
+    # Generate ensemble EnDA
     ud = {
         'aux' : '%s' %enda_sfx,
         # Do blending for initial time-step
@@ -162,15 +192,26 @@
     pi = p**0.4
     v = 4.0 / u_ref
 
+
+    sd_rho = 0.025
+    sd_rhou = 0.0325
+    sd_rhov = 0.0325
+    sd_rhoY = 0.00025
+    sd_pi = 0.0004
+    sds = np.array([sd_rho, sd_rhou, sd_rhov, sd_rhoY, sd_pi])
+    sds *= 2.0 # let's make it 10%
+    sds = np.sqrt(sds)
     dap = {
         'da_times' : da_times,
         # Assimilate all fields
         'obs_attrs' : ['rho', 'rhou', 'rhov', 'rhoY', 'p2_nodes'],
         # Path to the generated observation
-        'obs_path' : path_to_obs + 'output_travelling_vortex/output_travelling_vortex_ensemble=1_64_64_3.0_obs.h5',
-        'loc_setter' : (11, 11),
+        'obs_path' : path_to_obs + 'output_travelling_vortex/output_travelling_vortex_ensemble=1_64_64_10.0_obs%s.h5' %ref_aux,
+        'loc_setter' : (21, 21),
         # using variance: 1K for temperature, 4ms^-1 for velocity, 50Pa for pressure
         # 'sd_setter' : [rho**0.5, (rho*v)**0.5, (rho*v)**0.5, p**0.5, pi**0.5],
+        # 'sd_setter' :  [0.05**2, 0.05**2, 0.05**2, 0.1, 0.0004],
+        # 'sd_setter' :  list(sds),
     }
 
     # run simulation
@@ -222,7 +263,7 @@
         # Assimilate the momentum fields
         'obs_attrs' : ['rhou', 'rhov'],
         # Path to the generated observation
-        'obs_path' : path_to_obs + 'output_travelling_vortex/output_travelling_vortex_ensemble=1_64_64_3.0_obs.h5'
+        'obs_path' : path_to_obs + 'output_travelling_vortex/output_travelling_vortex_ensemble=1_64_64_3.0_obs%s.h5' %ref_aux
     }
 
     # run simulation
@@ -382,7 +423,7 @@
     # run simulation
     rp.ud = json.dumps(ud)
     rp.dap = json.dumps(dap)
-    # rp.queue_run()
+    rp.queue_run()
 
     ud = {
         'aux' : 'wda_7',
@@ -394,7 +435,7 @@
     # run simulation
     rp.ud = json.dumps(ud)
     rp.dap = json.dumps(dap)
-    # rp.queue_run()
+    rp.queue_run()