-
Notifications
You must be signed in to change notification settings - Fork 0
/
initial_conditions.py
74 lines (45 loc) · 3.32 KB
/
initial_conditions.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
import numpy as np
from filters import Filters, Grids
from functions_for_solver import compute_VgradW_hat_np, get_psi_hat
import h5py
def get_initial_conditions(from_restart_file:bool, filters:Filters, grid:Grids, fname:str="none", USE_HDF5:bool=False, filter_restart_state=False, transpose = False):
if from_restart_file:
return initial_conditions_from_restart( fname, filters, grid, USE_HDF5, filter_restart_state, transpose)
else:
return create_initial_conditions(grid, filters)
def initial_conditions_from_restart( fname, filters:Filters, grid:Grids, USE_HDF5, filter_restart_state=False, transpose = False):
##################### initial conditions ################################
#fname = 'C:/Users/rik/Documents/Python Scripts/CAF_reduced_SGS-master/source/ref_data/restart_t_' + str(np.around(t_by_day,1))+file_extension
# define high fidelity arrays
VgradW_hat_nm1_HF = np.zeros((grid.N_HF,grid.N_HF),dtype=np.complex128)
VgradW_hat_nm10_HF = np.zeros((grid.N_HF,grid.N_HF),dtype=np.complex128)
w_hat_n_HF = np.zeros((grid.N_HF,grid.N_HF),dtype=np.complex128)
w_hat_nm1_HF = np.zeros((grid.N_HF,grid.N_HF),dtype=np.complex128)
w_hat_nm10_HF = np.zeros((grid.N_HF,grid.N_HF),dtype=np.complex128)
#create HDF5 file
h5f = h5py.File(fname, 'r')
VgradW_hat_nm1_HF[filters.P_HF==1]= np.array(h5f['VgradW_hat_nm1_HF']).flatten()*(filters.N_HF/filters.N_HF_resolved)**2
VgradW_hat_nm10_HF[filters.P_HF==1]= np.array(h5f['VgradW_hat_nm10_HF']).flatten()*(filters.N_HF/filters.N_HF_resolved)**2
w_hat_n_HF[filters.P_HF==1]= np.array(h5f['w_hat_n_HF']).flatten()*(filters.N_HF/filters.N_HF_resolved)**2
w_hat_nm1_HF[filters.P_HF==1]= np.array(h5f['w_hat_nm1_HF']).flatten()*(filters.N_HF/filters.N_HF_resolved)**2
w_hat_nm10_HF[filters.P_HF==1]= np.array(h5f['w_hat_nm10_HF']).flatten()*(filters.N_HF/filters.N_HF_resolved)**2
h5f.close()
psi_hat_n_HF = get_psi_hat(w_hat_n_HF, grid.k_squared_nonzero_HF)
psi_hat_nm1_HF = get_psi_hat(w_hat_nm1_HF, grid.k_squared_nonzero_HF)
return w_hat_n_HF, w_hat_nm1_HF, VgradW_hat_nm1_HF, psi_hat_n_HF, psi_hat_nm1_HF
def create_initial_conditions(grid:Grids, filters:Filters):
#initial condition
w_HF = np.sin(4.0*grid.x_HF)*np.sin(4.0*grid.y_HF) + 0.4*np.cos(3.0*grid.x_HF)*np.cos(3.0*grid.y_HF) + \
0.3*np.cos(5.0*grid.x_HF)*np.cos(5.0*grid.y_HF) + 0.02*np.sin(grid.x_HF) + 0.02*np.cos(grid.y_HF)
w_hat_n_HF = filters.P_HF*np.fft.fft2(w_HF)
psi_hat_n_HF = get_psi_hat(w_hat_n_HF, grid.k_squared_nonzero_HF)
#initial Fourier coefficients at time n-1
w_hat_nm1_HF = np.copy(w_hat_n_HF)
#initial Fourier coefficients of the jacobian at time n and n-1
VgradW_hat_n_HF = compute_VgradW_hat_np(w_hat_n_HF, psi_hat_n_HF, filters.P_HF, grid.k_x_HF, grid.k_y_HF)
VgradW_hat_nm1_HF = np.copy(VgradW_hat_n_HF)
psi_hat_nm1_HF = get_psi_hat(w_hat_nm1_HF, grid.k_squared_nonzero_HF)
return w_hat_n_HF, w_hat_nm1_HF, VgradW_hat_nm1_HF, psi_hat_n_HF, psi_hat_nm1_HF
#V_i for Q_i = (1, omega)
#V_hat_w1 = P_LF*np.fft.fft2(np.ones([N_LF, N_LF]))
####################################################################################################################