CarpetPETSc: Beginning of an elliptic solver for mesh refinement

8 files changed, 384 insertions, 0 deletions

diff --git a/CarpetDev/CarpetPETSc/README b/CarpetDev/CarpetPETSc/README
new file mode 100644
index 000000000..c58d6e089
--- /dev/null
+++ b/CarpetDev/CarpetPETSc/README
@@ -0,0 +1,8 @@
+Cactus Code Thorn CarpetPETSc
+Thorn Author(s)     : Erik Schnetter <schnetter@cct.lsu.edu>
+Thorn Maintainer(s) : Erik Schnetter <schnetter@cct.lsu.edu>
+--------------------------------------------------------------------------
+
+Purpose of the thorn:
+
+Solve an elliptic equation with PETSc
diff --git a/CarpetDev/CarpetPETSc/configuration.ccl b/CarpetDev/CarpetPETSc/configuration.ccl
new file mode 100644
index 000000000..561fd590d
--- /dev/null
+++ b/CarpetDev/CarpetPETSc/configuration.ccl
@@ -0,0 +1,3 @@
+# Configuration definition for thorn CarpetPETSc
+
+REQUIRES PETSc
diff --git a/CarpetDev/CarpetPETSc/doc/documentation.tex b/CarpetDev/CarpetPETSc/doc/documentation.tex
new file mode 100644
index 000000000..d2681fe96
--- /dev/null
+++ b/CarpetDev/CarpetPETSc/doc/documentation.tex
@@ -0,0 +1,144 @@
+% *======================================================================*
+%  Cactus Thorn template for ThornGuide documentation
+%  Author: Ian Kelley
+%  Date: Sun Jun 02, 2002
+%  $Header: /cactusdevcvs/Cactus/doc/ThornGuide/template.tex,v 1.12 2004/01/07 20:12:39 rideout Exp $
+%
+%  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: /cactusdevcvs/Cactus/doc/ThornGuide/template.tex,v 1.12 2004/01/07 20:12:39 rideout Exp $
+
+\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@cct.lsu.edu\textgreater}
+
+% The title of the document (not necessarily the name of the Thorn)
+\title{CarpetPETSc}
+
+% the date your document was last changed, if your document is in CVS,
+% please use:
+%    \date{$ $Date: 2004/01/07 20:12:39 $ $}
+\date{February 22 2006}
+
+\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}
+
+\end{abstract}
+
+% The following sections are suggestive only.
+% Remove them or add your own.
+
+\section{Introduction}
+
+\section{Physical System}
+
+\section{Numerical Implementation}
+
+\section{Using This Thorn}
+
+\subsection{Obtaining This Thorn}
+
+\subsection{Basic Usage}
+
+\subsection{Special Behaviour}
+
+\subsection{Interaction With Other Thorns}
+
+\subsection{Examples}
+
+\subsection{Support and Feedback}
+
+\section{History}
+
+\subsection{Thorn Source Code}
+
+\subsection{Thorn Documentation}
+
+\subsection{Acknowledgements}
+
+
+\begin{thebibliography}{9}
+
+\end{thebibliography}
+
+% Do not delete next line
+% END CACTUS THORNGUIDE
+
+\end{document}
diff --git a/CarpetDev/CarpetPETSc/interface.ccl b/CarpetDev/CarpetPETSc/interface.ccl
new file mode 100644
index 000000000..c35818af2
--- /dev/null
+++ b/CarpetDev/CarpetPETSc/interface.ccl
@@ -0,0 +1,4 @@
+# Interface definition for thorn CarpetPETSc
+
+IMPLEMENTS: CarpetPETSc
+INHERITS: grid
diff --git a/CarpetDev/CarpetPETSc/param.ccl b/CarpetDev/CarpetPETSc/param.ccl
new file mode 100644
index 000000000..409df90b3
--- /dev/null
+++ b/CarpetDev/CarpetPETSc/param.ccl
@@ -0,0 +1 @@
+# Parameter definitions for thorn CarpetPETSc
diff --git a/CarpetDev/CarpetPETSc/schedule.ccl b/CarpetDev/CarpetPETSc/schedule.ccl
new file mode 100644
index 000000000..6774187a9
--- /dev/null
+++ b/CarpetDev/CarpetPETSc/schedule.ccl
@@ -0,0 +1 @@
+# Schedule definitions for thorn CarpetPETSc
diff --git a/CarpetDev/CarpetPETSc/src/make.code.defn b/CarpetDev/CarpetPETSc/src/make.code.defn
new file mode 100644
index 000000000..307fd8802
--- /dev/null
+++ b/CarpetDev/CarpetPETSc/src/make.code.defn
@@ -0,0 +1,7 @@
+# Main make.code.defn file for thorn CarpetPETSc
+
+# Source files in this directory
+SRCS = solve.cc
+
+# Subdirectories containing source files
+SUBDIRS = 
diff --git a/CarpetDev/CarpetPETSc/src/solve.cc b/CarpetDev/CarpetPETSc/src/solve.cc
new file mode 100644
index 000000000..64778e7a0
--- /dev/null
+++ b/CarpetDev/CarpetPETSc/src/solve.cc
@@ -0,0 +1,216 @@
+#include <cassert>
+#include <vector>
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+#include "util_ErrorCodes.h"
+#include "util_Table.h"
+
+#include "carpet.hh"
+#include "TATelliptic.h"
+
+#include "dist.hh"
+
+#include <petscvec.h>
+
+
+
+#ifdef CCTK_CXX_RESTRICT
+#  define restrict CCTK_CXX_RESTRICT
+#endif
+
+
+
+namespace CarpetPETSc {
+  
+  using namespace std;
+  using namespace Carpet;
+  
+  
+  
+  char const * const solver_name = "CarpetPETSc";
+  
+  
+  
+  struct options_t {
+    
+    // Width of outer boundary
+    CCTK_INT bndwidth;
+    
+    // Symmetry information
+    CCTK_INT symmetry_handle[2*dim];
+    CCTK_INT symmetry_zone_width[2*dim];
+  };
+  
+  
+  
+  void
+  copy_in (cGH const * restrict const cctkGH,
+           Vec dst,
+           vector<CCTK_INT> const & src,
+           options_t const & options);
+  
+  void
+  copy_out (cGH const * restrict const cctkGH,
+            vector<CCTK_INT> const & dst,
+            Vec const src,
+            options_t const & options);
+  
+  
+  
+  int
+  solve (cGH const * const cctkGH,
+         int const * const var,
+         int const * const val,
+         int const nvars,
+         int const options_table,
+         int (* const fun) (cGH const * cctkGH,
+                            int options_table,
+                            void * data),
+         int (* const bnd) (cGH const * cctkGH,
+                            int options_table,
+                            void * data),
+         void * const data)
+  {
+    DECLARE_CCTK_PARAMETERS;
+  }
+  
+  
+  
+  void
+  copy_in (cGH const * restrict const cctkGH,
+           Vec dst,
+           vector<CCTK_INT> const & src,
+           options_t const & options)
+  {
+    double * restrict dstptr;
+    VecGetArray (dst, & dstptr);
+    PetscInt dstsize;
+    VecGetSize (dst, & dstsize);
+    
+    int dstind = 0;
+    
+    BEGIN_MAP_LOOP (cctkGH, CCTK_GF) {
+      BEGIN_LOCAL_COMPONENT_LOOP (cctkGH, CCTK_GF) {
+        DECLARE_CCTK_ARGUMENTS;
+        
+        vector <CCTK_REAL const *> srcptrs (src.size());
+        for (int n = 0; n < src.size(); ++ n) {
+          srcptrs.at(n) =
+            static_cast<CCTK_REAL const *>
+            (CCTK_VarDataPtrI (cctkGH, 0, src.at(n)));
+          assert (srcptrs.at(n));
+        }
+        
+        int imin[3], imax[3];
+        interior_shape (cctkGH, options, imin, imax);
+        for (int k=imin[2]; k<imax[2]; ++k) {
+          for (int j=imin[1]; j<imax[1]; ++j) {
+            for (int i=imin[0]; i<imax[0]; ++i) {
+              int const srcind = CCTK_GFINDEX3D (cctkGH, i, j, k);
+              
+              for (int n = 0; n < src.size(); ++ n) {
+                assert (dstind < dstsize);
+                dstptr[dstind] = srcptr[srcind];
+                ++ dstind;
+              }
+              
+            }
+          }
+        }
+        
+      } END_LOCAL_COMPONENT_LOOP;
+    } END_MAP_LOOP;
+  }
+  
+  
+  
+  void
+  copy_out (cGH const * restrict const cctkGH,
+            vector<CCTK_INT> const & dst,
+            Vec const src,
+            options_t const & options)
+  {
+    double const * restrict srcptr;
+    VecGetArray (src, & srcptr);
+    PetscInt srcsize;
+    VecGetSize (dst, & srcsize);
+    
+    int srcind = 0;
+    
+    BEGIN_MAP_LOOP (cctkGH, CCTK_GF) {
+      BEGIN_LOCAL_COMPONENT_LOOP (cctkGH, CCTK_GF) {
+        DECLARE_CCTK_ARGUMENTS;
+        
+        vector <CCTK_REAL *> dstptrs (dst.size());
+        for (int n = 0; n < dst.size(); ++ n) {
+          dstptrs.at(n) =
+            static_cast<CCTK_REAL *> (CCTK_VarDataPtrI (cctkGH, 0, dst.at(n)));
+          assert (dstptrs.at(n));
+        }
+        
+        int imin[3], imax[3];
+        interior_shape (cctkGH, options, imin, imax);
+        for (int k=imin[2]; k<imax[2]; ++k) {
+          for (int j=imin[1]; j<imax[1]; ++j) {
+            for (int i=imin[0]; i<imax[0]; ++i) {
+              int const dstind = CCTK_GFINDEX3D (cctkGH, i, j, k);
+              
+              for (int n = 0; n < src.size(); ++ n) {
+                assert (srcind < srcsize);
+                dstptr[dstind] = srcptr[srcind];
+                ++ srcind;
+              }
+              
+            }
+          }
+        }
+        
+      } END_LOCAL_COMPONENT_LOOP;
+    } END_MAP_LOOP;
+  }
+  
+  
+  
+  // Calculate the shape of the interior
+  void
+  interior_shape (cGH const * restrict const cctkGH,
+                  options_t const & options,
+                  int * restrict const imin,
+                  int * restrict const imax)
+  {
+    DECLARE_CCTK_ARGUMENTS;
+    
+    for (int d=0; d<dim; ++d) {
+      
+      imin[d] = 0;
+      if (cctk_bbox[2*d]) {
+        if (options.symmetry_handle[2*d] >= 0) {
+          imin[d] += options.symmetry_zone_width[2*d];
+        } else {
+          imin[d] += options.bndwidth;
+        }
+      } else {
+        imin[d] += cctk_nghostzones[d];
+      }
+      
+      imax[d] = cctk_lsh[d];
+      if (cctk_bbox[2*d+1]) {
+        if (options.symmetry_handle[2*d+1] >= 0) {
+          imax[d] -= options.symmetry_zone_width[2*d+1];
+        } else {
+          imax[d] -= options.bndwidth;
+        }
+      } else {
+        imax[d] -= cctk_nghostzones[d];
+      }
+      
+      assert (imin[d] <= imax[d]);
+      
+    } // for d
+  }
+  
+  
+  
+} // namespace CarpetPETSc