SUBROUTINE SEPX4 (IORDER,A,B,M,MBDCND,BDA,ALPHA,BDB,BETA,C,D,N, 1 NBDCND,BDC,BDD,COFX,GRHS,USOL,IDMN,W,PERTRB, 2 IERROR) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * F I S H P A C K * C * * C * * C * A PACKAGE OF FORTRAN SUBPROGRAMS FOR THE SOLUTION OF * C * * C * SEPARABLE ELLIPTIC PARTIAL DIFFERENTIAL EQUATIONS * C * * C * (VERSION 4.0 , JUNE 1999) * C * * C * BY * C * * C * JOHN ADAMS, PAUL SWARZTRAUBER AND ROLAND SWEET * C * * C * OF * C * * C * THE NATIONAL CENTER FOR ATMOSPHERIC RESEARCH * C * * C * BOULDER, COLORADO (80307) U.S.A. * C * * C * WHICH IS SPONSORED BY * C * * C * THE NATIONAL SCIENCE FOUNDATION * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C C C DIMENSION OF BDA(N+1), BDB(N+1), BDC(M+1), BDD(M+1), C ARGUMENTS USOL(IDMN,N+1), GRHS(IDMN,N+1), C W (SEE ARGUMENT LIST) C C LATEST REVISION NOVEMBER 1988 C C PURPOSE SEPX4 SOLVES FOR EITHER THE SECOND-ORDER C FINITE DIFFERENCE APPROXIMATION OR A C FOURTH-ORDER APPROXIMATION TO A SEPARABLE C ELLIPTIC EQUATION C C AF(X)*UXX+BF(X)*UX+CF(X)*U+UYY = G(X,Y) C C ON A RECTANGLE (X GREATER THAN OR EQUAL TO C A AND LESS THAN OR EQUAL TO B, Y GREATER THAN C OR EQUAL TO C AND LESS THAN OR EQUAL TO D). C ANY COMBINATION OF PERIODIC OR MIXED BOUNDARY C CONDITIONS IS ALLOWED. IF BOUNDARY C CONDITIONS IN THE X DIRECTION ARE PERIODIC C (SEE MBDCND=0 BELOW) THEN THE COEFFICIENTS C MUST SATISFY C C AF(X)=C1,BF(X)=0,CF(X)=C2 FOR ALL X. C C HERE C1,C2 ARE CONSTANTS, C1.GT.0. C C THE POSSIBLE BOUNDARY CONDITIONS ARE: C IN THE X-DIRECTION: C (0) PERIODIC, U(X+B-A,Y)=U(X,Y) FOR C ALL Y,X C (1) U(A,Y), U(B,Y) ARE SPECIFIED FOR ALL Y C (2) U(A,Y), DU(B,Y)/DX+BETA*U(B,Y) ARE C SPECIFIED FOR ALL Y C (3) DU(A,Y)/DX+ALPHA*U(A,Y),DU(B,Y)/DX+ C BETA*U(B,Y) ARE SPECIFIED FOR ALL Y C (4) DU(A,Y)/DX+ALPHA*U(A,Y),U(B,Y) ARE C SPECIFIED FOR ALL Y C C IN THE Y-DIRECTION: C (0) PERIODIC, U(X,Y+D-C)=U(X,Y) FOR ALL X,Y C (1) U(X,C),U(X,D) ARE SPECIFIED FOR ALL X C (2) U(X,C),DU(X,D)/DY ARE SPECIFIED FOR C ALL X C (3) DU(X,C)/DY,DU(X,D)/DY ARE SPECIFIED FOR C ALL X C (4) DU(X,C)/DY,U(X,D) ARE SPECIFIED FOR C ALL X C C USAGE CALL SEPX4(IORDER,A,B,M,MBDCND,BDA,ALPHA,BDB, C BETA,C,D,N,NBDCND,BDC,BDD,COFX, C GRHS,USOL,IDMN,W,PERTRB,IERROR) C C ARGUMENTS C ON INPUT IORDER C = 2 IF A SECOND-ORDER APPROXIMATION IS C SOUGHT C = 4 IF A FOURTH-ORDER APPROXIMATION IS C SOUGHT C C A,B C THE RANGE OF THE X-INDEPENDENT VARIABLE, C I.E., X IS GREATER THAN OR EQUAL TO A C AND LESS THAN OR EQUAL TO B. A MUST BE C LESS THAN B. C C M C THE NUMBER OF PANELS INTO WHICH THE C INTERVAL (A,B) IS SUBDIVIDED. HENCE, C THERE WILL BE M+1 GRID POINTS IN THE X- C DIRECTION GIVEN BY XI=A+(I-1)*DLX C FOR I=1,2,...,M+1 WHERE DLX=(B-A)/M IS C THE PANEL WIDTH. M MUST BE LESS THAN C IDMN AND GREATER THAN 5. C C MBDCND C INDICATES THE TYPE OF BOUNDARY CONDITION C AT X=A AND X=B C = 0 IF THE SOLUTION IS PERIODIC IN X, I.E., C U(X+B-A,Y)=U(X,Y) FOR ALL Y,X C = 1 IF THE SOLUTION IS SPECIFIED AT X=A C AND X=B, I.E., U(A,Y) AND U(B,Y) ARE C SPECIFIED FOR ALL Y C = 2 IF THE SOLUTION IS SPECIFIED AT X=A C AND THE BOUNDARY CONDITION IS MIXED AT C X=B, I.E., U(A,Y) AND C DU(B,Y)/DX+BETA*U(B,Y) ARE SPECIFIED C FOR ALL Y C = 3 IF THE BOUNDARY CONDITIONS AT X=A AND C X=B ARE MIXED, I.E., C DU(A,Y)/DX+ALPHA*U(A,Y) AND C DU(B,Y)/DX+BETA*U(B,Y) ARE SPECIFIED C FOR ALL Y C = 4 IF THE BOUNDARY CONDITION AT X=A IS C MIXED AND THE SOLUTION IS SPECIFIED C AT X=B, I.E., DU(A,Y)/DX+ALPHA*U(A,Y) C AND U(B,Y) ARE SPECIFIED FOR ALL Y C C BDA C A ONE-DIMENSIONAL ARRAY OF LENGTH N+1 THAT C SPECIFIES THE VALUES OF C DU(A,Y)/DX+ ALPHA*U(A,Y) AT X=A, WHEN C MBDCND=3 OR 4. C BDA(J) = DU(A,YJ)/DX+ALPHA*U(A,YJ), C J=1,2,...,N+1 C WHEN MBDCND HAS ANY OTHER VALUE, BDA IS C A DUMMY PARAMETER. C C ALPHA C THE SCALAR MULTIPLYING THE SOLUTION IN CASE C OF A MIXED BOUNDARY CONDITION AT X=A C (SEE ARGUMENT BDA). IF MBDCND IS NOT EQUAL C TO EITHER 3 OR 4, THEN ALPHA IS A DUMMY C PARAMETER. C C BDB C A ONE-DIMENSIONAL ARRAY OF LENGTH N+1 THAT C SPECIFIES THE VALUES OF C DU(B,Y)/DX+ BETA*U(B,Y) AT X=B. C WHEN MBDCND=2 OR 3 C BDB(J) = DU(B,YJ)/DX+BETA*U(B,YJ), C J=1,2,...,N+1 C WHEN MBDCND HAS ANY OTHER VALUE, BDB IS C A DUMMY PARAMETER. C C BETA C THE SCALAR MULTIPLYING THE SOLUTION IN C CASE OF A MIXED BOUNDARY CONDITION AT X=B C (SEE ARGUMENT BDB). IF MBDCND IS NOT EQUAL C TO 2 OR 3, THEN BETA IS A DUMMY PARAMETER. C C C,D C THE RANGE OF THE Y-INDEPENDENT VARIABLE, C I.E., Y IS GREATER THAN OR EQUAL TO C AND C LESS THAN OR EQUAL TO D. C MUST BE LESS C THAN D. C C N C THE NUMBER OF PANELS INTO WHICH THE C INTERVAL (C,D) IS SUBDIVIDED. HENCE, C THERE WILL BE N+1 GRID POINTS IN THE Y- C DIRECTION GIVEN BY YJ=C+(J-1)*DLY FOR C J=1,2,...,N+1 WHERE DLY=(D-C)/N IS THE C PANEL WIDTH. IN ADDITION, N MUST BE C GREATER THAN 4. C C NBDCND C INDICATES THE TYPES OF BOUNDARY CONDITIONS C AT Y=C AND Y=D C = 0 IF THE SOLUTION IS PERIODIC IN Y, C I.E., U(X,Y+D-C)=U(X,Y) FOR ALL X,Y C = 1 IF THE SOLUTION IS SPECIFIED AT Y=C C AND Y = D, I.E., U(X,C) AND U(X,D) C ARE SPECIFIED FOR ALL X C = 2 IF THE SOLUTION IS SPECIFIED AT Y=C C AND THE BOUNDARY CONDITION IS MIXED C AT Y=D, I.E., DU(X,C)/DY AND U(X,D) C ARE SPECIFIED FOR ALL X C = 3 IF THE BOUNDARY CONDITIONS ARE MIXED C AT Y=CAND Y=D I.E., C DU(X,D)/DY AND DU(X,D)/DY ARE C SPECIFIED FOR ALL X C = 4 IF THE BOUNDARY CONDITION IS MIXED C AT Y=C AND THE SOLUTION IS SPECIFIED C AT Y=D, I.E. DU(X,C)/DY+GAMA*U(X,C) C AND U(X,D) ARE SPECIFIED FOR ALL X C C BDC C A ONE-DIMENSIONAL ARRAY OF LENGTH M+1 THAT C SPECIFIES THE VALUE DU(X,C)/DY AT Y=C. C C WHEN NBDCND=3 OR 4 C BDC(I) = DU(XI,C)/DY I=1,2,...,M+1. C C WHEN NBDCND HAS ANY OTHER VALUE, BDC IS C A DUMMY PARAMETER. C C BDD C A ONE-DIMENSIONAL ARRAY OF LENGTH M+1 THAT C SPECIFIED THE VALUE OF DU(X,D)/DY AT Y=D. C C WHEN NBDCND=2 OR 3 C BDD(I)=DU(XI,D)/DY I=1,2,...,M+1. C C WHEN NBDCND HAS ANY OTHER VALUE, BDD IS C A DUMMY PARAMETER. C C COFX C A USER-SUPPLIED SUBPROGRAM WITH PARAMETERS C X, AFUN, BFUN, CFUN WHICH RETURNS THE C VALUES OF THE X-DEPENDENT COEFFICIENTS C AF(X), BF(X), CF(X) IN THE ELLIPTIC C EQUATION AT X. IF BOUNDARY CONDITIONS IN C THE X DIRECTION ARE PERIODIC THEN THE C COEFFICIENTS MUST SATISFY AF(X)=C1,BF(X)=0, C CF(X)=C2 FOR ALL X. HERE C1.GT.0 C AND C2 ARE CONSTANTS. C C NOTE THAT COFX MUST BE DECLARED EXTERNAL C IN THE CALLING ROUTINE. C C GRHS C A TWO-DIMENSIONAL ARRAY THAT SPECIFIES THE C VALUES OF THE RIGHT-HAND SIDE OF THE C ELLIPTIC EQUATION, I.E.,GRHS(I,J)=G(XI,YI), C FOR I=2,...,M, J=2,...,N. AT THE C BOUNDARIES, GRHS IS DEFINED BY C C MBDCND GRHS(1,J) GRHS(M+1,J) C ------ --------- ----------- C 0 G(A,YJ) G(B,YJ) C 1 * * C 2 * G(B,YJ) J=1,2,...,N+1 C 3 G(A,YJ) G(B,YJ) C 4 G(A,YJ) * C C NBDCND GRHS(I,1) GRHS(I,N+1) C ------ --------- ----------- C 0 G(XI,C) G(XI,D) C 1 * * C 2 * G(XI,D) I=1,2,...,M+1 C 3 G(XI,C) G(XI,D) C 4 G(XI,C) * C C WHERE * MEANS THESE QUANTITES ARE NOT USED. C GRHS SHOULD BE DIMENSIONED IDMN BY AT LEAST C N+1 IN THE CALLING ROUTINE. C C USOL C A TWO-DIMENSIONAL ARRAY THAT SPECIFIES THE C VALUES OF THE SOLUTION ALONG THE BOUNDARIES. C AT THE BOUNDARIES, USOL IS DEFINED BY C C MBDCND USOL(1,J) USOL(M+1,J) C ------ --------- ----------- C 0 * * C 1 U(A,YJ) U(B,YJ) C 2 U(A,YJ) * J=1,2,...,N+1 C 3 * * C 4 * U(B,YJ) C C NBDCND USOL(I,1) USOL(I,N+1) C ------ --------- ----------- C 0 * * C 1 U(XI,C) U(XI,D) C 2 U(XI,C) * I=1,2,...,M+1 C 3 * * C 4 * U(XI,D) C C WHERE * MEANS THE QUANTITES ARE NOT USED C IN THE SOLUTION. C C IF IORDER=2, THE USER MAY EQUIVALENCE GRHS C AND USOL TO SAVE SPACE. NOTE THAT IN THIS C CASE THE TABLES SPECIFYING THE BOUNDARIES C OF THE GRHS AND USOL ARRAYS DETERMINE THE C BOUNDARIES UNIQUELY EXCEPT AT THE CORNERS. C IF THE TABLES CALL FOR BOTH G(X,Y) AND C U(X,Y) AT A CORNER THEN THE SOLUTION MUST C BE CHOSEN. C FOR EXAMPLE, IF MBDCND=2 AND NBDCND=4, C THEN U(A,C), U(A,D),U(B,D) MUST BE CHOSEN C AT THE CORNERS IN ADDITION TO G(B,C). C C IF IORDER=4, THEN THE TWO ARRAYS, USOL AND C GRHS, MUST BE DISTINCT. C C USOL SHOULD BE DIMENSIONED IDMN BY AT LEAST C N+1 IN THE CALLING ROUTINE. C C IDMN C THE ROW (OR FIRST) DIMENSION OF THE ARRAYS C GRHS AND USOL AS IT APPEARS IN THE PROGRAM C CALLING SEPELI. THIS PARAMETER IS USED C TO SPECIFY THE VARIABLE DIMENSION OF GRHS C AND USOL. IDMN MUST BE AT LEAST 7 AND C GREATER THAN OR EQUAL TO M+1. C C W C A ONE-DIMENSIONAL ARRAY THAT MUST BE C PROVIDED BY THE USER FOR WORK SPACE. C 10*N+(16+INT(LOG2(N+1)))*(M+1)+11 WILL C SUFFICE AS A LENGTH FOR W. THE ACTUAL C LENGTH OF W IN THE CALLING ROUTINE C MUST BE SET IN W(1) (SEE IERROR=11). C C C ON OUTPUT USOL C CONTAINS THE APPROXIMATE SOLUTION TO THE C ELLIPTIC EQUATION. USOL(I,J) IS THE C APPROXIMATION TO U(XI,YJ) FOR I=1,2...,M+1 C AND J=1,2,...,N+1. THE APPROXIMATION HAS C ERROR O(DLX**2+DLY**2) IF CALLED WITH C IORDER=2 AND O(DLX**4+DLY**4) IF CALLED C WITH IORDER=4. C C W C CONTAINS INTERMEDIATE VALUES THAT MUST NOT C BE DESTROYED IF SEPELI IS CALLED AGAIN C WITH INTL=1. IN ADDITION W(1) CONTAINS C THE EXACT MINIMAL LENGTH (IN FLOATING POINT) C REQUIRED FOR THE WORK SPACE (SEE IERROR=11). C C PERTRB C IF A COMBINATION OF PERIODIC OR DERIVATIVE C BOUNDARY CONDITIONS (I.E., ALPHA=BETA=0 IF C MBDCND=3) IS SPECIFIED AND IF CF(X)=0 FOR C ALL X THEN A SOLUTION TO THE DISCRETIZED C MATRIX EQUATION MAY NOT EXIST C (REFLECTING THE NON-UNIQUENESS OF SOLUTIONS C TO THE PDE). C PERTRB IS A CONSTANT CALCULATED AND C SUBTRACTED FROM THE RIGHT HAND SIDE OF THE C MATRIX EQUATION INSURING THE EXISTENCE OF A C SOLUTION. SEPX4 COMPUTES THIS SOLUTION C WHICH IS A WEIGHTED MINIMAL LEAST SQUARES C SOLUTION TO THE ORIGINAL PROBLEM. IF C SINGULARITY IS NOT DETECTED PERTRB=0.0 IS C RETURNED BY SEPX4. C C IERROR C AN ERROR FLAG THAT INDICATES INVALID INPUT C PARAMETERS OR FAILURE TO FIND A SOLUTION C C = 0 NO ERROR C = 1 IF A GREATER THAN B OR C GREATER C THAN D C = 2 IF MBDCND LESS THAN 0 OR MBDCND C GREATER THAN 4 C = 3 IF NBDCND LESS THAN 0 OR NBDCND C GREATER THAN 4 C = 4 IF ATTEMPT TO FIND A SOLUTION FAILS. C (THE LINEAR SYSTEM GENERATED IS NOT C DIAGONALLY DOMINANT.) C = 5 IF IDMN IS TOO SMALL (SEE DISCUSSION C OF IDMN) C = 6 IF M IS TOO SMALL OR TOO LARGE C (SEE DISCUSSION OF M) C = 7 IF N IS TOO SMALL (SEE DISCUSSION OF N) C = 8 IF IORDER IS NOT 2 OR 4 C = 9 IF INTL IS NOT 0 OR 1 C = 10 IF AFUN IS LESS THAN OR EQUAL TO ZERO C FOR SOME INTERIOR MESH POINT XI SOME C INTERIOR MESH POINT (XI,YJ) C = 11 IF THE WORK SPACE LENGTH INPUT IN W(1) C IS LESS THAN THE EXACT MINIMAL WORK C SPACE LENGTH REQUIRED OUTPUT IN W(1). C = 12 IF MBDCND=0 AND AF(X)=CF(X)=CONSTANT C OR BF(X)=0 FOR ALL X IS NOT TRUE. C C SPECIAL CONDITIONS NONE C C I/O NONE C C REQUIRED LIBRARY COMF, GENBUN, GNBNAUX, AND SEPAUX C FILES FROM FISHPAK C C C PRECISION SINGLE C C REQUIRED LIBRARY NONE C FILES C C LANGUAGE FORTRAN C C HISTORY SEPX4 WAS DEVELOPED AT NCAR BY JOHN C. C ADAMS OF THE SCIENTIFIC COMPUTING DIVISION C IN OCTOBER 1978. THE BASIS OF THIS CODE IS C NCAR ROUTINE SEPELI. BOTH PACKAGES WERE C RELEASED ON NCAR'S PUBLIC LIBRARIES IN C JANUARY 1980. C C PORTABILITY FORTRAN 77 C C ALGORITHM SEPX4 AUTOMATICALLY DISCRETIZES THE SEPARABLE C ELLIPTIC EQUATION WHICH IS THEN SOLVED BY A C GENERALIZED CYCLIC REDUCTION ALGORITHM IN THE C SUBROUTINE POIS. THE FOURTH ORDER SOLUTION C IS OBTAINED USING THE TECHNIQUE OF DEFFERRED C CORRECTIONS REFERENCED BELOW. C C TIMING WHEN POSSIBLE, SEPX4 SHOULD BE USED INSTEAD C OF PACKAGE SEPELI. THE INCREASE IN SPEED C IS AT LEAST A FACTOR OF THREE. C C REFERENCES KELLER, H.B., NUMERICAL METHODS FOR TWO-POINT C BOUNDARY-VALUE PROBLEMS, BLAISDEL (1968), C WALTHAM, MASS. C C SWARZTRAUBER, P., AND R. SWEET (1975): C EFFICIENT FORTRAN SUBPROGRAMS FOR THE C SOLUTION OF ELLIPTIC PARTIAL DIFFERENTIAL C EQUATIONS. NCAR TECHNICAL NOTE C NCAR-TN/IA-109, PP. 135-137. C***********************************************************************