Kreiss-Oliger artificial dissipation for MoL.

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

9 files changed, 326 insertions, 0 deletions

diff --git a/README b/README
new file mode 100644
index 0000000..502b638
--- /dev/null
+++ b/README
@@ -0,0 +1,19 @@
+CVS info   : $Header$
+
+Cactus Code Thorn Dissipation
+Thorn Author(s)     : Erik Schnetter <schnetter@aei.mpg.de>
+Thorn Maintainer(s) : Erik Schnetter <schnetter@aei.mpg.de>
+--------------------------------------------------------------------------
+
+Purpose of the thorn:
+
+Add fourth order Kreiss-Oliger dissipation to the right hand side of
+evolution equations.  This thorn works only for time evolutions that
+use MoL.  I thank Scott Hawley who wrote a very similar thorn
+HawleyThorns/Dissipation for evolutions that do not use MoL; this
+thorn here is modelled after his.
+
+For a description of the artificial dissipation, see H. Kreiss and
+J. Oliger, "Methods for the Approximate Solution of Time Dependent
+Problems", vol. 10 of Global Atmospheric Research Programme (GARP):
+GARP Publication Series (GARP Publication, 1973).
diff --git a/doc/documentation.tex b/doc/documentation.tex
new file mode 100644
index 0000000..ecc3634
--- /dev/null
+++ b/doc/documentation.tex
@@ -0,0 +1,128 @@
+% *======================================================================*
+%  Cactus Thorn template for ThornGuide documentation
+%  Author: Ian Kelley
+%  Date: Sun Jun 02, 2002
+%  $Header$                                                             
+%
+%  Thorn documentation in the latex file doc/documentation.tex 
+%  will be included in ThornGuides built with the Cactus make system.
+%  The scripts employed by the make system automatically include 
+%  pages about variables, parameters and scheduling parsed from the 
+%  relevant thorn CCL files.
+%  
+%  This template contains guidelines which help to assure that your     
+%  documentation will be correctly added to ThornGuides. More 
+%  information is available in the Cactus UsersGuide.
+%                                                    
+%  Guidelines:
+%   - Do not change anything before the line
+%       % START CACTUS THORNGUIDE",
+%     except for filling in the title, author, date, etc. fields.
+%        - Each of these fields should only be on ONE line.
+%        - Author names should be separated with a \\ or a comma.
+%   - You can define your own macros, but they must appear after
+%     the START CACTUS THORNGUIDE line, and must not redefine standard 
+%     latex commands.
+%   - To avoid name clashes with other thorns, 'labels', 'citations', 
+%     'references', and 'image' names should conform to the following 
+%     convention:          
+%       ARRANGEMENT_THORN_LABEL
+%     For example, an image wave.eps in the arrangement CactusWave and 
+%     thorn WaveToyC should be renamed to CactusWave_WaveToyC_wave.eps
+%   - Graphics should only be included using the graphicx package. 
+%     More specifically, with the "\includegraphics" command.  Do
+%     not specify any graphic file extensions in your .tex file. This 
+%     will allow us to create a PDF version of the ThornGuide
+%     via pdflatex.
+%   - References should be included with the latex "\bibitem" command. 
+%   - Use \begin{abstract}...\end{abstract} instead of \abstract{...}
+%   - Do not use \appendix, instead include any appendices you need as 
+%     standard sections. 
+%   - For the benefit of our Perl scripts, and for future extensions, 
+%     please use simple latex.     
+%
+% *======================================================================* 
+% 
+% Example of including a graphic image:
+%    \begin{figure}[ht]
+% 	\begin{center}
+%    	   \includegraphics[width=6cm]{MyArrangement_MyThorn_MyFigure}
+% 	\end{center}
+% 	\caption{Illustration of this and that}
+% 	\label{MyArrangement_MyThorn_MyLabel}
+%    \end{figure}
+%
+% Example of using a label:
+%   \label{MyArrangement_MyThorn_MyLabel}
+%
+% Example of a citation:
+%    \cite{MyArrangement_MyThorn_Author99}
+%
+% Example of including a reference
+%   \bibitem{MyArrangement_MyThorn_Author99}
+%   {J. Author, {\em The Title of the Book, Journal, or periodical}, 1 (1999), 
+%   1--16. {\tt http://www.nowhere.com/}}
+%
+% *======================================================================* 
+
+% If you are using CVS use this line to give version information
+% $Header$
+
+\documentclass{article}
+
+% Use the Cactus ThornGuide style file
+% (Automatically used from Cactus distribution, if you have a 
+%  thorn without the Cactus Flesh download this from the Cactus
+%  homepage at www.cactuscode.org)
+\usepackage{../../../../doc/latex/cactus}
+
+\begin{document}
+
+% The author of the documentation
+\author{Erik Schnetter \textless schnetter@aei.mpg.de\textgreater}
+
+% The title of the document (not necessarily the name of the Thorn)
+\title{Dissipation}
+
+% the date your document was last changed, if your document is in CVS, 
+% please use:
+\date{$ $Date$ $}
+
+\maketitle
+
+% Do not delete next line
+% START CACTUS THORNGUIDE
+
+% Add all definitions used in this documentation here 
+%   \def\mydef etc
+
+% Add an abstract for this thorn's documentation
+\begin{abstract}
+Add fourth order Kreiss-Oliger dissipation to the right hand side of
+evolution equations.  This thorn is intended for time evolutions that
+use MoL.
+\end{abstract}
+
+% The following sections are suggestive only.
+% Remove them or add your own.
+
+\section{Physical System}
+For a description of the artificial dissipation, see
+\cite{kreiss-oliger}.
+
+\subsection{Acknowledgements}
+I thank Scott Hawley who wrote a very similar thorn
+\texttt{HawleyThorns/Dissipation} for evolutions that do not use MoL;
+this thorn here is modelled after his.
+
+\begin{thebibliography}{9}
+\bibitem{kreiss-oliger}
+H. Kreiss and J. Oliger, \emph{Methods for the Approximate Solution of
+Time Dependent Problems}, vol.\ 10 of Global Atmospheric Research
+Programme (GARP): GARP Publication Series (GARP Publication, 1973)
+\end{thebibliography}
+
+% Do not delete next line
+% END CACTUS THORNGUIDE
+
+\end{document}
diff --git a/interface.ccl b/interface.ccl
new file mode 100644
index 0000000..80a3e75
--- /dev/null
+++ b/interface.ccl
@@ -0,0 +1,7 @@
+# Interface definition for thorn Dissipation
+# $Header$
+
+IMPLEMENTS: Dissipation
+
+CCTK_INT FUNCTION MoLQueryEvolvedRHS (CCTK_INT IN EvolvedIndex)
+USES FUNCTION MoLQueryEvolvedRHS
diff --git a/param.ccl b/param.ccl
new file mode 100644
index 0000000..b50b7da
--- /dev/null
+++ b/param.ccl
@@ -0,0 +1,16 @@
+# Parameter definitions for thorn Dissipation
+# $Header$
+
+BOOLEAN verbose "produce log output" STEERABLE=always
+{
+} "no"
+
+REAL epsdis "Dissipation strength" STEERABLE=always
+{
+  *:* :: "0 for no dissipation.  Unstable for epsdis<0 and epsdis>1/"
+} 0.2
+
+STRING vars "List of evolved grid functions that should have dissipation added" STEERABLE=always
+{
+  .* :: "must be a valid list of grid functions"
+} ""
diff --git a/schedule.ccl b/schedule.ccl
new file mode 100644
index 0000000..759bf04
--- /dev/null
+++ b/schedule.ccl
@@ -0,0 +1,7 @@
+# Schedule definitions for thorn Dissipation
+# $Header$
+
+SCHEDULE dissipation_add IN MoL_PostRHS
+{
+  LANG: C
+} "Add Kreiss-Oliger dissipation to the right hand sides"
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");
+  }
+}