From cab79658b8e19947c3b4058f75961f01a055ef2e Mon Sep 17 00:00:00 2001 From: allen Date: Mon, 1 Nov 1999 11:27:16 +0000 Subject: 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 --- src/AxiBrillBHIVP.F | 339 ++++++++++++ src/IDAxiBrillBH.F | 339 ++++++++++++ src/bhbrill.m | 126 +++++ src/bhbrill.x | 60 +++ src/gij.x | 133 +++++ src/interp2.F | 213 ++++++++ src/make.code.defn | 9 + src/mg59p.F | 896 +++++++++++++++++++++++++++++++ src/psi_1st_deriv.x | 18 + src/psi_2nd_deriv.x | 51 ++ src/shmgp.F77 | 1455 +++++++++++++++++++++++++++++++++++++++++++++++++++ 11 files changed, 3639 insertions(+) create mode 100644 src/AxiBrillBHIVP.F create mode 100644 src/IDAxiBrillBH.F create mode 100644 src/bhbrill.m create mode 100644 src/bhbrill.x create mode 100644 src/gij.x create mode 100644 src/interp2.F create mode 100644 src/make.code.defn create mode 100644 src/mg59p.F create mode 100644 src/psi_1st_deriv.x create mode 100644 src/psi_2nd_deriv.x create mode 100644 src/shmgp.F77 (limited to 'src') 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

+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 -- cgit v1.2.3