diff options
Diffstat (limited to 'src/sor_flat.F')
| -rw-r--r-- | src/sor_flat.F | 268 |
1 files changed, 0 insertions, 268 deletions
diff --git a/src/sor_flat.F b/src/sor_flat.F deleted file mode 100644 index 904f66d..0000000 --- a/src/sor_flat.F +++ /dev/null @@ -1,268 +0,0 @@ - /*@@ - @file sor.F - @date Thu Mar 13 07:46:55 1997 - @author Joan Masso + Paul Walker - @desc - The SOR solver - @enddesc - @version $Header$ - @@*/ - -#include "cctk.h" -#include "cctk_Arguments.h" -#include "cctk_Parameters.h" - -#include "EllBase.h" - - /*@@ - @routine sor - @date Thu Apr 24 13:29:52 1997 - @author Joan Masso - @desc - This is a standalone sor solver which does all the MPI itself. - It is a pretty good example of doing pugh-based communications - in FORTRAN. - @enddesc -@@*/ - - subroutine sor_flat_core3d( - & ierr,_CCTK_FARGUMENTS, - $ Mlinear_lsh,Mlinear, - $ Nsource_lsh,Nsource, - $ var,var_idx, - $ abstol,reltol) - - implicit none - - _DECLARE_CCTK_FARGUMENTS - DECLARE_CCTK_PARAMETERS - DECLARE_CCTK_FUNCTIONS - - integer ierr - INTEGER Mlinear_lsh(3) - CCTK_REAL Mlinear(Mlinear_lsh(1),Mlinear_lsh(2),Mlinear_lsh(3)) - - INTEGER Nsource_lsh(3) - CCTK_REAL Nsource(Nsource_lsh(1),Nsource_lsh(2),Nsource_lsh(3)) - - CCTK_REAL var(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) - INTEGER var_idx - - CCTK_REAL abstol(3),reltol(3) - - CCTK_REAL tol - INTEGER toltype - - CCTK_REAL dx,dy,dz - -c Temporaries - INTEGER sor_iteration - CCTK_REAL tmp - -c Numbers... - CCTK_REAL two, four - -c Total residual - CCTK_REAL resnorm, residual - -c Stencil - CCTK_REAL a(-1:1,-1:1,-1:1) - -c Loopers - INTEGER i,j,k - -c Acceleration factor - CCTK_REAL omega, rjacobian - -c conformal test - logical conformal - logical octant - -c Loop bounds. Starts, Ends, and deltas (steps) - INTEGER is, js, ks, ie, je, ke, di, dj, dk, kstep - -c Coeeficients for the solver: 19 point stencil... - CCTK_REAL ac,ac_orig,aw,ae,an,as,at,ab - CCTK_REAL ane,anw,ase,asw,ate,atw - CCTK_REAL abe,abw,atn,ats,abn,asb - - CCTK_REAL finf - CCTK_INT npow - - logical cheb, const, none, verb - integer Mlinear_storage,Nsource_storage - INTEGER sum_handle - -c stencil size - INTEGER sw(3) - - ierr = ELL_SUCCESS - - tol = AbsTol(1) - -c Get the reduction handel for the sum operation - call CCTK_ReductionArrayHandle(sum_handle,"sum"); - if (sum_handle.lt.0) then - call CCTK_WARN(1,"Cannot get reduction handle.") - endif - -c Set boundary related variables - if (CCTK_EQUALS(sor_bound,"robin")) then - sw(1)=1 - sw(2)=1 - sw(3)=1 - - call Ell_GetRealKey(ierr, finf, "EllLinFlat::Bnd::Robin::inf") - call Ell_GetIntKey (ierr, npow, "EllLinFlat::Bnd::Robin::falloff") - end if - -c We have no storage for M/N if they are of size one in each direction - if ((Mlinear_lsh(1).eq.1).and.(Mlinear_lsh(2).eq.1).and.(Mlinear_lsh(3).eq.1)) then - Mlinear_storage=0 - else - Mlinear_storage=1 - endif - - if ((Nsource_lsh(1).eq.1).and.(Nsource_lsh(2).eq.1).and.(Nsource_lsh(3).eq.1)) then - Nsource_storage=0 - else - Nsource_storage=1 - endif - -c Set up shorthand for the grid spacings - dx=cctk_delta_space(1) - dy=cctk_delta_space(2) - dz=cctk_delta_space(3) - - - verb = CCTK_Equals(elliptic_verbose,"yes").eq.1 - octant = CCTK_Equals(domain,"octant").eq.1 - -c cheb = contains("sor_accel","cheb").ne.0 -c const = contains("sor_accel","const").ne.0 -c none = contains("sor_accel","none").ne.0 - - verb = .true. - cheb = .false. - none = .false. - const = .false. - - if (verb .and. cheb) - $ print *,"Chebyshev Acceleration with radius of 1" - if (verb .and. const) - $ print *,"SOR with omega = 1.8" - if (verb .and. none) - $ print *,"Un-accelearted relaxation (omega = 1)" - - two = 2.0D0 - four = 4.0D0 - resnorm = 0 - - ae = 1.0d0/dx**2. - aw = 1.0d0/dx**2. - an = 1.0d0/dy**2. - as = 1.0d0/dy**2. - at = 1.0d0/dz**2. - ab = 1.0d0/dz**2. - ac_orig = -2.0d0/dx**2. - 2.0d0/dy**2. - 2.0d0/dz**2. - - do sor_iteration=1,maxit - -c We do not know the spectral radius of the Jacobi iteration, -c so we will take it to be one which empirically seems to be pretty good - rjacobian = 1.0D0 - -c Set up the omega factor - omega = 1.0D0 - if (cheb) then - do i=2,sor_iteration - omega = 1.0D0/(1.0D0 - .25D0*rjacobian**2*omega) - enddo - endif - if (const) then - omega = 1.8 - endif - -c Total norm of the residual zero - resnorm = 0. - -c Start loop with Red Black - ks = mod(sor_iteration,2)+2 - - if (cctk_lsh(3) .eq. 3) then - ks = 2 - endif - kstep = 2 - - do k=ks,cctk_lsh(3)-1,kstep - do j=2,cctk_lsh(2)-1 - do i=2,cctk_lsh(1)-1 - - ac = ac_orig - - if (Mlinear_storage.eq.1) then - ac = ac + Mlinear(i,j,k) - endif - - residual = ac*var(i,j,k) - & + ae*var(i+1,j,k) + aw*var(i-1,j,k) - & + an*var(i,j+1,k) + as*var(i,j-1,k) - & + at*var(i,j,k+1) + ab*var(i,j,k-1) - - if (Nsource_storage.eq.1) then - residual = residual + Nsource(i,j,k) - endif - -c Accumulate to the total Norm of the residual - resnorm = resnorm + abs(residual) -c Update - var(i,j,k) = var(i,j,k) - omega*residual/ac - end do - end do - end do - -c Reduce the norm - call CCTK_ReduceLocalScalar(ierr, cctkGH, -1, sum_handle, - $ resnorm, residual, CCTK_VARIABLE_REAL) - if (ierr.ne.0) then - call CCTK_WARN(1,"Reduction of norm failed!"); - endif - - residual = residual / (cctk_gsh(1)*cctk_gsh(2)*cctk_gsh(3)) - -c write (*,*) "Iteration/Norm",sor_iteration,residual - -c Synchronize the variables - call CCTK_SyncGroupWithVarI(cctkGH, var_idx) - - if (residual .lt. tol) then - goto 123 - endif - -c Apply boundary conditions -c call Ell_GetStringKey(nchar, mybound,"EllLinFlat::Bnd") - - if (CCTK_EQUALS(sor_bound,"robin")) then - call RobinBCVarI(ierr, cctkGH, finf, npow, sw, var_idx); - if (ierr.ne.0) then - call CCTK_WARN(1,"Could not apply Robin BC !") - endif - endif - -c Apply octant Symmetries - call CartSymVI(ierr, cctkGH, var_idx) - - enddo - - call CCTK_INFO("SOR solver did not converge") - ierr = ELL_NOCONVERGENCE - - 123 continue - - if (verb) write (*,*) "Iteration/Norm",sor_iteration,residual - - return - end - - - |