3 Scalar version of slope limiters. No periodic wrapping. Meant only for tests 4 and validation of the optimized slopes.py 6 Created on Sat Nov 10 11:27:58 2018 12 from numpy
import abs, zeros, min, product, roll
14 from random
import random
17 def slopes (xx, slope_type, timeit=False):
19 Computes slopes for hydro and radiative transfer 20 variables in 3 dimensions. 21 The different slope limiters are: 24 - 3 : Positivity preserving 25 - 4 : Harmonic average 26 - 5 : Van Albada & Van Leer 32 m1 = (3,m[0],m[1],m[2])
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):
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)
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)
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)
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):
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)
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)
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)
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):
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]
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]
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]
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]
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])
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))
130 dlim = np.minimum(1.0,np.minimum(abs(vmin),abs(vmax))/dff)
132 df[0,i,j,k] = dlim*dfx
133 df[1,i,j,k] = dlim*dfy
134 df[2,i,j,k] = dlim*dfz
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):
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)
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)
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)
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):
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 )
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)
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)
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):
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)
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)
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)
201 df = opt.slopes(xx,slope_type,timeit=timeit)
204 musec = (time()-start)/product(m[0:3])*1e6
205 print(
'slope_type={}: {:6.2f} microsec per cell'.format(slope_type,musec))
210 return 1.0
if b > 0.0
else -1.0
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
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)):
221 elif ((abs(pa) < abs(pb))
and (abs(pa) < abs(pc))):
223 elif ((abs(pb) < abs(pa))
and (abs(pb) < abs(pc))):
230 """ Harmonic averaging function """ 232 MMod = 2.*pa*pb/(pa+pb)
238 """ Van Albada & Van Leer averaging function """ 241 num = ( pa*pa + e*e )*pb + ( pb*pb + e*e )*pa
242 den = ( pa*pa + pb*pb + 2.0*e*e )
261 d=(d2-d1)[:,2:-2,2:-2,2:-2]
262 print(d.min(),d.max())
269 mus=(time()-start)*1e6/n**3
270 print(
'deriv by roll : {:7.4f} microsec per cell'.format(mus))
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))