CG  Version 25
sm/src/pml.h
Go to the documentation of this file.
1 ! ******** This file generated from pmlUpdate.h using pml.maple *****
2 
3 c ** NOTE: run pml.maple to take this file and generate the pml.h file
4 c restart; read "pml.maple";
5 c ====================================================================================================
6 c Fourth-order update on a side
7 c OPTION : fullUpdate or partialUpdate
8 c 2 : 2 or 3 space dimensions
9 c ====================================================================================================
10 #beginMacro update4x2d(va,wa, m,OPTION)
11 
12 
13  ! (u(n+1) - 2*u + u(n-1))/dt^2 = u_tt + (dt^4/12)*u_tttt + ...
14  !
15  ! [u(n)-u(n-1)]/dt = u_t + (dt/2)*u_tt + O(dt^2)
16  !
17  ! u_tt = Delta u - v_x - w
18  ! u_tttt = Delta u_tt - v_xtt - wtt
19  !
20 
21  v=va (i1,i2,i3,m)
22  vx = va x42r(i1,i2,i3,m)
23  vxx = va xx42r(i1,i2,i3,m)
24  vxxx= va xxx22r(i1,i2,i3,m)
25  vxyy= va xyy22r(i1,i2,i3,m)
26 
27  w=wa(i1,i2,i3,m)
28  wx = wa x42r(i1,i2,i3,m)
29  wxx = wa xx42r(i1,i2,i3,m)
30 
31  ux= ux42r(i1,i2,i3,m)
32  uxx= uxx42r(i1,i2,i3,m)
33  uxxx=uxxx22r(i1,i2,i3,m)
34  uxyy=uxyy22r(i1,i2,i3,m)
35 
36  uxxxx=uxxxx22r(i1,i2,i3,m)
37  uxxyy=uxxyy22r(i1,i2,i3,m)
38 
39  uLap = uLaplacian42r(i1,i2,i3,m)
40 
41  ! --- these change in 3D ---
42  #If "2" == "2"
43  uyyyy=uyyyy22r(i1,i2,i3,m)
44  uLapSq=uxxxx +2.*uxxyy +uyyyy
45  uLapx = uxxx+uxyy
46  uLapxx= uxxxx+uxxyy
47  vLapx=vxxx+vxyy
48  #Elif "2" == "3"
49  uLapSq=uLapSq23r(i1,i2,i3,m)
50  uxxx=uxxx23r(i1,i2,i3,m)
51  uxxxx=uxxx23r(i1,i2,i3,m)
52  vxxx= va xxx23r(i1,i2,i3,m)
53 
54  uLapx = uxxx+uxyy+uxxx
55  uLapxx= uxxxx+uxxyy+uxxxx
56  vLapx=vxxx+vxyy+vxxx
57  #Else
58  stop 111999
59  #End
60 
61  ut = (u(i1,i2,i3,m)-um(i1,i2,i3,m))/dt - (.5*dt*csq)*( uLap - vx - w )
62  uxt = ( ux-umx42r(i1,i2,i3,m))/dt - (.5*dt*csq)*( uLapx - vxx - wx )
63  uxxt= (uxx-umxx42r(i1,i2,i3,m))/dt - (.5*dt*csq)*( uLapxx - vxxx - wxx )
64  ! *** uxxxt= (uxxx-umxxx42r(i1,i2,i3,m))/dt ! only need to first order in dt
65  ! *** uxyyt= (uxyy-umxyy42r(i1,i2,i3,m))/dt ! only need to first order in dt
66  ! *** uxxxxt= (uxxxx-umxxxx42r(i1,i2,i3,m))/dt ! only need to first order in dt
67  ! *** uxxyyt= (uxxyy-umxxyy42r(i1,i2,i3,m))/dt ! only need to first order in dt
68 
69  vt = sigma1*( -v + ux )
70  vxt = sigma1*( -vx + uxx ) + sigma1x*( -v + ux )
71  vxtt = sigma1*( -vxt + uxxt ) + sigma1x*( -vt + uxt )
72 
73  wt = sigma1*( -w -vx + uxx )
74  wtt = sigma1*( -wt -vxt + uxxt )
75 
76  #If #OPTION == "fullUpdate"
77  un(i1,i2,i3,m)=2.*u(i1,i2,i3,m)-um(i1,i2,i3,m) \
78  + cdtsq*( uLap - vx -w ) \
79  + cdt4Over12*( uLapSq - vLapx - wa Laplacian42r(i1,i2,i3,m) - vxtt - wtt )
80  #Elif #OPTION == "partialUpdate"
81  ! on an edge just add the other terms
82  un(i1,i2,i3,m)=un(i1,i2,i3,m)\
83  + cdtsq*( - vx -w ) \
84  + cdt4Over12*( - vLapx - wa Laplacian42r(i1,i2,i3,m) - vxtt - wtt )
85  #Else
86  stop 88437
87  #End
88 
89  ! auxilliary variables
90  ! v_t = sigma1*( -v + u_x )
91  ! (v(n+1)-v(n-1))/(2*dt) = v_t + (dt^2/3)*vttt
92  ! vttt = sigma1*( -v_tt + u_xtt )
93 
94  uxtt = csq*( uLapx - vxx -wx )
95  uxxtt = csq*( uLapxx - vxxx -wxx )
96 
97  ! new:
98  ! *** uxttt = csq*( uxxxt +uxyyt - vxxt -wxt )
99  ! *** uxxttt = csq*( uxxxxt+uxxyyt - vxxxt -wxxt )
100 
101  vtt = sigma1*( -vt + uxt )
102  vttt = sigma1*( -vtt + uxtt )
103  vtttt = 0. ! *** sigma1*( -vttt + uxttt)
104 
105  ! (v(n+1)-v(n-1))/(2dt) = vt + (dt^2/6)*vttt
106  ! va n(i1,i2,i3,m)=va m(i1,i2,i3,m)+(2.*dt)*( vt + (dt**2/6.)*vttt )
107  va n(i1,i2,i3,m)=va(i1,i2,i3,m)+(dt)*( vt + dt*( .5*vtt + dt*( (1./6.)*vttt + dt*( (1./24.)*vtttt ) ) ) )
108  ! write(*,'(" i1,i2=",2i3," vt,vtt,vttt=",3e10.2," v,ux,uxt,uxtt=",4e10.2)') i1,i2,vt,vtt,vttt,v,ux,uxt,uxtt
109 
110  ! w_t = sigma1*( -w -vx + uxx )
111 
112  wttt = sigma1*( -wtt -vxtt + uxxtt )
113  wtttt = 0. ! **** sigma1*( -wttt -vxttt + uxxttt )
114 ! wan(i1,i2,i3,m)=wam(i1,i2,i3,m)+(2.*dt)*( wt + (dt**2/6.)*wttt )
115  wa n(i1,i2,i3,m)=wa(i1,i2,i3,m)+(dt)*( wt + dt*( .5*wtt + dt*( (1./6.)*wttt + dt*( (1./24.)*wtttt ) ) ) )
116 
117 
118 #endMacro
119 
120 ! ******** This file generated from pmlUpdate.h using pml.maple *****
121 
122 c ** NOTE: run pml.maple to take this file and generate the pml.h file
123 c restart; read "pml.maple";
124 c ====================================================================================================
125 c Fourth-order update on a side
126 c OPTION : fullUpdate or partialUpdate
127 c 2 : 2 or 3 space dimensions
128 c ====================================================================================================
129 #beginMacro update4y2d(va,wa, m,OPTION)
130 
131 
132  ! (u(n+1) - 2*u + u(n-1))/dt^2 = u_tt + (dt^4/12)*u_tttt + ...
133  !
134  ! [u(n)-u(n-1)]/dt = u_t + (dt/2)*u_tt + O(dt^2)
135  !
136  ! u_tt = Delta u - v_y - w
137  ! u_tttt = Delta u_tt - v_ytt - wtt
138  !
139 
140  v=va (i1,i2,i3,m)
141  vy = va y42r(i1,i2,i3,m)
142  vyy = va yy42r(i1,i2,i3,m)
143  vyyy= va yyy22r(i1,i2,i3,m)
144  vxxy= va xxy22r(i1,i2,i3,m)
145 
146  w=wa(i1,i2,i3,m)
147  wy = wa y42r(i1,i2,i3,m)
148  wyy = wa yy42r(i1,i2,i3,m)
149 
150  uy= uy42r(i1,i2,i3,m)
151  uyy= uyy42r(i1,i2,i3,m)
152  uyyy=uyyy22r(i1,i2,i3,m)
153  uxxy=uxxy22r(i1,i2,i3,m)
154 
155  uyyyy=uyyyy22r(i1,i2,i3,m)
156  uxxyy=uxxyy22r(i1,i2,i3,m)
157 
158  uLap = uLaplacian42r(i1,i2,i3,m)
159 
160  ! --- these change in 3D ---
161  #If "2" == "2"
162  uxxxx=uxxxx22r(i1,i2,i3,m)
163  uLapSq=uyyyy +2.*uxxyy +uxxxx
164  uLapy = uyyy+uxxy
165  uLapyy= uyyyy+uxxyy
166  vLapy=vyyy+vxxy
167  #Elif "2" == "3"
168  uLapSq=uLapSq23r(i1,i2,i3,m)
169  uyyy=uyyy23r(i1,i2,i3,m)
170  uyyyy=uyyy23r(i1,i2,i3,m)
171  vyyy= va yyy23r(i1,i2,i3,m)
172 
173  uLapy = uyyy+uxxy+uyyy
174  uLapyy= uyyyy+uxxyy+uyyyy
175  vLapy=vyyy+vxxy+vyyy
176  #Else
177  stop 111999
178  #End
179 
180  ut = (u(i1,i2,i3,m)-um(i1,i2,i3,m))/dt - (.5*dt*csq)*( uLap - vy - w )
181  uyt = ( uy-umy42r(i1,i2,i3,m))/dt - (.5*dt*csq)*( uLapy - vyy - wy )
182  uyyt= (uyy-umyy42r(i1,i2,i3,m))/dt - (.5*dt*csq)*( uLapyy - vyyy - wyy )
183  ! *** uyyyt= (uyyy-umyyy42r(i1,i2,i3,m))/dt ! onlx need to first order in dt
184  ! *** uxxyt= (uxxy-umxxy42r(i1,i2,i3,m))/dt ! onlx need to first order in dt
185  ! *** uyyyyt= (uyyyy-umyyyy42r(i1,i2,i3,m))/dt ! onlx need to first order in dt
186  ! *** uxxyyt= (uxxyy-umxxyy42r(i1,i2,i3,m))/dt ! onlx need to first order in dt
187 
188  vt = sigma2*( -v + uy )
189  vyt = sigma2*( -vy + uyy ) + sigma2y*( -v + uy )
190  vytt = sigma2*( -vyt + uyyt ) + sigma2y*( -vt + uyt )
191 
192  wt = sigma2*( -w -vy + uyy )
193  wtt = sigma2*( -wt -vyt + uyyt )
194 
195  #If #OPTION == "fullUpdate"
196  un(i1,i2,i3,m)=2.*u(i1,i2,i3,m)-um(i1,i2,i3,m) \
197  + cdtsq*( uLap - vy -w ) \
198  + cdt4Over12*( uLapSq - vLapy - wa Laplacian42r(i1,i2,i3,m) - vytt - wtt )
199  #Elif #OPTION == "partialUpdate"
200  ! on an edge just add the other terms
201  un(i1,i2,i3,m)=un(i1,i2,i3,m)\
202  + cdtsq*( - vy -w ) \
203  + cdt4Over12*( - vLapy - wa Laplacian42r(i1,i2,i3,m) - vytt - wtt )
204  #Else
205  stop 88437
206  #End
207 
208  ! auyilliarx variables
209  ! v_t = sigma2*( -v + u_y )
210  ! (v(n+1)-v(n-1))/(2*dt) = v_t + (dt^2/3)*vttt
211  ! vttt = sigma2*( -v_tt + u_ytt )
212 
213  uytt = csq*( uLapy - vyy -wy )
214  uyytt = csq*( uLapyy - vyyy -wyy )
215 
216  ! new:
217  ! *** uyttt = csq*( uyyyt +uxxyt - vyyt -wyt )
218  ! *** uyyttt = csq*( uyyyyt+uxxyyt - vyyyt -wyyt )
219 
220  vtt = sigma2*( -vt + uyt )
221  vttt = sigma2*( -vtt + uytt )
222  vtttt = 0. ! *** sigma2*( -vttt + uyttt)
223 
224  ! (v(n+1)-v(n-1))/(2dt) = vt + (dt^2/6)*vttt
225  ! va n(i1,i2,i3,m)=va m(i1,i2,i3,m)+(2.*dt)*( vt + (dt**2/6.)*vttt )
226  va n(i1,i2,i3,m)=va(i1,i2,i3,m)+(dt)*( vt + dt*( .5*vtt + dt*( (1./6.)*vttt + dt*( (1./24.)*vtttt ) ) ) )
227  ! write(*,'(" i1,i2=",2i3," vt,vtt,vttt=",3e10.2," v,uy,uyt,uytt=",4e10.2)') i1,i2,vt,vtt,vttt,v,uy,uyt,uytt
228 
229  ! w_t = sigma2*( -w -vy + uyy )
230 
231  wttt = sigma2*( -wtt -vytt + uyytt )
232  wtttt = 0. ! **** sigma2*( -wttt -vyttt + uyyttt )
233 ! wan(i1,i2,i3,m)=wam(i1,i2,i3,m)+(2.*dt)*( wt + (dt**2/6.)*wttt )
234  wa n(i1,i2,i3,m)=wa(i1,i2,i3,m)+(dt)*( wt + dt*( .5*wtt + dt*( (1./6.)*wttt + dt*( (1./24.)*wtttt ) ) ) )
235 
236 
237 #endMacro
238 
239 ! ******** This file generated from pmlUpdate.h using pml.maple *****
240 
241 c ** NOTE: run pml.maple to take this file and generate the pml.h file
242 c restart; read "pml.maple";
243 c ====================================================================================================
244 c Fourth-order update on a side
245 c OPTION : fullUpdate or partialUpdate
246 c 3 : 2 or 3 space dimensions
247 c ====================================================================================================
248 #beginMacro update4x3d(va,wa, m,OPTION)
249 
250 
251  ! (u(n+1) - 2*u + u(n-1))/dt^2 = u_tt + (dt^4/12)*u_tttt + ...
252  !
253  ! [u(n)-u(n-1)]/dt = u_t + (dt/2)*u_tt + O(dt^2)
254  !
255  ! u_tt = Delta u - v_x - w
256  ! u_tttt = Delta u_tt - v_xtt - wtt
257  !
258 
259  v=va (i1,i2,i3,m)
260  vx = va x43r(i1,i2,i3,m)
261  vxx = va xx43r(i1,i2,i3,m)
262  vxxx= va xxx23r(i1,i2,i3,m)
263  vxyy= va xyy23r(i1,i2,i3,m)
264 
265  w=wa(i1,i2,i3,m)
266  wx = wa x43r(i1,i2,i3,m)
267  wxx = wa xx43r(i1,i2,i3,m)
268 
269  ux= ux43r(i1,i2,i3,m)
270  uxx= uxx43r(i1,i2,i3,m)
271  uxxx=uxxx23r(i1,i2,i3,m)
272  uxyy=uxyy23r(i1,i2,i3,m)
273 
274  uxxxx=uxxxx23r(i1,i2,i3,m)
275  uxxyy=uxxyy23r(i1,i2,i3,m)
276 
277  uLap = uLaplacian43r(i1,i2,i3,m)
278 
279  ! --- these change in 3D ---
280  #If "3" == "2"
281  uyyyy=uyyyy23r(i1,i2,i3,m)
282  uLapSq=uxxxx +2.*uxxyy +uyyyy
283  uLapx = uxxx+uxyy
284  uLapxx= uxxxx+uxxyy
285  vLapx=vxxx+vxyy
286  #Elif "3" == "3"
287  uLapSq=uLapSq23r(i1,i2,i3,m)
288  uxzz=uxzz23r(i1,i2,i3,m)
289  uxxzz=uxzz23r(i1,i2,i3,m)
290  vxzz= va xzz23r(i1,i2,i3,m)
291 
292  uLapx = uxxx+uxyy+uxzz
293  uLapxx= uxxxx+uxxyy+uxxzz
294  vLapx=vxxx+vxyy+vxzz
295  #Else
296  stop 111999
297  #End
298 
299  ut = (u(i1,i2,i3,m)-um(i1,i2,i3,m))/dt - (.5*dt*csq)*( uLap - vx - w )
300  uxt = ( ux-umx43r(i1,i2,i3,m))/dt - (.5*dt*csq)*( uLapx - vxx - wx )
301  uxxt= (uxx-umxx43r(i1,i2,i3,m))/dt - (.5*dt*csq)*( uLapxx - vxxx - wxx )
302  ! *** uxxxt= (uxxx-umxxx43r(i1,i2,i3,m))/dt ! only need to first order in dt
303  ! *** uxyyt= (uxyy-umxyy43r(i1,i2,i3,m))/dt ! only need to first order in dt
304  ! *** uxxxxt= (uxxxx-umxxxx43r(i1,i2,i3,m))/dt ! only need to first order in dt
305  ! *** uxxyyt= (uxxyy-umxxyy43r(i1,i2,i3,m))/dt ! only need to first order in dt
306 
307  vt = sigma1*( -v + ux )
308  vxt = sigma1*( -vx + uxx ) + sigma1x*( -v + ux )
309  vxtt = sigma1*( -vxt + uxxt ) + sigma1x*( -vt + uxt )
310 
311  wt = sigma1*( -w -vx + uxx )
312  wtt = sigma1*( -wt -vxt + uxxt )
313 
314  #If #OPTION == "fullUpdate"
315  un(i1,i2,i3,m)=2.*u(i1,i2,i3,m)-um(i1,i2,i3,m) \
316  + cdtsq*( uLap - vx -w ) \
317  + cdt4Over12*( uLapSq - vLapx - wa Laplacian43r(i1,i2,i3,m) - vxtt - wtt )
318  #Elif #OPTION == "partialUpdate"
319  ! on an edge just add the other terms
320  un(i1,i2,i3,m)=un(i1,i2,i3,m)\
321  + cdtsq*( - vx -w ) \
322  + cdt4Over12*( - vLapx - wa Laplacian43r(i1,i2,i3,m) - vxtt - wtt )
323  #Else
324  stop 88437
325  #End
326 
327  ! auxilliary variables
328  ! v_t = sigma1*( -v + u_x )
329  ! (v(n+1)-v(n-1))/(2*dt) = v_t + (dt^2/3)*vttt
330  ! vttt = sigma1*( -v_tt + u_xtt )
331 
332  uxtt = csq*( uLapx - vxx -wx )
333  uxxtt = csq*( uLapxx - vxxx -wxx )
334 
335  ! new:
336  ! *** uxttt = csq*( uxxxt +uxyyt - vxxt -wxt )
337  ! *** uxxttt = csq*( uxxxxt+uxxyyt - vxxxt -wxxt )
338 
339  vtt = sigma1*( -vt + uxt )
340  vttt = sigma1*( -vtt + uxtt )
341  vtttt = 0. ! *** sigma1*( -vttt + uxttt)
342 
343  ! (v(n+1)-v(n-1))/(2dt) = vt + (dt^2/6)*vttt
344  ! va n(i1,i2,i3,m)=va m(i1,i2,i3,m)+(2.*dt)*( vt + (dt**2/6.)*vttt )
345  va n(i1,i2,i3,m)=va(i1,i2,i3,m)+(dt)*( vt + dt*( .5*vtt + dt*( (1./6.)*vttt + dt*( (1./24.)*vtttt ) ) ) )
346  ! write(*,'(" i1,i2=",2i3," vt,vtt,vttt=",3e10.2," v,ux,uxt,uxtt=",4e10.2)') i1,i2,vt,vtt,vttt,v,ux,uxt,uxtt
347 
348  ! w_t = sigma1*( -w -vx + uxx )
349 
350  wttt = sigma1*( -wtt -vxtt + uxxtt )
351  wtttt = 0. ! **** sigma1*( -wttt -vxttt + uxxttt )
352 ! wan(i1,i2,i3,m)=wam(i1,i2,i3,m)+(2.*dt)*( wt + (dt**2/6.)*wttt )
353  wa n(i1,i2,i3,m)=wa(i1,i2,i3,m)+(dt)*( wt + dt*( .5*wtt + dt*( (1./6.)*wttt + dt*( (1./24.)*wtttt ) ) ) )
354 
355 
356 #endMacro
357 
358 ! ******** This file generated from pmlUpdate.h using pml.maple *****
359 
360 c ** NOTE: run pml.maple to take this file and generate the pml.h file
361 c restart; read "pml.maple";
362 c ====================================================================================================
363 c Fourth-order update on a side
364 c OPTION : fullUpdate or partialUpdate
365 c 3 : 2 or 3 space dimensions
366 c ====================================================================================================
367 #beginMacro update4y3d(va,wa, m,OPTION)
368 
369 
370  ! (u(n+1) - 2*u + u(n-1))/dt^2 = u_tt + (dt^4/12)*u_tttt + ...
371  !
372  ! [u(n)-u(n-1)]/dt = u_t + (dt/2)*u_tt + O(dt^2)
373  !
374  ! u_tt = Delta u - v_y - w
375  ! u_tttt = Delta u_tt - v_ytt - wtt
376  !
377 
378  v=va (i1,i2,i3,m)
379  vy = va y43r(i1,i2,i3,m)
380  vyy = va yy43r(i1,i2,i3,m)
381  vyyy= va yyy23r(i1,i2,i3,m)
382  vyzz= va yzz23r(i1,i2,i3,m)
383 
384  w=wa(i1,i2,i3,m)
385  wy = wa y43r(i1,i2,i3,m)
386  wyy = wa yy43r(i1,i2,i3,m)
387 
388  uy= uy43r(i1,i2,i3,m)
389  uyy= uyy43r(i1,i2,i3,m)
390  uyyy=uyyy23r(i1,i2,i3,m)
391  uyzz=uyzz23r(i1,i2,i3,m)
392 
393  uyyyy=uyyyy23r(i1,i2,i3,m)
394  uyyzz=uyyzz23r(i1,i2,i3,m)
395 
396  uLap = uLaplacian43r(i1,i2,i3,m)
397 
398  ! --- these change in 3D ---
399  #If "3" == "2"
400  uzzzz=uzzzz23r(i1,i2,i3,m)
401  uLapSq=uyyyy +2.*uyyzz +uzzzz
402  uLapy = uyyy+uyzz
403  uLapyy= uyyyy+uyyzz
404  vLapy=vyyy+vyzz
405  #Elif "3" == "3"
406  uLapSq=uLapSq23r(i1,i2,i3,m)
407  uxxy=uxxy23r(i1,i2,i3,m)
408  uxxyy=uxxy23r(i1,i2,i3,m)
409  vxxy= va xxy23r(i1,i2,i3,m)
410 
411  uLapy = uyyy+uyzz+uxxy
412  uLapyy= uyyyy+uyyzz+uxxyy
413  vLapy=vyyy+vyzz+vxxy
414  #Else
415  stop 111999
416  #End
417 
418  ut = (u(i1,i2,i3,m)-um(i1,i2,i3,m))/dt - (.5*dt*csq)*( uLap - vy - w )
419  uyt = ( uy-umy43r(i1,i2,i3,m))/dt - (.5*dt*csq)*( uLapy - vyy - wy )
420  uyyt= (uyy-umyy43r(i1,i2,i3,m))/dt - (.5*dt*csq)*( uLapyy - vyyy - wyy )
421  ! *** uyyyt= (uyyy-umyyy43r(i1,i2,i3,m))/dt ! onlz need to first order in dt
422  ! *** uyzzt= (uyzz-umyzz43r(i1,i2,i3,m))/dt ! onlz need to first order in dt
423  ! *** uyyyyt= (uyyyy-umyyyy43r(i1,i2,i3,m))/dt ! onlz need to first order in dt
424  ! *** uyyzzt= (uyyzz-umyyzz43r(i1,i2,i3,m))/dt ! onlz need to first order in dt
425 
426  vt = sigma2*( -v + uy )
427  vyt = sigma2*( -vy + uyy ) + sigma2y*( -v + uy )
428  vytt = sigma2*( -vyt + uyyt ) + sigma2y*( -vt + uyt )
429 
430  wt = sigma2*( -w -vy + uyy )
431  wtt = sigma2*( -wt -vyt + uyyt )
432 
433  #If #OPTION == "fullUpdate"
434  un(i1,i2,i3,m)=2.*u(i1,i2,i3,m)-um(i1,i2,i3,m) \
435  + cdtsq*( uLap - vy -w ) \
436  + cdt4Over12*( uLapSq - vLapy - wa Laplacian43r(i1,i2,i3,m) - vytt - wtt )
437  #Elif #OPTION == "partialUpdate"
438  ! on an edge just add the other terms
439  un(i1,i2,i3,m)=un(i1,i2,i3,m)\
440  + cdtsq*( - vy -w ) \
441  + cdt4Over12*( - vLapy - wa Laplacian43r(i1,i2,i3,m) - vytt - wtt )
442  #Else
443  stop 88437
444  #End
445 
446  ! auyilliarz variables
447  ! v_t = sigma2*( -v + u_y )
448  ! (v(n+1)-v(n-1))/(2*dt) = v_t + (dt^2/3)*vttt
449  ! vttt = sigma2*( -v_tt + u_ytt )
450 
451  uytt = csq*( uLapy - vyy -wy )
452  uyytt = csq*( uLapyy - vyyy -wyy )
453 
454  ! new:
455  ! *** uyttt = csq*( uyyyt +uyzzt - vyyt -wyt )
456  ! *** uyyttt = csq*( uyyyyt+uyyzzt - vyyyt -wyyt )
457 
458  vtt = sigma2*( -vt + uyt )
459  vttt = sigma2*( -vtt + uytt )
460  vtttt = 0. ! *** sigma2*( -vttt + uyttt)
461 
462  ! (v(n+1)-v(n-1))/(2dt) = vt + (dt^2/6)*vttt
463  ! va n(i1,i2,i3,m)=va m(i1,i2,i3,m)+(2.*dt)*( vt + (dt**2/6.)*vttt )
464  va n(i1,i2,i3,m)=va(i1,i2,i3,m)+(dt)*( vt + dt*( .5*vtt + dt*( (1./6.)*vttt + dt*( (1./24.)*vtttt ) ) ) )
465  ! write(*,'(" i1,i2=",2i3," vt,vtt,vttt=",3e10.2," v,uy,uyt,uytt=",4e10.2)') i1,i2,vt,vtt,vttt,v,uy,uyt,uytt
466 
467  ! w_t = sigma2*( -w -vy + uyy )
468 
469  wttt = sigma2*( -wtt -vytt + uyytt )
470  wtttt = 0. ! **** sigma2*( -wttt -vyttt + uyyttt )
471 ! wan(i1,i2,i3,m)=wam(i1,i2,i3,m)+(2.*dt)*( wt + (dt**2/6.)*wttt )
472  wa n(i1,i2,i3,m)=wa(i1,i2,i3,m)+(dt)*( wt + dt*( .5*wtt + dt*( (1./6.)*wttt + dt*( (1./24.)*wtttt ) ) ) )
473 
474 
475 #endMacro
476 
477 ! ******** This file generated from pmlUpdate.h using pml.maple *****
478 
479 c ** NOTE: run pml.maple to take this file and generate the pml.h file
480 c restart; read "pml.maple";
481 c ====================================================================================================
482 c Fourth-order update on a side
483 c OPTION : fullUpdate or partialUpdate
484 c 3 : 2 or 3 space dimensions
485 c ====================================================================================================
486 #beginMacro update4z3d(va,wa, m,OPTION)
487 
488 
489  ! (u(n+1) - 2*u + u(n-1))/dt^2 = u_tt + (dt^4/12)*u_tttt + ...
490  !
491  ! [u(n)-u(n-1)]/dt = u_t + (dt/2)*u_tt + O(dt^2)
492  !
493  ! u_tt = Delta u - v_z - w
494  ! u_tttt = Delta u_tt - v_ztt - wtt
495  !
496 
497  v=va (i1,i2,i3,m)
498  vz = va z43r(i1,i2,i3,m)
499  vzz = va zz43r(i1,i2,i3,m)
500  vzzz= va zzz23r(i1,i2,i3,m)
501  vxxz= va xxz23r(i1,i2,i3,m)
502 
503  w=wa(i1,i2,i3,m)
504  wz = wa z43r(i1,i2,i3,m)
505  wzz = wa zz43r(i1,i2,i3,m)
506 
507  uz= uz43r(i1,i2,i3,m)
508  uzz= uzz43r(i1,i2,i3,m)
509  uzzz=uzzz23r(i1,i2,i3,m)
510  uxxz=uxxz23r(i1,i2,i3,m)
511 
512  uzzzz=uzzzz23r(i1,i2,i3,m)
513  uxxzz=uxxzz23r(i1,i2,i3,m)
514 
515  uLap = uLaplacian43r(i1,i2,i3,m)
516 
517  ! --- these change in 3D ---
518  #If "3" == "2"
519  uxxxx=uxxxx23r(i1,i2,i3,m)
520  uLapSq=uzzzz +2.*uxxzz +uxxxx
521  uLapz = uzzz+uxxz
522  uLapzz= uzzzz+uxxzz
523  vLapz=vzzz+vxxz
524  #Elif "3" == "3"
525  uLapSq=uLapSq23r(i1,i2,i3,m)
526  uyyz=uyyz23r(i1,i2,i3,m)
527  uyyzz=uyyz23r(i1,i2,i3,m)
528  vyyz= va yyz23r(i1,i2,i3,m)
529 
530  uLapz = uzzz+uxxz+uyyz
531  uLapzz= uzzzz+uxxzz+uyyzz
532  vLapz=vzzz+vxxz+vyyz
533  #Else
534  stop 111999
535  #End
536 
537  ut = (u(i1,i2,i3,m)-um(i1,i2,i3,m))/dt - (.5*dt*csq)*( uLap - vz - w )
538  uzt = ( uz-umz43r(i1,i2,i3,m))/dt - (.5*dt*csq)*( uLapz - vzz - wz )
539  uzzt= (uzz-umzz43r(i1,i2,i3,m))/dt - (.5*dt*csq)*( uLapzz - vzzz - wzz )
540  ! *** uzzzt= (uzzz-umzzz43r(i1,i2,i3,m))/dt ! onlx need to first order in dt
541  ! *** uxxzt= (uxxz-umxxz43r(i1,i2,i3,m))/dt ! onlx need to first order in dt
542  ! *** uzzzzt= (uzzzz-umzzzz43r(i1,i2,i3,m))/dt ! onlx need to first order in dt
543  ! *** uxxzzt= (uxxzz-umxxzz43r(i1,i2,i3,m))/dt ! onlx need to first order in dt
544 
545  vt = sigma3*( -v + uz )
546  vzt = sigma3*( -vz + uzz ) + sigma3z*( -v + uz )
547  vztt = sigma3*( -vzt + uzzt ) + sigma3z*( -vt + uzt )
548 
549  wt = sigma3*( -w -vz + uzz )
550  wtt = sigma3*( -wt -vzt + uzzt )
551 
552  #If #OPTION == "fullUpdate"
553  un(i1,i2,i3,m)=2.*u(i1,i2,i3,m)-um(i1,i2,i3,m) \
554  + cdtsq*( uLap - vz -w ) \
555  + cdt4Over12*( uLapSq - vLapz - wa Laplacian43r(i1,i2,i3,m) - vztt - wtt )
556  #Elif #OPTION == "partialUpdate"
557  ! on an edge just add the other terms
558  un(i1,i2,i3,m)=un(i1,i2,i3,m)\
559  + cdtsq*( - vz -w ) \
560  + cdt4Over12*( - vLapz - wa Laplacian43r(i1,i2,i3,m) - vztt - wtt )
561  #Else
562  stop 88437
563  #End
564 
565  ! auzilliarx variables
566  ! v_t = sigma3*( -v + u_z )
567  ! (v(n+1)-v(n-1))/(2*dt) = v_t + (dt^2/3)*vttt
568  ! vttt = sigma3*( -v_tt + u_ztt )
569 
570  uztt = csq*( uLapz - vzz -wz )
571  uzztt = csq*( uLapzz - vzzz -wzz )
572 
573  ! new:
574  ! *** uzttt = csq*( uzzzt +uxxzt - vzzt -wzt )
575  ! *** uzzttt = csq*( uzzzzt+uxxzzt - vzzzt -wzzt )
576 
577  vtt = sigma3*( -vt + uzt )
578  vttt = sigma3*( -vtt + uztt )
579  vtttt = 0. ! *** sigma3*( -vttt + uzttt)
580 
581  ! (v(n+1)-v(n-1))/(2dt) = vt + (dt^2/6)*vttt
582  ! va n(i1,i2,i3,m)=va m(i1,i2,i3,m)+(2.*dt)*( vt + (dt**2/6.)*vttt )
583  va n(i1,i2,i3,m)=va(i1,i2,i3,m)+(dt)*( vt + dt*( .5*vtt + dt*( (1./6.)*vttt + dt*( (1./24.)*vtttt ) ) ) )
584  ! write(*,'(" i1,i2=",2i3," vt,vtt,vttt=",3e10.2," v,uz,uzt,uztt=",4e10.2)') i1,i2,vt,vtt,vttt,v,uz,uzt,uztt
585 
586  ! w_t = sigma3*( -w -vz + uzz )
587 
588  wttt = sigma3*( -wtt -vztt + uzztt )
589  wtttt = 0. ! **** sigma3*( -wttt -vzttt + uzzttt )
590 ! wan(i1,i2,i3,m)=wam(i1,i2,i3,m)+(2.*dt)*( wt + (dt**2/6.)*wttt )
591  wa n(i1,i2,i3,m)=wa(i1,i2,i3,m)+(dt)*( wt + dt*( .5*wtt + dt*( (1./6.)*wttt + dt*( (1./24.)*wtttt ) ) ) )
592 
593 
594 #endMacro
595 
Generated on Fri Jan 4 2013 10:19:26 for CG by   doxygen 1.8.3