1 files changed, 144 insertions, 0 deletions
diff --git a/src/SourceData.F90 b/src/SourceData.F90
new file mode 100644
index 0000000..aee68bb
--- /dev/null
+++ b/src/SourceData.F90
@@ -0,0 +1,144 @@
+ /*@@
+   @file      SourceData.F90
+   @date      
+   @author    Gabrielle Allen
+   @desc 
+              Elliptic initial data for wave equation
+
+              Originally written in F77 fixed format. Converted to F90
+              free format by Peter Diener.
+   @enddesc 
+   @version $Header$
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
+
+#include "EllBase.h"
+
+subroutine UniformCharge(CCTK_ARGUMENTS)
+
+! Find static field for a uniformly charge sphere
+!
+! That is, solve Nabla^2 phi = - 4 pi rho
+! where rho = Q/(4/3 * Pi * R^3)
+! where Q is the total charge and R is the sphere radius
+
+  implicit none
+
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_FUNCTIONS
+
+  CCTK_REAL pi
+  CCTK_REAL AbsTol(3), RelTol(3)
+  CCTK_REAL charge_factor
+
+  integer iphi,iMcoeff,iNcoeff
+  integer i,j,k,ierr
+  CCTK_INT length
+
+  character*30 fsolver
+  character*200 infoline
+
+! Get variable indices for all grid functions
+
+  call CCTK_VarIndex (iMcoeff,  "idscalarwaveelliptic::Mcoeff")
+  if (iMcoeff .lt. 0) then
+    call CCTK_WARN(0,"Grid variable index for Mcoeff not found")
+  end if
+
+  call CCTK_VarIndex (iNcoeff,  "idscalarwaveelliptic::Ncoeff")
+  if (iNcoeff .lt. 0) then
+    call CCTK_WARN(0,"Grid variable index for Ncoeff not found")
+  end if
+
+  call CCTK_VarIndex (iphi,"wavetoy::phi")
+  if (iphi .lt. 0) then
+    call CCTK_WARN(0,"Grid variable index for iphi not found")
+  end if
+
+
+! Set up all coefficients and initial guess for solution
+
+  pi = 4.0d0*atan(1.0d0)
+  charge_factor = 4.0d0*pi*charge*3.0d0/(4.0d0*pi*radius**3)
+
+  phi = 0.0d0
+  Mcoeff = 0.0d0
+  
+  where ( r <= radius )
+    Ncoeff = charge_factor
+  elsewhere
+    Ncoeff = 0.0d0
+  end where
+
+! Set tolerance for stopping elliptic solve
+
+  AbsTol(1)=1.0d-5
+  AbsTol(2)=1.0d-5
+  AbsTol(3)=1.0d-5
+
+  RelTol(1)=-1
+  RelTol(2)=-1
+  RelTol(3)=-1
+
+! Set any options needed for the elliptic solve
+
+  call Ell_SetStrKey (ierr, "yes", "EllLinFlat::Bnd::Robin")
+  call Ell_SetRealKey(ierr, 0.0d0, "EllLinFlat::Bnd::Robin::inf")
+  call Ell_SetIntKey (ierr, 1,     "EllLinFlat::Bnd::Robin::falloff")
+
+! Set parameters for specific solvers
+  if (CCTK_EQUALS(solver,"sor")) then
+    call Ell_SetIntKey(ierr, sor_maxit,"Ell::SORmaxit")
+  end if
+     
+  call CCTK_FortranString(length,solver,fsolver)
+  write(infoline,'("Going into elliptic solver ",A)') fsolver
+  call CCTK_INFO(infoline)
+
+! Call elliptic solver to fill out phi
+
+  call Ell_LinFlatSolver( &
+         ierr, &
+         cctkGH, & 
+         iphi,  &
+         iMcoeff, iNcoeff, &
+         AbsTol, RelTol, &
+         fsolver)
+
+  if (ierr .eq. ELL_SUCCESS) then
+    call CCTK_INFO("Leaving elliptic solver: solve successful")
+  else if (ierr .eq. ELL_NOCONVERGENCE) then
+    call CCTK_INFO("Leaving elliptic solver: solver failed to converge")
+  else if (ierr .eq. ELL_NOSOLVER) then
+    call CCTK_INFO("Elliptic solver not found")
+  else 
+    write(infoline,'("Leaving elliptic solver: solve failed (Error ",I3,")")') ierr
+    call CCTK_INFO(infoline)
+  end if
+
+! Set up last timestep ... assume (first order) time symmetry
+
+  phi_p = phi
+
+! Output exact solution if required
+
+  if (output_tmp .eq. 1) then
+
+    where ( r .ge. radius )
+      temp = charge / r
+    elsewhere
+      temp = charge/(2.0d0*radius**3)*(3.0d0*radius**2-r**2)
+    end where
+
+    call CCTK_OutputVarAsByMethod(ierr,cctkGH, &
+            "idscalarwaveelliptic::temp", &
+            "IOASCII_1D","phi_exact")
+  end if      
+
+  return
+end