diff options
author | Anton Khirnov <anton@khirnov.net> | 2016-03-24 10:53:52 +0100 |
---|---|---|
committer | Anton Khirnov <anton@khirnov.net> | 2016-03-24 10:53:52 +0100 |
commit | 19df41ffde95df208c198ce3bf3041460280108e (patch) | |
tree | 94008c5ce2be74ce02b3cc466d028c970ebc9e40 /src/apply_dissipation.c | |
parent | 8c0efc3ebf2e8f5e8f55b52e7a977ffb39cd7520 (diff) |
Cartoonize, convert to C.
Diffstat (limited to 'src/apply_dissipation.c')
-rw-r--r-- | src/apply_dissipation.c | 40 |
1 files changed, 40 insertions, 0 deletions
diff --git a/src/apply_dissipation.c b/src/apply_dissipation.c new file mode 100644 index 0000000..6953757 --- /dev/null +++ b/src/apply_dissipation.c @@ -0,0 +1,40 @@ +#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 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; + default: + CCTK_WARN(0, "internal error"); + } +} |