ivan-pi
commented
1 year ago An unverified transcript can be found in the details box:
```fortran
SUBROUTINE PDEF1(X,T,U,NPDE,NPTS,NC,FEVAL,GEVAL,BNDRY,
* ALFA,BETA,GAMMA,UMEAN,UX,F,G)
DIMENSION X(NPTS),U(NPDE,NPTS),ALFA(NPDE),BETA(NPDE),
* GAMMA(NPDE),UX(NPDE),UMEAN(NPDE),F(NPDE),G(NPDE)
C
C THIS SUBROUTINE SERVES AS AN INTERFACE BETWEEN A
C PARABOLID PARTIAL DIFFERENTIAL EQUATION IN ONE
C SPACE VARIABLE AND AN INTEGRATOR OF INITIAL VALUE PROBLEMS.
C
C THE DIFFERNTIAL EUQATION SHOULD BE OF THE FROM
C
C (D/DT)U(I) = (D/DX)(X**NC*F(I,X,T,U,UX)))/X**NC
C + G(I,X,T,U,UX),I=1,...,NPDE;
C
C
C WITH BOUNDARY CONDITIONS
C
C
C ALFA(I)*U(I) + BETA(I)*UX(I) = GAMMA(I) ,I=1,...,NPDE;
C
C
C PDEF1 TRANSFORMS THE SYSTEM OF PARTIAL DIFFERENTIAL EQUATIONS
C INTO A SYSTEM OF ORDINARY EQUATIONS BY SEMI-DISCRETIZATION IN
C THE SPACE VARIABLE; THIS SEMI-DISCRETIZATION IS PERFORMED
C BY APPLICATION OF THE FINITE ELEMENT METHOD (SEE E.G.
C )
C
C
C DESCRIPTION OF THE PARAMETERS
C
C INPUT:
C
C X: DIMENSION X(NPTS);
C X(1),...,X(NPTS) IS A PARTITION OF [X(1),X(NPTS)]
C
C T: THE TIME VARIABLE;
C
C U: DIMENSION U(NPDE,NPTS);
C U(I,J) IS AN APPROXIMATION OF U(I,X(J),T),I=1,...,NPDE;
C J=1,...,NPTS;
C
C
C NPDE: THE NUMBER OF PARTIAL DIFFERENTIAL EQUATIONS;
C
C NPTS: THE NUMBER OF GRIDPOINTS:
C
C NC: A NUMBER DESIGNING THE KIND OF SPACE COORDINATES;
C NX = 0 : CARTESIAN COORDINATES;
C NX = 1 : CYLINDRICAL COORDINATES;
C NX = 2 : SPHERICAL COORDINATES;
C
C FEVAL: SUBROUTINE FEVAL(X,T,U,UX,F,NPDE)
C DIMENSION U(NPDE),UX(NPDE),F(NPDE)
C
C EXIT: THE VALUES OF F(I,X,T,U,UX) ARE ASSIGNED
C TO THE ARRAY F;
C
C GEVAL: SUBROUTINE GEVAL(X,T,U,UX,G,NPDE)
C DIMENSTION U(NPDE),UX(NPDE),G(NPDE)
C
C EXIT: THE VALUES OF G(I,X,T,U,UX) ARE ASSIGNED
C TO THE ARRAY G;
C
C BNDRY: SUBROUTINE BNDRY(T,ALPHA,BETA,GAMMA,U,NPDE,LEFT)
C DIMENSION ALFA(NPDE),BETA(NPDE),GAMMA(NPDE),U(NPDE)
C LOGICAL LEFT
C
C THIS SUBROUTINE DEFINES THE BOUNDARY CONDITIONS OF THE PDE
C IN X(1) OR X(NPTS): IF LEFT = .TRUE.,THE BOUNDARY CONDI-
C TIONS IN X(1) ARE DEFINED, OTHERWISE THE B.C. IN X(NPTS);
C ALFA,BETA AND GAMMA MAY DEPEND ON U AND T.
C
C ALFA,BETA,GAMMA,UMEAN,UX,F,G:
C
C WORK-ARRAYS OF DIMENSION NPDE.
C
C OUTPUT:
C
C U: DIMENSION U(NPDE,NPTS);
C THE SEMI-DISCRETIZED RIGHT HAND SIDE OF THE PDE OVERWRITTEN
C ON U;
C
WR = 0.
XR = X(1)
DO 150 L = 2,NPTS
C
C THE SEGMENTWISE ASSEMBLY OF THE RIGHT HAND SIDE BEGINS;
C AT FIRST, SOME NUMBERS CHARACTERISTIC FOR NC ARE COMPUTED.
C
XL = XR
L1 = L - 1
XR = X(L)
H = XR - XL
IF (NC .EQ. 1) GOTO 10
IF (NC .EQ. 2) GOTO 20
VL = 0.5
VR = 0.5
GOTO 30
10 VL = XL/3. + XR/6.
VR = XL/6. + XR/3.
GOTO 30
20 XLSQ = XL*XL/12.
XLXR = XL*XR/6.
XRSQ = XR*XR/12.
VL = 3.*XLSQ + XLXR + XRSQ
VR = XLSQ + XLXR + 3.*XRSQ
30 WL = H*VL + WR
WR = H*VR
VMEAN = VL + VR
WMEAN = H*VMEAN
PR = VR/VMEAN
PL = 1.0 - PR
XMEAN = XL + H*PR
IF (L .GT. 2) GOTO 90
C
C ON THE FIRST SEGMENT,AT FIRST THE BOUNDARY CONDITIONS ARE
C PROCESSED;IF THEY ARE OF DIRICHLET TYPE,THE VALUE OF U(X(1),T)
C IS UPDATED AND THE TIME DERIVATIVE IS PUT EQUAL TO ZERO;
C IF THEY ARE OF MIXED TYPE,THE SPATIAL DERIVATIVE OF U(X(1),T)
C OF U(X(1),T) IS EXPRESSED IN U AND T AND THAT EXPRESSION IS
C IS IMPLEMENTED IN THE STOCK TERM
C
DO 40 J = 1,NPDE
UMEAN(J) = U(J,1)
U(J,1) = 0.
40 CONTINUE
CALL BNDRY(T,ALFA,BETA,GAMMA,UMEAN,NPDE,.TRUE.)
ICT = 0
DO 60 J = 1,NPDE
IF (BETA(J) .NE. 0.) GOTO 60
ICT = ICT + 1
UMEAN(J) = GAMMA(J)/ALFA(J)
60 CONTINUE
XPOW = 1.
IF (NC .GT. 0) XPOW = XL**NC
IF (XPOW .EQ. 0. .OR. ICT .EQ. NPDE) GOTO 90
IF (ICT .GT. 0) CALL BNDRY(T,ALFA,BETA,GAMMA,UMEAN,NPDE,.TRUE.)
DO 70 J = 1,NPDE
IF (BETA(J) .NE. 0.)
* UX(J) = (GAMMA(J) - ALFA(J)*UMEAN(J))/BETA(J)
IF (BETA(J) .EQ. 0.) UX(J) = (U(J,2) - UMEAN(J))/H
70 CONTINUE
CALL FEVAL(XL,T,UMEAN,UX,F,NPDE)
C
C NOW,THE STOCK TERMS ARE IMPLEMENTED,AS FAR AS THE
C BOUNDARY CONDITIONS ARE OD MIXED TYPE.
C
DO 80 J = 1,NPDE
80 IF (BETA(J) .NE. 0.) U(J,1) = -XPOW*F(J)
C
C HERE THE IMPLEMENTATION OF THE LEFT BOUNDARY CONDITIONS
C ENDS
C
90 IF (L .LT. NPTS) GOTO 120
C
C AT THE LAST SEGMENT,THE RIGHT BOUNDARY CONDITIONS,AS FAR
C AS THEY ARE OF DIRICHLET TYPE ARE IMPLEMENTED BY UPDATING
C THE BOUNDARY VALUES;
C
DO 100 J = 1,NPDE
IF (BETA(J) .EQ. 0.) U(J,1) = 0.
UX(J) = U(J,NPTS)
100 CONTINUE
CALL BNDRY(T,ALFA,BETA,GAMMA,UX,NPDE,.FALSE.)
ICT= 0
DO 110 J = 1,NPDE
IF (BETA(J) .NE. 0.) GOTO 110
ICT = ICT + 1
U(J,NPTS) = GAMMA(J)/ALFA(J)
110 CONTINUE
C
C NOW, THE REAL ASSEMBLY BEGINS
C
120 DO 130 J = 1,NPDE
ULJ = UMEAN(J)
URJ = U(J,L)
UMEAN(J) = PL*ULJ + PR*URJ
UX(J) = (URJ - ULJ)/H
130 CONTINUE
CALL FEVAL(XMEAN,T,UMEAN,UX,F,NPDE)
CALL GEVAL(XMEAN,T,UMEAN,UX,G,NPDE)
DO 140 J = 1,NPDE
FMEAN = VMEAN*F(J)
GMEAN = WMEAN*G(J)
U(J,L1) = (U(J,L1) + FMEAN + PL*GMEAN)/WL
UMEAN(J) = U(J,L)
U(J,L) = - FMEAN + PR*GMEAN
140 CONTINUE
150 CONTINUE
C
C FINALLY,THE PROCESSING OF THE RIGHT BOUNDARY CONDITIONS
C IS PERFORMED.
C
DO 160 J = 1,NPDE
160 IF (BETA(J) .EQ. 0.) U(J,NPTS) = 0.
IF (ICT .EQ. NPDE) RETURN
IF (ICT .GT. 0) CALL BNDRY(T,ALFA,BETA,GAMMA,UMEAN,NPDE,.FALSE.)
XPOW = XR**NC
DO 170 J = 1,NPDE
170 IF (BETA(J) .NE. 0.) UX(J) = (GAMMA(J)-ALFA(J)*UMEAN(J))/BETA(J)
CALL FEVAL(XR,T,UMEAN,UX,F,NPDE)
DO 180 J = 1,NPDE
180 IF (BETA(J) .NE. 0.) U(J,NPTS) = (U(J,NPTS) + XPOW*F(J))/WR
RETURN
END
```
Would be nice to add the test problems from the work of Bakker: https://ir.cwi.nl/pub/9016 (this would include also a comparison with the PDEF1/PDEF2 codes)