C ALGORITHM 758, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 22, NO. 3, September, 1996, P. 302--328. C #! /bin/sh # This is a shell archive, meaning: # 1. Remove everything above the #! /bin/sh line. # 2. Save the resulting text in a file. # 3. Execute the file with /bin/sh (not csh) to create the files: # Doc # Drivers # Info # Src # This archive created: Wed Sep 25 11:40:52 1996 export PATH; PATH=/bin:$PATH if test ! -d 'Doc' then mkdir 'Doc' fi cd 'Doc' if test -f 'readme' then echo shar: will not over-write existing file "'readme'" else cat << \SHAR_EOF > 'readme' `VLUGR2: A Vectorizable Adaptive Grid Solver for PDEs in 2D' by J.G. Blom, R.A. Trompert, and J.G. Verwer. This code solves systems of PDEs of the type F(t,x,y,U,Ut,Ux,Uy,Uxx,Uxy,Uyy)=0 with boundary conditions B(t,x,y,U,Ut,Ux,Uy)=0 and initial values U(t0,x,y)=U0 on a 2D domain bounded by right-angled polygons. In space Local Uniform Grid Refinement is applied to resolve local sharp gradients in the solution. For the time integration the implicit BDF2 method is used with variable stepsizes. Description of contents of source code files -------------------------------------------- (Both single precision and double precision available) src.f Main module, contains documentation ilubsn.f ILU decomposition and backsolve for arbitrary number of PDE components ilubs1.f ILU decomposition and backsolve for optimal vector performance for PDE with 1 component ilubs2.f ILU decomposition and backsolve for optimal vector performance for PDE with 2 components ilubs3.f ILU decomposition and backsolve for optimal vector performance for PDE with 3 components user.f Default modules that can be replaced by user's own (see description in paper) blas.f BLAS modules exmpl.f Calling program for the first time interval of the example in the paper exmplr.f Calling program for the second time interval of the example in the paper probi.f Calling program for problem I in the paper probii.f Calling program for problem II in the paper prtsol.f Program to print out solution from file generated by the DUMP routine wrtuni.f Program that reads the file generated by the DUMP routine and writes the (interpolated) solution on a uniform grid of a specified grid level and the maximum used grid level in each point to file. plot.m Matlab plotting routine to plot the data generated by WRTUNI.f How to use the solver: ---------------------- Compile and link the modules in exmpl.f user.f (only the SUBROUTINE DERIVF) src.f ilubSn.f blas.f (if the BLAS library is not available on the platform) The module blas.f contains, a.o., the functions I1MACH and R/D1MACH which set machine-dependent values. These functions need to be adapted to the platform. The results and integration information can be found in the files exmpl_runinfo exmpl_output A file DUMP is created that contains all the necessary information to restart the integration on the second time interval. For the second run one should compile and link the modules in exmplr.f user.f (only (a dummy) FUNCTION INIDOM and SUBROUTINE CHSPCM) src.f ilubsn.f blas.f (if the BLAS library is not available on the platform) The results for this run is in the file exmplr_output To get an optimal vector performance for a small number of PDE components one should use the specific ilubs#.f code, in this case: ilubs2.f SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test ! -d 'Dp' then mkdir 'Dp' fi cd 'Dp' if test -f 'prtsol.f' then echo shar: will not over-write existing file "'prtsol.f'" else cat << \SHAR_EOF > 'prtsol.f' PROGRAM PRTSOL C C----------------------------------------------------------------------- C Ccc This program reads a file made by subroutine DUMP and prints the C solution on an output file. Both filenames are read from standard C input. C Ccc EXTERNALS USED: EXTERNAL PRSOL, RDDUMP C C Ccc INCLUDE 'CMNWRITEF' C C CMNWRITEF C C COMMON needed for continuation calls INTEGER MAXLVW, NPDEW, LRWKPS, LIWKPS, LRWKB, LIWKB LOGICAL FIRST, SECOND DOUBLE PRECISION T0, TW, TEW, DTW, XLW, YLW, XRW, YUW, DXB, DYB, + DTO COMMON /WRITIF/ MAXLVW, NPDEW, LRWKPS, LIWKPS, LRWKB, LIWKB COMMON /WRITLF/ FIRST, SECOND COMMON /WRITRF/ T0, TW, TEW, DTW, XLW,YLW, XRW,YUW, DXB, DYB, DTO SAVE /WRITIF/, /WRITLF/, /WRITRF/ C C end INCLUDE 'CMNWRITEF' C C C----------------------------------------------------------------------- C INTEGER MXLEV, NPD, NPTS, LENIWK, LENRWK PARAMETER (MXLEV=5, NPD=3, NPTS=10000) PARAMETER (LENIWK=NPTS*(7*MXLEV+20), + LENRWK=5*NPTS*NPD*MXLEV) C CHARACTER FILE*128 INTEGER IWK(LENIWK), + LSGNM1, LSGN, LSGNP1, LSUNM1, LSSN, LSUN DOUBLE PRECISION RWK(LENRWK) PRINT *, 'DUMP file?' READ '(A)', FILE C OPEN(UNIT=62,FILE=FILE,FORM='UNFORMATTED') CALL RDDUMP (62, RWK, LENRWK, IWK, LENIWK) CLOSE(62) C C Setup work storage LSGNM1 = 1 LSGN = LSGNM1 + MAXLVW+1 LSGNP1 = LSGN + MAXLVW+1 LSUNM1 = LSGNP1 + MAXLVW+1 LSSN = LSUNM1 + MAXLVW LSUN = LSSN + MAXLVW C C call print routine PRINT *, 'output file?' READ '(A)', FILE C OPEN(UNIT=61,FILE=FILE) CALL PRSOL (61, TW, NPDEW, XLW, YLW, DXB, DYB, + IWK(LSGN), IWK(LIWKPS), IWK(LSUN), RWK(LRWKPS)) CLOSE(61) END SUBROUTINE RDDUMP (LUNDMP, RWK, LENRWK, IWK, LENIWK) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER LENIWK INTEGER LUNDMP, LENRWK, IWK(LENIWK) DOUBLE PRECISION RWK(LENRWK) C Ccc PURPOSE: C Read all information necessary for a restart of VLUGR2 from file C Ccc PARAMETER DESCRIPTION: C LUNDMP : IN. Logical unit number of dumpfile. Should be opened as an C unformatted file. C RWK : OUT. Real workstorage intended to pass to VLUGR2 C LENRWK : IN. Dimension of RWK. C IWK : OUT. Integer workstorage intended to pass to VLUGR2 C LENIWK : IN. Dimension of IWK. C Ccc EXTERNALS USED: NONE C C Ccc INCLUDE 'CMNSTATS' C C CMNSTATS C C COMMON with integration statistics INTEGER MXCLEV, MXCNIT PARAMETER (MXCLEV = 10, MXCNIT = 20) INTEGER LUNPDS, LUNNLS, LUNLSS, LEVEL, NSTEPS, NREJS, + NJACS(MXCLEV), NRESID(MXCLEV), NNIT(MXCLEV), + NLSIT(MXCLEV,MXCNIT) COMMON /STATS/ LUNPDS, LUNNLS, LUNLSS, LEVEL, NSTEPS, NREJS, + NJACS, NRESID, NNIT, NLSIT SAVE /STATS/ C C end INCLUDE 'CMNSTATS' C C Ccc INCLUDE 'CMNWRITEF' C C CMNWRITEF C C COMMON needed for continuation calls INTEGER MAXLVW, NPDEW, LRWKPS, LIWKPS, LRWKB, LIWKB LOGICAL FIRST, SECOND DOUBLE PRECISION T0, TW, TEW, DTW, XLW, YLW, XRW, YUW, DXB, DYB, + DTO COMMON /WRITIF/ MAXLVW, NPDEW, LRWKPS, LIWKPS, LRWKB, LIWKB COMMON /WRITLF/ FIRST, SECOND COMMON /WRITRF/ T0, TW, TEW, DTW, XLW,YLW, XRW,YUW, DXB, DYB, DTO SAVE /WRITIF/, /WRITLF/, /WRITRF/ C C end INCLUDE 'CMNWRITEF' C C C----------------------------------------------------------------------- C INTEGER I, J READ(LUNDMP) MAXLVW, NPDEW, LRWKPS, LIWKPS, LRWKB, LIWKB, + FIRST, SECOND, + T0, TW, TEW, DTW, XLW, YLW, XRW, YUW, DXB, DYB, DTO IF (LENRWK .LT. LRWKPS+LRWKB .OR. LENIWK .LT. LIWKPS+LIWKB) THEN PRINT *, LENRWK, LRWKPS+LRWKB, LENIWK, LIWKPS+LIWKB STOP 'work space too small' ENDIF READ(LUNDMP) LUNPDS, LUNNLS, LUNLSS, LEVEL, NSTEPS, NREJS, + (NJACS(I), I=1,MXCLEV), (NRESID(I), I=1,MXCLEV), + (NNIT(I), I=1,MXCLEV), ((NLSIT(I,J), I=1,MXCLEV), J=1,MXCNIT) READ(LUNDMP) (RWK(I), I=1,LRWKPS+LRWKB) READ(LUNDMP) (IWK(I), I=1,LIWKPS+LIWKB) C RETURN END SUBROUTINE PRSOL (LUN, T, NPDE, XL, YL, DXB, DYB, LGRID, ISTRUC, + LSOL, SOL) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER LUN, NPDE, LGRID(0:*), ISTRUC(*), LSOL(*) DOUBLE PRECISION T, XL, YL, DXB, DYB, SOL(*) C Ccc PURPOSE: C Print solution and coordinate values at all grid levels. C Ccc PARAMETER DESCRIPTION: C LUN : IN. Logical unit number of print file C T : IN. Current value of time variable C NPDE : IN. # PDE components C XL : IN. X-coordinate of lowerleft corner of (virtual) domain C YL : IN. Y-coordinate of lowerleft corner of (virtual) domain C DXB : IN. Cell width in X-direction of base grid C DYB : IN. Cell width in Y-direction of base grid C LGRID : IN. (0:*) C LGRID(0) = max. grid level used at T C LGRID(1): pointer to base grid structure ISTRUC C LGRID(LEVEL): pointer to grid structure (LROW, IROW, ICOL) C of refinement level LEVEL for time T C ISTRUC : IN. (*) C ISTRUC(LGRID(LEVEL):.) contains (LROW,IROW,ICOL) of grid C level LEVEL, C LROW : (0:LROW(0)+1) C LROW(0) = NROWS: Actual # rows in grid C LROW(1:NROWS): pointers to the start of a row in the grid C LROW(NROWS+1) = NPTS+1: Actual # nodes in grid + 1 C IROW : (NROWS) C IROW(IR): row number of row IR in virtual rectangle C ICOL : (NPTS) C ICOL(IPT): column number of grid point IPT in virtual C rectangle C LSOL : IN. (*) C LSOL(LEVEL): pointer to (injected) solution at grid C of refinement level LEVEL for time T C SOL : IN. (*) C SOL(LSOL(LEVEL)+1:LSOL(LEVEL)+NPTS(LEVEL)*NPDE) contains C U_LEVEL(NPTS,NPDE) C Ccc EXTERNALS USED: EXTERNAL PRSOLL C C----------------------------------------------------------------------- C INTEGER MAXLEV, LEVEL, LLROW, NROWS, NPTS, LIROW, LICOL DOUBLE PRECISION DX, DY MAXLEV = LGRID(0) DX = DXB DY = DYB DO 10 LEVEL = 1, MAXLEV LLROW = LGRID(LEVEL) NROWS = ISTRUC(LLROW) NPTS = ISTRUC(LLROW+NROWS+1)-1 LIROW = LLROW+NROWS+2 LICOL = LIROW+NROWS CALL PRSOLL (LUN, LEVEL, T, NPTS, NPDE, XL, YL, DX, DY, + ISTRUC(LLROW), ISTRUC(LIROW), ISTRUC(LICOL), + SOL(LSOL(LEVEL)+1)) DX = DX/2 DY = DY/2 10 CONTINUE RETURN END SUBROUTINE PRSOLL (LUN, LEVEL, T, NPTS, NPDE, XL, YL, DX, DY, + LROW, IROW, ICOL, U) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER LUN, LEVEL, NPTS, NPDE, LROW(0:*), IROW(*), ICOL(*) DOUBLE PRECISION T, XL, YL, DX, DY, U(NPTS,NPDE) C Ccc PURPOSE: C Print solution and X- and Y-coordinates of gridlevel LEVEL. C Ccc PARAMETER DESCRIPTION: C LUN : IN. Logical unit number of print file C LEVEL : IN. Grid level corresponding with solution U. C T : IN. Current value of time variable C NPTS : IN. # grid points at this level C NPDE : IN. # PDE components C XL : IN. X-coordinate of lower-left point of virtual rectangle C YL : IN. Y-coordinate of lower-left point of virtual rectangle C DX : IN. Grid width in X-direction C DY : IN. Grid width in Y-direction C LROW : IN. (0:LROW(0)+1) C LROW(0) = NROWS: Actual # rows in grid C LROW(1:NROWS): pointers to the start of a row in the grid C LROW(NROWS+1) = NPTS+1: Actual # nodes in grid + 1 C IROW : IN. (NROWS) C IROW(IR): row number of row IR in virtual rectangle C ICOL : IN. (NPTS) C ICOL(IPT): column number of grid point IPT in virtual C rectangle C U : IN. Solution on this grid level C Ccc EXTERNALS USED: NONE C C----------------------------------------------------------------------- C INTEGER IC, IPT, IR, NROWS DOUBLE PRECISION X, Y C NROWS = LROW(0) WRITE(LUN,'(//// T10,A,T30,A,T46,A,T62,A,T71,A //)') + 'Level', 't', 'Y', 'X', 'Solution' IR = 1 Y = YL + IROW(IR)*DY IPT = LROW(IR) X = XL + ICOL(IPT)*DX WRITE(LUN, + '(T13,I2,T21,E12.5,T37,E12.5,T53,E12.5,T69,E12.5)') + LEVEL, T, Y, X, U(IPT,1) DO 10 IC = 2, NPDE WRITE(LUN,'(T69,E12.5)') U(IPT,IC) 10 CONTINUE DO 20 IPT = LROW(IR)+1, LROW(IR+1)-1 X = XL + ICOL(IPT)*DX WRITE(LUN,'(T53,E12.5,T69,E12.5)') X, U(IPT,1) DO 30 IC = 2, NPDE WRITE(LUN,'(T69,E12.5)') U(IPT,IC) 30 CONTINUE 20 CONTINUE DO 40 IR = 2, NROWS Y = YL + IROW(IR)*DY IPT = LROW(IR) X = XL + ICOL(IPT)*DX WRITE(LUN, + '(T21,E12.5,T37,E12.5,T53,E12.5,T69,E12.5)') + T, Y, X, U(IPT,1) DO 50 IC = 2, NPDE WRITE(LUN,'(T69,E12.5)') U(IPT,IC) 50 CONTINUE DO 60 IPT = LROW(IR)+1, LROW(IR+1)-1 X = XL + ICOL(IPT)*DX WRITE(LUN,'(T53,E12.5,T69,E12.5)') X, U(IPT,1) DO 70 IC = 2, NPDE WRITE(LUN,'(T69,E12.5)') U(IPT,IC) 70 CONTINUE 60 CONTINUE 40 CONTINUE RETURN END SHAR_EOF fi # end of overwriting check if test -f 'probi.f' then echo shar: will not over-write existing file "'probi.f'" else cat << \SHAR_EOF > 'probi.f' PROGRAM EXMPL C C Ccc INCLUDE 'PARNEWTON' C C PARNEWTON C C Parameters for Newton process C MAXNIT : Max. number of Newton iterations C MAXJAC : Max. number of Jacobian / preconditioner evaluations during C a Newton process C TOLNEW : Tolerance for Newton process: C rho/(1-rho)*|| corr.||_w < TOLNEW INTEGER MAXNIT, MAXJAC DOUBLE PRECISION TOLNEW PARAMETER (MAXNIT = 10, MAXJAC = 2, TOLNEW = 1.0) C C end INCLUDE 'PARNEWTON' C C Ccc INCLUDE 'PARGCRO' C C PARGCRO C C Parameters for linear system solver GCRO + (block-)diagonal C preconditioner C IDIAGP : 0: block-diagonal + first order derivatives C 1: block-diagonal neglecting first order derivatives C 2: diagonal + first order derivatives C 3: diagonal neglecting first order derivatives C NRRMAX : Max. number of restarts of outer loop C MAXLR : Max. number of iterations in outer loop C MAXL : Max. number of iterations in GMRES inner loop C TOLLSC : Tolerance for linear system solver INTEGER IDIAGP, NRRMAX, MAXLR, MAXL DOUBLE PRECISION TOLLSC PARAMETER (NRRMAX = 1, MAXLR = 5, MAXL = 20) C PARAMETER (NRRMAX = 1, MAXLR = 3, MAXL = 10) PARAMETER (TOLLSC = TOLNEW/10) COMMON /IGCRO/ IDIAGP SAVE /IGCRO/ C C end INCLUDE 'PARGCRO' C C INTEGER MXLEV, NPD, NPTS, LENIWK, LENRWK, LENLWK PARAMETER (MXLEV=2, NPD=3, NPTS=5000) PARAMETER (LENIWK=NPTS*(5*MXLEV+14), + LENRWK=NPTS*NPD*(5*MXLEV+9 + + 9*NPD+(2*MAXLR+MAXL+6+NPD)), + LENLWK=2*NPTS) C C----------------------------------------------------------------------- C INTEGER LUNDMP PARAMETER (LUNDMP = 89) C CHARACTER FILE*7 INTEGER NPDE, INFO(7), IWK(LENIWK), MNTR, I LOGICAL LWK(LENLWK) DOUBLE PRECISION T, TOUT(4), DT, XL, YL, XR, YU, DX, DY, + TOLS, TOLT, RINFO(2+3*NPD), RWK(LENRWK) C C First call of VLUGR2 MNTR = 0 NPDE = 3 T = 0.0 TOUT(1) = 500.0 TOUT(2) = 5000.0 TOUT(3) = 10000.0 TOUT(4) = 20000.0 DT = 0.1 XL = 0.0 XR = 1.0 YL = 0.0 YU = 1.0 DX = 0.05 DY = 0.05 TOLS = 0.1 TOLT = 0.1 INFO(1) = 1 C MAXLEV INFO(2) = 3 C Domain is a rectangle INFO(3) = 0 C Linear system solver PRINT *, 'Lin.sys.solver; BiCGStab, GCRO or matrix-free GCRO ?' PRINT *, ' (0 / 10,11,12,13 / 20,21,22,23 ) ?' READ *, INFO(4) OPEN (UNIT=61,FILE='RunInfo') C Write integration history to unit # 61 INFO(5) = 61 C Write Newton info to unit # 61 INFO(6) = 61 C Write Linear system solver info to unit # 61 INFO(7) = 61 C DTMIN = 1D-3 RINFO(1) = 1.0D-3 C DTMAX = 1.0 RINFO(2) = 10000.0 C UMAX RINFO(3) = 1.1D+5 RINFO(4) = 0.25 RINFO(5) = 292.0 C SPCWGT = 1.0 RINFO(6) = 1.0 RINFO(7) = 1.0 RINFO(8) = 1.0 C TIMWGT = 1.0 RINFO( 9) = 1.0 RINFO(10) = 1.0 RINFO(11) = 1.0 C C Call main routine FILE='DUMP' DO 10 I = 1, 4 CALL VLUGR2 (NPDE, T, TOUT(I), DT, XL, YL, XR, YU, DX, DY, + TOLS, TOLT, INFO, RINFO, RWK, LENRWK, IWK, LENIWK, + LWK, LENLWK, MNTR) C C Save info on file WRITE(FILE(5:7),'(I3.3)') I OPEN(UNIT=LUNDMP,FILE=FILE,FORM='UNFORMATTED') CALL DUMP (LUNDMP, RWK, IWK) CLOSE(LUNDMP) C Check MNTR value IF (MNTR .NE. 1) THEN PRINT *, 'VLUGR2 returned with MNTR=', MNTR STOP ENDIF 10 CONTINUE END SUBROUTINE PDEIV (T, X, Y, U, NPTS, NPDE) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER NPTS, NPDE DOUBLE PRECISION T, X(NPTS), Y(NPTS), U(NPTS,NPDE) C Ccc PURPOSE: C Define (initial) solution of PDE. C Ccc PARAMETER DESCRIPTION: C T : IN. Time at which (initial) solution should be given C X : IN. Array of X-coordinates for the gridpoints C Y : IN. Array of Y-coordinates for the gridpoints C U : OUT. Array of PDE component values for the gridpoints. C NPTS : IN. Number of gridpoints C NPDE : IN. # PDE components C C----------------------------------------------------------------------- C INTEGER I C DOUBLE PRECISION N, GAMMA, MU0, RHO0, P0, W0, G, DM, KAPPA, AL, + AT, CF, TKAPPA, LT, LL, CS, RHOS, T0, ALPHA, BETA, QC, TC COMMON /PROBLM/ N, GAMMA, MU0, RHO0, P0, W0, G, DM, KAPPA, AL, AT, + CF, TKAPPA, LT, LL, CS, RHOS, T0, ALPHA, BETA, QC, TC SAVE /PROBLM/ C Ccc Problem parameters N = 0.4 KAPPA = 1.0D-10 G = 9.81 DM = 0.0 AT = 0.002 AL = 0.01 CF = 4182.0 TKAPPA = 4.0 LT = 0.001 LL = 0.01 CS = 840.0 RHOS = 2500.0 RHO0 = 1.0D+3 T0 = 290.0 P0 = 1.0D+5 ALPHA = -3.0D-4 BETA = 4.45D-10 GAMMA = LOG(1.2) MU0 = 1.0D-3 W0 = 0.25 QC = 1.0D-4 TC = 292.0 C Ccc Initial solution DO 10 I = 1, NPTS U(I,1) = P0 + (1.0 - Y(I))*RHO0*G U(I,2) = 0.0 U(I,3) = T0 10 CONTINUE RETURN END SUBROUTINE PDEF (T, X, Y, U, UT, UX, UY, UXX, UXY, UYY, RES, + NPTS, NPDE) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER NPTS, NPDE DOUBLE PRECISION T, X(NPTS), Y(NPTS), U(NPTS,NPDE), + UT(NPTS,NPDE), UX(NPTS,NPDE), UY(NPTS,NPDE), + UXX(NPTS,NPDE), UXY(NPTS,NPDE), UYY(NPTS,NPDE), + RES(NPTS,NPDE) C Ccc PURPOSE: C Define residual of PDE on interior of domain. Boundary values will be C overwritten later on. C Ccc PARAMETER DESCRIPTION: C T : IN. Time at which residual should be evaluated C X : IN. Array of X-coordinates for the gridpoints C Y : IN. Array of Y-coordinates for the gridpoints C U : IN. Array of PDE components for the gridpoints. C UT : IN. Array of time derivative of PDE components C UX : IN. -I C UY : IN. I C UXX : IN. I Arrays containing space derivatives of PDE components C UXY : IN. I C UYY : IN. -I C RES : OUT. Array containg PDE residual at gridpoints in interior of C domain. The residual values at boundary points will be C overwritten by a call to PDEBC. C NPTS : IN. Number of gridpoints C NPDE : IN. Number of PDE components C C C Ccc INCLUDE 'CMNCMMACH' C C CMNCMMACH C C COMMON with `machine numbers' C LUNOUT : Logical unit # standard output -I C LUNERR : Logical unit # standard error I Set in the routine C UROUND : Smallest machine number such that I MACNUM C 1.0+UROUND > 1.0 and 1.0-UROUND < 1.0 I C XMIN : Smallest floating-point number -I INTEGER LUNOUT, LUNERR DOUBLE PRECISION UROUND, XMIN COMMON /IMACH/ LUNOUT, LUNERR COMMON /RMACH/ UROUND, XMIN SAVE /IMACH/, /RMACH/ C C end INCLUDE 'CMNCMMACH' C C----------------------------------------------------------------------- C DOUBLE PRECISION N, GAMMA, MU0, RHO0, P0, W0, G, DM, KAPPA, AL, + AT, CF, TKAPPA, LT, LL, CS, RHOS, T0, ALPHA, BETA, QC, TC COMMON /PROBLM/ N, GAMMA, MU0, RHO0, P0, W0, G, DM, KAPPA, AL, AT, + CF, TKAPPA, LT, LL, CS, RHOS, T0, ALPHA, BETA, QC, TC SAVE /PROBLM/ C INTEGER I DOUBLE PRECISION P, PT, PX, PY, W, WT, WX, WY, T1, TT, TX, TY, + RHO, RHOX, RHOY, + MU, MUX, MUY, KAPMU, KAPMU2, KAPMUX, KAPMUY, Q1, Q2, QL, + ND11, ND12, ND22, H11, H12, H22, + PXX, PXY, PYY, WXX, WXY, WYY, TXX, TXY, TYY, + ND11Q1, ND11Q2, ND12Q1, ND12Q2, ND22Q1, ND22Q2, + H11Q1, H11Q2, H12Q1, H12Q2, H22Q1, H22Q2, Q1X, Q1Y, Q2X, Q2Y, + ND11X, ND12X, ND12Y, ND22Y, JW1, JW2, JW1X, JW2Y, + H11X, H12X, H12Y, H22Y, JT1X, JT2Y C DO 10 I = 1, NPTS P = U(I,1) PT = UT(I,1) PX = UX(I,1) PY = UY(I,1) W = U(I,2) WT = UT(I,2) WX = UX(I,2) WY = UY(I,2) T1 = U(I,3) TT = UT(I,3) TX = UX(I,3) TY = UY(I,3) RHO = RHO0*EXP(ALPHA*(T1-T0)+BETA*(P-P0)+GAMMA*W) RHOX = RHO*(ALPHA*TX+BETA*PX+GAMMA*WX) RHOY = RHO*(ALPHA*TY+BETA*PY+GAMMA*WY) MU = MU0*(1+1.85*W-4.0*W*W) MUX = MU0*(1.85*WX-8.0*W*WX) MUY = MU0*(1.85*WY-8.0*W*WY) KAPMU = KAPPA/MU KAPMU2 = -KAPMU/MU KAPMUX = KAPMU2*MUX KAPMUY = KAPMU2*MUY Q1 = -KAPMU*PX Q2 = -KAPMU*(PY+RHO*G) QL = MAX(SQRT(Q1*Q1+Q2*Q2),UROUND) ND11 = N*DM + AT*QL + (AL-AT)*Q1*Q1/QL ND12 = (AL-AT)*Q1*Q2/QL ND22 = N*DM + AT*QL + (AL-AT)*Q2*Q2/QL H11 = TKAPPA + LT*QL + (LL-LT)*Q1*Q1/QL H12 = (LL-LT)*Q1*Q2/QL H22 = TKAPPA + LT*QL + (LL-LT)*Q2*Q2/QL PXX = UXX(I,1) PXY = UXY(I,1) PYY = UYY(I,1) WXX = UXX(I,2) WXY = UXY(I,2) WYY = UYY(I,2) TXX = UXX(I,3) TXY = UXY(I,3) TYY = UYY(I,3) ND11Q1 = (AT + (AL-AT)*(2-(Q1/QL)**2))*Q1/QL ND11Q2 = (AT - (AL-AT)*((Q1/QL)**2))*Q2/QL ND12Q1 = (AL-AT)*(Q2/QL)**3 ND12Q2 = (AL-AT)*(Q1/QL)**3 ND22Q1 = (AT - (AL-AT)*((Q2/QL)**2))*Q1/QL ND22Q2 = (AT + (AL-AT)*(2-(Q2/QL)**2))*Q2/QL H11Q1 = (LT + (LL-LT)*(2-(Q1/QL)**2))*Q1/QL H11Q2 = (LT - (LL-LT)*((Q1/QL)**2))*Q2/QL H12Q1 = (LL-LT)*(Q2/QL)**3 H12Q2 = (LL-LT)*(Q1/QL)**3 H22Q1 = (LT - (LL-LT)*((Q2/QL)**2))*Q1/QL H22Q2 = (LT + (LL-LT)*(2-(Q2/QL)**2))*Q2/QL Q1X = -(KAPMUX*PX+KAPMU*PXX) Q1Y = -(KAPMUY*PX+KAPMU*PXY) Q2X = -(KAPMUX*(PY+RHO*G)+KAPMU*(PXY+RHOX*G)) Q2Y = -(KAPMUY*(PY+RHO*G)+KAPMU*(PYY+RHOY*G)) ND11X = ND11Q1*Q1X + ND11Q2*Q2X ND12X = ND12Q1*Q1X + ND12Q2*Q2X ND12Y = ND12Q1*Q1Y + ND12Q2*Q2Y ND22Y = ND22Q1*Q1Y + ND22Q2*Q2Y JW1 = -(ND11*WX + ND12*WY) JW2 = -(ND12*WX + ND22*WY) JW1X = -(ND11X*WX+ND11*WXX + ND12X*WY+ND12*WXY) JW2Y = -(ND12Y*WX+ND12*WXY + ND22Y*WY+ND22*WYY) H11X = H11Q1*Q1X + H11Q2*Q2X H12X = H12Q1*Q1X + H12Q2*Q2X H12Y = H12Q1*Q1Y + H12Q2*Q2Y H22Y = H22Q1*Q1Y + H22Q2*Q2Y JT1X = -(H11X*TX+H11*TXX + H12X*TY+H12*TXY) JT2Y = -(H12Y*TX+H12*TXY + H22Y*TY+H22*TYY) C RES(I,1) = N*RHO*(BETA*PT+GAMMA*WT+ALPHA*TT) + + RHOX*Q1+RHO*Q1X + RHOY*Q2+RHO*Q2Y RES(I,2) = N*RHO*WT + + RHO*Q1*WX + RHO*Q2*WY + + RHOX*JW1+RHO*JW1X + RHOY*JW2+RHO*JW2Y RES(I,3) = (N*CF*RHO+(1-N)*RHOS*CS)*TT + + RHO*CF*Q1*TX + RHO*CF*Q2*TY + + JT1X + JT2Y 10 CONTINUE RETURN END SUBROUTINE PDEBC (T, X, Y, U, UT, UX, UY, RES, + NPTS, NPDE, LLBND, ILBND, LBND) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER NPTS, NPDE, LLBND(0:*), ILBND(*), LBND(*) DOUBLE PRECISION T, X(NPTS), Y(NPTS), U(NPTS,NPDE), + UT(NPTS,NPDE), UX(NPTS,NPDE), UY(NPTS,NPDE), + RES(NPTS,NPDE) C Ccc PURPOSE: C Define residual of boundary equations of PDE. The residual on interior C points has already been stored in RES. C Ccc PARAMETER DESCRIPTION: C T : IN. Time at which BC's should be evaluated C X : IN. Array of X-coordinates for the gridpoints C Y : IN. Array of Y-coordinates for the gridpoints C U : IN. Array of PDE components for the gridpoints. C UT : IN. Array of time derivative of PDE components C UX : IN. -I C UY : IN. -I Arrays containing space derivatives of PDE components C RES : INOUT. C IN: PDE residual for interior points (set by PDEF) C OUT: Array with PDE residual at physical boundary points C inserted C NPTS : IN. Number of grid components C NPDE : IN. Number of PDE components C LLBND : (0:LLBND(0)+2) C LLBND(0) = NBNDS: total # physical boundaries and corners in C actual domain. C NB. corners should be stored as an independent boundary C (cf. ILBND). The order in LLBND should be first the C boundaries and then the corners. C LLBND(1:NBNDS): pointers to a specific boundary or corner in C LBND C LLBND(NBNDS+1) = NBDPTS+1: total # physical boundary points C in LBND + 1 C LLBND(NBNDS+1): pointer to internal boundary in LBND C LLBND(NBNDS+2) = NBIPTS+1: total # points in LBND + 1 C ILBND : (NBNDS) C ILBND(IB): type of boundary: C 1: Lower boundary -I C 2: Left boundary I C 3: Upper boundary I max. first order derivative C 4: Right boundary -I C 12: Lowerleft corner -I C 23: Leftupper corner I corners of 90 degrees C 34: Upperright corner I (external corners) C 41: Rightlower corner -I max. first order deriv. C 21: Leftlower corner -I C 32: Upperleft corner I corners of 270 degrees C 43: Rightupper corner I (internal corners) C 14: Lowerright corner -I max. first order deriv. C LBND : IN. (NBDPTS) C LBND(LB): pointer to boundary point in actual grid C structure (as in X, Y, and U) C C Ccc INCLUDE 'CMNCMMACH' C C CMNCMMACH C C COMMON with `machine numbers' C LUNOUT : Logical unit # standard output -I C LUNERR : Logical unit # standard error I Set in the routine C UROUND : Smallest machine number such that I MACNUM C 1.0+UROUND > 1.0 and 1.0-UROUND < 1.0 I C XMIN : Smallest floating-point number -I INTEGER LUNOUT, LUNERR DOUBLE PRECISION UROUND, XMIN COMMON /IMACH/ LUNOUT, LUNERR COMMON /RMACH/ UROUND, XMIN SAVE /IMACH/, /RMACH/ C C end INCLUDE 'CMNCMMACH' C C C----------------------------------------------------------------------- C C DOUBLE PRECISION N, GAMMA, MU0, RHO0, P0, W0, G, DM, KAPPA, AL, + AT, CF, TKAPPA, LT, LL, CS, RHOS, T0, ALPHA, BETA, QC, TC COMMON /PROBLM/ N, GAMMA, MU0, RHO0, P0, W0, G, DM, KAPPA, AL, AT, + CF, TKAPPA, LT, LL, CS, RHOS, T0, ALPHA, BETA, QC, TC SAVE /PROBLM/ C INTEGER I, J, K, NBNDS DOUBLE PRECISION P, PY, W, T1, RHO, MU, KAPMU, Q2 C NBNDS = LLBND(0) DO 10 J = 1, NBNDS IF (ILBND(J) .EQ. 1) THEN C C yL boundary: dp/dy - rho.g2 = 0 C dw/dy = 0 0 'probii.f' PROGRAM EXMPL C C Ccc INCLUDE 'PARNEWTON' C C PARNEWTON C C Parameters for Newton process C MAXNIT : Max. number of Newton iterations C MAXJAC : Max. number of Jacobian / preconditioner evaluations during C a Newton process C TOLNEW : Tolerance for Newton process: C rho/(1-rho)*|| corr.||_w < TOLNEW INTEGER MAXNIT, MAXJAC DOUBLE PRECISION TOLNEW PARAMETER (MAXNIT = 10, MAXJAC = 2, TOLNEW = 1.0) C C end INCLUDE 'PARNEWTON' C C Ccc INCLUDE 'PARGCRO' C C PARGCRO C C Parameters for linear system solver GCRO + (block-)diagonal C preconditioner C IDIAGP : 0: block-diagonal + first order derivatives C 1: block-diagonal neglecting first order derivatives C 2: diagonal + first order derivatives C 3: diagonal neglecting first order derivatives C NRRMAX : Max. number of restarts of outer loop C MAXLR : Max. number of iterations in outer loop C MAXL : Max. number of iterations in GMRES inner loop C TOLLSC : Tolerance for linear system solver INTEGER IDIAGP, NRRMAX, MAXLR, MAXL DOUBLE PRECISION TOLLSC PARAMETER (NRRMAX = 1, MAXLR = 5, MAXL = 20) C PARAMETER (NRRMAX = 1, MAXLR = 3, MAXL = 10) PARAMETER (TOLLSC = TOLNEW/10) COMMON /IGCRO/ IDIAGP SAVE /IGCRO/ C C end INCLUDE 'PARGCRO' C INTEGER MXLEV, NPD, NPTS, LENIWK, LENRWK, LENLWK PARAMETER (MXLEV=2, NPD=2, NPTS=10000) PARAMETER (LENIWK=NPTS*(5*MXLEV+14), + LENRWK=NPTS*NPD*(5*MXLEV+9 + + 9*NPD+(2*MAXLR+MAXL+6+NPD)), + LENLWK=2*NPTS) C C----------------------------------------------------------------------- C INTEGER LUNDMP PARAMETER (LUNDMP = 89) C CHARACTER FILE*7 INTEGER NPDE, INFO(7), IWK(LENIWK), MNTR, I LOGICAL LWK(LENLWK) DOUBLE PRECISION T, TOUT(4), DT, XL, YL, XR, YU, DX, DY, + TOLS, TOLT, RINFO(2+3*NPD), RWK(LENRWK) C C First call of VLUGR2 MNTR = 0 NPDE = 2 T = 0.0 TOUT(1) = 10000.0 TOUT(2) = 20000.0 TOUT(3) = 30000.0 TOUT(4) = 100000.0 DT = 0.1 C Since domain is not a rectangle the domain parameters have not to be C specified here TOLS = 0.1 TOLT = 0.1 INFO(1) = 1 C MAXLEV INFO(2) = 3 C Domain is not a rectangle INFO(3) = 1 C Linear system solver PRINT *, 'Lin.sys.solver; BiCGStab, GCRO or matrix-free GCRO ?' PRINT *, ' (0 / 10,11,12,13 / 20,21,22,23 ) ?' READ *, INFO(4) OPEN (UNIT=61,FILE='RunInfo') C Write integration history to unit # 61 INFO(5) = 61 C Write Newton info to unit # 61 INFO(6) = 61 C Write Linear system solver info to unit # 61 INFO(7) = 61 C DTMIN RINFO(1) = 1.0D-3 C DTMAX RINFO(2) = 50000.0 C UMAX RINFO(3) = 1.1D+5 RINFO(4) = 0.25 C SPCWGT = 1.0 RINFO(5) = 1.0 RINFO(6) = 1.0 C TIMWGT = 1.0 RINFO(7) = 1.0 RINFO(8) = 1.0 C C Call main routine FILE='DUMP' DO 10 I = 1, 4 CALL VLUGR2 (NPDE, T, TOUT(I), DT, XL, YL, XR, YU, DX, DY, + TOLS, TOLT, INFO, RINFO, RWK, LENRWK, IWK, LENIWK, + LWK, LENLWK, MNTR) C C Save info on file WRITE(FILE(5:7),'(I3.3)') I OPEN(UNIT=LUNDMP,FILE=FILE,FORM='UNFORMATTED') CALL DUMP (LUNDMP, RWK, IWK) CLOSE(LUNDMP) C Check MNTR value IF (MNTR .NE. 1) THEN PRINT *, 'VLUGR2 returned with MNTR=', MNTR STOP ENDIF 10 CONTINUE END LOGICAL FUNCTION INIDOM (MAXPTS, XL, YL, XR, YU, DX, DY, + LROW, IROW, ICOL, LLBND, ILBND, LBND) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER MAXPTS, LROW(0:*), IROW(*), ICOL(*), + LLBND(0:*), ILBND(*), LBND(*) DOUBLE PRECISION XL, YL, XR, YU, DX, DY C Ccc PURPOSE: C Define initial domain. NB. Boundaries should consist of as many points C as are necessary to employ second order space discretization, i.e., C a boundary enclosing the internal part of the domain should not C include less than 3 grid points including the corners. If Neumann C boundaries are used the minimum is 4 since otherwise the Jacobian C matrix will be singular. C C A (virtual) rectangle is placed upon the (irregular) domain. The C lowerleft point of this rectangle is (XL,YL) in physical coordinates C and (0,0) in column, resp. row coordinates. The upperright point is C (XR,YU) resp. (Nx, Ny), where Nx = (XR-XL)/DX and Ny = (YU-YL)/DY. C Only real grid points are stored. C The coordinate values of the initial grid should be stored rowwise C in LROW, IROW, ICOL. C Pointers to the boundary points should be stored in a list together C with the type of the boundary. (LLBND, ILBND, LBND) C C On exit INIDOM = .FALSE. if the # grid points required is larger C than MAXPTS and MAXPTS is set to the required # points. C Ccc PARAMETER DESCRIPTION: C MAXPTS : INOUT. C IN: Max. # grid points allowed by the available workspace C OUT: # grid points required, if larger than # points allowed C XL : OUT. X-coordinate of lower-left point of virtual rectangle C YL : OUT. Y-coordinate of lower-left point of virtual rectangle C XR : OUT. X-coordinate of upper-right point of virtual rectangle C YU : OUT. Y-coordinate of upper-right point of virtual rectangle C DX : OUT. Grid width in X-direction C DY : OUT. Grid width in Y-direction C LROW : OUT. INTEGER array of dimension (0:LROW(0)+1) C LROW(0) = NROWS: Actual # rows in grid C LROW(1:NROWS): pointers to the start of a row in the grid C structure C LROW(NROWS+1) = NPTS+1: Actual # nodes in grid + 1 C IROW : OUT. INTEGER array of dimension (NROWS) C IROW(IR): row number of row IR in virtual rectangle C ICOL : OUT. INTEGER array of dimension (NPTS) C ICOL(IPT): column number of grid point IPT in virtual C rectangle C LLBND : (0:LLBND(0)+2) C LLBND(0) = NBNDS: total # physical boundaries and corners in C actual domain. C NB. corners should be stored as an independent boundary C (cf. ILBND). The order in LLBND should be first the C boundaries and then the corners. C LLBND(1:NBNDS): pointers to a specific boundary or corner in C LBND C LLBND(NBNDS+1) = NBDPTS+1: total # physical boundary points C in LBND + 1 C LLBND(NBNDS+1): pointer to internal boundary in LBND C LLBND(NBNDS+2) = NBIPTS+1: total # points in LBND + 1 C ILBND : (NBNDS) C ILBND(IB): type of boundary: C 1: Lower boundary -I C 2: Left boundary I C 3: Upper boundary I max. first order derivative C 4: Right boundary -I C 12: Lowerleft corner -I C 23: Leftupper corner I corners of 90 degrees C 34: Upperright corner I (external corners) C 41: Rightlower corner -I max. first order deriv. C 21: Leftlower corner -I C 32: Upperleft corner I corners of 270 degrees C 43: Rightupper corner I (internal corners) C 14: Lowerright corner -I max. first order deriv. C LBND : OUT. INTEGER array of dimension (NBDPTS) C LBND(IBPT): pointer to boundary point in actual grid C structure C Ccc EXTERNALS USED: NONE C C----------------------------------------------------------------------- C C NX should be even and NY a quintuple INTEGER NX, NY PARAMETER (NX = 20, NY = 20) C INTEGER I, IPT, J, NROWS, NPTS, NBNDS, NX1, NY1, NY2 C Ccc Make initial grid, check MAXPTS against rough estimate of NPTS IF (MAXPTS .LT. (NX+1)*(NY+1)) THEN INIDOM = .FALSE. MAXPTS = (NX+1)*(NY+1) RETURN ELSE INIDOM = .TRUE. ENDIF NROWS = NY+1 XL = 0.0 YL = 0.0 XR = 1.0 YU = 1.0 DX = (XR-XL)/NX DY = (YU-YL)/NY NX1 = NX/2 NY1 = NINT(NY*0.4) NY2 = NINT(NY*0.6) C C Make grid structure LROW(0) = NROWS IPT = 1 DO 10 I = 0, NY1 LROW(I+1) = IPT IROW(I+1) = I DO 20 J = 0, NX ICOL(IPT) = J IPT = IPT + 1 20 CONTINUE 10 CONTINUE DO 30 I = NY1+1, NY2-1 LROW(I+1) = IPT IROW(I+1) = I DO 40 J = NX/2, NX ICOL(IPT) = J IPT = IPT + 1 40 CONTINUE 30 CONTINUE DO 50 I = NY2, NY LROW(I+1) = IPT IROW(I+1) = I DO 60 J = 0, NX ICOL(IPT) = J IPT = IPT + 1 60 CONTINUE 50 CONTINUE LROW(NROWS+1) = IPT NPTS = IPT-1 C C Boundaries NBNDS = 16 ILBND(1) = 1 ILBND(2) = 2 ILBND(3) = 3 ILBND(4) = 2 ILBND(5) = 1 ILBND(6) = 2 ILBND(7) = 3 ILBND(8) = 4 ILBND( 9) = 12 ILBND(10) = 23 ILBND(11) = 32 ILBND(12) = 21 ILBND(13) = 12 ILBND(14) = 23 ILBND(15) = 34 ILBND(16) = 41 LLBND(0) = NBNDS LLBND(1) = 1 LLBND(2) = LLBND(1) + (NX-1) LLBND(3) = LLBND(2) + (NY1-1) LLBND(4) = LLBND(3) + (NX1-1) LLBND(5) = LLBND(4) + (NY2-NY1-1) LLBND(6) = LLBND(5) + (NX1-1) LLBND(7) = LLBND(6) + (NY1-1) LLBND(8) = LLBND(7) + (NX-1) LLBND( 9) = LLBND( 8) + (NY-1) LLBND(10) = LLBND( 9) + 1 LLBND(11) = LLBND(10) + 1 LLBND(12) = LLBND(11) + 1 LLBND(13) = LLBND(12) + 1 LLBND(14) = LLBND(13) + 1 LLBND(15) = LLBND(14) + 1 LLBND(16) = LLBND(15) + 1 LLBND(17) = LLBND(16) + 1 C Lower and upper boundary pointers DO 100 J = 1, NX-1 LBND(LLBND(1)+J-1) = J + 1 LBND(LLBND(7)+J-1) = NPTS - J 100 CONTINUE C Left boundary pointers DO 120 I = 1, NY1-1 LBND(LLBND(2)+I-1) = I*(NX+1) + 1 LBND(LLBND(6)+I-1) = NPTS - (I+1)*(NX+1) + 1 120 CONTINUE DO 130 I = 1, NY2-NY1-1 LBND(LLBND(4)+I-1) = NY1*(NX+1) + (I+1)*(NX1+1) 130 CONTINUE DO 140 I = 1, NX1-1 LBND(LLBND(3)+I-1) = NY1*(NX+1) + 1 + I LBND(LLBND(5)+I-1) = NPTS - (NY1+1)*(NX+1) + 1 + I 140 CONTINUE C Right boundary pointers DO 110 I = 1, NY1 LBND(LLBND(8)+I-1) = (I+1)*(NX+1) 110 CONTINUE J = LLBND(8)+NY1-1 DO 113 I = 1, NY2-NY1-1 LBND(J+I) = LBND(J) + I*(NX1+1) 113 CONTINUE J = LLBND(8)+NY2-1 DO 116 I = 0, NY1-1 LBND(J+I) = LBND(J-1) + (I+1)*(NX+1) 116 CONTINUE C Corners LBND(LLBND( 9)) = 1 LBND(LLBND(16)) = NX+1 LBND(LLBND(10)) = NY1*(NX+1) + 1 LBND(LLBND(11)) = LBND(LLBND(10)) + NX1 LBND(LLBND(13)) = (NY1+1)*(NX+1) + (NY2-NY1-1)*(NX1+1) + 1 LBND(LLBND(12)) = LBND(LLBND(13)) + NX1 LBND(LLBND(14)) = NPTS - NX LBND(LLBND(15)) = NPTS RETURN END SUBROUTINE PDEIV (T, X, Y, U, NPTS, NPDE) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER NPTS, NPDE DOUBLE PRECISION T, X(NPTS), Y(NPTS), U(NPTS,NPDE) C Ccc PURPOSE: C Define (initial) solution of PDE. C Ccc PARAMETER DESCRIPTION: C T : IN. Time at which (initial) solution should be given C X : IN. Array of X-coordinates for the gridpoints C Y : IN. Array of Y-coordinates for the gridpoints C U : OUT. Array of PDE component values for the gridpoints. C NPTS : IN. Number of gridpoints C NPDE : IN. # PDE components C C----------------------------------------------------------------------- C INTEGER I C DOUBLE PRECISION N, GAMMA, MU0, RHO0, P0, W0, G, DM, KAPPA, AL, + AT, CF, TKAPPA, LT, LL, CS, RHOS, T0, ALPHA, BETA, QC, TC COMMON /PROBLM/ N, GAMMA, MU0, RHO0, P0, W0, G, DM, KAPPA, AL, AT, + CF, TKAPPA, LT, LL, CS, RHOS, T0, ALPHA, BETA, QC, TC SAVE /PROBLM/ C Ccc Problem parameters N = 0.4 KAPPA = 1.0D-10 G = 9.81 DM = 0.0 AT = 0.002 AL = 0.01 RHO0 = 1.0D+3 P0 = 1.0D+5 BETA = 0.0 GAMMA = LOG(1.2) MU0 = 1.0D-3 W0 = 0.25 QC = 1.0D-4 C Ccc Initial solution DO 10 I = 1, NPTS U(I,1) = P0 + (1.0 - Y(I))*RHO0*G U(I,2) = 0.0 10 CONTINUE RETURN END SUBROUTINE PDEF (T, X, Y, U, UT, UX, UY, UXX, UXY, UYY, RES, + NPTS, NPDE) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER NPTS, NPDE DOUBLE PRECISION T, X(NPTS), Y(NPTS), U(NPTS,NPDE), + UT(NPTS,NPDE), UX(NPTS,NPDE), UY(NPTS,NPDE), + UXX(NPTS,NPDE), UXY(NPTS,NPDE), UYY(NPTS,NPDE), + RES(NPTS,NPDE) C Ccc PURPOSE: C Define residual of PDE on interior of domain. Boundary values will be C overwritten later on. C Ccc PARAMETER DESCRIPTION: C T : IN. Time at which residual should be evaluated C X : IN. Array of X-coordinates for the gridpoints C Y : IN. Array of Y-coordinates for the gridpoints C U : IN. Array of PDE components for the gridpoints. C UT : IN. Array of time derivative of PDE components C UX : IN. -I C UY : IN. I C UXX : IN. I Arrays containing space derivatives of PDE components C UXY : IN. I C UYY : IN. -I C RES : OUT. Array containg PDE residual at gridpoints in interior of C domain. The residual values at boundary points will be C overwritten by a call to PDEBC. C NPTS : IN. Number of gridpoints C NPDE : IN. Number of PDE components C C C Ccc INCLUDE 'CMNCMMACH' C C CMNCMMACH C C COMMON with `machine numbers' C LUNOUT : Logical unit # standard output -I C LUNERR : Logical unit # standard error I Set in the routine C UROUND : Smallest machine number such that I MACNUM C 1.0+UROUND > 1.0 and 1.0-UROUND < 1.0 I C XMIN : Smallest floating-point number -I INTEGER LUNOUT, LUNERR DOUBLE PRECISION UROUND, XMIN COMMON /IMACH/ LUNOUT, LUNERR COMMON /RMACH/ UROUND, XMIN SAVE /IMACH/, /RMACH/ C C end INCLUDE 'CMNCMMACH' C C----------------------------------------------------------------------- C DOUBLE PRECISION N, GAMMA, MU0, RHO0, P0, W0, G, DM, KAPPA, AL, + AT, CF, TKAPPA, LT, LL, CS, RHOS, T0, ALPHA, BETA, QC, TC COMMON /PROBLM/ N, GAMMA, MU0, RHO0, P0, W0, G, DM, KAPPA, AL, AT, + CF, TKAPPA, LT, LL, CS, RHOS, T0, ALPHA, BETA, QC, TC SAVE /PROBLM/ C INTEGER I DOUBLE PRECISION P, PT, PX, PY, W, WT, WX, WY, + RHO, RHOX, RHOY, + MU, MUX, MUY, KAPMU, KAPMU2, KAPMUX, KAPMUY, Q1, Q2, QL, + ND11, ND12, ND22, + PXX, PXY, PYY, WXX, WXY, WYY, + ND11Q1, ND11Q2, ND12Q1, ND12Q2, ND22Q1, ND22Q2, + Q1X, Q1Y, Q2X, Q2Y, + ND11X, ND12X, ND12Y, ND22Y, JW1, JW2, JW1X, JW2Y C DO 10 I = 1, NPTS P = U(I,1) PT = UT(I,1) PX = UX(I,1) PY = UY(I,1) W = U(I,2) WT = UT(I,2) WX = UX(I,2) WY = UY(I,2) RHO = RHO0*EXP(BETA*(P-P0)+GAMMA*W) RHOX = RHO*(BETA*PX+GAMMA*WX) RHOY = RHO*(BETA*PY+GAMMA*WY) MU = MU0*(1+1.85*W-4.0*W*W) MUX = MU0*(1.85*WX-8.0*W*WX) MUY = MU0*(1.85*WY-8.0*W*WY) KAPMU = KAPPA/MU KAPMU2 = -KAPMU/MU KAPMUX = KAPMU2*MUX KAPMUY = KAPMU2*MUY Q1 = -KAPMU*PX Q2 = -KAPMU*(PY+RHO*G) QL = MAX(SQRT(Q1*Q1+Q2*Q2),UROUND) ND11 = N*DM + AT*QL + (AL-AT)*Q1*Q1/QL ND12 = (AL-AT)*Q1*Q2/QL ND22 = N*DM + AT*QL + (AL-AT)*Q2*Q2/QL PXX = UXX(I,1) PXY = UXY(I,1) PYY = UYY(I,1) WXX = UXX(I,2) WXY = UXY(I,2) WYY = UYY(I,2) ND11Q1 = (AT + (AL-AT)*(2-(Q1/QL)**2))*Q1/QL ND11Q2 = (AT - (AL-AT)*((Q1/QL)**2))*Q2/QL ND12Q1 = (AL-AT)*(Q2/QL)**3 ND12Q2 = (AL-AT)*(Q1/QL)**3 ND22Q1 = (AT - (AL-AT)*((Q2/QL)**2))*Q1/QL ND22Q2 = (AT + (AL-AT)*(2-(Q2/QL)**2))*Q2/QL Q1X = -(KAPMUX*PX+KAPMU*PXX) Q1Y = -(KAPMUY*PX+KAPMU*PXY) Q2X = -(KAPMUX*(PY+RHO*G)+KAPMU*(PXY+RHOX*G)) Q2Y = -(KAPMUY*(PY+RHO*G)+KAPMU*(PYY+RHOY*G)) ND11X = ND11Q1*Q1X + ND11Q2*Q2X ND12X = ND12Q1*Q1X + ND12Q2*Q2X ND12Y = ND12Q1*Q1Y + ND12Q2*Q2Y ND22Y = ND22Q1*Q1Y + ND22Q2*Q2Y JW1 = -(ND11*WX + ND12*WY) JW2 = -(ND12*WX + ND22*WY) JW1X = -(ND11X*WX+ND11*WXX + ND12X*WY+ND12*WXY) JW2Y = -(ND12Y*WX+ND12*WXY + ND22Y*WY+ND22*WYY) C RES(I,1) = N*RHO*(BETA*PT+GAMMA*WT) + + RHOX*Q1+RHO*Q1X + RHOY*Q2+RHO*Q2Y RES(I,2) = N*RHO*WT + + RHO*Q1*WX + RHO*Q2*WY + + RHOX*JW1+RHO*JW1X + RHOY*JW2+RHO*JW2Y 10 CONTINUE RETURN END SUBROUTINE PDEBC (T, X, Y, U, UT, UX, UY, RES, + NPTS, NPDE, LLBND, ILBND, LBND) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER NPTS, NPDE, LLBND(0:*), ILBND(*), LBND(*) DOUBLE PRECISION T, X(NPTS), Y(NPTS), U(NPTS,NPDE), + UT(NPTS,NPDE), UX(NPTS,NPDE), UY(NPTS,NPDE), + RES(NPTS,NPDE) C Ccc PURPOSE: C Define residual of boundary equations of PDE. The residual on interior C points has already been stored in RES. C Ccc PARAMETER DESCRIPTION: C T : IN. Time at which BC's should be evaluated C X : IN. Array of X-coordinates for the gridpoints C Y : IN. Array of Y-coordinates for the gridpoints C U : IN. Array of PDE components for the gridpoints. C UT : IN. Array of time derivative of PDE components C UX : IN. -I C UY : IN. -I Arrays containing space derivatives of PDE components C RES : INOUT. C IN: PDE residual for interior points (set by PDEF) C OUT: Array with PDE residual at physical boundary points C inserted C NPTS : IN. Number of grid components C NPDE : IN. Number of PDE components C LLBND : IN. (0:LLBND(0)+1) C LLBND(0) = NBNDS: total # physical boundaries in actual grid C LLBND(1:NBNDS): pointers to a specific boundary in LBND C LLBND(NBNDS+1) = NBDPTS+1: total # physical boundary points C in the list + 1 C NB. corners with 2 different types of boundaries should be C pointed at twice. C ILBND : IN. (NBNDS) C ILBND(IB): type of boundary: C 1: Dirichlet C 2: Lower boundary -I C 3: Left boundary I C 4: Upper boundary I max. first order derivative C 5: Right boundary -I C LBND : IN. (NBDPTS) C LBND(LB): pointer to boundary point in actual grid C structure (as in X, Y, and U) C C Ccc INCLUDE 'CMNCMMACH' C C CMNCMMACH C C COMMON with `machine numbers' C LUNOUT : Logical unit # standard output -I C LUNERR : Logical unit # standard error I Set in the routine C UROUND : Smallest machine number such that I MACNUM C 1.0+UROUND > 1.0 and 1.0-UROUND < 1.0 I C XMIN : Smallest floating-point number -I INTEGER LUNOUT, LUNERR DOUBLE PRECISION UROUND, XMIN COMMON /IMACH/ LUNOUT, LUNERR COMMON /RMACH/ UROUND, XMIN SAVE /IMACH/, /RMACH/ C C end INCLUDE 'CMNCMMACH' C C C----------------------------------------------------------------------- C C DOUBLE PRECISION N, GAMMA, MU0, RHO0, P0, W0, G, DM, KAPPA, AL, + AT, CF, TKAPPA, LT, LL, CS, RHOS, T0, ALPHA, BETA, QC, TC COMMON /PROBLM/ N, GAMMA, MU0, RHO0, P0, W0, G, DM, KAPPA, AL, AT, + CF, TKAPPA, LT, LL, CS, RHOS, T0, ALPHA, BETA, QC, TC SAVE /PROBLM/ C INTEGER I, J, K DOUBLE PRECISION P, PY, W, RHO, MU, KAPMU, Q2 C J = 1 C C yL boundary: q2 = 0 C dw/dy = 0 0 'plot.m' % Plots solution and grid levels from data files sol.dat and grid.dat % generated by WRTUNI.f % NB. pcolor with default shading colors a cell with the lowerleft value and % ignores the last row and column. % nxb=input('nX base grid? '); nyb=input('nY base grid? '); load sol.dat load grid.dat [n,npde]=size(sol); unilev=floor(log(n/(nxb*nyb))/(log(2)*2)+1) nx=nxb*2^(unilev-1); ny=nyb*2^(unilev-1); for ic=1:npde Umin=input(['min. sol. value comp. ',int2str(ic),'? ']) Umax=input(['max. sol. value comp. ',int2str(ic),'? ']) U=zeros(ny+1,nx+1); for j = 0:ny for i = 0:nx U(j+1,i+1) = sol(j*(nx+1)+i+1,ic); end end figure(ic); colormap(jet); pcolor(U); shading('interp'); if Umin < Umax caxis([Umin Umax]) end; keyboard end; G=zeros(ny+2,nx+2); for j = 0:ny for i = 0:nx G(j+1,i+1) = grid(j*(nx+1)+i+1); end end figure(npde+1); grc=[1 1 1; 1 1 0; 0 1 0; 0 1 1; 0 0 1; 1 0 0; 1 0 1]; colormap(grc(1:unilev,:)); pcolor(G); caxis([0.99 unilev]) SHAR_EOF fi # end of overwriting check if test -f 'wrtuni.f' then echo shar: will not over-write existing file "'wrtuni.f'" else cat << \SHAR_EOF > 'wrtuni.f' PROGRAM WRTUNI C C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C !!! !!! C !!! In subroutine WRUNI the constant NONVAL should be adjusted to !!! C !!! the data (NONVAL = impossible value for the first componenent) !!! C !!! !!! C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C C----------------------------------------------------------------------- C Ccc This program reads a file made by subroutine DUMP and writes the C (interpolated) solution on a uniform grid of a specified grid level C to the output file sol.dat. The maximum grid level used in each point C is written to the file grid.dat. C NB. This program is not correct for a domain with holes in it with C a size of the width of the base grid, e.g. it will ignore some holes C in the domain of the example problem. C Ccc EXTERNALS USED: EXTERNAL WRUNI, RDDUMP C C Ccc INCLUDE 'CMNWRITEF' C C CMNWRITEF C C COMMON needed for continuation calls INTEGER MAXLVW, NPDEW, LRWKPS, LIWKPS, LRWKB, LIWKB LOGICAL FIRST, SECOND DOUBLE PRECISION T0, TW, TEW, DTW, XLW, YLW, XRW, YUW, DXB, DYB, + DTO COMMON /WRITIF/ MAXLVW, NPDEW, LRWKPS, LIWKPS, LRWKB, LIWKB COMMON /WRITLF/ FIRST, SECOND COMMON /WRITRF/ T0, TW, TEW, DTW, XLW,YLW, XRW,YUW, DXB, DYB, DTO SAVE /WRITIF/, /WRITLF/, /WRITRF/ C C end INCLUDE 'CMNWRITEF' C C C----------------------------------------------------------------------- C INTEGER MXLEV, NPD, NPTS, LENIWK, LENRWK PARAMETER (MXLEV=5, NPD=3, NPTS=10000) PARAMETER (LENIWK=NPTS*(7*MXLEV+20), + LENRWK=5*NPTS*NPD*MXLEV) C CHARACTER FILE*128 INTEGER IWK(LENIWK), + LSGNM1, LSGN, LSGNP1, LSUNM1, LSSN, LSUN, + LUNI, MAXLEV, NX, NXB, NY, NYB, UNILEV DOUBLE PRECISION RWK(LENRWK) PRINT *, 'DUMP file?' READ '(A)', FILE C OPEN(UNIT=62,FILE=FILE,FORM='UNFORMATTED') CALL RDDUMP (62, RWK, LENRWK, IWK, LENIWK) CLOSE(62) C C Setup work storage LSGNM1 = 1 LSGN = LSGNM1 + MAXLVW+1 LSGNP1 = LSGN + MAXLVW+1 LSUNM1 = LSGNP1 + MAXLVW+1 LSSN = LSUNM1 + MAXLVW LSUN = LSSN + MAXLVW C C Check workspace MAXLEV = IWK(LSGN) PRINT *, 'Max. grid level?' READ *, UNILEV UNILEV = MIN(UNILEV,MAXLEV) NXB = NINT((XRW - XLW)/DXB) NYB = NINT((YUW - YLW)/DYB) NX = NXB * 2**(UNILEV-1) NY = NYB * 2**(UNILEV-1) LUNI = LENRWK - (NX+1)*(NY+1)*NPDEW IF (LUNI .LT. IWK(LSUN+MAXLVW)) STOP 'workspace' C C Write problem info to standard output and write the interpolated C solution and grid levels to the files PRINT *, 'T, NPDE, XL, YL, DXB, DYB, NXB, NYB' PRINT *, TW, NPDEW, XLW, YLW, DXB, DYB, NXB, NYB FILE = 'sol.dat' OPEN(UNIT=61,FILE=FILE) FILE = 'grid.dat' OPEN(UNIT=63,FILE=FILE) CALL WRUNI (61, 63, UNILEV, + TW, NPDEW, XLW, YLW, DXB, DYB, NXB, NYB, + IWK(LSGN), IWK(LIWKPS), IWK(LSUN), RWK(LRWKPS), + RWK(LUNI), NX, NY) CLOSE(61) CLOSE(63) END SUBROUTINE RDDUMP (LUNDMP, RWK, LENRWK, IWK, LENIWK) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER LENIWK INTEGER LUNDMP, LENRWK, IWK(LENIWK) DOUBLE PRECISION RWK(LENRWK) C Ccc PURPOSE: C Read all information necessary for a restart of VLUGR2 from file C Ccc PARAMETER DESCRIPTION: C LUNDMP : IN. Logical unit number of dumpfile. Should be opened as an C unformatted file. C RWK : OUT. Real workstorage intended to pass to VLUGR2 C LENRWK : IN. Dimension of RWK. C IWK : OUT. Integer workstorage intended to pass to VLUGR2 C LENIWK : IN. Dimension of IWK. C Ccc EXTERNALS USED: NONE C C Ccc INCLUDE 'CMNSTATS' C C CMNSTATS C C COMMON with integration statistics INTEGER MXCLEV, MXCNIT PARAMETER (MXCLEV = 10, MXCNIT = 20) INTEGER LUNPDS, LUNNLS, LUNLSS, LEVEL, NSTEPS, NREJS, + NJACS(MXCLEV), NRESID(MXCLEV), NNIT(MXCLEV), + NLSIT(MXCLEV,MXCNIT) COMMON /STATS/ LUNPDS, LUNNLS, LUNLSS, LEVEL, NSTEPS, NREJS, + NJACS, NRESID, NNIT, NLSIT SAVE /STATS/ C C end INCLUDE 'CMNSTATS' C C Ccc INCLUDE 'CMNWRITEF' C C CMNWRITEF C C COMMON needed for continuation calls INTEGER MAXLVW, NPDEW, LRWKPS, LIWKPS, LRWKB, LIWKB LOGICAL FIRST, SECOND DOUBLE PRECISION T0, TW, TEW, DTW, XLW, YLW, XRW, YUW, DXB, DYB, + DTO COMMON /WRITIF/ MAXLVW, NPDEW, LRWKPS, LIWKPS, LRWKB, LIWKB COMMON /WRITLF/ FIRST, SECOND COMMON /WRITRF/ T0, TW, TEW, DTW, XLW,YLW, XRW,YUW, DXB, DYB, DTO SAVE /WRITIF/, /WRITLF/, /WRITRF/ C C end INCLUDE 'CMNWRITEF' C C C----------------------------------------------------------------------- C INTEGER I, J READ(LUNDMP) MAXLVW, NPDEW, LRWKPS, LIWKPS, LRWKB, LIWKB, + FIRST, SECOND, + T0, TW, TEW, DTW, XLW, YLW, XRW, YUW, DXB, DYB, DTO IF (LENRWK .LT. LRWKPS+LRWKB .OR. LENIWK .LT. LIWKPS+LIWKB) THEN PRINT *, LENRWK, LRWKPS+LRWKB, LENIWK, LIWKPS+LIWKB STOP 'work space too small' ENDIF READ(LUNDMP) LUNPDS, LUNNLS, LUNLSS, LEVEL, NSTEPS, NREJS, + (NJACS(I), I=1,MXCLEV), (NRESID(I), I=1,MXCLEV), + (NNIT(I), I=1,MXCLEV), ((NLSIT(I,J), I=1,MXCLEV), J=1,MXCNIT) READ(LUNDMP) (RWK(I), I=1,LRWKPS+LRWKB) READ(LUNDMP) (IWK(I), I=1,LIWKPS+LIWKB) C RETURN END SUBROUTINE WRUNI (LUNS, LUNG, UNILEV, + T, NPDE, XL, YL, DXB, DYB, NXB, NYB, + LGRID, ISTRUC, LSOL, SOL, UNIFRM, NX, NY) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER LUNS, LUNG, UNILEV, + NPDE, NXB, NYB, LGRID(0:*), ISTRUC(*), LSOL(*), NX, NY DOUBLE PRECISION T, XL, YL, DXB, DYB, SOL(*), + UNIFRM(0:NX,0:NY,NPDE) C Ccc PURPOSE: C Write (interpolated) solution values at grid level UNILEV to file C LUNS. C Write maximum gridlevel used in each point to file LUNG. C NB. The data will not be correct for a domain with holes in it with C a size of the width of the base grid, e.g. it will ignore some holes C in the domain of the example problem. C Ccc PARAMETER DESCRIPTION: C LUNS : IN. Logical unit number of solution file C LUNG : IN. Logical unit number of grid level file C UNILEV : IN. Maximum grid level to be used to generate the data C T : IN. Value of time variable C NPDE : IN. # PDE components C XL : IN. X-coordinate of lower left corner of (virtual) domain C YL : IN. Y-coordinate of lower left corner of (virtual) domain C DXB : IN. Cell width in X-direction of base grid C DYB : IN. Cell width in Y-direction of base grid C NXB,NYB: IN. # gridcells in X- and Y-direction, resp., on base grid C LGRID : IN. (0:*) C LGRID(0) = max. grid level used at T C LGRID(1): pointer to base grid structure ISTRUC C LGRID(LEVEL): pointer to grid structure (LROW, IROW, ICOL) C of refinement level LEVEL for time T C ISTRUC : IN. (*) C ISTRUC(LGRID(LEVEL):.) contains (LROW,IROW,ICOL) of grid C level LEVEL, C LROW : (0:LROW(0)+1) C LROW(0) = NROWS: Actual # rows in grid C LROW(1:NROWS): pointers to the start of a row in the grid C LROW(NROWS+1) = NPTS+1: Actual # nodes in grid + 1 C IROW : (NROWS) C IROW(IR): row number of row IR in virtual rectangle C ICOL : (NPTS) C ICOL(IPT): column number of grid point IPT in virtual C rectangle C LSOL : IN. (*) C LSOL(LEVEL): pointer to (injected) solution at grid C of refinement level LEVEL for time T C SOL : IN. (*) C SOL(LSOL(LEVEL)+1:LSOL(LEVEL)+NPTS(LEVEL)*NPDE) contains C U_LEVEL(NPTS,NPDE) C UNIFRM : WORK. (Interpolated) solution on level UNILEV / max. grid C level used. C NX, NY : IN. # gridcells in X- and Y-direction, resp., on grid of C of level UNILEV C C----------------------------------------------------------------------- C C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C !!! !!! C !!! In subroutine WRUNI the constant NONVAL should be adjusted to !!! C !!! the data (NONVAL = impossible value for the first componenent) !!! C !!! !!! DOUBLE PRECISION NONVAL PARAMETER (NONVAL = -999.999) C !!! !!! C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C C----------------------------------------------------------------------- C INTEGER I, IC, ICOL, IMUL, IPT, IR, IROW, J, + LEVEL, LLROW, LIROW, LICOL, MAXLEV, NROWS, NPTS DO 1 IC = 1, NPDE DO 1 IROW = 0, NY DO 1 ICOL = 0, NX UNIFRM(ICOL,IROW,IC) = NONVAL 1 CONTINUE MAXLEV = LGRID(0) DO 10 LEVEL = 1, UNILEV IMUL = 2**(UNILEV-LEVEL) LLROW = LGRID(LEVEL) NROWS = ISTRUC(LLROW) NPTS = ISTRUC(LLROW+NROWS+1)-1 LIROW = LLROW+NROWS+2 LICOL = LIROW+NROWS DO 20 IR= 1, NROWS IROW = ISTRUC(LIROW-1+IR)*IMUL DO 30 IPT = ISTRUC(LLROW+IR), ISTRUC(LLROW+IR+1)-1 ICOL = ISTRUC(LICOL-1+IPT)*IMUL DO 40 IC = 1, NPDE UNIFRM(ICOL,IROW,IC) = + SOL(LSOL(LEVEL)+(IC-1)*NPTS+IPT) 40 CONTINUE 30 CONTINUE 20 CONTINUE 10 CONTINUE DO 100 LEVEL = 2, UNILEV IMUL = 2**(UNILEV-LEVEL) DO 110 J = IMUL, NY, IMUL*2 DO 110 I = 0, NX, IMUL*2 IF (UNIFRM(I,J,1) .EQ. NONVAL) THEN DO 120 IC = 1, NPDE UNIFRM(I,J,IC) = + (UNIFRM(I,J-IMUL,IC)+UNIFRM(I,J+IMUL,IC))/2 120 CONTINUE ENDIF 110 CONTINUE DO 130 J = 0, NY, IMUL DO 130 I = IMUL, NX, IMUL*2 IF (UNIFRM(I,J,1) .EQ. NONVAL) THEN DO 140 IC = 1, NPDE UNIFRM(I,J,IC) = + (UNIFRM(I-IMUL,J,IC)+UNIFRM(I+IMUL,J,IC))/2 140 CONTINUE ENDIF 130 CONTINUE 100 CONTINUE DO 150 J = 0, NY DO 150 I = 0, NX WRITE(LUNS,'(100E13.3)') (UNIFRM(I,J,IC), IC = 1, NPDE) 150 CONTINUE C C Grids DO 201 IROW = 0, NY DO 201 ICOL = 0, NX UNIFRM(ICOL,IROW,1) = 0 201 CONTINUE DO 210 LEVEL = 1, UNILEV IMUL = 2**(UNILEV-LEVEL) LLROW = LGRID(LEVEL) NROWS = ISTRUC(LLROW) NPTS = ISTRUC(LLROW+NROWS+1)-1 LIROW = LLROW+NROWS+2 LICOL = LIROW+NROWS DO 220 IR= 1, NROWS IROW = ISTRUC(LIROW-1+IR)*IMUL DO 230 IPT = ISTRUC(LLROW+IR), ISTRUC(LLROW+IR+1)-1 ICOL = ISTRUC(LICOL-1+IPT)*IMUL UNIFRM(ICOL,IROW,1) = LEVEL 230 CONTINUE 220 CONTINUE 210 CONTINUE DO 300 LEVEL = 2, UNILEV IMUL = 2**(UNILEV-LEVEL) DO 310 J = IMUL, NY, IMUL*2 DO 310 I = 0, NX, IMUL*2 IF (UNIFRM(I,J,1) .LT. LEVEL) THEN UNIFRM(I,J,1) = + MIN(UNIFRM(I,J-IMUL,1),UNIFRM(I,J+IMUL,1)) ENDIF 310 CONTINUE DO 330 J = 0, NY, IMUL DO 330 I = IMUL, NX, IMUL*2 IF (UNIFRM(I,J,1) .LT. LEVEL) THEN UNIFRM(I,J,1) = + MIN(UNIFRM(I-IMUL,J,1),UNIFRM(I+IMUL,J,1)) ENDIF 330 CONTINUE 300 CONTINUE DO 350 J = 0, NY DO 350 I = 0, NX WRITE(LUNG,'(I2)') NINT(UNIFRM(I,J,1)) 350 CONTINUE RETURN END SHAR_EOF fi # end of overwriting check if test -f 'exmpl.f' then echo shar: will not over-write existing file "'exmpl.f'" else cat << \SHAR_EOF > 'exmpl.f' PROGRAM EXMPL C C Ccc INCLUDE 'PARNEWTON' C C PARNEWTON C C Parameters for Newton process C MAXNIT : Max. number of Newton iterations C MAXJAC : Max. number of Jacobian / preconditioner evaluations during C a Newton process C TOLNEW : Tolerance for Newton process: C rho/(1-rho)*|| corr.||_w < TOLNEW INTEGER MAXNIT, MAXJAC DOUBLE PRECISION TOLNEW PARAMETER (MAXNIT = 10, MAXJAC = 2, TOLNEW = 1.0) C C end INCLUDE 'PARNEWTON' C C Ccc INCLUDE 'PARGCRO' C C PARGCRO C C Parameters for linear system solver GCRO + (block-)diagonal C preconditioner C IDIAGP : 0: block-diagonal + first order derivatives C 1: block-diagonal neglecting first order derivatives C 2: diagonal + first order derivatives C 3: diagonal neglecting first order derivatives C NRRMAX : Max. number of restarts of outer loop C MAXLR : Max. number of iterations in outer loop C MAXL : Max. number of iterations in GMRES inner loop C TOLLSC : Tolerance for linear system solver INTEGER IDIAGP, NRRMAX, MAXLR, MAXL DOUBLE PRECISION TOLLSC PARAMETER (NRRMAX = 1, MAXLR = 5, MAXL = 20) C PARAMETER (NRRMAX = 1, MAXLR = 3, MAXL = 10) PARAMETER (TOLLSC = TOLNEW/10) COMMON /IGCRO/ IDIAGP SAVE /IGCRO/ C C end INCLUDE 'PARGCRO' C INTEGER MXLEV, NPD, NPTS, LENIWK, LENRWK, LENLWK PARAMETER (MXLEV=2, NPD=2, NPTS=2500) PARAMETER (LENIWK=NPTS*(5*MXLEV+11), + LENRWK=NPTS*NPD*(5*MXLEV+9 + + (9*NPD+2*MAXLR+MAXL+7)), + LENLWK=2*NPTS) C C----------------------------------------------------------------------- C INTEGER LUNDMP PARAMETER (LUNDMP = 89) C INTEGER NPDE, INFO(7), IWK(LENIWK), MNTR LOGICAL LWK(LENLWK) DOUBLE PRECISION T, TOUT, DT, XL, YL, XR, YU, DX, DY, + TOLS, TOLT, RINFO(2+3*NPD), RWK(LENRWK) C First call of VLUGR2 MNTR = 0 NPDE = 2 T = 0.0 TOUT = 1.0 DT = 0.001 C Since domain is not a rectangle the grid parameters need not to be C specified here (cf. INIDOM) TOLS = 0.1 TOLT = 0.05 INFO(1) = 1 C MAXLEV INFO(2) = 5 C Domain not a rectangle INFO(3) = 1 C Linear system solver: GCRO + Diagonal scaling C (no first order derivatives at the boundaries) INFO(4) = 13 OPEN (UNIT=61,FILE='RunInfo') C Write integration history to unit # 61 INFO(5) = 61 C Write Newton info to unit # 61 INFO(6) = 61 C Write GCRO info to unit # 61 INFO(7) = 61 C DTMIN = 1D-7 RINFO(1) = 1.0D-7 C DTMAX = 1.0 RINFO(2) = 1.0 C UMAX = 1.0 RINFO(3) = 1.0 RINFO(4) = 1.0 C SPCWGT = 1.0 RINFO(5) = 1.0 RINFO(6) = 1.0 C TIMWGT = 1.0 RINFO(7) = 1.0 RINFO(8) = 1.0 C C Call main routine CALL VLUGR2 (NPDE, T, TOUT, DT, XL, YL, XR, YU, DX, DY, + TOLS, TOLT, INFO, RINFO, RWK, LENRWK, IWK, LENIWK, LWK, LENLWK, + MNTR) PRINT *, 'VLUGR2 returned with MNTR=', MNTR C C Save info on file OPEN(UNIT=LUNDMP,FILE='DUMP',FORM='UNFORMATTED') CALL DUMP (LUNDMP, RWK, IWK) CLOSE(LUNDMP) END LOGICAL FUNCTION INIDOM (MAXPTS, XL, YL, XR, YU, DX, DY, + LROW, IROW, ICOL, LLBND, ILBND, LBND) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER MAXPTS, LROW(0:*), IROW(*), ICOL(*), + LLBND(0:*), ILBND(*), LBND(*) DOUBLE PRECISION XL, YL, XR, YU, DX, DY C Ccc PURPOSE: C Define initial domain. NB. Boundaries should consist of as many points C as are necessary to employ second order space discretization, i.e., C a boundary enclosing the internal part of the domain should not C include less than 3 grid points including the corners. If Neumann C boundaries are used the minimum is 4 since otherwise the Jacobian C matrix will be singular. C C A (virtual) rectangle is placed upon the (irregular) domain. The C lowerleft point of this rectangle is (XL,YL) in physical coordinates C and (0,0) in column, resp. row coordinates. The upperright point is C (XR,YU) resp. (Nx, Ny), where Nx = (XR-XL)/DX and Ny = (YU-YL)/DY. C Only real grid points are stored. C The coordinate values of the initial grid should be stored rowwise C in LROW, IROW, ICOL. C Pointers to the boundary points should be stored in a list together C with the type of the boundary. (LLBND, ILBND, LBND) C C On exit INIDOM = .FALSE. if the # grid points required is larger C than MAXPTS and MAXPTS is set to the required # points. C Ccc PARAMETER DESCRIPTION: C MAXPTS : INOUT. C IN: Max. # grid points allowed by the available workspace C OUT: # grid points required, if larger than # points allowed C XL : OUT. X-coordinate of lowerleft point of virtual rectangle C YL : OUT. Y-coordinate of lowerleft point of virtual rectangle C XR : OUT. X-coordinate of upperright point of virtual rectangle C YU : OUT. Y-coordinate of upperright point of virtual rectangle C DX : OUT. Grid width in X-direction C DY : OUT. Grid width in Y-direction C LROW : OUT. INTEGER array of dimension (0:LROW(0)+1) C LROW(0) = NROWS: Actual # rows in grid C LROW(1:NROWS): pointers to the start of a row in the grid C structure C LROW(NROWS+1) = NPTS+1: Actual # nodes in grid + 1 C IROW : OUT. INTEGER array of dimension (NROWS) C IROW(IR): row number of row IR in virtual rectangle C ICOL : OUT. INTEGER array of dimension (NPTS) C ICOL(IPT): column number of grid point IPT in virtual C rectangle C LLBND : (0:LLBND(0)+2) C LLBND(0) = NBNDS: total # physical boundaries and corners in C actual domain. C NB. corners should be stored as an independent boundary C (cf. ILBND). The order in LLBND should be first the C boundaries and then the corners. C LLBND(1:NBNDS): pointers to a specific boundary or corner in C LBND C LLBND(NBNDS+1) = NBDPTS+1: total # physical boundary points C in LBND + 1 C LLBND(NBNDS+1): pointer to internal boundary in LBND C LLBND(NBNDS+2) = NBIPTS+1: total # points in LBND + 1 C ILBND : (NBNDS) C ILBND(IB): type of boundary: C 1: Lower boundary -I C 2: Left boundary I C 3: Upper boundary I max. first order derivative C 4: Right boundary -I C 12: Lowerleft corner -I C 23: Leftupper corner I corners of 90 degrees C 34: Upperright corner I (external corners) C 41: Rightlower corner -I max. first order deriv. C 21: Leftlower corner -I C 32: Upperleft corner I corners of 270 degrees C 43: Rightupper corner I (internal corners) C 14: Lowerright corner -I max. first order deriv. C LBND : (NBIPTS) C LBND(IBPT): pointer to boundary point in actual grid C structure C Ccc EXTERNALS USED: NONE C C----------------------------------------------------------------------- C C Square domain [0,1]x[0,1] with holes. Dirichlet boundaries. C INTEGER NX, NY PARAMETER (NX = 10, NY = 10) INTEGER IDOM((NX+1)*(NY+1)) C INTEGER I, IPT, J, NROWS, NPTS, NBNDS NPTS = (NX+1)*(NY+1) - (NX-2) - 2*3 - 2 IF (MAXPTS .LT. NPTS) THEN INIDOM = .FALSE. MAXPTS = NPTS RETURN ELSE INIDOM = .TRUE. ENDIF NROWS = NY+1 XL = 0.0 YL = 0.0 XR = 1.0 YU = 1.0 DX = (XR-XL)/NX DY = (YU-YL)/NY C C Make grid structure LROW(0) = NROWS IPT = 1 DO 10 I = 0, 0 LROW(I+1) = IPT IROW(I+1) = I DO 15 J = 0, 2 ICOL(IPT) = J IPT = IPT + 1 15 CONTINUE 10 CONTINUE DO 20 I = 1, 3 LROW(I+1) = IPT IROW(I+1) = I DO 23 J = 0, 2 ICOL(IPT) = J IPT = IPT + 1 23 CONTINUE DO 26 J = 3, 10 ICOL(IPT) = J IPT = IPT + 1 26 CONTINUE 20 CONTINUE DO 30 I = 4, 4 LROW(I+1) = IPT IROW(I+1) = I DO 33 J = 0, 2 ICOL(IPT) = J IPT = IPT + 1 33 CONTINUE DO 36 J = 3, 5 ICOL(IPT) = J IPT = IPT + 1 36 CONTINUE DO 39 J = 8, 10 ICOL(IPT) = J IPT = IPT + 1 39 CONTINUE 30 CONTINUE DO 40 I = 5, 7 LROW(I+1) = IPT IROW(I+1) = I DO 43 J = 0, 2 ICOL(IPT) = J IPT = IPT + 1 43 CONTINUE DO 46 J = 3, 10 ICOL(IPT) = J IPT = IPT + 1 46 CONTINUE 40 CONTINUE DO 50 I = 8, 10 LROW(I+1) = IPT IROW(I+1) = I DO 56 J = 0, 8 ICOL(IPT) = J IPT = IPT + 1 56 CONTINUE 50 CONTINUE LROW(NROWS+1) = IPT C C Boundaries NBNDS = 28 ILBND(1) = 1 LLBND(1) = 1 IPT = 2 LBND(LLBND(1)) = IPT ILBND(2) = 2 LLBND(2) = LLBND(1) + 1 IPT = 4 LBND(LLBND(2) ) = IPT IPT = 15 LBND(LLBND(2)+1) = IPT IPT = 26 LBND(LLBND(2)+2) = IPT IPT = 37 LBND(LLBND(2)+3) = IPT IPT = 46 LBND(LLBND(2)+4) = IPT IPT = 57 LBND(LLBND(2)+5) = IPT IPT = 68 LBND(LLBND(2)+6) = IPT IPT = 79 LBND(LLBND(2)+7) = IPT IPT = 88 LBND(LLBND(2)+8) = IPT ILBND(3) = 3 LLBND(3) = LLBND(2) + 9 DO 130 J = 0, 6 IPT = 98+J LBND(LLBND(3)+J) = IPT 130 CONTINUE ILBND(4) = 4 LLBND(4) = LLBND(3) + 7 IPT = 96 LBND(LLBND(4)) = IPT ILBND(5) = 1 LLBND(5) = LLBND(4) + 1 DO 150 J = 0, 4 IPT = 86-J LBND(LLBND(5)+J) = IPT 150 CONTINUE ILBND(6) = 4 LLBND(6) = LLBND(5) + 5 DO 160 J = 0, 6 IPT = LBND(LLBND(2)+J) + 2 LBND(LLBND(6)+J) = IPT 160 CONTINUE ILBND(7) = 1 LLBND(7) = LLBND(6) + 7 DO 170 J = 0, 5 IPT = 8+J LBND(LLBND(7)+J) = IPT 170 CONTINUE ILBND(8) = 2 LLBND(8) = LLBND(7) + 6 DO 180 J = 0, 4 IPT = LBND(LLBND(6)+J+1) + 1 LBND(LLBND(8)+J) = IPT 180 CONTINUE ILBND(9) = 3 LLBND(9) = LLBND(8) + 5 DO 190 J = 0, 5 IPT = 72+J LBND(LLBND(9)+J) = IPT 190 CONTINUE ILBND(10) = 4 LLBND(10) = LLBND(9) + 6 IPT = 67 LBND(LLBND(10) ) = IPT IPT = 56 LBND(LLBND(10)+1) = IPT IPT = 45 LBND(LLBND(10)+2) = IPT IPT = 36 LBND(LLBND(10)+3) = IPT IPT = 25 LBND(LLBND(10)+4) = IPT ILBND(11) = 1 LLBND(11) = LLBND(10) + 5 IPT = 52 LBND(LLBND(11) ) = IPT IPT = 53 LBND(LLBND(11)+1) = IPT ILBND(12) = 2 LLBND(12) = LLBND(11) + 2 IPT = 43 LBND(LLBND(12) ) = IPT ILBND(13) = 3 LLBND(13) = LLBND(12) + 1 IPT = 32 LBND(LLBND(13) ) = IPT IPT = 33 LBND(LLBND(13)+1) = IPT ILBND(14) = 4 LLBND(14) = LLBND(13) + 2 IPT = 42 LBND(LLBND(14) ) = IPT C ILBND(15) = 12 LLBND(15) = LLBND(14) + 1 IPT = 1 LBND(LLBND(15)) = IPT ILBND(16) = 23 LLBND(16) = LLBND(15) + 1 IPT = 97 LBND(LLBND(16)) = IPT ILBND(17) = 34 LLBND(17) = LLBND(16) + 1 IPT = 105 LBND(LLBND(17)) = IPT ILBND(18) = 41 LLBND(18) = LLBND(17) + 1 IPT = 87 LBND(LLBND(18)) = IPT ILBND(19) = 14 LLBND(19) = LLBND(18) + 1 IPT = 81 LBND(LLBND(19)) = IPT ILBND(20) = 41 LLBND(20) = LLBND(19) + 1 IPT = 3 LBND(LLBND(20)) = IPT ILBND(21) = 12 LLBND(21) = LLBND(20) + 1 IPT = 7 LBND(LLBND(21)) = IPT ILBND(22) = 23 LLBND(22) = LLBND(21) + 1 IPT = 71 LBND(LLBND(22)) = IPT ILBND(23) = 34 LLBND(23) = LLBND(22) + 1 IPT = 78 LBND(LLBND(23)) = IPT ILBND(24) = 41 LLBND(24) = LLBND(23) + 1 IPT = 14 LBND(LLBND(24)) = IPT ILBND(25) = 14 LLBND(25) = LLBND(24) + 1 IPT = 51 LBND(LLBND(25)) = IPT ILBND(26) = 43 LLBND(26) = LLBND(25) + 1 IPT = 31 LBND(LLBND(26)) = IPT ILBND(27) = 32 LLBND(27) = LLBND(26) + 1 IPT = 34 LBND(LLBND(27)) = IPT ILBND(28) = 21 LLBND(28) = LLBND(27) + 1 IPT = 54 LBND(LLBND(28)) = IPT C LLBND(29) = LLBND(28) + 1 LLBND( 0) = NBNDS C No internal boundaries C (only necessary because we want to print the domain) LLBND(NBNDS+2) = LLBND(NBNDS+1) PRINT *, 'Input domain:' CALL PRDOM (LROW, IROW, ICOL, LLBND, ILBND, LBND, + IDOM, NX, NY) RETURN END SUBROUTINE CHSPCM (T, LEVEL, NPTS, X, Y, NPDE, U, SPCMON, TOL) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER LEVEL, NPTS, NPDE DOUBLE PRECISION T, X(NPTS), Y(NPTS), U(NPTS,NPDE), SPCMON(NPTS), + TOL C Ccc PURPOSE: C Force grid refinement. C If for a node IPT SPCMON(IPT) > TOL the 16 surrounding cells will be C refined. C Ccc PARAMETER DESCRIPTION: C T : IN. Current value of time variable C LEVEL : IN. Current grid level C NPTS : IN. Number of grid points at this level C X : IN. Array of X-coordinates for the gridpoints C Y : IN. Array of Y-coordinates for the gridpoints C NPDE : IN. Number of PDE components C U : IN. Array of PDE components for the gridpoints C SPCMON : INOUT. C IN: Space monitor values as determined by VLUGR2 C OUT: Changed to a value > TOL where refinement is required C TOL : IN. Tolerance with which SPCMON will be compared C C----------------------------------------------------------------------- C INTEGER I C IF (LEVEL .GE. 3) RETURN DO 10 I = 1, NPTS IF (ABS(X(I)-1.0) .LT. 0.0001 .AND. + ABS(Y(I)-0.1) .LT. 0.0001) THEN SPCMON(I) = 2*TOL ENDIF 10 CONTINUE C RETURN END SUBROUTINE MONITR (T, DT, DTNEW, XL, YL, DXB, DYB, + LGRID, ISTRUC, LSOL, SOL) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER LGRID(0:*), ISTRUC(*), LSOL(*) DOUBLE PRECISION T, DT, DTNEW, XL, YL, DXB, DYB, SOL(*) C Ccc PURPOSE: C Control after a successful time step. The solution can be printed, C plotted or compared with the exact solution. C Ccc PARAMETER DESCRIPTION: C T : IN. Current value of time variable C DT : IN. Current time step size C DTNEW : IN. Time step size for next time step C XL : IN. X-coordinate of lowerleft corner of (virtual) rectangle C YL : IN. Y-coordinate of lowerleft corner of (virtual) rectangle C DXB : IN. Cell width in X-direction of base grid C DYB : IN. Cell width in Y-direction of base grid C LGRID : IN. (0:*) C LGRID(0) = max. grid level used at T C LGRID(1): pointer to base grid structure ISTRUC C LGRID(LEVEL): pointer to grid structure (LROW, IROW, ICOL) C of refinement level LEVEL for time T C ISTRUC : IN. (*) C ISTRUC(LGRID(LEVEL):.) contains (LROW,IROW,ICOL) of grid C level LEVEL, C LROW : (0:LROW(0)+1) C LROW(0) = NROWS: Actual # rows in grid C LROW(1:NROWS): pointers to the start of a row in the grid C LROW(NROWS+1) = NPTS+1: Actual # nodes in grid + 1 C IROW : (NROWS) C IROW(IR): row number of row IR in virtual rectangle C ICOL : (NPTS) C ICOL(IPT): column number of grid point IPT in virtual C rectangle C LSOL : IN. (*) C LSOL(LEVEL): pointer to (injected) solution at grid C of refinement level LEVEL for time T C SOL : IN. (*) C SOL(LSOL(LEVEL)+1:LSOL(LEVEL)+NPTS(LEVEL)*NPDE) contains C U_LEVEL(NPTS,NPDE) C C C Local arrays: INTEGER MAXPTS, NPDE PARAMETER (MAXPTS=10000, NPDE=2) DOUBLE PRECISION X(MAXPTS), Y(MAXPTS), UEX(MAXPTS*NPDE) C C----------------------------------------------------------------------- C INTEGER MAXLEV, LEVEL, LLROW, NROWS, NPTS, LIROW, LICOL DOUBLE PRECISION DX, DY C C Loop over the grid levels from coarse to fine. C Get physical coordinates of grid points C Compute ||err||_max MAXLEV = LGRID(0) DX = DXB DY = DYB DO 10 LEVEL = 1, MAXLEV LLROW = LGRID(LEVEL) NROWS = ISTRUC(LLROW) NPTS = ISTRUC(LLROW+NROWS+1)-1 LIROW = LLROW+NROWS+2 LICOL = LIROW+NROWS CALL SETXY (XL, YL, DX, DY, + ISTRUC(LLROW), ISTRUC(LIROW), ISTRUC(LICOL), X, Y) DX = DX/2 DY = DY/2 CALL PRERR (LEVEL, T, NPTS, NPDE, X, Y, SOL(LSOL(LEVEL)+1), + UEX) 10 CONTINUE RETURN END SUBROUTINE PRERR (LEVEL, T, NPTS, NPDE, X, Y, U, UEX) INTEGER LEVEL, NPTS, NPDE DOUBLE PRECISION T, X(NPTS), Y(NPTS), U(NPTS,NPDE), UEX(NPTS,NPDE) INTEGER I,J DOUBLE PRECISION RMAX(2) CALL PDEIV (T, X, Y, UEX, NPTS, NPDE) RMAX(1) = 0.0 RMAX(2) = 0.0 DO 10 I = 1, NPTS J = 1 RMAX(J) = MAX(RMAX(J),ABS(UEX(I,J)-U(I,J))) J = 2 RMAX(J) = MAX(RMAX(J),ABS(UEX(I,J)-U(I,J))) 10 CONTINUE PRINT '(''Error at T='',E9.3,'', level='',I1,'' :'',2E11.3,I10)', + T, LEVEL, RMAX(1), RMAX(2), NPTS RETURN END SUBROUTINE PDEIV (T, X, Y, U, NPTS, NPDE) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER NPTS, NPDE DOUBLE PRECISION T, X(NPTS), Y(NPTS), U(NPTS,NPDE) C Ccc PURPOSE: C Define (initial) solution of PDE. C Ccc PARAMETER DESCRIPTION: C T : IN. Time at which (initial) solution should be given C X : IN. Array of X-coordinates for the gridpoints C Y : IN. Array of Y-coordinates for the gridpoints C U : OUT. Array of PDE component values for the gridpoints. C NPTS : IN. Number of gridpoints C NPDE : IN. # PDE components C C----------------------------------------------------------------------- C C Burgers' equation, solution wave front at y = x+0.25t, speed of C propagation sqrt(2)/8 perpendicular to wave front. C U = 3/4 - 1/4/(1+exp((-4x+4y-t)/(32*eps))) C V = 3/4 + 1/4/(1+exp((-4x+4y-t)/(32*eps))) C DOUBLE PRECISION EPS PARAMETER (EPS = 1D-3) INTEGER I DO 10 I = 1, NPTS U(I,1) = 0.75 - 0.25/(1+EXP((-4*X(I)+4*Y(I)-T)/(32*EPS))) U(I,2) = 0.75 + 0.25/(1+EXP((-4*X(I)+4*Y(I)-T)/(32*EPS))) 10 CONTINUE RETURN END SUBROUTINE PDEF (T, X, Y, U, UT, UX, UY, UXX, UXY, UYY, RES, + NPTS, NPDE) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER NPTS, NPDE DOUBLE PRECISION T, X(NPTS), Y(NPTS), U(NPTS,NPDE), + UT(NPTS,NPDE), UX(NPTS,NPDE), UY(NPTS,NPDE), + UXX(NPTS,NPDE), UXY(NPTS,NPDE), UYY(NPTS,NPDE), + RES(NPTS,NPDE) C Ccc PURPOSE: C Define residual of PDE on interior of domain. Boundary values will be C overwritten later on. C Ccc PARAMETER DESCRIPTION: C T : IN. Time at which residual should be evaluated C X : IN. Array of X-coordinates for the gridpoints C Y : IN. Array of Y-coordinates for the gridpoints C U : IN. Array of PDE components for the gridpoints. C UT : IN. Array of time derivative of PDE components C UX : IN. -I C UY : IN. I C UXX : IN. I Arrays containing space derivatives of PDE components C UXY : IN. I C UYY : IN. -I C RES : OUT. Array containg PDE residual at gridpoints in interior of C domain. The residual values at boundary points will be C overwritten by a call to PDEBC. C NPTS : IN. Number of gridpoints C NPDE : IN. Number of PDE components C C----------------------------------------------------------------------- C C Burgers' equation Ut = - U.Ux - V.Uy + eps.(Uxx + Uyy) C Vt = - U.Vx - V.Vy + eps.(Vxx + Vyy) C DOUBLE PRECISION EPS PARAMETER (EPS = 1D-3) INTEGER I DO 10 I = 2, NPTS-1 RES(I,1) = UT(I,1) - + (-U(I,1)*UX(I,1) - U(I,2)*UY(I,1) + EPS*(UXX(I,1)+UYY(I,1))) RES(I,2) = UT(I,2) - + (-U(I,1)*UX(I,2) - U(I,2)*UY(I,2) + EPS*(UXX(I,2)+UYY(I,2))) 10 CONTINUE RETURN END SUBROUTINE PDEBC (T, X, Y, U, UT, UX, UY, RES, + NPTS, NPDE, LLBND, ILBND, LBND) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER NPTS, NPDE, LLBND(0:*), ILBND(*), LBND(*) DOUBLE PRECISION T, X(NPTS), Y(NPTS), U(NPTS,NPDE), + UT(NPTS,NPDE), UX(NPTS,NPDE), UY(NPTS,NPDE), + RES(NPTS,NPDE) C Ccc PURPOSE: C Define residual of boundary equations of PDE. The residual on interior C points has already been stored in RES. C Ccc PARAMETER DESCRIPTION: C T : IN. Time at which BC's should be evaluated C X : IN. Array of X-coordinates for the gridpoints C Y : IN. Array of Y-coordinates for the gridpoints C U : IN. Array of PDE components for the gridpoints. C UT : IN. Array of time derivative of PDE components C UX : IN. -I C UY : IN. -I Arrays containing space derivatives of PDE components C RES : INOUT. C IN: PDE residual for interior points (set by PDEF) C OUT: Array with PDE residual at physical boundary points C inserted C NPTS : IN. Number of grid components C NPDE : IN. Number of PDE components C LLBND : IN. (0:LLBND(0)+1) C LLBND(0) = NBNDS: total # physical boundaries in actual grid C LLBND(1:NBNDS): pointers to a specific boundary in LBND C LLBND(NBNDS+1) = NBDPTS+1: total # physical boundary points C in the list + 1 C NB. corners with 2 different types of boundaries should be C pointed at twice. C ILBND : IN. (NBNDS) C ILBND(IB): type of boundary: C 1: Dirichlet C 2: Lower boundary -I C 3: Left boundary I C 4: Upper boundary I max. first order derivative C 5: Right boundary -I C LBND : IN. (NBDPTS) C LBND(LB): pointer to boundary point in actual grid C structure (as in X, Y, and U) C C----------------------------------------------------------------------- C C Burgers' equation, Dirichlet boundaries. C U = 3/4 - 1/4/(1+exp((-4x+4y-t)/(32*eps))) C V = 3/4 + 1/4/(1+exp((-4x+4y-t)/(32*eps))) C DOUBLE PRECISION EPS PARAMETER (EPS = 1D-3) INTEGER I, K, NBNDS NBNDS = LLBND(0) DO 10 K = LLBND(1), LLBND(NBNDS+1)-1 I = LBND(K) RES(I,1) = U(I,1) - + (0.75 - 0.25/(1+EXP((-4*X(I)+4*Y(I)-T)/(32*EPS)))) RES(I,2) = U(I,2) - + (0.75 + 0.25/(1+EXP((-4*X(I)+4*Y(I)-T)/(32*EPS)))) 10 CONTINUE RETURN END SHAR_EOF fi # end of overwriting check if test -f 'exmplr.f' then echo shar: will not over-write existing file "'exmplr.f'" else cat << \SHAR_EOF > 'exmplr.f' PROGRAM EXMPLR C C Restart of EXMPL, default values, Jacobian derivatives exact INTEGER MXLEV, NPD, NPTS, LENIWK, LENRWK, LENLWK PARAMETER (MXLEV=2, NPD=2, NPTS=2500) PARAMETER (LENIWK=NPTS*(5*MXLEV+14), + LENRWK=NPTS*NPD*(5*MXLEV + 9+18*NPD), + LENLWK=2*NPTS) C C----------------------------------------------------------------------- C INTEGER LUNDMP PARAMETER (LUNDMP = 89) C INTEGER NPDE, INFO(1), IWK(LENIWK), MNTR LOGICAL LWK(LENLWK) DOUBLE PRECISION T, TOUT, DT, XL, YL, XR, YU, DX, DY, + TOLS, TOLT, RINFO(1), RWK(LENRWK) C Continuation call of VLUGR2 MNTR = 1 TOUT = 3.0 TOLS = 0.1 TOLT = 0.05 C Default choices INFO(1) = 0 C OPEN(UNIT=LUNDMP,FILE='DUMP',FORM='UNFORMATTED') CALL RDDUMP (LUNDMP, RWK, LENRWK, IWK, LENIWK) CLOSE(LUNDMP) C C call main routine CALL VLUGR2 (NPDE, T, TOUT, DT, XL, YL, XR, YU, DX, DY, + TOLS, TOLT, INFO, RINFO, RWK, LENRWK, IWK, LENIWK, LWK, LENLWK, + MNTR) PRINT *, 'VLUGR2 returned with MNTR=', MNTR C OPEN(UNIT=LUNDMP,FILE='DUMP2',FORM='UNFORMATTED') CALL DUMP (LUNDMP, RWK, IWK) CLOSE(LUNDMP) END SUBROUTINE DERIVF (F, T, X, Y, NPTS, NPDE, U, A0, DT, DX, DY, + LLBND, ILBND, LBND, UIB, UT, UX, UY, UXX, UXY, UYY, + ABSTOL, DEL, WORK, + FU, FUX, FUY, FUXX, FUXY, FUYY) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER NPTS, NPDE, LLBND(0:*), ILBND(*), LBND(*) DOUBLE PRECISION F(NPTS,NPDE), T, X(NPTS), Y(NPTS), U(NPTS,NPDE), + A0, DT, DX, DY, UIB(*), + UT(NPTS,NPDE), UX(NPTS,NPDE), UY(NPTS,NPDE), + UXX(NPTS,NPDE), UXY(NPTS,NPDE), UYY(NPTS,NPDE), + ABSTOL(NPDE), DEL(NPTS), WORK(2*NPTS*NPDE), + FU(NPTS,NPDE,NPDE), FUX(NPTS,NPDE,NPDE), FUY(NPTS,NPDE,NPDE), + FUXX(NPTS,NPDE,NPDE),FUXY(NPTS,NPDE,NPDE),FUYY(NPTS,NPDE,NPDE) C Ccc PURPOSE: C Compute derivatives of residual wrt (derivatives of) U C C PARAMETER DESCRIPTION: C F : IN. Residual F(t,U,Ut) C T : IN. Current time C X,Y : IN. Physical coordinates of gridpoints C NPTS : IN. # grid points C NPDE : IN. # PDE components C U : IN. Solution at T on current grid C A0 : IN. Coefficient of U_n+1 in time derivative C DT : IN. Current time step size C DX : IN. Cell width in X-direction for current grid C DY : IN. Cell width in Y-direction for current grid C LLBND : (0:LLBND(0)+2) C LLBND(0) = NBNDS: total # physical boundaries and corners in C actual domain. C NB. corners should be stored as an independent boundary C (cf. ILBND). The order in LLBND should be first the C boundaries and then the corners. C LLBND(1:NBNDS): pointers to a specific boundary or corner in C LBND C LLBND(NBNDS+1) = NBDPTS+1: total # physical boundary points C in LBND + 1 C LLBND(NBNDS+1): pointer to internal boundary in LBND C LLBND(NBNDS+2) = NBIPTS+1: total # points in LBND + 1 C ILBND : (NBNDS) C ILBND(IB): type of boundary: C 1: Lower boundary -I C 2: Left boundary I C 3: Upper boundary I max. first order derivative C 4: Right boundary -I C 12: Lowerleft corner -I C 23: Leftupper corner I corners of 90 degrees C 34: Upperright corner I (external corners) C 41: Rightlower corner -I max. first order deriv. C 21: Leftlower corner -I C 32: Upperleft corner I corners of 270 degrees C 43: Rightupper corner I (internal corners) C 14: Lowerright corner -I max. first order deriv. C LBND : IN. (NBIPTS) C LBND(IBPT): pointer to boundary point in actual grid C UIB : IN. Solution at T on internal boundaries C UT : IN. Time derivative of U on current grid C UX : IN. -I C UY : IN. I C UXX : IN. I Space derivatives of U on current grid C UXY : IN. I C UYY : IN. -I C ABSTOL : IN. Absolute tolerance for Newton process C DEL : WORK. (NPTS) C WORK : WORK. (2.LENU) C FU : OUT. dF(U,Ut)dU C FUX : OUT. dF(Ux)dUx C FUY : OUT. dF(Uy)dUy C FUXX : OUT. dF(Uxx)dUxx C FUXY : OUT. dF(Uxy)dUxy C FUYY : OUT. dF(Uyy)dUyy C C----------------------------------------------------------------------- C DOUBLE PRECISION EPS PARAMETER (EPS = 1D-3) C INTEGER IC, IPT, LB, NBNDS C Ccc Loop over the components of the (derivatives of) U IC = 1 C C dF(U,Ut)/dU_ic DO 20 IPT = 1, NPTS FU(IPT,1,IC) = A0 - (-UX(IPT,1)) FU(IPT,2,IC) = - (-UX(IPT,2)) 20 CONTINUE C C dF(Ux)/dUx_ic DO 40 IPT = 1, NPTS FUX(IPT,1,IC) = - (-U(IPT,1)) FUX(IPT,2,IC) = 0.0 40 CONTINUE C C dF(Uy)/dUy_ic DO 50 IPT = 1, NPTS FUY(IPT,1,IC) = - (-U(IPT,2)) FUY(IPT,2,IC) = 0.0 50 CONTINUE C C dF(Uxx)/dUxx_ic DO 60 IPT = 1, NPTS FUXX(IPT,1,IC) = - (EPS) FUXX(IPT,2,IC) = 0.0 60 CONTINUE C C dF(Uxy)/dUxy_ic DO 70 IPT = 1, NPTS FUXY(IPT,1,IC) = 0.0 FUXY(IPT,2,IC) = 0.0 70 CONTINUE C C dF(Uyy)/dUyy_ic DO 80 IPT = 1, NPTS FUYY(IPT,1,IC) = - (EPS) FUYY(IPT,2,IC) = 0.0 80 CONTINUE IC = 2 C C dF(U,Ut)/dU_ic DO 120 IPT = 1, NPTS FU(IPT,1,IC) = - (-UY(IPT,1)) FU(IPT,2,IC) = A0 - (-UY(IPT,2)) 120 CONTINUE C C dF(Ux)/dUx_ic DO 140 IPT = 1, NPTS FUX(IPT,1,IC) = 0.0 FUX(IPT,2,IC) = - (-U(IPT,1)) 140 CONTINUE C C dF(Uy)/dUy_ic DO 150 IPT = 1, NPTS FUY(IPT,1,IC) = 0.0 FUY(IPT,2,IC) = - (-U(IPT,2)) 150 CONTINUE C C dF(Uxx)/dUxx_ic DO 160 IPT = 1, NPTS FUXX(IPT,1,IC) = 0.0 FUXX(IPT,2,IC) = - (EPS) 160 CONTINUE C C dF(Uxy)/dUxy_ic DO 170 IPT = 1, NPTS FUXY(IPT,1,IC) = 0.0 FUXY(IPT,2,IC) = 0.0 170 CONTINUE C C dF(Uyy)/dUyy_ic DO 180 IPT = 1, NPTS FUYY(IPT,1,IC) = 0.0 FUYY(IPT,2,IC) = - (EPS) 180 CONTINUE C C Correct boundaries (incl. the internal) NBNDS = LLBND(0) DO 100 LB = LLBND(1), LLBND(NBNDS+2)-1 IPT = LBND(LB) FU(IPT,1,1) = 1.0 FU(IPT,1,2) = 0.0 FU(IPT,2,1) = 0.0 FU(IPT,2,2) = 1.0 FUX(IPT,1,1) = 0.0 FUX(IPT,1,2) = 0.0 FUX(IPT,2,1) = 0.0 FUX(IPT,2,2) = 0.0 FUY(IPT,1,1) = 0.0 FUY(IPT,1,2) = 0.0 FUY(IPT,2,1) = 0.0 FUY(IPT,2,2) = 0.0 FUXX(IPT,1,1) = 0.0 FUXX(IPT,1,2) = 0.0 FUXX(IPT,2,1) = 0.0 FUXX(IPT,2,2) = 0.0 FUXY(IPT,1,1) = 0.0 FUXY(IPT,1,2) = 0.0 FUXY(IPT,2,1) = 0.0 FUXY(IPT,2,2) = 0.0 FUYY(IPT,1,1) = 0.0 FUYY(IPT,1,2) = 0.0 FUYY(IPT,2,1) = 0.0 FUYY(IPT,2,2) = 0.0 100 CONTINUE RETURN END SUBROUTINE MONITR (T, DT, DTNEW, XL, YL, DXB, DYB, + LGRID, ISTRUC, LSOL, SOL) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER LGRID(0:*), ISTRUC(*), LSOL(*) DOUBLE PRECISION T, DT, DTNEW, XL, YL, DXB, DYB, SOL(*) C Ccc PURPOSE: C Control after a successful time step. The solution can be printed, C plotted or compared with the exact solution. C Ccc PARAMETER DESCRIPTION: C T : IN. Current value of time variable C DT : IN. Current time step size C DTNEW : IN. Time step size for next time step C XL : IN. X-coordinate of lowerleft corner of (virtual) rectangle C YL : IN. Y-coordinate of lowerleft corner of (virtual) rectangle C DXB : IN. Cell width in X-direction of base grid C DYB : IN. Cell width in Y-direction of base grid C LGRID : IN. (0:*) C LGRID(0) = max. grid level used at T C LGRID(1): pointer to base grid structure ISTRUC C LGRID(LEVEL): pointer to grid structure (LROW, IROW, ICOL) C of refinement level LEVEL for time T C ISTRUC : IN. (*) C ISTRUC(LGRID(LEVEL):.) contains (LROW,IROW,ICOL) of grid C level LEVEL, C LROW : (0:LROW(0)+1) C LROW(0) = NROWS: Actual # rows in grid C LROW(1:NROWS): pointers to the start of a row in the grid C LROW(NROWS+1) = NPTS+1: Actual # nodes in grid + 1 C IROW : (NROWS) C IROW(IR): row number of row IR in virtual rectangle C ICOL : (NPTS) C ICOL(IPT): column number of grid point IPT in virtual C rectangle C LSOL : IN. (*) C LSOL(LEVEL): pointer to (injected) solution at grid C of refinement level LEVEL for time T C SOL : IN. (*) C SOL(LSOL(LEVEL)+1:LSOL(LEVEL)+NPTS(LEVEL)*NPDE) contains C U_LEVEL(NPTS,NPDE) C C C Local arrays: INTEGER MAXPTS, NPDE PARAMETER (MAXPTS=10000, NPDE=2) DOUBLE PRECISION X(MAXPTS), Y(MAXPTS), UEX(MAXPTS*NPDE) C C----------------------------------------------------------------------- C INTEGER MAXLEV, LEVEL, LLROW, NROWS, NPTS, LIROW, LICOL DOUBLE PRECISION DX, DY C C Loop over the grid levels from coarse to fine. C Get physical coordinates of grid points C Compute ||err||_max MAXLEV = LGRID(0) DX = DXB DY = DYB DO 10 LEVEL = 1, MAXLEV LLROW = LGRID(LEVEL) NROWS = ISTRUC(LLROW) NPTS = ISTRUC(LLROW+NROWS+1)-1 LIROW = LLROW+NROWS+2 LICOL = LIROW+NROWS CALL SETXY (XL, YL, DX, DY, + ISTRUC(LLROW), ISTRUC(LIROW), ISTRUC(LICOL), X, Y) DX = DX/2 DY = DY/2 CALL PRERR (LEVEL, T, NPTS, NPDE, X, Y, SOL(LSOL(LEVEL)+1), + UEX) 10 CONTINUE RETURN END SUBROUTINE PRERR (LEVEL, T, NPTS, NPDE, X, Y, U, UEX) INTEGER LEVEL, NPTS, NPDE DOUBLE PRECISION T, X(NPTS), Y(NPTS), U(NPTS,NPDE), UEX(NPTS,NPDE) INTEGER I,J DOUBLE PRECISION RMAX(2) CALL PDEIV (T, X, Y, UEX, NPTS, NPDE) RMAX(1) = 0.0 RMAX(2) = 0.0 DO 10 I = 1, NPTS J = 1 RMAX(J) = MAX(RMAX(J),ABS(UEX(I,J)-U(I,J))) J = 2 RMAX(J) = MAX(RMAX(J),ABS(UEX(I,J)-U(I,J))) 10 CONTINUE PRINT '(''Error at T='',E9.3,'', level='',I1,'' :'',2E11.3,I10)', + T, LEVEL, RMAX(1), RMAX(2), NPTS RETURN END SUBROUTINE PDEIV (T, X, Y, U, NPTS, NPDE) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER NPTS, NPDE DOUBLE PRECISION T, X(NPTS), Y(NPTS), U(NPTS,NPDE) C Ccc PURPOSE: C Define (initial) solution of PDE. C Ccc PARAMETER DESCRIPTION: C T : IN. Time at which (initial) solution should be given C X : IN. Array of X-coordinates for the gridpoints C Y : IN. Array of Y-coordinates for the gridpoints C U : OUT. Array of PDE component values for the gridpoints. C NPTS : IN. Number of gridpoints C NPDE : IN. # PDE components C C----------------------------------------------------------------------- C C Burgers' equation, solution wave front at y = x+0.25t, speed of C propagation sqrt(2)/8 perpendicular to wave front. C U = 3/4 - 1/4/(1+exp((-4x+4y-t)/(32*eps))) C V = 3/4 + 1/4/(1+exp((-4x+4y-t)/(32*eps))) C DOUBLE PRECISION EPS PARAMETER (EPS = 1D-3) INTEGER I DO 10 I = 1, NPTS U(I,1) = 0.75 - 0.25/(1+EXP((-4*X(I)+4*Y(I)-T)/(32*EPS))) U(I,2) = 0.75 + 0.25/(1+EXP((-4*X(I)+4*Y(I)-T)/(32*EPS))) 10 CONTINUE RETURN END SUBROUTINE PDEF (T, X, Y, U, UT, UX, UY, UXX, UXY, UYY, RES, + NPTS, NPDE) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER NPTS, NPDE DOUBLE PRECISION T, X(NPTS), Y(NPTS), U(NPTS,NPDE), + UT(NPTS,NPDE), UX(NPTS,NPDE), UY(NPTS,NPDE), + UXX(NPTS,NPDE), UXY(NPTS,NPDE), UYY(NPTS,NPDE), + RES(NPTS,NPDE) C Ccc PURPOSE: C Define residual of PDE on interior of domain. Boundary values will be C overwritten later on. C Ccc PARAMETER DESCRIPTION: C T : IN. Time at which residual should be evaluated C X : IN. Array of X-coordinates for the gridpoints C Y : IN. Array of Y-coordinates for the gridpoints C U : IN. Array of PDE components for the gridpoints. C UT : IN. Array of time derivative of PDE components C UX : IN. -I C UY : IN. I C UXX : IN. I Arrays containing space derivatives of PDE components C UXY : IN. I C UYY : IN. -I C RES : OUT. Array containg PDE residual at gridpoints in interior of C domain. The residual values at boundary points will be C overwritten by a call to PDEBC. C NPTS : IN. Number of gridpoints C NPDE : IN. Number of PDE components C C----------------------------------------------------------------------- C C Burgers' equation Ut = - U.Ux - V.Uy + eps.(Uxx + Uyy) C Vt = - U.Vx - V.Vy + eps.(Vxx + Vyy) C DOUBLE PRECISION EPS PARAMETER (EPS = 1D-3) INTEGER I DO 10 I = 2, NPTS-1 RES(I,1) = UT(I,1) - + (-U(I,1)*UX(I,1) - U(I,2)*UY(I,1) + EPS*(UXX(I,1)+UYY(I,1))) RES(I,2) = UT(I,2) - + (-U(I,1)*UX(I,2) - U(I,2)*UY(I,2) + EPS*(UXX(I,2)+UYY(I,2))) 10 CONTINUE RETURN END SUBROUTINE PDEBC (T, X, Y, U, UT, UX, UY, RES, + NPTS, NPDE, LLBND, ILBND, LBND) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER NPTS, NPDE, LLBND(0:*), ILBND(*), LBND(*) DOUBLE PRECISION T, X(NPTS), Y(NPTS), U(NPTS,NPDE), + UT(NPTS,NPDE), UX(NPTS,NPDE), UY(NPTS,NPDE), + RES(NPTS,NPDE) C Ccc PURPOSE: C Define residual of boundary equations of PDE. The residual on interior C points has already been stored in RES. C Ccc PARAMETER DESCRIPTION: C T : IN. Time at which BC's should be evaluated C X : IN. Array of X-coordinates for the gridpoints C Y : IN. Array of Y-coordinates for the gridpoints C U : IN. Array of PDE components for the gridpoints. C UT : IN. Array of time derivative of PDE components C UX : IN. -I C UY : IN. -I Arrays containing space derivatives of PDE components C RES : INOUT. C IN: PDE residual for interior points (set by PDEF) C OUT: Array with PDE residual at physical boundary points C inserted C NPTS : IN. Number of grid components C NPDE : IN. Number of PDE components C LLBND : IN. (0:LLBND(0)+1) C LLBND(0) = NBNDS: total # physical boundaries in actual grid C LLBND(1:NBNDS): pointers to a specific boundary in LBND C LLBND(NBNDS+1) = NBDPTS+1: total # physical boundary points C in the list + 1 C NB. corners with 2 different types of boundaries should be C pointed at twice. C ILBND : IN. (NBNDS) C ILBND(IB): type of boundary: C 1: Dirichlet C 2: Lower boundary -I C 3: Left boundary I C 4: Upper boundary I max. first order derivative C 5: Right boundary -I C LBND : IN. (NBDPTS) C LBND(LB): pointer to boundary point in actual grid C structure (as in X, Y, and U) C C----------------------------------------------------------------------- C C Burgers' equation, Dirichlet boundaries. C U = 3/4 - 1/4/(1+exp((-4x+4y-t)/(32*eps))) C V = 3/4 + 1/4/(1+exp((-4x+4y-t)/(32*eps))) C DOUBLE PRECISION EPS PARAMETER (EPS = 1D-3) INTEGER I, K, NBNDS NBNDS = LLBND(0) DO 10 K = LLBND(1), LLBND(NBNDS+1)-1 I = LBND(K) RES(I,1) = U(I,1) - + (0.75 - 0.25/(1+EXP((-4*X(I)+4*Y(I)-T)/(32*EPS)))) RES(I,2) = U(I,2) - + (0.75 + 0.25/(1+EXP((-4*X(I)+4*Y(I)-T)/(32*EPS)))) 10 CONTINUE RETURN END SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Sp' then mkdir 'Sp' fi cd 'Sp' if test -f 'prtsol.f' then echo shar: will not over-write existing file "'prtsol.f'" else cat << \SHAR_EOF > 'prtsol.f' PROGRAM PRTSOL C C----------------------------------------------------------------------- C Ccc This program reads a file made by subroutine DUMP and prints the C solution on an output file. Both filenames are read from standard C input. C Ccc EXTERNALS USED: EXTERNAL PRSOL, RDDUMP C C Ccc INCLUDE 'CMNWRITEF' C C CMNWRITEF C C COMMON needed for continuation calls INTEGER MAXLVW, NPDEW, LRWKPS, LIWKPS, LRWKB, LIWKB LOGICAL FIRST, SECOND REAL T0, TW, TEW, DTW, XLW, YLW, XRW, YUW, DXB, DYB, DTO COMMON /WRITIF/ MAXLVW, NPDEW, LRWKPS, LIWKPS, LRWKB, LIWKB COMMON /WRITLF/ FIRST, SECOND COMMON /WRITRF/ T0, TW, TEW, DTW, XLW,YLW, XRW,YUW, DXB, DYB, DTO SAVE /WRITIF/, /WRITLF/, /WRITRF/ C C end INCLUDE 'CMNWRITEF' C C C----------------------------------------------------------------------- C INTEGER MXLEV, NPD, NPTS, LENIWK, LENRWK PARAMETER (MXLEV=5, NPD=3, NPTS=10000) PARAMETER (LENIWK=NPTS*(7*MXLEV+20), + LENRWK=5*NPTS*NPD*MXLEV) C CHARACTER FILE*128 INTEGER IWK(LENIWK), + LSGNM1, LSGN, LSGNP1, LSUNM1, LSSN, LSUN REAL RWK(LENRWK) PRINT *, 'DUMP file?' READ '(A)', FILE C OPEN(UNIT=62,FILE=FILE,FORM='UNFORMATTED') CALL RDDUMP (62, RWK, LENRWK, IWK, LENIWK) CLOSE(62) C C Setup work storage LSGNM1 = 1 LSGN = LSGNM1 + MAXLVW+1 LSGNP1 = LSGN + MAXLVW+1 LSUNM1 = LSGNP1 + MAXLVW+1 LSSN = LSUNM1 + MAXLVW LSUN = LSSN + MAXLVW C C call print routine PRINT *, 'output file?' READ '(A)', FILE C OPEN(UNIT=61,FILE=FILE) CALL PRSOL (61, TW, NPDEW, XLW, YLW, DXB, DYB, + IWK(LSGN), IWK(LIWKPS), IWK(LSUN), RWK(LRWKPS)) CLOSE(61) END SUBROUTINE RDDUMP (LUNDMP, RWK, LENRWK, IWK, LENIWK) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER LENIWK INTEGER LUNDMP, LENRWK, IWK(LENIWK) REAL RWK(LENRWK) C Ccc PURPOSE: C Read all information necessary for a restart of VLUGR2 from file C Ccc PARAMETER DESCRIPTION: C LUNDMP : IN. Logical unit number of dumpfile. Should be opened as an C unformatted file. C RWK : OUT. Real workstorage intended to pass to VLUGR2 C LENRWK : IN. Dimension of RWK. C IWK : OUT. Integer workstorage intended to pass to VLUGR2 C LENIWK : IN. Dimension of IWK. C Ccc EXTERNALS USED: NONE C C Ccc INCLUDE 'CMNSTATS' C C CMNSTATS C C COMMON with integration statistics INTEGER MXCLEV, MXCNIT PARAMETER (MXCLEV = 10, MXCNIT = 20) INTEGER LUNPDS, LUNNLS, LUNLSS, LEVEL, NSTEPS, NREJS, + NJACS(MXCLEV), NRESID(MXCLEV), NNIT(MXCLEV), + NLSIT(MXCLEV,MXCNIT) COMMON /STATS/ LUNPDS, LUNNLS, LUNLSS, LEVEL, NSTEPS, NREJS, + NJACS, NRESID, NNIT, NLSIT SAVE /STATS/ C C end INCLUDE 'CMNSTATS' C C Ccc INCLUDE 'CMNWRITEF' C C CMNWRITEF C C COMMON needed for continuation calls INTEGER MAXLVW, NPDEW, LRWKPS, LIWKPS, LRWKB, LIWKB LOGICAL FIRST, SECOND REAL T0, TW, TEW, DTW, XLW, YLW, XRW, YUW, DXB, DYB, DTO COMMON /WRITIF/ MAXLVW, NPDEW, LRWKPS, LIWKPS, LRWKB, LIWKB COMMON /WRITLF/ FIRST, SECOND COMMON /WRITRF/ T0, TW, TEW, DTW, XLW,YLW, XRW,YUW, DXB, DYB, DTO SAVE /WRITIF/, /WRITLF/, /WRITRF/ C C end INCLUDE 'CMNWRITEF' C C C----------------------------------------------------------------------- C INTEGER I, J READ(LUNDMP) MAXLVW, NPDEW, LRWKPS, LIWKPS, LRWKB, LIWKB, + FIRST, SECOND, + T0, TW, TEW, DTW, XLW, YLW, XRW, YUW, DXB, DYB, DTO IF (LENRWK .LT. LRWKPS+LRWKB .OR. LENIWK .LT. LIWKPS+LIWKB) THEN PRINT *, LENRWK, LRWKPS+LRWKB, LENIWK, LIWKPS+LIWKB STOP 'work space too small' ENDIF READ(LUNDMP) LUNPDS, LUNNLS, LUNLSS, LEVEL, NSTEPS, NREJS, + (NJACS(I), I=1,MXCLEV), (NRESID(I), I=1,MXCLEV), + (NNIT(I), I=1,MXCLEV), ((NLSIT(I,J), I=1,MXCLEV), J=1,MXCNIT) READ(LUNDMP) (RWK(I), I=1,LRWKPS+LRWKB) READ(LUNDMP) (IWK(I), I=1,LIWKPS+LIWKB) C RETURN END SUBROUTINE PRSOL (LUN, T, NPDE, XL, YL, DXB, DYB, LGRID, ISTRUC, + LSOL, SOL) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER LUN, NPDE, LGRID(0:*), ISTRUC(*), LSOL(*) REAL T, XL, YL, DXB, DYB, SOL(*) C Ccc PURPOSE: C Print solution and coordinate values at all grid levels. C Ccc PARAMETER DESCRIPTION: C LUN : IN. Logical unit number of print file C T : IN. Current value of time variable C NPDE : IN. # PDE components C XL : IN. X-coordinate of lowerleft corner of (virtual) domain C YL : IN. Y-coordinate of lowerleft corner of (virtual) domain C DXB : IN. Cell width in X-direction of base grid C DYB : IN. Cell width in Y-direction of base grid C LGRID : IN. (0:*) C LGRID(0) = max. grid level used at T C LGRID(1): pointer to base grid structure ISTRUC C LGRID(LEVEL): pointer to grid structure (LROW, IROW, ICOL) C of refinement level LEVEL for time T C ISTRUC : IN. (*) C ISTRUC(LGRID(LEVEL):.) contains (LROW,IROW,ICOL) of grid C level LEVEL, C LROW : (0:LROW(0)+1) C LROW(0) = NROWS: Actual # rows in grid C LROW(1:NROWS): pointers to the start of a row in the grid C LROW(NROWS+1) = NPTS+1: Actual # nodes in grid + 1 C IROW : (NROWS) C IROW(IR): row number of row IR in virtual rectangle C ICOL : (NPTS) C ICOL(IPT): column number of grid point IPT in virtual C rectangle C LSOL : IN. (*) C LSOL(LEVEL): pointer to (injected) solution at grid C of refinement level LEVEL for time T C SOL : IN. (*) C SOL(LSOL(LEVEL)+1:LSOL(LEVEL)+NPTS(LEVEL)*NPDE) contains C U_LEVEL(NPTS,NPDE) C Ccc EXTERNALS USED: EXTERNAL PRSOLL C C----------------------------------------------------------------------- C INTEGER MAXLEV, LEVEL, LLROW, NROWS, NPTS, LIROW, LICOL REAL DX, DY MAXLEV = LGRID(0) DX = DXB DY = DYB DO 10 LEVEL = 1, MAXLEV LLROW = LGRID(LEVEL) NROWS = ISTRUC(LLROW) NPTS = ISTRUC(LLROW+NROWS+1)-1 LIROW = LLROW+NROWS+2 LICOL = LIROW+NROWS CALL PRSOLL (LUN, LEVEL, T, NPTS, NPDE, XL, YL, DX, DY, + ISTRUC(LLROW), ISTRUC(LIROW), ISTRUC(LICOL), + SOL(LSOL(LEVEL)+1)) DX = DX/2 DY = DY/2 10 CONTINUE RETURN END SUBROUTINE PRSOLL (LUN, LEVEL, T, NPTS, NPDE, XL, YL, DX, DY, + LROW, IROW, ICOL, U) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER LUN, LEVEL, NPTS, NPDE, LROW(0:*), IROW(*), ICOL(*) REAL T, XL, YL, DX, DY, U(NPTS,NPDE) C Ccc PURPOSE: C Print solution and X- and Y-coordinates of gridlevel LEVEL. C Ccc PARAMETER DESCRIPTION: C LUN : IN. Logical unit number of print file C LEVEL : IN. Grid level corresponding with solution U. C T : IN. Current value of time variable C NPTS : IN. # grid points at this level C NPDE : IN. # PDE components C XL : IN. X-coordinate of lower-left point of virtual rectangle C YL : IN. Y-coordinate of lower-left point of virtual rectangle C DX : IN. Grid width in X-direction C DY : IN. Grid width in Y-direction C LROW : IN. (0:LROW(0)+1) C LROW(0) = NROWS: Actual # rows in grid C LROW(1:NROWS): pointers to the start of a row in the grid C LROW(NROWS+1) = NPTS+1: Actual # nodes in grid + 1 C IROW : IN. (NROWS) C IROW(IR): row number of row IR in virtual rectangle C ICOL : IN. (NPTS) C ICOL(IPT): column number of grid point IPT in virtual C rectangle C U : IN. Solution on this grid level C Ccc EXTERNALS USED: NONE C C----------------------------------------------------------------------- C INTEGER IC, IPT, IR, NROWS REAL X, Y C NROWS = LROW(0) WRITE(LUN,'(//// T10,A,T30,A,T46,A,T62,A,T71,A //)') + 'Level', 't', 'Y', 'X', 'Solution' IR = 1 Y = YL + IROW(IR)*DY IPT = LROW(IR) X = XL + ICOL(IPT)*DX WRITE(LUN, + '(T13,I2,T21,E12.5,T37,E12.5,T53,E12.5,T69,E12.5)') + LEVEL, T, Y, X, U(IPT,1) DO 10 IC = 2, NPDE WRITE(LUN,'(T69,E12.5)') U(IPT,IC) 10 CONTINUE DO 20 IPT = LROW(IR)+1, LROW(IR+1)-1 X = XL + ICOL(IPT)*DX WRITE(LUN,'(T53,E12.5,T69,E12.5)') X, U(IPT,1) DO 30 IC = 2, NPDE WRITE(LUN,'(T69,E12.5)') U(IPT,IC) 30 CONTINUE 20 CONTINUE DO 40 IR = 2, NROWS Y = YL + IROW(IR)*DY IPT = LROW(IR) X = XL + ICOL(IPT)*DX WRITE(LUN, + '(T21,E12.5,T37,E12.5,T53,E12.5,T69,E12.5)') + T, Y, X, U(IPT,1) DO 50 IC = 2, NPDE WRITE(LUN,'(T69,E12.5)') U(IPT,IC) 50 CONTINUE DO 60 IPT = LROW(IR)+1, LROW(IR+1)-1 X = XL + ICOL(IPT)*DX WRITE(LUN,'(T53,E12.5,T69,E12.5)') X, U(IPT,1) DO 70 IC = 2, NPDE WRITE(LUN,'(T69,E12.5)') U(IPT,IC) 70 CONTINUE 60 CONTINUE 40 CONTINUE RETURN END SHAR_EOF fi # end of overwriting check if test -f 'probi.f' then echo shar: will not over-write existing file "'probi.f'" else cat << \SHAR_EOF > 'probi.f' PROGRAM EXMPL C C Ccc INCLUDE 'PARNEWTON' C C PARNEWTON C C Parameters for Newton process C MAXNIT : Max. number of Newton iterations C MAXJAC : Max. number of Jacobian / preconditioner evaluations during C a Newton process C TOLNEW : Tolerance for Newton process: C rho/(1-rho)*|| corr.||_w < TOLNEW INTEGER MAXNIT, MAXJAC REAL TOLNEW PARAMETER (MAXNIT = 10, MAXJAC = 2, TOLNEW = 1.0) C C end INCLUDE 'PARNEWTON' C C Ccc INCLUDE 'PARGCRO' C C PARGCRO C C Parameters for linear system solver GCRO + (block-)diagonal C preconditioner C IDIAGP : 0: block-diagonal + first order derivatives C 1: block-diagonal neglecting first order derivatives C 2: diagonal + first order derivatives C 3: diagonal neglecting first order derivatives C NRRMAX : Max. number of restarts of outer loop C MAXLR : Max. number of iterations in outer loop C MAXL : Max. number of iterations in GMRES inner loop C TOLLSC : Tolerance for linear system solver INTEGER IDIAGP, NRRMAX, MAXLR, MAXL REAL TOLLSC PARAMETER (NRRMAX = 1, MAXLR = 5, MAXL = 20) C PARAMETER (NRRMAX = 1, MAXLR = 3, MAXL = 10) PARAMETER (TOLLSC = TOLNEW/10) COMMON /IGCRO/ IDIAGP SAVE /IGCRO/ C C end INCLUDE 'PARGCRO' C C INTEGER MXLEV, NPD, NPTS, LENIWK, LENRWK, LENLWK PARAMETER (MXLEV=2, NPD=3, NPTS=5000) PARAMETER (LENIWK=NPTS*(5*MXLEV+14), + LENRWK=NPTS*NPD*(5*MXLEV+9 + + 9*NPD+(2*MAXLR+MAXL+6+NPD)), + LENLWK=2*NPTS) C C----------------------------------------------------------------------- C INTEGER LUNDMP PARAMETER (LUNDMP = 89) C CHARACTER FILE*7 INTEGER NPDE, INFO(7), IWK(LENIWK), MNTR, I LOGICAL LWK(LENLWK) REAL T, TOUT(4), DT, XL, YL, XR, YU, DX, DY, + TOLS, TOLT, RINFO(2+3*NPD), RWK(LENRWK) C C First call of VLUGR2 MNTR = 0 NPDE = 3 T = 0.0 TOUT(1) = 500.0 TOUT(2) = 5000.0 TOUT(3) = 10000.0 TOUT(4) = 20000.0 DT = 0.1 XL = 0.0 XR = 1.0 YL = 0.0 YU = 1.0 DX = 0.05 DY = 0.05 TOLS = 0.1 TOLT = 0.1 INFO(1) = 1 C MAXLEV INFO(2) = 3 C Domain is a rectangle INFO(3) = 0 C Linear system solver PRINT *, 'Lin.sys.solver; BiCGStab, GCRO or matrix-free GCRO ?' PRINT *, ' (0 / 10,11,12,13 / 20,21,22,23 ) ?' READ *, INFO(4) OPEN (UNIT=61,FILE='RunInfo') C Write integration history to unit # 61 INFO(5) = 61 C Write Newton info to unit # 61 INFO(6) = 61 C Write Linear system solver info to unit # 61 INFO(7) = 61 C DTMIN = 1E-3 RINFO(1) = 1.0E-3 C DTMAX = 1.0 RINFO(2) = 10000.0 C UMAX RINFO(3) = 1.1E+5 RINFO(4) = 0.25 RINFO(5) = 292.0 C SPCWGT = 1.0 RINFO(6) = 1.0 RINFO(7) = 1.0 RINFO(8) = 1.0 C TIMWGT = 1.0 RINFO( 9) = 1.0 RINFO(10) = 1.0 RINFO(11) = 1.0 C C Call main routine FILE='DUMP' DO 10 I = 1, 4 CALL VLUGR2 (NPDE, T, TOUT(I), DT, XL, YL, XR, YU, DX, DY, + TOLS, TOLT, INFO, RINFO, RWK, LENRWK, IWK, LENIWK, + LWK, LENLWK, MNTR) C C Save info on file WRITE(FILE(5:7),'(I3.3)') I OPEN(UNIT=LUNDMP,FILE=FILE,FORM='UNFORMATTED') CALL DUMP (LUNDMP, RWK, IWK) CLOSE(LUNDMP) C Check MNTR value IF (MNTR .NE. 1) THEN PRINT *, 'VLUGR2 returned with MNTR=', MNTR STOP ENDIF 10 CONTINUE END SUBROUTINE PDEIV (T, X, Y, U, NPTS, NPDE) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER NPTS, NPDE REAL T, X(NPTS), Y(NPTS), U(NPTS,NPDE) C Ccc PURPOSE: C Define (initial) solution of PDE. C Ccc PARAMETER DESCRIPTION: C T : IN. Time at which (initial) solution should be given C X : IN. Array of X-coordinates for the gridpoints C Y : IN. Array of Y-coordinates for the gridpoints C U : OUT. Array of PDE component values for the gridpoints. C NPTS : IN. Number of gridpoints C NPDE : IN. # PDE components C C----------------------------------------------------------------------- C INTEGER I C REAL N, GAMMA, MU0, RHO0, P0, W0, G, DM, KAPPA, AL, AT, + CF, TKAPPA, LT, LL, CS, RHOS, T0, ALPHA, BETA, QC, TC COMMON /PROBLM/ N, GAMMA, MU0, RHO0, P0, W0, G, DM, KAPPA, AL, AT, + CF, TKAPPA, LT, LL, CS, RHOS, T0, ALPHA, BETA, QC, TC SAVE /PROBLM/ C Ccc Problem parameters N = 0.4 KAPPA = 1.0E-10 G = 9.81 DM = 0.0 AT = 0.002 AL = 0.01 CF = 4182.0 TKAPPA = 4.0 LT = 0.001 LL = 0.01 CS = 840.0 RHOS = 2500.0 RHO0 = 1.0E+3 T0 = 290.0 P0 = 1.0E+5 ALPHA = -3.0E-4 BETA = 4.45E-10 GAMMA = LOG(1.2) MU0 = 1.0E-3 W0 = 0.25 QC = 1.0E-4 TC = 292.0 C Ccc Initial solution DO 10 I = 1, NPTS U(I,1) = P0 + (1.0 - Y(I))*RHO0*G U(I,2) = 0.0 U(I,3) = T0 10 CONTINUE RETURN END SUBROUTINE PDEF (T, X, Y, U, UT, UX, UY, UXX, UXY, UYY, RES, + NPTS, NPDE) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER NPTS, NPDE REAL T, X(NPTS), Y(NPTS), U(NPTS,NPDE), + UT(NPTS,NPDE), UX(NPTS,NPDE), UY(NPTS,NPDE), + UXX(NPTS,NPDE), UXY(NPTS,NPDE), UYY(NPTS,NPDE), + RES(NPTS,NPDE) C Ccc PURPOSE: C Define residual of PDE on interior of domain. Boundary values will be C overwritten later on. C Ccc PARAMETER DESCRIPTION: C T : IN. Time at which residual should be evaluated C X : IN. Array of X-coordinates for the gridpoints C Y : IN. Array of Y-coordinates for the gridpoints C U : IN. Array of PDE components for the gridpoints. C UT : IN. Array of time derivative of PDE components C UX : IN. -I C UY : IN. I C UXX : IN. I Arrays containing space derivatives of PDE components C UXY : IN. I C UYY : IN. -I C RES : OUT. Array containg PDE residual at gridpoints in interior of C domain. The residual values at boundary points will be C overwritten by a call to PDEBC. C NPTS : IN. Number of gridpoints C NPDE : IN. Number of PDE components C C C Ccc INCLUDE 'CMNCMMACH' C C CMNCMMACH C C COMMON with `machine numbers' C LUNOUT : Logical unit # standard output -I C LUNERR : Logical unit # standard error I Set in the routine C UROUND : Smallest machine number such that I MACNUM C 1.0+UROUND > 1.0 and 1.0-UROUND < 1.0 I C XMIN : Smallest floating-point number -I INTEGER LUNOUT, LUNERR REAL UROUND, XMIN COMMON /IMACH/ LUNOUT, LUNERR COMMON /RMACH/ UROUND, XMIN SAVE /IMACH/, /RMACH/ C C end INCLUDE 'CMNCMMACH' C C----------------------------------------------------------------------- C REAL N, GAMMA, MU0, RHO0, P0, W0, G, DM, KAPPA, AL, AT, + CF, TKAPPA, LT, LL, CS, RHOS, T0, ALPHA, BETA, QC, TC COMMON /PROBLM/ N, GAMMA, MU0, RHO0, P0, W0, G, DM, KAPPA, AL, AT, + CF, TKAPPA, LT, LL, CS, RHOS, T0, ALPHA, BETA, QC, TC SAVE /PROBLM/ C INTEGER I REAL P, PT, PX, PY, W, WT, WX, WY, T1, TT, TX, TY, + RHO, RHOX, RHOY, + MU, MUX, MUY, KAPMU, KAPMU2, KAPMUX, KAPMUY, Q1, Q2, QL, + ND11, ND12, ND22, H11, H12, H22, + PXX, PXY, PYY, WXX, WXY, WYY, TXX, TXY, TYY, + ND11Q1, ND11Q2, ND12Q1, ND12Q2, ND22Q1, ND22Q2, + H11Q1, H11Q2, H12Q1, H12Q2, H22Q1, H22Q2, Q1X, Q1Y, Q2X, Q2Y, + ND11X, ND12X, ND12Y, ND22Y, JW1, JW2, JW1X, JW2Y, + H11X, H12X, H12Y, H22Y, JT1X, JT2Y C DO 10 I = 1, NPTS P = U(I,1) PT = UT(I,1) PX = UX(I,1) PY = UY(I,1) W = U(I,2) WT = UT(I,2) WX = UX(I,2) WY = UY(I,2) T1 = U(I,3) TT = UT(I,3) TX = UX(I,3) TY = UY(I,3) RHO = RHO0*EXP(ALPHA*(T1-T0)+BETA*(P-P0)+GAMMA*W) RHOX = RHO*(ALPHA*TX+BETA*PX+GAMMA*WX) RHOY = RHO*(ALPHA*TY+BETA*PY+GAMMA*WY) MU = MU0*(1+1.85*W-4.0*W*W) MUX = MU0*(1.85*WX-8.0*W*WX) MUY = MU0*(1.85*WY-8.0*W*WY) KAPMU = KAPPA/MU KAPMU2 = -KAPMU/MU KAPMUX = KAPMU2*MUX KAPMUY = KAPMU2*MUY Q1 = -KAPMU*PX Q2 = -KAPMU*(PY+RHO*G) QL = MAX(SQRT(Q1*Q1+Q2*Q2),UROUND) ND11 = N*DM + AT*QL + (AL-AT)*Q1*Q1/QL ND12 = (AL-AT)*Q1*Q2/QL ND22 = N*DM + AT*QL + (AL-AT)*Q2*Q2/QL H11 = TKAPPA + LT*QL + (LL-LT)*Q1*Q1/QL H12 = (LL-LT)*Q1*Q2/QL H22 = TKAPPA + LT*QL + (LL-LT)*Q2*Q2/QL PXX = UXX(I,1) PXY = UXY(I,1) PYY = UYY(I,1) WXX = UXX(I,2) WXY = UXY(I,2) WYY = UYY(I,2) TXX = UXX(I,3) TXY = UXY(I,3) TYY = UYY(I,3) ND11Q1 = (AT + (AL-AT)*(2-(Q1/QL)**2))*Q1/QL ND11Q2 = (AT - (AL-AT)*((Q1/QL)**2))*Q2/QL ND12Q1 = (AL-AT)*(Q2/QL)**3 ND12Q2 = (AL-AT)*(Q1/QL)**3 ND22Q1 = (AT - (AL-AT)*((Q2/QL)**2))*Q1/QL ND22Q2 = (AT + (AL-AT)*(2-(Q2/QL)**2))*Q2/QL H11Q1 = (LT + (LL-LT)*(2-(Q1/QL)**2))*Q1/QL H11Q2 = (LT - (LL-LT)*((Q1/QL)**2))*Q2/QL H12Q1 = (LL-LT)*(Q2/QL)**3 H12Q2 = (LL-LT)*(Q1/QL)**3 H22Q1 = (LT - (LL-LT)*((Q2/QL)**2))*Q1/QL H22Q2 = (LT + (LL-LT)*(2-(Q2/QL)**2))*Q2/QL Q1X = -(KAPMUX*PX+KAPMU*PXX) Q1Y = -(KAPMUY*PX+KAPMU*PXY) Q2X = -(KAPMUX*(PY+RHO*G)+KAPMU*(PXY+RHOX*G)) Q2Y = -(KAPMUY*(PY+RHO*G)+KAPMU*(PYY+RHOY*G)) ND11X = ND11Q1*Q1X + ND11Q2*Q2X ND12X = ND12Q1*Q1X + ND12Q2*Q2X ND12Y = ND12Q1*Q1Y + ND12Q2*Q2Y ND22Y = ND22Q1*Q1Y + ND22Q2*Q2Y JW1 = -(ND11*WX + ND12*WY) JW2 = -(ND12*WX + ND22*WY) JW1X = -(ND11X*WX+ND11*WXX + ND12X*WY+ND12*WXY) JW2Y = -(ND12Y*WX+ND12*WXY + ND22Y*WY+ND22*WYY) H11X = H11Q1*Q1X + H11Q2*Q2X H12X = H12Q1*Q1X + H12Q2*Q2X H12Y = H12Q1*Q1Y + H12Q2*Q2Y H22Y = H22Q1*Q1Y + H22Q2*Q2Y JT1X = -(H11X*TX+H11*TXX + H12X*TY+H12*TXY) JT2Y = -(H12Y*TX+H12*TXY + H22Y*TY+H22*TYY) C RES(I,1) = N*RHO*(BETA*PT+GAMMA*WT+ALPHA*TT) + + RHOX*Q1+RHO*Q1X + RHOY*Q2+RHO*Q2Y RES(I,2) = N*RHO*WT + + RHO*Q1*WX + RHO*Q2*WY + + RHOX*JW1+RHO*JW1X + RHOY*JW2+RHO*JW2Y RES(I,3) = (N*CF*RHO+(1-N)*RHOS*CS)*TT + + RHO*CF*Q1*TX + RHO*CF*Q2*TY + + JT1X + JT2Y 10 CONTINUE RETURN END SUBROUTINE PDEBC (T, X, Y, U, UT, UX, UY, RES, + NPTS, NPDE, LLBND, ILBND, LBND) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER NPTS, NPDE, LLBND(0:*), ILBND(*), LBND(*) REAL T, X(NPTS), Y(NPTS), U(NPTS,NPDE), + UT(NPTS,NPDE), UX(NPTS,NPDE), UY(NPTS,NPDE), + RES(NPTS,NPDE) C Ccc PURPOSE: C Define residual of boundary equations of PDE. The residual on interior C points has already been stored in RES. C Ccc PARAMETER DESCRIPTION: C T : IN. Time at which BC's should be evaluated C X : IN. Array of X-coordinates for the gridpoints C Y : IN. Array of Y-coordinates for the gridpoints C U : IN. Array of PDE components for the gridpoints. C UT : IN. Array of time derivative of PDE components C UX : IN. -I C UY : IN. -I Arrays containing space derivatives of PDE components C RES : INOUT. C IN: PDE residual for interior points (set by PDEF) C OUT: Array with PDE residual at physical boundary points C inserted C NPTS : IN. Number of grid components C NPDE : IN. Number of PDE components C LLBND : (0:LLBND(0)+2) C LLBND(0) = NBNDS: total # physical boundaries and corners in C actual domain. C NB. corners should be stored as an independent boundary C (cf. ILBND). The order in LLBND should be first the C boundaries and then the corners. C LLBND(1:NBNDS): pointers to a specific boundary or corner in C LBND C LLBND(NBNDS+1) = NBDPTS+1: total # physical boundary points C in LBND + 1 C LLBND(NBNDS+1): pointer to internal boundary in LBND C LLBND(NBNDS+2) = NBIPTS+1: total # points in LBND + 1 C ILBND : (NBNDS) C ILBND(IB): type of boundary: C 1: Lower boun