triple_decomp/phase_average.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
import sys
from pathlib import Path
import numpy as np
import pyvista as pv
import vtk
import h5py

class Dataset():
    """ Load BSL-specific data and common ops.
    """
    def __init__(self, folder):
        self.folder = Path(folder)

        file_glob_key = '*_curcyc_*up.h5'

        mesh_glob_key = '*h5'
        self.mesh_file = sorted(folder.glob(mesh_glob_key), key=lambda x: len(x.stem))[0]
        self.ts = '_ts='
        self.tssplit = '_'

        self.up_files = sorted(folder.glob(file_glob_key), key=self._get_ts)
        self.tsteps = len(sorted(folder.glob(file_glob_key), key=self._get_ts))
        self.times = sorted(folder.glob(file_glob_key), key=self._get_time)

    def __call__(self, idx=None, array='u', file=None):
        """ Return velocity in u_file. """
        if idx == None:
            h5_file = file
            with h5py.File(h5_file, 'r') as hf:
                val = np.array(hf[array])
        else:
            if array in ['u', 'p']:
                h5_file = self.up_files[idx]
                with h5py.File(h5_file, 'r') as hf:
                    val = np.array(hf['Solution'][array])
        return val

    def _get_ts(self, h5_file):
        """ Given a simulation h5_file, get ts. """
        return int(h5_file.stem.split(self.ts)[1].split(self.tssplit)[0])

    def _get_time(self, h5_file):
        """ Given a simulation h5_file, get time. """
        return float(h5_file.stem.split('_t=')[1].split('_')[0]) / 1000.0

    def assemble_mesh(self):
        """ Create UnstructuredGrid from h5 mesh file. """
        assert self.mesh_file.exists(), 'mesh_file does not exist.'

        with h5py.File(self.mesh_file, 'r') as hf:
            points = np.array(hf['Mesh']['coordinates'])*(10**-3)
            cells = np.array(hf['Mesh']['topology'])

            celltypes = np.empty(cells.shape[0], dtype=np.uint8)
            celltypes[:] = vtk.VTK_TETRA
            cell_type = np.ones((cells.shape[0], 1), dtype=int) * 4
            cells = np.concatenate([cell_type, cells], axis = 1)
            self.mesh = pv.UnstructuredGrid(cells.ravel(), celltypes, points)
            self.surf = self.mesh.extract_surface()
        return self

def phase_average(dd, idx, n_cycles, ts_pc, outfolder):
    u_c = dd(idx,array='u')
    u_pavg=np.zeros(u_c.shape)
    n_cycles_c = n_cycles
    for n in range(n_cycles):
        ndx = n*ts_pc+idx
        u_n = dd(ndx,array='u')
        if (n>0) and (np.linalg.norm(u_c-u_n)/np.linalg.norm(u_c)/u_c.shape[0]>0.1):
            print("Velocity greater than 10 percent different, skipping tstep={}({}), cycle={}".format(idx,ndx,n))
            n_cycles_c -= 1
        else:
            u_pavg += u_n
    u_pavg /=n_cycles_c
    f = h5py.File(outfolder + '/{}.h5'.format(idx), 'w')
    f.create_dataset('u', data = u_pavg)

if __name__ == "__main__":
    results_folder = Path(sys.argv[1])
    n_cycles = int(sys.argv[2])
    ts_pc = int(sys.argv[3])
    idx = int(sys.argv[4])
    outfolder = './triple_decomp_data/phase_avg_{}_cycles'.format(n_cycles)
    if not Path(outfolder).exists():
        Path(outfolder).mkdir(parents=True, exist_ok=True)
    if not Path(outfolder + '/{}.h5'.format(idx)).exists():
        dd = Dataset(results_folder)
        main_folder = Path(results_folder).parents[1]
        dd=dd.assemble_mesh()
        phase_average(dd,idx,n_cycles, ts_pc, outfolder)