DISPATCH
legacy_reader.py
1 # Pythn 2/3 compatibility
2 from __future__ import print_function
3 import sys
4 if sys.version_info[0]!=2:
5  from importlib import reload
6 
7 import numpy as np
8 
9 """
10  Support functions for reading legacy format DISPATCH2 data.
11 """
12 
13 mesh_types = {1: 'Cartesian',
14  2: 'spherical',
15  3: 'cylindrical',}
16 mesh_types_r = {v: k for k, v in mesh_types.items()}
17 
18 def get_mesh_type(current_mesh_type):
19  """
20  A quick and hard-wired way to turn the integer mesh type index into a
21  string. The following must mirror the `mesh_types` variable in the
22  header of `mesh_mod.f90`.
23  """
24  assert isinstance(current_mesh_type, int), "expecting an integer for `mesh_type`"
25  try:
26  mesh_string = mesh_types[current_mesh_type]
27  except:
28  raise Exception("Unrecognised mesh type!")
29  return mesh_string
30 
31 def info (self):
32  print(" patch id:", self.id)
33  print(" no. of steps:", self.istep)
34  print(" coords.:", self.mesh_type)
35  print(" full patch:", self.m)
36  print(" unigrid patch:", self.n)
37  print(" total cells:", self.gn)
38  print(" lower indices:", self.li)
39  print(" upper indices:", self.ui)
40  print("no. of variables:", self.nvar)
41  print(" ghost cells:", self.nghost)
42  print(" time:", self.time)
43  print(" position:", self.pos)
44  print(" delta s:", self.ds)
45  print(" velocity:", self.vel)
46  print(" origin:", self.origin)
47  print(" l.l.c.:", self.llc)
48  print(" vol. centre:", self.centre)
49  print(" level:", self.level)
50  print(" quality:", self.quality)
51  print(" task type:", self.kind)
52 
53 def read_version0 (self, fp, line, verbose):
54  """
55  An concrete method for reading information for a single patch.
56  The following must mirror what is in `patch_mod.f90::output_legacy`.
57  """
58  self.kind = 'ramses'
59  self.eos = 'ideal'
60  self.opacity = 'none'
61  self.gamma = 1.4
62  # line 1: time, pos(1:3), ds(1:3), vel(1:3), origin(1:3), llc(1:3), cent(1:3)
63  self.time = float(line[0])
64  self.level = int (line[19])
65  self.quality = float(line[20])
66  self.pos = np.array(line[ 1: 4],dtype=np.float64)
67  self.ds = np.array(line[ 4: 7],dtype=np.float64)
68  self.vel = np.array(line[ 7:10],dtype=np.float64)
69  self.origin = np.array(line[10:13],dtype=np.float64)
70  self.llc = np.array(line[13:16],dtype=np.float64)
71  self.centre = np.array(line[16:19],dtype=np.float64)
72  # line 2: m(1:3), n(1:3), li(1:3), ui(1:3), gn(1:3), nghost(1:3), nvar
73  line = fp.readline().split()
74  self.m = np.array(line[ 0: 3],dtype=np.int32)
75  self.nwrite = (int(line[3]),int(line[4]), int(line[5]))
76  self.li = np.array(line[ 6: 9],dtype=np.int32)-1
77  self.ui = np.array(line[ 9:12],dtype=np.int32)-1
78  self.n = np.array(line[12:15],dtype=np.int32)
79  self.gn = np.array(line[15:18],dtype=np.int32)
80  self.nghost = np.array(line[18:21],dtype=np.int32)
81  self.nvar = int(line[21])
82  # line 3: id, iout, istep, mesh_type
83  line = fp.readline().split()
84  self.id = int(line[0])
85  self.iout = int(line[1])
86  self.istep = int(line[2])
87  self.dt = float(line[3])
88  self.mesh_type = get_mesh_type(int(line[4]))
89  if (verbose):
90  info (self)
91 
92 def read_version1 (self, fp, verbose):
93  """
94  An concrete method for reading information for a single patch.
95  The following must mirror what is in `patch_mod.f90::output_legacy`.
96  """
97  # lines 1-3: solver, eos, opacity
98  self.kind = fp.readline().strip()
99  self.eos = fp.readline().strip()
100  self.opacity = fp.readline().strip()
101  # line 4: time, dtime, pos(1:3), ds(1:3), vel(1:3), origin(1:3), llc(1:3), cent(1:3)
102  line = fp.readline().split()
103  self.time,self.dt = np.array(line[ 0: 2],dtype=np.float64)
104  self.pos = np.array(line[ 2: 5],dtype=np.float64)
105  self.ds = np.array(line[ 5: 8],dtype=np.float64)
106  self.vel = np.array(line[ 8:11],dtype=np.float64)
107  self.origin = np.array(line[11:14],dtype=np.float64)
108  self.llc = np.array(line[14:17],dtype=np.float64)
109  self.centre = np.array(line[17:20],dtype=np.float64)
110  self.quality,self.gamma = np.array(line[20:22],dtype=np.float64)
111  # line 5: m(1:3), n(1:3), li(1:3), ui(1:3), gn(1:3), nghost(1:3), nvar
112  line = fp.readline().split()
113  self.m = np.array(line[ 0: 3],dtype=np.int32)
114  self.nwrite = (int(line[3]),int(line[4]), int(line[5]))
115  self.li = np.array(line[ 6: 9],dtype=np.int32)-1
116  self.ui = np.array(line[ 9:12],dtype=np.int32)-1
117  self.n = np.array(line[12:15],dtype=np.int32)
118  self.gn = np.array(line[15:18],dtype=np.int32)
119  self.nghost = np.array(line[18:21],dtype=np.int32)
120  self.nvar = np.int (line[21 ])
121  self.level = np.int (line[22 ])
122  # line 3: id, iout, istep, mesh_type
123  line = fp.readline().split()
124  self.id = np.int(line[0])
125  self.iout = np.int(line[1])
126  self.istep = np.int(line[2])
127  self.mesh_type = get_mesh_type(int(line[3]))
128 
129 def single_patch (self, fp, fd=None, verbose=False):
130  if fd:
131  self.ioformat=3
132  else:
133  line=fp.readline()
134  chars=line[0:2]
135  if chars==' ':
136  self.ioformat=0
137  line=line.split()
138  read_version0 (self, fp, line, verbose)
139  else:
140  self.ioformat=np.int(chars)
141  read_version1 (self, fp, verbose)
142  derived (self, verbose)
143 
144 def derived(self, verbose=False):
145  self.dx = self.ds
146  self.corners = self._determine_corners()
147  self.active_corners = self._determine_corners(active_only=True)
148  self.corners = np.array(self.corners)
149  self.active_corners = np.array(self.active_corners)
150  self.size = self.ds*self.n
151  self.x=self.corners[0][0]+(0.5+np.arange(self.gn[0]))*self.ds[0]
152  self.y=self.corners[0][1]+(0.5+np.arange(self.gn[1]))*self.ds[1]
153  self.z=self.corners[0][2]+(0.5+np.arange(self.gn[2]))*self.ds[2]
154  if self.ioformat==2 or self.ioformat==4:
155  self.size = self.ds*(self.n-1)
156  self.active_corners[1] = self.active_corners[1] - self.ds
157  self.corners[1] = self.corners[1] - self.ds
158  self.x = self.x - 0.5*self.ds[0]
159  self.y = self.y - 0.5*self.ds[1]
160  self.z = self.z - 0.5*self.ds[2]
161  self.xs = self.x - 0.5*self.ds[0]
162  self.ys = self.y - 0.5*self.ds[1]
163  self.zs = self.z - 0.5*self.ds[2]
164  a = self.active_corners
165  self.extent=np.array([[a[0,1],a[1,1],a[0,2],a[1,2]],\
166  [a[0,0],a[1,0],a[0,2],a[1,2]],\
167  [a[0,0],a[1,0],a[0,1],a[1,1]]])
168  if (verbose):
169  info (self)