triple_decomp/stats.py at main · Biomedical-Simulation-Lab/triple_decomp · GitHub

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
119
120
121
122
123
import sys
from pathlib import Path
import numpy as np
import pyvista as pv
import vtk
import h5py
import phase_average as PA

def time_avg(dd,n_cycles,ts_pc, outfolder):
    u_A=np.zeros((dd(0, array='u').shape))
    for iidx in range(ts_pc):
        f = h5py.File(str(outfolder) + '/phase_avg_{}_cycles/{}.h5'.format(n_cycles, iidx), 'r')
        u_A += np.array(f['u'])/ts_pc
    f_ii = h5py.File(str(outfolder) + 'phase_avg_{}_cycles/time_avg.h5'.format(n_cycles), 'w')
    f_ii.create_dataset('u_A', data = u_A)

def stats(dd, idx, n_cycles, ts_pc, outfolder, nu):
    f = h5py.File(str(outfolder) + '/phase_avg_{}_cycles/{}.h5'.format(n_cycles, idx), 'r')
    f_ii = h5py.File(str(outfolder) + 'phase_avg_{}_cycles/time_avg.h5'.format(n_cycles), 'r')
    #get fluctuations of the first cycle
    dd.mesh.point_data['u_p1'] = dd(idx, array='u')-np.array(f['u'])
    u_p1 = dd.mesh.point_data['u_p1']
    u_pp1 = dd(idx, array='u')-np.array(f_ii['u_A'])
    Re_sp1 = np.zeros((len(u_p1),3,3))
    Re_spp1 = np.zeros((len(u_pp1),3,3))
    for i in range(3):
        for j in range(3):
            Re_sp1[:,i,j] = np.array((u_p1[:,i]*u_p1[:,j]))
            Re_spp1[:,i,j] = np.array((u_pp1[:,i]*u_pp1[:,j]))
    TKE_p1 = np.sum(u_p1**2, axis=1)
    TKE_pp1 = np.sum(u_pp1**2, axis=1)
    grad = dd.mesh.compute_derivative(scalars="u_p1", gradient=True, qcriterion=False, faster=False)
    J_1 = grad.point_data['gradient'].reshape(-1, 3, 3)
    s_p1 = 0.5*(J_1+np.transpose(J_1, axes=(0,2,1)))
    eps_p1 = 2*nu*np.sum(np.sum(s_p1*s_p1, axis=2), axis=1)

    #get fluctuations of the last cycle
    dd.mesh.point_data['u_pL'] = dd(idx+n_cycles*ts_pc, array='u')-np.array(f['u'])
    u_pL = dd.mesh.point_data['u_pL']
    u_ppL = dd(idx+n_cycles*ts_pc, array='u')-np.array(f_ii['u_A'])
    Re_spL = np.zeros((len(u_pL),3,3))
    Re_sppL = np.zeros((len(u_ppL),3,3))
    for i in range(3):
        for j in range(3):
            Re_spL[:,i,j] = np.array((u_pL[:,i]*u_pL[:,j]))
            Re_sppL[:,i,j] = np.array((u_ppL[:,i]*u_ppL[:,j]))
    TKE_pL = np.sum(u_pL**2, axis=1)
    TKE_ppL = np.sum(u_ppL**2, axis=1)
    grad = dd.mesh.compute_derivative(scalars="u_pL", gradient=True, qcriterion=False, faster=False)
    J_l = grad.point_data['gradient'].reshape(-1, 3, 3)
    s_l = 0.5*(J_l+np.transpose(J_l, axes=(0,2,1)))
    eps_pL = 2*nu*np.sum(np.sum(s_l*s_l, axis=2), axis=1)

    #get average fluctuations for all cycles
    dd.mesh.point_data['u_p_avg']=np.zeros((dd(idx, array='u').shape))
    u_p = np.zeros((dd.mesh.point_data['u_p_avg'].shape))
    u_pp = np.zeros((dd.mesh.point_data['u_p_avg'].shape))
    Re_s = np.zeros((len(dd.mesh.point_data['u_p_avg']),3,3))
    for n in range(n_cycles):
        u_p += (dd(n*ts_pc+idx, array='u')-np.array(f['u']))
        u_pp += (dd(n*ts_pc+idx, array='u')-np.array(f_ii['u_A']))

    Re_sp = np.zeros((len(dd.mesh.point_data['u_p_avg']),3,3))
    Re_spp = np.zeros((len(dd.mesh.point_data['u_p_avg']),3,3))
    for i in range(3):
        for j in range(3):
            Re_sp[:,i,j] = np.array((u_p[:,i]*u_p[:,j]))
            Re_spp[:,i,j] = np.array((u_pp[:,i]*u_pp[:,j]))
    TKE_p = np.sum(u_p**2, axis=1)
    TKE_pp = np.sum(u_pp**2, axis=1)
    #compute average dissipation rate
    dd.mesh.point_data['u_p_avg'] = u_p
    grad = dd.mesh.compute_derivative(scalars="u_p_avg", gradient=True, qcriterion=False, faster=False)
    J = grad.point_data['gradient'].reshape(-1, 3, 3)
    s_p = 0.5*(J+np.transpose(J, axes=(0,2,1)))
    eps = 2*nu*np.sum(np.sum(s_p*s_p, axis=2), axis=1)
    #get averages
    u_p /= n_cycles
    Re_s /= n_cycles
    eps /= n_cycles
    f_f = h5py.File(str(outfolder) + 'phase_avg_{}_cycles/stats/stats_{}.h5'.format(n_cycles,idx), 'w')
    #average data
    f_f.create_dataset('u_pAvg', data = u_p)
    f_f.create_dataset('Re_sp', data = Re_sp)
    f_f.create_dataset('Re_spp', data = Re_spp)
    f_f.create_dataset('TKE_p', data = TKE_p)
    f_f.create_dataset('TKE_pp', data = TKE_pp)
    f_f.create_dataset('eps', data = eps)
    #first cycle data
    f_f.create_dataset('u_p1', data = u_p1)
    f_f.create_dataset('u_pp1', data = u_pp1)
    f_f.create_dataset('Re_sp1', data = Re_sp1)
    f_f.create_dataset('Re_spp1', data = Re_spp1)
    f_f.create_dataset('TKE_p1', data = TKE_p1)
    f_f.create_dataset('TKE_pp1', data = TKE_pp1)
    f_f.create_dataset('eps_1', data = eps_p1)
    #last cycle data
    f_f.create_dataset('u_pL', data = u_pL)
    f_f.create_dataset('u_ppL', data = u_ppL)
    f_f.create_dataset('Re_spL', data = Re_spL)
    f_f.create_dataset('Re_sppL', data = Re_sppL)
    f_f.create_dataset('TKE_pL', data = TKE_pL)
    f_f.create_dataset('TKE_ppL', data = TKE_ppL)
    f_f.create_dataset('eps_L', data = eps_pL)

if __name__ == "__main__":
    results_folder = Path(sys.argv[1])
    n_cycles = int(sys.argv[2])
    ts_pc = int(sys.argv[3])
    outfolder = 'triple_decomp_data/'
    nu = 0.0000035
    if not Path(outfolder+'/phase_avg_{}_cycles/stats/'.format(n_cycles)).exists():
        Path(outfolder+'/phase_avg_{}_cycles/stats/'.format(n_cycles)).mkdir(parents=True, exist_ok=True)
    if len(sys.argv)>4:
        idx = int(sys.argv[4])
        if not Path(outfolder + 'phase_avg_{}_cycles/stats/stats_{}.h5'.format(n_cycles,idx)).exists():
            dd = PA.Dataset(results_folder)
            dd=dd.assemble_mesh()
            stats(dd, idx, n_cycles, ts_pc, outfolder, nu)
            print('Completed index {}, {} cycles'.format(idx, n_cycles))
    else:
        dd = PA.Dataset(results_folder)
        time_avg(dd, n_cycles, ts_pc, outfolder)