Kreiss-Oliger artificial dissipation for MoL.

git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/Dissipation/trunk@2 850bcc8b-0e4f-0410-8c26-8d28fbf1eda9

4 files changed, 149 insertions, 0 deletions

diff --git a/src/apply_dissipation.F77 b/src/apply_dissipation.F77
new file mode 100644
index 0000000..9363824
--- /dev/null
+++ b/src/apply_dissipation.F77
@@ -0,0 +1,27 @@
+c     $Header$
+
+#include "cctk.h"
+
+      subroutine apply_dissipation (var, rhs, ni, nj, nk, epsdis)
+      implicit none
+      
+      integer   ni, nj, nk
+      CCTK_REAL var(ni,nj,nk), rhs(ni,nj,nk)
+      CCTK_REAL epsdis
+      
+      integer   i, j, k
+      
+      do k = 3, nk-2
+         do j = 3, nj-2
+            do i = 3, ni-2
+               
+               rhs(i,j,k) = rhs(i,j,k) - epsdis / 16
+     $              * (  var(i-2,j,k) - 4*var(i-1,j,k) + 6*var(i,j,k) - 4*var(i+1,j,k) + var(i+2,j,k)
+     $                 + var(i,j-2,k) - 4*var(i,j-1,k) + 6*var(i,j,k) - 4*var(i,j+1,k) + var(i,j+2,k)
+     $                 + var(i,j,k-2) - 4*var(i,j,k-1) + 6*var(i,j,k) - 4*var(i,j,k+1) + var(i,j,k+2))
+               
+            end do
+         end do
+      end do
+      
+      end
diff --git a/src/dissipation.c b/src/dissipation.c
new file mode 100644
index 0000000..362bfbb
--- /dev/null
+++ b/src/dissipation.c
@@ -0,0 +1,93 @@
+/* $Header$ */
+
+#include <assert.h>
+#include <stdlib.h>
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+void CCTK_FCALL
+CCTK_FNAME(apply_dissipation) (CCTK_REAL const * const var,
+                               CCTK_REAL       * const rhs,
+                               int       const * const ni,
+                               int       const * const nj,
+                               int       const * const nk,
+                               CCTK_REAL const * const epsdis);
+
+static void
+apply (int const varindex, char const * const optstring, void * const arg);
+
+void
+dissipation_add (CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  CCTK_TraverseString (vars, apply, cctkGH, CCTK_GROUP_OR_VAR);
+}
+
+void
+apply (int const varindex, char const * const optstring, void * const arg)
+{
+  cGH const * const cctkGH = (cGH const *) arg;
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  int rhsindex;
+  int vargroup, rhsgroup;
+  cGroup vardata, rhsdata;
+  CCTK_REAL const * varptr;
+  CCTK_REAL       * rhsptr;
+  int n;
+  int ierr;
+  
+  assert (varindex >= 0);
+  
+  rhsindex = MoLQueryEvolvedRHS (varindex);
+  if (rhsindex < 0) {
+    char * const fullvarname = CCTK_FullName (varindex);
+    assert (fullvarname);
+    CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "There is no RHS variable registered with MoL for the evolved variable \"%s\"",
+                fullvarname);
+    free (fullvarname);
+  }
+  assert (rhsindex >= 0);
+  
+  if (verbose) {
+    char * const fullvarname = CCTK_FullName (varindex);
+    char * const fullrhsname = CCTK_FullName (rhsindex);
+    assert (fullvarname);
+    assert (fullrhsname);
+    CCTK_VInfo (CCTK_THORNSTRING,
+                "Applying dissipation to \"%s\" (RHS \"%s\")",
+                fullvarname, fullrhsname);
+    free (fullvarname);
+    free (fullrhsname);
+  }
+  
+  vargroup = CCTK_GroupIndexFromVarI (varindex);
+  assert (vargroup >= 0);
+  rhsgroup = CCTK_GroupIndexFromVarI (rhsindex);
+  assert (rhsgroup >= 0);
+  
+  ierr = CCTK_GroupData (vargroup, &vardata);
+  assert (!ierr);
+  ierr = CCTK_GroupData (rhsgroup, &rhsdata);
+  assert (!ierr);
+  
+  assert (vardata.grouptype == CCTK_GF);
+  assert (vardata.vartype == CCTK_VARIABLE_REAL);
+  assert (vardata.dim == cctk_dim);
+  assert (rhsdata.grouptype == CCTK_GF);
+  assert (rhsdata.vartype == CCTK_VARIABLE_REAL);
+  assert (rhsdata.dim == cctk_dim);
+  
+  varptr = CCTK_VarDataPtrI (cctkGH, 0, varindex);
+  assert (varptr);
+  rhsptr = CCTK_VarDataPtrI (cctkGH, 0, rhsindex);
+  assert (rhsptr);
+  
+  CCTK_FNAME(apply_dissipation)
+    (varptr, rhsptr, &cctk_lsh[0], &cctk_lsh[1], &cctk_lsh[2], &epsdis);
+}
diff --git a/src/make.code.defn b/src/make.code.defn
new file mode 100644
index 0000000..a36a5cf
--- /dev/null
+++ b/src/make.code.defn
@@ -0,0 +1,8 @@
+# Main make.code.defn file for thorn Dissipation
+# $Header$
+
+# Source files in this directory
+SRCS = apply_dissipation.F77 dissipation.c paramcheck.c
+
+# Subdirectories containing source files
+SUBDIRS = 
diff --git a/src/paramcheck.c b/src/paramcheck.c
new file mode 100644
index 0000000..3cce974
--- /dev/null
+++ b/src/paramcheck.c
@@ -0,0 +1,21 @@
+/* $Header$ */
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+void dissipation_paramcheck (CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  if (cctk_dim != 3) {
+    CCTK_PARAMWARN ("This thorn supports only dim=3");
+    CCTK_WARN (0, "This thorn supports only dim=3");
+  }
+  
+  if (cctk_nghostzones[0] < 2 || cctk_nghostzones[1] < 2
+      || cctk_nghostzones[2] < 2) {
+    CCTK_PARAMWARN ("This thorn requires at least 2 ghost zones");
+  }
+}