aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorallen <allen@0a4070d5-58f5-498f-b6c0-2693e757fa0f>1999-11-01 11:27:16 +0000
committerallen <allen@0a4070d5-58f5-498f-b6c0-2693e757fa0f>1999-11-01 11:27:16 +0000
commitcab79658b8e19947c3b4058f75961f01a055ef2e (patch)
tree0532e9dac09c72f6e1aa2942186ac65dba2b5e66 /src
parentf6ebcf616a79e8a6ed1683187dc8d58e33980c21 (diff)
This commit was generated by cvs2svn to compensate for changes in r2, which
included commits to RCS files with non-trunk default branches. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/IDAxiBrillBH/trunk@3 0a4070d5-58f5-498f-b6c0-2693e757fa0f
Diffstat (limited to 'src')
-rw-r--r--src/AxiBrillBHIVP.F339
-rw-r--r--src/IDAxiBrillBH.F339
-rw-r--r--src/bhbrill.m126
-rw-r--r--src/bhbrill.x60
-rw-r--r--src/gij.x133
-rw-r--r--src/interp2.F213
-rw-r--r--src/make.code.defn9
-rw-r--r--src/mg59p.F896
-rw-r--r--src/psi_1st_deriv.x18
-rw-r--r--src/psi_2nd_deriv.x51
-rw-r--r--src/shmgp.F771455
11 files changed, 3639 insertions, 0 deletions
diff --git a/src/AxiBrillBHIVP.F b/src/AxiBrillBHIVP.F
new file mode 100644
index 0000000..35faded
--- /dev/null
+++ b/src/AxiBrillBHIVP.F
@@ -0,0 +1,339 @@
+c/*@@
+c @file AxiBrillBHIVP.F
+c @date
+c @author
+c @desc
+c
+c @enddesc
+c@@*/
+
+#include "cctk.h"
+#include "cctk_parameters.h"
+#include "cctk_arguments.h"
+
+c/*@@
+c @routine AxiBrillBHIVP
+c @date
+c @author
+c @desc
+c
+c @enddesc
+c @calls
+c @calledby
+c @history
+c
+c @endhistory
+c@@*/
+
+c Need include file from Einstein
+#include "CactusEinstein/Einstein/src/Einstein.h"
+
+ subroutine AxiBrillBHIVP(CCTK_FARGUMENTS)
+
+ implicit none
+
+ DECLARE_CCTK_FARGUMENTS
+ DECLARE_CCTK_PARAMETERS
+
+c Perhaps this and others should go into cctk.h
+ integer CCTK_Equals
+
+
+ real*8 axibheps, rmax, dq, deta
+ integer levels,id5,id9,idi,idg,ier
+ real*8, allocatable :: cc(:,:),ce(:,:),cw(:,:),cn(:,:),cs(:,:),
+ $ rhs(:,:),psi2d(:,:),detapsi2d(:,:),dqpsi2d(:,:),
+ $ detaetapsi2d(:,:),detaqpsi2d(:,:),dqqpsi2d(:,:)
+ real*8, allocatable :: etagrd(:),qgrd(:)
+ real*8, allocatable :: eta(:,:,:),abseta(:,:,:),sign_eta(:,:,:),
+ $ q(:,:,:),phi(:,:,:)
+ real*8, allocatable :: psi2dv(:,:,:),detapsi2dv(:,:,:),
+ $ dqpsi2dv(:,:,:),detaetapsi2dv(:,:,:),detaqpsi2dv(:,:,:),
+ $ dqqpsi2dv(:,:,:)
+ parameter(axibheps = 1.0d-12)
+ real*8 ep1,ep2
+ real*8 o1,o2,o3,o4,o5,o6,o7,o8,o9
+ real*8 o10,o11,o12,o13,o14,o15,o16,o17,o18,o19
+ real*8 o20,o21,o22,o23,o24,o25,o26,o27,o28,o29
+ real*8 o30,o31,o32,o33,o34,o35,o36,o37,o38,o39
+ real*8 o40,o41,o42,o43,o44,o45,o46,o47,o48,o49
+ real*8 o50,o51,o52,o53,o54,o55,o56,o57,o58,o59
+ real*8 o60,o61,o62,o63,o64,o65,o66,o67,o68,o69
+ real*8 o70,o71,o72,o73,o74,o75,o76,o77,o78,o79
+ real*8 o80,o81,o82,o83,o84,o85,o86,o87,o88,o89
+ real*8 o90,o91,o92,o93,o94,o95,o96,o97,o98,o99
+ real*8 pi
+ real*8 adm
+ integer :: nx,ny,nz
+ integer i,j,k,nquads
+ integer npoints,handle,ierror
+
+ pi = 4.0d0*atan(1.0d0)
+
+c Set up the grid spacings
+ nx = cctk_lsh(1)
+ ny = cctk_lsh(2)
+ nz = cctk_lsh(3)
+
+c Brill wave parameters
+
+ print *,"Brill wave + BH Axisymmetric solve"
+ write (*,123)amp,eta0,sigma,n
+ print *,'etamax=',etamax
+ 123 format(1x,'Pars : Amp',f8.5,' eta0',f8.5,' sigma',f8.5,' n ',i3)
+
+c Solve on this sized cartesian grid
+c 2D grid size NE x NQ
+c Add 2 zones for boundaries...
+ ne = ne+2
+ nq = nq+2
+ ! do I need to call free?
+ allocate(cc(ne,nq),ce(ne,nq),cw(ne,nq),cn(ne,nq),cs(ne,nq),
+ $ rhs(ne,nq),psi2d(ne,nq),detapsi2d(ne,nq),dqpsi2d(ne,nq),
+ $ detaetapsi2d(ne,nq),detaqpsi2d(ne,nq),dqqpsi2d(ne,nq),
+ $ etagrd(ne),qgrd(nq))
+ allocate(eta(nx,ny,nz),abseta(nx,ny,nz),sign_eta(nx,ny,nz),
+ $ q(nx,ny,nz),phi(nx,ny,nz),
+ $ psi2dv(nx,ny,nz),detapsi2dv(nx,ny,nz),dqpsi2dv(nx,ny,nz),
+ $ detaetapsi2dv(nx,ny,nz),detaqpsi2dv(nx,ny,nz),
+ $ dqqpsi2dv(nx,ny,nz))
+c Initialize some arrays
+ psi2d = 1.0d0
+ detapsi2d = 0.0d0
+
+ nquads = 2
+ dq = nquads*0.5*pi/(nq-2)
+ deta = etamax/(ne-3)
+
+ do j=1,nq
+ qgrd(j) = (j-1.5)*dq
+ do i=1,ne
+ etagrd(i) = (i-2)*deta
+#include "Development/IDAxiBrillBH/src/bhbrill.x"
+ enddo
+ enddo
+c Boundary conditions
+ do j=1,nq
+ ce(2,j)=ce(2,j)+cw(2,j)
+ cw(2,j)=0.0
+
+ cw(ne-1,j)=cw(ne-1,j)+ce(ne-1,j)
+ cc(ne-1,j)=cc(ne-1,j)-deta*ce(ne-1,j)
+ ce(ne-1,j)=0.0
+
+ enddo
+ do i=1,ne
+ cc(i,2)=cc(i,2)+cs(i,2)
+ cs(i,2)=0.0
+ cc(i,nq-1)=cc(i,nq-1)+cn(i,nq-1)
+ cn(i,nq-1)=0.0
+ enddo
+
+c Do the solve
+ print *, " Calling axisymmetric solver"
+ call mgparm (levels,5,id5,id9,idi,idg,ne,nq)
+ call mg5 (ne,2,ne-1,nq,2,nq-1,
+ $ cc,cn,cs,cw,ce,psi2d,rhs,
+ $ id5,id9,idi,idg,1,axibheps,rmax,ier)
+ print *, " Solve complete"
+c The solution is now available.
+c Debugging is needed, a stop statement should
+c be called if the IVP solve is not successful
+
+ if(ier .ne. 0) stop 'bad solution to brill wave problem'
+ print *,'rmax = ',rmax
+ print *,'axibheps = ',axibheps
+ print *,'psi2d = ',maxval(psi2d),' ',minval(psi2d)
+
+ ep2 = 0.0
+ do j=2,nq-1
+ do i=2,ne-1
+ ep1 = rhs(i,j)-psi2d(i,j)*cc(i,j)-psi2d(i,j+1)*cn(i,j)-psi2d(i,j-1)*cs(i,j)-
+ & psi2d(i+1,j)*ce(i,j)-psi2d(i-1,j)*cw(i,j)
+ ep2 = max(abs(ep1),ep2)
+ enddo
+ enddo
+ print *,'Resulting eps =',ep2
+
+ ! what a pain all this is....
+ do j=1,nq
+ psi2d(1,j)=psi2d(3,j)
+ psi2d(ne,j)=-deta*psi2d(ne-1,j)+psi2d(ne-2,j)
+ enddo
+ do i=1,ne
+ psi2d(i,1)=psi2d(i,2)
+ psi2d(i,nq)=psi2d(i,nq-1)
+ enddo
+c goto 111
+ do j=2,nq-1
+ do i=2,ne-1
+ dqpsi2d(i,j)=0.5*(psi2d(i,j+1)-psi2d(i,j-1))/dq
+ dqqpsi2d(i,j)=(psi2d(i,j+1)+psi2d(i,j-1)-2.*psi2d(i,j))/dq**2
+ detapsi2d(i,j)=sinh(0.5*etagrd(i))+0.5*(psi2d(i+1,j)-psi2d(i-1,j))/deta
+ detaetapsi2d(i,j)=0.5*cosh(0.5*etagrd(i))+
+ $ (psi2d(i+1,j)+psi2d(i-1,j)-2.*psi2d(i,j))/deta**2
+ enddo
+ enddo
+ do j=1,nq
+ detapsi2d(1,j)=-detapsi2d(3,j)
+ detapsi2d(ne,j)=detapsi2d(ne-2,j) ! simplified
+
+ detaetapsi2d(1,j)=detaetapsi2d(3,j)
+ detaetapsi2d(ne,j)=detaetapsi2d(ne-2,j) ! simplified...
+
+ dqqpsi2d(1,j)=dqqpsi2d(3,j)
+ dqqpsi2d(ne,j)=dqqpsi2d(ne-2,j) ! simplified
+
+ dqpsi2d(1,j)=dqpsi2d(3,j)
+ dqpsi2d(ne,j)=-dq*dqpsi2d(ne-1,j)+dqpsi2d(ne-2,j)
+ enddo
+ do i=1,ne
+ detapsi2d(i,1)=detapsi2d(i,2)
+ detapsi2d(i,nq)=detapsi2d(i,nq-1)
+
+ detaetapsi2d(i,1)=detaetapsi2d(i,2)
+ detaetapsi2d(i,nq)=detaetapsi2d(i,nq-1)
+
+ dqqpsi2d(i,1)=dqqpsi2d(i,2)
+ dqqpsi2d(i,nq)=dqqpsi2d(i,nq-1)
+
+ dqpsi2d(i,1)=-dqpsi2d(i,2)
+ dqpsi2d(i,nq)=-dqpsi2d(i,nq-1)
+ enddo
+ do j=2,nq-1
+ do i=2,ne-1
+ detaqpsi2d(i,j)=0.5*(detapsi2d(i,j+1)-detapsi2d(i,j-1))/dq
+ enddo
+ enddo
+ do j=1,nq
+ detaqpsi2d(1,j)=-detaqpsi2d(3,j)
+ detaqpsi2d(ne,j)=detaqpsi2d(ne-2,j) ! simplified
+ enddo
+ do i=1,ne
+ detaqpsi2d(i,1)=-detaqpsi2d(i,2)
+ detaqpsi2d(i,nq)=-detaqpsi2d(i,nq-1)
+ enddo
+ do j=1,nq
+ psi2d(:,j)=psi2d(:,j)+2.*cosh(0.5*etagrd)
+ enddo
+
+c Now evaluate each of the following at x(i,j,k), y(i,j,k) and
+c z(i,j,k) where i,j,k go from 1 to nx,ny,nz
+
+c Conformal factor
+
+ eta = 0.5d0 * dlog (x**2 + y**2 + z**2)
+ abseta = abs (eta)
+ q = datan2 (sqrt (x**2 + y**2), z)
+ phi = datan2 (y, x)
+
+ do k=1,nz
+ do j=1,ny
+ do i=1,nx
+c eta(i,j,k) = 0.5d0*dlog(x(i,j,k)**2+y(i,j,k)**2+z(i,j,k)**2)
+c abseta(i,j,k) = abs(eta(i,j,k))
+ if(eta(i,j,k) .lt. 0)then
+ sign_eta(i,j,k) = -1
+ else
+ sign_eta(i,j,k) = 1
+ endif
+c q(i,j,k) = datan2(sqrt(x(i,j,k)**2+y(i,j,k)**2),z(i,j,k))
+c TYPO HERE ???????????
+c |
+c |
+c phi(i,j,k)= datan2(y(i,j,k),x(i,j,k))
+ enddo
+ enddo
+ enddo
+
+ call CCTK_GetInterpHandle (handle, "simple_local")
+
+ npoints = nx*ny*nz
+
+ call CCTK_Interp (ierror,cctkGH,handle,npoints,2,6,6,
+ $ ne,nq,abseta,q,
+ $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
+ $ etagrd(1),qgrd(1),deta,dq,
+ $ psi2d,detapsi2d,dqpsi2d,detaetapsi2d,detaqpsi2d,dqqpsi2d,
+ $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
+ $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
+ $ psi2dv,detapsi2dv,dqpsi2dv,detaetapsi2dv,detaqpsi2dv,
+ $ dqqpsi2dv,
+ $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
+ $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL)
+
+ psi = psi2dv * exp (-0.5 * eta)
+ detapsi2dv = sign_eta * detapsi2dv
+ detaqpsi2dv = sign_eta * detaqpsi2dv
+
+ do k=1,nz
+ do j=1,ny
+ do i=1,nx
+
+c psi(i,j,k) = psi2dv(i,j,k)*exp(-0.5*eta(i,j,k))
+c detapsi2dv(i,j,k) = sign_eta(i,j,k)*detapsi2dv(i,j,k)
+
+c psix = \partial psi / \partial x / psi
+#include "Development/IDAxiBrillBH/src/psi_1st_deriv.x"
+
+c detaqpsi2dv(i,j,k) = sign_eta(i,j,k)*detaqpsi2dv(i,j,k)
+
+c psixx = \partial^2\psi / \partial x^2 / psi
+#include "Development/IDAxiBrillBH/src/psi_2nd_deriv.x"
+ enddo
+ enddo
+ enddo
+
+c Conformal metric
+c gxx = ...
+
+c Derivatives of the metric
+c dxgxx = 1/2 \partial gxx / \partial x
+
+ do k=1,nz
+ do j=1,ny
+ do i=1,nx
+c THESE WERE ALREADY CALCULATED ABOVE !!!
+c eta(i,j,k) = 0.5d0*dlog(x(i,j,k)**2+y(i,j,k)**2+z(i,j,k)**2)
+c q(i,j,k) = datan2(sqrt(x(i,j,k)**2+y(i,j,k)**2),z(i,j,k))
+c phi(i,j,k) = datan2(y(i,j,k),x(i,j,k))
+#include "Development/IDAxiBrillBH/src/gij.x"
+ enddo
+ enddo
+ enddo
+
+c Curvature
+ kxx = 0.0D0
+ kxy = 0.0D0
+ kxz = 0.0D0
+ kyy = 0.0D0
+ kyz = 0.0D0
+ kzz = 0.0D0
+
+ 111 continue
+c Set ADM mass
+ i = ne-15
+ adm = 0.0
+ do j=2,nq-1
+ adm=adm+(psi2d(i,j)-(psi2d(i+1,j)-psi2d(i-1,j))/deta)*exp(0.5*
+ $ etagrd(i))
+ enddo
+ adm=adm/(nq-2)
+ print *,'ADM mass: ',adm
+ if (CCTK_Equals(initial_lapse,"schwarz")==1) then
+ write (*,*)"Initial with schwarzschild-like lapse"
+ write (*,*)"using alp = (2.*r - adm)/(2.*r+adm)."
+ alp = (2.*r - adm)/(2.*r+adm)
+ endif
+
+ conformal_state = CONFORMAL_METRIC
+
+ deallocate(cc,ce,cw,cn,cs,rhs,psi2d,detapsi2d,dqpsi2d,
+ $ detaetapsi2d,detaqpsi2d,dqqpsi2d,
+ $ etagrd,qgrd,
+ $ eta,abseta,sign_eta,q,phi,psi2dv,detapsi2dv,dqpsi2dv,
+ $ detaetapsi2dv,detaqpsi2dv,dqqpsi2dv)
+
+ return
+ end
+
diff --git a/src/IDAxiBrillBH.F b/src/IDAxiBrillBH.F
new file mode 100644
index 0000000..218e43e
--- /dev/null
+++ b/src/IDAxiBrillBH.F
@@ -0,0 +1,339 @@
+c/*@@
+c @file IDAxiBrillBH.F
+c @date
+c @author
+c @desc
+c
+c @enddesc
+c@@*/
+
+#include "cctk.h"
+#include "cctk_parameters.h"
+#include "cctk_arguments.h"
+
+c/*@@
+c @routine IDAxiBrillBH
+c @date
+c @author
+c @desc
+c
+c @enddesc
+c @calls
+c @calledby
+c @history
+c
+c @endhistory
+c@@*/
+
+c Need include file from Einstein
+#include "CactusEinstein/Einstein/src/Einstein.h"
+
+ subroutine IDAxiBrillBH(CCTK_FARGUMENTS)
+
+ implicit none
+
+ DECLARE_CCTK_FARGUMENTS
+ DECLARE_CCTK_PARAMETERS
+
+c Perhaps this and others should go into cctk.h
+ integer CCTK_Equals
+
+
+ real*8 axibheps, rmax, dq, deta
+ integer levels,id5,id9,idi,idg,ier
+ real*8, allocatable :: cc(:,:),ce(:,:),cw(:,:),cn(:,:),cs(:,:),
+ $ rhs(:,:),psi2d(:,:),detapsi2d(:,:),dqpsi2d(:,:),
+ $ detaetapsi2d(:,:),detaqpsi2d(:,:),dqqpsi2d(:,:)
+ real*8, allocatable :: etagrd(:),qgrd(:)
+ real*8, allocatable :: eta(:,:,:),abseta(:,:,:),sign_eta(:,:,:),
+ $ q(:,:,:),phi(:,:,:)
+ real*8, allocatable :: psi2dv(:,:,:),detapsi2dv(:,:,:),
+ $ dqpsi2dv(:,:,:),detaetapsi2dv(:,:,:),detaqpsi2dv(:,:,:),
+ $ dqqpsi2dv(:,:,:)
+ parameter(axibheps = 1.0d-12)
+ real*8 ep1,ep2
+ real*8 o1,o2,o3,o4,o5,o6,o7,o8,o9
+ real*8 o10,o11,o12,o13,o14,o15,o16,o17,o18,o19
+ real*8 o20,o21,o22,o23,o24,o25,o26,o27,o28,o29
+ real*8 o30,o31,o32,o33,o34,o35,o36,o37,o38,o39
+ real*8 o40,o41,o42,o43,o44,o45,o46,o47,o48,o49
+ real*8 o50,o51,o52,o53,o54,o55,o56,o57,o58,o59
+ real*8 o60,o61,o62,o63,o64,o65,o66,o67,o68,o69
+ real*8 o70,o71,o72,o73,o74,o75,o76,o77,o78,o79
+ real*8 o80,o81,o82,o83,o84,o85,o86,o87,o88,o89
+ real*8 o90,o91,o92,o93,o94,o95,o96,o97,o98,o99
+ real*8 pi
+ real*8 adm
+ integer :: nx,ny,nz
+ integer i,j,k,nquads
+ integer npoints,handle,ierror
+
+ pi = 4.0d0*atan(1.0d0)
+
+c Set up the grid spacings
+ nx = cctk_lsh(1)
+ ny = cctk_lsh(2)
+ nz = cctk_lsh(3)
+
+c Brill wave parameters
+
+ print *,"Brill wave + BH Axisymmetric solve"
+ write (*,123)amp,eta0,sigma,n
+ print *,'etamax=',etamax
+ 123 format(1x,'Pars : Amp',f8.5,' eta0',f8.5,' sigma',f8.5,' n ',i3)
+
+c Solve on this sized cartesian grid
+c 2D grid size NE x NQ
+c Add 2 zones for boundaries...
+ ne = ne+2
+ nq = nq+2
+ ! do I need to call free?
+ allocate(cc(ne,nq),ce(ne,nq),cw(ne,nq),cn(ne,nq),cs(ne,nq),
+ $ rhs(ne,nq),psi2d(ne,nq),detapsi2d(ne,nq),dqpsi2d(ne,nq),
+ $ detaetapsi2d(ne,nq),detaqpsi2d(ne,nq),dqqpsi2d(ne,nq),
+ $ etagrd(ne),qgrd(nq))
+ allocate(eta(nx,ny,nz),abseta(nx,ny,nz),sign_eta(nx,ny,nz),
+ $ q(nx,ny,nz),phi(nx,ny,nz),
+ $ psi2dv(nx,ny,nz),detapsi2dv(nx,ny,nz),dqpsi2dv(nx,ny,nz),
+ $ detaetapsi2dv(nx,ny,nz),detaqpsi2dv(nx,ny,nz),
+ $ dqqpsi2dv(nx,ny,nz))
+c Initialize some arrays
+ psi2d = 1.0d0
+ detapsi2d = 0.0d0
+
+ nquads = 2
+ dq = nquads*0.5*pi/(nq-2)
+ deta = etamax/(ne-3)
+
+ do j=1,nq
+ qgrd(j) = (j-1.5)*dq
+ do i=1,ne
+ etagrd(i) = (i-2)*deta
+#include "Development/IDAxiBrillBH/src/bhbrill.x"
+ enddo
+ enddo
+c Boundary conditions
+ do j=1,nq
+ ce(2,j)=ce(2,j)+cw(2,j)
+ cw(2,j)=0.0
+
+ cw(ne-1,j)=cw(ne-1,j)+ce(ne-1,j)
+ cc(ne-1,j)=cc(ne-1,j)-deta*ce(ne-1,j)
+ ce(ne-1,j)=0.0
+
+ enddo
+ do i=1,ne
+ cc(i,2)=cc(i,2)+cs(i,2)
+ cs(i,2)=0.0
+ cc(i,nq-1)=cc(i,nq-1)+cn(i,nq-1)
+ cn(i,nq-1)=0.0
+ enddo
+
+c Do the solve
+ print *, " Calling axisymmetric solver"
+ call mgparm (levels,5,id5,id9,idi,idg,ne,nq)
+ call mg5 (ne,2,ne-1,nq,2,nq-1,
+ $ cc,cn,cs,cw,ce,psi2d,rhs,
+ $ id5,id9,idi,idg,1,axibheps,rmax,ier)
+ print *, " Solve complete"
+c The solution is now available.
+c Debugging is needed, a stop statement should
+c be called if the IVP solve is not successful
+
+ if(ier .ne. 0) stop 'bad solution to brill wave problem'
+ print *,'rmax = ',rmax
+ print *,'axibheps = ',axibheps
+ print *,'psi2d = ',maxval(psi2d),' ',minval(psi2d)
+
+ ep2 = 0.0
+ do j=2,nq-1
+ do i=2,ne-1
+ ep1 = rhs(i,j)-psi2d(i,j)*cc(i,j)-psi2d(i,j+1)*cn(i,j)-psi2d(i,j-1)*cs(i,j)-
+ & psi2d(i+1,j)*ce(i,j)-psi2d(i-1,j)*cw(i,j)
+ ep2 = max(abs(ep1),ep2)
+ enddo
+ enddo
+ print *,'Resulting eps =',ep2
+
+ ! what a pain all this is....
+ do j=1,nq
+ psi2d(1,j)=psi2d(3,j)
+ psi2d(ne,j)=-deta*psi2d(ne-1,j)+psi2d(ne-2,j)
+ enddo
+ do i=1,ne
+ psi2d(i,1)=psi2d(i,2)
+ psi2d(i,nq)=psi2d(i,nq-1)
+ enddo
+c goto 111
+ do j=2,nq-1
+ do i=2,ne-1
+ dqpsi2d(i,j)=0.5*(psi2d(i,j+1)-psi2d(i,j-1))/dq
+ dqqpsi2d(i,j)=(psi2d(i,j+1)+psi2d(i,j-1)-2.*psi2d(i,j))/dq**2
+ detapsi2d(i,j)=sinh(0.5*etagrd(i))+0.5*(psi2d(i+1,j)-psi2d(i-1,j))/deta
+ detaetapsi2d(i,j)=0.5*cosh(0.5*etagrd(i))+
+ $ (psi2d(i+1,j)+psi2d(i-1,j)-2.*psi2d(i,j))/deta**2
+ enddo
+ enddo
+ do j=1,nq
+ detapsi2d(1,j)=-detapsi2d(3,j)
+ detapsi2d(ne,j)=detapsi2d(ne-2,j) ! simplified
+
+ detaetapsi2d(1,j)=detaetapsi2d(3,j)
+ detaetapsi2d(ne,j)=detaetapsi2d(ne-2,j) ! simplified...
+
+ dqqpsi2d(1,j)=dqqpsi2d(3,j)
+ dqqpsi2d(ne,j)=dqqpsi2d(ne-2,j) ! simplified
+
+ dqpsi2d(1,j)=dqpsi2d(3,j)
+ dqpsi2d(ne,j)=-dq*dqpsi2d(ne-1,j)+dqpsi2d(ne-2,j)
+ enddo
+ do i=1,ne
+ detapsi2d(i,1)=detapsi2d(i,2)
+ detapsi2d(i,nq)=detapsi2d(i,nq-1)
+
+ detaetapsi2d(i,1)=detaetapsi2d(i,2)
+ detaetapsi2d(i,nq)=detaetapsi2d(i,nq-1)
+
+ dqqpsi2d(i,1)=dqqpsi2d(i,2)
+ dqqpsi2d(i,nq)=dqqpsi2d(i,nq-1)
+
+ dqpsi2d(i,1)=-dqpsi2d(i,2)
+ dqpsi2d(i,nq)=-dqpsi2d(i,nq-1)
+ enddo
+ do j=2,nq-1
+ do i=2,ne-1
+ detaqpsi2d(i,j)=0.5*(detapsi2d(i,j+1)-detapsi2d(i,j-1))/dq
+ enddo
+ enddo
+ do j=1,nq
+ detaqpsi2d(1,j)=-detaqpsi2d(3,j)
+ detaqpsi2d(ne,j)=detaqpsi2d(ne-2,j) ! simplified
+ enddo
+ do i=1,ne
+ detaqpsi2d(i,1)=-detaqpsi2d(i,2)
+ detaqpsi2d(i,nq)=-detaqpsi2d(i,nq-1)
+ enddo
+ do j=1,nq
+ psi2d(:,j)=psi2d(:,j)+2.*cosh(0.5*etagrd)
+ enddo
+
+c Now evaluate each of the following at x(i,j,k), y(i,j,k) and
+c z(i,j,k) where i,j,k go from 1 to nx,ny,nz
+
+c Conformal factor
+
+ eta = 0.5d0 * dlog (x**2 + y**2 + z**2)
+ abseta = abs (eta)
+ q = datan2 (sqrt (x**2 + y**2), z)
+ phi = datan2 (y, x)
+
+ do k=1,nz
+ do j=1,ny
+ do i=1,nx
+c eta(i,j,k) = 0.5d0*dlog(x(i,j,k)**2+y(i,j,k)**2+z(i,j,k)**2)
+c abseta(i,j,k) = abs(eta(i,j,k))
+ if(eta(i,j,k) .lt. 0)then
+ sign_eta(i,j,k) = -1
+ else
+ sign_eta(i,j,k) = 1
+ endif
+c q(i,j,k) = datan2(sqrt(x(i,j,k)**2+y(i,j,k)**2),z(i,j,k))
+c TYPO HERE ???????????
+c |
+c |
+c phi(i,j,k)= datan2(y(i,j,k),x(i,j,k))
+ enddo
+ enddo
+ enddo
+
+ call CCTK_GetInterpHandle (handle, "simple_local")
+
+ npoints = nx*ny*nz
+
+ call CCTK_Interp (ierror,cctkGH,handle,npoints,2,6,6,
+ $ ne,nq,abseta,q,
+ $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
+ $ etagrd(1),qgrd(1),deta,dq,
+ $ psi2d,detapsi2d,dqpsi2d,detaetapsi2d,detaqpsi2d,dqqpsi2d,
+ $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
+ $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
+ $ psi2dv,detapsi2dv,dqpsi2dv,detaetapsi2dv,detaqpsi2dv,
+ $ dqqpsi2dv,
+ $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
+ $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL)
+
+ psi = psi2dv * exp (-0.5 * eta)
+ detapsi2dv = sign_eta * detapsi2dv
+ detaqpsi2dv = sign_eta * detaqpsi2dv
+
+ do k=1,nz
+ do j=1,ny
+ do i=1,nx
+
+c psi(i,j,k) = psi2dv(i,j,k)*exp(-0.5*eta(i,j,k))
+c detapsi2dv(i,j,k) = sign_eta(i,j,k)*detapsi2dv(i,j,k)
+
+c psix = \partial psi / \partial x / psi
+#include "Development/IDAxiBrillBH/src/psi_1st_deriv.x"
+
+c detaqpsi2dv(i,j,k) = sign_eta(i,j,k)*detaqpsi2dv(i,j,k)
+
+c psixx = \partial^2\psi / \partial x^2 / psi
+#include "Development/IDAxiBrillBH/src/psi_2nd_deriv.x"
+ enddo
+ enddo
+ enddo
+
+c Conformal metric
+c gxx = ...
+
+c Derivatives of the metric
+c dxgxx = 1/2 \partial gxx / \partial x
+
+ do k=1,nz
+ do j=1,ny
+ do i=1,nx
+c THESE WERE ALREADY CALCULATED ABOVE !!!
+c eta(i,j,k) = 0.5d0*dlog(x(i,j,k)**2+y(i,j,k)**2+z(i,j,k)**2)
+c q(i,j,k) = datan2(sqrt(x(i,j,k)**2+y(i,j,k)**2),z(i,j,k))
+c phi(i,j,k) = datan2(y(i,j,k),x(i,j,k))
+#include "Development/IDAxiBrillBH/src/gij.x"
+ enddo
+ enddo
+ enddo
+
+c Curvature
+ kxx = 0.0D0
+ kxy = 0.0D0
+ kxz = 0.0D0
+ kyy = 0.0D0
+ kyz = 0.0D0
+ kzz = 0.0D0
+
+ 111 continue
+c Set ADM mass
+ i = ne-15
+ adm = 0.0
+ do j=2,nq-1
+ adm=adm+(psi2d(i,j)-(psi2d(i+1,j)-psi2d(i-1,j))/deta)*exp(0.5*
+ $ etagrd(i))
+ enddo
+ adm=adm/(nq-2)
+ print *,'ADM mass: ',adm
+ if (CCTK_Equals(initial_lapse,"schwarz")==1) then
+ write (*,*)"Initial with schwarzschild-like lapse"
+ write (*,*)"using alp = (2.*r - adm)/(2.*r+adm)."
+ alp = (2.*r - adm)/(2.*r+adm)
+ endif
+
+ conformal_state = CONFORMAL_METRIC
+
+ deallocate(cc,ce,cw,cn,cs,rhs,psi2d,detapsi2d,dqpsi2d,
+ $ detaetapsi2d,detaqpsi2d,dqqpsi2d,
+ $ etagrd,qgrd,
+ $ eta,abseta,sign_eta,q,phi,psi2dv,detapsi2dv,dqpsi2dv,
+ $ detaetapsi2dv,detaqpsi2dv,dqqpsi2dv)
+
+ return
+ end
+
diff --git a/src/bhbrill.m b/src/bhbrill.m
new file mode 100644
index 0000000..6c78a1f
--- /dev/null
+++ b/src/bhbrill.m
@@ -0,0 +1,126 @@
+$Path = Union[$Path,{"~/SetTensor"}];
+Needs["SetTensor`"];
+
+Dimension = 3;
+x[1] = eta; x[2] = q; x[3] = phi;
+qf[eta_,q_] := amp (Exp[-(eta-eta0)^2/sigma^2]+Exp[-(eta+eta0)^2/sigma^2]) Sin[q]^n
+
+md = {
+{Exp[2 qf[eta,q]],0,0},
+{0,Exp[2 qf[eta,q]],0},
+{0,0,Sin[q]^2}} psi2d[eta,q]^4;
+InitializeMetric[md];
+
+Clear[exc];
+DefineTensor[exc];
+SetTensor[exc[la,lb],{{0,0,0},{0,0,0},{0,0,0}}];
+
+tmp = RicciR[la,lb] Metricg[ua,ub]+exc[la,lb] Metricg[ua,ub] exc[lc,ld] Metricg[uc,ud]-
+ exc[la,lb] exc[lc,ld] Metricg[ua,uc] Metricg[ub,ud];
+tmp = RicciToAffine[tmp];
+tmp = EvalMT[tmp];
+tmp = ExpandAll[-Exp[2 qf[eta,q]] psi2d[eta,q]^5/8 tmp]
+sav=tmp
+tmp = SubFun[sav,psi2d[eta,q],2 Cosh[eta/2]+psi2d[eta,q]]
+
+(* Make the stencil... *)
+
+stencil = ExpandAll[tmp /. {
+ D[psi2d[eta,q],eta]->(psi2d[i+1,j]-psi2d[i-1,j])/(2 deta),
+ D[psi2d[eta,q],eta,eta]->(psi2d[i+1,j]+psi2d[i-1,j]-2 psi2d[i,j])/(deta deta),
+ D[psi2d[eta,q],q]->(psi2d[i,j+1]-psi2d[i,j-1])/(2 dq),
+ D[psi2d[eta,q],q,q]->(psi2d[i,j+1]+psi2d[i,j-1]-2 psi2d[i,j])/(dq dq),
+ psi2d[eta,q]->psi2d[i,j]
+ }];
+
+cn = Coefficient[stencil,psi2d[i,j+1]]
+cs = Coefficient[stencil,psi2d[i,j-1]]
+ce = Coefficient[stencil,psi2d[i+1,j]]
+cw = Coefficient[stencil,psi2d[i-1,j]]
+cc = Coefficient[stencil,psi2d[i,j]]
+rhs = -SubFun[tmp,psi2d[eta,q],0]
+
+FortranOutputOfDepList = "(i,j)";
+$FortranReplace = Union[{
+ "UND"->"_",
+ "(eta,q)"->"(i,j)"
+}];
+fd = FortranOpen["bhbrill.x"];
+FortranWrite[fd,Cn[i,j],cn ];
+FortranWrite[fd,Cs[i,j],cs ];
+FortranWrite[fd,Cw[i,j],cw ];
+FortranWrite[fd,Cc[i,j],cc ];
+FortranWrite[fd,Ce[i,j],ce ];
+FortranWrite[fd,Rhs[i,j],rhs ];
+FortranClose[fd];
+
+(* Next part, write out conformal g's and d's *)
+
+
+xv = Exp[eta] Sin[q] Cos[phi];
+yv = Exp[eta] Sin[q] Sin[phi];
+zv = Exp[eta] Cos[q];
+
+mc = Table[ D[ {xv,yv,zv}[[i]], {eta,q,phi}[[j]] ],{i,1,3},{j,1,3}];
+mci = Simplify[Inverse[mc]];
+
+Clear[mct];
+DefineTensor[mct,{{1,2},1}];
+Iter[mct[ua,lb],
+ mct[ua,lb]=mc[[ua,-lb]];
+];
+
+Clear[mcti];
+DefineTensor[mcti,{{1,2},1}];
+Iter[mcti[ua,lb],
+ mcti[ua,lb]=mci[[ua,-lb]];
+];
+
+gijtmp = Exp[2 eta]/psi2d[eta,q]^4 Metricg[lc,ld] mcti[uc,la] mcti[ud,lb]
+
+Clear[i2];
+DefineTensor[i2,{{2,1},1}];
+
+fd = FortranOpen["gij.x"];
+Iter[i2[ua,ub],
+ v1 = {x,y,z}[[ua]];
+ v2 = {x,y,z}[[ub]];
+ metv = ToExpression["g"<>ToString[v1]<>ToString[v2]<>"[i,j,k]"];
+ gg[v1,v2]=Simplify[EvalMT[gijtmp,{la->-ua,lb->-ub}]];
+ FortranWrite[fd,metv,gg[v1,v2]];
+ For[ii=1,ii<=3,ii++,
+ v3 = {x,y,z}[[ii]];
+ dmetv = ToExpression["d"<>ToString[v3]<>ToString[metv]];
+ res = OD[gg[v1,v2],lc] mcti[uc,ld]/2;
+ res = EvalMT[res,ld-> -ii];
+ res = Simplify[res];
+ FortranWrite[fd,dmetv,res];
+ ];
+];
+FortranClose[fd];
+
+$FortranReplace = {
+ "UND"->"_",
+ "(eta,q)"->""
+};
+
+fd = FortranOpen["psi_1st_deriv.x"];
+For[ii=1,ii<=3,ii++,
+ v1 = {x,y,z}[[ii]];
+ psv =ToExpression["psi"<>ToString[v1]<>"[i,j,k]"];
+ rhs = CD[Exp[-eta/2] psi2dv[eta,q],lc] mcti[uc,la];
+ rhs = EvalMT[rhs,{la->-ii}]/(Exp[-eta/2] psi2dv[eta,q]);
+ FortranWrite[fd,psv,rhs];
+];
+FortranClose[fd];
+
+fd = FortranOpen["psi_2nd_deriv.x"];
+Iter[i2[ua,ub],
+ v1 = {x,y,z}[[ua]];
+ v2 = {x,y,z}[[ub]];
+ psv = ToExpression["psi"<>ToString[v1]<>ToString[v2]<>"[i,j,k]"];
+ rhs = OD[OD[Exp[-eta/2] psi2dv[eta,q],lc] mcti[uc,la],ld] mcti[ud,lb];
+ rhs = EvalMT[rhs,{la->-ua,lb->-ub}]/(Exp[-eta/2] psi2dv[eta,q]);
+ FortranWrite[fd,psv,rhs];
+];
+FortranClose[fd];
diff --git a/src/bhbrill.x b/src/bhbrill.x
new file mode 100644
index 0000000..b7721c9
--- /dev/null
+++ b/src/bhbrill.x
@@ -0,0 +1,60 @@
+ o1 = dq**2
+ o2 = 1/o1
+ o3 = 1/dq
+ o4 = tan(qgrd(j))
+ o5 = 1/o4
+ o6 = deta**2
+ o7 = 1/o6
+ o8 = sin(qgrd(j))
+ o9 = o8**2
+ o10 = 1/o9
+ o11 = cos(qgrd(j))
+ o12 = o11**2
+ o13 = etagrd(i)**2
+ o14 = sigma**2
+ o15 = 1/o14
+ o16 = -(o13*o15)
+ o17 = -2.00000000000000d0*etagrd(i)*eta0*o15
+ o18 = eta0**2
+ o19 = -(o15*o18)
+ o20 = o16 + o17 + o19
+ o21 = exp(o20)
+ o22 = -2.00000000000000d0 + n
+ o23 = o8**o22
+ o24 = n**2
+ o25 = 2.00000000000000d0*etagrd(i)*eta0*o15
+ o26 = o16 + o19 + o25
+ o27 = exp(o26)
+ o28 = o8**n
+ o29 = o14**2
+ o30 = 1/o29
+ o31 = o4**2
+ o32 = 1/o31
+ o33 = 5.0000000000000d-1*etagrd(i)
+ o34 = cosh(o33)
+ Cn(i,j) = o2 + 5.0000000000000d-1*o3*o5
+ Cs(i,j) = o2 - 5.0000000000000d-1*o3*o5
+ Cw(i,j) = o7
+ Cc(i,j) = -1.25000000000000d-1 - 1.25000000000000d-1*o10 - 2.000
+ & 00000000000d0*o2 - 2.50000000000000d-1*amp*n*o12*o21*o23 + 2.500
+ & 00000000000d-1*amp*o12*o21*o23*o24 - 2.50000000000000d-1*amp*n*o
+ & 12*o23*o27 + 2.50000000000000d-1*amp*o12*o23*o24*o27 - 2.5000000
+ & 0000000d-1*amp*n*o21*o28 - 5.0000000000000d-1*amp*o15*o21*o28 -
+ & 2.50000000000000d-1*amp*n*o27*o28 - 5.0000000000000d-1*amp*o15*o
+ & 27*o28 + 2.00000000000000d0*amp*etagrd(i)*eta0*o21*o28*o30 + amp*o13*o
+ & 21*o28*o30 + amp*o18*o21*o28*o30 - 2.00000000000000d0*amp*etagrd(i)*et
+ & a0*o27*o28*o30 + amp*o13*o27*o28*o30 + amp*o18*o27*o28*o30 + 1.2
+ & 5000000000000d-1*o32 - 2.00000000000000d0*o7
+ Ce(i,j) = o7
+ Rhs(i,j) = -2.50000000000000d-1*o34 + 2.50000000000000d-1*o10*o3
+ & 4 + 5.0000000000000d-1*amp*n*o12*o21*o23*o34 - 5.0000000000000d-
+ & 1*amp*o12*o21*o23*o24*o34 + 5.0000000000000d-1*amp*n*o12*o23*o27
+ & *o34 - 5.0000000000000d-1*amp*o12*o23*o24*o27*o34 + 5.0000000000
+ & 000d-1*amp*n*o21*o28*o34 + amp*o15*o21*o28*o34 + 5.0000000000000
+ & d-1*amp*n*o27*o28*o34 + amp*o15*o27*o28*o34 - 4.0000000000000d0*
+ & amp*etagrd(i)*eta0*o21*o28*o30*o34 - 2.00000000000000d0*amp*o13*o21*o2
+ & 8*o30*o34 - 2.00000000000000d0*amp*o18*o21*o28*o30*o34 + 4.00000
+ & 00000000d0*amp*etagrd(i)*eta0*o27*o28*o30*o34 - 2.00000000000000d0*amp
+ & *o13*o27*o28*o30*o34 - 2.00000000000000d0*amp*o18*o27*o28*o30*o3
+ & 4 - 2.50000000000000d-1*o32*o34
+
diff --git a/src/gij.x b/src/gij.x
new file mode 100644
index 0000000..8b644f2
--- /dev/null
+++ b/src/gij.x
@@ -0,0 +1,133 @@
+ o1 = 2.00000000000000d0*phi(i,j,k)
+ o2 = cos(o1)
+ o3 = -eta0
+ o4 = eta(i,j,k) + o3
+ o5 = o4**2
+ o6 = sigma**2
+ o7 = 1/o6
+ o8 = -(o5*o7)
+ o9 = exp(o8)
+ o10 = eta(i,j,k) + eta0
+ o11 = o10**2
+ o12 = -(o11*o7)
+ o13 = exp(o12)
+ o14 = o13 + o9
+ o15 = sin(q(i,j,k))
+ o16 = o15**n
+ o17 = 2.00000000000000d0*amp*o14*o16
+ o18 = exp(o17)
+ o19 = -1.00000000000000d0 + o18
+ o20 = eta(i,j,k)**2
+ o21 = 2.00000000000000d0*o20
+ o22 = eta0**2
+ o23 = 2.00000000000000d0*o22
+ o24 = eta(i,j,k)*o6
+ o25 = o21 + o23 + o24
+ o26 = -(o25*o7)
+ o27 = exp(o26)
+ o28 = 1/o15
+ o29 = o20 + o22
+ o30 = 2.00000000000000d0*o29*o7
+ o31 = exp(o30)
+ o32 = sin(phi(i,j,k))
+ o33 = sin(o1)
+ o34 = o11*o7
+ o35 = exp(o34)
+ o36 = cos(phi(i,j,k))
+ o37 = o32**2
+ o38 = cos(q(i,j,k))
+ o39 = o38**2
+ o40 = o5*o7
+ o41 = exp(o40)
+ o42 = o15**2
+ o43 = o36**2
+ o44 = -2.00000000000000d0*o19*o31*o6
+ o45 = -2.00000000000000d0*eta(i,j,k)
+ o46 = -2.00000000000000d0*eta0
+ o47 = n*o6
+ o48 = 4.0000000000000d0*eta(i,j,k)*eta0*o7
+ o49 = exp(o48)
+ o50 = -2.00000000000000d0*eta(i,j,k)*o49
+ o51 = 2.00000000000000d0*eta0*o49
+ o52 = n*o49*o6
+ o53 = 2.00000000000000d0*q(i,j,k)
+ o54 = cos(o53)
+ o55 = 2.00000000000000d0*eta(i,j,k)
+ o56 = 2.00000000000000d0*eta0
+ o57 = 2.00000000000000d0*eta(i,j,k)*o49
+ o58 = -2.00000000000000d0*eta0*o49
+ o59 = o47 + o52 + o55 + o56 + o57 + o58
+ o60 = o54*o59
+ o61 = o45 + o46 + o47 + o50 + o51 + o52 + o60
+ o62 = o17 + o40
+ o63 = exp(o62)
+ o64 = amp*o16*o61*o63
+ o65 = o44 + o64
+ o66 = -1.00000000000000d0 + o49
+ o67 = -2.00000000000000d0*eta0*o66
+ o68 = 1.00000000000000d0 + o49
+ o69 = 2.00000000000000d0*eta(i,j,k)*o68
+ o70 = n*o6*o68
+ o71 = o67 + o69 + o70
+ o72 = o56 + o6
+ o73 = eta(i,j,k)*o72
+ o74 = o20 + o22 + o73
+ o75 = -(o7*o74)
+ o76 = o17 + o75
+ o77 = exp(o76)
+ o78 = -eta(i,j,k)
+ o79 = exp(o78)
+ o80 = -2.00000000000000d0*o29*o7
+ o81 = o17 + o80
+ o82 = exp(o81)
+ o83 = -1.00000000000000d0 + n
+ o84 = o15**o83
+ o85 = o35 + o41
+ o86 = -(n*o39*o6*o85)
+ o87 = eta(i,j,k)*o49
+ o88 = -(eta0*o49)
+ o89 = eta(i,j,k) + eta0 + o87 + o88
+ o90 = 2.00000000000000d0*o41*o42*o89
+ o91 = o86 + o90
+ o92 = o17 + o26
+ o93 = exp(o92)
+ o94 = n*o39*o6*o85
+ o95 = -2.00000000000000d0*o41*o42*o89
+ o96 = o94 + o95
+ gxx(i,j,k) = 5.0000000000000d-1*(1.00000000000000d0 + o18 + o19*
+ & o2)
+c dxgxx(i,j,k) = 5.0000000000000d-1*o27*o28*o7*(o36*(-2.0000000000
+c & 0000d0*o31*o37*o6 - amp*o16*o18*(2.00000000000000d0*eta(i,j,k)*o35*o42
+c & - 3.00000000000000d0*eta0*o35*o42 + 2.00000000000000d0*eta(i,j,k)*o2*o3
+c & 5*o42 + 2.00000000000000d0*eta(i,j,k)*o41*o42 + 2.00000000000000d0*eta0
+c & *o41*o42 + 2.00000000000000d0*eta(i,j,k)*o2*o41*o42 + 2.00000000000000d
+c & 0*eta0*o2*o41*o42 - n*o35*o39*o6 - n*o2*o35*o39*o6 - n*o39*o41*o
+c & 6 - n*o2*o39*o41*o6)) + o18*(o31*o32*o33*o6 + amp*eta0*o15**(2.0
+c & 0000000000000d0 + n)*o35*cos(3.00000000000000d0*phi(i,j,k))))
+c dygxx(i,j,k) = 5.0000000000000d-1*o27*o28*o32*o43*o65*o7
+c dzgxx(i,j,k) = -(amp*o16*o38*o43*o7*o71*o77)
+ gxy(i,j,k) = 5.0000000000000d-1*o19*o33
+c dxgxy(i,j,k) = 5.0000000000000d-1*o7*o79*(-(o19*o2*o28*o32*o6) -
+c & amp*o33*o36*o82*o84*o91)
+c dygxy(i,j,k) = 5.0000000000000d-1*o7*o79*(1.00000000000000d0*o19
+c & *o2*o28*o36*o6 - amp*o32*o33*o82*o84*o91)
+c dzgxy(i,j,k) = -5.0000000000000d-1*amp*o16*o33*o38*o7*o71*o77
+ gxz(i,j,k) = 0
+c dxgxz(i,j,k) = 0
+c dygxz(i,j,k) = 0
+c dzgxz(i,j,k) = 0
+ gyy(i,j,k) = 5.0000000000000d-1*(1.00000000000000d0 + o18 - o19*
+ & o2)
+c dxgyy(i,j,k) = 5.0000000000000d-1*o27*o28*o36*o37*o65*o7
+c dygyy(i,j,k) = 5.0000000000000d-1*o27*o28*o7*(-2.00000000000000d
+c & 0*o31*o32*o43*o6 - 2.00000000000000d0*amp*o16*o18*o32*o37*o91 +
+c & o33*o36*o6*exp(2.00000000000000d0*(amp*o14*o16 + o29*o7)))
+c dzgyy(i,j,k) = -(amp*o16*o37*o38*o7*o71*o77)
+ gyz(i,j,k) = 0
+c dxgyz(i,j,k) = 0
+c dygyz(i,j,k) = 0
+c dzgyz(i,j,k) = 0
+ gzz(i,j,k) = o18
+c dxgzz(i,j,k) = amp*o36*o7*o84*o93*o96
+c dygzz(i,j,k) = amp*o32*o7*o84*o93*o96
+c dzgzz(i,j,k) = -(amp*o16*o38*o7*o71*o77)
diff --git a/src/interp2.F b/src/interp2.F
new file mode 100644
index 0000000..006287f
--- /dev/null
+++ b/src/interp2.F
@@ -0,0 +1,213 @@
+c /*@@
+c@routine interp2d
+c@date Fri Feb 14 08:46:53 1997
+c@author Paul Walker
+c@desc
+c Interpolates from 2D data var with coordinates x and y and
+c sizes nx and ny onto 1D data out with position outx and outy
+c and nout points.
+c <p>
+c This has linear, quadratic and cubic interpolators in it.
+c Or will one day very soon.
+c@enddesc
+c@calls
+c@calledby numerical_axig
+c@@*/
+
+ subroutine interp2d(var,x,y,nx,ny,out,outx,outy,nout,order)
+ implicit none
+ integer nx,ny,nout
+ real*8 var(nx,ny), x(nx), y(ny)
+ real*8 out(nout),outx(nout),outy(nout)
+ integer order
+c Interpolation goes from ibelow to ibelow+[1,2,3] depending on order
+ integer i,j,ibelow,jbelow,pt
+ real*8 xsym,ysym,findx,findy,frac
+ real*8 ydir(order+1)
+ real*8 ft(10), xt(10)
+ real*8 poly2inter, quad_2d, cubic_2d
+ real*8 dx, dy, PI
+ integer twobhjsad
+
+ PI = 3.14159265
+
+c Set up the grid spacings
+ dx = x(2) - x(1)
+ dy = y(2) - y(1)
+
+c Loop over all out points
+ do pt=1,nout
+ ysym = 1.0D0
+ xsym = 1.0D0
+c Check bounds
+ findx = outx(pt)
+ if (findx .lt. x(1)) then
+ write (*,*) "Below inner bound at ",pt,outx(pt)
+ STOP
+ endif
+ if (findx .gt. x(nx)) then
+ write (*,*) "Above x bounds at ",pt,outx(pt),x(nx)
+ STOP
+ endif
+ findy = outy(pt)
+ if (findy .lt. y(1)) then
+ write (*,*) "Below y inner bound at ",pt,outy(pt)
+ STOP
+ endif
+ if (findy .gt. y(ny)) then
+ write (*,*) "Below y inner bound at ",pt,outy(pt)
+ STOP
+ endif
+
+c Locate ourselves in i,j space
+c do i=1,nx
+c if (x(i) .lt. findx) then
+c ibelow = i
+c endif
+c enddo
+c Assume a regular grid
+ ibelow = (findx-x(1))/dx+1
+
+ if (order .eq. 3 .and. ibelow .gt. 1) then
+ ibelow = ibelow - 1
+ endif
+ if (ibelow + order .gt. nx) then
+ ibelow = nx - order
+ endif
+
+c do i=1,ny
+c if (y(i) .lt. findy) then
+c jbelow = i
+c endif
+c enddo
+c Assume a regular grid
+ jbelow = (findy-y(1))/dy+1
+ if (order .eq. 3 .and. jbelow .gt. 1) then
+ jbelow = jbelow - 1
+ endif
+ if (jbelow + order .gt. ny) then
+ jbelow = ny - order
+ endif
+
+c write (*,*) "PT :",findx,findy
+c write (*,*) "SYM:",sym
+c write (*,*) "BOUND X ",ibelow,x(ibelow),x(ibelow+1)
+c write (*,*) "BOUND Y ",jbelow,y(jbelow),y(jbelow+1)
+
+
+c So do the interpolation
+ if (order .eq. 1) then
+c Interp in the x direction
+ frac = (findx-x(ibelow))/(x(ibelow+1)-x(ibelow))
+ ydir(1) = frac * var(ibelow+1,jbelow) +
+ $ (1.0 - frac)*var(ibelow,jbelow)
+ ydir(2) = frac * var(ibelow+1,jbelow+1) +
+ $ (1.0 - frac)*var(ibelow,jbelow+1)
+c And now in the y direction
+ frac = (findy-y(jbelow))/(y(jbelow+1)-y(jbelow))
+ out(pt) = xsym * ysym *
+ $ (frac * ydir(2) + (1.0 - frac) * ydir(1))
+ else if (order .eq. 2) then
+c Load up for calls to poly2inter
+ do j=1,3
+ do i=1,3
+ ft(i) = var(ibelow+i-1,jbelow+j-1)
+ xt(i) = x(ibelow+i-1)
+ enddo
+ ydir(j) = quad_2d(ft,xt(1),dx,findx)
+ enddo
+ do j=1,3
+ xt(j) = y(jbelow+j-1)
+ enddo
+ out(pt) = xsym * ysym*quad_2d(ydir,xt(1),dy,findy)
+ else if (order .eq. 3) then
+c Load up for calls to cubic
+ do j=1,4
+ do i=1,4
+ ft(i) = var(ibelow+i-1,jbelow+j-1)
+ xt(i) = x(ibelow+i-1)
+ enddo
+ ydir(j) = cubic_2d(ft,xt(1),dx,findx)
+ enddo
+ do j=1,4
+ xt(j) = y(jbelow+j-1)
+ enddo
+ out(pt) = xsym * ysym*cubic_2d(ydir,xt(1),dy,findy)
+ else
+ write (*,*) "ORDER set wrong in interp2d",order
+ stop
+ endif
+ enddo
+
+ return
+ end
+
+ real*8 function linear_2d(f, x0, dx, findx)
+ implicit none
+ real*8 f(2),x0,dx,findx
+ real*8 frac
+
+ frac = (findx-x0)/dx
+ linear_2d = (frac)*f(2) + (1.0-frac)*f(1)
+
+ return
+ end
+
+ real*8 function quad_2d(f, x0, dx, findx)
+ implicit none
+ real*8 f(3),x0, dx, findx
+ real*8 f0,f1,f2
+ real*8 a,b,c, dx2, x02, o2dx2
+c Mathematica tells us
+c - List(List(Rule(c,(2*dx**2*f0 + 3*dx*f0*x0 - 4*dx*f1*x0
+c - + dx*f2*x0 + f0*x0**2 - 2*f1*x0**2 + f2*x0**2)
+c - /(2*dx**2)), Rule(b,(-3*dx*f0 + 4*dx*f1 - dx*f2 - 2*f0*x0
+c - + 4*f1*x0 - 2*f2*x0)/(2*dx**2)),Rule(a,(f0 - 2*f1 +
+c - f2)/(2*dx**2))))
+
+ f0 = f(1)
+ f1 = f(2)
+ f2 = f(3)
+ dx2 = dx**2
+ x02 = x0**2
+ o2dx2 = 1.0D0/(2.0D0*dx2)
+
+ c = (2.0D0*dx2*f0 + dx*x0*(3.0D0*f0 - 4.0D0*f1 + f2) +
+ $ x02*(f0 - 2.0D0*f1 + f2))*o2dx2
+ b = (dx * (-3.0D0*f0 + 4.0D0*f1 - f2) + x0 * (- 2.0D0*f0 +
+ $ 4.0D0*f1 - 2.0D0*f2))*o2dx2
+
+ a = (f0 - 2.0D0*f1 + f2)*o2dx2
+
+ quad_2d = a*findx**2 + b*findx + c
+
+ return
+ end
+
+ real*8 function cubic_2d(f, x0, dx, findx)
+ implicit none
+ real*8 a,b,c,d
+ real*8 f(4),x0,dx,findx
+
+ a = -(f(1)-3.0*f(2)+3.0*f(3)-f(4)) / (6.0*(dx**3))
+
+ b = (f(1)-2.0*f(2)+f(3))/(2.0*(dx**2)) +
+ $ (f(1)-3.0*f(2)+3.0*f(3)-f(4))*(dx+x0)/(2.0*(dx**3))
+
+ c = ((dx**2)*(-11.0*f(1) + 18.0*f(2) - 9.0*f(3) + 2.0*f(4)) +
+ $ dx*x0* (-12.0*f(1) + 30.0*f(2) - 24.0*f(3) + 6.0*f(4)) +
+ $ (x0**2)*( -3.0*f(1) + 9.0*f(2) - 9.0*f(3) + 3.0*f(4))) /
+ $ (6.0*(dx**3))
+
+ d = ((dx**3)* ( 6.0*f(1) ) +
+ $ (dx**2)*x0*( 11.0*f(1) - 18.0*f(2) + 9.0*f(3) - 2.0*f(4)) +
+ $ (x0**2)*dx*( 6.0*f(1) - 15.0*f(2) + 12.0*f(3) - 3.0*f(4)) +
+ $ (x0**3)* ( 1.0*f(1) - 3.0*f(2) + 3.0*f(3) - 1.0*f(4)))/
+ $ (6.0*(dx**3))
+
+ cubic_2d = ((a*findx + b)*findx + c)*findx + d
+
+ return
+ end
+
+
diff --git a/src/make.code.defn b/src/make.code.defn
new file mode 100644
index 0000000..7bd9c1e
--- /dev/null
+++ b/src/make.code.defn
@@ -0,0 +1,9 @@
+# Main make.code.defn file for thorn AxiBrillBHIVP
+# $Header$
+
+# Source files in this directory
+SRCS = IDAxiBrillBH.F mg59p.F shmgp.F77
+
+# Subdirectories containing source files
+SUBDIRS =
+
diff --git a/src/mg59p.F b/src/mg59p.F
new file mode 100644
index 0000000..63dcb4a
--- /dev/null
+++ b/src/mg59p.F
@@ -0,0 +1,896 @@
+c----------------------------------------------------------------------
+ subroutine mgparm(m,ifd59,id5,id9,idi,idg,imx,jmx)
+ implicit real*8 (a-h,o-z)
+c----------------------------------------------------------------------
+c Given imx, jmx and ifd59 (See comments in mgsu2), mgparm calculates
+c the number of grids that will be needed, and the dimensions
+* needed for the coefficient, right hand side and solution arrays
+* to store values on all grid levels.
+c .....................................................................
+*
+**** Parameters
+*
+* INTEGERS:
+* m :
+* This is the number of grid levels that the multigrid
+* routine will use.
+*
+* ifmg :
+* ifmg = 0 - The full multigrid algorithm is not used to
+* obtain a good initial guess on the fine grid.
+* (use this if you can provide a good initial guess)
+* ifmg = 1 - The full multigrid algorithm is used to obtain a
+* good initial guess on the fine grid.
+*
+* id5 :
+* Dimension of the arrays ac,aw,as,ae,an,q and f. id5 is the
+* total number of grid points on the finest grid and all
+* coarser grids.
+*
+* id9 :
+* Dimension of the arrays asw,ase,ane,anw. If ifd59=5 then
+* id9=idi. If ifd59=9 then id9=id5.
+* (NOTE: This routine specifically written for 9-point)
+*
+* idi :
+* Dimension of the work arrays pu and pd. idi is the total
+* number of grid points on all of the coarser grids.
+*
+* idg :
+* Dimension of the work array gam. It is set to the value im,
+* the number of grid points in the i-direction on the finest
+* grid.
+*
+* imx,jmx :
+* These are the number of points in the i and j directions
+* including the two ficticious points.
+*
+ parameter(n5=1,n9=2,ni=3,jm=4,i9=5,j9=6,ifd=7,jred=8)
+ dimension np2(20,8)
+ iq5=1
+ iq9=1
+ iqi=1
+ m=1
+ np2(m,1)=jmx
+ np2(m,2)=3
+ 10 if(np2(m,1).le.3) go to 20
+ m=m+1
+ np2(m,1)=np2(m-1,1)/2+1
+ if(np2(m-1,2).eq.2.and.mod(np2(m-1,1),2).eq.1)
+ + np2(m,1)=np2(m,1)+1
+ np2(m,2)=2
+ go to 10
+ 20 do 30 k=1,m
+ np2(m-k+1,jm)=np2(k,1)
+ 30 np2(m-k+1,jred)=np2(k,2)
+ do 40 k=m,1,-1
+ ktot=imx*np2(k,jm)
+ np2(k,n5)=iq5
+ iq5=iq5+ktot
+ np2(k,n9)=iq9
+ if(k.lt.m.or.ifd59.eq.9) iq9=iq9+ktot
+ np2(k,ni)=iqi
+ 40 if(k.lt.m) iqi=iqi+ktot
+ do 50 k=1,m
+ np2(k,i9)=imx
+ np2(k,j9)=np2(k,jm)
+ 50 np2(k,ifd)=9
+ if(ifd59.eq.5) then
+ np2(m,i9)=1
+ np2(m,j9)=1
+ np2(m,ifd)=5
+ endif
+ id5=iq5-1
+ id9=iq9-1
+ idi=iqi-1
+ idg=imx
+ return
+ end
+*
+*
+*
+************************************************************************
+*
+ subroutine mg9 (idim,ilower,iupper,jdim,jlower,jupper,
+ & cc,cn,cs,cw,ce,cnw,cne,csw,cse,u,rhs,
+ & id5,id9,idi,idg,ifmg,eps,rmax,ier)
+ implicit real*8 (a-h,o-z)
+*
+************************************************************************
+*
+*
+* This routine is a wrapper for the multigrid solver. The input
+* arrays are the finite difference coefficients on the 2D grid with
+* the boundary conditions absorbed into them.
+*
+**** Parameters
+*
+* INTEGERS:
+* idim,jdim :
+* These define sizes of the coefficient arrays passed to this
+* routine.
+*
+* ilower,iupper,
+* jlower,jupper :
+* These are the indices of the computational grid that
+* correspond to upper and lower index limits for points on
+* which computations are actually done.
+*
+* id5 :
+* Dimension of the arrays ac,aw,as,ae,an,q and f. id5 is the
+* total number of grid points on the finest grid and all
+* coarser grids.
+*
+* id9 :
+* Dimension of the arrays asw,ase,ane,anw. If ifd59=5 then
+* id9=idi. If ifd59=9 then id9=id5.
+* (NOTE: This routine specifically written for 9-point)
+*
+* idi :
+* Dimension of the work arrays pu and pd. idi is the total
+* number of grid points on all of the coarser grids.
+*
+* idg :
+* Dimension of the work array gam. It is set to the value im,
+* the number of grid points in the i-direction on the finest
+* grid.
+*
+* ifmg :
+* ifmg = 0 - The full multigrid algorithm is not used to
+* obtain a good initial guess on the fine grid.
+* (use this if you can provide a good initial guess)
+* ifmg = 1 - The full multigrid algorithm is used to obtain a
+* good initial guess on the fine grid.
+*
+* ier :
+* This is an error flag with the following possible return
+* values:
+* ier = 0 => solver converged without error
+* ier = -1 => solver did not converge
+*
+* REALS:
+* cc(idim,jdim),cn(idim,jdim),
+* cs(idim,jdim),cw(idim,jdim),ce(idim,jdim),
+* cnw(idim,jdim),cne(idim,jdim),
+* csw(idim,jdim),cse(idim,jdim) :
+* These are the finite difference coefficient arrays for a
+* nine-point stencil on the two dimensional grid as follows:
+*
+*
+* cnw cn cne
+* ^
+* | cw cc ce
+* increasing |
+* j values | csw cs cse
+* (theta) |
+* --------> increasing i values (eta)
+*
+* Ie. the coresspondence is : nw : i-1,j+1
+* n : i,j+1
+* ne : i+1,j+1
+* w : i-1,j
+* c : i,j
+* e : i+1,j
+* sw : i-1,j-1
+* s : i,j-1
+* se : i+1,j-1
+*
+* u :
+* Input: this contains the initial guess to the solution of the
+* equation
+* Output: This contains the final approximation to the solution
+* determined by the multigrid solver.
+*
+* rhs :
+* This array contains the values of the right hand side of the
+* equation at every point on the rwo dimensional grid.
+*
+* eps :
+* eps > 0. The maximum norm of the residual is calculated at
+* the end of each multigrid cycle. The algorithm is
+* terminated when this maximum becomes less than eps
+* or when the maximum number of iterations (see ncyc)
+* is exceeded. It is up to the user to provide a
+* meaningfull tolerance criteria for the particular
+* problem being solved.
+*
+* rmax:
+* This is the final value of the residual calcu;ated.
+*
+************************************************************************
+*
+*
+
+ integer idim,ilower,iupper,jdim,jlower,jupper,ifmg
+ integer id5,id9,idi,idg,ier
+ real*8 cc(idim,jdim),cn(idim,jdim),cs(idim,jdim),cw(idim,jdim),
+ & ce(idim,jdim),cnw(idim,jdim),cne(idim,jdim),csw(idim,jdim),
+ & cse(idim,jdim),u(idim,jdim),rhs(idim,jdim)
+ real*8 eps,rmax
+
+*
+************************************************************************
+*
+* Variable definitions:
+*
+* Integers:
+* ifd59 :
+* ifd59 = 5 - means a 5-point finite difference stencil
+* (ac,an,as,aw,ae) is defined on the finest grid.
+* ifd59 = 9 - means a 9-point finite difference stencil
+* (ac,an,as,aw,ae,anw,ane,asw,ase) is defined on
+* the finest grid by the user.
+* (NOTE: This routine specifically written
+* for 9-point)
+*
+* ncyc :
+* The maximum number of multigrid "v"-cycles to be used. If
+* the maximum norm of the residual is not less than tol at the
+* end of ncyc cycles, the algorithm is terminated.
+* (NOTE: ncyc <= 40 )
+*
+* np(20,8) :
+* Input: When the iskip=1,-1 or -2 option is used, np2 is
+* assumed to contain the grid information for umgs2.
+* Output: When the iskip=0 option is used, the grid
+* information for umgs2 is returned in np2.
+* (NOTE: This is only useful for multiple instance problems)
+*
+* iskip :
+* iskip = 0 - The coarse grid information, coarse grid
+* operators and interpolation coefficients are
+* calculated by umgs2. This information is stored
+* in the arrays ac, aw, as, asw, ase, pu, pd, np2
+* and the variable.
+* iskip = 1 - The calculation of the coarse grid information,
+* coarse grid operators and interpolation
+* coefficients is skipped. This option would be
+* used when umgs2 has been called with iskip=0 and
+* is being called again to solve a system of
+* equations with the same matrix. This would be
+* the case in, say, parabolic problems with time
+* independent coefficients.
+* iskip =-1 - The set up of pointers (ugrdfn) is skipped.
+* Coarse grid operators and interpolation
+* coefficients are calculated and the given matrix
+* equation is solved. This option would be used
+* when umgs2 has been called with iskip=0 and is
+* being called again to solve a system of equations
+* with a different matrix of the same dimensions.
+* This would be the case for, say, parabolic
+* problems with time dependent coefficients.
+* iskip =-2 - The set up of pointers (ugrdfn) is skipped.
+* Coarse grid operators and interpolation
+* coefficients are calculated and returned.
+* No matrix solve.
+*
+* ipc :
+* ipc = 0 or 1.
+* ipc is a multigrid parameter which determines the type of
+* interpolation to be used. Usually ipc=1 is best. However,
+* if the boundary condition equations have been absorbed into
+* the interior equations then ipc=0 can be used which results
+* in a slightly more efficient algorithm.
+*
+* nman :
+* nman = 0 usually.
+* nman = 1 signals that the fine grid equations are singulari
+* for the case when homogeneous Neumann boundary conditions are
+* applied along the entire boundary. In this case, the
+* difference equations are singular and the condition that the
+* integral of q over the domain be zero is added to the set of
+* difference equations. This condition is satisfied by adding
+* the appropriate constant vector to q on the fine grid. It is
+* assumed, in this case, that a well-defined problem has been
+* given to mgss2, i.e. the integral of f over the domain is
+* zero.
+*
+* im :
+* The number of grid points in the x-direction (including two
+* ficticious points)
+* jm :
+* The number of grid points in the y-direction (including two
+* ficticious points)
+*
+* linp :
+* This is a dummy argument left over from the authors
+* development of the code
+* Use: common /io/ linp,lout
+*
+* lout :
+* lout = unit number of output file into which the maximum norm
+* of the residual after each multigrid v-cycle" is printed.
+* Use: common /io/ linp,lout
+*
+* iscale :
+* Flag to indicate whether problem can be diagonally scaled to
+* speed convergence of the multigrid solver.
+*
+* REALS:
+*
+* ac(id5),an(id5),as(id5),aw(id5),ae(id5),
+* anw(id9),ane(id9),asw(id9),ase(id9) :
+* Input: ac, an, as, aw, ae, anw, ane, asw and ase contain the
+* stencil coefficients for the difference operator on
+* the finest grid. When the iskip=1 option is used,
+* these arrays also are assumed to contain the coarse
+* grid difference stencil coeficients.
+* Output: when the iskip=0 option is used, the coarse grid
+* stencil coeficients are returned in ac, an, as, aw,
+* ae, anw, ane, asw and ase.
+*
+* ru(idi),rd(idi),rc(idi) :
+* Real work arrays.
+*
+* pu(idi),pd(idi),pc(idi) :
+* Real work arrays.
+* Input: when the iskip=1 option is used, these arrays are
+* assumed to contain the interpolation coefficients used
+* in the semi-coarsening multigrid algorithm.
+* Output: when the iskip=0 option is used, the interpolation
+* coeficients are returned in pu and pd.
+*
+* f(id5) :
+* f contains the right hand side vector of the matrix
+* equation to be solved by umgs2.
+*
+* q(id5) :
+* If ifmg=0, q contains the initial guess on the fine grid.
+* If ifmg=1, the initial guess on the fine grid is determined
+* by the full multigrid process and the value of
+* q on input to umgs2 not used.
+*
+* tol :
+* tol > 0. The maximum norm of the residual is calculated at
+* the end of each multigrid cycle. The algorithm is
+* terminated when this maximum becomes less than tol
+* or when the maximum number of iterations (see ncyc)
+* is exceeded. It is up to the user to provide a
+* meaningfull tolerance criteria for the particular
+* problem being solved.
+* tol = 0. Perform ncyc multigrid cycles. Calculate and print
+* the maximum norm of the residual after each cycle.
+* tol =-1. Perform ncyc multigrid cycles. The maximum norm of
+* the final residual is calculated and returned in
+* the variable rmax in the calling list of umgs2.
+* tol =-2. Perform ncyc multigrid cycles. The maximum norm of
+* the residual is not calculated.
+*
+* rmax :
+* If tol.ge.-1., the final residual norm is returned in rmax.
+*
+************************************************************************
+*
+*
+
+ integer ncyc,ifd59
+ parameter (ncyc=40,ifd59=9)
+
+ integer np2(20,8)
+ integer iskip,ipc,nman
+ parameter (iskip=0,ipc=1,nman=0)
+ integer irc,irurd,im,jm
+ integer linp,lout
+ common /io/ linp,lout
+
+ real*8 ac(id5),an(id5),as(id5),aw(id5),ae(id5),
+ & anw(id9),ane(id9),asw(id9),ase(id9),
+ & q(id5),f(id5),gam(idg)
+
+ real*8 ru(idi),rd(idi),rc(idi),pu(idi),pd(idi),pc(idi)
+
+ real*8 tol
+
+ integer iscale,itry
+
+
+ ier=0
+
+*
+* Set some parameters for multigrid solver
+*
+ lout=6
+c rewind(unit=lout)
+ irc=0
+ irurd=0
+ im=iupper-ilower+3
+ jm=jupper-jlower+3
+ tol=eps
+
+*
+* Set up coefficients into vectors with correct indexing
+*
+ do 110 j=jlower,jupper
+ do 100 i=ilower,iupper
+ n=(j+1-jlower)*im + i+1-(ilower-1)
+ ac(n)=cc(i,j)
+ an(n)=cn(i,j)
+ as(n)=cs(i,j)
+ aw(n)=cw(i,j)
+ ae(n)=ce(i,j)
+ anw(n)=cnw(i,j)
+ ane(n)=cne(i,j)
+ asw(n)=csw(i,j)
+ ase(n)=cse(i,j)
+ q(n)=u(i,j)
+ f(n)=rhs(i,j)
+ 100 continue
+ 110 continue
+
+
+*
+* Determine whether we can diagonal scale the problem to speed
+* convergence. Can only be done if there are no zeros on the main
+* diagonal (ie. central difference coefficient).
+*
+ iscale=1
+ do 200 j=jlower,jupper
+ do 205 i=ilower,iupper
+ n=(j+1-jlower)*im + i+1-(ilower-1)
+ if (ac(n) .eq. 0.) then
+ iscale=0
+ endif
+ 205 continue
+ 200 continue
+
+*
+* Do the diagonal scaling if we can.
+*
+ if (iscale.eq.1) then
+
+ do 210 j=jlower,jupper
+ do 215 i=ilower,iupper
+ n=(j+1-jlower)*im + i+1-(ilower-1)
+ f(n)=f(n)/ac(n)
+ ase(n)=ase(n)/ac(n)
+ asw(n)=asw(n)/ac(n)
+ ane(n)=ane(n)/ac(n)
+ anw(n)=anw(n)/ac(n)
+ ae(n)=ae(n)/ac(n)
+ aw(n)=aw(n)/ac(n)
+ as(n)=as(n)/ac(n)
+ an(n)=an(n)/ac(n)
+ ac(n)=1.
+ 215 continue
+ 210 continue
+
+ endif
+
+c
+c Now call the multigrid routine
+
+ itry=0
+
+ 1122 call umgs2(
+ . ac,aw,as,ae,an,asw,ase,ane,anw,q,f,pu,pd,pc,ru,rd,rc,gam,np2,
+ . ifd59,ifmg,ncyc,tol,nman,im,jm,id5,id9,idi,m,iskip,rmax,
+ . ipc,irc,irurd)
+
+ if ((rmax.gt.tol).and.(itry.le.5)) then
+ itry=itry+1
+ print*,"Retry #",itry
+ goto 1122
+ endif
+
+ if (rmax.gt.tol) then
+ ier = -1
+ print*,"Did not converge"
+ print*," maximum residual = ",rmax
+ print*," tolerance = ",tol
+ endif
+
+*
+* Convert the solution back to the 2D array form
+*
+ do 510 j=jlower,jupper
+ do 500 i=ilower,iupper
+ n=(j+1-jlower)*im + i+1-(ilower-1)
+ u(i,j)=q(n)
+ 500 continue
+ 510 continue
+
+ return
+ end
+
+
+
+*
+*
+*
+************************************************************************
+*
+ subroutine mg5 (idim,ilower,iupper,jdim,jlower,jupper,
+ & cc,cn,cs,cw,ce,u,rhs,
+ & id5,id9,idi,idg,ifmg,eps,rmax,ier)
+ implicit real*8 (a-h,o-z)
+*
+************************************************************************
+*
+*
+* This routine is a wrapper for the multigrid solver. The input
+* arrays are the finite difference coefficients on the 2D grid with
+* the boundary conditions absorbed into them.
+*
+**** Parameters
+*
+* INTEGERS:
+* idim,jdim :
+* These define sizes of the coefficient arrays passed to this
+* routine.
+*
+* ilower,iupper,
+* jlower,jupper :
+* These are the indices of the computational grid that
+* correspond to upper and lower index limits for points on
+* which computations are actually done.
+*
+* id5 :
+* Dimension of the arrays ac,aw,as,ae,an,q and f. id5 is the
+* total number of grid points on the finest grid and all
+* coarser grids.
+*
+* id9 :
+* Dimension of the arrays asw,ase,ane,anw. If ifd59=5 then
+* id9=idi. If ifd59=9 then id9=id5.
+* (NOTE: This routine specifically written for 5-point)
+*
+* idi :
+* Dimension of the work arrays pu and pd. idi is the total
+* number of grid points on all of the coarser grids.
+*
+* idg :
+* Dimension of the work array gam. It is set to the value im,
+* the number of grid points in the i-direction on the finest
+* grid.
+*
+* ifmg :
+* ifmg = 0 - The full multigrid algorithm is not used to
+* obtain a good initial guess on the fine grid.
+* (use this if you can provide a good initial guess)
+* ifmg = 1 - The full multigrid algorithm is used to obtain a
+* good initial guess on the fine grid.
+*
+* ier :
+* This is an error flag with the following possible return
+* values:
+* ier = 0 => solver converged without error
+* ier = -1 => solver did not converge
+*
+* REALS:
+* cc(idim,jdim),cn(idim,jdim),
+* cs(idim,jdim),cw(idim,jdim),ce(idim,jdim):
+* These are the finite difference coefficient arrays for a
+* five-point stencil on the two dimensional grid as follows:
+*
+*
+* cn
+* ^
+* | cw cc ce
+* increasing |
+* j values | cs
+* (theta) |
+* --------> increasing i values (eta)
+*
+* Ie. the coresspondence is : n : i,j+1
+* ne : i+1,j+1
+* c : i,j
+* e : i+1,j
+* s : i,j-1
+*
+* u :
+* Input: this contains the initial guess to the solution of the
+* equation
+* Output: This contains the final approximation to the solution
+* determined by the multigrid solver.
+*
+* rhs :
+* This array contains the values of the right hand side of the
+* equation at every point on the rwo dimensional grid.
+*
+* eps :
+* eps > 0. The maximum norm of the residual is calculated at
+* the end of each multigrid cycle. The algorithm is
+* terminated when this maximum becomes less than tol
+* or when the maximum number of iterations (see ncyc)
+* is exceeded. It is up to the user to provide a
+* meaningfull tolerance criteria for the particular
+* problem being solved.
+* rmax:
+* This is the final value of the residual calcu;ated.
+*
+************************************************************************
+*
+*
+
+ integer idim,ilower,iupper,jdim,jlower,jupper,ifmg
+ integer id5,id9,idi,idg,ier
+ real*8 cc(idim,jdim),cn(idim,jdim),cs(idim,jdim),cw(idim,jdim),
+ & ce(idim,jdim),u(idim,jdim),rhs(idim,jdim)
+ real*8 eps,rmax
+
+*
+************************************************************************
+*
+* Variable definitions:
+*
+* Integers:
+* ifd59 :
+* ifd59 = 5 - means a 5-point finite difference stencil
+* (ac,an,as,aw,ae) is defined on the finest grid.
+* ifd59 = 9 - means a 9-point finite difference stencil
+* (ac,an,as,aw,ae,anw,ane,asw,ase) is defined on
+* the finest grid by the user.
+* (NOTE: This routine specifically written
+* for 5-point)
+*
+* ncyc :
+* The maximum number of multigrid "v"-cycles to be used. If
+* the maximum norm of the residual is not less than tol at the
+* end of ncyc cycles, the algorithm is terminated.
+* (NOTE: ncyc <= 40 )
+*
+* np(20,8) :
+* Input: When the iskip=1,-1 or -2 option is used, np2 is
+* assumed to contain the grid information for umgs2.
+* Output: When the iskip=0 option is used, the grid
+* information for umgs2 is returned in np2.
+* (NOTE: This is only useful for multiple instance problems)
+*
+* iskip :
+* iskip = 0 - The coarse grid information, coarse grid
+* operators and interpolation coefficients are
+* calculated by umgs2. This information is stored
+* in the arrays ac, aw, as, asw, ase, pu, pd, np2
+* and the variable.
+* iskip = 1 - The calculation of the coarse grid information,
+* coarse grid operators and interpolation
+* coefficients is skipped. This option would be
+* used when umgs2 has been called with iskip=0 and
+* is being called again to solve a system of
+* equations with the same matrix. This would be
+* the case in, say, parabolic problems with time
+* independent coefficients.
+* iskip =-1 - The set up of pointers (ugrdfn) is skipped.
+* Coarse grid operators and interpolation
+* coefficients are calculated and the given matrix
+* equation is solved. This option would be used
+* when umgs2 has been called with iskip=0 and is
+* being called again to solve a system of equations
+* with a different matrix of the same dimensions.
+* This would be the case for, say, parabolic
+* problems with time dependent coefficients.
+* iskip =-2 - The set up of pointers (ugrdfn) is skipped.
+* Coarse grid operators and interpolation
+* coefficients are calculated and returned.
+* No matrix solve.
+*
+* ipc :
+* ipc = 0 or 1.
+* ipc is a multigrid parameter which determines the type of
+* interpolation to be used. Usually ipc=1 is best. However,
+* if the boundary condition equations have been absorbed into
+* the interior equations then ipc=0 can be used which results
+* in a slightly more efficient algorithm.
+*
+* nman :
+* nman = 0 usually.
+* nman =1 signals that the fine grid equations are singular for
+* the case when homogeneous Neumann boundary conditions are
+* applied along the entire boundary. In this case, the
+* difference equations are singular and the condition that the
+* integral of q over the domain be zero is added to the set of
+* difference equations. This condition is satisfied by adding
+* the appropriate constant vector to q on the fine grid. It is
+* assumed, in this case, that a well-defined problem has been
+* given to mgss2, i.e. the integral of f over the domain is
+* zero.
+*
+* im :
+* The number of grid points in the x-direction (including two
+* ficticious points)
+* jm :
+* The number of grid points in the y-direction (including two
+* ficticious points)
+*
+* linp :
+* This is a dummy argument left over from the authors
+* development of the code
+* Use: common /io/ linp,lout
+*
+* lout :
+* lout = unit number of output file into which the maximum norm
+* of the residual after each multigrid v-cycle" is printed.
+* Use: common /io/ linp,lout
+*
+* iscale :
+* Flag to indicate whether problem can be diagonally scaled to
+* speed convergence of the multigrid solver.
+*
+* REALS:
+*
+* ac(id5),an(id5),as(id5),aw(id5),ae(id5),
+* anw(id9),ane(id9),asw(id9),ase(id9) :
+* Input: ac, an, as, aw, ae, anw, ane, asw and ase contain the
+* stencil coefficients for the difference operator on
+* the finest grid. When the iskip=1 option is used,
+* these arrays also are assumed to contain the coarse
+* grid difference stencil coeficients.
+* Output: when the iskip=0 option is used, the coarse grid
+* stencil coeficients are returned in ac, an, as, aw,
+* ae, anw, ane, asw and ase.
+*
+* ru(idi),rd(idi),rc(idi) :
+* Real work arrays.
+*
+* pu(idi),pd(idi),pc(idi) :
+* Real work arrays.
+* Input: when the iskip=1 option is used, these arrays are
+* assumed to contain the interpolation coefficients used
+* in the semi-coarsening multigrid algorithm.
+* Output: when the iskip=0 option is used, the interpolation
+* coeficients are returned in pu and pd.
+*
+* f(id5) :
+* f contains the right hand side vector of the matrix
+* equation to be solved by umgs2.
+*
+* q(id5) :
+* If ifmg=0, q contains the initial guess on the fine grid.
+* If ifmg=1, the initial guess on the fine grid is determined
+* by the full multigrid process and the value of
+* q on input to umgs2 not used.
+*
+* tol :
+* tol > 0. The maximum norm of the residual is calculated at
+* the end of each multigrid cycle. The algorithm is
+* terminated when this maximum becomes less than tol
+* or when the maximum number of iterations (see ncyc)
+* is exceeded. It is up to the user to provide a
+* meaningfull tolerance criteria for the particular
+* problem being solved.
+* tol = 0. Perform ncyc multigrid cycles. Calculate and print
+* the maximum norm of the residual after each cycle.
+* tol =-1. Perform ncyc multigrid cycles. The maximum norm of
+* the final residual is calculated and returned in
+* the variable rmax in the calling list of umgs2.
+* tol =-2. Perform ncyc multigrid cycles. The maximum norm of
+* the residual is not calculated.
+*
+************************************************************************
+*
+*
+
+ integer ncyc,ifd59
+ parameter (ncyc=40,ifd59=5)
+
+* integer id5,id9,idi,idg
+* This is for a 103x28 grid
+* parameter(id5=6695,id9=3811,idi=3811,idg=103)
+* This is for a 203x56 grid
+* parameter(id5=24969,id9=13601,idi=13601,idg=203)
+* This is for a 403x118 grid
+* parameter(id5=100750,id9=52793,idi=52793,idg=403)
+
+ integer np2(20,8)
+ integer iskip,ipc,nman
+ parameter (iskip=0,ipc=1,nman=0)
+ integer irc,irurd,im,jm
+ integer linp,lout
+ common /io/ linp,lout
+
+ real*8 ac(id5),an(id5),as(id5),aw(id5),ae(id5),
+ & anw(id9),ane(id9),asw(id9),ase(id9),
+ & q(id5),f(id5),gam(idg)
+
+ real*8 ru(idi),rd(idi),rc(idi),pu(idi),pd(idi),pc(idi)
+
+ real*8 tol
+
+ integer iscale,itry
+
+
+
+ ier=0
+*
+* Set some parameters for multigrid solver
+*
+ lout=6
+c rewind(unit=lout)
+ irc=0
+ irurd=0
+ im=iupper-ilower+3
+ jm=jupper-jlower+3
+ tol=eps
+
+
+*
+* Set up coefficients into vectors with correct indexing
+*
+ do 110 j=jlower,jupper
+ do 100 i=ilower,iupper
+ n=(j+1-jlower)*im + i+1-(ilower-1)
+ ac(n)=cc(i,j)
+ an(n)=cn(i,j)
+ as(n)=cs(i,j)
+ aw(n)=cw(i,j)
+ ae(n)=ce(i,j)
+ anw(n)=0.
+ ane(n)=0.
+ asw(n)=0.
+ ase(n)=0.
+ q(n)=u(i,j)
+ f(n)=rhs(i,j)
+ 100 continue
+ 110 continue
+
+
+*
+* Determine whether we can diagonal scale the problem to speed
+* convergence. Can only be done if there are no zeros on the main
+* diagonal (ie. central difference coefficient).
+*
+ iscale=1
+ do 200 j=jlower,jupper
+ do 205 i=ilower,iupper
+ n=(j+1-jlower)*im + i+1-(ilower-1)
+ if (ac(n) .eq. 0.) then
+ iscale=0
+ endif
+ 205 continue
+ 200 continue
+
+*
+* Do the diagonal scaling if we can.
+*
+ if (iscale.eq.1) then
+
+ do 210 j=jlower,jupper
+ do 215 i=ilower,iupper
+ n=(j+1-jlower)*im + i+1-(ilower-1)
+ f(n)=f(n)/ac(n)
+ ae(n)=ae(n)/ac(n)
+ aw(n)=aw(n)/ac(n)
+ as(n)=as(n)/ac(n)
+ an(n)=an(n)/ac(n)
+ ac(n)=1.
+ 215 continue
+ 210 continue
+
+ endif
+
+c
+c Now call the multigrid routine
+
+ itry=0
+
+ 1122 call umgs2(
+ + ac,aw,as,ae,an,asw,ase,ane,anw,q,f,pu,pd,pc,ru,rd,rc,gam,np2,
+ + ifd59,ifmg,ncyc,tol,nman,im,jm,id5,id9,idi,m,iskip,rmax,
+ + ipc,irc,irurd)
+
+
+ if ((rmax.gt.tol).and.(itry.le.5)) then
+ itry=itry+1
+ print*,"Retry #",itry
+ goto 1122
+ endif
+
+ if (rmax.gt.tol) then
+ ier = -1
+ print*,"Did not converge"
+ print*," maximum residual = ",rmax
+ print*," tolerance = ",tol
+ endif
+
+*
+* Convert the solution back to the 2D array form
+*
+ do 510 j=jlower,jupper
+ do 500 i=ilower,iupper
+ n=(j+1-jlower)*im + i+1-(ilower-1)
+ u(i,j)=q(n)
+ 500 continue
+ 510 continue
+
+ return
+ end
diff --git a/src/psi_1st_deriv.x b/src/psi_1st_deriv.x
new file mode 100644
index 0000000..3c3335e
--- /dev/null
+++ b/src/psi_1st_deriv.x
@@ -0,0 +1,18 @@
+ o1 = 5.0000000000000d-1*eta(i,j,k)
+ o2 = exp(o1)
+ o3 = psi2dv(i,j,k)
+ o4 = 1/o3
+ o5 = cos(phi(i,j,k))
+ o6 = cos(q(i,j,k))
+ o7 = dqpsi2dv(i,j,k)
+ o8 = -1.50000000000000d0*eta(i,j,k)
+ o9 = exp(o8)
+ o10 = detapsi2dv(i,j,k)
+ o11 = sin(q(i,j,k))
+ o12 = sin(phi(i,j,k))
+ psix(i,j,k) = o2*o4*(o10*o11*o5*o9 - 5.0000000000000d-1*o11*o3*o
+ & 5*o9 + o5*o6*o7*o9)
+ psiy(i,j,k) = o2*o4*(o10*o11*o12*o9 - 5.0000000000000d-1*o11*o12
+ & *o3*o9 + o12*o6*o7*o9)
+ psiz(i,j,k) = o2*o4*(o10*o6*o9 - 5.0000000000000d-1*o3*o6*o9 - o
+ & 11*o7*o9)
diff --git a/src/psi_2nd_deriv.x b/src/psi_2nd_deriv.x
new file mode 100644
index 0000000..0cfc82a
--- /dev/null
+++ b/src/psi_2nd_deriv.x
@@ -0,0 +1,51 @@
+ o1 = 5.0000000000000d-1*eta(i,j,k)
+ o2 = exp(o1)
+ o3 = psi2dv(i,j,k)
+ o4 = 1/o3
+ o5 = cos(phi(i,j,k))
+ o6 = o5**2
+ o7 = cos(q(i,j,k))
+ o8 = o7**2
+ o9 = detapsi2dv(i,j,k)
+ o10 = -2.50000000000000d0*eta(i,j,k)
+ o11 = exp(o10)
+ o12 = dqqpsi2dv(i,j,k)
+ o13 = sin(phi(i,j,k))
+ o14 = o13**2
+ o15 = detaqpsi2dv(i,j,k)
+ o16 = sin(q(i,j,k))
+ o17 = dqpsi2dv(i,j,k)
+ o18 = detaetapsi2dv(i,j,k)
+ o19 = o16**2
+ o20 = tan(q(i,j,k))
+ o21 = 1/o20
+ psixx(i,j,k) = o2*o4*(1.00000000000000d0*o11*o14*o17*o21 - 5.000
+ & 0000000000d-1*o11*o14*o3 + o11*o18*o19*o6 + 7.5000000000000d-1*o
+ & 11*o19*o3*o6 + 2.00000000000000d0*o11*o15*o16*o6*o7 - 3.00000000
+ & 000000d0*o11*o16*o17*o6*o7 + o11*o12*o6*o8 - 5.0000000000000d-1*
+ & o11*o3*o6*o8 + o11*o14*o9 - 2.00000000000000d0*o11*o19*o6*o9 + o
+ & 11*o6*o8*o9)
+ psixy(i,j,k) = o2*o4*(o11*o13*o18*o19*o5 - o11*o13*o17*o21*o5 +
+ & 5.0000000000000d-1*o11*o13*o3*o5 + 7.5000000000000d-1*o11*o13*o1
+ & 9*o3*o5 + 2.00000000000000d0*o11*o13*o15*o16*o5*o7 - 3.000000000
+ & 00000d0*o11*o13*o16*o17*o5*o7 + o11*o12*o13*o5*o8 - 5.0000000000
+ & 000d-1*o11*o13*o3*o5*o8 - o11*o13*o5*o9 - 2.00000000000000d0*o11
+ & *o13*o19*o5*o9 + o11*o13*o5*o8*o9)
+ psixz(i,j,k) = o2*o4*(-(o11*o15*o19*o5) + 1.50000000000000d0*o11
+ & *o17*o19*o5 - o11*o12*o16*o5*o7 + o11*o16*o18*o5*o7 + 1.25000000
+ & 000000d0*o11*o16*o3*o5*o7 + o11*o15*o5*o8 - 1.50000000000000d0*o
+ & 11*o17*o5*o8 - 3.00000000000000d0*o11*o16*o5*o7*o9)
+ psiyy(i,j,k) = o2*o4*(o11*o14*o18*o19 + 7.5000000000000d-1*o11*o
+ & 14*o19*o3 + 1.00000000000000d0*o11*o17*o21*o6 - 5.0000000000000d
+ & -1*o11*o3*o6 + 2.00000000000000d0*o11*o14*o15*o16*o7 - 3.0000000
+ & 0000000d0*o11*o14*o16*o17*o7 + o11*o12*o14*o8 - 5.0000000000000d
+ & -1*o11*o14*o3*o8 - 2.00000000000000d0*o11*o14*o19*o9 + o11*o6*o9
+ & + o11*o14*o8*o9)
+ psiyz(i,j,k) = o2*o4*(-(o11*o13*o15*o19) + 1.50000000000000d0*o1
+ & 1*o13*o17*o19 - o11*o12*o13*o16*o7 + o11*o13*o16*o18*o7 + 1.2500
+ & 0000000000d0*o11*o13*o16*o3*o7 + o11*o13*o15*o8 - 1.500000000000
+ & 00d0*o11*o13*o17*o8 - 3.00000000000000d0*o11*o13*o16*o7*o9)
+ psizz(i,j,k) = o2*o4*(o11*o12*o19 - 5.0000000000000d-1*o11*o19*o
+ & 3 - 2.00000000000000d0*o11*o15*o16*o7 + 3.00000000000000d0*o11*o
+ & 16*o17*o7 + o11*o18*o8 + 7.5000000000000d-1*o11*o3*o8 + o11*o19*
+ & o9 - 2.00000000000000d0*o11*o8*o9)
diff --git a/src/shmgp.F77 b/src/shmgp.F77
new file mode 100644
index 0000000..1982591
--- /dev/null
+++ b/src/shmgp.F77
@@ -0,0 +1,1455 @@
+
+c----------------------------------------------------------------------
+ subroutine umgs2(
+ + ac,aw,as,ae,an,asw,ase,ane,anw,q,f,pu,pd,pc,ru,rd,rc,gam,np2,
+ + ifd59,ifmg,ncyc,tol,nman,im,jm,id5,id9,idi,m,iskip,rmax,
+ + ipc,irc,irurd)
+ implicit real*8(a-h,o-z)
+c----------------------------------------------------------------------
+cdir$ noinline
+c** SUBROUTINE UMGS2
+c
+c** COPYRIGHT: Ecodynamics Research Associates, Inc.
+c
+c** Date written: June, 1990
+c** Author: Steve Schaffer
+c Mathematics Department
+c New Mexico Tech
+c Socorro, NM 87801
+c 505-835-5811
+c
+c** DESCRIPTION:
+c umgs2 is a black box symmetric matrix solver. It is written
+c in unsymmetric storage mode and can be used to solve mildly
+c nonsymmetric problems. The user provides a matrix and right hand
+c side vector corresponding to a 5 or 9 point finite difference/
+c finite volume discretization of a symmetric second order PDE.
+c umgs2 will construct a sequence of coarse grids and coarse
+c grid operators and then solve the matrix equation using a
+c y-direction semi-coarsening multigrid algorithm and return
+c the solution vector. If a sequence of matrix problems are
+c to be solved using the same matrix, computational time can
+c be saved by skipping the construction of the coarse grid
+c information in subsequent calls to umgs2.
+c
+c The matrix on the finest grid is stored in the arrays ac,aw,as,
+c an,asw,ase,ane and anw. The difference stencil at the point
+c (i,j) given by
+c
+c nw n ne anw(i,j) an(i,j) ane(i,j)
+c w c e = aw(i,j) ac(i,j) aw(i,j)
+c sw s se asw(i,j) as(i,j) ase(i,j)
+c
+c If the difference stencil on the fine grid is a 5 point stencil
+c then the arrays asw,ase,ane,anw are not used and the
+c stencil is given by
+c
+c n an(i,j)
+c w c e = aw(i,j) ac(i,j) ae(i,j)
+c s as(i,j)
+c
+c However, asw,ase,ane,anw still need to be dimensioned (by id9)
+c in the calling program as they are used in the coarse grid
+c calculations.
+c
+c** STORAGE
+c It is assumed that a set of ficticious points have been defined
+c along the entire boundary. These points have nothing to do with
+c the solution and are used for programming convenience and
+c vectorization purposes. Storage is allocated for the stencil
+c elements ac,aw,as,asw,ase, the solution vector, q, and the
+c right hand side vector, f, at these ficticious points. The
+c stencils at these ficticious points and all stencil connections
+c to them are set to zero in the subroutine useta which is called
+c by umgs2. The computational grid is depicted by
+c
+c x x x x x x x x
+c
+c x * * * * * * x
+c
+c x * * * * * * x
+c .
+c .
+c .
+c x * * * * * * x
+c
+c x * * * * * * x
+c
+c x x x x x x x x
+c
+c where x depicts the ficticious points and * depicts the interior
+c points. The total storage requirements for the fine grid problem
+c is then 5*im*jm for 5 point stencils and 7*im*jm for 9 point
+c stencils. The total storage requirements for the multigrid
+c solution is approximately 2 to 3 times that of the storage
+c requirements of the fine grid problem. (See DIMENSION PARAMETERS).
+c Note: The first im*jm elements of the arrays ac,aw,as,[asw,ase],
+c q and f correspond to the finest grid.
+c
+c** DIMENSION PARAMETERS
+c The arrays ac,aw,ae,asw,ase,q,f,pu,pd and pc are dimensioned as one
+c dimensional arrays in the calling program. They are dimensioned
+c as two dimensional arrays in the working subroutines. The one
+c dimensional storage of the arrays, say q, follows: n=(j-1)*jm+i,
+c where n is the element location in the one dimensional storage of
+c q corresponding to the (i,j)th element of the two dimensional
+c storage of q and jm is the number of grid points in the j
+c direction (including the two ficticious points).
+c
+c The dimension parameters are id5, id9, idi and idg. They can be
+c determined by running the companion program MSS2DIM.F.
+c id5 - Integer variable.
+c Dimension of the arrays ac,aw,as,ae,an,q and f in the
+c calling program. id5 is the total number of grid points
+c on the finest grid and all coarser grids.
+c id9 - Integer variable.
+c Dimension of the arrays asw,ase,ane,anw in the calling
+c program. If ifd59=5 then id9=idi. If ifd59=9 then
+c id9=id5.
+c idi - Integer variable.
+c Dimension of the work arrays pu and pd in the calling
+c program. idi is the total number of grid points on all
+c of the coarser grids.
+c idg - Integer variable.
+c Dimension of the work array gam in the calling program.
+c It is set to the value im, the number of grid points
+c in the i-direction on the finest grid.
+c
+c** INPUT
+c (Note: all variable types are set implicitly)
+c ac,aw,as
+c ae,an - Real arrays. Dimensioned (id5) in calling program.
+c See comments in DESCRIPTION and DIMENSION PARAMETERS.
+c asw,ase
+c ane,anw - Real arrays. Dimensioned (id9) in calling program.
+c See comments in DESCRIPTION and DIMENSION PARAMETERS.
+c f - Real array. Dimensioned (id5) in calling program.
+c f contains the right hand side vector of the matrix
+c equation to be solved by umgs2.
+c q - Real array. Dimensioned (id5) in calling program.
+c If ifmg=0, q contains the initial guess on the fine
+c grid. If ifmg=1, the initial guess on the fine grid
+c is determined by the full multigrid process and the
+c value of q on input to umgs2 not used.
+c ifd59 - Integer variable.
+c =5 - means a 5-point finite difference stencil (ac,aw and
+c as) is defined on the finest grid by the user.
+c =9 - means a 9-point finite difference stencil (ac,aw,as,
+c asw, ase) is defined on the finest grid by the user.
+c ifmg - Integer variable.
+c =0 - The full multigrid algorithm is not used to obtain a
+c good initial guess on the fine grid.
+c =1 - The full multigrid algorithm is used to obtain a good
+c initial guess on the fine grid.
+c ncyc - Integer variable.
+c The maximum number of multigrid v-cycles to be used.
+c If the maximum norm of the residual is not less than tol
+c at the end of ncyc cycles, the algorithm is terminated.
+c tol - Real variable.
+c >0 - The maximum norm of the residual is calculated at the
+c end of each multigrid cycle. The algorithm is terminated
+c when this maximum becomes less than tol or when the maximum
+c number of iterations (see ncyc) is exceeded. It is up to
+c the user to provide a meaningfull tolerance criteria for
+c the particular problem being solved.
+c =0 - Perform ncyc multigrid cycles. Calculate and print
+c the maximum norm of the residual after each cycle.
+c =-1. - Perform ncyc multigrid cycles. The maximum norm of
+c the final residual is calculated and returned in the
+c variable rmax in the calling list of umgs2.
+c =-2. - Perform ncyc multigrid cycles. The maximum norm of
+c the residual is not calculated.
+c iskip - Integer variable.
+c =0 - The coarse grid information, coarse grid operators
+c and interpolation coefficients are calculated by
+c umgs2. This information is stored in the arrays
+c ac, aw, as, asw, ase, pu, pd, np2 and the variable m
+c and returned to the calling program.
+c =1 - The calculation of the coarse grid information, coarse
+c grid operators and interpolation coefficients is
+c skipped. This option would be used when umgs2 has
+c been called with iskip=0 and is being called again
+c to solve a system of equations with the same matrix.
+c This would be the case in, say, parabolic problems
+c with time independent coefficients.
+c =-1 -The set up of pointers (ugrdfn) is skipped. Coarse grid
+c operators and interpolation coefficients are calculated
+c and the given matrix equation is solved. This option
+c would be used when umgs2 has been called with iskip=0
+c and is being called again to solve a system of
+c equations with a different matrix of the same
+c dimensions. This would be the case for, say,
+c parabolic problems with time dependent coefficients.
+c =-2 -The set up of pointers (ugrdfn) is skipped. Coarse grid
+c operators and interpolation coefficients are calculated
+c and returned to the calling program. No matrix solve.
+c ipc - Integer variable.
+c =0 or 1.
+c ipc is a multigrid parameter which determines the type of
+c interpolation to be used. Usually ipc=1 is best. However, if
+c the boundary contition equations have been absorbed into the
+c interior equations then ipc=0 can be used which results in a
+c slightly more efficient algorithm.
+c nman - Integer variable.
+c =0 usually.
+c =1 signals that the fine grid equations are singular for
+c the case when homogeneous Neumann boundary conditions are
+c applied along the entire boundary. In this case, the
+c difference equations are singular and the condition that
+c the integral of q over the domain be zero is added to the
+c set of difference equations. This condition is satisfied
+c by adding the appropriate constant vector to q on the fine
+c grid. It is assumed, in this case, that a well-defined
+c problem has been given to mgss2, i.e. the integral of f
+c over the domain is zero.
+c im - Integer variable.
+c The number of grid points in the x-direction (including two
+c ficticious points)
+c jm - Integer variable.
+c The number of grid points in the y-direction (including two
+c ficticious points)
+c lout - Integer variable.
+c = unit number of output file into which the maximum norm
+c of the residual after each multigrid v-cycle is printed.
+c Use: common /iout/ lout
+c
+c** INPUT/OUTPUT
+c q - Real array. Dimensioned (id5)
+c On input, if ifmg=0, q contains the initial guess on the
+c finest grid for umgs2. On output, q contains the final
+c solution on the finest grid.
+c ac-anw - Real arrays. See DIMENSION.
+c On input, ac, aw, as, [asw and ase] contain the stencil
+c coefficients for the difference operator on the finest
+c grid. When the iskip=1 option is used, these arrays
+c also are assumed to contain the coarse grid difference
+c stencil coeficients.
+c On output, when the iskip=0 option is used, the coarse
+c grid stencil coeficients are returned in ac - ase.
+c
+c ru,rd,rc - Real work arrays. Dimensioned (idi)
+c
+c pu,pd,pc - Real work arrays. Dimensioned (idi).
+c On input, when the iskip=1 option is used, these arrays
+c are assumed to contain the interpolation coefficients
+c used in the semi-coarsening multigrid algorithm.
+c On output, when the iskip=0 option is used, the
+c interpolation coeficients are returned in pu and pd.
+c np2 - Integer work array. Dimensioned np2(20,8).
+c On input, when the iskip=1,-1 or -2 option is used, np2 is
+c assumed to contain the grid information for umgs2.
+c On output, when the iskip=0 option is used, the grid
+c information for umgs2 is returned in np2.
+c** OUTPUT
+c rmax - If tol.ge.-1., the final residual norm is returned in rmax.
+c
+c** SUBROUTINES CALLED BY UMGS2
+c
+c - ugrdfn, ukey, uintad, urelax, urscal, ursrhs, useta
+c
+c** END OF DESCRIPTION OF UMGS2
+c .....................................................................
+ real*8 ac(id5),aw(id5),as(id5),ae(id5),an(id5),asw(id9),
+ + ase(id9),ane(id9),anw(id9),q(id5),f(id5)
+ real*8 pu(idi),pd(idi),pc(idi),gam(im)
+ integer np2(20,8)
+ real*8 ru(idi),rd(idi),rc(idi)
+ real*8 resid(0:40),confac(0:40)
+ common /io/ linp,lout
+c
+c-time tsu0=second()
+ if(iskip.eq.0) then
+ call ugrdfn(m,ifd59,is5,is9,isi,np2,im,jm)
+ iquit=0
+ if(m.gt.20) then
+ iquit=1
+ write(lout,*) ' m=',m,' > 20 - np2 is dimensioned np2(m=20,8)'
+ endif
+ if(is5.gt.id5) then
+ iquit=1
+ write(lout,*) ' id5=',id5,' too small. Should be set to',is5
+ endif
+ if(is9.gt.id9) then
+ iquit=1
+ write(lout,*) ' id9=',id9,' too small. Should be set to',is9
+ endif
+ if(isi.gt.idi) then
+ iquit=1
+ write(lout,*) ' idi=',idi,' too small. Should be set to',isi
+ endif
+ if(is5.lt.2*im*jm) then
+ iquit=1
+ write(lout,*) ' id5.lt.2*im*jm can cause problems in useta'
+ write(lout,*) ' this can be remedied by setting id5 larger'
+ endif
+ if(iquit.eq.1) return
+ endif
+ if(iskip.le.0) then
+c ---------- interpolation and coarse grid operators -----------
+ do 5 k=m-1,1,-1
+cdir$ inline
+ call ukey(k+1,np2,n5,n9,ni,jm,i9,j9,ifd,jr)
+ call ukey(k,np2,n5c,n9c,nic,jmc,i9c,j9c,ifdc,jrc)
+cdir$ noinline
+ if(k.eq.m-1) n5cqf=n5c
+ 5 call useta(
+ + ac(n5),aw(n5),as(n5),ae(n5),an(n5),asw(n9),ase(n9),
+ + ane(n9),anw(n9),ac(n5c),aw(n5c),as(n5c),ae(n5c),an(n5c),asw(n9c),
+ + ase(n9c),ane(n9c),anw(n9c),pu(nic),pd(nic),pc(nic),ru(nic),
+ + rd(nic),rc(nic),q(n5cqf),f(n5cqf),gam,
+ + im,jm,jmc,ifd,i9,j9,nman,k+1,m,jr,ipc,irc,irurd)
+ endif
+ if(iskip.eq.-2) return
+c
+ if(ifmg.ge.1) then
+ do 6 k=m-1,1,-1
+cdir$ inline
+ call ukey(k+1,np2,n5,n9,ni,jm,i9,j9,ifd,jr)
+ call ukey(k,np2,n5c,n9c,nic,jmc,i9c,j9c,ifdc,jrc)
+cdir$ noinline
+ 6 call ursrhs(f(n5),f(n5c),pu(nic),pd(nic),pc(nic),ru(nic),
+ + rd(nic),rc(nic),im,jm,jmc,m,k+1,jr,irc)
+ endif
+c-time tsu1=second()
+c-time write(lout,*) ' time for setup =',tsu1-tsu0
+ l=1
+ if(ifmg.eq.0) l=m
+ k=l
+ mcyc=0
+ rmaxo=1.
+c ---------- begin multigrid cycling ----------------------------
+c
+ if(l.eq.1) go to 20
+cdir$ inline
+ 10 call ukey(k,np2,n5,n9,ni,jm,i9,j9,ifd,jr)
+cdir$ noinline
+ call urelax(
+ + ac(n5),aw(n5),as(n5),ae(n5),an(n5),asw(n9),ase(n9),
+ + ane(n9),anw(n9),f(n5),q(n5),gam,
+ + im,jm,i9,j9,ifd,nman,k,m,jr,0,0)
+cdir$ inline
+ call ukey(k-1,np2,n5c,n9c,nic,jmc,i9c,j9c,ifdc,jrc)
+cdir$ noinline
+ call urscal(
+ + ac(n5),aw(n5),as(n5),ae(n5),an(n5),asw(n9),ase(n9),
+ + ane(n9),anw(n9),q(n5),f(n5),f(n5c),q(n5c),rc(nic),
+ + im,jm,jmc,ifd,i9,j9,k,m,jr,tol,rmax,ipc,irc)
+ if(k.eq.m.and.rmax.lt.tol) go to 60
+ if(k.eq.m.and.tol.ge.-.5) then
+ if(rmaxo.ne.0.) rate=rmax/rmaxo
+ rmaxo=rmax
+ if(mcyc.eq.0) rmax0=rmax
+ resid(mcyc)=rmax
+ confac(mcyc)=rate
+ endif
+ if(tol.eq.-.5) write(lout,*) ' down ',k,rmax
+ k=k-1
+ if(k.gt.1) go to 10
+c --------- solve coarsest grid ----------------------------------
+c
+cdir$ inline
+ 20 call ukey(1,np2,n5,n9,ni,jm,i9,j9,ifd,jr)
+cdir$ noinline
+ call urelax(
+ + ac(n5),aw(n5),as(n5),ae(n5),an(n5),asw(n9),ase(n9),
+ + ane(n9),anw(n9),f(n5),q(n5),gam,
+ + im,jm,i9,j9,ifd,nman,k,m,jr,0,0)
+ if(l.eq.1) go to 40
+c ---------- interpolate correction to next finer grid -----------
+c
+ 30 k=k+1
+cdir$ inline
+ call ukey(k,np2,n5,n9,ni,jm,i9,j9,ifd,jr)
+ call ukey(k-1,np2,n5c,n9c,nic,jmc,i9c,j9c,ifdc,jrc)
+cdir$ noinline
+ call uintad(
+ + q(n5),q(n5c),pu(nic),pd(nic),im,jm,jmc,1,jr,ipc)
+ call urelax(
+ + ac(n5),aw(n5),as(n5),ae(n5),an(n5),asw(n9),ase(n9),
+ + ane(n9),anw(n9),f(n5),q(n5),gam,
+ + im,jm,i9,j9,ifd,nman,k,m,jr,0,0)
+ if(tol.eq.-.5) then
+ call urscal(
+ + ac(n5),aw(n5),as(n5),ae(n5),an(n5),asw(n9),ase(n9),
+ + ane(n9),anw(n9),q(n5),f(n5),f(n5c),q(n5c),rc(nic),
+ + im,jm,jmc,ifd,i9,j9,k,m,jr,tol,rmax,ipc,irc)
+ write(lout,*) ' up ',k,rmax
+ endif
+ if(k.lt.l) go to 30
+ if(l.eq.m) go to 50
+c ---------- interpolate solution to new finest grid l+1 in fmg ----
+c
+ 40 l=l+1
+ k=l
+cdir$ inline
+ call ukey(l,np2,n5,n9,ni,jm,i9,j9,ifd,jr)
+ call ukey(l-1,np2,n5c,n9c,nic,jmc,i9c,j9c,ifdc,jrc)
+cdir$ noinline
+ call uintad(
+ + q(n5),q(n5c),pu(nic),pd(nic),im,jm,jmc,0,jr,0)
+ go to 10
+c
+ 50 if(nman.eq.1) call uneuman(q(n5),im,jm)
+ mcyc=mcyc+1
+c ---------- Cycle ncyc times on grid m ----------------------------
+ if(mcyc.lt.ncyc) go to 10
+c-time tmg1=second()
+c-time write(lout,*) ' time in ',ncyc,' cycles =',tmg1-tsu1
+c
+c ---------- print out final residual and work units ---------------
+ if(tol.ge.-1.) then
+cdir$ inline
+ call ukey(m,np2,n5,n9,ni,jm,i9,j9,ifd,jr)
+ call ukey(m-1,np2,n5c,n9c,nic,jmc,i9c,j9c,ifdc,jrc)
+cdir$ noinline
+ call urscal(
+ + ac(n5),aw(n5),as(n5),ae(n5),an(n5),asw(n9),ase(n9),
+ + ane(n9),anw(n9),q(n5),f(n5),f(n5c),q(n5c),rc(nic),
+ + im,jm,jmc,ifd,i9,j9,k,m,jr,1.,rmax,ipc,irc)
+ resid(mcyc)=rmax
+ confac(mcyc)=rmax/rmaxo
+ nb=0
+ ne=min0(6,mcyc)
+ 2029 write(lout,2033) (mc,mc=nb,ne)
+ write(lout,2032) (resid(mc),mc=nb,ne)
+ write(lout,2031) (confac(mc),mc=nb,ne)
+ nb=ne+1
+ ne=ne+min0(6,mcyc-ne)
+ if(nb.le.ne) go to 2029
+ fconfac=(rmax/rmax0)**(1./float(mcyc))
+ write(lout,2034) fconfac
+ 2034 format(30x,6(1h*)/,' average convergence factor =',f7.3,/,
+ + 30x,6(1h*))
+ 2033 format(7(4x,i2,4x))
+ 2031 format(7(1x,f9.3))
+ 2032 format(7(1x,e9.3))
+ endif
+ return
+ 60 write(lout,1003) mcyc,resid(mcyc),tol
+ return
+ 1003 format(' cyc=',i2,' max(res)=',1pe8.2/
+ + ' tolerance condition tol=',1pe8.2,' satisfied')
+ end
+c----------------------------------------------------------------------
+ subroutine ugrdfn(m,ifd59,is5,is9,isi,np2,imx,jmx)
+ implicit real*8(a-h,o-z)
+c----------------------------------------------------------------------
+cdir$ noinline
+c Given imx, jmx and ifd59 (See comments in mgss2), ugrdfn calculates
+c the number of grids that will be needed. Pointers into the arrays
+c ac, aw, as, asw, ase, q, f, pu, pd, pc, ru, rd and rc and the size
+c of each grid is calculated and stored in the array np2. The
+c subroutine ukey is called to retrieve the grid information.
+c .....................................................................
+ parameter(n5=1,n9=2,ni=3,jm=4,i9=5,j9=6,ifd=7,jred=8)
+ integer np2(20,8)
+ common /cs/ icorstr,iprint
+ iq5=1
+ iq9=1
+ iqi=1
+ m=1
+ np2(m,1)=jmx
+ np2(m,2)=3
+ 10 if(np2(m,1).le.3) go to 20
+ m=m+1
+ np2(m,1)=np2(m-1,1)/2+1
+ if(np2(m-1,2).eq.2.and.mod(np2(m-1,1),2).eq.1)
+ + np2(m,1)=np2(m,1)+1
+ np2(m,2)=2
+ go to 10
+ 20 do 30 k=1,m
+ np2(m-k+1,jm)=np2(k,1)
+ 30 np2(m-k+1,jred)=np2(k,2)
+ do 40 k=m,1,-1
+ ktot=imx*np2(k,jm)
+ np2(k,n5)=iq5
+ iq5=iq5+ktot
+ np2(k,n9)=iq9
+ if(k.lt.m.or.ifd59.eq.9) iq9=iq9+ktot
+ np2(k,ni)=iqi
+ 40 if(k.lt.m) iqi=iqi+ktot
+ do 50 k=1,m
+ np2(k,i9)=imx
+ np2(k,j9)=np2(k,jm)
+ 50 np2(k,ifd)=9
+ if(ifd59.eq.5) then
+ np2(m,i9)=1
+ np2(m,j9)=1
+ np2(m,ifd)=5
+ endif
+ is5=iq5-1
+ is9=iq9-1
+ isi=iqi-1
+ return
+ end
+c----------------------------------------------------------------------
+ subroutine ukey(k,np2,nn5,nn9,nni,jjm,ii9,jj9,iifd,jjred)
+ implicit real*8(a-h,o-z)
+c----------------------------------------------------------------------
+c Returns the grid pointers and dimension variables for grid k. The
+c information is stored in the array np2.
+c......................................................................
+ parameter(n5=1,n9=2,ni=3,jm=4,i9=5,j9=6,ifd=7,jred=8)
+ integer np2(20,8)
+ nn5=np2(k,n5)
+ nn9=np2(k,n9)
+ nni=np2(k,ni)
+ jjm=np2(k,jm)
+ ii9=np2(k,i9)
+ jj9=np2(k,j9)
+ iifd=np2(k,ifd)
+ jjred=np2(k,jred)
+ return
+ end
+c----------------------------------------------------------------------
+ subroutine uintad(q,qc,pu,pd,im,jm,jmc,iadd,jred,ipc)
+ implicit real*8 (a-h,o-z)
+c----------------------------------------------------------------------
+c iadd=1:
+c Interpolates and adds the coarse grid (kf-1) correction, qc, to the
+c fine grid (kf) approximation, q, at the black y-lines.
+c iadd=0:
+c In the full multigrid algorithm, the solution to the coarse grid
+c (kf-1) difference equation is interpolated to the fine grid (kf)
+c to be used as the initial guess vector for kf=2,3,...,m.
+c Interpolation is at black y-lines only.
+c .....................................................................
+ real*8 q(im,jm),qc(im,jmc),pu(im,jmc),pd(im,jmc)
+ im1=im-1
+ jm1=jm-1
+ jblack=5-jred
+c add correction to next finer grid
+ 1000 if(iadd.eq.1) then
+ jc=3-jred
+ do 10 j=jblack,jm1,2
+ jc=jc+1
+ do 10 i=2,im1
+ 10 q(i,j)=q(i,j)+pd(i,jc)*qc(i,jc)+pu(i,jc)*qc(i,jc+1)
+c
+c interpolate solution to next finer grid in fmg
+ 1001 else
+ jc=3-jred
+ do 40 j=jblack,jm1,2
+ jc=jc+1
+ do 40 i=2,im1
+ 40 q(i,j)=pd(i,jc)*qc(i,jc)+pu(i,jc)*qc(i,jc+1)
+ 1002 endif
+ return
+ end
+c----------------------------------------------------------------------
+ subroutine uneuman(q,im,jm)
+ implicit real*8 (a-h,o-z)
+c----------------------------------------------------------------------
+c For problems with homogeneous Neumann boundary contitions, the
+c condition that the integral of q over the domain be zero is added
+c to the set of difference equations in order to obtain a unique
+c solution.
+c......................................................................
+ real*8 q(im,jm)
+ im1=im-1
+ jm1=jm-1
+ con=0.
+ do 10 j=2,jm1
+ do 10 i=2,im1
+ 10 con=con+q(i,j)
+ con=con/((im-2)*(jm-2))
+ do 20 j=2,jm1
+ do 20 i=2,im1
+ 20 q(i,j)=q(i,j)-con
+ return
+ end
+c----------------------------------------------------------------------
+ subroutine urelax(ac,aw,as,ae,an,asw,ase,ane,anw,f,q,gam,
+ + im,jm,i9,j9,ifd,nman,k,m,jred,ipc,iprcud)
+ implicit real*8 (a-h,o-z)
+c----------------------------------------------------------------------
+c Performs red/black x-line relaxation. The Thomas algorithm is used
+c to solve the tridiagonal matrices.
+c** INPUT -
+c ac-anw= finite difference operator coeficients
+c q= initial approximation
+c f= right hand side vector
+c im,jm= the number of grid points in the x,y-directions
+c i9,j9= the i,j-dimensions of the arrays asw,ase
+c ifd= 5 or 9 - the size of the stencil
+c nman- =0 usually.
+c =1 signals that the fine grid equations are singular
+c for the case when Neumann boundary conditions are
+c applied along the entire boundary. In this case, the
+c equations on the coarsest grid (consisting of a single
+c line of unknowns) is a singular tridiagonal system
+c and the Thomas algorithm is modified on this grid to
+c obtain a solution with an arbitrary constant vector
+c component. This constant vector is removed on the
+c finest grid by the call to subroutine uneuman.
+c** OUTPUT -
+c q= final approximation after a red/black relaxation sweep
+c .....................................................................
+ real*8 ac(im,jm),aw(im,jm),as(im,jm),ae(im,jm),an(im,jm),
+ + asw(i9,j9),ase(i9,j9),ane(i9,j9),anw(i9,j9)
+ real*8 f(im,jm),q(im,jm),gam(im)
+ jm1=jm-1
+ im1=im-1
+ im2=im-2
+ jblack=5-jred
+c usual red/black relaxatio
+ nrel=2
+ jrb=jred
+c ipc ..brbr relaxation swee
+ 1000 if(iprcud.eq.1) then
+ nrel=ipc
+ if(mod(ipc,2).eq.0) jrb=jblack
+c 1 black relax for calc g pu,pd,ru,
+ 1001 elseif(iprcud.eq.2) then
+ nrel=1
+ jrb=jblack
+ 1002 endif
+c
+c
+ do 109 nrr=1,nrel
+ 5000 if(jrb.eq.jblack) then
+c black rela
+ 6000 if(jblack.le.jm1) then
+ 1400 if(iprcud.ne.2) then
+c
+ do 110 j=jblack,jm1,2
+ do 110 i=2,im1
+ 110 q(i,j)=f(i,j)-as(i,j)*q(i,j-1)-an(i,j)*q(i,j+1)
+ 7000 if(ifd.eq.9) then
+ do 120 j=jblack,jm1,2
+ do 120 i=2,im1
+ 120 q(i,j)=q(i,j)-asw(i,j)*q(i-1,j-1)-ase(i,j)*q(i+1,j-1)-
+ + anw(i,j)*q(i-1,j+1)-ane(i,j)*q(i+1,j+1)
+ 7001 endif
+ 1401 endif
+c black tridiagonal solve
+c**
+c** Moved calculation of loop 129 from loop 130 for vectorization
+c** on vector machines (ie. Cray)
+c** By: John Towns 2/6/92
+c**
+ do 129 j=jblack,jm1,2
+ 129 q(2,j)=q(2,j)/ac(2,j)
+
+c**
+c** Changed bet=(quantity) to bet=1./(quantity) to trade two divisions
+c** for one division and two multiplies (more efficient on all
+c** machines)
+c** By: John Towns 2/6/92
+c**
+
+c**
+c** Cray compiler directives to parallelize tridiagonal solve.
+c** By: John Towns 4/13/92
+c**
+
+cmic$ parallel private(bet,gam,i)
+cmic$1shared(ac,ae,aw,q,jblack,jm1,im1,im2)
+cmic$ do parallel
+ do 130 j=jblack,jm1,2
+ bet=1./ac(2,j)
+ do 140 i=3,im1
+ gam(i)=ae(i-1,j)*bet
+ bet=1./(ac(i,j)-aw(i,j)*gam(i))
+ 140 q(i,j)=(q(i,j)-aw(i,j)*q(i-1,j))*bet
+ do 150 i=im2,2,-1
+ 150 q(i,j)=q(i,j)-gam(i+1)*q(i+1,j)
+ 130 continue
+cmic$ end do
+cmic$ end parallel
+ 6001 endif
+c red relax
+ 5001 else
+c
+ do 210 j=jred,jm1,2
+ do 210 i=2,im1
+ 210 q(i,j)=f(i,j)-as(i,j)*q(i,j-1)-an(i,j)*q(i,j+1)
+ 1100 if(ifd.eq.9) then
+ do 220 j=jred,jm1,2
+ do 220 i=2,im1
+ 220 q(i,j)=q(i,j)-asw(i,j)*q(i-1,j-1)-ase(i,j)*q(i+1,j-1)-
+ + anw(i,j)*q(i-1,j+1)-ane(i,j)*q(i+1,j+1)
+ 1101 endif
+c tridiagonal solve
+c nman=1 ==> avoid singularity on coarsest grid
+ imm=im1
+ if(nman.eq.1.and.k.eq.1) then
+ imm=im-2
+ q(im1,2)=0.
+ gam(im1)=0.
+ endif
+c
+c**
+c** Moved calculation of loop 229 from loop 230 for vectorization
+c** on vector machines (ie. Cray)
+c** By: John Towns 2/6/92
+c**
+ do 229 j=jred,jm1,2
+ 229 q(2,j)=q(2,j)/ac(2,j)
+
+c**
+c** Changed bet=(quantity) to bet=1./(quantity) to trade two divisions
+c** for one division and two multiplies (more efficient on all
+c** machines)
+c** By: John Towns 2/6/92
+c**
+
+c**
+c** Cray compiler directives to parallelize tridiagonal solve.
+c** By: John Towns 4/13/92
+c**
+
+cmic$ parallel private(bet,gam,i)
+cmic$1shared(ac,ae,aw,q,jred,jm1,im2,imm)
+cmic$ do parallel
+ do 230 j=jred,jm1,2
+ bet=1./ac(2,j)
+ do 240 i=3,imm
+ gam(i)=ae(i-1,j)*bet
+ bet=1./(ac(i,j)-aw(i,j)*gam(i))
+ 240 q(i,j)=(q(i,j)-aw(i,j)*q(i-1,j))*bet
+ do 250 i=im2,2,-1
+ 250 q(i,j)=q(i,j)-gam(i+1)*q(i+1,j)
+ 230 continue
+cmic$ end do
+cmic$ end parallel
+ 5002 endif
+ jrb=5-jrb
+ 109 continue
+ return
+ end
+c----------------------------------------------------------------------
+ subroutine urscal(
+ + ac,aw,as,ae,an,asw,ase,ane,anw,q,f,fc,qc,rc,
+ + im,jm,jmc,ifd,i9,j9,kf,m,jred,tol,rmax,ipc,irc)
+ implicit real*8 (a-h,o-z)
+c----------------------------------------------------------------------
+c Defines the grid kf-1 right hand side, fc, as the restriction of the
+c grid kf residual. The restriction operator is the transpose of the
+c interpolation operator. Note: The grid kf residual is zero at the
+c black lines (j-direction) as a result of red/black relaxation.
+c Thus, the restriction is simple injection. The initial guess, qc,
+c for the coarse grid correction equation is set to zero. The
+c maximum norm of the residual is calculated and returned in rmax.
+c......................................................................
+ real*8 ac(im,jm),aw(im,jm),as(im,jm),ae(im,jm),an(im,jm),
+ + asw(i9,j9),ase(i9,j9),ane(i9,j9),anw(i9,j9)
+ real*8 f(im,jm),q(im,jm),fc(im,jmc),qc(im,jmc)
+ real*8 rc(im,jmc)
+ rmax=0.
+ im1=im-1
+ jm1=jm-1
+ jmc1=jmc-1
+ jc=1
+ do 10 j=jred,jm1,2
+ jc=jc+1
+ do 10 i=2,im1
+ 10 fc(i,jc)=f(i,j)-as(i,j)*q(i,j-1)-an(i,j)*q(i,j+1)-
+ + aw(i,j)*q(i-1,j)-ae(i,j)*q(i+1,j)-ac(i,j)*q(i,j)
+ 1000 if(ifd.eq.9) then
+ jc=1
+ do 20 j=jred,jm1,2
+ jc=jc+1
+ do 20 i=2,im1
+ 20 fc(i,jc)=fc(i,jc)-asw(i,j)*q(i-1,j-1)-ane(i,j)*q(i+1,j+1)-
+ + ase(i,j)*q(i+1,j-1)-anw(i,j)*q(i-1,j+1)
+ 1001 endif
+c zero out qc as initial guess
+ do 25 jc=1,jmc
+ do 25 i=1,im
+ 25 qc(i,jc)=0.
+c if kf=m calculate residual norm
+ 2000 if((kf.eq.m.and.tol.ge.0.).or.tol.eq.-.5) then
+ do 30 jc=2,jmc1
+ do 30 i=2,im1
+ resmax=abs(fc(i,jc))
+ 30 if(resmax.gt.rmax) rmax=resmax
+ 2001 endif
+c weight rhs if irc.ge.1
+ 3000 if(irc.eq.1.and.ipc.ge.1) then
+ do 40 jc=2,jmc1
+ do 40 i=2,im1
+ 40 fc(i,jc)=rc(i,jc)*fc(i,jc)
+ 3001 endif
+c
+ return
+ end
+c----------------------------------------------------------------------
+ subroutine ursrhs(f,fc,ru,rd,rc,im,jm,jmc,m,kf,jred,irc)
+ implicit real*8 (a-h,o-z)
+c----------------------------------------------------------------------
+c Restricts the right hand side vector on grid kf onto grid kf-1 when
+c the full multigrid (ifmg>0) option is used. The restriction operator
+c is NOT necessarily the transpose of the interpolation operator.
+c......................................................................
+ real*8 f(im,jm),fc(im,jmc),ru(im,jmc),rd(im,jmc),rc(im,jmc)
+ jm1=jm-1
+ im1=im-1
+ jc=1
+ 1000 if(irc.eq.0) then
+ do 10 j=jred,jm1,2
+ jc=jc+1
+ do 10 i=2,im1
+ 10 fc(i,jc)=ru(i,jc-1)*f(i,j-1)+rd(i,jc)*f(i,j+1)+f(i,j)
+ 1001 else
+ do 20 j=jred,jm1,2
+ jc=jc+1
+ do 20 i=2,im1
+ 20 fc(i,jc)=ru(i,jc-1)*f(i,j-1)+rd(i,jc)*f(i,j+1)+
+ + rc(i,jc)*f(i,j)
+ 1002 endif
+ return
+ end
+c----------------------------------------------------------------------
+ subroutine useta(
+ + ac,aw,as,ae,an,asw,ase,ane,anw,acc,awc,asc,aec,
+ + anc,aswc,asec,anec,anwc,pu,pd,pc,ru,rd,rc,qw,fw,gam,
+ + im,jm,jmc,ifd,i9,j9,nman,kf,m,jred,ipc,irc,irurd)
+ implicit real*8 (a-h,o-z)
+c----------------------------------------------------------------------
+cdir$ noinline
+c Calculates the interpolation coefficients from grid kf-1 to
+c grid kf and the coarse grid operator on grid kf-1.
+c** INPUT -
+c ac - anw = fine grid (kf) array stencil coeficients
+c m= total number of grids
+c kf= grid number of the fine grid
+c ifd= the size of the fine grid stencil (= 5 or 9)
+c i9,j9= the i,j-dimensions of the arrays asw,ase
+c qw,fw= coarse grid portions of q and f used for work arrays here
+c (See comments in MGSS2 for details)
+c** OUTPUT -
+c acc - anwc = coarse grid (kf-1) array stencil coeficients
+c pu,pd= arrays of interpolation coefficients from grid kf-1
+c to grid kf
+c .....................................................................
+ real*8 ac(im,jm),aw(im,jm),as(im,jm),ae(im,jm),an(im,jm),
+ + asw(i9,j9),ase(i9,j9),ane(i9,j9),anw(i9,j9),
+ + ru(im,jmc),rd(im,jmc),rc(im,jmc),
+ + pu(im,jmc),pd(im,jmc),pc(im,jmc),gam(im)
+ real*8 acc(im,jmc),awc(im,jmc),asc(im,jmc),aec(im,jmc),
+ + anc(im,jmc),aswc(im,jmc),asec(im,jmc),anec(im,jmc),anwc(im,jmc)
+ real*8 qw(im,jm),fw(im,jm)
+ common /io/ linp,lout
+ common /prsol/ iprsol
+c
+ pcscale=.001
+c
+ im1=im-1
+ jm1=jm-1
+ jmc1=jmc-1
+ jblack=5-jred
+c zeroing out connections to fictitious points
+ do 1 j=1,jm
+ do 2 i=1,im,im1
+ ac(i,j)=0.
+ aw(i,j)=0.
+ as(i,j)=0.
+ ae(i,j)=0.
+ 2 an(i,j)=0.
+ aw(2,j)=0.
+ 1 ae(im1,j)=0.
+ do 3 i=1,im
+ do 4 j=1,jm,jm1
+ ac(i,j)=0.
+ aw(i,j)=0.
+ as(i,j)=0.
+ ae(i,j)=0.
+ 4 an(i,j)=0.
+ as(i,2)=0.
+ 3 an(i,jm1)=0.
+ 1000 if(ifd.eq.9) then
+ do 5 j=1,jm
+ do 6 i=1,im,im1
+ asw(i,j)=0.
+ ase(i,j)=0.
+ ane(i,j)=0.
+ 6 anw(i,j)=0.
+ asw(2,j)=0.
+ anw(2,j)=0.
+ ase(im1,j)=0.
+ 5 ane(im1,j)=0.
+ do 7 i=1,im
+ do 8 j=1,jm,jm1
+ asw(i,j)=0.
+ ase(i,j)=0.
+ ane(i,j)=0.
+ 8 anw(i,j)=0.
+ ase(i,2)=0.
+ asw(i,2)=0.
+ ane(i,jm1)=0.
+ 7 anw(i,jm1)=0.
+ 1001 endif
+c
+ do 9 jc=1,jmc
+ do 9 i=1,im
+ pc(i,jc)=0.
+ pu(i,jc)=0.
+ pd(i,jc)=0.
+ rc(i,jc)=0.
+ ru(i,jc)=0.
+ 9 rd(i,jc)=0.
+c
+c calculation of interpolation coeficients
+c
+
+c define pc
+ 2000 if(ipc.ge.1) then
+ do 20 j=2,jm1
+ do 20 i=2,im1
+ fw(i,j)=0.
+ 20 qw(i,j)=1.
+c
+ call urelax(ac,aw,as,ae,an,asw,ase,ane,anw,fw,qw,gam,im,jm,
+ + i9,j9,ifd,nman,kf,m,jred,ipc,1)
+c scale pc
+ pcmax=0.
+ jc=1
+ do 40 j=jred,jm1,2
+ jc=jc+1
+ do 40 i=2,im1
+ pc(i,jc)=qw(i,j)
+ 40 pcmax=max(pcmax,abs(qw(i,j)))
+ do 50 jc=2,jmc1
+ do 50 i=2,im1
+ if(pc(i,jc).eq.0.) pc(i,jc)=pcscale
+ 50 pc(i,jc)=pc(i,jc)/pcmax
+c
+ 2001 else
+ do 55 jc=2,jmc1
+ do 55 i=2,im1
+ 55 pc(i,jc)=1.
+ 2002 endif
+c
+c define pu
+ jc=3-jred
+ do 60 j=jblack,jm1,2
+ jc=jc+1
+ 4000 if(ipc.eq.0) then
+ do 70 i=2,im1
+ 70 qw(i,j)=-an(i,j)
+ 5000 if(ifd.eq.9) then
+ do 80 i=2,im1
+ 80 qw(i,j)=qw(i,j)-ane(i,j)-anw(i,j)
+ 5001 endif
+ 4001 else
+ do 90 i=2,im1
+ 90 qw(i,j)=-an(i,j)*pc(i,jc+1)
+ 6000 if(ifd.eq.9) then
+ do 100 i=2,im1
+ 100 qw(i,j)=qw(i,j)-ane(i,j)*pc(i+1,jc+1)-anw(i,j)*pc(i-1,jc+1)
+ 6001 endif
+ 4002 endif
+ 60 continue
+c solve for pu
+ call urelax(ac,aw,as,ae,an,asw,ase,ane,anw,fw,qw,gam,im,jm,
+ + i9,j9,ifd,nman,kf,m,jred,ipc,2)
+c
+
+ jc=3-jred
+ do 102 j=jblack,jm1,2
+ jc=jc+1
+ 3020 if(j.lt.jm1) then
+ do 103 i=2,im1
+ 103 pu(i,jc)=qw(i,j)
+ 3021 endif
+ 102 continue
+c
+c define pd
+ jc=3-jred
+ do 106 j=jblack,jm1,2
+ jc=jc+1
+ 8000 if(ipc.eq.0) then
+ do 130 i=2,im1
+ 130 qw(i,j)=-as(i,j)
+ 9000 if(ifd.eq.9) then
+ do 140 i=2,im1
+ 140 qw(i,j)=qw(i,j)-ase(i,j)-asw(i,j)
+ 9001 endif
+c
+ 8001 else
+c
+ do 150 i=2,im1
+ 150 qw(i,j)=-as(i,j)*pc(i,jc)
+ 1100 if(ifd.eq.9) then
+ do 160 i=2,im1
+ 160 qw(i,j)=qw(i,j)-ase(i,j)*pc(i+1,jc)-asw(i,j)*pc(i-1,jc)
+ 1101 endif
+ 8002 endif
+ 106 continue
+c solve for pd
+ call urelax(ac,aw,as,ae,an,asw,ase,ane,anw,fw,qw,gam,im,jm,
+ + i9,j9,ifd,nman,kf,m,jred,ipc,2)
+c
+ jc=3-jred
+ do 105 j=jblack,jm1,2
+ jc=jc+1
+ 7010 if(j.gt.2) then
+ do 104 i=2,im1
+ 104 pd(i,jc)=qw(i,j)
+ 7011 endif
+ 105 continue
+c
+c define restriction operator
+c
+c define rc
+ 1200 if(irc.eq.1) then
+ do 500 jc=2,jmc1
+ do 500 i=2,im1
+ 500 rc(i,jc)=pc(i,jc)
+ else
+ do 502 jc=2,jmc1
+ do 502 i=2,im1
+ 502 rc(i,jc)=1.
+ 1201 endif
+c
+c compute qw = -Cb(inv) * eb*
+ 1300 if(irurd.ge.1) then
+ jc=3-jred
+ 3300 if(irurd.eq.1) then
+ do 560 j=jblack,jm1,2
+ jc=jc+1
+ do 560 i=2,im1
+ 560 qw(i,j)=1.
+ 3301 elseif(irurd.eq.2) then
+ do 561 j=jblack,jm1,2
+ jc=jc+1
+ do 561 i=2,im1
+ 561 qw(i,j)=(pd(i,jc)*pc(i,jc)+pu(i,jc)*pc(i,jc+1))
+ 3302 endif
+c
+ call urelax(ac,aw,as,ae,an,asw,ase,ane,anw,fw,qw,gam,im,jm,
+ + i9,j9,ifd,nman,kf,m,jred,ipc,2)
+c
+ jc=3-jred
+ do 566 j=jblack,jm1,2
+ jc=jc+1
+c compute ru = -b(j+1) * qw
+ 1400 if(j.lt.jm1) then
+ do 570 i=2,im1
+ 570 ru(i,jc)=-as(i,j+1)*qw(i,j)
+ 1500 if(ifd.eq.9) then
+ do 580 i=2,im1
+ 580 ru(i,jc)=ru(i,jc)-ase(i,j+1)*qw(i+1,j)-asw(i,j+1)*qw(i-1,j)
+ 1501 endif
+ 1401 endif
+c compute rd = -a(j-1) * c(j)(inv) * qw
+ 1600 if(j.gt.2) then
+ do 650 i=2,im1
+ 650 rd(i,jc)=-an(i,j-1)*qw(i,j)
+ 1700 if(ifd.eq.9) then
+ do 660 i=2,im1
+ 660 rd(i,jc)=rd(i,jc)-ane(i,j-1)*qw(i+1,j)-anw(i,j-1)*qw(i-1,j)
+ 1701 endif
+ 1601 endif
+ 566 continue
+c
+ 1301 else
+c else set ru=pu and rd=pd
+ jc=3-jred
+ do 670 j=jblack,jm1,2
+ jc=jc+1
+ do 670 i=2,im1
+ ru(i,jc)=pu(i,jc)
+ 670 rd(i,jc)=pd(i,jc)
+ 1303 endif
+c
+c calculating the coarse grid operator
+c
+ 1800 if(ipc+irc+irurd.eq.0) then
+ j=jred-2
+ do 200 jc=2,jmc1
+ j=j+2
+ do 200 i=2,im1
+ acc(i,jc)=ac(i,j)+an(i,j-1)*pu(i,jc-1)+as(i,j+1)*pd(i,jc)+
+ + pu(i,jc-1)*(as(i,j)+ac(i,j-1)*pu(i,jc-1))+
+ + pd(i,jc)*(an(i,j)+ac(i,j+1)*pd(i,jc))
+ awc(i,jc)=aw(i,j)+pd(i-1,jc)*aw(i,j+1)*pd(i,jc)+pu(i-1,jc-1)*
+ + aw(i,j-1)*pu(i,jc-1)
+ asc(i,jc)=as(i,j)*pd(i,jc-1)+pu(i,jc-1)*(as(i,j-1)+
+ + ac(i,j-1)*pd(i,jc-1))
+ aec(i,jc)=ae(i,j)+pd(i+1,jc)*ae(i,j+1)*pd(i,jc)+pu(i+1,jc-1)*
+ + ae(i,j-1)*pu(i,jc-1)
+ anc(i,jc)=an(i,j)*pu(i,jc)+pd(i,jc)*(an(i,j+1)+
+ + ac(i,j+1)*pu(i,jc))
+ aswc(i,jc)=pd(i-1,jc-1)*aw(i,j-1)*pu(i,jc-1)
+ asec(i,jc)=pd(i+1,jc-1)*ae(i,j-1)*pu(i,jc-1)
+ anec(i,jc)=pu(i+1,jc)*ae(i,j+1)*pd(i,jc)
+ 200 anwc(i,jc)=pu(i-1,jc)*aw(i,j+1)*pd(i,jc)
+ 1900 if(ifd.eq.9) then
+ j=jred-2
+ do 210 jc=2,jmc1
+ j=j+2
+ do 210 i=2,im1
+ awc(i,jc)=awc(i,jc)+asw(i,j+1)*pd(i,jc)+anw(i,j-1)*pu(i,jc-1)+
+ + pd(i-1,jc)*anw(i,j)+pu(i-1,jc-1)*asw(i,j)
+ aec(i,jc)=aec(i,jc)+ase(i,j+1)*pd(i,jc)+ane(i,j-1)*pu(i,jc-1)+
+ + pd(i+1,jc)*ane(i,j)+pu(i+1,jc-1)*ase(i,j)
+ aswc(i,jc)=aswc(i,jc)+asw(i,j-1)*pu(i,jc-1)+pd(i-1,jc-1)*asw(i,j)
+ asec(i,jc)=asec(i,jc)+ase(i,j-1)*pu(i,jc-1)+pd(i+1,jc-1)*ase(i,j)
+ anec(i,jc)=anec(i,jc)+ane(i,j+1)*pd(i,jc)+pu(i+1,jc)*ane(i,j)
+ 210 anwc(i,jc)=anwc(i,jc)+anw(i,j+1)*pd(i,jc)+pu(i-1,jc)*anw(i,j)
+ 1901 endif
+c
+ 1801 else
+c
+ j=jred-2
+ do 300 jc=2,jmc1
+ j=j+2
+ do 300 i=2,im1
+ acc(i,jc)=rc(i,jc)*(ac(i,j)*pc(i,jc)+
+ + as(i,j)*pu(i,jc-1)+
+ + an(i,j)*pd(i,jc))+
+ + ru(i,jc-1)*(ac(i,j-1)*pu(i,jc-1)+
+ + an(i,j-1)*pc(i,jc))+
+ + rd(i,jc)*(ac(i,j+1)*pd(i,jc)+
+ + as(i,j+1)*pc(i,jc))
+ awc(i,jc)=rc(i,jc)*aw(i,j)*pc(i-1,jc)+
+ + rd(i,jc)*aw(i,j+1)*pd(i-1,jc)+
+ + ru(i,jc-1)*aw(i,j-1)*pu(i-1,jc-1)
+ asc(i,jc)=rc(i,jc)*as(i,j)*pd(i,jc-1)+
+ + ru(i,jc-1)*(ac(i,j-1)*pd(i,jc-1)+
+ + as(i,j-1)*pc(i,jc-1))
+ aec(i,jc)=rc(i,jc)*ae(i,j)*pc(i+1,jc)+
+ + rd(i,jc)*ae(i,j+1)*pd(i+1,jc)+
+ + ru(i,jc-1)*ae(i,j-1)*pu(i+1,jc-1)
+ anc(i,jc)=rc(i,jc)*an(i,j)*pu(i,jc)+
+ + rd(i,jc)*(ac(i,j+1)*pu(i,jc)+
+ + an(i,j+1)*pc(i,jc+1))
+ aswc(i,jc)=ru(i,jc-1)*aw(i,j-1)*pd(i-1,jc-1)
+ asec(i,jc)=ru(i,jc-1)*ae(i,j-1)*pd(i+1,jc-1)
+ anec(i,jc)=rd(i,jc)*ae(i,j+1)*pu(i+1,jc)
+ 300 anwc(i,jc)=rd(i,jc)*aw(i,j+1)*pu(i-1,jc)
+ 2100 if(ifd.eq.9) then
+ j=jred-2
+ do 310 jc=2,jmc1
+ j=j+2
+ do 310 i=2,im1
+ awc(i,jc)=awc(i,jc)+(rd(i,jc)*asw(i,j+1)+
+ + ru(i,jc-1)*anw(i,j-1))*pc(i-1,jc)+
+ + rc(i,jc)*(anw(i,j)*pd(i-1,jc)+
+ + asw(i,j)*pu(i-1,jc-1))
+ aec(i,jc)=aec(i,jc)+(rd(i,jc)*ase(i,j+1)+
+ + ru(i,jc-1)*ane(i,j-1))*pc(i+1,jc)+
+ + rc(i,jc)*(ane(i,j)*pd(i+1,jc)+
+ + ase(i,j)*pu(i+1,jc-1))
+ aswc(i,jc)=aswc(i,jc)+ru(i,jc-1)*asw(i,j-1)*pc(i-1,jc-1)+
+ + rc(i,jc)*asw(i,j)*pd(i-1,jc-1)
+ asec(i,jc)=asec(i,jc)+ru(i,jc-1)*ase(i,j-1)*pc(i+1,jc-1)+
+ + rc(i,jc)*ase(i,j)*pd(i+1,jc-1)
+ anec(i,jc)=anec(i,jc)+rd(i,jc)*ane(i,j+1)*pc(i+1,jc+1)+
+ + rc(i,jc)*ane(i,j)*pu(i+1,jc)
+ 310 anwc(i,jc)=anwc(i,jc)+rd(i,jc)*anw(i,j+1)*pc(i-1,jc+1)+
+ + rc(i,jc)*anw(i,j)*pu(i-1,jc)
+ 2101 endif
+ 1802 endif
+cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
+ if(iprsol.eq.2.and.kf.eq.m.and.ifd.eq.5) then
+ do 111 j=2,jm1
+ do 111 i=2,im1
+ write(lout,1010) im,jm,an(i,j)
+ write(lout,1012) i,j,aw(i,j),ac(i,j),ae(i,j)
+ 111 write(lout,1011) as(i,j)
+ 1010 format(2(1x,i2),14x,f12.5)
+ 1011 format(20x,f12.5)
+ 1012 format(2(1x,i2),3(1x,f12.5))
+ endif
+ if(iprsol.eq.2.and.kf.eq.m.and.ifd.eq.9) then
+ do 115 j=2,jm1
+ do 115 i=2,im1
+ write(lout,1017) im,jm,anw(i,j),an(i,j),ane(i,j)
+ write(lout,1017) i,j,aw(i,j),ac(i,j),ae(i,j)
+ 115 write(lout,1016) asw(i,j),as(i,j),ase(i,j)
+ 1016 format(6x,3(1x,f12.5))
+ 1017 format(2(1x,i2),3(1x,f12.5))
+ endif
+cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
+ if(iprsol.eq.2.or.iprsol.eq.3) then
+ write(lout,*) ' coarse grid kf-1 ==',kf-1
+ do 211 jc=2,jmc1
+ do 211 i=2,im1
+ write(lout,2015) im,jmc,anwc(i,jc),anc(i,jc),anec(i,jc),pd(i,jc),
+ + rd(i,jc)
+ write(lout,2013) i,jc,awc(i,jc),acc(i,jc),aec(i,jc),pc(i,jc),
+ + rc(i,jc)
+ write(lout,2014) aswc(i,jc),asc(i,jc),asec(i,jc),pu(i,jc-1),
+ + ru(i,jc-1)
+ 211 write(lout,*) ' m, kf-1=',m,kf-1
+ 2014 format(9x,3(1x,f11.5),3x,2(f11.5,1x))
+ 2013 format(3x,2(1x,i2),3(1x,f11.5),3x,2(f11.5,1x))
+ 2015 format(3x,2(1x,i2),3(1x,f11.5),3x,2(f11.5,1x))
+ endif
+cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
+ return
+ end
+c**********************************************************************
+ subroutine uoutpt(q,im,jm)
+ implicit real*8 (a-h,o-z)
+c**********************************************************************
+c Sample output subroutine. Prints out the values of q at the
+c interior points of the finest grid.
+c**********************************************************************
+ common /io/ linp,lout
+ real*8 q(im,jm)
+ im1=im-1
+ jm1=jm-1
+ ie=1
+ 20 ib=ie+1
+ ie=ib+min0(5,im1-ib)
+ do 10 j=jm1,2,-1
+ 10 write(lout,100) j,(q(i,j),i=ib,ie)
+ if(ie.lt.im1) go to 20
+ 100 format(1x,i2,1x,6(1x,f10.4))
+ return
+ end
+c**********************************************************************
+ subroutine uputf(ac,aw,as,ae,an,f,nx,ny,
+ + lo,nxd,nyd,i32su)
+ implicit real*8 (a-h,o-z)
+c**********************************************************************
+ real*8 ac(lo:nxd,lo:nyd),aw(lo:nxd,lo:nyd),
+ + as(lo:nxd,lo:nyd),ae(lo:nxd,lo:nyd),
+ + an(lo:nxd,lo:nyd),f(lo:nxd,lo:nyd)
+ real*8 a(20),b(20),ab(20)
+ integer il(20),ir(20),jb(20),jt(20)
+ integer ibc(4)
+ common /io/ linp,lout
+c
+c dcell is the value assigned to the diagonal element of a dead
+c cell, i.e. a cell that has 0 conections to all its neighbors.
+c icendif determines the differencing scheme for the first order terms
+c icendif=0 - central differencing, =1 - forward differencing.
+c
+ dcell = 1.
+ do 321 j=lo,nyd
+ do 321 i=lo,nxd
+ ac(i,j)=0.
+ aw(i,j)=0.
+ as(i,j)=0.
+ ae(i,j)=0.
+ 321 an(i,j)=0.
+c
+ nx1=nx+1
+ ny1=ny+1
+ read(linp,*) icendif
+ read(linp,*) ibc(1),ibc(2),ibc(3),ibc(4)
+ read(linp,*) nreg
+ hx=1./nx
+ hy=1./ny
+ hx2=hx*hx
+ hy2=hy*hy
+ hxy2=hx2*hy2
+ dcell=hxy2*dcell
+ write(lout,1011) ibc(1),ibc(2),ibc(3),ibc(4),icendif,dcell,hx,hy
+ 1011 format(' ibc_l ibc_b ibc_r ibc_t icendiff dcell hx hy'/
+ + 4x,i1,5x,i1,5x,i1,5x,i1,7x,i1,6x,e8.2,2x,f6.4,2x,f6.4/)
+c
+ do 10 irg=1,nreg
+ read(linp,*) il(irg),ir(irg),jb(irg),jt(irg)
+ read(linp,*) xk,yk,sreg,freg
+ read(linp,*) a(irg),b(irg),ab(irg)
+ write(lout,1000) il(irg),ir(irg),jb(irg),jt(irg),xk,yk,sreg,freg
+ 1000 format(1x,i3,',',i3,' X ',i3,',',i3,2x,1pe8.2,1x,1pe8.2,2x,
+ + 1pe8.1,1x,1pe8.1)
+ write(lout,1001) a(irg),b(irg),ab(irg)
+ 1001 format(17x,' a=',1pe10.3,' b=',1pe10.3,' ab=',1pe10.3)
+ xk=xk*hy2
+ yk=yk*hx2
+ sreg=sreg*hxy2
+ freg=freg*hxy2
+ a(irg)=a(irg)*hx2*hy
+ b(irg)=b(irg)*hx*hy2
+ ab(irg)=ab(irg)*hx*hy
+ if(icendif.eq.0) then
+ a(irg)=a(irg)/2.
+ b(irg)=b(irg)/2.
+ ab(irg)=ab(irg)/4.
+ endif
+ if(il(irg).eq.1) il(irg)=0
+ if(ir(irg).eq.nx) ir(irg)=nx1
+ if(jb(irg).eq.1) jb(irg)=0
+ if(jt(irg).eq.ny) jt(irg)=ny1
+ do 20 i=il(irg),ir(irg)
+ do 20 j=jb(irg),jt(irg)
+ aw(i,j)=xk
+ as(i,j)=yk
+ ac(i,j)=sreg
+ 20 f(i,j)=freg
+ 10 continue
+ write(lout,*) ' - - - - - - - - - - - - - - - - - - - - - - - - -'
+c defining coeficients by harmonic averaging
+ do 30 i=1,nx
+ asio=as(i,0)
+ do 30 j=1,ny1
+ aa=as(i,j)*asio
+ if(aa.gt.0.) then
+ t=2.*aa/(as(i,j)+asio)
+ asio=as(i,j)
+ as(i,j)=t
+ else
+ asio=as(i,j)
+ as(i,j)=0.
+ endif
+ 30 continue
+ do 40 j=1,ny
+ awoj=aw(0,j)
+ do 40 i=1,nx1
+ aa=aw(i,j)*awoj
+ if(aa.gt.0.) then
+ t=2.*aa/(aw(i,j)+awoj)
+ awoj=aw(i,j)
+ aw(i,j)=t
+ else
+ awoj=aw(i,j)
+ aw(i,j)=0.
+ endif
+ 40 continue
+ do 45 i=0,nx
+ do 45 j=0,ny
+ ae(i,j)=aw(i+1,j)
+ 45 an(i,j)=as(i,j+1)
+ do 50 i=1,nx
+ do 50 j=1,ny
+ ac(i,j)=ac(i,j)-aw(i,j)-as(i,j)-ae(i,j)-
+ + an(i,j)
+ 50 if(ac(i,j).eq.0.) ac(i,j)=dcell
+c adding on the unsymmetric terms
+ do 51 irg=1,nreg
+c icendif=0 ==> central diff g
+ if(icendif.eq.0) then
+ do 52 i=il(irg),ir(irg)
+ do 52 j=jb(irg),jt(irg)
+ aw(i,j)=aw(i,j)-a(irg)
+ ae(i,j)=ae(i,j)+a(irg)
+ an(i,j)=an(i,j)+b(irg)
+ 52 as(i,j)=as(i,j)-b(irg)
+c icendif=1 ==> upstream diff s
+ elseif(icendif.eq.1) then
+ do 54 i=il(irg),ir(irg)
+ do 54 j=jb(irg),jt(irg)
+ ac(i,j)=ac(i,j)-a(irg)
+ ae(i,j)=ae(i,j)+a(irg)
+ an(i,j)=an(i,j)+b(irg)
+ 54 ac(i,j)=ac(i,j)-b(irg)
+ endif
+ 51 continue
+c set boundary conditions for 5 point operator
+ do 60 j=1,ny
+c left boundary
+ ae(0,j)=aw(1,j)
+ if(ibc(1).eq.1) then
+ ac(0,j)=aw(1,j)
+ f(0,j)=2.*aw(1,j)*0.0
+ elseif(ibc(1).eq.2) then
+ ac(0,j)=-aw(1,j)
+ f(0,j)=hx*aw(1,j)*0.0
+ endif
+ if(ac(0,j).eq.0.) then
+ ae(0,j)=0.
+ ac(0,j)=dcell
+ f(0,j)=0.
+ endif
+c right boundary
+ aw(nx1,j)=ae(nx,j)
+ if(ibc(3).eq.1) then
+ ac(nx1,j)=ae(nx,j)
+ f(nx1,j)=2.*ae(nx,j)*0.
+ elseif(ibc(3).eq.2) then
+ ac(nx1,j)=-ae(nx,j)
+ f(nx1,j)=hx*ae(nx,j)*0.
+ endif
+ if(ac(nx1,j).eq.0.) then
+ aw(nx1,j)=0.
+ ac(nx1,j)=dcell
+ f(nx1,j)=0.
+ endif
+ 60 continue
+c
+ do 80 i=1,nx
+c lower boundary
+ an(i,0)=as(i,1)
+ if(ibc(2).eq.1) then
+ ac(i,0)=as(i,1)
+ f(i,0)=2.*as(i,1)*0.
+ elseif(ibc(2).eq.2) then
+ ac(i,0)=-as(i,1)
+ f(i,0)=hy*as(i,1)*0.
+ endif
+ if(ac(i,0).eq.0.) then
+ an(i,0)=0.
+ ac(i,0)=dcell
+ f(i,0)=0.
+ endif
+c upper boundary
+ as(i,ny1)=an(i,ny)
+ if(ibc(4).eq.1) then
+ ac(i,ny1)=an(i,ny)
+ f(i,ny1)=2.*an(i,ny)*0.
+ elseif(ibc(4).eq.2) then
+ ac(i,ny1)=-an(i,ny)
+ f(i,ny1)=2.*an(i,ny)*0.
+ endif
+ if(ac(i,ny1).eq.0.) then
+ as(i,ny1)=0.
+ ac(i,ny1)=dcell
+ f(i,ny1)=0.
+ endif
+ 80 continue
+c connections between ghost boundary points zeroed
+ do 83 j=1,ny1
+ as(0,j)=0.
+ an(0,j-1)=0.
+ as(nx1,j)=0.
+ 83 an(nx1,j-1)=0.
+ do 86 i=1,nx1
+ aw(i,0)=0.
+ ae(i-1,0)=0.
+ aw(i,ny1)=0.
+ 86 ae(i-1,ny1)=0.
+c corner stencils and rhs defined
+ if(i32su.eq.32) then
+ do 90 j=0,ny1,ny1
+ do 90 i=0,nx1,nx1
+ ac(i,j)=dcell
+ aw(i,j)=0.
+ ae(i,j)=0.
+ as(i,j)=0.
+ an(i,j)=0.
+ 90 f(i,j)=0.
+ endif
+c i32su=22 - boundary conditions absorbed
+ if(i32su.eq.22) then
+ do 100 j=1,ny
+ awac=aw(1,j)/ac(0,j)
+ ac(1,j)=ac(1,j)-awac*ae(0,j)
+ aw(1,j)=0.
+ f(1,j)=f(1,j)-awac*f(0,j)
+ ac(0,j)=0.
+ ae(0,j)=0.
+ f(0,j)=0
+ awac=aw(nx1,j)/ac(nx1,j)
+ ac(nx,j)=ac(nx,j)-awac*ae(nx1,j)
+ ae(nx,j)=0.
+ f(nx,j)=f(nx,j)-awac*f(nx1,j)
+ ac(nx1,j)=0.
+ aw(nx1,j)=0.
+ 100 f(nx1,j)=0.
+c
+ do 110 i=1,nx
+ asac=as(i,1)/ac(i,0)
+ ac(i,1)=ac(i,1)-asac*an(i,0)
+ as(i,1)=0.
+ f(i,1)=f(i,1)-asac*f(i,0)
+ ac(i,0)=0.
+ an(i,0)=0.
+ f(i,0)=0.
+ anac=an(i,ny)/ac(i,ny1)
+ ac(i,ny)=ac(i,ny)-anac*as(i,ny1)
+ an(i,ny)=0.
+ f(i,ny)=f(i,ny)-anac*f(i,ny1)
+ ac(i,ny1)=0.
+ as(i,ny1)=0.
+ 110 f(i,ny1)=0.
+ endif
+ return
+ end
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-21 21:26:37 +0000