1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
|
#include <math.h>
#include "cctk.h"
void apply_dissipation(const cGH *gh, const double *var, double *rhs, double dx[3],
int order, double *epsdis)
{
double inv_dx[3] = { 1.0 / dx[0], 1.0 / dx[1], 1.0 / dx[2] };
int ni = gh->cctk_lsh[0];
int nj = gh->cctk_lsh[1];
int nk = gh->cctk_lsh[2];
double *y;
int y_idx = -1;
y = CCTK_VarDataPtr(gh, 0, "grid::y");
for (int j = 0; j < nj; j++) {
if (fabs(y[CCTK_GFINDEX3D(gh, 0, j, 0)]) < 1e-8) {
y_idx = j;
break;
}
}
switch (order) {
case 1:
#pragma omp parallel for
for (int k = 1; k < nk - 1; k++)
for (int i = 1; i < nk - 1; i++) {
int idx = CCTK_GFINDEX3D(gh, i, y_idx, k);
rhs[idx] += epsdis[idx] *
((var[idx - 1] - 2.0 * var[idx] + var[idx + 1]) * inv_dx[0] +
(var[idx - 1 * ni * nj] - 2.0 * var[idx] + var[idx + 1 * ni * nj]) * inv_dx[2]);
}
break;
case 3:
#pragma omp parallel for
for (int k = 2; k < nk - 2; k++)
for (int i = 2; i < nk - 2; i++) {
int idx = CCTK_GFINDEX3D(gh, i, y_idx, k);
rhs[idx] -= epsdis[idx] / 16.0 *
((var[idx - 2] - 4.0 * var[idx - 1] + 6.0 * var[idx] - 4.0 * var[idx + 1] + var[idx + 2]) * inv_dx[0] +
(var[idx - 2 * ni * nj] - 4.0 * var[idx - 1 * ni * nj] + 6.0 * var[idx] - 4.0 * var[idx + 1 * ni * nj] + var[idx + 2 * ni * nj]) * inv_dx[2]);
}
break;
case 5:
#pragma omp parallel for
for (int k = 3; k < nk - 3; k++)
for (int i = 3; i < nk - 3; i++) {
int idx = CCTK_GFINDEX3D(gh, i, y_idx, k);
rhs[idx] += epsdis[idx] / 64 *
((var[idx - 3] - 6 * var[idx - 2] + 15 * var[idx - 1] - 20 * var[idx] + 15 * var[idx + 1] - 6 * var[idx + 2] + var[idx + 3]) * inv_dx[0] +
(var[idx - 3 * ni * nj] - 6 * var[idx - 2 * ni * nj] + 15 * var[idx - 1 * ni * nj] - 20 * var[idx] + 15 * var[idx + 1 * ni * nj] - 6 * var[idx + 2 * ni * nj] + var[idx + 3 * ni * nj]) * inv_dx[2]);
// + (var(i,j-3,k) - 6*var(i,j-2,k) + 15*var(i,j-1,k) - 20*var(i,j,k) + 15*var(i,j+1,k) - 6*var(i,j+2,k) + var(i,j+3,k) ) * idx(2)
}
break;
case 9:
#pragma omp parallel for
for (int k = 5; k < nk - 5; k++)
for (int i = 5; i < nk - 5; i++) {
int idx = CCTK_GFINDEX3D(gh, i, y_idx, k);
rhs[idx] += epsdis[idx] / 1024 *
((var[idx - 5] - 10 * var[idx - 4] + 45 * var[idx - 3] - 120 * var[idx - 2] + 210 * var[idx - 1] - 252 * var[idx] +
var[idx + 5] - 10 * var[idx + 4] + 45 * var[idx + 3] - 120 * var[idx + 2] + 210 * var[idx + 1]) * inv_dx[0] +
(var[idx - 5 * ni * nj] - 10 * var[idx - 4 * ni * nj] + 45 * var[idx - 3 * ni * nj] - 120 * var[idx - 2 * ni * nj] + 210 * var[idx - 1 * ni * nj] - 252 * var[idx] +
var[idx + 5 * ni * nj] - 10 * var[idx + 4 * ni * nj] + 45 * var[idx + 3 * ni * nj] - 120 * var[idx + 2 * ni * nj] + 210 * var[idx + 1 * ni * nj]) * inv_dx[2]);
//$ + (var(i,j-5,k) - 10*var(i,j-4,k) + 45*var(i,j-3,k) - 120*var(i,j-2,k) + 210*var(i,j-1,k) - 252*var(i,j,k) + 210*var(i,j+1,k) - 120*var(i,j+2,k) + 45*var(i,j+3,k) - 10*var(i,j+4,k) + var(i,j+5,k)) * idx(2)
}
break;
default:
CCTK_WARN(0, "internal error");
}
}
|
