This commit was generated by cvs2svn to compensate for changes in r2, which

included commits to RCS files with non-trunk default branches.

git-svn-id: http://svn.cactuscode.org/arrangements/CactusWave/IDScalarWaveElliptic/trunk@3 41e88fdd-2190-4c69-9c84-4659c8cf322e

7 files changed, 203 insertions, 0 deletions

diff --git a/README b/README
new file mode 100644
index 0000000..faeb579
--- /dev/null
+++ b/README
@@ -0,0 +1,7 @@
+Cactus Code Thorn IDScalarWaveElliptic
+Authors    : ...
+CVS info   : $Header$
+--------------------------------------------------------------------------
+
+Purpose of the thorn:
+
diff --git a/interface.ccl b/interface.ccl
new file mode 100644
index 0000000..1067c00
--- /dev/null
+++ b/interface.ccl
@@ -0,0 +1,18 @@
+# Interface definition for thorn IDScalarWaveElliptic
+# $Header$
+
+implements: idscalarwaveelliptic
+inherits:grid wavetoy
+friend: wavetoy
+ 
+private:
+
+cctk_real ellcoeffs type = GF
+{
+  Mcoeff,Ncoeff
+} "Coefficients of the elliptic equation"
+
+ 
+ 
+
+
diff --git a/par/charge.par b/par/charge.par
new file mode 100644
index 0000000..ec465d1
--- /dev/null
+++ b/par/charge.par
@@ -0,0 +1,45 @@
+# wavetoyf77_rad.par - wavetoy evolution with radiation boundaries
+# $Id$
+#
+# /*@@
+#   @file      wavetoyf77_rad.par
+#   @date      Sunday 24th July
+#   @author    Gabrielle Allen
+#   @desc
+#   Wavetoy parameter file demonstrating radiation boundaries in Octant
+#   @enddesc
+# @@*/          
+
+ActiveThorns = "idscalarwaveelliptic idscalarwave time wavetoyf77 pugh cartgrid3d ioutil ioascii"
+ 
+time::dtfac = 0.5
+
+idscalarwave::initial_data = "charge"
+idscalarwaveelliptic::radius = 0.5
+idscalarwaveelliptic::charge = 1
+
+wavetoyf77::bound = "radiation"
+
+grid::type = "BySpacing"
+grid::domain = "octant"
+grid::dxyz = 0.6
+
+driver::global_nx = 60
+driver::global_ny = 60
+driver::global_nz = 60
+
+cctk_itlast = 120
+
+IOASCII::out0D_every = 2
+IOASCII::out0D_vars = "wavetoy::phi"
+
+IOASCII::out1D_every = 2
+IOASCII::out1D_z = "no"
+IOASCII::out1D_y = "no"
+IOASCII::out1D_vars = "wavetoy::phi "
+
+IOASCII::outinfo_every = 10
+IOASCII::outinfo_vars = "wavetoy::phi"
+ 
+IO::outdir = "charge"
+
diff --git a/param.ccl b/param.ccl
new file mode 100644
index 0000000..0b40989
--- /dev/null
+++ b/param.ccl
@@ -0,0 +1,24 @@
+# Parameter definitions for thorn IDScalarWaveElliptic
+# $Header$
+
+shares: idscalarwave
+
+EXTENDS KEYWORD initial_data ""
+{
+  "charge"      :: "Uniformly charged sphere"
+}
+
+
+private:
+
+REAL radius "Radius of uniformly charged sphere"
+{
+ 0:* ::
+} 1.0
+
+REAL charge "Charge of uniformly charged sphere"
+{
+ 0:* ::
+} 1.0
+
+
diff --git a/schedule.ccl b/schedule.ccl
new file mode 100644
index 0000000..1205472
--- /dev/null
+++ b/schedule.ccl
@@ -0,0 +1,9 @@
+# Schedule definitions for thorn IDScalarWaveElliptic
+# $Header$
+
+schedule UniformCharge at CCTK_INITIAL
+{
+  STORAGE:       wavetoy::scalarevolve,ellcoeffs
+  COMMUNICATION: wavetoy::scalarevolve
+  LANG:          Fortran
+} "Initial data (Uniform Charge) for 3D wave equation"
diff --git a/src/SourceData.F b/src/SourceData.F
new file mode 100644
index 0000000..d16c0de
--- /dev/null
+++ b/src/SourceData.F
@@ -0,0 +1,91 @@
+
+#include "cctk.h"
+#include "cctk_arguments.h"
+#include "cctk_parameters.h"
+
+      subroutine UniformCharge(CCTK_FARGUMENTS)
+
+c     Find static field for a uniformly charge sphere
+c
+c     That is, solve Nabla^2 phi = - 4 pi rho
+c     where rho = Q/(4/3 * Pi * R^3)
+c     where Q is the total charge and R is the sphere radius
+
+      implicit none
+
+      DECLARE_CCTK_FARGUMENTS
+      DECLARE_CCTK_PARAMETERS
+
+      CCTK_REAL tolerance,pi
+      CCTK_REAL AbsTol(3), RelTol(3)
+
+      integer iphi,iMcoeff,iNcoeff
+      integer i,j,k
+      CCTK_REAL charge_factor
+
+      pi = 4.0*atan(1.0)
+
+      charge_factor = 4.0d0*pi*charge*3.0d0/(4.0d0*pi*radius**3)
+
+      tolerance = 1.0d-5
+
+      do k=1, cctk_lsh(3)
+         do j=1, cctk_lsh(2)
+            do i=1, cctk_lsh(1)
+                     
+               Mcoeff(i,j,k) = 0.0d0
+
+               if (r(i,j,k) <= radius) then
+                  Ncoeff(i,j,k) = charge_factor
+               else
+                  Ncoeff(i,j,k) = 0d0
+               end if
+
+            end do
+         end do
+      end do
+
+      phi(i,j,k) = 0.0d0
+
+      call CCTK_VarIndex (iMcoeff,  "idscalarwaveelliptic::Mcoeff")
+      call CCTK_VarIndex (iNcoeff,  "idscalarwaveelliptic::Ncoeff")
+      call CCTK_VarIndex (iphi,"wavetoy::phi")
+      
+      write (*,*) iMcoeff, iNcoeff, iphi
+
+      AbsTol(1)=1.0d-5
+      AbsTol(2)=1.0d-5
+      AbsTol(3)=1.0d-5
+
+      RelTol(1)=-1
+      RelTol(2)=-1
+      RelTol(3)=-1
+
+c     Call to elliptic solver will go here
+c      call Ell_LinConfMetricSolver(cctkGH, 
+c     $     metpsi_index, afield_index, Mlin_index, Nsrc_index, AbsTol, RelTol,
+c     $     "petsc")
+
+c     Cheat, and use exact solution
+      do k=1, cctk_lsh(3)
+         do j=1, cctk_lsh(2)
+            do i=1, cctk_lsh(1)
+                     
+               if (r(i,j,k) >= radius) then
+                  phi(i,j,k) = charge/r(i,j,k)
+               else
+                  phi(i,j,k) = charge/(2.0d0*radius**3)*
+     &                 (3.0d0*radius**2-r(i,j,k)**2)
+               end if
+
+c              Should put in an option for time symmetric first step
+               phi_old(i,j,k) = phi(i,j,k)
+
+            end do
+         end do
+      end do
+
+      
+
+      return
+      end
diff --git a/src/make.code.defn b/src/make.code.defn
new file mode 100644
index 0000000..0f9d963
--- /dev/null
+++ b/src/make.code.defn
@@ -0,0 +1,9 @@
+# Main make.code.defn file for thorn IDScalarWaveElliptic
+# $Header$
+
+# Source files in this directory
+SRCS = SourceData.F
+
+# Subdirectories containing source files
+SUBDIRS = 
+