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. 6 Created on Sat Nov 10 11:27:58 2018 10 from numpy
import abs, zeros, ones, product, roll, where, min, minimum, maximum, asarray
12 from random
import random
13 import slopes_scalar
as scalar
17 d[1:,:,:]=f[1:,:,:]-f[0:-1,:,:]
18 d[0,:,:]=f[0,:,:]-f[-1,:,:]
23 d[0:-1,:,:]=f[1:,:,:]-f[0:-1,:,:]
24 d[-1,:,:]=f[0,:,:]-f[-1,:,:]
29 d[:,1:,:]=f[:,1:,:]-f[:,0:-1,:]
30 d[:,0,:]=f[:,0,:]-f[:,-1,:]
35 d[:,0:-1,:]=f[:,1:,:]-f[:,0:-1,:]
36 d[:,-1,:]=f[:,0,:]-f[:,-1,:]
41 d[:,:,1:]=f[:,:,1:]-f[:,:,0:-1]
42 d[:,:,0]=f[:,:,0]-f[:,:,-1]
47 d[:,:,0:-1]=f[:,:,1:]-f[:,:,0:-1]
48 d[:,:,-1]=f[:,:,0]-f[:,:,-1]
67 def slopes (xx, slope_type, timeit=False):
69 Computes slopes for hydro and radiative transfer 70 variables in 3 dimensions. 71 The different slope limiters are: 74 - 3 : Positivity preserving 75 - 4 : Harmonic average 76 - 5 : Van Albada & Van Leer 82 m1 = (3,m[0],m[1],m[2])
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
101 for i
in range(ndim):
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
112 dm = min([abs(dn),abs(up),abs(cn)],axis=0)
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
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
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
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
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])
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))
163 dlim = minimum(abs(vmin),abs(vmax))/dff
164 dlim = minimum(1.0,dlim)
171 for i
in range(ndim):
174 df[i] = maximum(dn*up,0.0)*2./(dn+up)
177 for i
in range(ndim):
181 num = (pa*pa + e*e)*pb + (pb*pb + e*e)*pa
182 den = (pa*pa + pb*pb + 2.0*e*e )
184 vv[where(pa*pb < 0.0)] = 0.0
188 df = scalar.slopes(xx,slope_type-10,timeit=timeit)
191 musec = (time()-start)/product(m[0:3])*1e6
192 print(
'slope_type={}: {:6.2f} microsec per cell'.format(slope_type,musec))
209 d=(d2-d1)[:,2:-2,2:-2,2:-2]
210 print(d.min(),d.max())