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");
+    }
+}