-
Notifications
You must be signed in to change notification settings - Fork 24
/
Copy pathpostProc.py
100 lines (82 loc) · 4.21 KB
/
postProc.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
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
import numpy as np
import json
#import argparse
import matplotlib.pyplot as plt
import matplotlib.ticker
from matplotlib.legend_handler import HandlerTuple
import sys
import os
import subprocess
import h5py
import pandas
# load HTR modules
sys.path.insert(0, os.path.expandvars("$HTR_DIR/scripts/modules"))
import MulticomponentMix
import HTRrestart
dir_name = os.path.join(os.environ['HTR_DIR'], 'testcases/Speelman')
input_file = os.path.join(dir_name, 'Speelman.json')
ref_file = os.path.join(dir_name, 'Ref.dat')
##############################################################################
# Read HTR input file #
##############################################################################
with open(input_file) as f:
data = json.load(f)
xNum = data["Grid"]["xNum"]
yNum = data["Grid"]["yNum"]
zNum = data["Grid"]["zNum"]
mix = MulticomponentMix.Mix(data["Flow"]["mixture"])
##############################################################################
# Read reference solution #
##############################################################################
Ref = pandas.read_csv(ref_file, delim_whitespace=True, encoding= "unicode_escape")
##############################################################################
# Process result file #
##############################################################################
def process(nstep):
##############################################################################
# Read HTR output data #
##############################################################################
restart = HTRrestart.HTRrestart(data)
restart.attach(sampleDir=os.path.join(dir_name, 'sample0'), step=nstep)
# Get the data
x = restart.load("centerCoordinates")[0,:,0,1]
pressure = restart.load("pressure" )[0,:,0]
temperature = restart.load("temperature" )[0,:,0]
density = restart.load("rho" )[0,:,0]
velocity = restart.load("velocity" )[0,:,0,:]
molarFracs = restart.load("MolarFracs" )[0,:,0,:]
return x*mix.LRef, pressure*mix.PRef, temperature*mix.TRef, density*mix.rhoRef, velocity*mix.uRef, molarFracs, f.attrs.get("SpeciesNames")
##############################################################################
# Plot #
##############################################################################
x, pressure, temperature, density, velocity, molarFracs, specieNames = process(3700000)
plt.rcParams.update({'font.size': 16})
plt.figure(1)
plt.plot(x, temperature, '-k', label="HTR solver")
plt.plot(Ref["y"], Ref["T"], 'ok', label="Cantera", markevery=5e-2)
plt.xlabel(r'$x [m]$', fontsize = 20)
plt.ylabel(r'$T [K]$', fontsize = 20)
plt.ticklabel_format(axis="x", style="sci", scilimits=(0,0), useMathText=True)
plt.gca().set_xlim(0, 5e-3)
plt.gca().set_ylim(350, 2400)
plt.legend()
plt.savefig('Temperature.eps', bbox_inches='tight')
plt.figure(2)
CH4my, = plt.semilogy(x, molarFracs[:, mix.FindSpecies("CH4", specieNames)], '-k')
CO2my, = plt.semilogy(x, molarFracs[:, mix.FindSpecies("CO2", specieNames)], '-r')
CHmy, = plt.semilogy(x, molarFracs[:, mix.FindSpecies("CH", specieNames)], '-b')
H2my, = plt.semilogy(x, molarFracs[:, mix.FindSpecies("H2", specieNames)], '-g')
CH4ca, = plt.semilogy(Ref["y"], Ref["X_CH4"], 'ok', markevery=0.05)
CO2ca, = plt.semilogy(Ref["y"], Ref["X_CO2"], 'or', markevery=0.05)
CHca, = plt.semilogy(Ref["y"], Ref["X_CH"], 'ob', markevery=0.05)
H2ca, = plt.semilogy(Ref["y"], Ref["X_H2"], 'og', markevery=0.05)
plt.xlabel(r'$x [m]$', fontsize = 20)
plt.ylabel(r'$X_i$', fontsize = 20)
plt.ticklabel_format(axis="x", style="sci", scilimits=(0,0), useMathText=True)
plt.gca().set_xlim(0, 5e-3)
plt.gca().set_ylim(1e-8, 0.5)
plt.legend([(CH4my, CH4ca), (CO2my, CO2ca), (CHmy, CHca), (H2my, H2ca),],
[r"$X_{CH_4}$", r"$X_{CO_2}$", r"$X_{CH}$", r"$X_{H2}$"],
handler_map={tuple: HandlerTuple(ndivide=None)}, handlelength=8.0)
plt.savefig('MolarFractions.eps', bbox_inches='tight')
plt.show()