c c file sepx4.txt (documentation for the FISHPACK solver SEPX4) c c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . c . . c . copyright (c) 2004 by UCAR . c . . c . UNIVERSITY CORPORATION for ATMOSPHERIC RESEARCH . c . . c . all rights reserved . c . . c . . c . FISHPACK version 5.0 . c . . c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 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 5.0 , JUNE 2004) * 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 SUBROUTINE SEPX4 (IORDER,A,B,M,MBDCND,BDA,ALPHA,BDB,BETA,C,D,N, C + NBDCND,BDC,BDD,COFX,GRHS,USOL,IDMN,PERTRB, C + IERROR) 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 C C LATEST REVISION June 2004 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 *** caution *** GRHS SHOULD BE RESET IF SEPX4 WAS FIRST CALLED C WITH IORDER=2 AND WILL BE CALLED AGAIN WITH C IORDER=4. VALUES IN GRHS ARE DESTROYED BY THE C IORDER=2 CALL. C 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 *** caution GRHS SHOULD BE RESET IF SEPX4 WAS FIRST CALLED C WITH IORDER=2 AND WILL BE CALLED AGAIN WITH C IORDER=4. VALUES IN GRHS ARE DESTROYED BY THE C IORDER=2 CALL. 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 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 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 = 12 IF MBDCND=0 AND AF(X)=CF(X)=CONSTANT C OR BF(X)=0 FOR ALL X IS NOT TRUE. C = 20 If the dynamic allocation of real and C complex work space required for solution C fails (for example if N,M are too large C for your computer) C C SPECIAL CONDITIONS NONE C C I/O NONE C C REQUIRED files fish.f,comf.f,genbun.f,gnbnaux.f,sepaux.f C C C PRECISION SINGLE C C C LANGUAGE FORTRAN 90 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. SEPX4 was modified in June 2004 c incorporating fortran 90 dynamical storage c allocation for work space requirements C C PORTABILITY FORTRAN 90 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***********************************************************************