From 7108b952a8fb7bc64c193a924d09db29f592b86d Mon Sep 17 00:00:00 2001
From: Philipp Griewank <philipp.griewank@uni-koeln.de>
Date: Tue, 25 May 2021 09:14:54 +0200
Subject: [PATCH] Added sine wave initial conditions
Also fixed a bug where the truth was one time step out of sync.
---
da_functions.py | 25 +++++++++++++++++++++----
model_functions.py | 6 ++++--
2 files changed, 25 insertions(+), 6 deletions(-)
diff --git a/da_functions.py b/da_functions.py
index 06e0e50..0b3d725 100644
--- a/da_functions.py
+++ b/da_functions.py
@@ -7,7 +7,7 @@ from model_functions import *
def set_da_constants(ncyc=10,nt=1,dt=500,u_std=0.5,dhdt_std=1e-4,True_std_obs=0.1,used_std_obs=0.1,pert_std_obs=0.00,obs_loc_h=np.arange(5,101,15),nens=100,nexp=1,init_noise=0.,fixed_seed=True,
loc=None,init_spread = False, init_spread_h=0.5,init_spread_x = 500.,
- loc_length = 1000,loc_type='gaussian'):
+ loc_length = 1000,loc_type='gaspari_cohn'):
"""
Sets constants needed for data assimilation and stores them in a dictionary.
There is some confusting misnaming going on, e.g. "u_std_ens" = u_std, but nothing that needs to be changed immediately
@@ -56,7 +56,21 @@ def set_da_constants(ncyc=10,nt=1,dt=500,u_std=0.5,dhdt_std=1e-4,True_std_obs=0.
DA_const["loc_type"] = loc_type # localization yes or no, no for now
return DA_const
+
+
+def sine_initial_condition(x,n):
+ """
+ Generates a sin curve on the x grid with overlapping waves defined by n.
+ """
+ mu = 0
+ y = np.zeros_like(x)
+ for ni in n:
+ y = y+ np.sin(x/np.max(x)*2*np.pi*ni)
+ return y
+
+
+
def create_states_dict(j,states,m_const,da_const):
"""
Generates the initial analysis and truth.
@@ -79,7 +93,8 @@ def create_states_dict(j,states,m_const,da_const):
nx = m_const["nx"]
#initial conditions
- h_initial = gaussian_initial_condition(m_const["x_grid"],m_const["h_init_std"])
+ if m_const['init_func']=='gaus': h_initial = gaussian_initial_condition(m_const["x_grid"],m_const["h_init_std"])
+ if m_const['init_func']=='sine': h_initial = sine_initial_condition( m_const["x_grid"],m_const["sine_init"])
#Generate truth
@@ -139,12 +154,14 @@ def generate_obs(truth,states,m_const,da_const,j,t):
u_truth = np.random.normal(m_const["u_ref"],m_const["u_std_truth"])
dhdt_truth = np.random.normal(m_const["dhdt_ref"],m_const["dhdt_std_truth"])
if da_const["fixed_seed"]==True:
- truth = linear_advection_model(states['truth'][0],u_truth,dhdt_truth,m_const["dx"],da_const["dt"]*t,da_const["nt"])
+ #t starts at zero, that is why we need a plus 1
+ truth = linear_advection_model(states['truth'][0],u_truth,dhdt_truth,m_const["dx"],da_const["dt"]*(t+1),da_const["nt"])
+
else:
truth = linear_advection_model(truth,u_truth,dhdt_truth,m_const["dx"],da_const["dt"],da_const["nt"])
#make obs by adding some noise, a fixed seed should change over time so that the differences are not always the same for each measurement location
- if da_const["fixed_seed"]==True: np.random.seed((j+1)*t)
+ if da_const["fixed_seed"]==True: np.random.seed((j+1)*(t+1))
obs = truth + np.random.normal(0,da_const["True_std_obs"],m_const["nx"])
states["truth"] = states["truth"]+[truth]
states["obs"] = states["obs"]+[obs]
diff --git a/model_functions.py b/model_functions.py
index d3b7c2b..ad05856 100644
--- a/model_functions.py
+++ b/model_functions.py
@@ -3,7 +3,7 @@
import numpy as np
-def set_model_constants(nx=101,dx=100,u=2.33333333,u_std=0,dhdt=0,dhdt_std=0,h_init_std = 500):
+def set_model_constants(nx=101,dx=100,u=2.33333333,u_std=0,dhdt=0,dhdt_std=0,h_init_std = 500, init_func='gaus',sine_init=[1]):
"""
Sets all the constants used to describe the model.
Currently has only gaussian initial conditions.
@@ -18,7 +18,9 @@ def set_model_constants(nx=101,dx=100,u=2.33333333,u_std=0,dhdt=0,dhdt_std=0,h_i
const["x_grid"] = np.arange(nx)*dx
#properties of truth wave
- const["h_init_std"] = h_init_std # width of the initial gaussian distribution
+ const["h_init_std"] = h_init_std # width of the initial gaussian distribution, or number of sine waves
+ const["sine_init"] = sine_init # width of the initial gaussian distribution, or number of sine waves
+ const["init_func"] = init_func # currently a gaussian ('gaus') or a sine wave ('sine') are available
const["u_ref"] = u # reference u around which ensemble is generated, truth can
const["u_std_truth"] = u_std # standard deviation of u for truth
const["dhdt_ref"] = dhdt # reference dhdt around which ensemble is generated
--
GitLab