Test problems from Bakker

[ivan-pi]

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)

[ivan-pi]

An unverified transcript can be found in the details box:

Details

      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

[ivan-pi]

For Problem C of Bakker, the current code provides a good match:

[ivan-pi]

Same for Problem B: