DISPATCH
slopes.py
1 # -*- coding: utf-8 -*-
2 """
3 Optimized version of slope limiters, with periodic wrapping. The slope limiters
4 are some of the most common ones used in Godunov type Riemann solvers.
5 
6 Created on Sat Nov 10 11:27:58 2018
7 
8 @author: Aake
9 """
10 from numpy import abs, zeros, ones, product, roll, where, min, minimum, maximum, asarray
11 from time import time
12 from random import random
13 import slopes_scalar as scalar
14 
15 def dxdn(f):
16  d=zeros(f.shape)
17  d[1:,:,:]=f[1:,:,:]-f[0:-1,:,:]
18  d[0,:,:]=f[0,:,:]-f[-1,:,:]
19  return d
20 
21 def dxup(f):
22  d=zeros(f.shape)
23  d[0:-1,:,:]=f[1:,:,:]-f[0:-1,:,:]
24  d[-1,:,:]=f[0,:,:]-f[-1,:,:]
25  return d
26 
27 def dydn(f):
28  d=zeros(f.shape)
29  d[:,1:,:]=f[:,1:,:]-f[:,0:-1,:]
30  d[:,0,:]=f[:,0,:]-f[:,-1,:]
31  return d
32 
33 def dyup(f):
34  d=zeros(f.shape)
35  d[:,0:-1,:]=f[:,1:,:]-f[:,0:-1,:]
36  d[:,-1,:]=f[:,0,:]-f[:,-1,:]
37  return d
38 
39 def dzdn(f):
40  d=zeros(f.shape)
41  d[:,:,1:]=f[:,:,1:]-f[:,:,0:-1]
42  d[:,:,0]=f[:,:,0]-f[:,:,-1]
43  return d
44 
45 def dzup(f):
46  d=zeros(f.shape)
47  d[:,:,0:-1]=f[:,:,1:]-f[:,:,0:-1]
48  d[:,:,-1]=f[:,:,0]-f[:,:,-1]
49  return d
50 
51 def ddn(f,i):
52  if i==0:
53  return dxdn(f)
54  if i==1:
55  return dydn(f)
56  if i==2:
57  return dzdn(f)
58 
59 def dup(f,i):
60  if i==0:
61  return dxup(f)
62  if i==1:
63  return dyup(f)
64  if i==2:
65  return dzup(f)
66 
67 def slopes (xx, slope_type, timeit=False):
68  """
69  Computes slopes for hydro and radiative transfer
70  variables in 3 dimensions.
71  The different slope limiters are:
72  - 1 : Minmod
73  - 2 : Moncen
74  - 3 : Positivity preserving
75  - 4 : Harmonic average
76  - 5 : Van Albada & Van Leer
77  - 6 : Superbee
78  """
79  m = xx.shape
80  ndim = 3-m.count(1)
81  m = list(m)
82  m1 = (3,m[0],m[1],m[2])
83  df = zeros(m1)
84 
85  if timeit:
86  start=time()
87 
88  if(slope_type == 0):
89  return zeros(m1)
90 
91  elif slope_type==1: # MINMOD
92  for i in range(ndim):
93  dn = ddn(xx,i)
94  up = dup(xx,i)
95  s = ones(xx.shape)
96  s[where(dn < 0.0)] -= 1.0
97  s[where(up < 0.0)] -= 1.0
98  df[i] = minimum(abs(dn),abs(up))*s
99 
100  elif slope_type==2: # MONCEN
101  for i in range(ndim):
102  # Down, up, and center slopes
103  dn = ddn(xx,i)*2.0
104  up = dup(xx,i)*2.0
105  cn = (dn+up)*0.25
106  # Switch, vanishes when one of three slopes differs
107  sw = ones(xx.shape)
108  sw[where(dn*up < 0.0)] = 0.0
109  sw[where(dn*cn < 0.0)] = 0.0
110  sw[where(up*cn < 0.0)] = 0.0
111  # Slope is smallest of the three
112  dm = min([abs(dn),abs(up),abs(cn)],axis=0)
113  d = cn
114  w = where(dm==abs(dn)); d[w] = dn[w]
115  w = where(dm==abs(up)); d[w] = up[w]
116  d[where(cn*up < 0.0)] = 0.0
117  d[where(cn*dn < 0.0)] = 0.0
118  d[where(dn*up < 0.0)] = 0.0
119  df[i] = d
120 
121  elif slope_type==3: # positivity preserving
122  dflll = roll(xx,( 1, 1, 1),(0,1,2))-xx
123  dflml = roll(xx,( 1, 0, 1),(0,1,2))-xx
124  dflrl = roll(xx,( 1,-1, 1),(0,1,2))-xx
125  dfmll = roll(xx,( 0, 1, 1),(0,1,2))-xx
126  dfmml = roll(xx,( 0, 0, 1),(0,1,2))-xx
127  dfmrl = roll(xx,( 0,-1, 1),(0,1,2))-xx
128  dfrll = roll(xx,(-1, 1, 1),(0,1,2))-xx
129  dfrml = roll(xx,(-1, 0, 1),(0,1,2))-xx
130  dfrrl = roll(xx,(-1,-1, 1),(0,1,2))-xx
131 
132  dfllm = roll(xx,( 1, 1, 0),(0,1,2))-xx
133  dflmm = roll(xx,( 1, 0, 0),(0,1,2))-xx
134  dflrm = roll(xx,( 1,-1, 0),(0,1,2))-xx
135  dfmlm = roll(xx,( 0, 1, 0),(0,1,2))-xx
136  dfmmm = roll(xx,( 0, 0, 0),(0,1,2))-xx
137  dfmrm = roll(xx,( 0,-1, 0),(0,1,2))-xx
138  dfrlm = roll(xx,(-1, 1, 0),(0,1,2))-xx
139  dfrmm = roll(xx,(-1, 0, 0),(0,1,2))-xx
140  dfrrm = roll(xx,(-1,-1, 0),(0,1,2))-xx
141 
142  dfllr = roll(xx,( 1, 1,-1),(0,1,2))-xx
143  dflmr = roll(xx,( 1, 0,-1),(0,1,2))-xx
144  dflrr = roll(xx,( 1,-1,-1),(0,1,2))-xx
145  dfmlr = roll(xx,( 0, 1,-1),(0,1,2))-xx
146  dfmmr = roll(xx,( 0, 0,-1),(0,1,2))-xx
147  dfmrr = roll(xx,( 0,-1,-1),(0,1,2))-xx
148  dfrlr = roll(xx,(-1, 1,-1),(0,1,2))-xx
149  dfrmr = roll(xx,(-1, 0,-1),(0,1,2))-xx
150  dfrrr = roll(xx,(-1,-1,-1),(0,1,2))-xx
151 
152  V = asarray([dflll,dflml,dflrl,dfmll,dfmml,dfmrl,dfrll,dfrml,dfrrl,
153  dfllm,dflmm,dflrm,dfmlm,dfmmm,dfmrm,dfrlm,dfrmm,dfrrm,
154  dfllr,dflmr,dflrr,dfmlr,dfmmr,dfmrr,dfrlr,dfrmr,dfrrr])
155  vmin=V.min(0)
156  vmax=V.max(0)
157 
158  dfx = 0.5*(dfrmm-dflmm)
159  dfy = 0.5*(dfmrm-dfmlm)
160  dfz = 0.5*(dfmmr-dfmml)
161  dff = 0.5*(abs(dfx)+abs(dfy)+abs(dfz))
162 
163  dlim = minimum(abs(vmin),abs(vmax))/dff
164  dlim = minimum(1.0,dlim)
165 
166  df[0] = dlim*dfx
167  df[1] = dlim*dfy
168  df[2] = dlim*dfz
169 
170  elif slope_type==4: # HARMONIC AVERAGE
171  for i in range(ndim):
172  dn = ddn(xx,i)
173  up = dup(xx,i)
174  df[i] = maximum(dn*up,0.0)*2./(dn+up)
175 
176  elif slope_type==5: # Van Albada Van Leer
177  for i in range(ndim):
178  pa = dup(xx,i)
179  pb = ddn(xx,i)
180  e = 1.0e-8
181  num = (pa*pa + e*e)*pb + (pb*pb + e*e)*pa
182  den = (pa*pa + pb*pb + 2.0*e*e )
183  vv = num/den
184  vv[where(pa*pb < 0.0)] = 0.0
185  df[i] = vv
186 
187  elif slope_type>10: # otimized
188  df = scalar.slopes(xx,slope_type-10,timeit=timeit)
189 
190  if timeit:
191  musec = (time()-start)/product(m[0:3])*1e6
192  print('slope_type={}: {:6.2f} microsec per cell'.format(slope_type,musec))
193 
194  return df
195 
196 def make_random(n):
197  f=zeros((n,n,n))
198  for i in range(n):
199  for j in range(n):
200  for k in range(n):
201  f[i,j,k]=random()
202  return f
203 
204 def test(n=32):
205  f=make_random(n)
206  for i in range(1,6):
207  d1=slopes(f,i)
208  d2=slopes(f,i+10)
209  d=(d2-d1)[:,2:-2,2:-2,2:-2]
210  print(d.min(),d.max())
Definition: slopes.py:1