Castro/Exec/hydro_tests/Sod_stellar/testsuite_analysis/test1-helm.py at eeb529b86c250bd7fa27dd73b757eaabdb696aaa · zingale/Castro · 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
124
125
126
#!/usr/bin/env python

# run as: ./test1-helm.py castro_exec_dir plotfle

# note: this relies on fextract.XXXX.ex being in your path somewhere

import sys
import os
import shutil
import numpy as np

import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt

def process(castro_exec_dir, plotfile):

    run_dir = os.getcwd()

    # 1. find the fextract tool by looking through the User's path
    path = os.environ["PATH"].split(":")

    for d in path:
        full_dir = os.path.expanduser(d)
        if not os.path.isdir(full_dir):
            continue
        for f in os.listdir(full_dir):
            if (os.path.isfile(full_dir+"/"+f) and
                f.startswith("fextract") and f.endswith(".ex")):
                analysis_routine = full_dir+"/"+f
                break

    print(f"analysis_routine = {analysis_routine}")

    shutil.copy(analysis_routine, run_dir)


    # 2. analyze the data


    # output the average profile
    os.system("./{} -s {} {}".format(os.path.basename(analysis_routine), "test1-helm.out", plotfile))


    analytic = castro_exec_dir + "Exec/hydro_tests/Sod_stellar/Verification/test1.exact.128.out"
    analytic_data = np.loadtxt(analytic)

    # need to be more flexible with the data from the simulations, as the
    # columns shift depending on the dimensionality.  This gets the column
    # names from the header
    data = np.genfromtxt("test1-helm.out", skip_header=2, names=True)

    # 3. make the plot
    plt.subplot(411)

    plt.plot(analytic_data[:,1], analytic_data[:,2])
    plt.scatter(data["x"], data["density"], marker="+", color="r")

    plt.xlabel("x")
    plt.ylabel("density")
    plt.xlim(0,1.e6)


    plt.subplot(412)

    # figure out which dimensions are present
    d = 1
    for n in data.dtype.names:
        if n == "ymom" or n == "zmom": d += 1

    dim = "xmom"
    if d >= 2 and np.ptp(data["ymom"]) > np.ptp(data["xmom"]):
        dim = "ymom"
    if d == 3 and np.ptp(data["zmom"]) > np.ptp(data[dim]):
        dim = "zmom"

    plt.plot(analytic_data[:,1], analytic_data[:,3])
    plt.scatter(data["x"], data[dim]/data["density"], marker="+", color="r")

    plt.xlabel("x")
    plt.ylabel("velocity")
    plt.xlim(0,1.e6)


    plt.subplot(413)

    plt.plot(analytic_data[:,1], analytic_data[:,4])
    plt.scatter(data["x"], data["pressure"], marker="+", color="r")

    plt.xlabel("x")
    plt.ylabel("pressure")
    plt.xlim(0,1.e6)


    plt.subplot(414)

    plt.plot(analytic_data[:,1], analytic_data[:,5])
    plt.scatter(data["x"], data["Temp"], marker="+", color="r")

    plt.xlabel("x")
    plt.ylabel("temperature")
    plt.xlim(0,1.e6)


    ax = plt.gca()
    ax.set_yscale("log")

    f = plt.gcf()
    f.set_size_inches(6.0, 9.0)

    plt.tight_layout()

    index = plotfile.rfind("_plt")
    if (index > 0):
        plt.savefig(plotfile[:index] + ".png")
    else:
        plt.savefig("test1-helm.png")


if __name__ == "__main__":

    castro_exec_dir = str(sys.argv[1])
    plotfile = str(sys.argv[2])

    process(castro_exec_dir, plotfile)