From 69cfc43e5c87e89e1bb1c44ef4c26689a1ebe708 Mon Sep 17 00:00:00 2001 From: lanfer Date: Tue, 26 Sep 2000 10:59:21 +0000 Subject: removing Fortran files git-svn-id: http://svn.cactuscode.org/arrangements/CactusElliptic/EllSOR/trunk@62 fa3da13c-9f13-4301-a575-cf5b8c5e1907 --- src/sor_confmetric.F | 329 --------------------------------------------------- src/sor_flat.F | 268 ----------------------------------------- 2 files changed, 597 deletions(-) delete mode 100644 src/sor_confmetric.F delete mode 100644 src/sor_flat.F diff --git a/src/sor_confmetric.F b/src/sor_confmetric.F deleted file mode 100644 index bee8e8e..0000000 --- a/src/sor_confmetric.F +++ /dev/null @@ -1,329 +0,0 @@ - /*@@ - @file sor.F - @date Thu Mar 13 07:46:55 1997 - @author Joan Masso + Paul Walker - @desc - The SOR solver - @enddesc - @@*/ - -#include "cctk.h" -#include "cctk_Arguments.h" -#include "cctk_Parameters.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_confmetric_core3d(_CCTK_FARGUMENTS, - $ Mlinear_lsh,Mlinear, - $ Nsource_lsh,Nsource, - $ gxx,gxy,gxz,gyy,gyz,gzz, - & psi,var, var_idx, - $ abstol,reltol) - - implicit none - - _DECLARE_CCTK_FARGUMENTS - DECLARE_CCTK_PARAMETERS - INTEGER CCTK_Equals - - 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 gxx(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)), gxy(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) - CCTK_REAL gxz(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)), gyy(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) - CCTK_REAL gyz(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)), gzz(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) - CCTK_REAL psi(cctk_lsh(1),cctk_lsh(2),cctk_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 pm4, tmp - CCTK_REAL psixp, psiyp, psizp - -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 Upper metric and determinant. What a pig. - CCTK_REAL uxx(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)), uyy(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) - CCTK_REAL uzz(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)), uxy(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) - CCTK_REAL uxz(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)), uyz(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) - CCTK_REAL det(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) - -c Coeeficients for the solver: 19 point stencil... - CCTK_REAL ac,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 reduction_handle,ierr - -c stencil size - INTEGER sw(3) - - tol = AbsTol(1) - -c Get the reduction handel for the sum operation - call CCTK_ReductionArrayHandle(reduction_handle,"sum"); - if (reduction_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(2)=1 - - call Ell_GetRealKey(ierr, finf, "EllLinConfMetric::Bnd::Robin::inf") - call Ell_GetIntKey (ierr, npow, "EllLinConfMetric::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 - - none = .false. - const = .false. - cheb = .false. - - if (CCTK_EQUALS(sor_accel,"const")) then - const = .true. - else if (CCTK_EQUALS(sor_accel,"cheb")) then - cheb = .true. - else - none = .true. - endif - - - 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.0d0 - -c compute determinant - det= -(gxz**2*gyy) + 2*gxy*gxz*gyz - gxx*gyz**2 - & - gxy**2*gzz + gxx*gyy*gzz - -c invert metric - uxx=(-gyz**2 + gyy*gzz)/det/psi**4 - uxy=(gxz*gyz - gxy*gzz)/det/psi**4 - uyy=(-gxz**2 + gxx*gzz)/det/psi**4 - uxz=(-gxz*gyy + gxy*gyz)/det/psi**4 - uyz=(gxy*gxz - gxx*gyz)/det/psi**4 - uzz=(-gxy**2 + gxx*gyy)/det/psi**4 - - det = det * psi**12 - - if (Mlinear_storage.eq.1) then - Mlinear = Mlinear*sqrt(det) - endif - if (Nsource_storage.eq.1) then - Nsource = Nsource*sqrt(det) - endif - -c Re-scaled uppwemetric. PATCHY but effective.... watch out: we get it back after -c the solve. - uxx=uxx/(2.*dx*dx)*sqrt(det) - uyy=uyy/(2.*dy*dy)*sqrt(det) - uzz=uzz/(2.*dz*dz)*sqrt(det) - uxy=uxy/(4.*dx*dy)*sqrt(det) - uxz=uxz/(4.*dx*dz)*sqrt(det) - uyz=uyz/(4.*dy*dz)*sqrt(det) - - - 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 -c which empirically seems to be pretty good - rjacobian = 1. - -c Set up the omega factor - omega = 1. - if (cheb) then - do i=2,sor_iteration - omega = 1./(1. - .25*rjacobian**2*omega) - enddo - endif - if (const) then - omega = 1.8 - endif - -c Total norm of the residual zero - resnorm = 0. -c End of first-iteration stuff - -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 = -uxx(i+1,j,k) -2.*uxx(i,j,k) -uxx(i-1,j,k) - & -uyy(i,j+1,k) -2.*uyy(i,j,k) -uyy(i,j-1,k) - & -uzz(i,j,k+1) -2.*uzz(i,j,k) -uzz(i,j,k-1) - if (Mlinear_storage.eq.1) then - ac = ac + Mlinear(i,j,k) - endif - - ae = uxx(i+1,j,k)+uxx(i,j,k) - aw = uxx(i-1,j,k)+uxx(i,j,k) - an = uyy(i,j+1,k)+uyy(i,j,k) - as = uyy(i,j-1,k)+uyy(i,j,k) - at = uzz(i,j,k+1)+uzz(i,j,k) - ab = uzz(i,j,k-1)+uzz(i,j,k) - ane = uxy(i,j+1,k)+uxy(i+1,j,k) - anw = - uxy(i-1,j,k)-uxy(i,j+1,k) - ase = - uxy(i+1,j,k)-uxy(i,j-1,k) - asw = uxy(i-1,j,k)+uxy(i,j-1,k) - ate = uxz(i,j,k+1)+uxz(i+1,j,k) - atw = - uxz(i-1,j,k)-uxz(i,j,k+1) - abe = - uxz(i+1,j,k)-uxz(i,j,k-1) - abw = uxz(i-1,j,k)+uxz(i,j,k-1) - atn = uyz(i,j+1,k)+uyz(i,j,k+1) - ats = - uyz(i,j,k+1)-uyz(i,j-1,k) - abn = - uyz(i,j,k-1)-uyz(i,j+1,k) - asb = uyz(i,j,k-1)+uyz(i,j-1,k) - - 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) - & + ane*var(i+1,j+1,k) + anw*var(i-1,j+1,k) - & + ase*var(i+1,j-1,k) + asw*var(i-1,j-1,k) - & + ate*var(i+1,j,k+1) + atw*var(i-1,j,k+1) - & + abe*var(i+1,j,k-1) + abw*var(i-1,j,k-1) - & + atn*var(i,j+1,k+1) + ats*var(i,j-1,k+1) - & + abn*var(i,j+1,k-1) + asb*var(i,j-1,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 - -c Reduce the norm - call CCTK_ReduceLocalScalar(ierr, cctkGH, -1, reduction_handle, - $ resnorm, residual, CCTK_VARIABLE_REAL) - if (ierr.ne.0) then - call CCTK_WARN(1,"Reduction of norm failed!"); - endif - residual = resnorm - residual = residual / (cctk_gsh(1)*cctk_gsh(2)*cctk_gsh(3)) - - call CCTK_SyncGroupWithVarI(cctkGH, var_idx) - - if (residual .lt. tol) then - goto 123 - endif - -c Apply Robin boundary - 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 - - write (*,*) "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" - write (*,*) "!! WARNING: SOR SOLVER DID NOT CONVERGE !!" - write (*,*) "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" - - 123 continue - if (Mlinear_storage.eq.1) then - Mlinear = Mlinear/sqrt(det) - endif - if (Nsource_storage.eq.1) then - Nsource = Nsource/sqrt(det) - endif - if (verb) write (*,*) "Iteration/Norm",sor_iteration,residual - - return - end - - - 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 - - - -- cgit v1.2.3