DISPATCH
_dispatch_grid.py
1 # -*- coding: utf-8 -*-
2 
3 """
4  Support functions for reading or calculating DISPATCH2 grid and geometry information.
5 
6 """
7 
8 import numpy as np
9 
10 class GeometricFactors(dict):
11  """Calculate and store the geometric factors used by curvilinear grids."""
12 
13  def __init__(self, patch):
14  """Constructor."""
15 
16  # Define geometric factors with the same notation as in `mesh_mod` ("c"
17  # for zone-centred and "f" for face-centred).
18  self['h2c'] = None
19  self['h2f'] = None
20  self['h31c'] = None
21  self['h31f'] = None
22  self['h32c'] = None
23  self['h32f'] = None
24  self['dx1c'] = None
25  self['dx1f'] = None
26  self['dx2c'] = None
27  self['dx2f'] = None
28  self['dx3c'] = None
29  self['dx3f'] = None
30  self['dvol1c'] = None
31  self['dvol1f'] = None
32  self['dvol2c'] = None
33  self['dvol2f'] = None
34  self['dvol3c'] = None
35  self['dvol3f'] = None
36  self['dar1c'] = None
37  self['dar1f'] = None
38  self['dar2c'] = None
39  self['dar2f'] = None
40  self['dar31c'] = None
41  self['dar31f'] = None
42  self['dar32c'] = None
43  self['dar32f'] = None
44  # initialize the grid
45  self.init_grid(patch)
46 
47  def init_grid(self, p):
48  """Initialise geometric factors based on coord. type."""
49 
50  if p.mesh_type == 'Cartesian': self.init_Cartesian(p)
51  elif p.mesh_type == 'cylindrical': self.init_cylindrical(p)
52  elif p.mesh_type == 'spherical': self.init_spherical(p)
53 
54  def init_Cartesian(self, p):
55  """Initialise geometric factors for a Cartesian coord. system."""
56 
57  n1, n2, n3 = p.ncell
58 
59  # 1-direction
60  self['h2c'] = np.ones(n1)
61  self['h2f'] = np.ones(n1)
62  self['h31c'] = self['h2c'].view()
63  self['h31f'] = self['h2f'].view()
64 
65  # 2-direction
66  self['h32c'] = np.ones(n2)
67  self['h32f'] = self['h32c'].view()
68 
69  # linear size elements
70  self['dx1c'] = np.ones(n1) * p.ds[0]
71  self['dx1f'] = np.ones(n1) * p.ds[0]
72  self['dx2c'] = np.ones(n2) * p.ds[1]
73  self['dx2f'] = np.ones(n2) * p.ds[1]
74  self['dx3c'] = np.ones(n3) * p.ds[2]
75  self['dx3f'] = np.ones(n3) * p.ds[2]
76 
77  # volume elements
78  self['dvol1c'] = np.ones(n1) * p.ds[0]
79  self['dvol1f'] = np.ones(n1) * p.ds[0]
80  self['dvol2c'] = np.ones(n2) * p.ds[1]
81  self['dvol2f'] = np.ones(n2) * p.ds[1]
82  self['dvol3c'] = np.ones(n3) * p.ds[2]
83  self['dvol3f'] = np.ones(n3) * p.ds[2]
84 
85  # area elements
86  self['dar1c'] = self['h31c'] * self['h2c']
87  self['dar1f'] = self['h31f'] * self['h2f']
88  self['dar2c'] = self['h31f'] * p.ds[0] / self['dvol1c']
89  self['dar2f'] = self['h31c'] * p.ds[0] / self['dvol1f']
90  self['dar31c'] = self['h2f'] * p.ds[0] / self['dvol1c']
91  self['dar31f'] = self['h2c'] * p.ds[0] / self['dvol1f']
92  self['dar32c'] = p.ds[1] / self['dvol2c']
93  self['dar32f'] = p.ds[1] / self['dvol2f']
94 
95  def init_cylindrical(self, p):
96  """Initialise geometric factors for a cylindrical coord. system."""
97 
98  n1, n2, n3 = p.ncell
99 
100  # 1-direction
101  self['h2c'] = np.ones(n1)
102  self['h2f'] = np.ones(n1)
103  self['h31c'] = self['h2c'].view()
104  self['h31f'] = self['h2f'].view()
105 
106  # 2-direction
107  pos_c = np.array(p.y )
108  pos_f = np.array(p.ys)
109  self['h32c'] = abs(pos_c)
110  self['h32f'] = abs(pos_f)
111 
112  # linear size elements
113  self['dx1c'] = np.ones(n1) * p.ds[0]
114  self['dx1f'] = np.ones(n1) * p.ds[0]
115  self['dx2c'] = np.ones(n2) * p.ds[1]
116  self['dx2f'] = np.ones(n2) * p.ds[1]
117  self['dx3c'] = np.ones(n3) * p.ds[2]
118  self['dx3f'] = np.ones(n3) * p.ds[2]
119 
120  # volume elements
121  self['dvol1c'] = np.ones(n1) * p.ds[0]
122  self['dvol1f'] = np.ones(n1) * p.ds[0]
123  self['dvol2c'] = np.empty_like(pos_c)
124  self['dvol2f'] = np.empty_like(pos_f)
125  for j in xrange(len(pos_c)-1):
126  self['dvol2c'][j] = 0.5 * abs( self['h32f'][j+1] * pos_f[j+1]
127  - self['h32f'][j ] * pos_f[j ] )
128  self['dvol2f'][j] = 0.5 * abs( self['h32c'][j ] * pos_c[j ]
129  - self['h32c'][j-1] * pos_c[j-1] )
130  self['dvol3c'] = np.ones(n3) * p.ds[2]
131  self['dvol3f'] = np.ones(n3) * p.ds[2]
132 
133  # area elements
134  self['dar1c'] = self['h31c'] * self['h2c']
135  self['dar1f'] = self['h31f'] * self['h2f']
136  self['dar2c'] = self['h31f'] * p.ds[0] / self['dvol1c']
137  self['dar2f'] = self['h31c'] * p.ds[0] / self['dvol1f']
138  self['dar31c'] = self['h2f'] * p.ds[0] / self['dvol1c']
139  self['dar31f'] = self['h2c'] * p.ds[0] / self['dvol1f']
140  self['dar32c'] = p.ds[1] / self['dvol2c']
141  self['dar32f'] = p.ds[1] / self['dvol2f']
142 
143  def init_spherical(self, p):
144  """Initialise geometric factors for a spherical coord. system."""
145 
146  n1, n2, n3 = p.ncell
147 
148  # 1-direction
149  rpos_c = np.array(p.x)
150  rpos_f = np.array(p.xs)
151  self['h2c'] = abs(rpos_c)
152  self['h2f'] = abs(rpos_f)
153  self['h31c'] = self['h2c'].view()
154  self['h31f'] = self['h2f'].view()
155 
156  # 2-direction
157  tpos_c = np.array(p.y)
158  tpos_f = np.array(p.ys)
159  self['h32c'] = abs(np.sin(tpos_c))
160  self['h32f'] = abs(np.sin(tpos_f))
161 
162  # linear size elements
163  self['dx1c'] = np.ones(n1) * p.ds[0]
164  self['dx1f'] = np.ones(n1) * p.ds[0]
165  self['dx2c'] = np.ones(n2) * p.ds[1]
166  self['dx2f'] = np.ones(n2) * p.ds[1]
167  self['dx3c'] = np.ones(n3) * p.ds[2]
168  self['dx3f'] = np.ones(n3) * p.ds[2]
169 
170  # volume elements
171  self['dvol1c'] = np.empty_like(rpos_c)
172  self['dvol1f'] = np.empty_like(rpos_f)
173  for i in xrange(len(rpos_c)-1):
174  self['dvol1c'] = (self['h2f'][i+1] * self['h31f'][i+1] * rpos_f[i+1] / 3.0
175  - self['h2f'][i ] * self['h31f'][i ] * rpos_f[i ] / 3.0)
176  self['dvol1f'] = (self['h2c'][i ] * self['h31c'][i ] * rpos_c[i ] / 3.0
177  - self['h2c'][i-1] * self['h31c'][i-1] * rpos_c[i-1] / 3.0)
178  self['dvol2c'] = np.empty_like(tpos_c)
179  self['dvol2f'] = np.empty_like(tpos_f)
180  for j in xrange(len(tpos_c)-1):
181  self['dvol2c'] = -np.cos(tpos_c[j+1]) + np.cos(tpos_c[j ])
182  self['dvol2f'] = -np.cos(tpos_f[j ]) + np.cos(tpos_f[j-1])
183  self['dvol3c'] = np.ones(n3) * p.ds[2]
184  self['dvol3f'] = np.ones(n3) * p.ds[2]
185 
186  # area elements
187  self['dar1c'] = self['h31c'] * self['h2c']
188  self['dar1f'] = self['h31f'] * self['h2f']
189  self['dar2c'] = self['h31f'] * p.ds[0] / self['dvol1c']
190  self['dar2f'] = self['h31c'] * p.ds[0] / self['dvol1f']
191  self['dar31c'] = self['h2f'] * p.ds[0] / self['dvol1c']
192  self['dar31f'] = self['h2c'] * p.ds[0] / self['dvol1f']
193  self['dar32c'] = p.ds[1] / self['dvol2c']
194  self['dar32f'] = p.ds[1] / self['dvol2f']