From 2671f9023f4c4ad1a08e2aa6dff1e96b4601816a Mon Sep 17 00:00:00 2001 From: allen Date: Fri, 24 Sep 1999 09:15:26 +0000 Subject: Developing git-svn-id: http://svn.cactuscode.org/arrangements/CactusWave/IDScalarWaveElliptic/trunk@6 41e88fdd-2190-4c69-9c84-4659c8cf322e --- src/SourceData.F | 49 +++++++++++++++++++++++++++++++++---------------- 1 file changed, 33 insertions(+), 16 deletions(-) (limited to 'src') diff --git a/src/SourceData.F b/src/SourceData.F index d16c0de..2913efb 100644 --- a/src/SourceData.F +++ b/src/SourceData.F @@ -16,19 +16,17 @@ c where Q is the total charge and R is the sphere radius DECLARE_CCTK_FARGUMENTS DECLARE_CCTK_PARAMETERS - CCTK_REAL tolerance,pi + CCTK_REAL pi CCTK_REAL AbsTol(3), RelTol(3) integer iphi,iMcoeff,iNcoeff - integer i,j,k + integer i,j,k,ierr 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) @@ -51,8 +49,6 @@ c where Q is the total charge and R is the sphere radius 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 @@ -62,30 +58,51 @@ c where Q is the total charge and R is the sphere radius 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") + call Ell_LinFlatSolver(cctkGH, + & iphi, + & iMcoeff, iNcoeff, + & AbsTol, RelTol, + & "sor") + +c Set up last timestep + do k=1, cctk_lsh(3) + do j=1, cctk_lsh(2) + do i=1, cctk_lsh(1) + + phi_old(i,j,k) = phi(i,j,k) + + end do + end do + end do + -c Cheat, and use exact solution +c Output exact solution if required 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) + temp(i,j,k) = charge/r(i,j,k) else - phi(i,j,k) = charge/(2.0d0*radius**3)* + temp(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 - + if (output_tmp==1) then + call CCTK_OutputVarAsByMethod(ierr,cctkGH, + & "idscalarwaveelliptic::temp", + & "IOASCII_1D","phi_exact") + end if return end + + + + + + -- cgit v1.2.3