














































 







































































 


































































































































































 























































































































































































 



















































 










 








































































































































































































































 










                                
                                

                                
                                

                                
                                

                                
























































































 
















subroutine UniformCharge(cctk_dim,cctk_gsh,cctk_lsh,cctk_lbnd,cctk_ubnd,cctk_tile_min,cctk_tile_max,cctk_ash,cctk_alignment,cctk_al&
  &ignment_offset,cctk_from,cctk_to,cctk_bbox,cctk_delta_time,cctk_time,cctk_delta_space,cctk_origin_space,cctk_levfac,cctk_levoff,&
  &cctk_levoffdenom,cctk_timefac,cctk_convlevel,cctk_convfac,cctk_nghostzones,cctk_iteration,cctkGH,cctk_ash1,cctk_ash2,cctk_ash3, &
  &Mcoeff,Ncoeff,temp, coarse_dx,coarse_dy,coarse_dz,phi,phi_p,phi_p_p,r,x,y,z)







  implicit none


  INTEGER cctk_dim
  integer, parameter :: cctki_use_cctk_dim = kind(cctk_dim)
  INTEGER cctk_gsh (cctk_dim)
  integer, parameter :: cctki_use_cctk_gsh = kind(cctk_gsh)
  INTEGER cctk_lsh (cctk_dim)
  integer, parameter :: cctki_use_cctk_lsh = kind(cctk_lsh)
  INTEGER cctk_lbnd (cctk_dim)
  integer, parameter :: cctki_use_cctk_lbnd = kind(cctk_lbnd)
  INTEGER cctk_ubnd (cctk_dim)
  integer, parameter :: cctki_use_cctk_ubnd = kind(cctk_ubnd)
  INTEGER cctk_tile_min (cctk_dim)
  integer, parameter :: cctki_use_cctk_tile_min = kind(cctk_tile_min)
  INTEGER cctk_tile_max (cctk_dim)
  integer, parameter :: cctki_use_cctk_tile_max = kind(cctk_tile_max)
  INTEGER cctk_ash (cctk_dim)
  integer, parameter :: cctki_use_cctk_ash = kind(cctk_ash)
  INTEGER cctk_alignment
  integer, parameter :: cctki_use_cctk_alignment = kind(cctk_alignment)
  INTEGER cctk_alignment_offset
  integer, parameter :: cctki_use_cctk_alignment_offset = kind(cctk_alignment_offset)
  INTEGER cctk_from (cctk_dim)
  integer, parameter :: cctki_use_cctk_from = kind(cctk_from)
  INTEGER cctk_to (cctk_dim)
  integer, parameter :: cctki_use_cctk_to = kind(cctk_to)
  INTEGER cctk_bbox (2*cctk_dim)
  integer, parameter :: cctki_use_cctk_bbox = kind(cctk_bbox)
  REAL*8 cctk_delta_time
  integer, parameter :: cctki_use_cctk_delta_time = kind(cctk_delta_time)
  REAL*8 cctk_time
  integer, parameter :: cctki_use_cctk_time = kind(cctk_time)
  REAL*8 cctk_delta_space (cctk_dim)
  integer, parameter :: cctki_use_cctk_delta_space = kind(cctk_delta_space)
  REAL*8 cctk_origin_space (cctk_dim)
  integer, parameter :: cctki_use_cctk_origin_space = kind(cctk_origin_space)
  INTEGER cctk_levfac (cctk_dim)
  integer, parameter :: cctki_use_cctk_levfac = kind(cctk_levfac)
  INTEGER cctk_levoff (cctk_dim)
  integer, parameter :: cctki_use_cctk_levoff = kind(cctk_levoff)
  INTEGER cctk_levoffdenom (cctk_dim)
  integer, parameter :: cctki_use_cctk_levoffdenom = kind(cctk_levoffdenom)
  INTEGER cctk_timefac
  integer, parameter :: cctki_use_cctk_timefac = kind(cctk_timefac)
  INTEGER cctk_convlevel
  integer, parameter :: cctki_use_cctk_convlevel = kind(cctk_convlevel)
  INTEGER cctk_convfac
  integer, parameter :: cctki_use_cctk_convfac = kind(cctk_convfac)
  INTEGER cctk_nghostzones (cctk_dim)
  integer, parameter :: cctki_use_cctk_nghostzones = kind(cctk_nghostzones)
  INTEGER cctk_iteration
  integer, parameter :: cctki_use_cctk_iteration = kind(cctk_iteration)
  integer*8 cctkGH
  integer, parameter :: cctki_use_cctkGH = kind(cctkGH)
  INTEGER cctk_ash1
  integer, parameter :: cctki_use_cctk_ash1 = kind(cctk_ash1)
  INTEGER cctk_ash2
  integer, parameter :: cctki_use_cctk_ash2 = kind(cctk_ash2)
  INTEGER cctk_ash3
  integer, parameter :: cctki_use_cctk_ash3 = kind(cctk_ash3)
  REAL*8 Mcoeff (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Mcoeff = kind(Mcoeff)
  REAL*8 Ncoeff (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Ncoeff = kind(Ncoeff)
  REAL*8 temp (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_temp = kind(temp)
  REAL*8 coarse_dx
  integer, parameter :: cctki_use_coarse_dx = kind(coarse_dx)
  REAL*8 coarse_dy
  integer, parameter :: cctki_use_coarse_dy = kind(coarse_dy)
  REAL*8 coarse_dz
  integer, parameter :: cctki_use_coarse_dz = kind(coarse_dz)
  REAL*8 phi (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_phi = kind(phi)
  REAL*8 phi_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_phi_p = kind(phi_p)
  REAL*8 phi_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_phi_p_p = kind(phi_p_p)
  REAL*8 r (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_r = kind(r)
  REAL*8 x (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_x = kind(x)
  REAL*8 y (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_y = kind(y)
  REAL*8 z (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_z = kind(z)
  REAL*8 charge
  integer, parameter :: cctki_use_charge = kind(charge)
  REAL*8 radius
  integer, parameter :: cctki_use_radius = kind(radius)
  integer*8 solver
  integer, parameter :: cctki_use_solver = kind(solver)
  INTEGER*4 output_tmp
  integer, parameter :: cctki_use_output_tmp = kind(output_tmp)
  INTEGER*4 sor_maxit
  integer, parameter :: cctki_use_sor_maxit = kind(sor_maxit)
  COMMON /IDScalarWaveEllipticpriv/charge, radius, solver, output_tmp, sor_maxit
  integer*8 initial_data
  integer, parameter :: cctki_use_initial_data = kind(initial_data)
  COMMON /IDSCALARWAVErest/initial_data
  interface
  integer*8 function CCTK_PointerTo (var)
  implicit none
  type(*),dimension(..),TARGET :: var
  end function CCTK_PointerTo
  end interface
  interface
  integer function CCTK_Equals (arg1, arg2)
  implicit none
  integer*8 arg1
  character(*) arg2
  end function CCTK_Equals
  integer function CCTK_MyProc (cctkGH)
  implicit none
  integer*8 cctkGH
  end function CCTK_MyProc
  integer function CCTK_nProcs (cctkGH)
  implicit none
  integer*8 cctkGH
  end function CCTK_nProcs
  integer function CCTK_IsThornActive (name)
  implicit none
  character(*) name
  end function CCTK_IsThornActive
  integer*8 function CCTK_NullPointer ()
  implicit none
  end function CCTK_NullPointer
  end interface

  REAL*8 pi
  REAL*8 AbsTol(3), RelTol(3)
  REAL*8 charge_factor

  integer iphi,iMcoeff,iNcoeff
  integer ierr
  INTEGER*4 length

  character*30 fsolver
  character*200 infoline



  call CCTK_VarIndex (iMcoeff,  "idscalarwaveelliptic::Mcoeff")
  if (iMcoeff .lt. 0) then
    call CCTK_Warn(0,50,"SourceData.F90","IDScalarWaveElliptic","Grid variable index for Mcoeff not found")
  end if

  call CCTK_VarIndex (iNcoeff,  "idscalarwaveelliptic::Ncoeff")
  if (iNcoeff .lt. 0) then
    call CCTK_Warn(0,55,"SourceData.F90","IDScalarWaveElliptic","Grid variable index for Ncoeff not found")
  end if

  call CCTK_VarIndex (iphi,"wavetoy::phi")
  if (iphi .lt. 0) then
    call CCTK_Warn(0,60,"SourceData.F90","IDScalarWaveElliptic","Grid variable index for iphi not found")
  end if




  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



  AbsTol(1)=1.0d-5
  AbsTol(2)=1.0d-5
  AbsTol(3)=1.0d-5

  RelTol(1)=-1
  RelTol(2)=-1
  RelTol(3)=-1



  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")


  if ((CCTK_Equals(solver,"sor") .ne. 0)) 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("IDScalarWaveElliptic",infoline)



  call Ell_LinFlatSolver( &
         ierr, &
         cctkGH, & 
         iphi,  &
         iMcoeff, iNcoeff, &
         AbsTol, RelTol, &
         fsolver)

  if (ierr .eq. 0) then
    call CCTK_Info("IDScalarWaveElliptic","Leaving elliptic solver: solve successful")
  else if (ierr .eq. -2) then
    call CCTK_Info("IDScalarWaveElliptic","Leaving elliptic solver: solver failed to converge")
  else if (ierr .eq. -1) then
    call CCTK_Info("IDScalarWaveElliptic","Elliptic solver not found")
  else 
    write(infoline,'("Leaving elliptic solver: solve failed (Error ",I3,")")') ierr
    call CCTK_Info("IDScalarWaveElliptic",infoline)
  end if



  phi_p = phi



  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
