plotSolution.py 4.14 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
################################################################################
# This file is part of SWIFT.
# Copyright (c) 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
################################################################################

##
# This script plots the Disc-Patch_*.hdf5 snapshots.
# It takes two (optional) parameters: the counter value of the first and last
23
# snapshot to plot (default: 0 21).
24
25
26
27
28
29
30
31
32
33
34
35
36
##

import numpy as np
import h5py
import matplotlib
matplotlib.use("Agg")
import pylab as pl
import glob
import sys

# Parameters
surface_density = 10.
scale_height = 100.
37
z_disc = 400.
38
39
z_trunc = 300.
z_max = 350.
40
utherm = 20.2678457288
41
42
43
gamma = 5. / 3.

start = 0
44
stop = 21
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
if len(sys.argv) > 1:
  start = int(sys.argv[1])
if len(sys.argv) > 2:
  stop = int(sys.argv[2])

# Get the analytic solution for the density
def get_analytic_density(x):
  return 0.5 * surface_density / scale_height / \
           np.cosh( (x - z_disc) / scale_height )**2

# Get the analytic solution for the (isothermal) pressure
def get_analytic_pressure(x):
  return (gamma - 1.) * utherm * get_analytic_density(x)

# Get the data fields to plot from the snapshot file with the given name:
#  snapshot time, z-coord, density, pressure, velocity norm
def get_data(name):
  file = h5py.File(name, "r")
  coords = np.array(file["/PartType0/Coordinates"])
  rho = np.array(file["/PartType0/Density"])
  u = np.array(file["/PartType0/InternalEnergy"])
  v = np.array(file["/PartType0/Velocities"])

  P = (gamma - 1.) * rho * u

  vtot = np.sqrt( v[:,0]**2 + v[:,1]**2 + v[:,2]**2 )

  return float(file["/Header"].attrs["Time"]), coords[:,2], rho, P, vtot

# scan the folder for snapshot files and plot all of them (within the requested
# range)
for f in sorted(glob.glob("Disc-Patch_*.hdf5")):
  num = int(f[-8:-5])
  if num < start or num > stop:
    continue

  print "processing", f, "..."

83
  zrange = np.linspace(0., 2. * z_disc, 1000)
84
85
86
87
88
89
  time, z, rho, P, v = get_data(f)

  fig, ax = pl.subplots(3, 1, sharex = True)

  ax[0].plot(z, rho, "r.")
  ax[0].plot(zrange, get_analytic_density(zrange), "k-")
90
91
92
93
94
  ax[0].plot([z_disc - z_max, z_disc - z_max], [0, 10], "k--", alpha=0.5)
  ax[0].plot([z_disc + z_max, z_disc + z_max], [0, 10], "k--", alpha=0.5)
  ax[0].plot([z_disc - z_trunc, z_disc - z_trunc], [0, 10], "k--", alpha=0.5)
  ax[0].plot([z_disc + z_trunc, z_disc + z_trunc], [0, 10], "k--", alpha=0.5)
  ax[0].set_ylim(0., 1.2 * get_analytic_density(z_disc))
95
96
97
98
  ax[0].set_ylabel("density")

  ax[1].plot(z, v, "r.")
  ax[1].plot(zrange, np.zeros(len(zrange)), "k-")
99
100
101
102
103
  ax[1].plot([z_disc - z_max, z_disc - z_max], [0, 10], "k--", alpha=0.5)
  ax[1].plot([z_disc + z_max, z_disc + z_max], [0, 10], "k--", alpha=0.5)
  ax[1].plot([z_disc - z_trunc, z_disc - z_trunc], [0, 10], "k--", alpha=0.5)
  ax[1].plot([z_disc + z_trunc, z_disc + z_trunc], [0, 10], "k--", alpha=0.5)
  ax[1].set_ylim(-0.5, 10.)
104
105
106
107
  ax[1].set_ylabel("velocity norm")

  ax[2].plot(z, P, "r.")
  ax[2].plot(zrange, get_analytic_pressure(zrange), "k-")
108
109
110
111
112
113
  ax[2].plot([z_disc - z_max, z_disc - z_max], [0, 10], "k--", alpha=0.5)
  ax[2].plot([z_disc + z_max, z_disc + z_max], [0, 10], "k--", alpha=0.5)
  ax[2].plot([z_disc - z_trunc, z_disc - z_trunc], [0, 10], "k--", alpha=0.5)
  ax[2].plot([z_disc + z_trunc, z_disc + z_trunc], [0, 10], "k--", alpha=0.5)
  ax[2].set_xlim(0., 2. * z_disc)
  ax[2].set_ylim(0., 1.2 * get_analytic_pressure(z_disc))
114
115
116
117
118
119
120
  ax[2].set_xlabel("z")
  ax[2].set_ylabel("pressure")

  pl.suptitle("t = {0:.2f}".format(time))

  pl.savefig("{name}.png".format(name = f[:-5]))
  pl.close()