Skip to content
Snippets Groups Projects
Select Git revision
  • master default protected
1 result
Blame Permalink
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
solver.py 2.76 KiB
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 






import numpy as np
from pyTurbulence.syntheticTurbulence import compute_solenoidal_velocities, compute_solenoidal_pressure, compute_thermodynamic_fields
from pyTurbulence.spectrum import compute_tke_spectrum, plot_spectrum, energy_spectrum
import time
import os

def solve(params, results_dir):
    # Read the params dictionary
    nbModes = params["nbModes"]
    Mt = params["Mt"]
    mean_pressure = params["Pressure"]
    mean_temperature = params["Temperature"]
    mean_density = mean_pressure / mean_temperature
    gamma = params["gamma"]
    celerity = np.sqrt(gamma * mean_pressure / mean_density)
    Ek = 0.5 * (Mt * celerity) ** 2
    urms = np.sqrt(2. * Ek / 3.)
    k0 = params["k0"]
    domain_size = params["domain_size"]
    domain_resolution = params["domain_resolution"]
    xi = gamma * Mt ** 2
    spectrum = params["Spectrum"]

    # Generate the solenoidal velocity field
    start = time.time()
    u, v, w, wmax = compute_solenoidal_velocities(spectrum, nbModes, urms, k0, domain_size, domain_resolution, seed=params["seed"])
    TKE = np.mean(0.5 * (u.reshape(-1) ** 2 + v.reshape(-1) ** 2 + w.reshape(-1) ** 2))
    end = time.time()
    solenoidal_time = end - start

    # Generate the incompressible pressure fluctuations
    start = time.time()
    incompressible_pressure_fluctuations = compute_solenoidal_pressure(u, v, w, domain_resolution, domain_size)
    end = time.time()
    pressure_time = end - start

    # Compute the thermodynamic fields
    start = time.time()
    density, pressure, temperature = compute_thermodynamic_fields(mean_density, mean_pressure, mean_temperature, incompressible_pressure_fluctuations, gamma, Mt, params["case"])
    end = time.time()
    thermodynamic_time = end - start

    # Plot the spectrum
    start = time.time()
    knyquist, wave_numbers, tke_spectrum = compute_tke_spectrum(u, v, w, domain_size[0], domain_size[1], domain_size[2])
    real_spectrum = energy_spectrum(spectrum, wave_numbers, urms, k0)
    spectrum_path = os.path.join(results_dir, "spectrum.png")
    plot_spectrum(wmax, wave_numbers, tke_spectrum, real_spectrum=real_spectrum, save_path=spectrum_path)
    end = time.time()
    spectrum_time = end - start

    # Save the data
    start = time.time()
    u.tofile(os.path.join(results_dir, "u.dat"))
    v.tofile(os.path.join(results_dir, "v.dat"))
    w.tofile(os.path.join(results_dir, "w.dat"))
    pressure.tofile(os.path.join(results_dir, "p.dat"))
    temperature.tofile(os.path.join(results_dir, "T.dat"))
    end = time.time()
    save_time = end - start

    return {
        "TKE": TKE,
        "Ek": Ek,
        "solenoidal_time": solenoidal_time,
        "pressure_time": pressure_time,
        "thermodynamic_time": thermodynamic_time,
        "spectrum_time": spectrum_time,
        "save_time": save_time
    }