CG  Version 25
secondOrderCurvilinear2D.h
Go to the documentation of this file.
1  t2 = cc ** 2
2  t5 = rx(i + 1,j,0,0) * rx(i + 1,j,1,1) - rx(i + 1,j,1,0) * rx(i +
3  #1,j,0,1)
4  t6 = 0.1D1 / t5
5  t7 = rx(i + 1,j,0,0) ** 2
6  t8 = rx(i + 1,j,0,1) ** 2
7  t10 = t6 * (t7 + t8)
8  t13 = rx(i,j,0,0) * rx(i,j,1,1) - rx(i,j,1,0) * rx(i,j,0,1)
9  t14 = 0.1D1 / t13
10  t15 = rx(i,j,0,0) ** 2
11  t16 = rx(i,j,0,1) ** 2
12  t18 = t14 * (t15 + t16)
13  t19 = t10 + t18
14  t20 = t2 * t19 / 0.2D1
15  t22 = 0.1D1 / dx
16  t23 = (u(i + 1,j,n) - u(i,j,n)) * t22
17  t24 = t20 * t23
18  t25 = ut(i + 1,j,n) - ut(i,j,n)
19  t34 = 0.1D1 / (rx(i + 2,j,0,0) * rx(i + 2,j,1,1) - rx(i + 2,j,1,0)
20  # * rx(i + 2,j,0,1))
21  t35 = rx(i + 2,j,0,0) ** 2
22  t36 = rx(i + 2,j,0,1) ** 2
23  t42 = (u(i + 2,j,n) - u(i + 1,j,n)) * t22
24  t51 = 0.1D1 / dy
25  t58 = t2 * t6
26  t61 = rx(i + 1,j,0,0) * rx(i + 1,j,1,0) + rx(i + 1,j,0,1) * rx(i +
27  # 1,j,1,1)
28  t63 = (u(i + 1,j + 1,n) - u(i + 1,j,n)) * t51
29  t65 = (u(i + 1,j,n) - u(i + 1,j - 1,n)) * t51
30  t68 = t58 * t61 * (t63 + t65) / 0.2D1
31  t72 = t2 * t14
32  t75 = rx(i,j,0,0) * rx(i,j,1,0) + rx(i,j,0,1) * rx(i,j,1,1)
33  t77 = (u(i,j + 1,n) - u(i,j,n)) * t51
34  t79 = (u(i,j,n) - u(i,j - 1,n)) * t51
35  t82 = t72 * t75 * (t77 + t79) / 0.2D1
36  t85 = (t68 - t82) * t22 / 0.2D1
37  t89 = 0.1D1 / (rx(i + 1,j + 1,0,0) * rx(i + 1,j + 1,1,1) - rx(i +
38  #1,j + 1,1,0) * rx(i + 1,j + 1,0,1))
39  t90 = t2 * t89
40  t93 = rx(i + 1,j + 1,0,0) * rx(i + 1,j + 1,1,0) + rx(i + 1,j + 1,0
41  #,1) * rx(i + 1,j + 1,1,1)
42  t97 = (u(i + 1,j + 1,n) - u(i,j + 1,n)) * t22
43  t103 = t58 * t61 * (t42 + t23) / 0.2D1
44  t110 = 0.1D1 / (rx(i + 1,j - 1,0,0) * rx(i + 1,j - 1,1,1) - rx(i +
45  # 1,j - 1,1,0) * rx(i + 1,j - 1,0,1))
46  t111 = t2 * t110
47  t114 = rx(i + 1,j - 1,0,0) * rx(i + 1,j - 1,1,0) + rx(i + 1,j - 1,
48  #0,1) * rx(i + 1,j - 1,1,1)
49  t118 = (u(i + 1,j - 1,n) - u(i,j - 1,n)) * t22
50  t125 = rx(i + 1,j + 1,1,0) ** 2
51  t126 = rx(i + 1,j + 1,1,1) ** 2
52  t129 = rx(i + 1,j,1,0) ** 2
53  t130 = rx(i + 1,j,1,1) ** 2
54  t132 = t6 * (t129 + t130)
55  t136 = rx(i + 1,j - 1,1,0) ** 2
56  t137 = rx(i + 1,j - 1,1,1) ** 2
57  t151 = t25 * t22
58  t157 = rx(i - 1,j,0,0) * rx(i - 1,j,1,1) - rx(i - 1,j,1,0) * rx(i
59  #- 1,j,0,1)
60  t158 = 0.1D1 / t157
61  t159 = rx(i - 1,j,0,0) ** 2
62  t160 = rx(i - 1,j,0,1) ** 2
63  t162 = t158 * (t159 + t160)
64  t163 = t18 + t162
65  t164 = t2 * t163 / 0.2D1
66  t166 = (u(i,j,n) - u(i - 1,j,n)) * t22
67  t167 = t164 * t166
68  t170 = t2 * t158
69  t173 = rx(i - 1,j,0,0) * rx(i - 1,j,1,0) + rx(i - 1,j,0,1) * rx(i
70  #- 1,j,1,1)
71  t175 = (u(i - 1,j + 1,n) - u(i - 1,j,n)) * t51
72  t177 = (u(i - 1,j,n) - u(i - 1,j - 1,n)) * t51
73  t180 = t170 * t173 * (t175 + t177) / 0.2D1
74  t183 = (t82 - t180) * t22 / 0.2D1
75  t186 = rx(i,j + 1,0,0) * rx(i,j + 1,1,1) - rx(i,j + 1,1,0) * rx(i,
76  #j + 1,0,1)
77  t187 = 0.1D1 / t186
78  t188 = t2 * t187
79  t191 = rx(i,j + 1,0,0) * rx(i,j + 1,1,0) + rx(i,j + 1,0,1) * rx(i,
80  #j + 1,1,1)
81  t193 = (u(i,j + 1,n) - u(i - 1,j + 1,n)) * t22
82  t196 = t188 * t191 * (t97 + t193) / 0.2D1
83  t199 = t72 * t75 * (t23 + t166) / 0.2D1
84  t202 = (t196 - t199) * t51 / 0.2D1
85  t205 = rx(i,j - 1,0,0) * rx(i,j - 1,1,1) - rx(i,j - 1,1,0) * rx(i,
86  #j - 1,0,1)
87  t206 = 0.1D1 / t205
88  t207 = t2 * t206
89  t210 = rx(i,j - 1,0,0) * rx(i,j - 1,1,0) + rx(i,j - 1,0,1) * rx(i,
90  #j - 1,1,1)
91  t212 = (u(i,j - 1,n) - u(i - 1,j - 1,n)) * t22
92  t215 = t207 * t210 * (t118 + t212) / 0.2D1
93  t218 = (t199 - t215) * t51 / 0.2D1
94  t219 = rx(i,j + 1,1,0) ** 2
95  t220 = rx(i,j + 1,1,1) ** 2
96  t222 = t187 * (t219 + t220)
97  t223 = rx(i,j,1,0) ** 2
98  t224 = rx(i,j,1,1) ** 2
99  t226 = t14 * (t223 + t224)
100  t227 = t222 + t226
101  t228 = t2 * t227 / 0.2D1
102  t229 = t228 * t77
103  t230 = rx(i,j - 1,1,0) ** 2
104  t231 = rx(i,j - 1,1,1) ** 2
105  t233 = t206 * (t230 + t231)
106  t234 = t226 + t233
107  t235 = t2 * t234 / 0.2D1
108  t236 = t235 * t79
109  t242 = dt * ((t24 - t167) * t22 + t85 + t183 + t202 + t218 + (t229
110  # - t236) * t51) * t13 / 0.2D1
111  t243 = ut(i,j,n) - ut(i - 1,j,n)
112  t244 = t243 * t22
113  t247 = dx * (t151 + t244) / 0.4D1
114  t254 = t24 + t20 * dt * t25 * t22 / 0.2D1 + cc * t19 * (ut(i + 1,j
115  #,n) + dt * ((t2 * (t34 * (t35 + t36) + t10) * t42 / 0.2D1 - t24) *
116  # t22 + (t2 * t34 * (rx(i + 2,j,0,0) * rx(i + 2,j,1,0) + rx(i + 2,j
117  #,0,1) * rx(i + 2,j,1,1)) * ((u(i + 2,j + 1,n) - u(i + 2,j,n)) * t5
118  #1 + (u(i + 2,j,n) - u(i + 2,j - 1,n)) * t51) / 0.2D1 - t68) * t22
119  #/ 0.2D1 + t85 + (t90 * t93 * ((u(i + 2,j + 1,n) - u(i + 1,j + 1,n)
120  #) * t22 + t97) / 0.2D1 - t103) * t51 / 0.2D1 + (t103 - t111 * t114
121  # * ((u(i + 2,j - 1,n) - u(i + 1,j - 1,n)) * t22 + t118) / 0.2D1) *
122  # t51 / 0.2D1 + (t2 * (t89 * (t125 + t126) + t132) * t63 / 0.2D1 -
123  #t2 * (t132 + t110 * (t136 + t137)) * t65 / 0.2D1) * t51) * t5 / 0.
124  #2D1 - dx * ((ut(i + 2,j,n) - ut(i + 1,j,n)) * t22 + t151) / 0.4D1
125  #- ut(i,j,n) - t242 - t247) * sqrt(0.2D1) * (t7 + t8 + t15 + t16) *
126  #* (-0.1D1 / 0.2D1) / 0.4D1
127  t255 = dt ** 2
128  t258 = t14 * t75
129  t260 = t2 * (t6 * t61 + t258) / 0.2D1
130  t267 = ut(i,j + 1,n) - ut(i,j,n)
131  t268 = t267 * t51
132  t269 = ut(i,j,n) - ut(i,j - 1,n)
133  t270 = t269 * t51
134  t275 = t260 * (t63 + t65 + t77 + t79) / 0.4D1 + t260 * dt * ((ut(i
135  # + 1,j + 1,n) - ut(i + 1,j,n)) * t51 + (ut(i + 1,j,n) - ut(i + 1,j
136  # - 1,n)) * t51 + t268 + t270) / 0.8D1
137  t285 = 0.1D1 / (rx(i - 2,j,0,0) * rx(i - 2,j,1,1) - rx(i - 2,j,1,0
138  #) * rx(i - 2,j,0,1))
139  t286 = rx(i - 2,j,0,0) ** 2
140  t287 = rx(i - 2,j,0,1) ** 2
141  t293 = (u(i - 1,j,n) - u(i - 2,j,n)) * t22
142  t314 = 0.1D1 / (rx(i - 1,j + 1,0,0) * rx(i - 1,j + 1,1,1) - rx(i -
143  # 1,j + 1,1,0) * rx(i - 1,j + 1,0,1))
144  t315 = t2 * t314
145  t318 = rx(i - 1,j + 1,0,0) * rx(i - 1,j + 1,1,0) + rx(i - 1,j + 1,
146  #0,1) * rx(i - 1,j + 1,1,1)
147  t326 = t170 * t173 * (t166 + t293) / 0.2D1
148  t333 = 0.1D1 / (rx(i - 1,j - 1,0,0) * rx(i - 1,j - 1,1,1) - rx(i -
149  # 1,j - 1,1,0) * rx(i - 1,j - 1,0,1))
150  t334 = t2 * t333
151  t337 = rx(i - 1,j - 1,0,0) * rx(i - 1,j - 1,1,0) + rx(i - 1,j - 1,
152  #0,1) * rx(i - 1,j - 1,1,1)
153  t346 = rx(i - 1,j + 1,1,0) ** 2
154  t347 = rx(i - 1,j + 1,1,1) ** 2
155  t350 = rx(i - 1,j,1,0) ** 2
156  t351 = rx(i - 1,j,1,1) ** 2
157  t353 = t158 * (t350 + t351)
158  t357 = rx(i - 1,j - 1,1,0) ** 2
159  t358 = rx(i - 1,j - 1,1,1) ** 2
160  t381 = t167 + t164 * dt * t243 * t22 / 0.2D1 + cc * t163 * (ut(i,j
161  #,n) + t242 - t247 - ut(i - 1,j,n) - dt * ((t167 - t2 * (t162 + t28
162  #5 * (t286 + t287)) * t293 / 0.2D1) * t22 + t183 + (t180 - t2 * t28
163  #5 * (rx(i - 2,j,0,0) * rx(i - 2,j,1,0) + rx(i - 2,j,0,1) * rx(i -
164  #2,j,1,1)) * ((u(i - 2,j + 1,n) - u(i - 2,j,n)) * t51 + (u(i - 2,j,
165  #n) - u(i - 2,j - 1,n)) * t51) / 0.2D1) * t22 / 0.2D1 + (t315 * t31
166  #8 * (t193 + (u(i - 1,j + 1,n) - u(i - 2,j + 1,n)) * t22) / 0.2D1 -
167  # t326) * t51 / 0.2D1 + (t326 - t334 * t337 * (t212 + (u(i - 1,j -
168  #1,n) - u(i - 2,j - 1,n)) * t22) / 0.2D1) * t51 / 0.2D1 + (t2 * (t3
169  #14 * (t346 + t347) + t353) * t175 / 0.2D1 - t2 * (t353 + t333 * (t
170  #357 + t358)) * t177 / 0.2D1) * t51) * t157 / 0.2D1 - dx * (t244 +
171  #(ut(i - 1,j,n) - ut(i - 2,j,n)) * t22) / 0.4D1) * sqrt(0.2D1) * (t
172  #15 + t16 + t159 + t160) ** (-0.1D1 / 0.2D1) / 0.4D1
173  t385 = t2 * (t258 + t158 * t173) / 0.2D1
174  t396 = t385 * (t77 + t79 + t175 + t177) / 0.4D1 + t385 * dt * (t26
175  #8 + t270 + (ut(i - 1,j + 1,n) - ut(i - 1,j,n)) * t51 + (ut(i - 1,j
176  #,n) - ut(i - 1,j - 1,n)) * t51) / 0.8D1
177  t403 = t2 * (t187 * t191 + t258) / 0.2D1
178  t414 = t403 * (t97 + t193 + t23 + t166) / 0.4D1 + t403 * dt * ((ut
179  #(i + 1,j + 1,n) - ut(i,j + 1,n)) * t22 + (ut(i,j + 1,n) - ut(i - 1
180  #,j + 1,n)) * t22 + t151 + t244) / 0.8D1
181  t421 = rx(i + 1,j + 1,0,0) ** 2
182  t422 = rx(i + 1,j + 1,0,1) ** 2
183  t425 = rx(i,j + 1,0,0) ** 2
184  t426 = rx(i,j + 1,0,1) ** 2
185  t428 = t187 * (t425 + t426)
186  t432 = rx(i - 1,j + 1,0,0) ** 2
187  t433 = rx(i - 1,j + 1,0,1) ** 2
188  t447 = (u(i,j + 2,n) - u(i,j + 1,n)) * t51
189  t450 = t188 * t191 * (t447 + t77) / 0.2D1
190  t465 = 0.1D1 / (rx(i,j + 2,0,0) * rx(i,j + 2,1,1) - rx(i,j + 2,1,0
191  #) * rx(i,j + 2,0,1))
192  t480 = rx(i,j + 2,1,0) ** 2
193  t481 = rx(i,j + 2,1,1) ** 2
194  t500 = dy * (t268 + t270) / 0.4D1
195  t507 = t229 + t228 * dt * t267 * t51 / 0.2D1 + cc * t227 * (ut(i,j
196  # + 1,n) + dt * ((t2 * (t89 * (t421 + t422) + t428) * t97 / 0.2D1 -
197  # t2 * (t428 + t314 * (t432 + t433)) * t193 / 0.2D1) * t22 + (t90 *
198  # t93 * ((u(i + 1,j + 2,n) - u(i + 1,j + 1,n)) * t51 + t63) / 0.2D1
199  # - t450) * t22 / 0.2D1 + (t450 - t315 * t318 * ((u(i - 1,j + 2,n)
200  #- u(i - 1,j + 1,n)) * t51 + t175) / 0.2D1) * t22 / 0.2D1 + (t2 * t
201  #465 * (rx(i,j + 2,0,0) * rx(i,j + 2,1,0) + rx(i,j + 2,0,1) * rx(i,
202  #j + 2,1,1)) * ((u(i + 1,j + 2,n) - u(i,j + 2,n)) * t22 + (u(i,j +
203  #2,n) - u(i - 1,j + 2,n)) * t22) / 0.2D1 - t196) * t51 / 0.2D1 + t2
204  #02 + (t2 * (t465 * (t480 + t481) + t222) * t447 / 0.2D1 - t229) *
205  #t51) * t186 / 0.2D1 - dy * ((ut(i,j + 2,n) - ut(i,j + 1,n)) * t51
206  #+ t268) / 0.4D1 - ut(i,j,n) - t242 - t500) * sqrt(0.2D1) * (t219 +
207  # t220 + t223 + t224) ** (-0.1D1 / 0.2D1) / 0.4D1
208  t511 = t2 * (t258 + t206 * t210) / 0.2D1
209  t522 = t511 * (t23 + t166 + t118 + t212) / 0.4D1 + t511 * dt * (t1
210  #51 + t244 + (ut(i + 1,j - 1,n) - ut(i,j - 1,n)) * t22 + (ut(i,j -
211  #1,n) - ut(i - 1,j - 1,n)) * t22) / 0.8D1
212  t529 = rx(i + 1,j - 1,0,0) ** 2
213  t530 = rx(i + 1,j - 1,0,1) ** 2
214  t533 = rx(i,j - 1,0,0) ** 2
215  t534 = rx(i,j - 1,0,1) ** 2
216  t536 = t206 * (t533 + t534)
217  t540 = rx(i - 1,j - 1,0,0) ** 2
218  t541 = rx(i - 1,j - 1,0,1) ** 2
219  t555 = (u(i,j - 1,n) - u(i,j - 2,n)) * t51
220  t558 = t207 * t210 * (t79 + t555) / 0.2D1
221  t573 = 0.1D1 / (rx(i,j - 2,0,0) * rx(i,j - 2,1,1) - rx(i,j - 2,1,0
222  #) * rx(i,j - 2,0,1))
223  t588 = rx(i,j - 2,1,0) ** 2
224  t589 = rx(i,j - 2,1,1) ** 2
225  t612 = t236 + t235 * dt * t269 * t51 / 0.2D1 + cc * t234 * (ut(i,j
226  #,n) + t242 - t500 - ut(i,j - 1,n) - dt * ((t2 * (t110 * (t529 + t5
227  #30) + t536) * t118 / 0.2D1 - t2 * (t536 + t333 * (t540 + t541)) *
228  #t212 / 0.2D1) * t22 + (t111 * t114 * (t65 + (u(i + 1,j - 1,n) - u(
229  #i + 1,j - 2,n)) * t51) / 0.2D1 - t558) * t22 / 0.2D1 + (t558 - t33
230  #4 * t337 * (t177 + (u(i - 1,j - 1,n) - u(i - 1,j - 2,n)) * t51) /
231  #0.2D1) * t22 / 0.2D1 + t218 + (t215 - t2 * t573 * (rx(i,j - 2,0,0)
232  # * rx(i,j - 2,1,0) + rx(i,j - 2,0,1) * rx(i,j - 2,1,1)) * ((u(i +
233  #1,j - 2,n) - u(i,j - 2,n)) * t22 + (u(i,j - 2,n) - u(i - 1,j - 2,n
234  #)) * t22) / 0.2D1) * t51 / 0.2D1 + (t236 - t2 * (t233 + t573 * (t5
235  #88 + t589)) * t555 / 0.2D1) * t51) * t205 / 0.2D1 - dy * (t270 + (
236  #ut(i,j - 1,n) - ut(i,j - 2,n)) * t51) / 0.4D1) * sqrt(0.2D1) * (t2
237  #23 + t224 + t230 + t231) ** (-0.1D1 / 0.2D1) / 0.4D1
238  cg(1,1) = u(i,j,n) + dt * ut(i,j,n) + (t254 * t255 + t275 * t255 -
239  # t381 * t255 - t396 * t255) * t13 * t22 / 0.2D1 + (t414 * t255 + t
240  #507 * t255 - t522 * t255 - t612 * t255) * t13 * t51 / 0.2D1
241  cg(2,1) = ut(i,j,n) + (t254 * dt + t275 * dt - t381 * dt - t396 *
242  #dt) * t13 * t22 + (t414 * dt + t507 * dt - t522 * dt - t612 * dt)
243  #* t13 * t51