DISPATCH
slopes_scalar.py
1 # -*- coding: utf-8 -*-
2 """
3 Scalar version of slope limiters. No periodic wrapping. Meant only for tests
4 and validation of the optimized slopes.py
5 
6 Created on Sat Nov 10 11:27:58 2018
7 
8 @author: Aake
9 """
10 
11 import numpy as np
12 from numpy import abs, zeros, min, product, roll
13 from time import time
14 from random import random
15 import slopes as opt
16 
17 def slopes (xx, slope_type, timeit=False):
18  """
19  Computes slopes for hydro and radiative transfer
20  variables in 3 dimensions.
21  The different slope limiters are:
22  - 1 : Minmod
23  - 2 : Moncen
24  - 3 : Positivity preserving
25  - 4 : Harmonic average
26  - 5 : Van Albada & Van Leer
27  - 6 : Superbee
28  """
29  m = xx.shape
30  ndim = 3-m.count(1)
31  m = list(m)
32  m1 = (3,m[0],m[1],m[2])
33  df = np.zeros(m1)
34 
35  if timeit:
36  start=time()
37 
38  if(slope_type == 0):
39  return zeros(m1)
40 
41  elif slope_type==1: # MINMOD
42  for i in range(1,m[0]-1):
43  for j in range(1,m[1]-1):
44  for k in range(1,m[2]-1):
45 
46  dflmm = ( xx[i-1,j ,k ]-xx[i,j,k] )
47  dfrmm = ( xx[i+1,j ,k ]-xx[i,j,k] )
48  df[0,i,j,k] = MinMod(dfrmm, -dflmm)
49 
50  if (ndim > 1): # slopes along y direction
51  dfmlm = ( xx[i ,j-1,k ]-xx[i,j,k] )
52  dfmrm = ( xx[i ,j+1,k ]-xx[i,j,k] )
53  df[1,i,j,k] = MinMod(dfmrm, -dfmlm)
54 
55  if (ndim > 2): # slopes along z direction
56  dfmml = ( xx[i ,j ,k-1]-xx[i,j,k] )
57  dfmmr = ( xx[i ,j ,k+1]-xx[i,j,k] )
58  df[2,i,j,k] = MinMod(dfmmr, -dfmml)
59 
60  elif slope_type==2: # MONCEN
61 
62  for i in range(1,m[0]-1):
63  for j in range(1,m[1]-1):
64  for k in range(1,m[2]-1):
65 
66  dflmm = 2.0 * ( xx[i, j, k ]-xx[i-1,j ,k ] )
67  dfrmm = 2.0 * ( xx[i+1,j ,k ]-xx[i ,j ,k ] )
68  dfcen = ( xx[i+1,j ,k ]-xx[i-1,j ,k ] ) / 2.0
69  df[0,i,j,k] = MonCen(dflmm,dfrmm,dfcen)
70 
71  if (ndim > 1): # slopes along y direction
72  dfmlm = 2.0 * ( xx[i, j, k ]-xx[i ,j-1,k ] )
73  dfmrm = 2.0 * ( xx[i ,j+1,k ]-xx[i ,j ,k ] )
74  dfcen = ( xx[i ,j+1,k ]-xx[i ,j-1,k ] ) / 2.0
75  df[1,i,j,k] = MonCen(dfmlm,dfmrm,dfcen)
76 
77  if (ndim > 2): # slopes along z direction
78  dfmml = 2.0 * ( xx[i, j, k ]-xx[i ,j ,k-1] )
79  dfmmr = 2.0 * ( xx[i ,j ,k+1]-xx[i ,j ,k ] )
80  dfcen = ( xx[i ,j ,k+1]-xx[i ,j ,k-1] ) / 2.0
81  df[2,i,j,k] = MonCen(dfmml,dfmmr,dfcen)
82 
83  elif slope_type==3: # positivity preserving
84  for i in range(1,m[0]-1):
85  for j in range(1,m[1]-1):
86  for k in range(1,m[2]-1):
87 
88  dflll = xx[i-1,j-1,k-1]-xx[i,j,k]
89  dflml = xx[i-1,j ,k-1]-xx[i,j,k]
90  dflrl = xx[i-1,j+1,k-1]-xx[i,j,k]
91  dfmll = xx[i ,j-1,k-1]-xx[i,j,k]
92  dfmml = xx[i ,j ,k-1]-xx[i,j,k]
93  dfmrl = xx[i ,j+1,k-1]-xx[i,j,k]
94  dfrll = xx[i+1,j-1,k-1]-xx[i,j,k]
95  dfrml = xx[i+1,j ,k-1]-xx[i,j,k]
96  dfrrl = xx[i+1,j+1,k-1]-xx[i,j,k]
97 
98  dfllm = xx[i-1,j-1,k ]-xx[i,j,k]
99  dflmm = xx[i-1,j ,k ]-xx[i,j,k]
100  dflrm = xx[i-1,j+1,k ]-xx[i,j,k]
101  dfmlm = xx[i ,j-1,k ]-xx[i,j,k]
102  dfmmm = 0.0
103  dfmrm = xx[i ,j+1,k ]-xx[i,j,k]
104  dfrlm = xx[i+1,j-1,k ]-xx[i,j,k]
105  dfrmm = xx[i+1,j ,k ]-xx[i,j,k]
106  dfrrm = xx[i+1,j+1,k ]-xx[i,j,k]
107 
108  dfllr = xx[i-1,j-1,k+1]-xx[i,j,k]
109  dflmr = xx[i-1,j ,k+1]-xx[i,j,k]
110  dflrr = xx[i-1,j+1,k+1]-xx[i,j,k]
111  dfmlr = xx[i ,j-1,k+1]-xx[i,j,k]
112  dfmmr = xx[i ,j ,k+1]-xx[i,j,k]
113  dfmrr = xx[i ,j+1,k+1]-xx[i,j,k]
114  dfrlr = xx[i+1,j-1,k+1]-xx[i,j,k]
115  dfrmr = xx[i+1,j ,k+1]-xx[i,j,k]
116  dfrrr = xx[i+1,j+1,k+1]-xx[i,j,k]
117 
118  vmin = min([dflll,dflml,dflrl,dfmll,dfmml,dfmrl,dfrll,dfrml,dfrrl,
119  dfllm,dflmm,dflrm,dfmlm,dfmmm,dfmrm,dfrlm,dfrmm,dfrrm,
120  dfllr,dflmr,dflrr,dfmlr,dfmmr,dfmrr,dfrlr,dfrmr,dfrrr])
121  vmax = max([dflll,dflml,dflrl,dfmll,dfmml,dfmrl,dfrll,dfrml,dfrrl,
122  dfllm,dflmm,dflrm,dfmlm,dfmmm,dfmrm,dfrlm,dfrmm,dfrrm,
123  dfllr,dflmr,dflrr,dfmlr,dfmmr,dfmrr,dfrlr,dfrmr,dfrrr])
124 
125  dfx = 0.5*(xx[i+1,j,k]-xx[i-1,j,k])
126  dfy = 0.5*(xx[i,j+1,k]-xx[i,j-1,k])
127  dfz = 0.5*(xx[i,j,k+1]-xx[i,j,k-1])
128  dff = 0.5*(abs(dfx)+abs(dfy)+abs(dfz))
129 
130  dlim = np.minimum(1.0,np.minimum(abs(vmin),abs(vmax))/dff)
131 
132  df[0,i,j,k] = dlim*dfx
133  df[1,i,j,k] = dlim*dfy
134  df[2,i,j,k] = dlim*dfz
135 
136  elif slope_type==4: # HARMONIC AVERAGE
137 
138  for i in range(1,m[0]-1):
139  for j in range(1,m[1]-1):
140  for k in range(1,m[2]-1):
141 
142  dflmm = ( xx[i-1,j ,k ]-xx[i,j,k] )
143  dfrmm = ( xx[i+1,j ,k ]-xx[i,j,k] )
144  df[0,i,j,k] = MoyHarm(dfrmm, -dflmm)
145 
146  if (ndim > 1): # slopes along y direction
147  dfmlm = ( xx[i ,j-1,k ]-xx[i,j,k] )
148  dfmrm = ( xx[i ,j+1,k ]-xx[i,j,k] )
149  df[1,i,j,k] = MoyHarm(dfmrm, -dfmlm)
150 
151  if (ndim > 2): # slopes along z direction
152  dfmml = ( xx[i ,j ,k-1]-xx[i,j,k] )
153  dfmmr = ( xx[i ,j ,k+1]-xx[i,j,k] )
154  df[2,i,j,k] = MoyHarm(dfmmr, -dfmml)
155 
156  elif slope_type==5: # Van Albada Van Leer
157  for i in range(1,m[0]-1):
158  for j in range(1,m[1]-1):
159  for k in range(1,m[2]-1):
160 
161  dflmm = ( xx[i-1,j ,k ]-xx[i,j,k] )
162  dfrmm = ( xx[i+1,j ,k ]-xx[i,j,k] )
163  df[0,i,j,k] = VAVL(dfrmm, -dflmm )
164 
165  if(ndim > 1): # slopes along y direction
166  dfmlm = ( xx[i ,j-1,k ]-xx[i,j,k] )
167  dfmrm = ( xx[i ,j+1,k ]-xx[i,j,k] )
168  df[1,i,j,k] = VAVL(dfmrm, -dfmlm)
169 
170  if(ndim > 2): # slopes along z direction
171  dfmml = ( xx[i ,j ,k-1]-xx[i,j,k] )
172  dfmmr = ( xx[i ,j ,k+1]-xx[i,j,k] )
173  df[2,i,j,k] = VAVL(dfmmr, -dfmml)
174 
175  elif slope_type==6: # Superbee
176  pxxD = np.zeros(m1)
177  pxxG = np.zeros(m1)
178  for i in range(1,m[0]-1):
179  for j in range(1,m[1]-1):
180  for k in range(1,m[2]-1):
181 
182  dflmm = ( xx[i-1,j ,k ]-xx[i,j,k] )
183  dfrmm = ( xx[i+1,j ,k ]-xx[i,j,k] )
184  pxxD[i,j,k,0] = MinMod(dfrmm, -2.0*dflmm)
185  pxxG[i,j,k,0] = MinMod(2.0*dfrmm, -dflmm)
186 
187  if (ndim > 1): # slopes along y direction
188  dfmlm = ( xx[i ,j-1,k ]-xx[i,j,k] )
189  dfmrm = ( xx[i ,j+1,k ]-xx[i,j,k] )
190  pxxD[i,j,k,1] = MinMod(dfmrm, -2.0*dfmlm)
191  pxxG[i,j,k,1] = MinMod(2.0*dfmrm, -dfmlm)
192 
193  if (ndim > 2): # slopes along z direction
194  dfmml = ( xx[i ,j ,k-1]-xx[i,j,k] )
195  dfmmr = ( xx[i ,j ,k+1]-xx[i,j,k] )
196  pxxD[i,j,k,2] = MinMod(dfmmr, -2.0*dfmml)
197  pxxG[i,j,k,2] = MinMod(2.0*dfmmr, -dfmml)
198  return pxxD, pxxG
199 
200  elif slope_type>10: # otimized
201  df = opt.slopes(xx,slope_type,timeit=timeit)
202 
203  if timeit:
204  musec = (time()-start)/product(m[0:3])*1e6
205  print('slope_type={}: {:6.2f} microsec per cell'.format(slope_type,musec))
206 
207  return df
208 
209 def SIGN(b):
210  return 1.0 if b > 0.0 else -1.0
211 
212 def MinMod(pa,pb):
213  """ Minmod averaging function """
214  s = (0.5 if pa>0.0 else -0.5) + (0.5 if pb>0.0 else -0.5)
215  return min([abs(pa),abs(pb)])*s
216 
217 def MonCen (pa,pb,pc):
218  """ Monotonised centred averaging function """
219  if ((pa*pb < 0.0) or (pa*pc < 0.0) or (pb*pc < 0.0)):
220  mcen = 0.0
221  elif ((abs(pa) < abs(pb)) and (abs(pa) < abs(pc))):
222  mcen = pa
223  elif ((abs(pb) < abs(pa)) and (abs(pb) < abs(pc))):
224  mcen = pb
225  else:
226  mcen = pc
227  return mcen
228 
229 def MoyHarm(pa,pb):
230  """ Harmonic averaging function """
231  if (pa*pb > 0.0):
232  MMod = 2.*pa*pb/(pa+pb)
233  else:
234  MMod = 0.
235  return MMod
236 
237 def VAVL(pa,pb):
238  """ Van Albada & Van Leer averaging function """
239  if(pa*pb > 0.0):
240  e = 1.0e-8
241  num = ( pa*pa + e*e )*pb + ( pb*pb + e*e )*pa
242  den = ( pa*pa + pb*pb + 2.0*e*e )
243  vv = num/den
244  else:
245  vv = 0.0
246  return vv
247 
248 def make_random(n):
249  f=zeros((n,n,n))
250  for i in range(n):
251  for j in range(n):
252  for k in range(n):
253  f[i,j,k]=random()
254  return f
255 
256 def test(n=32):
257  f=make_random(n)
258  for i in range(1,6):
259  d1=slopes(f,i)
260  d2=slopes(f,i+10)
261  d=(d2-d1)[:,2:-2,2:-2,2:-2]
262  print(d.min(),d.max())
263  test2(n=4*n)
264 
265 def test2(n=64):
266  f=make_random(n)
267  start=time()
268  f-roll(f,1,1)
269  mus=(time()-start)*1e6/n**3
270  print('deriv by roll : {:7.4f} microsec per cell'.format(mus))
271  start=time()
272  f[0:-1,:,:]=f[0:-1,:,:]-f[1:,:,:]
273  mus=(time()-start)*1e6/n**3
274  print('deriv by index: {:7.4f} microsec per cell'.format(mus))
Definition: slopes.py:1