SUBROUTINE SEPELI (INTL,IORDER,A,B,M,MBDCND,BDA,ALPHA,BDB,BETA,C, 1 D,N,NBDCND,BDC,GAMA,BDD,XNU,COFX,COFY,GRHS, 2 USOL,IDMN,W,PERTRB,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 SEPELI SOLVES FOR EITHER THE SECOND-ORDER C FINITE DIFFERENCE APPROXIMATION OR A C FOURTH-ORDER APPROXIMATION TO A SEPARABLE C ELLIPTIC EQUATION C C 2 2 C AF(X)*D U/DX + BF(X)*DU/DX + CF(X)*U + C 2 2 C DF(Y)*D U/DY + EF(Y)*DU/DY + FF(Y)*U C C = G(X,Y) C C ON A RECTANGLE (X GREATER THAN OR EQUAL TO A C 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. C C THE POSSIBLE BOUNDARY CONDITIONS ARE: C IN THE X-DIRECTION: C (0) PERIODIC, U(X+B-A,Y)=U(X,Y) FOR ALL C Y,X (1) U(A,Y), U(B,Y) ARE SPECIFIED FOR C 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+XNU*U(X,D) ARE C SPECIFIED FOR ALL X C (3) DU(X,C)/DY+GAMA*U(X,C),DU(X,D)/DY+ C XNU*U(X,D) ARE SPECIFIED FOR ALL X C (4) DU(X,C)/DY+GAMA*U(X,C),U(X,D) ARE C SPECIFIED FOR ALL X C C USAGE CALL SEPELI (INTL,IORDER,A,B,M,MBDCND,BDA, C ALPHA,BDB,BETA,C,D,N,NBDCND,BDC, C GAMA,BDD,XNU,COFX,COFY,GRHS,USOL, C IDMN,W,PERTRB,IERROR) C C ARGUMENTS C ON INPUT INTL C = 0 ON INITIAL ENTRY TO SEPELI OR IF ANY C OF THE ARGUMENTS C,D, N, NBDCND, COFY C ARE CHANGED FROM A PREVIOUS CALL C = 1 IF C, D, N, NBDCND, COFY ARE UNCHANGED C FROM THE PREVIOUS CALL. C C IORDER C = 2 IF A SECOND-ORDER APPROXIMATION C IS SOUGHT C = 4 IF A FOURTH-ORDER APPROXIMATION C IS 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 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 AND C THE BOUNDARY CONDITION IS MIXED AT X=B, C I.E., U(A,Y) AND DU(B,Y)/DX+BETA*U(B,Y) C ARE SPECIFIED 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 C THAT 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. WHEN MBDCND HAS ANY OTHER C OTHER VALUE, BDA IS A DUMMY PARAMETER. C C ALPHA C THE SCALAR MULTIPLYING THE SOLUTION IN C CASE OF A MIXED BOUNDARY CONDITION AT X=A C (SEE ARGUMENT BDA). IF MBDCND IS NOT C EQUAL TO 3 OR 4 THEN ALPHA IS A DUMMY C PARAMETER. C C BDB C A ONE-DIMENSIONAL ARRAY OF LENGTH N+1 C THAT 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. WHEN MBDCND HAS ANY OTHER C OTHER VALUE, BDB IS A DUMMY PARAMETER. C C BETA C THE SCALAR MULTIPLYING THE SOLUTION IN C CASE OF A MIXED BOUNDARY CONDITION AT C X=B (SEE ARGUMENT BDB). IF MBDCND IS C NOT EQUAL TO 2 OR 3 THEN BETA IS A DUMMY C PARAMETER. C C C,D C THE RANGE OF THE Y-INDEPENDENT VARIABLE, C I.E., Y IS GREATER THAN OR EQUAL TO C C AND LESS THAN OR EQUAL TO D. C MUST BE C LESS THAN D. C C N C THE NUMBER OF PANELS INTO WHICH THE C INTERVAL [C,D] IS SUBDIVIDED. C HENCE, THERE WILL BE N+1 GRID POINTS C IN THE Y-DIRECTION GIVEN BY C YJ=C+(J-1)*DLY FOR J=1,2,...,N+1 WHERE C DLY=(D-C)/N IS THE PANEL WIDTH. C IN ADDITION, N MUST BE GREATER THAN 4. C C NBDCND C INDICATES THE TYPES OF BOUNDARY CONDITIONS C AT Y=C AND Y=D C 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., U(X,C) AND C DU(X,D)/DY+XNU*U(X,D) ARE SPECIFIED C FOR ALL X C = 3 IF THE BOUNDARY CONDITIONS ARE MIXED C AT Y=C AND Y=D, I.E., C DU(X,D)/DY+GAMA*U(X,C) AND C DU(X,D)/DY+XNU*U(X,D) ARE SPECIFIED C 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 C THAT SPECIFIES THE VALUE OF C DU(X,C)/DY+GAMA*U(X,C) AT Y=C. C WHEN NBDCND=3 OR 4 BDC(I) = DU(XI,C)/DY + C GAMA*U(XI,C), I=1,2,...,M+1. C WHEN NBDCND HAS ANY OTHER VALUE, BDC C IS A DUMMY PARAMETER. C C GAMA C THE SCALAR MULTIPLYING THE SOLUTION IN C CASE OF A MIXED BOUNDARY CONDITION AT C Y=C (SEE ARGUMENT BDC). IF NBDCND IS C NOT EQUAL TO 3 OR 4 THEN GAMA IS A DUMMY C PARAMETER. C C BDD C A ONE-DIMENSIONAL ARRAY OF LENGTH M+1 C THAT SPECIFIES THE VALUE OF C DU(X,D)/DY + XNU*U(X,D) AT Y=C. C WHEN NBDCND=2 OR 3 BDD(I) = DU(XI,D)/DY + C XNU*U(XI,D), I=1,2,...,M+1. C WHEN NBDCND HAS ANY OTHER VALUE, BDD C IS A DUMMY PARAMETER. C C XNU C THE SCALAR MULTIPLYING THE SOLUTION IN C CASE OF A MIXED BOUNDARY CONDITION AT C Y=D (SEE ARGUMENT BDD). IF NBDCND IS C NOT EQUAL TO 2 OR 3 THEN XNU IS A C DUMMY PARAMETER. C C COFX C A USER-SUPPLIED SUBPROGRAM WITH C PARAMETERS X, AFUN, BFUN, CFUN WHICH C RETURNS THE VALUES OF THE X-DEPENDENT C COEFFICIENTS AF(X), BF(X), CF(X) IN THE C ELLIPTIC EQUATION AT X. C C COFY C A USER-SUPPLIED SUBPROGRAM WITH PARAMETERS C Y, DFUN, EFUN, FFUN WHICH RETURNS THE C VALUES OF THE Y-DEPENDENT COEFFICIENTS C DF(Y), EF(Y), FF(Y) IN THE ELLIPTIC C EQUATION AT Y. C C NOTE: COFX AND COFY MUST BE DECLARED C EXTERNAL IN THE CALLING ROUTINE. C THE VALUES RETURNED IN AFUN AND DFUN C MUST SATISFY AFUN*DFUN GREATER THAN 0 C FOR A LESS THAN X LESS THAN B, C LESS C THAN Y LESS THAN D (SEE IERROR=10). C THE COEFFICIENTS PROVIDED MAY LEAD TO A C MATRIX EQUATION WHICH IS NOT DIAGONALLY C DOMINANT IN WHICH CASE SOLUTION MAY FAIL C (SEE IERROR=4). 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., C GRHS(I,J)=G(XI,YI), FOR I=2,...,M, C J=2,...,N. AT THE BOUNDARIES, GRHS IS C 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 QUANTITIES 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 QUANTITIES 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. FOR EXAMPLE, IF MBDCND=2 AND C NBDCND=4, THEN U(A,C), U(A,D), U(B,D) MUST C BE CHOSEN AT THE CORNERS IN ADDITION C 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 LET K=INT(LOG2(N+1))+1 AND SET L=2**(K+1). C THEN (K-2)*L+K+10*N+12*M+27 WILL SUFFICE C AS A LENGTH OF W. THE ACTUAL LENGTH OF W C IN THE CALLING ROUTINE MUST BE SET IN W(1) C (SEE IERROR=11). C C ON OUTPUT USOL C CONTAINS THE APPROXIMATE SOLUTION TO THE C ELLIPTIC EQUATION. C USOL(I,J) IS THE APPROXIMATION TO U(XI,YJ) C FOR I=1,2...,M+1 AND J=1,2,...,N+1. C THE APPROXIMATION HAS ERROR C O(DLX**2+DLY**2) IF CALLED WITH IORDER=2 C AND O(DLX**4+DLY**4) IF CALLED WITH C IORDER=4. C C W C CONTAINS INTERMEDIATE VALUES THAT MUST NOT C BE DESTROYED IF SEPELI IS CALLED AGAIN WITH C INTL=1. IN ADDITION W(1) CONTAINS THE C 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 C (I.E., ALPHA=BETA=0 IF MBDCND=3; C GAMA=XNU=0 IF NBDCND=3) IS SPECIFIED C AND IF THE COEFFICIENTS OF U(X,Y) IN THE C SEPARABLE ELLIPTIC EQUATION ARE ZERO C (I.E., CF(X)=0 FOR X GREATER THAN OR EQUAL C TO A AND LESS THAN OR EQUAL TO B; C FF(Y)=0 FOR Y GREATER THAN OR EQUAL TO C C AND LESS THAN OR EQUAL TO D) THEN A C SOLUTION MAY NOT EXIST. PERTRB IS A C CONSTANT CALCULATED AND SUBTRACTED FROM C THE RIGHT-HAND SIDE OF THE MATRIX EQUATIONS C GENERATED BY SEPELI WHICH INSURES THAT A C SOLUTION EXISTS. SEPELI THEN COMPUTES THIS C SOLUTION WHICH IS A WEIGHTED MINIMAL LEAST C SQUARES SOLUTION TO THE ORIGINAL PROBLEM. C C IERROR C AN ERROR FLAG THAT INDICATES INVALID INPUT C PARAMETERS OR FAILURE TO FIND A SOLUTION C = 0 NO ERROR C = 1 IF A GREATER THAN B OR C GREATER THAN D C = 2 IF MBDCND LESS THAN 0 OR MBDCND GREATER C THAN 4 C = 3 IF NBDCND LESS THAN 0 OR NBDCND GREATER C 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 C (SEE DISCUSSION 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*DFUN LESS THAN OR EQUAL TO 0 C FOR SOME 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 C NOTE (CONCERNING IERROR=4): FOR THE C COEFFICIENTS INPUT THROUGH COFX, COFY, C THE DISCRETIZATION MAY LEAD TO A BLOCK C TRIDIAGONAL LINEAR SYSTEM WHICH IS NOT C DIAGONALLY DOMINANT (FOR EXAMPLE, THIS C HAPPENS IF CFUN=0 AND BFUN/(2.*DLX) GREATER C THAN AFUN/DLX**2). IN THIS CASE SOLUTION C MAY FAIL. THIS CANNOT HAPPEN IN THE LIMIT C AS DLX, DLY APPROACH ZERO. HENCE, THE C CONDITION MAY BE REMEDIED BY TAKING LARGER C VALUES FOR M OR N. C C SPECIAL CONDITIONS SEE COFX, COFY ARGUMENT DESCRIPTIONS ABOVE. C C I/O NONE C C PRECISION SINGLE C C REQUIRED LIBRARY BLKTRI, COMF, AND SEPAUX C FILES FROM FISHPAK C C LANGUAGE FORTRAN C C HISTORY DEVELOPED AT NCAR DURING 1975-76 BY C JOHN C. ADAMS OF THE SCIENTIFIC COMPUTING C DIVISION. RELEASED ON NCAR'S PUBLIC SOFTWARE C LIBRARIES IN JANUARY 1980. C C PORTABILITY FORTRAN 77 C C ALGORITHM SEPELI AUTOMATICALLY DISCRETIZES THE C SEPARABLE ELLIPTIC EQUATION WHICH IS THEN C SOLVED BY A GENERALIZED CYCLIC REDUCTION C ALGORITHM IN THE SUBROUTINE, BLKTRI. THE C FOURTH-ORDER SOLUTION IS OBTAINED USING C 'DEFERRED CORRECTIONS' WHICH IS DESCRIBED C AND REFERENCED IN SECTIONS, REFERENCES AND C METHOD. C C TIMING THE OPERATIONAL COUNT IS PROPORTIONAL TO C M*N*LOG2(N). C C ACCURACY THE FOLLOWING ACCURACY RESULTS WERE OBTAINED C ON A CDC 7600. NOTE THAT THE FOURTH-ORDER C ACCURACY IS NOT REALIZED UNTIL THE MESH IS C SUFFICIENTLY REFINED. C C SECOND-ORDER FOURTH-ORDER C M N ERROR ERROR C C 6 6 6.8E-1 1.2E0 C 14 14 1.4E-1 1.8E-1 C 30 30 3.2E-2 9.7E-3 C 62 62 7.5E-3 3.0E-4 C 126 126 1.8E-3 3.5E-6 C 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***********************************************************************