diff options
Diffstat (limited to 'src/SourceData.F77')
| -rw-r--r-- | src/SourceData.F77 | 175 |
1 files changed, 0 insertions, 175 deletions
diff --git a/src/SourceData.F77 b/src/SourceData.F77 deleted file mode 100644 index 1cacee5..0000000 --- a/src/SourceData.F77 +++ /dev/null @@ -1,175 +0,0 @@ - /*@@ - @file SourceData.F77 - @date - @author Gabrielle Allen - @desc - Elliptic initial data for wave equation - @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) - -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_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 - -c 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 - - -c 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) - - do k=1, cctk_lsh(3) - do j=1, cctk_lsh(2) - do i=1, cctk_lsh(1) - - phi(i,j,k) = 0.0d0 - - Mcoeff(i,j,k) = 0.0d0 - - if (r(i,j,k) .le. radius) then - Ncoeff(i,j,k) = charge_factor - else - Ncoeff(i,j,k) = 0.0d0 - end if - - end do - end do - end do - - -c 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 - - -c 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") - -c 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) - -c 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 - -c Set up last timestep ... assume (first order) time symmetry - - do k=1, cctk_lsh(3) - do j=1, cctk_lsh(2) - do i=1, cctk_lsh(1) - - phi_p(i,j,k) = phi(i,j,k) - - end do - end do - end do - - -c Output exact solution if required - - if (output_tmp .eq. 1) then - do k=1, cctk_lsh(3) - do j=1, cctk_lsh(2) - do i=1, cctk_lsh(1) - - if (r(i,j,k) .ge. radius) then - temp(i,j,k) = charge/r(i,j,k) - else - temp(i,j,k) = charge/(2.0d0*radius**3)* - & (3.0d0*radius**2-r(i,j,k)**2) - end if - - end do - end do - end do - - call CCTK_OutputVarAsByMethod(ierr,cctkGH, - & "idscalarwaveelliptic::temp", - & "IOASCII_1D","phi_exact") - end if - - return - end - - - - - - |