From cb37bf615da4d77d1114b2cbca9c9f4b74f28ff4 Mon Sep 17 00:00:00 2001
From: Anton Khirnov <anton@khirnov.net>
Date: Wed, 20 Feb 2019 12:13:34 +0100
Subject: Disable dissipation close to the boundary.

---
 schedule.ccl                 | 12 ++++++++++++
 src/maximal_slicing_axi_mg.c | 30 ++++++++++++++++++++++++++++++
 2 files changed, 42 insertions(+)

diff --git a/schedule.ccl b/schedule.ccl
index 0b992b9..7297aa7 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -38,6 +38,18 @@ if (CCTK_Equals(lapse_evolution_method, "maximal_axi_mg")) {
         LANG: C
     } ""
 
+    SCHEDULE maximal_slicing_axi_mg_modify_diss IN CCTK_POSTSTEP AFTER setup_epsdis {
+        LANG: C
+    } ""
+
+    SCHEDULE maximal_slicing_axi_mg_modify_diss IN CCTK_POSTREGRID AFTER setup_epsdis {
+        LANG: C
+    } ""
+
+    SCHEDULE maximal_slicing_axi_mg_modify_diss IN CCTK_POSTREGRIDINITIAL AFTER setup_epsdis {
+        LANG: C
+    } ""
+
     STORAGE: lapse_mg[2]
     STORAGE: lapse_mg_eval
     STORAGE: lapse_prev0
diff --git a/src/maximal_slicing_axi_mg.c b/src/maximal_slicing_axi_mg.c
index c37cbef..14bd2c8 100644
--- a/src/maximal_slicing_axi_mg.c
+++ b/src/maximal_slicing_axi_mg.c
@@ -871,3 +871,33 @@ void msa_mg_inithist(CCTK_ARGUMENTS)
         lapse_prev1_time[i] = DBL_MAX;
     }
 }
+
+void maximal_slicing_axi_mg_modify_diss(CCTK_ARGUMENTS)
+{
+    DECLARE_CCTK_ARGUMENTS;
+    DECLARE_CCTK_PARAMETERS;
+
+    CoordPatch *cp;
+
+    const int reflevel = ctz(cctkGH->cctk_levfac[0]);
+
+    double *epsdis;
+
+    epsdis = CCTK_VarDataPtr(cctkGH, 0, "Dissipation::epsdisA");
+    if (!epsdis)
+        abort();
+
+    cp = get_coord_patch(ms, reflevel);
+
+    for (int idx_z = 0; idx_z < cp->grid_size[2]; idx_z++)
+        for (int idx_x = cp->grid_size[0] - cp->offset_right[0] - (ms->fd_stencil + 1); idx_x < cp->grid_size[0]; idx_x++) {
+            const ptrdiff_t idx_dst = CPINDEX(cp, idx_x, cp->y_idx, idx_z);
+            epsdis[idx_dst] = 0.0;
+        }
+
+    for (int idx_x = 0; idx_x < cp->grid_size[0]; idx_x++)
+        for (int idx_z = cp->grid_size[2] - cp->offset_right[1] - (ms->fd_stencil + 1); idx_z < cp->grid_size[2]; idx_z++) {
+            const ptrdiff_t idx_dst = CPINDEX(cp, idx_x, cp->y_idx, idx_z);
+            epsdis[idx_dst] = 0.0;
+        }
+}
-- 
cgit v1.2.3