C C * * * * * * * * * * * * * C * * C * HS3CRT * C * * C * * * * * * * * * * * * * C C C (VERSION 1, AUGUST 1985) C C A VECTORIZED PACKAGE OF FORTRAN SUBPROGRAMS FOR THE SOLUTION C OF A THREE-DIMENSIONAL HELMHOLTZ EQUATION ON A STAGGERED GRID C C BY C C ROLAND A. SWEET C SCIENTIFIC COMPUTING DIVISION C NATIONAL BUREAU OF STANDARDS C BOULDER, COLORADO 80303 C C C * * * * * * * * * * * * * C * * C * PURPOSE * C * * C * * * * * * * * * * * * * C C C THIS PACKAGE SOLVES THE STANDARD SEVEN-POINT FINITE DIFFERENCE C APPROXIMATION ON A STAGGERED GRID TO THE HELMHOLTZ EQUATION IN C CARTESIAN COORDINATES C C (D/DX)(DU/DX)+(D/DY)(DU/DY)+(D/DZ)(DU/DZ) + LAMBDA*U = F(X,Y,Z) C C WITH A VARIETY OF POSSIBLE BOUNDARY CONDITIONS. C C * * * * * * * * * * * * * * * * * * * * C * * C * DESCRIPTION OF PACKAGE * C * * C * * * * * * * * * * * * * * * * * * * * C C C THIS PACKAGE CONSISTS OF TWO USER-CALLABLE FORTRAN SUBROUTINES C FOR THE SOLUTION OF FINITE DIFFERENCE APPROXIMATIONS ON A STAGGERED C (CELL-CENTERED) GRID TO THREE-DIMENSIONAL HELMHOLTZ EQUATIONS. C THE PACKAGE IS USED TO SOLVE A PARTICULAR PROBLEM BY CALLING: C C 1) SUBROUTINE HS3CRI(XS,XF,L,LBDCND,YS,YF,M,MBDCND,ZS,ZF,N,NBDCND, C ELMBDA,LDIMF,MDIMF,IERROR,W) C C THAT DEFINES THE LINEAR OPERATOR REPRESENTING THE FINITE DIFFERENCE C APPROXIMATION, AND C C 2) SUBROUTINE HS3CRT(BDXS,BDXF,BDYS,BDYF,BDZS,BDZF, C LDIMF,MDIMF,F,PERTRB,W) C C THAT DEFINES THE DATA FOR THE LINEAR SYSTEM OF EQUATIONS AND SOLVES C THE SYSTEM BY INVOKING THE SUBROUTINE PACKAGE PSTG3D. C C THE USER MUST CALL SUBROUTINE HS3CRI TO DEFINE THE APPROXIMATION C AND FOLLOW THAT BY A CALL TO SUBROUTINE HS3CRT TO SOLVE A PARTICULAR C PROBLEM. SUBSEQUENT SOLUTIONS TO THE SAME FINITE DIFFERENCE C APPROXIMATION BUT WITH DIFFERENT DATA MAY BE OBTAINED BY ONLY CALLING C SUBROUTINE HS3CRT. SUCH SEPARATION OF TASKS RESULTS IN FASTER C REPEATED SOLUTIONS. C C C C * * * * * * * * * * * * * * * * * * * * * C * * C * PARAMETER DESCRIPTION * C * * C * * * * * * * * * * * * * * * * * * * * * C C C INPUT PARAMETERS C C XS,XF C THE RANGE OF X, I.E. XS .LE. X .LE. XF. XS MUST BE LESS THAN XF. C C L C THE NUMBER OF GRID POINTS IN THE INTERVAL (XS,XF). THE GRID C POINTS IN THE X-DIRECTION ARE GIVEN BY X(I) = XS + (I-0.5)DX FOR C I=1,2,...,L WHERE DX =(XF-XS)/L. L MUST BE AT LEAST 3. C C LBDCND C INDICATES THE TYPE OF BOUNDARY CONDITIONS AT X = XS AND X = XF. C C = 0 IF THE SOLUTION IS PERIODIC IN X, C U(L+I,J,K) = U(I,J,K). C C = 1 IF THE SOLUTION IS SPECIFIED AT X = XS AND X = XF. C C = 2 IF THE SOLUTION IS SPECIFIED AT X = XS AND THE DERIVATIVE C OF THE SOLUTION WITH RESPECT TO X IS SPECIFIED AT X = XF. C C = 3 IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO X IS C SPECIFIED AT X = XS AND X = XF. C C = 4 IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO X IS C SPECIFIED AT X = XS AND THE SOLUTION IS SPECIFIED AT X = XF. C C BDXS C A TWO-DIMENSIONAL ARRAY THAT SPECIFIES THE BOUNDARY VALUES C (IF ANY) OF THE SOLUTION AT X = XS. WHEN LBDCND = 1 OR 2, C C BDXS(J,K) = U(XS,Y(J),Z(K)) , J=1,2,...,M, C K=1,2,...,N. C C WHEN LBDCND = 3 OR 4, C C BDXS(J,K) = (D/DX)U(XS,Y(J),Z(K)) , J=1,2,...,M, C K=1,2,...,N. C C WHEN LBDCND = 0, BDXS IS A DUMMY VARIABLE. BDXS MUST BE C DIMENSIONED MDIMF X N. C C BDXF C A TWO-DIMENSIONAL ARRAY THAT SPECIFIES THE BOUNDARY VALUES OF C THE SOLUTION AT X = XF. WHEN LBDCND = 1 OR 4 C C BDXF(J,K) = U(XF,Y(J),Z(K)) , J=1,2,...,M, C K=1,2,...,N. C C WHEN LBDCND = 2 OR 3 C C BDXF(J,K) = (D/DX)U(XF,Y(J),Z(K)) , J=1,2,...,M, C K=1,2,...,N. C C WHEN LBDCND = 0, BDXF IS A DUMMY VARIABLE. BDXF MUST BE C DIMENSIONED MDIMF X N. C C YS,YF C THE RANGE OF Y, I.E. YS. LE. Y .LE. YF. YS MUST BE LESS C THAN YF. C C M C THE NUMBER OF UNKNOWNS IN THE INTERVAL (YS,YF). THE UNKNOWNS IN C THE Y-DIRECTION ARE GIVEN BY Y(J) = YS + (J-0.5)DY, C J=1,2,...,M, WHERE DY = (YF-YS)/M. M MUST BE AT LEAST 3. C C MBDCND C INDICATES THE TYPE OF BOUNDARY CONDITIONS AT Y = YS C AND Y = YF. C C = 0 IF THE SOLUTION IS PERIODIC IN Y, I.E. C U(I,M+J,K) = U(I,J,K). C C = 1 IF THE SOLUTION IS SPECIFIED AT Y = YS AND Y = YF. C C = 2 IF THE SOLUTION IS SPECIFIED AT Y = YS AND THE DERIVATIVE C OF THE SOLUTION WITH RESPECT TO Y IS SPECIFIED AT Y = YF. C C = 3 IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO Y IS C SPECIFIED AT Y = YS AND Y = YF. C C = 4 IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO Y IS C SPECIFIED AT Y = YS AND THE SOLUTION IS SPECIFIED AT Y = YF. C C BDYS C A TWO-DIMENSIONAL ARRAY THAT SPECIFIES THE BOUNDARY VALUES OF THE C SOLUTION AT Y = YS. WHEN MBDCND = 1 OR 2, C C BDYS(I,K) = U(X(I),YS,Z(K)) , I=1,2,...,L, C K=1,2,...,N. C C WHEN MBDCND = 3 OR 4, C C BDYS(I,K) = (D/DY)U(X(I),YS,Z(K)), I=1,2,...,L, C K=1,2,...,N. C C WHEN MBDCND = 0, BDYS IS A DUMMY VARIABLE. BDYS MUST BE C DIMENSIONED LDIMF X N. C C BDYF C A TWO-DIMENSIONAL ARRAY THAT SPECIFIES THE BOUNDARY VALUES OF THE C SOLUTION AT Y = YF. WHEN MBDCND = 1 OR 4, C C BDYF(I,K) = U(X(I),YF,Z(K)) , I=1,2,...,L, C K=1,2,...,N. C C WHEN MBDCND = 2 OR 3, C C BDYF(I,K) = (D/DY)U(X(I),YF,Z(K)) , I=1,2,...,L, C K=1,2,...,N. C C WHEN MBDCND = 0, BDYF IS A DUMMY VARIABLE. BDYF MUST BE C DIMENSIONED LDIMF X N. C C ZS,ZF C THE RANGE OF X, I.E. ZS. LE. Z .LE. ZF. ZS MUST BE LESS C THAN ZF. C C N C THE NUMBER OF UNKNOWNS IN THE INTERVAL (ZS,ZF). THE UNKNOWNS IN C THE Z-DIRECTION ARE GIVEN BY Z(J) = ZS + (K-0.5)DZ, C K=1,2,...,N, WHERE DZ = (ZF-ZS)/N. N MUST BE AT LEAST 3. C C NBDCND C INDICATES THE TYPE OF BOUNDARY CONDITIONS AT Z = ZS C AND Z = ZF. C C = 0 IF THE SOLUTION IS PERIODIC IN Z, I.E. C U(I,J,N+K) = U(I,J,K). C C = 1 IF THE SOLUTION IS SPECIFIED AT Z = ZS AND Z = ZF. C C = 2 IF THE SOLUTION IS SPECIFIED AT Z = ZS AND THE DERIVATIVE C OF THE SOLUTION WITH RESPECT TO Z IS SPECIFIED AT Z = ZF. C C = 3 IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO Z IS C SPECIFIED AT Z = ZS AND Z = ZF. C C = 4 IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO Z IS C SPECIFIED AT Z = ZS AND THE SOLUTION IS SPECIFIED AT Z = ZF. C C BDZS C A TWO-DIMENSIONAL ARRAY THAT SPECIFIES THE BOUNDARY VALUES OF THE C SOLUTION AT Z = ZS. WHEN NBDCND = 1 OR 2, C C BDZS(I,J) = U(X(I),Y(J),ZS) , I=1,2,...,L, C J=1,2,...,M. C C WHEN NBDCND = 3 OR 4, C C BDZS(I,J) = (D/DZ)U(X(I),Y(J),ZS), I=1,2,...,L, C J=1,2,...,M. C C WHEN NBDCND = 0, BDZS IS A DUMMY VARIABLE. BDZS MUST BE C DIMENSIONED LDIMF X M. C C BDZF C A TWO-DIMENSIONAL ARRAY THAT SPECIFIES THE BOUNDARY VALUES OF THE C SOLUTION AT Z = ZF. WHEN NBDCND = 1 OR 4, C C BDZF(I,J) = U(X(I),Y(J),ZF) , I=1,2,...,L, C J=1,2,...,M. C C WHEN NBDCND = 2 OR 3, C C BDZF(I,J) = (D/DZ)U(X(I),Y(J),ZF) , I=1,2,...,L, C J=1,2,...,M. C C WHEN NBDCND = 0, BDZF IS A DUMMY VARIABLE. BDZF MUST BE C DIMENSIONED LDIMF X M. C C ELMBDA C THE CONSTANT LAMBDA IN THE HELMHOLTZ EQUATION. IF LAMBDA IS C GREATER THAN 0, A SOLUTION MAY NOT EXIST. HOWEVER, HS3CRT WILL C ATTEMPT TO FIND A SOLUTION. C C F C A THREE-DIMENSIONAL ARRAY THAT SPECIFIES THE VALUES OF THE RIGHT C SIDE OF THE HELMHOLTZ EQUATION. FOR I=1,2,...,L, J=1,2,...,M, C AND K=1,2,...,N C C F(I,J,K) = F(X(I),Y(J),Z(K)) . C C F MUST BE DIMENSIONED LDIMF X MDIMF X N. C C LDIMF C THE ROW (OR FIRST) DIMENSION OF THE ARRAYS F, BDYS, BDYF, BDZS, C AND BDZF AS THEY APPEAR IN THE PROGRAM CALLING HS3CRT. THIS C PARAMETER IS USED TO SPECIFY THE VARIABLE DIMENSION OF THESE C ARRAYS. LDIMF MUST BE AT LEAST L. C C MDIMF C THE COLUMN (OR SECOND) DIMENSION OF THE ARRAY F AND THE ROW C (OR FIRST) DIMENSION OF THE ARRAYS BDXS AND BDXF AS THEY APPEAR C IN THE PROGRAM CALLING HS3CRT. THIS PARAMETER IS USED TO SPECIFY C THE VARIABLE DIMENSION OF THESE ARRAYS. MDIMF MUST BE AT C LEAST M. C C W C A ONE-DIMENSIONAL ARRAY THAT MUST BE PROVIDED BY THE USER FOR C WORK SPACE. THE LENGTH OF W MUST BE AT LEAST 2*L*M*N + 6*L + C 5*(M+N) + 46. C C C OUTPUT PARAMETERS C C F C CONTAINS THE SOLUTION U(I,J,K) OF THE FINITE DIFFERENCE C APPROXIMATION FOR THE GRID POINT (X(I),Y(J),Z(K)) FOR C I=1,2,...,L, J=1,2,...,M, AND K=1,2,...,N. C C PERTRB C IF A COMBINATION OF PERIODIC OR DERIVATIVE BOUNDARY CONDITIONS IS C SPECIFIED FOR A POISSON EQUATION (LAMBDA = 0), A SOLUTION MAY NOT C EXIST. PERTRB IS A CONSTANT, CALCULATED AND SUBTRACTED FROM F, C WHICH ENSURES THAT A SOLUTION EXISTS. HS3CRT THEN COMPUTES THIS C SOLUTION, WHICH IS A LEAST SQUARES SOLUTION TO THE ORIGINAL C APPROXIMATION. THIS SOLUTION PLUS ANY CONSTANT IS ALSO A C SOLUTION; HENCE, THE SOLUTION IS NOT UNIQUE. THE VALUE OF PERTRB C SHOULD BE SMALL COMPARED TO THE RIGHT SIDE F. OTHERWISE, A C SOLUTION IS OBTAINED TO AN ESSENTIALLY DIFFERENT PROBLEM. THIS C COMPARISON SHOULD ALWAYS BE MADE TO INSURE THAT A MEANINGFUL C SOLUTION HAS BEEN OBTAINED. C C IERROR C AN ERROR FLAG THAT INDICATES INVALID INPUT PARAMETERS. C EXCEPT FOR NUMBERS 0 AND 10, A SOLUTION IS NOT ATTEMPTED. C C = 0 NO ERROR C C = 1 XS .GT. XF C C = 2 L .LT. 3 C C = 3 LBDCND .NE. 0, 1, 2, 3, OR 4 C C = 4 YS .GT. YF C C = 5 M .LT. 3 C C = 6 MBDCND .NE. 0, 1, 2, 3, OR 4 C C = 7 ZS .GT. ZF C C = 8 N .LT. 3 C C = 9 NBDCND .NE. 0, 1, 2, 3, OR 4 C C = 10 ELMBDA .GT. 0. C C = 11 LDIMF .LT. L C C = 12 MDIMF .LT. M C C SINCE THIS IS THE ONLY MEANS OF INDICATING A POSSIBLY C INCORRECT CALL TO HS3CRI, THE USER SHOULD TEST IERROR AFTER C THE CALL. C C W C CONTAINS INTERMEDIATE QUANTITIES (CALCULATED IN ROUTINE HS3CRI) C THAT MUST NOT BE DESTROYED IF THIS SUBROUTINE WILL BE USED FOR C REPEAT SOLUTIONS. C C C * * * * * * * * * * * * * * * * * * * * * C * * C * PROGRAM SPECIFICATIONS * C * * C * * * * * * * * * * * * * * * * * * * * * C C C DIMENSION OF BDXS(MDIMF,N), BDXF(MDIMF,N), BDYS(LDIMF,N), C ARGUMENTS BDYF(LDIMF,N), BDZS(LDIMF,M), BDZF(LDIMF,M), C F(LDIMF,MDIMF,N), W(2*L*M*N+6*L+5*(M+N)+46) C C LATEST OCTOBER 1, 1984 C REVISION C C SUBPROGRAMS HS3CRT,HS3CRI,PSTG3D(PACKAGE) C REQUIRED VSFFTI(PACKAGE),VSFFT(PACKAGE), C VRFFTI(PACKAGE),VRFFT(PACKAGE) C C SPECIAL NONE C CONDITIONS C C COMMON NONE C BLOCKS C C I/O NONE C C PRECISION SINGLE C C SPECIALIST ROLAND SWEET C C LANGUAGE FORTRAN C C HISTORY WRITTEN BY ROLAND SWEET AT THE NATIONAL BUREAU OF C STANDARDS (BOULDER). C C ALGORITHM THIS SUBROUTINE DEFINES THE FINITE-DIFFERENCE C EQUATIONS (AS DEFINED BELOW), INCORPORATES BOUNDARY C DATA, ADJUSTS THE RIGHT SIDE WHEN THE SYSTEM IS C SINGULAR AND CALLS PSTG3D WHICH SOLVES THE LINEAR C SYSTEM OF EQUATIONS. C C PORTABILITY AMERICAN NATIONAL STANDARDS INSTITUTE FORTRAN 77. C ALL MACHINE DEPENDENT CONSTANTS ARE LOCATED IN THE C FUNCTION PIMACH. C C REQUIRED COS,SIN C RESIDENT C ROUTINES C C * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * FINITE DIFFERENCE APPROXIMATIONS * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * C C AT INTERIOR GRIDPOINTS (X(I),Y(J),Z(K)) THE SECOND DERIVATIVES IN C THE POISSON EQUATION ARE REPLACED BY STANDARD SECOND-ORDER CENTRAL C DIFFERENCE APPROXIMATIONS TO YIELD THE LINEAR EQUATIONS C C (U(I-1,J,K) - 2*U(I,J,K) + U(I+1,J,K))/(DX**2) C + (U(I,J-1,K) - 2*U(I,J,K) + U(I,J+1,K))/(DY**2) C + (U(I,J,K-1) - 2*U(I,J,K) + U(I,J,K+1))/(DZ**2) C + ELMBDA*U(I,J,K) = F(I,J,K) C C FOR I=2,...,L-1, J=2,...,M-1, AND K=2,...,N-1. NEAR THE C BOUNDARIES THE APPROXIMATIONS TO THE BOUNDARY CONDITIONS ARE USED C TO ELIMINATE FICTIOUS POINTS. THE APPROXIMATIONS TO THE C BOUNDARY CONDITIONS ARE (FOR, SAY, THE X-BOUNDARIES): C C 1) IF U(XS,Y(J),Z(K)) IS SPECIFIED, THEN BY AVERAGING C C (U(0,J,K) + U(1,J,K))/2 = U(XS,Y(J),Z(K)). C C 2) IF (D/DX)U(XS,Y(J),Z(K)) IS SPECIFIED, THEN BY A CENTRAL C DIFFERENCE APPROXIMATION C C (U(1,J,K) - U(0,J,K))/DX = (D/DX)U(XS,Y(J),Z(K)). C C C