diff options
Diffstat (limited to 'src/SourceData.F77')
| -rw-r--r-- | src/SourceData.F77 | 46 |
1 files changed, 35 insertions, 11 deletions
diff --git a/src/SourceData.F77 b/src/SourceData.F77 index bbd786c..5ab3157 100644 --- a/src/SourceData.F77 +++ b/src/SourceData.F77 @@ -26,18 +26,40 @@ c where Q is the total charge and R is the sphere radius CCTK_REAL pi CCTK_REAL AbsTol(3), RelTol(3) + CCTK_REAL charge_factor integer iphi,iMcoeff,iNcoeff integer i,j,k,ierr - CCTK_REAL charge_factor - pi = 4.0*atan(1.0) +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.0*atan(1.0) 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 @@ -50,13 +72,10 @@ c where Q is the total charge and R is the sphere radius end do end do end do + - phi(i,j,k) = 0.0d0 +c Set tolerance for stopping elliptic solve - call CCTK_VarIndex (iMcoeff, "idscalarwaveelliptic::Mcoeff") - call CCTK_VarIndex (iNcoeff, "idscalarwaveelliptic::Ncoeff") - call CCTK_VarIndex (iphi,"wavetoy::phi") - AbsTol(1)=1.0d-5 AbsTol(2)=1.0d-5 AbsTol(3)=1.0d-5 @@ -65,22 +84,27 @@ c where Q is the total charge and R is the sphere radius RelTol(2)=-1 RelTol(3)=-1 -c Call elliptic solver to fill out phi + +c Set any options needed for the elliptic solve c call Ell_SetStringKey(ierr, "yes", "EllLinFlat::Bnd::Robin") c write (*,*) "SourceData: SetString ",ierr - call Ell_SetRealKey(ierr, 1.0d0, "EllLinFlat::Bnd::Robin::inf") call Ell_SetIntKey (ierr, 1, "EllLinFlat::Bnd::Robin::falloff") +c Call elliptic solver to fill out phi + write (*,*) "charge: Going into elliptic solver" - call Ell_LinFlatSolver(ierr,cctkGH, + call Ell_LinFlatSolver( + & ierr, + & cctkGH, & iphi, & iMcoeff, iNcoeff, & AbsTol, RelTol, & "sor") - write (*,*) "charge: Exit elliptic solver" + write (*,*) "charge: Exit elliptic solver ierr =",ierr + c Set up last timestep ... assume time symmetry |