CosmoGRaPH v0.0
math.h
1 #ifndef COSMO_UTILS_MATH_H
2 #define COSMO_UTILS_MATH_H
3 
4 #include "../cosmo_macros.h"
5 #include "../cosmo_types.h"
6 #include "../cosmo_includes.h"
7 #include "../cosmo_globals.h"
8 
9 namespace cosmo
10 {
11 
36 inline real_t KO_dissipation_Q(idx_t i, idx_t j, idx_t k, arr_t & field, real_t ko_coeff)
37 {
38  if(ko_coeff == 0)
39  return 0;
40 
41 # if STENCIL_ORDER == 2
42  real_t stencil = (
43  1.0*field[INDEX(i-2,j,k)] + 1.0*field[INDEX(i,j-2,k)] + 1.0*field[INDEX(i,j,k-2)]
44  - 4.0*field[INDEX(i-1,j,k)] - 4.0*field[INDEX(i,j-1,k)] - 4.0*field[INDEX(i,j,k-1)]
45  + 6.0*field[INDEX(i ,j,k)] + 6.0*field[INDEX(i,j ,k)] + 6.0*field[INDEX(i,j,k )]
46  - 4.0*field[INDEX(i+1,j,k)] - 4.0*field[INDEX(i,j+1,k)] - 4.0*field[INDEX(i,j,k+1)]
47  + 1.0*field[INDEX(i+2,j,k)] + 1.0*field[INDEX(i,j+2,k)] + 1.0*field[INDEX(i,j,k+2)]
48  )/pow(dx, 4.0);
49  real_t dissipation = ko_coeff*pow(dx, 3.0)/64.0*stencil;
50  return dissipation;
51 # endif
52 
53 # if STENCIL_ORDER == 8
54  real_t stencil = (
55  1.0*field[INDEX(i-5,j,k)] + 1.0*field[INDEX(i,j-5,k)] + 1.0*field[INDEX(i,j,k-5)]
56  - 10.0*field[INDEX(i-4,j,k)] - 10.0*field[INDEX(i,j-4,k)] - 10.0*field[INDEX(i,j,k-4)]
57  + 45.0*field[INDEX(i-3,j,k)] + 45.0*field[INDEX(i,j-3,k)] + 45.0*field[INDEX(i,j,k-3)]
58  - 120.0*field[INDEX(i-2,j,k)] - 120.0*field[INDEX(i,j-2,k)] - 120.0*field[INDEX(i,j,k-2)]
59  + 210.0*field[INDEX(i-1,j,k)] + 210.0*field[INDEX(i,j-1,k)] + 210.0*field[INDEX(i,j,k-1)]
60  - 252.0*field[INDEX(i ,j,k)] - 252.0*field[INDEX(i,j ,k)] - 252.0*field[INDEX(i,j,k )]
61  + 210.0*field[INDEX(i+1,j,k)] + 210.0*field[INDEX(i,j+1,k)] + 210.0*field[INDEX(i,j,k+1)]
62  - 120.0*field[INDEX(i+2,j,k)] - 120.0*field[INDEX(i,j+2,k)] - 120.0*field[INDEX(i,j,k+2)]
63  + 45.0*field[INDEX(i+3,j,k)] + 45.0*field[INDEX(i,j+3,k)] + 45.0*field[INDEX(i,j,k+3)]
64  - 10.0*field[INDEX(i+4,j,k)] - 10.0*field[INDEX(i,j+4,k)] - 10.0*field[INDEX(i,j,k+4)]
65  + 1.0*field[INDEX(i+5,j,k)] + 1.0*field[INDEX(i,j+5,k)] + 1.0*field[INDEX(i,j,k+5)]
66  )/pow(dx, 10.0);
67  real_t dissipation = -ko_coeff*pow(dx, 9.0)/1024.0*stencil;
68  return dissipation;
69 # endif
70 
71  return 0.0;
72 }
73 
74 inline real_t derivative_Odx2(idx_t i, idx_t j, idx_t k, int d,
75  arr_t & field)
76 {
77  switch (d) {
78  case 1:
79  return ((
80  - 1.0/2.0*field[INDEX(i-1,j,k)]
81  + 1.0/2.0*field[INDEX(i+1,j,k)]
82  )/dx);
83  break;
84  case 2:
85  return ((
86  - 1.0/2.0*field[INDEX(i,j-1,k)]
87  + 1.0/2.0*field[INDEX(i,j+1,k)]
88  )/dx);
89  break;
90  case 3:
91  return ((
92  - 1.0/2.0*field[INDEX(i,j,k-1)]
93  + 1.0/2.0*field[INDEX(i,j,k+1)]
94  )/dx);
95  break;
96  }
97 
98  /* XXX */
99  return 0;
100 }
101 
102 inline real_t forward_derivative_Odx2(idx_t i, idx_t j, idx_t k, int d,
103  arr_t & field)
104 {
105  switch (d) {
106  case 1:
107  return ((
108  - 1.0/2.0*field[INDEX(i+2,j,k)]
109  + 2.0*field[INDEX(i+1,j,k)]
110  - 3.0/2.0*field[INDEX(i,j,k)]
111  )/dx);
112  break;
113  case 2:
114  return ((
115  - 1.0/2.0*field[INDEX(i,j+2,k)]
116  + 2.0*field[INDEX(i,j+1,k)]
117  - 3.0/2.0*field[INDEX(i,j,k)]
118  )/dx);
119  break;
120  case 3:
121  return ((
122  - 1.0/2.0*field[INDEX(i,j,k+2)]
123  + 2.0*field[INDEX(i,j,k+1)]
124  - 3.0/2.0*field[INDEX(i,j,k)]
125  )/dx);
126  break;
127  }
128 
129  /* XXX */
130  return 0;
131 }
132 
133 
134 
135 inline real_t backward_derivative_Odx2(idx_t i, idx_t j, idx_t k, int d,
136  arr_t & field)
137 {
138  switch (d) {
139  case 1:
140  return ((
141  + 1.0/2.0*field[INDEX(i-2,j,k)]
142  - 2.0*field[INDEX(i-1,j,k)]
143  + 3.0/2.0*field[INDEX(i,j,k)]
144  )/dx);
145  break;
146  case 2:
147  return ((
148  + 1.0/2.0*field[INDEX(i,j-2,k)]
149  - 2.0*field[INDEX(i,j-1,k)]
150  + 3.0/2.0*field[INDEX(i,j,k)]
151  )/dx);
152  break;
153  case 3:
154  return ((
155  + 1.0/2.0*field[INDEX(i,j,k-2)]
156  - 2.0*field[INDEX(i,j,k-1)]
157  + 3.0/2.0*field[INDEX(i,j,k)]
158  )/dx);
159  break;
160  }
161 
162  /* XXX */
163  return 0;
164 }
165 
166 inline real_t lop_forward_derivative_Odx2(idx_t i, idx_t j, idx_t k, int d,
167  arr_t & field)
168 {
169  switch (d) {
170  case 1:
171  return ((
172  - 1.0/2.0*field[INDEX(i+2,j,k)]
173  + 2.0*field[INDEX(i+1,j,k)]
174  - 3.0/2.0*field[INDEX(i,j,k)]
175  )/dx);
176  break;
177  case 2:
178  return ((
179  - 1.0/2.0*field[INDEX(i,j+2,k)]
180  + 2.0*field[INDEX(i,j+1,k)]
181  - 3.0/2.0*field[INDEX(i,j,k)]
182  )/dx);
183  break;
184  case 3:
185  return ((
186  - 1.0/2.0*field[INDEX(i,j,k+2)]
187  + 2.0*field[INDEX(i,j,k+1)]
188  - 3.0/2.0*field[INDEX(i,j,k)]
189  )/dx);
190  break;
191  }
192 
193  /* XXX */
194  return 0;
195 }
196 
197 
198 
199 inline real_t lop_backward_derivative_Odx2(idx_t i, idx_t j, idx_t k, int d,
200  arr_t & field)
201 {
202  switch (d) {
203  case 1:
204  return ((
205  + 1.0/2.0*field[INDEX(i-2,j,k)]
206  - 2.0*field[INDEX(i-1,j,k)]
207  + 3.0/2.0*field[INDEX(i,j,k)]
208  )/dx);
209  break;
210  case 2:
211  return ((
212  + 1.0/2.0*field[INDEX(i,j-2,k)]
213  - 2.0*field[INDEX(i,j-1,k)]
214  + 3.0/2.0*field[INDEX(i,j,k)]
215  )/dx);
216  break;
217  case 3:
218  return ((
219  + 1.0/2.0*field[INDEX(i,j,k-2)]
220  - 2.0*field[INDEX(i,j,k-1)]
221  + 3.0/2.0*field[INDEX(i,j,k)]
222  )/dx);
223  break;
224  }
225 
226  /* XXX */
227  return 0;
228 }
229 
230 
231 inline real_t derivative_Odx4(idx_t i, idx_t j, idx_t k, int d,
232  arr_t & field)
233 {
234  switch (d) {
235  case 1:
236  return (
237  + 1.0/12.0*field[INDEX(i-2,j,k)]
238  - 2.0/3.0*field[INDEX(i-1,j,k)]
239  + 2.0/3.0*field[INDEX(i+1,j,k)]
240  - 1.0/12.0*field[INDEX(i+2,j,k)]
241  )/dx;
242  break;
243  case 2:
244  return (
245  + 1.0/12.0*field[INDEX(i,j-2,k)]
246  - 2.0/3.0*field[INDEX(i,j-1,k)]
247  + 2.0/3.0*field[INDEX(i,j+1,k)]
248  - 1.0/12.0*field[INDEX(i,j+2,k)]
249  )/dx;
250  break;
251  case 3:
252  return (
253  + 1.0/12.0*field[INDEX(i,j,k-2)]
254  - 2.0/3.0*field[INDEX(i,j,k-1)]
255  + 2.0/3.0*field[INDEX(i,j,k+1)]
256  - 1.0/12.0*field[INDEX(i,j,k+2)]
257  )/dx;
258  break;
259  }
260 
261  /* XXX */
262  return 0;
263 }
264 
265 inline real_t lop_forward_derivative_Odx4(idx_t i, idx_t j, idx_t k, int d,
266  arr_t & field)
267 {
268  switch (d) {
269  case 1:
270  return (
271  + 1.0/12.0*field[INDEX(i+3,j,k)]
272  - 1.0/2.0*field[INDEX(i+2,j,k)]
273  + 3.0/2.0*field[INDEX(i+1,j,k)]
274  - 5.0/6.0*field[INDEX(i,j,k)]
275  - 1.0/4.0*field[INDEX(i-1,j,k)]
276  )/dx;
277  break;
278  case 2:
279  return (
280  + 1.0/12.0*field[INDEX(i,j+3,k)]
281  - 1.0/2.0*field[INDEX(i,j+2,k)]
282  + 3.0/2.0*field[INDEX(i,j+1,k)]
283  - 5.0/6.0*field[INDEX(i,j,k)]
284  - 1.0/4.0*field[INDEX(i,j-1,k)]
285  )/dx;
286  break;
287  case 3:
288  return (
289  + 1.0/12.0*field[INDEX(i,j,k+3)]
290  - 1.0/2.0*field[INDEX(i,j,k+2)]
291  + 3.0/2.0*field[INDEX(i,j,k+1)]
292  - 5.0/6.0*field[INDEX(i,j,k)]
293  - 1.0/4.0*field[INDEX(i,j,k-1)]
294  )/dx;
295  break;
296  }
297 
298  /* XXX */
299  return 0;
300 }
301 
302 inline real_t lop_backward_derivative_Odx4(idx_t i, idx_t j, idx_t k, int d,
303  arr_t & field)
304 {
305  switch (d) {
306  case 1:
307  return (
308  - 1.0/12.0*field[INDEX(i-3,j,k)]
309  + 1.0/2.0*field[INDEX(i-2,j,k)]
310  - 3.0/2.0*field[INDEX(i-1,j,k)]
311  + 5.0/6.0*field[INDEX(i,j,k)]
312  + 1.0/4.0*field[INDEX(i+1,j,k)]
313  )/dx;
314  break;
315  case 2:
316  return (
317  - 1.0/12.0*field[INDEX(i,j-3,k)]
318  + 1.0/2.0*field[INDEX(i,j-2,k)]
319  - 3.0/2.0*field[INDEX(i,j-1,k)]
320  + 5.0/6.0*field[INDEX(i,j,k)]
321  + 1.0/4.0*field[INDEX(i,j+1,k)]
322  )/dx;
323  break;
324  case 3:
325  return (
326  - 1.0/12.0*field[INDEX(i,j,k-3)]
327  + 1.0/2.0*field[INDEX(i,j,k-2)]
328  - 3.0/2.0*field[INDEX(i,j,k-1)]
329  + 5.0/6.0*field[INDEX(i,j,k)]
330  + 1.0/4.0*field[INDEX(i,j,k+1)]
331  )/dx;
332  break;
333  }
334 
335  /* XXX */
336  return 0;
337 
338 }
339 inline real_t forward_derivative_Odx4(idx_t i, idx_t j, idx_t k, int d,
340  arr_t & field)
341 {
342  switch (d) {
343  case 1:
344  return (
345  - 1.0/4.0*field[INDEX(i+4,j,k)]
346  + 4.0/3.0*field[INDEX(i+3,j,k)]
347  - 3.0*field[INDEX(i+2,j,k)]
348  + 4.0*field[INDEX(i+1,j,k)]
349  - 25.0/12.0*field[INDEX(i,j,k)]
350  )/dx;
351  break;
352  case 2:
353  return (
354  - 1.0/4.0*field[INDEX(i,j+4,k)]
355  + 4.0/3.0*field[INDEX(i,j+3,k)]
356  - 3.0*field[INDEX(i,j+2,k)]
357  + 4.0*field[INDEX(i,j+1,k)]
358  - 25.0/12.0*field[INDEX(i,j,k)]
359  )/dx;
360  break;
361  case 3:
362  return (
363  - 1.0/4.0*field[INDEX(i,j,k+4)]
364  + 4.0/3.0*field[INDEX(i,j,k+3)]
365  - 3.0*field[INDEX(i,j,k+2)]
366  + 4.0*field[INDEX(i,j,k+1)]
367  - 25.0/12.0*field[INDEX(i,j,k)]
368  )/dx;
369  break;
370  }
371 
372  /* XXX */
373  return 0;
374 }
375 
376 inline real_t backward_derivative_Odx4(idx_t i, idx_t j, idx_t k, int d,
377  arr_t & field)
378 {
379  switch (d) {
380  case 1:
381  return (
382  + 1.0/4.0*field[INDEX(i-4,j,k)]
383  - 4.0/3.0*field[INDEX(i-3,j,k)]
384  + 3.0*field[INDEX(i-2,j,k)]
385  - 4.0*field[INDEX(i-1,j,k)]
386  + 25.0/12.0*field[INDEX(i,j,k)]
387  )/dx;
388  break;
389  case 2:
390  return (
391  + 1.0/4.0*field[INDEX(i,j-4,k)]
392  - 4.0/3.0*field[INDEX(i,j-3,k)]
393  + 3.0*field[INDEX(i,j-2,k)]
394  - 4.0*field[INDEX(i,j-1,k)]
395  + 25.0/12.0*field[INDEX(i,j,k)]
396  )/dx;
397  break;
398  case 3:
399  return (
400  + 1.0/4.0*field[INDEX(i,j,k-4)]
401  - 4.0/3.0*field[INDEX(i,j,k-3)]
402  + 3.0*field[INDEX(i,j,k-2)]
403  - 4.0*field[INDEX(i,j,k-1)]
404  + 25.0/12.0*field[INDEX(i,j,k)]
405  )/dx;
406  break;
407  }
408  return 0;
409 }
410 
411 
412 inline real_t derivative_Odx6(idx_t i, idx_t j, idx_t k, int d,
413  arr_t & field)
414 {
415  switch (d) {
416  case 1:
417  return (
418  - 1.0/60.0*field[INDEX(i-3,j,k)]
419  + 3.0/20.0*field[INDEX(i-2,j,k)]
420  - 3.0/4.0*field[INDEX(i-1,j,k)]
421  + 3.0/4.0*field[INDEX(i+1,j,k)]
422  - 3.0/20.0*field[INDEX(i+2,j,k)]
423  + 1.0/60.0*field[INDEX(i+3,j,k)]
424  )/dx;
425  break;
426  case 2:
427  return (
428  - 1.0/60.0*field[INDEX(i,j-3,k)]
429  + 3.0/20.0*field[INDEX(i,j-2,k)]
430  - 3.0/4.0*field[INDEX(i,j-1,k)]
431  + 3.0/4.0*field[INDEX(i,j+1,k)]
432  - 3.0/20.0*field[INDEX(i,j+2,k)]
433  + 1.0/60.0*field[INDEX(i,j+3,k)]
434  )/dx;
435  break;
436  case 3:
437  return (
438  - 1.0/60.0*field[INDEX(i,j,k-3)]
439  + 3.0/20.0*field[INDEX(i,j,k-2)]
440  - 3.0/4.0*field[INDEX(i,j,k-1)]
441  + 3.0/4.0*field[INDEX(i,j,k+1)]
442  - 3.0/20.0*field[INDEX(i,j,k+2)]
443  + 1.0/60.0*field[INDEX(i,j,k+3)]
444  )/dx;
445  break;
446  }
447 
448  /* XXX */
449  return 0;
450 }
451 
452 
453 inline real_t forward_derivative_Odx6(idx_t i, idx_t j, idx_t k, int d,
454  arr_t & field)
455 {
456  switch (d) {
457  case 1:
458  return (
459  - 49.0/20.0*field[INDEX(i,j,k)]
460  + 6.0*field[INDEX(i+1,j,k)]
461  - 15.0/2.0*field[INDEX(i+2,j,k)]
462  + 20.0/3.0*field[INDEX(i+3,j,k)]
463  - 15.0/4.0*field[INDEX(i+4,j,k)]
464  + 6.0/5.0*field[INDEX(i+5,j,k)]
465  - 1.0/6.0*field[INDEX(i+6,j,k)]
466  )/dx;
467  break;
468  case 2:
469  return (
470  - 49.0/20.0*field[INDEX(i,j,k)]
471  + 6.0*field[INDEX(i,j+1,k)]
472  - 15.0/2.0*field[INDEX(i,j+2,k)]
473  + 20.0/3.0*field[INDEX(i,j+3,k)]
474  - 15.0/4.0*field[INDEX(i,j+4,k)]
475  + 6.0/5.0*field[INDEX(i,j+5,k)]
476  - 1.0/6.0*field[INDEX(i,j+6,k)]
477  )/dx;
478  break;
479  case 3:
480  return (
481  - 49.0/20.0*field[INDEX(i,j,k)]
482  + 6.0*field[INDEX(i,j,k+1)]
483  - 15.0/2.0*field[INDEX(i,j,k+2)]
484  + 20.0/3.0*field[INDEX(i,j,k+3)]
485  - 15.0/4.0*field[INDEX(i,j,k+4)]
486  + 6.0/5.0*field[INDEX(i,j,k+5)]
487  - 1.0/6.0*field[INDEX(i,j,k+6)]
488  )/dx;
489  break;
490  }
491 
492  /* XXX */
493  return 0;
494 }
495 
496 inline real_t backward_derivative_Odx6(idx_t i, idx_t j, idx_t k, int d,
497  arr_t & field)
498 {
499  switch (d) {
500  case 1:
501  return (
502  + 49.0/20.0*field[INDEX(i,j,k)]
503  - 6.0*field[INDEX(i-1,j,k)]
504  + 15.0/2.0*field[INDEX(i-2,j,k)]
505  - 20.0/3.0*field[INDEX(i-3,j,k)]
506  + 15.0/4.0*field[INDEX(i-4,j,k)]
507  - 6.0/5.0*field[INDEX(i-5,j,k)]
508  + 1.0/6.0*field[INDEX(i-6,j,k)]
509  )/dx;
510  break;
511  case 2:
512  return (
513  + 49.0/20.0*field[INDEX(i,j,k)]
514  - 6.0*field[INDEX(i,j-1,k)]
515  + 15.0/2.0*field[INDEX(i,j-2,k)]
516  - 20.0/3.0*field[INDEX(i,j-3,k)]
517  + 15.0/4.0*field[INDEX(i,j-4,k)]
518  - 6.0/5.0*field[INDEX(i,j-5,k)]
519  + 1.0/6.0*field[INDEX(i,j-6,k)]
520  )/dx;
521  break;
522  case 3:
523  return (
524  + 49.0/20.0*field[INDEX(i,j,k)]
525  - 6.0*field[INDEX(i,j,k-1)]
526  + 15.0/2.0*field[INDEX(i,j,k-2)]
527  - 20.0/3.0*field[INDEX(i,j,k-3)]
528  + 15.0/4.0*field[INDEX(i,j,k-4)]
529  - 6.0/5.0*field[INDEX(i,j,k-5)]
530  + 1.0/6.0*field[INDEX(i,j,k-6)]
531  )/dx;
532  break;
533  }
534 
535  /* XXX */
536  return 0;
537 }
538 
539 inline real_t lop_forward_derivative_Odx6(idx_t i, idx_t j, idx_t k, int d,
540  arr_t & field)
541 {
542  /* XXX */
543  return 0;
544 }
545 inline real_t lop_backward_derivative_Odx6(idx_t i, idx_t j, idx_t k, int d,
546  arr_t & field)
547 {
548  return 0;
549 }
550 
551 
552 
553 
554 inline real_t derivative_Odx8(idx_t i, idx_t j, idx_t k, int d,
555  arr_t & field)
556 {
557  switch (d) {
558  case 1:
559  return (
560  ( 1.0/280.0*field[INDEX(i-4,j,k)] - 1.0/280.0*field[INDEX(i+4,j,k)] )
561  - ( 4.0/105.0*field[INDEX(i-3,j,k)] - 4.0/105.0*field[INDEX(i+3,j,k)] )
562  + ( 1.0/5.0*field[INDEX(i-2,j,k)] - 1.0/5.0*field[INDEX(i+2,j,k)] )
563  - ( 4.0/5.0*field[INDEX(i-1,j,k)] - 4.0/5.0*field[INDEX(i+1,j,k)] )
564  )/dx;
565  break;
566  case 2:
567  return (
568  ( 1.0/280.0*field[INDEX(i,j-4,k)] - 1.0/280.0*field[INDEX(i,j+4,k)] )
569  - ( 4.0/105.0*field[INDEX(i,j-3,k)] - 4.0/105.0*field[INDEX(i,j+3,k)] )
570  + ( 1.0/5.0*field[INDEX(i,j-2,k)] - 1.0/5.0*field[INDEX(i,j+2,k)] )
571  - ( 4.0/5.0*field[INDEX(i,j-1,k)] - 4.0/5.0*field[INDEX(i,j+1,k)] )
572  )/dx;
573  break;
574  case 3:
575  return (
576  ( 1.0/280.0*field[INDEX(i,j,k-4)] - 1.0/280.0*field[INDEX(i,j,k+4)] )
577  - ( 4.0/105.0*field[INDEX(i,j,k-3)] - 4.0/105.0*field[INDEX(i,j,k+3)] )
578  + ( 1.0/5.0*field[INDEX(i,j,k-2)] - 1.0/5.0*field[INDEX(i,j,k+2)] )
579  - ( 4.0/5.0*field[INDEX(i,j,k-1)] - 4.0/5.0*field[INDEX(i,j,k+1)] )
580  )/dx;
581  break;
582  }
583 
584  /* XXX */
585  return 0;
586 }
587 
588 inline real_t forward_derivative_Odx8(idx_t i, idx_t j, idx_t k, int d,
589  arr_t & field)
590 {
591  switch (d) {
592  case 1:
593  return (
594  - 761.0/280.0*field[INDEX(i,j,k)] + 8.0*field[INDEX(i+1,j,k)]
595  - 14.0*field[INDEX(i+2,j,k)] + 56.0/3.0*field[INDEX(i+3,j,k)]
596  - 35.0/2.0*field[INDEX(i+4,j,k)] + 56.0/5.0*field[INDEX(i+5,j,k)]
597  - 14.0/3.0*field[INDEX(i+6,j,k)] + 8.0/7.0*field[INDEX(i+7,j,k)]
598  - 1.0/8.0*field[INDEX(i+8,j,k)]
599  )/dx;
600  break;
601  case 2:
602  return (
603  - 761.0/280.0*field[INDEX(i,j,k)] + 8.0*field[INDEX(i,j+1,k)]
604  - 14.0*field[INDEX(i,j+2,k)] + 56.0/3.0*field[INDEX(i,j+3,k)]
605  - 35.0/2.0*field[INDEX(i,j+4,k)] + 56.0/5.0*field[INDEX(i,j+5,k)]
606  - 14.0/3.0*field[INDEX(i,j+6,k)] + 8.0/7.0*field[INDEX(i,j+7,k)]
607  - 1.0/8.0*field[INDEX(i,j+8,k)]
608  )/dx;
609  break;
610  case 3:
611  return (
612  - 761.0/280.0*field[INDEX(i,j,k)] + 8.0*field[INDEX(i,j,k+1)]
613  - 14.0*field[INDEX(i,j,k+2)] + 56.0/3.0*field[INDEX(i,j,k+3)]
614  - 35.0/2.0*field[INDEX(i,j,k+4)] + 56.0/5.0*field[INDEX(i,j,k+5)]
615  - 14.0/3.0*field[INDEX(i,j,k+6)] + 8.0/7.0*field[INDEX(i,j,k+7)]
616  - 1.0/8.0*field[INDEX(i,j,k+8)]
617  )/dx;
618  break;
619  }
620 
621  /* XXX */
622  return 0;
623 }
624 
625 inline real_t backward_derivative_Odx8(idx_t i, idx_t j, idx_t k, int d,
626  arr_t & field)
627 {
628  switch (d) {
629  case 1:
630  return (
631  + 761.0/280.0*field[INDEX(i,j,k)] - 8.0*field[INDEX(i-1,j,k)]
632  + 14.0*field[INDEX(i-2,j,k)] - 56.0/3.0*field[INDEX(i-3,j,k)]
633  + 35.0/2.0*field[INDEX(i-4,j,k)] - 56.0/5.0*field[INDEX(i-5,j,k)]
634  + 14.0/3.0*field[INDEX(i-6,j,k)] - 8.0/7.0*field[INDEX(i-7,j,k)]
635  + 1.0/8.0*field[INDEX(i-8,j,k)]
636  )/dx;
637  break;
638  case 2:
639  return (
640  + 761.0/280.0*field[INDEX(i,j,k)] - 8.0*field[INDEX(i,j-1,k)]
641  + 14.0*field[INDEX(i,j-2,k)] - 56.0/3.0*field[INDEX(i,j-3,k)]
642  + 35.0/2.0*field[INDEX(i,j-4,k)] - 56.0/5.0*field[INDEX(i,j-5,k)]
643  + 14.0/3.0*field[INDEX(i,j-6,k)] - 8.0/7.0*field[INDEX(i,j-7,k)]
644  + 1.0/8.0*field[INDEX(i,j-8,k)]
645  )/dx;
646  break;
647  case 3:
648  return (
649  + 761.0/280.0*field[INDEX(i,j,k)] - 8.0*field[INDEX(i,j,k-1)]
650  + 14.0*field[INDEX(i,j,k-2)] - 56.0/3.0*field[INDEX(i,j,k-3)]
651  + 35.0/2.0*field[INDEX(i,j,k-4)] - 56.0/5.0*field[INDEX(i,j,k-5)]
652  + 14.0/3.0*field[INDEX(i,j,k-6)] - 8.0/7.0*field[INDEX(i,j,k-7)]
653  + 1.0/8.0*field[INDEX(i,j,k-8)]
654  )/dx;
655  break;
656  }
657 
658  /* XXX */
659  return 0;
660 }
661 inline real_t lop_forward_derivative_Odx8(idx_t i, idx_t j, idx_t k, int d,
662  arr_t & field)
663 {
664 
665  /* XXX */
666  return 0;
667 }
668 inline real_t lop_backward_derivative_Odx8(idx_t i, idx_t j, idx_t k, int d,
669  arr_t & field)
670 {
671 
672  return 0;
673 }
674 
675 
676 inline real_t mixed_derivative_stencil_Odx2(idx_t i, idx_t j, idx_t k, int d1, int d2, arr_t & field)
677 {
678  if( (d1 == 1 && d2 == 2) || (d1 == 2 && d2 == 1) ) {
679  return (
680  - field[INDEX(i+1,j-1,k)] + field[INDEX(i+1,j+1,k)]
681  + field[INDEX(i-1,j-1,k)] - field[INDEX(i-1,j+1,k)]
682  )/4.0/dx/dx;
683  }
684 
685  if( (d1 == 1 && d2 == 3) || (d1 == 3 && d2 == 1) ) {
686  return (
687  - field[INDEX(i+1,j,k-1)] + field[INDEX(i+1,j,k+1)]
688  + field[INDEX(i-1,j,k-1)] - field[INDEX(i-1,j,k+1)]
689  )/4.0/dx/dx;
690  }
691 
692  if( (d1 == 3 && d2 == 2) || (d1 == 2 && d2 == 3) ) {
693  return (
694  - field[INDEX(i,j+1,k-1)] + field[INDEX(i,j+1,k+1)]
695  + field[INDEX(i,j-1,k-1)] - field[INDEX(i,j-1,k+1)]
696  )/4.0/dx/dx;
697  }
698 
699  /* XXX */
700  return 0;
701 }
702 
703 inline real_t mixed_derivative_stencil_Odx4(idx_t i, idx_t j, idx_t k, int d1, int d2, arr_t & field)
704 {
705  if( (d1 == 1 && d2 == 2) || (d1 == 2 && d2 == 1) ) {
706  return (
707  (
708  - field[INDEX(i+1,j-1,k)] + field[INDEX(i+1,j+1,k)]
709  + field[INDEX(i-1,j-1,k)] - field[INDEX(i-1,j+1,k)]
710  ) - 1.0/16.0*(
711  - field[INDEX(i+2,j-2,k)] + field[INDEX(i+2,j+2,k)]
712  + field[INDEX(i-2,j-2,k)] - field[INDEX(i-2,j+2,k)]
713  )
714  )/3.0/dx/dx;
715  }
716 
717  if( (d1 == 1 && d2 == 3) || (d1 == 3 && d2 == 1) ) {
718  return (
719  (
720  - field[INDEX(i+1,j,k-1)] + field[INDEX(i+1,j,k+1)]
721  + field[INDEX(i-1,j,k-1)] - field[INDEX(i-1,j,k+1)]
722  ) - 1.0/16.0*(
723  - field[INDEX(i+2,j,k-2)] + field[INDEX(i+2,j,k+2)]
724  + field[INDEX(i-2,j,k-2)] - field[INDEX(i-2,j,k+2)]
725  )
726  )/3.0/dx/dx;
727  }
728 
729  if( (d1 == 3 && d2 == 2) || (d1 == 2 && d2 == 3) ) {
730  return (
731  (
732  - field[INDEX(i,j+1,k-1)] + field[INDEX(i,j+1,k+1)]
733  + field[INDEX(i,j-1,k-1)] - field[INDEX(i,j-1,k+1)]
734  ) - 1.0/16.0*(
735  - field[INDEX(i,j+2,k-2)] + field[INDEX(i,j+2,k+2)]
736  + field[INDEX(i,j-2,k-2)] - field[INDEX(i,j-2,k+2)]
737  )
738  )/3.0/dx/dx;
739  }
740 
741  /* XXX */
742  return 0;
743 }
744 
745 inline real_t mixed_derivative_stencil_Odx6(idx_t i, idx_t j, idx_t k, int d1, int d2, arr_t & field)
746 {
747  if( (d1 == 1 && d2 == 2) || (d1 == 2 && d2 == 1) ) {
748  return (
749  135.0*(
750  - field[INDEX(i+1,j-1,k)] + field[INDEX(i+1,j+1,k)]
751  + field[INDEX(i-1,j-1,k)] - field[INDEX(i-1,j+1,k)]
752  ) - 27.0/2.0*(
753  - field[INDEX(i+2,j-2,k)] + field[INDEX(i+2,j+2,k)]
754  + field[INDEX(i-2,j-2,k)] - field[INDEX(i-2,j+2,k)]
755  ) + (
756  - field[INDEX(i+3,j-3,k)] + field[INDEX(i+3,j+3,k)]
757  + field[INDEX(i-3,j-3,k)] - field[INDEX(i-3,j+3,k)]
758  )
759  )/360.0/dx/dx;
760  }
761 
762  if( (d1 == 1 && d2 == 3) || (d1 == 3 && d2 == 1) ) {
763  return (
764  135.0*(
765  - field[INDEX(i+1,j,k-1)] + field[INDEX(i+1,j,k+1)]
766  + field[INDEX(i-1,j,k-1)] - field[INDEX(i-1,j,k+1)]
767  ) - 27.0/2.0*(
768  - field[INDEX(i+2,j,k-2)] + field[INDEX(i+2,j,k+2)]
769  + field[INDEX(i-2,j,k-2)] - field[INDEX(i-2,j,k+2)]
770  ) + (
771  - field[INDEX(i+3,j,k-3)] + field[INDEX(i+3,j,k+3)]
772  + field[INDEX(i-3,j,k-3)] - field[INDEX(i-3,j,k+3)]
773  )
774  )/360.0/dx/dx;
775  }
776 
777  if( (d1 == 3 && d2 == 2) || (d1 == 2 && d2 == 3) ) {
778  return (
779  135.0*(
780  - field[INDEX(i,j+1,k-1)] + field[INDEX(i,j+1,k+1)]
781  + field[INDEX(i,j-1,k-1)] - field[INDEX(i,j-1,k+1)]
782  ) - 27.0/2.0*(
783  - field[INDEX(i,j+2,k-2)] + field[INDEX(i,j+2,k+2)]
784  + field[INDEX(i,j-2,k-2)] - field[INDEX(i,j-2,k+2)]
785  ) + (
786  - field[INDEX(i,j+3,k-3)] + field[INDEX(i,j+3,k+3)]
787  + field[INDEX(i,j-3,k-3)] - field[INDEX(i,j-3,k+3)]
788  )
789  )/360.0/dx/dx;
790  }
791 
792  /* XXX */
793  return 0;
794 }
795 
796 inline real_t mixed_derivative_stencil_Odx8(idx_t i, idx_t j, idx_t k, int d1,
797  int d2, arr_t & field)
798 {
799  if( (d1 == 1 && d2 == 2) || (d1 == 2 && d2 == 1) ) {
800  return (
801  2.0/5.0*(
802  - field[INDEX(i+1,j-1,k)] + field[INDEX(i+1,j+1,k)]
803  + field[INDEX(i-1,j-1,k)] - field[INDEX(i-1,j+1,k)]
804  ) - 1.0/20.0*(
805  - field[INDEX(i+2,j-2,k)] + field[INDEX(i+2,j+2,k)]
806  + field[INDEX(i-2,j-2,k)] - field[INDEX(i-2,j+2,k)]
807  ) + 2.0/315.0*(
808  - field[INDEX(i+3,j-3,k)] + field[INDEX(i+3,j+3,k)]
809  + field[INDEX(i-3,j-3,k)] - field[INDEX(i-3,j+3,k)]
810  ) - 1.0/2240.0*(
811  - field[INDEX(i+4,j-4,k)] + field[INDEX(i+4,j+4,k)]
812  + field[INDEX(i-4,j-4,k)] - field[INDEX(i-4,j+4,k)]
813  )
814  )/dx/dx;
815  }
816 
817  if( (d1 == 1 && d2 == 3) || (d1 == 3 && d2 == 1) ) {
818  return (
819  2.0/5.0*(
820  - field[INDEX(i+1,j,k-1)] + field[INDEX(i+1,j,k+1)]
821  + field[INDEX(i-1,j,k-1)] - field[INDEX(i-1,j,k+1)]
822  ) - 1.0/20.0*(
823  - field[INDEX(i+2,j,k-2)] + field[INDEX(i+2,j,k+2)]
824  + field[INDEX(i-2,j,k-2)] - field[INDEX(i-2,j,k+2)]
825  ) + 2.0/315.0*(
826  - field[INDEX(i+3,j,k-3)] + field[INDEX(i+3,j,k+3)]
827  + field[INDEX(i-3,j,k-3)] - field[INDEX(i-3,j,k+3)]
828  ) - 1.0/2240.0*(
829  - field[INDEX(i+4,j,k-4)] + field[INDEX(i+4,j,k+4)]
830  + field[INDEX(i-4,j,k-4)] - field[INDEX(i-4,j,k+4)]
831  )
832  )/dx/dx;
833  }
834 
835  if( (d1 == 3 && d2 == 2) || (d1 == 2 && d2 == 3) ) {
836  return (
837  2.0/5.0*(
838  - field[INDEX(i,j+1,k-1)] + field[INDEX(i,j+1,k+1)]
839  + field[INDEX(i,j-1,k-1)] - field[INDEX(i,j-1,k+1)]
840  ) - 1.0/20.0*(
841  - field[INDEX(i,j+2,k-2)] + field[INDEX(i,j+2,k+2)]
842  + field[INDEX(i,j-2,k-2)] - field[INDEX(i,j-2,k+2)]
843  ) + 2.0/315.0*(
844  - field[INDEX(i,j+3,k-3)] + field[INDEX(i,j+3,k+3)]
845  + field[INDEX(i,j-3,k-3)] - field[INDEX(i,j-3,k+3)]
846  ) - 1.0/2240.0*(
847  - field[INDEX(i,j+4,k-4)] + field[INDEX(i,j+4,k+4)]
848  + field[INDEX(i,j-4,k-4)] - field[INDEX(i,j-4,k+4)]
849  )
850  )/dx/dx;
851  }
852 
853  /* XXX */
854  return 0;
855 }
856 
857 inline real_t double_derivative_stencil_Odx2(idx_t i, idx_t j, idx_t k, int d,
858  arr_t & field)
859 {
860  switch (d) {
861  case 1:
862  return (
863  field[INDEX(i-1,j,k)]
864  - 2.0*field[INDEX(i-0,j,k)]
865  + field[INDEX(i+1,j,k)]
866  )/dx/dx;
867  break;
868  case 2:
869  return (
870  field[INDEX(i,j-1,k)]
871  - 2.0*field[INDEX(i,j-0,k)]
872  + field[INDEX(i,j+1,k)]
873  )/dx/dx;
874  break;
875  case 3:
876  return (
877  field[INDEX(i,j,k-1)]
878  - 2.0*field[INDEX(i,j,k-0)]
879  + field[INDEX(i,j,k+1)]
880  )/dx/dx;
881  break;
882  }
883 
884  /* XXX */
885  return 0;
886 }
887 
888 inline real_t forward_double_derivative_stencil_Odx2(idx_t i, idx_t j, idx_t k, int d,
889  arr_t & field)
890 {
891  switch (d) {
892  case 1:
893  return (
894  2.0*field[INDEX(i,j,k)]
895  - 5.0*field[INDEX(i+1,j,k)]
896  + 4.0*field[INDEX(i+2,j,k)]
897  - 1.0*field[INDEX(i+3,j,k)]
898  )/dx/dx;
899  break;
900  case 2:
901  return (
902  2.0*field[INDEX(i,j,k)]
903  - 5.0*field[INDEX(i,j+1,k)]
904  + 4.0*field[INDEX(i,j+2,k)]
905  - 1.0*field[INDEX(i,j+3,k)]
906  )/dx/dx;
907  break;
908  case 3:
909  return (
910  2.0*field[INDEX(i,j,k)]
911  - 5.0*field[INDEX(i,j,k+1)]
912  + 4.0*field[INDEX(i,j,k+2)]
913  - 1.0*field[INDEX(i,j,k+3)]
914  )/dx/dx;
915  break;
916  }
917 
918  /* XXX */
919  return 0;
920 }
921 
922 
923 inline real_t backward_double_derivative_stencil_Odx2(idx_t i, idx_t j, idx_t k, int d,
924  arr_t & field)
925 {
926  switch (d) {
927  case 1:
928  return (
929  2.0*field[INDEX(i,j,k)]
930  - 5.0*field[INDEX(i-1,j,k)]
931  + 4.0*field[INDEX(i-2,j,k)]
932  - 1.0*field[INDEX(i-3,j,k)]
933  )/dx/dx;
934  break;
935  case 2:
936  return (
937  2.0*field[INDEX(i,j,k)]
938  - 5.0*field[INDEX(i,j-1,k)]
939  + 4.0*field[INDEX(i,j-2,k)]
940  - 1.0*field[INDEX(i,j-3,k)]
941  )/dx/dx;
942  break;
943  case 3:
944  return (
945  2.0*field[INDEX(i,j,k)]
946  - 5.0*field[INDEX(i,j,k-1)]
947  + 4.0*field[INDEX(i,j,k-2)]
948  - 1.0*field[INDEX(i,j,k-3)]
949  )/dx/dx;
950  break;
951  }
952 
953  /* XXX */
954  return 0;
955 }
956 
957 
958 inline real_t double_derivative_stencil_Odx4(idx_t i, idx_t j, idx_t k, int d,
959  arr_t & field)
960 {
961  switch (d) {
962  case 1:
963  return (
964  - 1.0/12.0*field[INDEX(i-2,j,k)]
965  + 4.0/3.0*field[INDEX(i-1,j,k)]
966  - 5.0/2.0*field[INDEX(i-0,j,k)]
967  + 4.0/3.0*field[INDEX(i+1,j,k)]
968  - 1.0/12.0*field[INDEX(i+2,j,k)]
969  )/dx/dx;
970  break;
971  case 2:
972  return (
973  - 1.0/12.0*field[INDEX(i,j-2,k)]
974  + 4.0/3.0*field[INDEX(i,j-1,k)]
975  - 5.0/2.0*field[INDEX(i,j-0,k)]
976  + 4.0/3.0*field[INDEX(i,j+1,k)]
977  - 1.0/12.0*field[INDEX(i,j+2,k)]
978  )/dx/dx;
979  break;
980  case 3:
981  return (
982  - 1.0/12.0*field[INDEX(i,j,k-2)]
983  + 4.0/3.0*field[INDEX(i,j,k-1)]
984  - 5.0/2.0*field[INDEX(i,j,k-0)]
985  + 4.0/3.0*field[INDEX(i,j,k+1)]
986  - 1.0/12.0*field[INDEX(i,j,k+2)]
987  )/dx/dx;
988  break;
989  }
990 
991  /* XXX */
992  return 0;
993 }
994 
995 inline real_t double_derivative_stencil_Odx6(idx_t i, idx_t j, idx_t k, int d,
996  arr_t & field)
997 {
998  switch (d) {
999  case 1:
1000  return (
1001  1.0/90.0*field[INDEX(i-3,j,k)]
1002  - 3.0/20.0*field[INDEX(i-2,j,k)]
1003  + 3.0/2.0*field[INDEX(i-1,j,k)]
1004  - 49.0/18.0*field[INDEX(i-0,j,k)]
1005  + 3.0/2.0*field[INDEX(i+1,j,k)]
1006  - 3.0/20.0*field[INDEX(i+2,j,k)]
1007  + 1.0/90.0*field[INDEX(i+3,j,k)]
1008  )/dx/dx;
1009  break;
1010  case 2:
1011  return (
1012  1.0/90.0*field[INDEX(i,j-3,k)]
1013  - 3.0/20.0*field[INDEX(i,j-2,k)]
1014  + 3.0/2.0*field[INDEX(i,j-1,k)]
1015  - 49.0/18.0*field[INDEX(i,j-0,k)]
1016  + 3.0/2.0*field[INDEX(i,j+1,k)]
1017  - 3.0/20.0*field[INDEX(i,j+2,k)]
1018  + 1.0/90.0*field[INDEX(i,j+3,k)]
1019  )/dx/dx;
1020  break;
1021  case 3:
1022  return (
1023  1.0/90.0*field[INDEX(i,j,k-3)]
1024  - 3.0/20.0*field[INDEX(i,j,k-2)]
1025  + 3.0/2.0*field[INDEX(i,j,k-1)]
1026  - 49.0/18.0*field[INDEX(i,j,k-0)]
1027  + 3.0/2.0*field[INDEX(i,j,k+1)]
1028  - 3.0/20.0*field[INDEX(i,j,k+2)]
1029  + 1.0/90.0*field[INDEX(i,j,k+3)]
1030  )/dx/dx;
1031  break;
1032  }
1033 
1034  /* XXX */
1035  return 0;
1036 }
1037 
1038 inline real_t double_derivative_stencil_Odx8(idx_t i, idx_t j, idx_t k, int d,
1039  arr_t & field)
1040 {
1041  switch (d) {
1042  case 1:
1043  return (
1044  - 1.0/560.0*field[INDEX(i-4,j,k)]
1045  + 8.0/315.0*field[INDEX(i-3,j,k)]
1046  - 1.0/5.0*field[INDEX(i-2,j,k)]
1047  + 8.0/5.0*field[INDEX(i-1,j,k)]
1048  - 205.0/72.0*field[INDEX(i-0,j,k)]
1049  + 8.0/5.0*field[INDEX(i+1,j,k)]
1050  - 1.0/5.0*field[INDEX(i+2,j,k)]
1051  + 8.0/315.0*field[INDEX(i+3,j,k)]
1052  - 1.0/560.0*field[INDEX(i+4,j,k)]
1053  )/dx/dx;
1054  break;
1055  case 2:
1056  return (
1057  - 1.0/560.0*field[INDEX(i,j-4,k)]
1058  + 8.0/315.0*field[INDEX(i,j-3,k)]
1059  - 1.0/5.0*field[INDEX(i,j-2,k)]
1060  + 8.0/5.0*field[INDEX(i,j-1,k)]
1061  - 205.0/72.0*field[INDEX(i,j-0,k)]
1062  + 8.0/5.0*field[INDEX(i,j+1,k)]
1063  - 1.0/5.0*field[INDEX(i,j+2,k)]
1064  + 8.0/315.0*field[INDEX(i,j+3,k)]
1065  - 1.0/560.0*field[INDEX(i,j+4,k)]
1066  )/dx/dx;
1067  break;
1068  case 3:
1069  return (
1070  - 1.0/560.0*field[INDEX(i,j,k-4)]
1071  + 8.0/315.0*field[INDEX(i,j,k-3)]
1072  - 1.0/5.0*field[INDEX(i,j,k-2)]
1073  + 8.0/5.0*field[INDEX(i,j,k-1)]
1074  - 205.0/72.0*field[INDEX(i,j,k-0)]
1075  + 8.0/5.0*field[INDEX(i,j,k+1)]
1076  - 1.0/5.0*field[INDEX(i,j,k+2)]
1077  + 8.0/315.0*field[INDEX(i,j,k+3)]
1078  - 1.0/560.0*field[INDEX(i,j,k+4)]
1079  )/dx/dx;
1080  break;
1081  }
1082 
1083  /* XXX */
1084  return 0;
1085 }
1086 
1087 inline real_t forward_dissipation_stencil_Odx2(idx_t i, idx_t j, idx_t k, int d,
1088  arr_t & field)
1089 {
1090  switch (d) {
1091  case 1:
1092  return -1.0/2.0*dx*
1093  forward_double_derivative_stencil_Odx2(i,j,k,d,field);
1094  break;
1095  case 2:
1096  return -1.0/2.0*dx*
1097  forward_double_derivative_stencil_Odx2(i,j,k,d,field);
1098  break;
1099  case 3:
1100  return -1.0/2.0*dx*
1101  forward_double_derivative_stencil_Odx2(i,j,k,d,field);
1102  break;
1103  }
1104 
1105  /* XXX */
1106  return 0;
1107 }
1108 
1109 
1110 inline real_t backward_dissipation_stencil_Odx2(idx_t i, idx_t j, idx_t k, int d,
1111  arr_t & field)
1112 {
1113  switch (d) {
1114  case 1:
1115  return +1.0/2.0*dx*
1116  backward_double_derivative_stencil_Odx2(i,j,k,d,field);
1117  break;
1118  case 2:
1119  return +1.0/2.0*dx*
1120  backward_double_derivative_stencil_Odx2(i,j,k,d,field);
1121  break;
1122  case 3:
1123  return +1.0/2.0*dx*
1124  backward_double_derivative_stencil_Odx2(i,j,k,d,field);
1125  break;
1126  }
1127 
1128  /* XXX */
1129  return 0;
1130 }
1131 
1132 inline real_t forward_dissipation_stencil_Odx4(idx_t i, idx_t j, idx_t k, int d,
1133  arr_t & field)
1134 {
1135  /* XXX */
1136  return 0;
1137 }
1138 
1139 
1140 inline real_t backward_dissipation_stencil_Odx4(idx_t i, idx_t j, idx_t k, int d,
1141  arr_t & field)
1142 {
1143  /* XXX */
1144  return 0;
1145 }
1146 
1147 
1148 inline real_t forward_dissipation_stencil_Odx6(idx_t i, idx_t j, idx_t k, int d,
1149  arr_t & field)
1150 {
1151  /* XXX */
1152  return 0;
1153 }
1154 
1155 
1156 inline real_t backward_dissipation_stencil_Odx6(idx_t i, idx_t j, idx_t k, int d,
1157  arr_t & field)
1158 {
1159  /* XXX */
1160  return 0;
1161 }
1162 
1163 
1164 inline real_t forward_dissipation_stencil_Odx8(idx_t i, idx_t j, idx_t k, int d,
1165  arr_t & field)
1166 {
1167  /* XXX */
1168  return 0;
1169 }
1170 
1171 
1172 inline real_t backward_dissipation_stencil_Odx8(idx_t i, idx_t j, idx_t k, int d,
1173  arr_t & field)
1174 {
1175  /* XXX */
1176  return 0;
1177 }
1178 
1179 
1180 inline real_t forward_dissipation(idx_t i, idx_t j, idx_t k, int d,
1181  arr_t & field)
1182 {
1183  return STENCIL_ORDER_FUNCTION(forward_dissipation_stencil_Odx)(i, j, k, d, field);
1184 }
1185 
1186 
1187 inline real_t backward_dissipation(idx_t i, idx_t j, idx_t k, int d,
1188  arr_t & field)
1189 {
1190  return STENCIL_ORDER_FUNCTION(backward_dissipation_stencil_Odx)(i, j, k, d, field);
1191 }
1192 
1193 
1205 inline real_t derivative(idx_t i, idx_t j, idx_t k, int d,
1206  arr_t & field)
1207 {
1208 
1209 #if NY == 1
1210  if(d == 2) return 0;
1211 #endif
1212 
1213 #if NZ == 1
1214  if(d == 3) return 0;
1215 #endif
1216 
1217  return STENCIL_ORDER_FUNCTION(derivative_Odx)(i, j, k, d, field);
1218 }
1219 
1220 inline real_t lop_forward_derivative(idx_t i, idx_t j, idx_t k, int d,
1221  arr_t & field)
1222 {
1223  return STENCIL_ORDER_FUNCTION(lop_forward_derivative_Odx)(i, j, k, d, field)
1224  + STENCIL_ORDER_FUNCTION(forward_dissipation_stencil_Odx)(i, j, k, d, field);
1225 }
1226 
1227 inline real_t lop_backward_derivative(idx_t i, idx_t j, idx_t k, int d,
1228  arr_t & field)
1229 {
1230  return STENCIL_ORDER_FUNCTION(lop_backward_derivative_Odx)(i, j, k, d, field)
1231  + STENCIL_ORDER_FUNCTION(backward_dissipation_stencil_Odx)(i, j, k, d, field);
1232 }
1233 inline real_t forward_derivative(idx_t i, idx_t j, idx_t k, int d,
1234  arr_t & field)
1235 {
1236  return STENCIL_ORDER_FUNCTION(forward_derivative_Odx)(i, j, k, d, field)
1237  + STENCIL_ORDER_FUNCTION(forward_dissipation_stencil_Odx)(i, j, k, d, field);
1238 }
1239 
1240 inline real_t backward_derivative(idx_t i, idx_t j, idx_t k, int d,
1241  arr_t & field)
1242 {
1243  return STENCIL_ORDER_FUNCTION(backward_derivative_Odx)(i, j, k, d, field)
1244  + STENCIL_ORDER_FUNCTION(backward_dissipation_stencil_Odx)(i, j, k, d, field);
1245 }
1246 
1247 inline real_t upwind_derivative(idx_t i, idx_t j, idx_t k, int d,
1248  arr_t & field, real_t c)
1249 {
1250  if( c > 0) return c * lop_forward_derivative(i,j,k,d,field);
1251  else return c * lop_backward_derivative(i,j,k,d,field);
1252 }
1253 
1254 
1267 inline real_t mixed_derivative_stencil(idx_t i, idx_t j, idx_t k, int d1, int d2, arr_t & field)
1268 {
1269 
1270 #if NY == 1
1271  if(d1 == 2 || d2 == 2) return 0;
1272 #endif
1273 
1274 #if NZ == 1
1275  if(d1 == 3 || d2 == 3) return 0;
1276 #endif
1277 
1278  return STENCIL_ORDER_FUNCTION(mixed_derivative_stencil_Odx)(i, j, k, d1, d2, field);
1279 }
1280 
1292 inline real_t double_derivative_stencil(idx_t i, idx_t j, idx_t k, int d,
1293  arr_t & field)
1294 {
1295 
1296 #if NY == 1
1297  if(d == 2) return 0;
1298 #endif
1299 
1300 #if NZ == 1
1301  if(d == 3) return 0;
1302 #endif
1303 
1304  return STENCIL_ORDER_FUNCTION(double_derivative_stencil_Odx)(i, j, k, d, field);
1305 }
1306 
1319 inline real_t double_derivative(idx_t i, idx_t j, idx_t k, int d1, int d2,
1320  arr_t & field)
1321 {
1322  if(d1 == d2) {
1323  return double_derivative_stencil(i, j, k, d1, field);
1324  } else {
1325  return mixed_derivative_stencil(i, j, k, d1, d2, field);
1326  }
1327 
1328  /* XXX */
1329  return 0;
1330 }
1331 
1344 inline real_t laplacian(idx_t i, idx_t j, idx_t k, arr_t & field)
1345 {
1346  return (
1347  double_derivative(i, j, k, 1, 1, field)
1348  + double_derivative(i, j, k, 2, 2, field)
1349  + double_derivative(i, j, k, 3, 3, field)
1350  );
1351 }
1352 
1359 inline real_t average(arr_t & field)
1360 {
1361  // note this may have poor precision for large datasets
1362  real_t sum = 0.0;
1363  idx_t i=0, j=0, k=0;
1364 
1365  #pragma omp parallel for default(shared) private(i, j, k) reduction(+:sum)
1366  LOOP3(i, j, k)
1367  {
1368  sum += field[NP_INDEX(i,j,k)];
1369  }
1370  return sum/POINTS;
1371 }
1372 
1373 
1382 inline real_t volume_average(arr_t & DIFFphi, real_t phi_FRW)
1383 {
1384  real_t sum = 0.0;
1385  idx_t i=0, j=0, k=0;
1386 
1387  #pragma omp parallel for default(shared) private(i, j, k) reduction(+:sum)
1388  LOOP3(i, j, k)
1389  {
1390  sum += exp(6.0*(DIFFphi[NP_INDEX(i,j,k)] + phi_FRW));
1391  }
1392  return sum/POINTS;
1393 }
1394 
1403 inline real_t conformal_average(arr_t & field, arr_t & DIFFphi, real_t phi_FRW)
1404 {
1405  real_t sum = 0.0;
1406  idx_t i=0, j=0, k=0;
1407 
1408  #pragma omp parallel for default(shared) private(i, j, k) reduction(+:sum)
1409  LOOP3(i, j, k)
1410  {
1411  sum += exp(6.0*(DIFFphi[NP_INDEX(i,j,k)] + phi_FRW))*field[NP_INDEX(i,j,k)];
1412  }
1413  real_t vol = volume_average(DIFFphi, phi_FRW);
1414  return sum/POINTS/vol;
1415 }
1416 
1424 inline real_t standard_deviation(arr_t & field, real_t avg)
1425 {
1426  // note this may have poor precision for large datasets
1427  idx_t i=0, j=0, k=0;
1428  real_t sum = 0.0;
1429  #pragma omp parallel for default(shared) private(i, j, k) reduction(+:sum)
1430  LOOP3(i, j, k)
1431  {
1432  sum += pw2(avg - field[NP_INDEX(i,j,k)]);
1433  }
1434  return sqrt(sum/(POINTS-1));
1435 }
1436 
1443 inline real_t standard_deviation(arr_t & field)
1444 {
1445  real_t avg = average(field);
1446  return standard_deviation(field, avg);
1447 }
1448 
1458 inline real_t conformal_standard_deviation(arr_t & field, arr_t & DIFFphi, real_t phi_FRW, real_t avg)
1459 {
1460  real_t sum = 0.0;
1461  idx_t i=0, j=0, k=0;
1462 
1463  #pragma omp parallel for default(shared) private(i, j, k) reduction(+:sum)
1464  LOOP3(i, j, k)
1465  {
1466  sum += exp(6.0*(DIFFphi[NP_INDEX(i,j,k)] + phi_FRW))*pw2(avg - field[NP_INDEX(i,j,k)]);
1467  }
1468  real_t vol = volume_average(DIFFphi, phi_FRW);
1469  return sqrt(sum/(POINTS-1)/vol);
1470 }
1471 
1480 inline real_t conformal_standard_deviation(arr_t & field, arr_t & DIFFphi, real_t phi_FRW)
1481 {
1482  real_t avg = conformal_average(field, DIFFphi, phi_FRW);
1483  return conformal_standard_deviation(field, DIFFphi, phi_FRW, avg);
1484 }
1485 
1489 inline real_t max(arr_t & field)
1490 {
1491  idx_t i=0, j=0, k=0;
1492  real_t max_val = field[0];
1493 #pragma omp parallel for default(shared) private(i, j, k) reduction(max:max_val)
1494  LOOP3(i, j, k)
1495  {
1496  if(field[INDEX(i,j,k)] > max_val) {
1497  max_val = field[INDEX(i,j,k)];
1498  }
1499  }
1500  return max_val;
1501 }
1502 
1506 inline real_t min(arr_t & field)
1507 {
1508  idx_t i=0, j=0, k=0;
1509  real_t min_val = field[0];
1510 #pragma omp parallel for default(shared) private(i, j, k) reduction(min:min_val)
1511  LOOP3(i, j, k)
1512  {
1513  if(field[INDEX(i,j,k)] < min_val) {
1514  min_val = field[INDEX(i,j,k)];
1515  }
1516  }
1517  return min_val;
1518 }
1519 
1524 inline idx_t numNaNs(arr_t & field)
1525 {
1526  idx_t i=0, j=0, k=0;
1527  idx_t NaNs = 0;
1528 
1529  #pragma omp parallel for default(shared) private(i, j, k) reduction(+:NaNs)
1530  LOOP3(i,j,k)
1531  {
1532  real_t val = field[NP_INDEX(i,j,k)];
1533  union { float val; uint32_t x; } u = { (float) val };
1534  if((u.x << 1) > 0xff000000u)
1535  {
1536  NaNs += 1;
1537  }
1538  }
1539 
1540  return NaNs;
1541 }
1542 
1543 
1544 inline idx_t idx_t_mod(idx_t n, idx_t d)
1545 {
1546  idx_t mod = n % d;
1547  if(mod < 0)
1548  mod += d;
1549  return mod;
1550 }
1551 
1552 inline real_t real_t_mod(real_t n, real_t d)
1553 {
1554  // is there a (small) chance of this not actually working due to roundoff?
1555  return n - d*std::floor(n/d);
1556 }
1557 
1558 
1559 
1560 
1561 /******************************************************************************************/
1562 
1563 
1564 inline real_t derivative_Odx2(idx_t i, idx_t j, idx_t k, idx_t nx, idx_t ny, idx_t nz, int d, arr_t & field)
1565 {
1566  real_t dx = H_LEN_FRAC / nx, dy = H_LEN_FRAC / ny, dz = H_LEN_FRAC / nz;
1567  switch (d) {
1568  case 1:
1569  return ((
1570  - 1.0/2.0*field[H_INDEX(i-1,j,k,nx,ny,nz)]
1571  + 1.0/2.0*field[H_INDEX(i+1,j,k,nx,ny,nz)]
1572  )/dx);
1573  break;
1574  case 2:
1575  return ((
1576  - 1.0/2.0*field[H_INDEX(i,j-1,k,nx,ny,nz)]
1577  + 1.0/2.0*field[H_INDEX(i,j+1,k,nx,ny,nz)]
1578  )/dy);
1579  break;
1580  case 3:
1581  return ((
1582  - 1.0/2.0*field[H_INDEX(i,j,k-1,nx,ny,nz)]
1583  + 1.0/2.0*field[H_INDEX(i,j,k+1,nx,ny,nz)]
1584  )/dz);
1585  break;
1586  }
1587 
1588  /* XXX */
1589  return 0;
1590 }
1591 
1592 inline real_t derivative_Odx4(idx_t i, idx_t j, idx_t k, idx_t nx, idx_t ny, idx_t nz, int d, arr_t & field)
1593 {
1594  real_t dx = H_LEN_FRAC / nx, dy = H_LEN_FRAC / ny, dz = H_LEN_FRAC / nz;
1595  switch (d) {
1596  case 1:
1597  return (
1598  + 1.0/12.0*field[H_INDEX(i-2,j,k,nx,ny,nz)]
1599  - 2.0/3.0*field[H_INDEX(i-1,j,k,nx,ny,nz)]
1600  + 2.0/3.0*field[H_INDEX(i+1,j,k,nx,ny,nz)]
1601  - 1.0/12.0*field[H_INDEX(i+2,j,k,nx,ny,nz)]
1602  )/dx;
1603  break;
1604  case 2:
1605  return (
1606  + 1.0/12.0*field[H_INDEX(i,j-2,k,nx,ny,nz)]
1607  - 2.0/3.0*field[H_INDEX(i,j-1,k,nx,ny,nz)]
1608  + 2.0/3.0*field[H_INDEX(i,j+1,k,nx,ny,nz)]
1609  - 1.0/12.0*field[H_INDEX(i,j+2,k,nx,ny,nz)]
1610  )/dy;
1611  break;
1612  case 3:
1613  return (
1614  + 1.0/12.0*field[H_INDEX(i,j,k-2,nx,ny,nz)]
1615  - 2.0/3.0*field[H_INDEX(i,j,k-1,nx,ny,nz)]
1616  + 2.0/3.0*field[H_INDEX(i,j,k+1,nx,ny,nz)]
1617  - 1.0/12.0*field[H_INDEX(i,j,k+2,nx,ny,nz)]
1618  )/dz;
1619  break;
1620  }
1621 
1622  /* XXX */
1623  return 0;
1624 }
1625 
1626 inline real_t derivative_Odx6(idx_t i, idx_t j, idx_t k,idx_t nx, idx_t ny, idx_t nz, int d,
1627 arr_t & field)
1628 {
1629  real_t dx = H_LEN_FRAC / nx, dy = H_LEN_FRAC / ny, dz = H_LEN_FRAC / nz;
1630  switch (d) {
1631  case 1:
1632  return (
1633  - 1.0/60.0*field[H_INDEX(i-3,j,k,nx,ny,nz)]
1634  + 3.0/20.0*field[H_INDEX(i-2,j,k,nx,ny,nz)]
1635  - 3.0/4.0*field[H_INDEX(i-1,j,k,nx,ny,nz)]
1636  + 3.0/4.0*field[H_INDEX(i+1,j,k,nx,ny,nz)]
1637  - 3.0/20.0*field[H_INDEX(i+2,j,k,nx,ny,nz)]
1638  + 1.0/60.0*field[H_INDEX(i+3,j,k,nx,ny,nz)]
1639  )/dx;
1640  break;
1641  case 2:
1642  return (
1643  - 1.0/60.0*field[H_INDEX(i,j-3,k,nx,ny,nz)]
1644  + 3.0/20.0*field[H_INDEX(i,j-2,k,nx,ny,nz)]
1645  - 3.0/4.0*field[H_INDEX(i,j-1,k,nx,ny,nz)]
1646  + 3.0/4.0*field[H_INDEX(i,j+1,k,nx,ny,nz)]
1647  - 3.0/20.0*field[H_INDEX(i,j+2,k,nx,ny,nz)]
1648  + 1.0/60.0*field[H_INDEX(i,j+3,k,nx,ny,nz)]
1649  )/dy;
1650  break;
1651  case 3:
1652  return (
1653  - 1.0/60.0*field[H_INDEX(i,j,k-3,nx,ny,nz)]
1654  + 3.0/20.0*field[H_INDEX(i,j,k-2,nx,ny,nz)]
1655  - 3.0/4.0*field[H_INDEX(i,j,k-1,nx,ny,nz)]
1656  + 3.0/4.0*field[H_INDEX(i,j,k+1,nx,ny,nz)]
1657  - 3.0/20.0*field[H_INDEX(i,j,k+2,nx,ny,nz)]
1658  + 1.0/60.0*field[H_INDEX(i,j,k+3,nx,ny,nz)]
1659  )/dz;
1660  break;
1661  }
1662 
1663  /* XXX */
1664  return 0;
1665 }
1666 
1667 inline real_t derivative_Odx8(idx_t i, idx_t j, idx_t k,idx_t nx, idx_t ny, idx_t nz, int d,arr_t & field)
1668 {
1669  real_t dx = H_LEN_FRAC / nx, dy = H_LEN_FRAC / ny, dz = H_LEN_FRAC / nz;
1670  switch (d) {
1671  case 1:
1672  return (
1673  ( 1.0/280.0*field[H_INDEX(i-4,j,k,nx,ny,nz)] - 1.0/280.0*field[H_INDEX(i+4,j,k,nx,ny,nz)] )
1674  - ( 4.0/105.0*field[H_INDEX(i-3,j,k,nx,ny,nz)] - 4.0/105.0*field[H_INDEX(i+3,j,k,nx,ny,nz)] )
1675  + ( 1.0/5.0*field[H_INDEX(i-2,j,k,nx,ny,nz)] - 1.0/5.0*field[H_INDEX(i+2,j,k,nx,ny,nz)] )
1676  - ( 4.0/5.0*field[H_INDEX(i-1,j,k,nx,ny,nz)] - 4.0/5.0*field[H_INDEX(i+1,j,k,nx,ny,nz)] )
1677  )/dx;
1678  break;
1679  case 2:
1680  return (
1681  ( 1.0/280.0*field[H_INDEX(i,j-4,k,nx,ny,nz)] - 1.0/280.0*field[H_INDEX(i,j+4,k,nx,ny,nz)] )
1682  - ( 4.0/105.0*field[H_INDEX(i,j-3,k,nx,ny,nz)] - 4.0/105.0*field[H_INDEX(i,j+3,k,nx,ny,nz)] )
1683  + ( 1.0/5.0*field[H_INDEX(i,j-2,k,nx,ny,nz)] - 1.0/5.0*field[H_INDEX(i,j+2,k,nx,ny,nz)] )
1684  - ( 4.0/5.0*field[H_INDEX(i,j-1,k,nx,ny,nz)] - 4.0/5.0*field[H_INDEX(i,j+1,k,nx,ny,nz)] )
1685  )/dy;
1686  break;
1687  case 3:
1688  return (
1689  ( 1.0/280.0*field[H_INDEX(i,j,k-4,nx,ny,nz)] - 1.0/280.0*field[H_INDEX(i,j,k+4,nx,ny,nz)] )
1690  - ( 4.0/105.0*field[H_INDEX(i,j,k-3,nx,ny,nz)] - 4.0/105.0*field[H_INDEX(i,j,k+3,nx,ny,nz)] )
1691  + ( 1.0/5.0*field[H_INDEX(i,j,k-2,nx,ny,nz)] - 1.0/5.0*field[H_INDEX(i,j,k+2,nx,ny,nz)] )
1692  - ( 4.0/5.0*field[H_INDEX(i,j,k-1,nx,ny,nz)] - 4.0/5.0*field[H_INDEX(i,j,k+1,nx,ny,nz)] )
1693  )/dz;
1694  break;
1695  }
1696 
1697  /* XXX */
1698  return 0;
1699 }
1700 
1701 inline real_t mixed_derivative_stencil_Odx2(idx_t i, idx_t j, idx_t k,idx_t nx, idx_t ny, idx_t nz, int d1, int d2, arr_t & field)
1702 {
1703  real_t dx = H_LEN_FRAC / nx, dy = H_LEN_FRAC / ny, dz = H_LEN_FRAC / nz;
1704  if( (d1 == 1 && d2 == 2) || (d1 == 2 && d2 == 1) ) {
1705  return (
1706  - field[H_INDEX(i+1,j-1,k,nx,ny,nz)] + field[H_INDEX(i+1,j+1,k,nx,ny,nz)]
1707  + field[H_INDEX(i-1,j-1,k,nx,ny,nz)] - field[H_INDEX(i-1,j+1,k,nx,ny,nz)]
1708  )/4.0/dx/dy;
1709  }
1710 
1711  if( (d1 == 1 && d2 == 3) || (d1 == 3 && d2 == 1) ) {
1712  return (
1713  - field[H_INDEX(i+1,j,k-1,nx,ny,nz)] + field[H_INDEX(i+1,j,k+1,nx,ny,nz)]
1714  + field[H_INDEX(i-1,j,k-1,nx,ny,nz)] - field[H_INDEX(i-1,j,k+1,nx,ny,nz)]
1715  )/4.0/dx/dz;
1716  }
1717 
1718  if( (d1 == 3 && d2 == 2) || (d1 == 2 && d2 == 3) ) {
1719  return (
1720  - field[H_INDEX(i,j+1,k-1,nx,ny,nz)] + field[H_INDEX(i,j+1,k+1,nx,ny,nz)]
1721  + field[H_INDEX(i,j-1,k-1,nx,ny,nz)] - field[H_INDEX(i,j-1,k+1,nx,ny,nz)]
1722  )/4.0/dy/dz;
1723  }
1724 
1725  /* XXX */
1726  return 0;
1727 }
1728 
1729 inline real_t mixed_derivative_stencil_Odx4(idx_t i, idx_t j, idx_t k,idx_t nx, idx_t ny, idx_t nz, int d1, int d2, arr_t & field)
1730 {
1731  real_t dx = H_LEN_FRAC / nx, dy = H_LEN_FRAC / ny, dz = H_LEN_FRAC / nz;
1732  if( (d1 == 1 && d2 == 2) || (d1 == 2 && d2 == 1) ) {
1733  return (
1734  (
1735  - field[H_INDEX(i+1,j-1,k,nx,ny,nz)] + field[H_INDEX(i+1,j+1,k,nx,ny,nz)]
1736  + field[H_INDEX(i-1,j-1,k,nx,ny,nz)] - field[H_INDEX(i-1,j+1,k,nx,ny,nz)]
1737  ) - 1.0/16.0*(
1738  - field[H_INDEX(i+2,j-2,k,nx,ny,nz)] + field[H_INDEX(i+2,j+2,k,nx,ny,nz)]
1739  + field[H_INDEX(i-2,j-2,k,nx,ny,nz)] - field[H_INDEX(i-2,j+2,k,nx,ny,nz)]
1740  )
1741  )/3.0/dx/dy;
1742  }
1743 
1744  if( (d1 == 1 && d2 == 3) || (d1 == 3 && d2 == 1) ) {
1745  return (
1746  (
1747  - field[H_INDEX(i+1,j,k-1,nx,ny,nz)] + field[H_INDEX(i+1,j,k+1,nx,ny,nz)]
1748  + field[H_INDEX(i-1,j,k-1,nx,ny,nz)] - field[H_INDEX(i-1,j,k+1,nx,ny,nz)]
1749  ) - 1.0/16.0*(
1750  - field[H_INDEX(i+2,j,k-2,nx,ny,nz)] + field[H_INDEX(i+2,j,k+2,nx,ny,nz)]
1751  + field[H_INDEX(i-2,j,k-2,nx,ny,nz)] - field[H_INDEX(i-2,j,k+2,nx,ny,nz)]
1752  )
1753  )/3.0/dx/dz;
1754  }
1755 
1756  if( (d1 == 3 && d2 == 2) || (d1 == 2 && d2 == 3) ) {
1757  return (
1758  (
1759  - field[H_INDEX(i,j+1,k-1,nx,ny,nz)] + field[H_INDEX(i,j+1,k+1,nx,ny,nz)]
1760  + field[H_INDEX(i,j-1,k-1,nx,ny,nz)] - field[H_INDEX(i,j-1,k+1,nx,ny,nz)]
1761  ) - 1.0/16.0*(
1762  - field[H_INDEX(i,j+2,k-2,nx,ny,nz)] + field[H_INDEX(i,j+2,k+2,nx,ny,nz)]
1763  + field[H_INDEX(i,j-2,k-2,nx,ny,nz)] - field[H_INDEX(i,j-2,k+2,nx,ny,nz)]
1764  )
1765  )/3.0/dy/dz;
1766  }
1767 
1768  /* XXX */
1769  return 0;
1770 }
1771 
1772 inline real_t mixed_derivative_stencil_Odx6(idx_t i, idx_t j, idx_t k, idx_t nx, idx_t ny, idx_t nz, int d1, int d2, arr_t & field)
1773 {
1774  real_t dx = H_LEN_FRAC / nx, dy = H_LEN_FRAC / ny, dz = H_LEN_FRAC / nz;
1775  if( (d1 == 1 && d2 == 2) || (d1 == 2 && d2 == 1) ) {
1776  return (
1777  135.0*(
1778  - field[H_INDEX(i+1,j-1,k,nx,ny,nz)] + field[H_INDEX(i+1,j+1,k,nx,ny,nz)]
1779  + field[H_INDEX(i-1,j-1,k,nx,ny,nz)] - field[H_INDEX(i-1,j+1,k,nx,ny,nz)]
1780  ) - 27.0/2.0*(
1781  - field[H_INDEX(i+2,j-2,k,nx,ny,nz)] + field[H_INDEX(i+2,j+2,k,nx,ny,nz)]
1782  + field[H_INDEX(i-2,j-2,k,nx,ny,nz)] - field[H_INDEX(i-2,j+2,k,nx,ny,nz)]
1783  ) + (
1784  - field[H_INDEX(i+3,j-3,k,nx,ny,nz)] + field[H_INDEX(i+3,j+3,k,nx,ny,nz)]
1785  + field[H_INDEX(i-3,j-3,k,nx,ny,nz)] - field[H_INDEX(i-3,j+3,k,nx,ny,nz)]
1786  )
1787  )/360.0/dx/dy;
1788  }
1789 
1790  if( (d1 == 1 && d2 == 3) || (d1 == 3 && d2 == 1) ) {
1791  return (
1792  135.0*(
1793  - field[H_INDEX(i+1,j,k-1,nx,ny,nz)] + field[H_INDEX(i+1,j,k+1,nx,ny,nz)]
1794  + field[H_INDEX(i-1,j,k-1,nx,ny,nz)] - field[H_INDEX(i-1,j,k+1,nx,ny,nz)]
1795  ) - 27.0/2.0*(
1796  - field[H_INDEX(i+2,j,k-2,nx,ny,nz)] + field[H_INDEX(i+2,j,k+2,nx,ny,nz)]
1797  + field[H_INDEX(i-2,j,k-2,nx,ny,nz)] - field[H_INDEX(i-2,j,k+2,nx,ny,nz)]
1798  ) + (
1799  - field[H_INDEX(i+3,j,k-3,nx,ny,nz)] + field[H_INDEX(i+3,j,k+3,nx,ny,nz)]
1800  + field[H_INDEX(i-3,j,k-3,nx,ny,nz)] - field[H_INDEX(i-3,j,k+3,nx,ny,nz)]
1801  )
1802  )/360.0/dx/dz;
1803  }
1804 
1805  if( (d1 == 3 && d2 == 2) || (d1 == 2 && d2 == 3) ) {
1806  return (
1807  135.0*(
1808  - field[H_INDEX(i,j+1,k-1,nx,ny,nz)] + field[H_INDEX(i,j+1,k+1,nx,ny,nz)]
1809  + field[H_INDEX(i,j-1,k-1,nx,ny,nz)] - field[H_INDEX(i,j-1,k+1,nx,ny,nz)]
1810  ) - 27.0/2.0*(
1811  - field[H_INDEX(i,j+2,k-2,nx,ny,nz)] + field[H_INDEX(i,j+2,k+2,nx,ny,nz)]
1812  + field[H_INDEX(i,j-2,k-2,nx,ny,nz)] - field[H_INDEX(i,j-2,k+2,nx,ny,nz)]
1813  ) + (
1814  - field[H_INDEX(i,j+3,k-3,nx,ny,nz)] + field[H_INDEX(i,j+3,k+3,nx,ny,nz)]
1815  + field[H_INDEX(i,j-3,k-3,nx,ny,nz)] - field[H_INDEX(i,j-3,k+3,nx,ny,nz)]
1816  )
1817  )/360.0/dy/dz;
1818  }
1819 
1820  /* XXX */
1821  return 0;
1822 }
1823 
1824 inline real_t mixed_derivative_stencil_Odx8(idx_t i, idx_t j, idx_t k, idx_t nx, idx_t ny, idx_t nz, int d1, int d2, arr_t & field)
1825 {
1826  real_t dx = H_LEN_FRAC / nx, dy = H_LEN_FRAC / ny, dz = H_LEN_FRAC / nz;
1827  if( (d1 == 1 && d2 == 2) || (d1 == 2 && d2 == 1) ) {
1828  return (
1829  2.0/5.0*(
1830  - field[H_INDEX(i+1,j-1,k,nx,ny,nz)] + field[H_INDEX(i+1,j+1,k,nx,ny,nz)]
1831  + field[H_INDEX(i-1,j-1,k,nx,ny,nz)] - field[H_INDEX(i-1,j+1,k,nx,ny,nz)]
1832  ) - 1.0/20.0*(
1833  - field[H_INDEX(i+2,j-2,k,nx,ny,nz)] + field[H_INDEX(i+2,j+2,k,nx,ny,nz)]
1834  + field[H_INDEX(i-2,j-2,k,nx,ny,nz)] - field[H_INDEX(i-2,j+2,k,nx,ny,nz)]
1835  ) + 2.0/315.0*(
1836  - field[H_INDEX(i+3,j-3,k,nx,ny,nz)] + field[H_INDEX(i+3,j+3,k,nx,ny,nz)]
1837  + field[H_INDEX(i-3,j-3,k,nx,ny,nz)] - field[H_INDEX(i-3,j+3,k,nx,ny,nz)]
1838  ) - 1.0/2240.0*(
1839  - field[H_INDEX(i+4,j-4,k,nx,ny,nz)] + field[H_INDEX(i+4,j+4,k,nx,ny,nz)]
1840  + field[H_INDEX(i-4,j-4,k,nx,ny,nz)] - field[H_INDEX(i-4,j+4,k,nx,ny,nz)]
1841  )
1842  )/dx/dy;
1843  }
1844 
1845  if( (d1 == 1 && d2 == 3) || (d1 == 3 && d2 == 1) ) {
1846  return (
1847  2.0/5.0*(
1848  - field[H_INDEX(i+1,j,k-1,nx,ny,nz)] + field[H_INDEX(i+1,j,k+1,nx,ny,nz)]
1849  + field[H_INDEX(i-1,j,k-1,nx,ny,nz)] - field[H_INDEX(i-1,j,k+1,nx,ny,nz)]
1850  ) - 1.0/20.0*(
1851  - field[H_INDEX(i+2,j,k-2,nx,ny,nz)] + field[H_INDEX(i+2,j,k+2,nx,ny,nz)]
1852  + field[H_INDEX(i-2,j,k-2,nx,ny,nz)] - field[H_INDEX(i-2,j,k+2,nx,ny,nz)]
1853  ) + 2.0/315.0*(
1854  - field[H_INDEX(i+3,j,k-3,nx,ny,nz)] + field[H_INDEX(i+3,j,k+3,nx,ny,nz)]
1855  + field[H_INDEX(i-3,j,k-3,nx,ny,nz)] - field[H_INDEX(i-3,j,k+3,nx,ny,nz)]
1856  ) - 1.0/2240.0*(
1857  - field[H_INDEX(i+4,j,k-4,nx,ny,nz)] + field[H_INDEX(i+4,j,k+4,nx,ny,nz)]
1858  + field[H_INDEX(i-4,j,k-4,nx,ny,nz)] - field[H_INDEX(i-4,j,k+4,nx,ny,nz)]
1859  )
1860  )/dx/dz;
1861  }
1862 
1863  if( (d1 == 3 && d2 == 2) || (d1 == 2 && d2 == 3) ) {
1864  return (
1865  2.0/5.0*(
1866  - field[H_INDEX(i,j+1,k-1,nx,ny,nz)] + field[H_INDEX(i,j+1,k+1,nx,ny,nz)]
1867  + field[H_INDEX(i,j-1,k-1,nx,ny,nz)] - field[H_INDEX(i,j-1,k+1,nx,ny,nz)]
1868  ) - 1.0/20.0*(
1869  - field[H_INDEX(i,j+2,k-2,nx,ny,nz)] + field[H_INDEX(i,j+2,k+2,nx,ny,nz)]
1870  + field[H_INDEX(i,j-2,k-2,nx,ny,nz)] - field[H_INDEX(i,j-2,k+2,nx,ny,nz)]
1871  ) + 2.0/315.0*(
1872  - field[H_INDEX(i,j+3,k-3,nx,ny,nz)] + field[H_INDEX(i,j+3,k+3,nx,ny,nz)]
1873  + field[H_INDEX(i,j-3,k-3,nx,ny,nz)] - field[H_INDEX(i,j-3,k+3,nx,ny,nz)]
1874  ) - 1.0/2240.0*(
1875  - field[H_INDEX(i,j+4,k-4,nx,ny,nz)] + field[H_INDEX(i,j+4,k+4,nx,ny,nz)]
1876  + field[H_INDEX(i,j-4,k-4,nx,ny,nz)] - field[H_INDEX(i,j-4,k+4,nx,ny,nz)]
1877  )
1878  )/dy/dz;
1879  }
1880 
1881  /* XXX */
1882  return 0;
1883 }
1884 
1885 inline real_t double_derivative_stencil_Odx2(idx_t i, idx_t j, idx_t k, idx_t nx, idx_t ny, idx_t nz, int d, arr_t & field)
1886 {
1887  real_t dx = H_LEN_FRAC / nx, dy = H_LEN_FRAC / ny, dz = H_LEN_FRAC / nz;
1888  switch (d) {
1889  case 1:
1890  return (
1891  field[H_INDEX(i-1,j,k,nx,ny,nz)]
1892  - 2.0*field[H_INDEX(i-0,j,k,nx,ny,nz)]
1893  + field[H_INDEX(i+1,j,k,nx,ny,nz)]
1894  )/dx/dx;
1895  break;
1896  case 2:
1897  return (
1898  field[H_INDEX(i,j-1,k,nx,ny,nz)]
1899  - 2.0*field[H_INDEX(i,j-0,k,nx,ny,nz)]
1900  + field[H_INDEX(i,j+1,k,nx,ny,nz)]
1901  )/dy/dy;
1902  break;
1903  case 3:
1904  return (
1905  field[H_INDEX(i,j,k-1,nx,ny,nz)]
1906  - 2.0*field[H_INDEX(i,j,k-0,nx,ny,nz)]
1907  + field[H_INDEX(i,j,k+1,nx,ny,nz)]
1908  )/dz/dz;
1909  break;
1910  }
1911 
1912  /* XXX */
1913  return 0;
1914 }
1915 
1916 inline real_t double_derivative_stencil_Odx4(idx_t i, idx_t j, idx_t k, idx_t nx, idx_t ny, idx_t nz, int d, arr_t & field)
1917 {
1918  real_t dx = H_LEN_FRAC / nx, dy = H_LEN_FRAC / ny, dz = H_LEN_FRAC / nz;
1919  switch (d) {
1920  case 1:
1921  return (
1922  - 1.0/12.0*field[H_INDEX(i-2,j,k,nx,ny,nz)]
1923  + 4.0/3.0*field[H_INDEX(i-1,j,k,nx,ny,nz)]
1924  - 5.0/2.0*field[H_INDEX(i-0,j,k,nx,ny,nz)]
1925  + 4.0/3.0*field[H_INDEX(i+1,j,k,nx,ny,nz)]
1926  - 1.0/12.0*field[H_INDEX(i+2,j,k,nx,ny,nz)]
1927  )/dx/dx;
1928  break;
1929  case 2:
1930  return (
1931  - 1.0/12.0*field[H_INDEX(i,j-2,k,nx,ny,nz)]
1932  + 4.0/3.0*field[H_INDEX(i,j-1,k,nx,ny,nz)]
1933  - 5.0/2.0*field[H_INDEX(i,j-0,k,nx,ny,nz)]
1934  + 4.0/3.0*field[H_INDEX(i,j+1,k,nx,ny,nz)]
1935  - 1.0/12.0*field[H_INDEX(i,j+2,k,nx,ny,nz)]
1936  )/dy/dy;
1937  break;
1938  case 3:
1939  return (
1940  - 1.0/12.0*field[H_INDEX(i,j,k-2,nx,ny,nz)]
1941  + 4.0/3.0*field[H_INDEX(i,j,k-1,nx,ny,nz)]
1942  - 5.0/2.0*field[H_INDEX(i,j,k-0,nx,ny,nz)]
1943  + 4.0/3.0*field[H_INDEX(i,j,k+1,nx,ny,nz)]
1944  - 1.0/12.0*field[H_INDEX(i,j,k+2,nx,ny,nz)]
1945  )/dz/dz;
1946  break;
1947  }
1948 
1949  /* XXX */
1950  return 0;
1951 }
1952 
1953 inline real_t double_derivative_stencil_Odx6(idx_t i, idx_t j, idx_t k, idx_t nx, idx_t ny, idx_t nz, int d, arr_t & field)
1954 {
1955  real_t dx = H_LEN_FRAC / nx, dy = H_LEN_FRAC / ny, dz = H_LEN_FRAC / nz;
1956  switch (d) {
1957  case 1:
1958  return (
1959  1.0/90.0*field[H_INDEX(i-3,j,k,nx,ny,nz)]
1960  - 3.0/20.0*field[H_INDEX(i-2,j,k,nx,ny,nz)]
1961  + 3.0/2.0*field[H_INDEX(i-1,j,k,nx,ny,nz)]
1962  - 49.0/18.0*field[H_INDEX(i-0,j,k,nx,ny,nz)]
1963  + 3.0/2.0*field[H_INDEX(i+1,j,k,nx,ny,nz)]
1964  - 3.0/20.0*field[H_INDEX(i+2,j,k,nx,ny,nz)]
1965  + 1.0/90.0*field[H_INDEX(i+3,j,k,nx,ny,nz)]
1966  )/dx/dx;
1967  break;
1968  case 2:
1969  return (
1970  1.0/90.0*field[H_INDEX(i,j-3,k,nx,ny,nz)]
1971  - 3.0/20.0*field[H_INDEX(i,j-2,k,nx,ny,nz)]
1972  + 3.0/2.0*field[H_INDEX(i,j-1,k,nx,ny,nz)]
1973  - 49.0/18.0*field[H_INDEX(i,j-0,k,nx,ny,nz)]
1974  + 3.0/2.0*field[H_INDEX(i,j+1,k,nx,ny,nz)]
1975  - 3.0/20.0*field[H_INDEX(i,j+2,k,nx,ny,nz)]
1976  + 1.0/90.0*field[H_INDEX(i,j+3,k,nx,ny,nz)]
1977  )/dy/dy;
1978  break;
1979  case 3:
1980  return (
1981  1.0/90.0*field[H_INDEX(i,j,k-3,nx,ny,nz)]
1982  - 3.0/20.0*field[H_INDEX(i,j,k-2,nx,ny,nz)]
1983  + 3.0/2.0*field[H_INDEX(i,j,k-1,nx,ny,nz)]
1984  - 49.0/18.0*field[H_INDEX(i,j,k-0,nx,ny,nz)]
1985  + 3.0/2.0*field[H_INDEX(i,j,k+1,nx,ny,nz)]
1986  - 3.0/20.0*field[H_INDEX(i,j,k+2,nx,ny,nz)]
1987  + 1.0/90.0*field[H_INDEX(i,j,k+3,nx,ny,nz)]
1988  )/dz/dz;
1989  break;
1990  }
1991 
1992  /* XXX */
1993  return 0;
1994 }
1995 
1996 inline real_t double_derivative_stencil_Odx8(idx_t i, idx_t j, idx_t k, idx_t nx, idx_t ny, idx_t nz, int d, arr_t & field)
1997 {
1998  real_t dx = H_LEN_FRAC / nx, dy = H_LEN_FRAC / ny, dz = H_LEN_FRAC / nz;
1999  switch (d) {
2000  case 1:
2001  return (
2002  - 1.0/560.0*field[H_INDEX(i-4,j,k,nx,ny,nz)]
2003  + 8.0/315.0*field[H_INDEX(i-3,j,k,nx,ny,nz)]
2004  - 1.0/5.0*field[H_INDEX(i-2,j,k,nx,ny,nz)]
2005  + 8.0/5.0*field[H_INDEX(i-1,j,k,nx,ny,nz)]
2006  - 205.0/72.0*field[H_INDEX(i-0,j,k,nx,ny,nz)]
2007  + 8.0/5.0*field[H_INDEX(i+1,j,k,nx,ny,nz)]
2008  - 1.0/5.0*field[H_INDEX(i+2,j,k,nx,ny,nz)]
2009  + 8.0/315.0*field[H_INDEX(i+3,j,k,nx,ny,nz)]
2010  - 1.0/560.0*field[H_INDEX(i+4,j,k,nx,ny,nz)]
2011  )/dx/dx;
2012  break;
2013  case 2:
2014  return (
2015  - 1.0/560.0*field[H_INDEX(i,j-4,k,nx,ny,nz)]
2016  + 8.0/315.0*field[H_INDEX(i,j-3,k,nx,ny,nz)]
2017  - 1.0/5.0*field[H_INDEX(i,j-2,k,nx,ny,nz)]
2018  + 8.0/5.0*field[H_INDEX(i,j-1,k,nx,ny,nz)]
2019  - 205.0/72.0*field[H_INDEX(i,j-0,k,nx,ny,nz)]
2020  + 8.0/5.0*field[H_INDEX(i,j+1,k,nx,ny,nz)]
2021  - 1.0/5.0*field[H_INDEX(i,j+2,k,nx,ny,nz)]
2022  + 8.0/315.0*field[H_INDEX(i,j+3,k,nx,ny,nz)]
2023  - 1.0/560.0*field[H_INDEX(i,j+4,k,nx,ny,nz)]
2024  )/dy/dy;
2025  break;
2026  case 3:
2027  return (
2028  - 1.0/560.0*field[H_INDEX(i,j,k-4,nx,ny,nz)]
2029  + 8.0/315.0*field[H_INDEX(i,j,k-3,nx,ny,nz)]
2030  - 1.0/5.0*field[H_INDEX(i,j,k-2,nx,ny,nz)]
2031  + 8.0/5.0*field[H_INDEX(i,j,k-1,nx,ny,nz)]
2032  - 205.0/72.0*field[H_INDEX(i,j,k-0,nx,ny,nz)]
2033  + 8.0/5.0*field[H_INDEX(i,j,k+1,nx,ny,nz)]
2034  - 1.0/5.0*field[H_INDEX(i,j,k+2,nx,ny,nz)]
2035  + 8.0/315.0*field[H_INDEX(i,j,k+3,nx,ny,nz)]
2036  - 1.0/560.0*field[H_INDEX(i,j,k+4,nx,ny,nz)]
2037  )/dz/dz;
2038  break;
2039  }
2040 
2041  /* XXX */
2042  return 0;
2043 }
2044 
2059 inline real_t derivative(idx_t i, idx_t j, idx_t k, idx_t nx, idx_t ny, idx_t nz, int d,
2060  arr_t & field)
2061 {
2062  return STENCIL_ORDER_FUNCTION(derivative_Odx)(i, j, k, nx, ny, nz, d, field);
2063 }
2064 
2080 inline real_t mixed_derivative_stencil(idx_t i, idx_t j, idx_t k, idx_t nx, idx_t ny, idx_t nz, int d1, int d2, arr_t & field)
2081 {
2082  return STENCIL_ORDER_FUNCTION(mixed_derivative_stencil_Odx)(i, j, k, nx, ny, nz, d1, d2, field);
2083 }
2084 
2099 inline real_t double_derivative_stencil(idx_t i, idx_t j, idx_t k, idx_t nx, idx_t ny, idx_t nz, int d, arr_t & field)
2100 {
2101  return STENCIL_ORDER_FUNCTION(double_derivative_stencil_Odx)(i, j, k, nx, ny, nz, d, field);
2102 }
2103 
2119 inline real_t double_derivative(idx_t i, idx_t j, idx_t k, idx_t nx, idx_t ny, idx_t nz, int d1, int d2, arr_t & field)
2120 {
2121  if(d1 == d2) {
2122  return double_derivative_stencil(i, j, k, nx, ny, nz, d1, field);
2123  } else {
2124  return mixed_derivative_stencil(i, j, k, nx, ny, nz, d1, d2, field);
2125  }
2126 
2127  /* XXX */
2128  return 0;
2129 }
2130 
2146 inline real_t laplacian(idx_t i, idx_t j, idx_t k, idx_t nx, idx_t ny, idx_t nz, arr_t & field)
2147 {
2148  return (
2149  double_derivative(i, j, k, nx, ny, nz, 1, 1, field)
2150  + double_derivative(i, j, k, nx, ny, nz, 2, 2, field)
2151  + double_derivative(i, j, k, nx, ny, nz, 3, 3, field)
2152  );
2153 }
2154 
2155 inline real_t interp(real_t i, real_t j, real_t k, idx_t nx, idx_t ny,
2156  idx_t nz, arr_t & field)
2157 {
2158  real_t C000 = field[INDEX( (idx_t) i, (idx_t) j, (idx_t) k )];
2159  real_t C001 = field[INDEX( (idx_t) i, (idx_t) j, (idx_t) k + 1 )];
2160  real_t C010 = field[INDEX( (idx_t) i, (idx_t) j + 1, (idx_t) k )];
2161  real_t C011 = field[INDEX( (idx_t) i, (idx_t) j + 1, (idx_t) k + 1 )];
2162  real_t C100 = field[INDEX( (idx_t) i + 1, (idx_t) j, (idx_t) k )];
2163  real_t C101 = field[INDEX( (idx_t) i + 1, (idx_t) j, (idx_t) k + 1 )];
2164  real_t C110 = field[INDEX( (idx_t) i + 1, (idx_t) j + 1, (idx_t) k )];
2165  real_t C111 = field[INDEX( (idx_t) i + 1, (idx_t) j + 1, (idx_t) k + 1 )];
2166 
2167  real_t i_frac = real_t_mod( i, 1.0 );
2168  real_t j_frac = real_t_mod( j, 1.0 );
2169  real_t k_frac = real_t_mod( k, 1.0 );
2170 
2171  real_t C00 = C000*(1.0 - i_frac) + C100*i_frac;
2172  real_t C01 = C001*(1.0 - i_frac) + C101*i_frac;
2173  real_t C10 = C010*(1.0 - i_frac) + C110*i_frac;
2174  real_t C11 = C011*(1.0 - i_frac) + C111*i_frac;
2175 
2176  real_t C0 = C00*(1.0 - j_frac) + C10*j_frac;
2177  real_t C1 = C01*(1.0 - j_frac) + C11*j_frac;
2178 
2179  // return "C"
2180  return C0*(1.0 - k_frac) + C1*k_frac;
2181 
2182  return 0.0;
2183 }
2184 
2185 }
2186 #endif
cosmo
Definition: bardeen.cc:5
Generated by Doxygen 1.8.11 | View cosmograph on GitHub