C ALGORITHM 759, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 22, NO. 3, September, 1996, P. 329--347. 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:41:19 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' `VLUGR3: A Vectorizable Adaptive Grid Solver for PDEs in 3D. Part II. Code Description' by J.G. Blom and J.G. Verwer. This code solves systems of PDEs of the type F(t,x,y,z,U,Ut,Ux,Uy,Uz,Uxx,Uyy,Uzz,Uxy,Uxz,Uyz)=0 with boundary conditions B(t,x,y,z,U,Ut,Ux,Uy,Uz)=0 and initial values U(t0,x,y,z)=U0 on a 3D domain bounded by right-angled polyhedrons. 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) VLUGR3.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 programs for the first two problems in the companion paper `VLUGR3: A Vectorizable Adaptive Grid Solver for PDEs in 3D. I. Algorithmic Aspects and Applications' 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 files EXMPL.f USER.f (only the SUBROUTINE DERIVF) VLUGR3.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_28 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 files EXMPLR.f USER.f (only (a dummy) FUNCTION INIDOM and SUBROUTINE CHSPCM) VLUGR3.f ILUBSn.f blas.f (if the BLAS library is not available on the platform) The results for this run is in the files EXMPLR_28 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: ILUBS3.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 '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 = 15) 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=3, NPTS=61000) PARAMETER (LENIWK=NPTS*(7*MXLEV+7), + LENRWK=NPTS*NPD*(5*MXLEV+13 + (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, YF, ZD, XR, YB, ZU, DX, DY, DZ, + TOLS, TOLT, RINFO(2+3*NPD), RWK(LENRWK) C First call of VLUGR3 MNTR = 0 NPDE = 3 T = 0.0 TOUT = 1.0 DT = 0.001 C Since domain is not a rectangular prism the grid parameters need not C to be specified here (cf. INIDOM) TOLS = 0.1 TOLT = 0.1 INFO(1) = 1 C MAXLEV INFO(2) = 4 C Domain not a rectangular prism INFO(3) = 1 C Linear system solver: matrix-free 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 RINFO(5) = 1.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 CALL VLUGR3 (NPDE, T, TOUT, DT, XL,YF,ZD, XR,YB,ZU, DX, DY, DZ, + TOLS, TOLT, INFO, RINFO, RWK, LENRWK, IWK, LENIWK, LWK, LENLWK, + MNTR) PRINT *, 'VLUGR3 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, YF, ZD, XR, YB, ZU, DX, DY, DZ, + LPLN, IPLN, LROW, IROW, ICOL, LLBND, ILBND, LBND) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER MAXPTS, LPLN(0:*), IPLN(*), LROW(*), IROW(*), ICOL(*), + LLBND(0:*), ILBND(*), LBND(*) DOUBLE PRECISION XL, YF, ZD, XR, YB, ZU, DX, DY, DZ 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 in any coordinate direction including C the corners. If Neumann boundaries are used the minimum is 4 since C otherwise the Jacobian matrix will be singular. C C A (virtual) box is placed upon the (irregular) domain. C The left/front/down point of this box is (XL,YF,ZD) in physical C coordinates and (0,0,0) in column, row, plane coordinates, resp.. C The right/back/upper point is (XR,YB,ZU) resp. (Nx,Ny,Nz), where C Nx = (XR-XL)/DX, Ny = (YB-YF)/DY, and Nz = (ZU-ZD)/DZ. C Only real grid points are stored. C The coordinate values of the initial grid should be stored plane C after plane and rowwise in LPLN, IPLN, 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 left/front/down point of virtual box C YF : OUT. Y-coordinate of left/front/down point of virtual box C ZD : OUT. Z-coordinate of left/front/down point of virtual box C XR : OUT. X-coordinate of right/back/upper point of virtual box C YB : OUT. Y-coordinate of right/back/upper point of virtual box C ZU : OUT. Z-coordinate of right/back/upper point of virtual box C DX : OUT. Grid width in X-direction C DY : OUT. Grid width in Y-direction C DZ : OUT. Grid width in Z-direction C LPLN : OUT. (0:LPLN(0)+1) C LPLN(0) = NPLNS: Actual # horizontal planes in LROW C LPLN(1:NPLNS): pointers to the start of a plane in LROW C LPLN(NPLNS+1) = NROWS+1: Total # rows in grid + 1 C IPLN : OUT. (NPLNS) C IPLN(IP): plane number of plane IP in virtual box C LROW : OUT. (NROWS+1) 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 : OUT. (NROWS) C IROW(IR): row number of row IR in virtual box C ICOL : OUT. (NPTS) C ICOL(IPT): column number of grid point IPT in virtual box C LLBND : OUT. (0:LLBND(0)+2) C LLBND(0) = NBNDS: total # physical planes in actual domain. C NB. edges and corners are stored for each plane they C belong to. C LLBND(1:NBNDS): pointers to a specific boundary in 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 : OUT. (NBNDS) C ILBND(IB): type of boundary: C 1: Left plane -I C 2: Down plane I C 3: Right plane I max. first order derivative C 4: Up plane I C 5: Front plane I C 6: Back plane -I C LBND : OUT. (NBIPTS) C LBND(IBPT): pointer to boundary point in actual grid C structure C Ccc EXTERNALS USED: NONE C C----------------------------------------------------------------------- C C Domain [0,1]x[0,1]x[0,1] with hole in the middle and projection C at [1,1.333]x[0,1]x[0.666,1]. C Virtual box: [0,1.333]x[0,1]x[0,1] C INTEGER NX, NY, NZ PARAMETER (NX = 8, NY = 6, NZ = 6) INTEGER IDOM(0:(NX+1)*(NY+1)*(NZ+1)) C INTEGER I, IPT, IR, J, K, NPLNS, NROWS, NPTS, NBNDS, + NPTSP1, NPTSP2 NPLNS = NZ+1 NROWS = (NY+1)*NPLNS NPTS = (NX+1)*NROWS-1-2*(NY+1)*(NZ-2) IF (MAXPTS .LT. NPTS) THEN INIDOM = .FALSE. MAXPTS = NPTS RETURN ELSE INIDOM = .TRUE. ENDIF XL = 0.0 YF = 0.0 ZD = 0.0 XR = 4.0/3.0 YB = 1.0 ZU = 1.0 DX = (XR-XL)/NX DY = (YB-YF)/NY DZ = (ZU-ZD)/NZ C C Make grid structure LPLN(0) = NPLNS IPT = 1 IR = 1 DO 10 K = 0, 2 LPLN(K+1) = IR IPLN(K+1) = K DO 20 I = 0, NY LROW(IR) = IPT IROW(IR) = I IR = IR+1 DO 30 J = 0, NX-2 ICOL(IPT) = J IPT = IPT + 1 30 CONTINUE 20 CONTINUE 10 CONTINUE K = 3 LPLN(K+1) = IR IPLN(K+1) = K DO 40 I = 0, 2 LROW(IR) = IPT IROW(IR) = I IR = IR+1 DO 50 J = 0, NX-2 ICOL(IPT) = J IPT = IPT + 1 50 CONTINUE 40 CONTINUE I = 3 LROW(IR) = IPT IROW(IR) = I IR = IR+1 DO 60 J = 0, 2 ICOL(IPT) = J IPT = IPT + 1 60 CONTINUE DO 70 J = 4, NX-2 ICOL(IPT) = J IPT = IPT + 1 70 CONTINUE DO 80 I = 4, NY LROW(IR) = IPT IROW(IR) = I IR = IR+1 DO 90 J = 0, NX-2 ICOL(IPT) = J IPT = IPT + 1 90 CONTINUE 80 CONTINUE DO 100 K = 4, NZ LPLN(K+1) = IR IPLN(K+1) = K DO 110 I = 0, NY LROW(IR) = IPT IROW(IR) = I IR = IR+1 DO 120 J = 0, NX ICOL(IPT) = J IPT = IPT + 1 120 CONTINUE 110 CONTINUE 100 CONTINUE LROW(NROWS+1) = NPTS+1 LPLN(NPLNS+1) = NROWS+1 C Ccc Boundaries NBNDS = 14 ILBND( 1) = 1 ILBND( 2) = 2 ILBND( 3) = 3 ILBND( 4) = 2 ILBND( 5) = 3 ILBND( 6) = 4 ILBND( 7) = 5 ILBND( 8) = 6 ILBND( 9) = 1 ILBND(10) = 2 ILBND(11) = 3 ILBND(12) = 4 ILBND(13) = 5 ILBND(14) = 6 LLBND( 0) = NBNDS LLBND( 1) = 1 LLBND( 2) = LLBND( 1) + (NY+1)*(NZ+1) LLBND( 3) = LLBND( 2) + (NX-1)*(NY+1) LLBND( 4) = LLBND( 3) + (NY+1)*(NZ-1) LLBND( 5) = LLBND( 4) + 3 *(NY+1) LLBND( 6) = LLBND( 5) + (NY+1)* 3 LLBND( 7) = LLBND( 6) + (NX+1)*(NY+1) LLBND( 8) = LLBND( 7) + (NX-1)*(NZ-2)+(NX+1)*3 LLBND( 9) = LLBND( 8) + (NX-1)*(NZ-2)+(NX+1)*3 LLBND(10) = LLBND( 9) + 9 LLBND(11) = LLBND(10) + 9 LLBND(12) = LLBND(11) + 9 LLBND(13) = LLBND(12) + 9 LLBND(14) = LLBND(13) + 9 LLBND(15) = LLBND(14) + 9 C Ccc Outer planes C Left boundary plane pointers NPTSP1 = (NX-1)*(NY+1) NPTSP2 = (NX+1)*(NY+1) DO 200 K = 0, 3 DO 201 I = 0, NY IPT = K*NPTSP1 + I*(NX-1) + 1 IF (K .EQ. 3 .AND. I .GT. 3) IPT = IPT-1 LBND(LLBND(1)+K*(NY+1)+I) = IPT 201 CONTINUE 200 CONTINUE DO 202 K = NZ-2, NZ DO 203 I = 0, NY IPT = (NZ-2)*NPTSP1+(K-NZ+2)*NPTSP2 + I*(NX+1) LBND(LLBND(1)+K*(NY+1)+I) = IPT 203 CONTINUE 202 CONTINUE C Right boundary plane pointers DO 210 K = 0, 3 DO 211 I = 0, NY IPT = (K+1)*NPTSP1 - I*(NX-1) IF (K .EQ. 3 .AND. I .LE. 3) IPT = IPT-1 LBND(LLBND(3)+K*(NY+1)+I) = IPT 211 CONTINUE 210 CONTINUE K = NZ-2 DO 209 I = 0, NY IPT = 4*NPTSP1 + NPTSP2-3 - I*(NX+1) LBND(LLBND(3)+K*(NY+1)+I) = IPT 209 CONTINUE DO 212 I = 0, NY DO 213 J = NX-2, NX IPT = NPTSP1*(NZ-2)+I*(NX+1) + J LBND(LLBND(4)+I*3+J-NX+2) = IPT 213 CONTINUE 212 CONTINUE DO 214 K = NZ-2, NZ DO 215 I = 0, NY IPT = NPTSP1*(NZ-2)+(K-NZ+3)*NPTSP2 - I*(NX+1) - 1 LBND(LLBND(5)+(K-NZ+2)*(NY+1)+I) = IPT 215 CONTINUE 214 CONTINUE C Down and up boundary plane pointers DO 220 I = 0, NY DO 221 J = 0, NX-2 IPT = I*(NX-1) + J + 1 LBND(LLBND(2)+I*(NX-1)+J) = IPT 221 CONTINUE 220 CONTINUE DO 230 I = 0, NY DO 231 J = 0, NX IPT = (NPTS - (I*(NX+1)+J)) LBND(LLBND(6)+I*(NX+1)+J) = IPT 231 CONTINUE 230 CONTINUE C Front and back boundary plane pointers DO 240 K = 0, 3 DO 241 J = 0, NX-2 IPT = K*NPTSP1 + J + 1 LBND(LLBND(7)+K*(NX-1)+J) = IPT IPT = (K+1)*NPTSP1 - J IF (K .EQ. 3) IPT = IPT-1 LBND(LLBND(8)+K*(NX-1)+J) = IPT 241 CONTINUE 240 CONTINUE DO 242 K = NZ-2, NZ DO 243 J = 0, NX IPT = (NZ-2)*NPTSP1+(K-NZ+2)*NPTSP2 + J LBND(LLBND(7)+4*(NX-1)+(K-NZ+2)*(NX+1)+J) = IPT IPT = 4*NPTSP1+(K-NZ+3)*NPTSP2 - J - 1 LBND(LLBND(8)+4*(NX-1)+(K-NZ+2)*(NX+1)+J) = IPT 243 CONTINUE 242 CONTINUE C Ccc Inner planes C Left and right boundary plane pointers LBND(LLBND( 9) ) = 117 LBND(LLBND( 9)+1) = 124 LBND(LLBND( 9)+2) = 131 LBND(LLBND( 9)+3) = 166 LBND(LLBND( 9)+4) = 172 LBND(LLBND( 9)+5) = 179 LBND(LLBND( 9)+6) = 218 LBND(LLBND( 9)+7) = 227 LBND(LLBND( 9)+8) = 236 LBND(LLBND(11) ) = 115 LBND(LLBND(11)+1) = 122 LBND(LLBND(11)+2) = 129 LBND(LLBND(11)+3) = 164 LBND(LLBND(11)+4) = 171 LBND(LLBND(11)+5) = 177 LBND(LLBND(11)+6) = 216 LBND(LLBND(11)+7) = 225 LBND(LLBND(11)+8) = 234 C Down and up boundary plane pointers DO 260 I = 0, 2 DO 270 J = 0, 2 LBND(LLBND(10)+I*3+J) = 236 - (I*(NX+1)+J) LBND(LLBND(12)+I*3+J) = 115 + I*(NX-1) + J 270 CONTINUE 260 CONTINUE C Front and back boundary plane pointers LBND(LLBND(13) ) = 129 LBND(LLBND(13)+1) = 130 LBND(LLBND(13)+2) = 131 LBND(LLBND(13)+3) = 177 LBND(LLBND(13)+4) = 178 LBND(LLBND(13)+5) = 179 LBND(LLBND(13)+6) = 234 LBND(LLBND(13)+7) = 235 LBND(LLBND(13)+8) = 236 LBND(LLBND(14) ) = 115 LBND(LLBND(14)+1) = 116 LBND(LLBND(14)+2) = 117 LBND(LLBND(14)+3) = 164 LBND(LLBND(14)+4) = 165 LBND(LLBND(14)+5) = 166 LBND(LLBND(14)+6) = 216 LBND(LLBND(14)+7) = 217 LBND(LLBND(14)+8) = 218 C LLBND(NBNDS+2) = LLBND(NBNDS+1) PRINT *, 'Input domain:' CALL PRDOM (LPLN, IPLN, LROW, IROW, ICOL, LLBND, ILBND, LBND, + IDOM, NX, NY, NZ) RETURN END SUBROUTINE CHSPCM (T, LEVEL, NPTS, X, Y, Z, NPDE, U, SPCMON, TOL) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER LEVEL, NPTS, NPDE DOUBLE PRECISION T, X(NPTS), Y(NPTS), Z(NPTS), U(NPTS,NPDE), + SPCMON(NPTS), TOL C Ccc PURPOSE: C Force grid refinement. C If for a node IPT SPCMON(IPT) > TOL the 64 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,Y,Z : IN. Arrays of physical 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 VLUGR3 C OUT: Changed to a value > TOL where refinement is required C TOL : IN. Tolerance with which SPCMON will be compared 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.5) .LT. 0.0001 .AND. + ABS(Z(I)) .LT. 0.0001) THEN SPCMON(I) = 2*TOL ENDIF 10 CONTINUE C RETURN END SUBROUTINE MONITR (T, DT, DTNEW, XL, YF, ZD, DXB, DYB, DZB, + LGRID, ISTRUC, LSOL, SOL) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER LGRID(0:*), ISTRUC(*), LSOL(*) DOUBLE PRECISION T, DT, DTNEW, XL, YF, ZD, DXB, DYB, DZB, 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 left/front/down point of (virtual) box C YF : IN. Y-coordinate of left/front/down point of (virtual) box C ZD : IN. Z-coordinate of left/front/down point of (virtual) box C DXB : IN. Cell width in X-direction of base grid C DYB : IN. Cell width in Y-direction of base grid C DZB : IN. Cell width in Z-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 C (LPLN,IPLN,LROW,IROW,ICOL) of refinement C level LEVEL for time T C ISTRUC : IN. (*) C ISTRUC(LGRID(LEVEL):.) contains (LPLN,IPLN,LROW,IROW,ICOL) C of grid level LEVEL, C LPLN : (0:LPLN(0)+1) C LPLN(0) = NPLNS: Actual # planes in LROW C LPLN(1:NPLNS): pointers to the start of a plane in LROW C LPLN(NPLNS+1) = NROWS+1: Total # rows in grid + 1 C IPLN : (NPLNS) C IPLN(IP): plane number of plane IP in virtual box C LROW : (NROWS+1) 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 box C ICOL : (NPTS) C ICOL(IPT): column number of grid point IPT in virtual box 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=250000, NPDE=3) DOUBLE PRECISION X(MAXPTS), Y(MAXPTS), Z(MAXPTS), UEX(MAXPTS*NPDE) C C----------------------------------------------------------------------- C INTEGER MAXLEV, LEVEL, LLPLN, LIPLN, LLROW, LIROW, LICOL, + NPLNS, NROWS, NPTS DOUBLE PRECISION DX, DY, DZ 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 DZ = DZB DO 10 LEVEL = 1, MAXLEV LLPLN = LGRID(LEVEL) NPLNS = ISTRUC(LLPLN) NROWS = ISTRUC(LLPLN+NPLNS+1)-1 LIPLN = LLPLN+NPLNS+2 LLROW = LIPLN+NPLNS NPTS = ISTRUC(LLROW+NROWS)-1 IF (NPTS .GT. MAXPTS) stop 'MONITR' LIROW = LLROW+NROWS+1 LICOL = LIROW+NROWS CALL SETXYZ (XL, YF, ZD, DX, DY, DZ, + ISTRUC(LLPLN), ISTRUC(LIPLN), + ISTRUC(LLROW), ISTRUC(LIROW), ISTRUC(LICOL), X, Y, Z) DX = DX/2 DY = DY/2 DZ = DZ/2 CALL PRERR (LEVEL, T, NPTS, NPDE, X, Y, Z, SOL(LSOL(LEVEL)+1), + UEX) 10 CONTINUE RETURN END SUBROUTINE PRERR (LEVEL, T, NPTS, NPDE, X, Y, Z, U, UEX) INTEGER LEVEL, NPTS, NPDE DOUBLE PRECISION T, X(NPTS), Y(NPTS), Z(NPTS), U(NPTS,NPDE), + UEX(NPTS,NPDE) INTEGER I,J DOUBLE PRECISION RMAX(3) CALL PDEIV (T, X, Y, Z, UEX, NPTS, NPDE) DO 1 J = 1,NPDE RMAX(J) = 0.0 DO 10 I = 1, NPTS RMAX(J) = MAX(RMAX(J),ABS(UEX(I,J)-U(I,J))) 10 CONTINUE 1 CONTINUE WRITE(28,'(''Error at T='',E9.3,'', level='',I1,'' :'', + I10,3E12.3)') + T, LEVEL, NPTS, (RMAX(J), J=1, NPDE) RETURN END SUBROUTINE PDEIV (T, X, Y, Z, U, NPTS, NPDE) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER NPTS, NPDE DOUBLE PRECISION T, X(NPTS), Y(NPTS), Z(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,Y,Z : IN. Arrays of physical 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 INTEGER I DOUBLE PRECISION EPS PARAMETER (EPS = 5.0D-3) DO 10 I = 1, NPTS U(I,1) = 1-0.5/(1+EXP((-X(I)+Y(I)+Z(I)-0.75*T)/(4*EPS))) U(I,2) = 1.5-U(I,1) U(I,3) = 1.5-U(I,1) 10 CONTINUE RETURN END SUBROUTINE PDEF (T, X, Y, Z, U, + UT, UX, UY, UZ, UXX, UYY, UZZ, UXY, UXZ, UYZ, RES, NPTS, NPDE) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER NPTS, NPDE DOUBLE PRECISION T, X(NPTS), Y(NPTS), Z(NPTS), U(NPTS,NPDE), + UT(NPTS,NPDE), UX(NPTS,NPDE), UY(NPTS,NPDE), UZ(NPTS,NPDE), + UXX(NPTS,NPDE), UYY(NPTS,NPDE), UZZ(NPTS,NPDE), + UXY(NPTS,NPDE), UXZ(NPTS,NPDE), UYZ(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,Y,Z : IN. Physical coordinates of 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 UZ : IN. I C UXX : IN. I C UYY : IN. I Space derivatives of U on current grid C UZZ : IN. I C UXY : IN. I C UXZ : IN. I C UYZ : 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 INTEGER I DOUBLE PRECISION EPS PARAMETER (EPS = 5.0D-3) DO 10 I = 1, NPTS RES(I,1) = UT(I,1) + U(I,1)*UX(I,1) + + U(I,2)*UY(I,1) + U(I,3)*UZ(I,1) - + EPS*(UXX(I,1)+UYY(I,1)+UZZ(I,1)) RES(I,2) = UT(I,2) + U(I,1)*UX(I,2) + + U(I,2)*UY(I,2) + U(I,3)*UZ(I,2) - + EPS*(UXX(I,2)+UYY(I,2)+UZZ(I,2)) RES(I,3) = UT(I,3) + U(I,1)*UX(I,3) + + U(I,2)*UY(I,3) + U(I,3)*UZ(I,3) - + EPS*(UXX(I,3)+UYY(I,3)+UZZ(I,3)) 10 CONTINUE RETURN END SUBROUTINE PDEBC (T, X, Y, Z, U, UT, UX, UY, UZ, 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), Z(NPTS), U(NPTS,NPDE), + UT(NPTS,NPDE), UX(NPTS,NPDE), UY(NPTS,NPDE), UZ(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,Y,Z : IN. Physical coordinates of 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 UZ : IN. -I 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)+2) C LLBND(0) = NBNDS: total # physical planes in actual domain. C NB. edges and corners are stored for each plane they C belong to. C LLBND(1:NBNDS): pointers to a specific boundary in 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 : IN. (NBNDS) C ILBND(IB): type of boundary: C 1: Left plane -I C 2: Down plane I C 3: Right plane I max. first order derivative C 4: Up plane I C 5: Front plane I C 6: Back plane -I C LBND : IN. (NBIPTS) C LBND(IBPT): pointer to boundary point in actual grid C C----------------------------------------------------------------------- C INTEGER I, J, K, NBNDS DOUBLE PRECISION EPS, UI PARAMETER (EPS = 5.0D-3) NBNDS = LLBND(0) DO 10 J = 1, NBNDS DO 20 K = LLBND(J), LLBND(J+1)-1 I = LBND(K) UI = 1-0.5/(1+EXP((-X(I)+Y(I)+Z(I)-0.75*T)/(4*EPS))) RES(I,1) = U(I,1) - UI RES(I,2) = U(I,2) - (1.5-UI) RES(I,3) = U(I,3) - (1.5-UI) 20 CONTINUE 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 INTEGER MXLEV, NPD, NPTS, LENIWK, LENRWK, LENLWK PARAMETER (MXLEV=2, NPD=3, NPTS=60000) PARAMETER (LENIWK=NPTS*(7*MXLEV+25), + LENRWK=NPTS*NPD*(5*MXLEV+13+38*NPD), + 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, YF, ZD, XR, YB, ZU, DX, DY, DZ, + TOLS, TOLT, RINFO(2+3*NPD), RWK(LENRWK) C Continuation call of VLUGR3 MNTR = 1 TOUT = 2.0 TOLS = 0.1 TOLT = 0.1 C Default choices INFO(1) = 0 C C Read info from file OPEN(UNIT=LUNDMP,FILE='DUMP',FORM='UNFORMATTED') CALL RDDUMP (LUNDMP, RWK, LENRWK, IWK, LENIWK) CLOSE(LUNDMP) C C Call main routine CALL VLUGR3 (NPDE, T, TOUT, DT, XL,YF,ZD, XR,YB,ZU, DX, DY, DZ, + TOLS, TOLT, INFO, RINFO, RWK, LENRWK, IWK, LENIWK, LWK, LENLWK, + MNTR) PRINT *, 'VLUGR3 returned with MNTR=', MNTR C C Save info on file OPEN(UNIT=LUNDMP,FILE='DUMP2',FORM='UNFORMATTED') CALL DUMP (LUNDMP, RWK, IWK) CLOSE(LUNDMP) END SUBROUTINE DERIVF (F, T, X, Y, Z, NPTS, NPDE, U, + A0, DT, DX, DY, DZ, + LLBND, ILBND, LBND, UIB, UT, UX, UY, UZ, UXX, UYY, UZZ, + UXY, UXZ, UYZ, ATOL, DEL, WORK, + FU, FUX, FUY, FUZ, FUXX, FUYY, FUZZ, FUXY, FUXZ, FUYZ) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER NPTS, NPDE, LLBND(0:*), ILBND(*), LBND(*) DOUBLE PRECISION F(NPTS,NPDE), T, X(NPTS), Y(NPTS), Z(NPTS), + U(NPTS,NPDE), + A0, DT, DX, DY, DZ, UIB(*), + UT(NPTS,NPDE), UX(NPTS,NPDE), UY(NPTS,NPDE), UZ(NPTS,NPDE), + UXX(NPTS,NPDE), UYY(NPTS,NPDE), UZZ(NPTS,NPDE), + UXY(NPTS,NPDE), UXZ(NPTS,NPDE), UYZ(NPTS,NPDE), + ATOL(NPDE), DEL(NPTS), WORK(2*NPTS*NPDE), + FU(NPTS,NPDE,NPDE), + FUX(NPTS,NPDE,NPDE), FUY(NPTS,NPDE,NPDE), FUZ(NPTS,NPDE,NPDE), + FUXX(NPTS,NPDE,NPDE),FUYY(NPTS,NPDE,NPDE),FUZZ(NPTS,NPDE,NPDE), + FUXY(NPTS,NPDE,NPDE),FUXZ(NPTS,NPDE,NPDE),FUYZ(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,Z : 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 DZ : IN. Cell width in Z-direction for current grid C LLBND : IN. (0:LLBND(0)+2) C LLBND(0) = NBNDS: total # physical planes in actual domain. C NB. edges and corners are stored for each plane they C belong to. C LLBND(1:NBNDS): pointers to a specific boundary in 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 : IN. (NBNDS) C ILBND(IB): type of boundary: C 1: Left plane -I C 2: Down plane I C 3: Right plane I max. first order derivative C 4: Up plane I C 5: Front plane I C 6: Back plane -I 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 UZ : IN. I C UXX : IN. I C UYY : IN. I Space derivatives of U on current grid C UZZ : IN. I C UXY : IN. I C UXZ : IN. I C UYZ : IN. -I C ATOL : 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 FUZ : OUT. dF(Uz)dUz C FUXX : OUT. dF(Uxx)dUxx C FUYY : OUT. dF(Uyy)dUyy C FUZZ : OUT. dF(Uzz)dUzz C FUXY : OUT. dF(Uxy)dUxy C FUXZ : OUT. dF(Uxz)dUxz C FUYZ : OUT. dF(Uyz)dUyz C Ccc EXTERNALS USED: EXTERNAL ZERO C C----------------------------------------------------------------------- C DOUBLE PRECISION EPS PARAMETER (EPS = 5D-3) C INTEGER IC, IPT, JC, LB, NBNDS C CALL ZERO (NPTS*NPDE*NPDE, FUX) CALL ZERO (NPTS*NPDE*NPDE, FUY) CALL ZERO (NPTS*NPDE*NPDE, FUZ) CALL ZERO (NPTS*NPDE*NPDE, FUXX) CALL ZERO (NPTS*NPDE*NPDE, FUYY) CALL ZERO (NPTS*NPDE*NPDE, FUZZ) CALL ZERO (NPTS*NPDE*NPDE, FUXY) CALL ZERO (NPTS*NPDE*NPDE, FUXZ) CALL ZERO (NPTS*NPDE*NPDE, FUYZ) C DO 10 IPT = 1, NPTS C dF_1(U,Ut)/dU_ic FU(IPT,1,1) = UX(IPT,1) + A0 FU(IPT,1,2) = UY(IPT,1) FU(IPT,1,3) = UZ(IPT,1) C dF_1(Up)/dUp_ic FUX(IPT,1,1) = U(IPT,1) FUY(IPT,1,1) = U(IPT,2) FUZ(IPT,1,1) = U(IPT,3) C dF_1(Upp)/dUpp_ic FUXX(IPT,1,1) = -EPS FUYY(IPT,1,1) = -EPS FUZZ(IPT,1,1) = -EPS C dF_2(U,Ut)/dU_ic FU(IPT,2,1) = UX(IPT,2) FU(IPT,2,2) = UY(IPT,2) + A0 FU(IPT,2,3) = UZ(IPT,2) C dF_2(Up)/dUp_ic FUX(IPT,2,2) = U(IPT,1) FUY(IPT,2,2) = U(IPT,2) FUZ(IPT,2,2) = U(IPT,3) C dF_2(Upp)/dUpp_ic FUXX(IPT,2,2) = -EPS FUYY(IPT,2,2) = -EPS FUZZ(IPT,2,2) = -EPS C dF_3(U,Ut)/dU_ic FU(IPT,3,1) = UX(IPT,3) FU(IPT,3,2) = UY(IPT,3) FU(IPT,3,3) = UZ(IPT,3) + A0 C dF_3(Up)/dUp_ic FUX(IPT,3,3) = U(IPT,1) FUY(IPT,3,3) = U(IPT,2) FUZ(IPT,3,3) = U(IPT,3) C dF_3(Upp)/dUpp_ic FUXX(IPT,3,3) = -EPS FUYY(IPT,3,3) = -EPS FUZZ(IPT,3,3) = -EPS 10 CONTINUE C C Correct boundaries (incl. the internal) NBNDS = LLBND(0) DO 100 LB = LLBND(1), LLBND(NBNDS+2)-1 IPT = LBND(LB) DO 110 IC = 1, NPDE DO 120 JC = 1, NPDE FU (IPT,IC,JC) = 0.0 FUX (IPT,IC,JC) = 0.0 FUY (IPT,IC,JC) = 0.0 FUZ (IPT,IC,JC) = 0.0 FUXX(IPT,IC,JC) = 0.0 FUYY(IPT,IC,JC) = 0.0 FUZZ(IPT,IC,JC) = 0.0 120 CONTINUE FU(IPT,IC,IC) = 1.0 110 CONTINUE 100 CONTINUE RETURN END SUBROUTINE MONITR (T, DT, DTNEW, XL, YF, ZD, DXB, DYB, DZB, + LGRID, ISTRUC, LSOL, SOL) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER LGRID(0:*), ISTRUC(*), LSOL(*) DOUBLE PRECISION T, DT, DTNEW, XL, YF, ZD, DXB, DYB, DZB, 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 left/front/down point of (virtual) box C YF : IN. Y-coordinate of left/front/down point of (virtual) box C ZD : IN. Z-coordinate of left/front/down point of (virtual) box C DXB : IN. Cell width in X-direction of base grid C DYB : IN. Cell width in Y-direction of base grid C DZB : IN. Cell width in Z-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 C (LPLN,IPLN,LROW,IROW,ICOL) of refinement C level LEVEL for time T C ISTRUC : IN. (*) C ISTRUC(LGRID(LEVEL):.) contains (LPLN,IPLN,LROW,IROW,ICOL) C of grid level LEVEL, C LPLN : (0:LPLN(0)+1) C LPLN(0) = NPLNS: Actual # planes in LROW C LPLN(1:NPLNS): pointers to the start of a plane in LROW C LPLN(NPLNS+1) = NROWS+1: Total # rows in grid + 1 C IPLN : (NPLNS) C IPLN(IP): plane number of plane IP in virtual box C LROW : (NROWS+1) 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 box C ICOL : (NPTS) C ICOL(IPT): column number of grid point IPT in virtual box 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=75000, NPDE=3) DOUBLE PRECISION X(MAXPTS), Y(MAXPTS), Z(MAXPTS), UEX(MAXPTS*NPDE) C C----------------------------------------------------------------------- C INTEGER MAXLEV, LEVEL, LLPLN, LIPLN, LLROW, LIROW, LICOL, + NPLNS, NROWS, NPTS DOUBLE PRECISION DX, DY, DZ 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 DZ = DZB DO 10 LEVEL = 1, MAXLEV LLPLN = LGRID(LEVEL) NPLNS = ISTRUC(LLPLN) NROWS = ISTRUC(LLPLN+NPLNS+1)-1 LIPLN = LLPLN+NPLNS+2 LLROW = LIPLN+NPLNS NPTS = ISTRUC(LLROW+NROWS)-1 IF (NPTS .GT. MAXPTS) stop 'MONITR' LIROW = LLROW+NROWS+1 LICOL = LIROW+NROWS CALL SETXYZ (XL, YF, ZD, DX, DY, DZ, + ISTRUC(LLPLN), ISTRUC(LIPLN), + ISTRUC(LLROW), ISTRUC(LIROW), ISTRUC(LICOL), X, Y, Z) DX = DX/2 DY = DY/2 DZ = DZ/2 CALL PRERR (LEVEL, T, NPTS, NPDE, X, Y, Z, SOL(LSOL(LEVEL)+1), + UEX) 10 CONTINUE RETURN END SUBROUTINE PRERR (LEVEL, T, NPTS, NPDE, X, Y, Z, U, UEX) INTEGER LEVEL, NPTS, NPDE DOUBLE PRECISION T, X(NPTS), Y(NPTS), Z(NPTS), U(NPTS,NPDE), + UEX(NPTS,NPDE) INTEGER I,J DOUBLE PRECISION RMAX(3) CALL PDEIV (T, X, Y, Z, UEX, NPTS, NPDE) DO 1 J = 1, NPDE RMAX(J) = 0.0 DO 10 I = 1, NPTS RMAX(J) = MAX(RMAX(J),ABS(UEX(I,J)-U(I,J))) 10 CONTINUE 1 CONTINUE WRITE(28,'(''Error at T='',E9.3,'', level='',I1,'' :'', + I10,3E12.3)') + T, LEVEL, NPTS, (RMAX(J), J=1, NPDE) RETURN END SUBROUTINE PDEIV (T, X, Y, Z, U, NPTS, NPDE) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER NPTS, NPDE DOUBLE PRECISION T, X(NPTS), Y(NPTS), Z(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,Y,Z : IN. Arrays of physical 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 INTEGER I DOUBLE PRECISION EPS PARAMETER (EPS = 5.0D-3) DO 10 I = 1, NPTS U(I,1) = 1-0.5/(1+EXP((-X(I)+Y(I)+Z(I)-0.75*T)/(4*EPS))) U(I,2) = 1.5-U(I,1) U(I,3) = 1.5-U(I,1) 10 CONTINUE RETURN END SUBROUTINE PDEF (T, X, Y, Z, U, + UT, UX, UY, UZ, UXX, UYY, UZZ, UXY, UXZ, UYZ, RES, NPTS, NPDE) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER NPTS, NPDE DOUBLE PRECISION T, X(NPTS), Y(NPTS), Z(NPTS), U(NPTS,NPDE), + UT(NPTS,NPDE), UX(NPTS,NPDE), UY(NPTS,NPDE), UZ(NPTS,NPDE), + UXX(NPTS,NPDE), UYY(NPTS,NPDE), UZZ(NPTS,NPDE), + UXY(NPTS,NPDE), UXZ(NPTS,NPDE), UYZ(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,Y,Z : IN. Physical coordinates of 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 UZ : IN. I C UXX : IN. I C UYY : IN. I Space derivatives of U on current grid C UZZ : IN. I C UXY : IN. I C UXZ : IN. I C UYZ : 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 INTEGER I DOUBLE PRECISION EPS PARAMETER (EPS = 5.0D-3) DO 10 I = 1, NPTS RES(I,1) = UT(I,1) + U(I,1)*UX(I,1) + + U(I,2)*UY(I,1) + U(I,3)*UZ(I,1) - + EPS*(UXX(I,1)+UYY(I,1)+UZZ(I,1)) RES(I,2) = UT(I,2) + U(I,1)*UX(I,2) + + U(I,2)*UY(I,2) + U(I,3)*UZ(I,2) - + EPS*(UXX(I,2)+UYY(I,2)+UZZ(I,2)) RES(I,3) = UT(I,3) + U(I,1)*UX(I,3) + + U(I,2)*UY(I,3) + U(I,3)*UZ(I,3) - + EPS*(UXX(I,3)+UYY(I,3)+UZZ(I,3)) 10 CONTINUE RETURN END SUBROUTINE PDEBC (T, X, Y, Z, U, UT, UX, UY, UZ, 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), Z(NPTS), U(NPTS,NPDE), + UT(NPTS,NPDE), UX(NPTS,NPDE), UY(NPTS,NPDE), UZ(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,Y,Z : IN. Physical coordinates of 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 UZ : IN. -I 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)+2) C LLBND(0) = NBNDS: total # physical planes in actual domain. C NB. edges and corners are stored for each plane they C belong to. C LLBND(1:NBNDS): pointers to a specific boundary in 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 : IN. (NBNDS) C ILBND(IB): type of boundary: C 1: Left plane -I C 2: Down plane I C 3: Right plane I max. first order derivative C 4: Up plane I C 5: Front plane I C 6: Back plane -I C LBND : IN. (NBIPTS) C LBND(IBPT): pointer to boundary point in actual grid C C----------------------------------------------------------------------- C INTEGER I, J, K, NBNDS DOUBLE PRECISION EPS, UI PARAMETER (EPS = 5.0D-3) NBNDS = LLBND(0) DO 10 J = 1, NBNDS DO 20 K = LLBND(J), LLBND(J+1)-1 I = LBND(K) UI = 1-0.5/(1+EXP((-X(I)+Y(I)+Z(I)-0.75*T)/(4*EPS))) RES(I,1) = U(I,1) - UI RES(I,2) = U(I,2) - (1.5-UI) RES(I,3) = U(I,3) - (1.5-UI) 20 CONTINUE 10 CONTINUE RETURN END SHAR_EOF fi # end of overwriting check if test -f 'probia.f' then echo shar: will not over-write existing file "'probia.f'" else cat << \SHAR_EOF > 'probia.f' PROGRAM EXIA C INTEGER MXLEV, NPD, NPTS, LENIWK, LENRWK, LENLWK PARAMETER (MXLEV=2, NPD=1, NPTS=170000) PARAMETER (LENIWK=NPTS*(7*MXLEV+25), + LENRWK=NPTS*NPD*(5*MXLEV+38*NPD+13), + 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, YF, ZD, XR, YB, ZU, DX, DY, DZ, + TOLS, TOLT, RINFO(2+3*NPD), RWK(LENRWK) C First call of VLUGR3 MNTR = 0 NPDE = 1 T = 0.0 TOUT = 1.0 DT = 0.001 XL = 0.0 YF = 0.0 ZD = 0.0 XR = 1.0 YB = 1.0 ZU = 1.0 DX = 0.1 DY = 0.1 DZ = 0.1 TOLS = 0.1 TOLT = 0.1 INFO(1) = 1 C MAXLEV INFO(2) = 4 C Domain a rectangular prism INFO(3) = 0 C Linear system solver print *, 'Lin.sys.solver BiCGStab or GCRO (0 / 10,11,12,13) ?' 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 solver 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 C SPCWGT = 1.0 RINFO(4) = 1.0 C TIMWGT = 1.0 RINFO(5) = 1.0 C C Call main routine CALL VLUGR3 (NPDE, T, TOUT, DT, XL,YF,ZD, XR,YB,ZU, DX, DY, DZ, + TOLS, TOLT, INFO, RINFO, RWK, LENRWK, IWK, LENIWK, LWK, LENLWK, + MNTR) PRINT *, 'VLUGR3 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 SUBROUTINE MONITR (T, DT, DTNEW, XL, YF, ZD, DXB, DYB, DZB, + LGRID, ISTRUC, LSOL, SOL) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER LGRID(0:*), ISTRUC(*), LSOL(*) DOUBLE PRECISION T, DT, DTNEW, XL, YF, ZD, DXB, DYB, DZB, 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 left/front/down point of (virtual) box C YF : IN. Y-coordinate of left/front/down point of (virtual) box C ZD : IN. Z-coordinate of left/front/down point of (virtual) box C DXB : IN. Cell width in X-direction of base grid C DYB : IN. Cell width in Y-direction of base grid C DZB : IN. Cell width in Z-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 C (LPLN,IPLN,LROW,IROW,ICOL) of refinement C level LEVEL for time T C ISTRUC : IN. (*) C ISTRUC(LGRID(LEVEL):.) contains (LPLN,IPLN,LROW,IROW,ICOL) C of grid level LEVEL, C LPLN : (0:LPLN(0)+1) C LPLN(0) = NPLNS: Actual # planes in LROW C LPLN(1:NPLNS): pointers to the start of a plane in LROW C LPLN(NPLNS+1) = NROWS+1: Total # rows in grid + 1 C IPLN : (NPLNS) C IPLN(IP): plane number of plane IP in virtual box C LROW : (NROWS+1) 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 box C ICOL : (NPTS) C ICOL(IPT): column number of grid point IPT in virtual box 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=250000, NPDE=1) DOUBLE PRECISION X(MAXPTS), Y(MAXPTS), Z(MAXPTS), UEX(MAXPTS*NPDE) C C----------------------------------------------------------------------- C INTEGER MAXLEV, LEVEL, LLPLN, LIPLN, LLROW, LIROW, LICOL, + NPLNS, NROWS, NPTS DOUBLE PRECISION DX, DY, DZ 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 DZ = DZB DO 10 LEVEL = 1, MAXLEV LLPLN = LGRID(LEVEL) NPLNS = ISTRUC(LLPLN) NROWS = ISTRUC(LLPLN+NPLNS+1)-1 LIPLN = LLPLN+NPLNS+2 LLROW = LIPLN+NPLNS NPTS = ISTRUC(LLROW+NROWS)-1 IF (NPTS .GT. MAXPTS) stop 'MONITR' LIROW = LLROW+NROWS+1 LICOL = LIROW+NROWS CALL SETXYZ (XL, YF, ZD, DX, DY, DZ, + ISTRUC(LLPLN), ISTRUC(LIPLN), + ISTRUC(LLROW), ISTRUC(LIROW), ISTRUC(LICOL), X, Y, Z) DX = DX/2 DY = DY/2 DZ = DZ/2 CALL PRERR (LEVEL, T, NPTS, NPDE, X, Y, Z, SOL(LSOL(LEVEL)+1), + UEX) 10 CONTINUE RETURN END SUBROUTINE PRERR (LEVEL, T, NPTS, NPDE, X, Y, Z, U, UEX) INTEGER LEVEL, NPTS, NPDE DOUBLE PRECISION T, X(NPTS), Y(NPTS), Z(NPTS), U(NPTS,NPDE), + UEX(NPTS,NPDE) INTEGER I,J DOUBLE PRECISION RMAX CALL PDEIV (T, X, Y, Z, UEX, NPTS, NPDE) RMAX = 0.0 J = 1 DO 10 I = 1, NPTS RMAX = MAX(RMAX,ABS(UEX(I,J)-U(I,J))) 10 CONTINUE WRITE(28,'(''Error at T='',E9.3,'', level='',I1,'' :'',E9.3,I10)') + T, LEVEL, RMAX, NPTS RETURN END SUBROUTINE PDEIV (T, X, Y, Z, U, NPTS, NPDE) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER NPTS, NPDE DOUBLE PRECISION T, X(NPTS), Y(NPTS), Z(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,Y,Z : IN. Arrays of physical 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 DOUBLE PRECISION EPS PARAMETER (EPS = 2.0D-3) DO 10 I = 1, NPTS U(I,1) = 1-0.5/(1+EXP((-X(I)+Y(I)+Z(I)-0.75*T)/(4*EPS))) 10 CONTINUE RETURN END SUBROUTINE PDEF (T, X, Y, Z, U, + UT, UX, UY, UZ, UXX, UYY, UZZ, UXY, UXZ, UYZ, RES, NPTS, NPDE) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER NPTS, NPDE DOUBLE PRECISION T, X(NPTS), Y(NPTS), Z(NPTS), U(NPTS,NPDE), + UT(NPTS,NPDE), UX(NPTS,NPDE), UY(NPTS,NPDE), UZ(NPTS,NPDE), + UXX(NPTS,NPDE), UYY(NPTS,NPDE), UZZ(NPTS,NPDE), + UXY(NPTS,NPDE), UXZ(NPTS,NPDE), UYZ(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,Y,Z : IN. Physical coordinates of 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 UZ : IN. I C UXX : IN. I C UYY : IN. I Space derivatives of U on current grid C UZZ : IN. I C UXY : IN. I C UXZ : IN. I C UYZ : 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 INTEGER I DOUBLE PRECISION EPS PARAMETER (EPS = 2.0D-3) DO 10 I = 1, NPTS RES(I,1) = UT(I,1) + U(I,1)*UX(I,1) + + (1.5-U(I,1))*(UY(I,1)+UZ(I,1)) - + EPS*(UXX(I,1)+UYY(I,1)+UZZ(I,1)) 10 CONTINUE RETURN END SUBROUTINE PDEBC (T, X, Y, Z, U, UT, UX, UY, UZ, 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), Z(NPTS), U(NPTS,NPDE), + UT(NPTS,NPDE), UX(NPTS,NPDE), UY(NPTS,NPDE), UZ(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,Y,Z : IN. Physical coordinates of 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 UZ : IN. -I 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)+2) C LLBND(0) = NBNDS: total # physical planes in actual domain. C NB. edges and corners are stored for each plane they C belong to. C LLBND(1:NBNDS): pointers to a specific boundary in 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 : IN. (NBNDS) C ILBND(IB): type of boundary: C 1: Left plane -I C 2: Down plane I C 3: Right plane I max. first order derivative C 4: Up plane I C 5: Front plane I C 6: Back plane -I C LBND : IN. (NBIPTS) C LBND(IBPT): pointer to boundary point in actual grid C C----------------------------------------------------------------------- C INTEGER I, J, K, NBNDS DOUBLE PRECISION EPS, UI PARAMETER (EPS = 2.0D-3) NBNDS = LLBND(0) DO 10 J = 1, NBNDS DO 20 K = LLBND(J), LLBND(J+1)-1 I = LBND(K) UI = 1-0.5/(1+EXP((-X(I)+Y(I)+Z(I)-0.75*T)/(4*EPS))) RES(I,1) = U(I,1) - UI 20 CONTINUE 10 CONTINUE RETURN END SHAR_EOF fi # end of overwriting check if test -f 'probib.f' then echo shar: will not over-write existing file "'probib.f'" else cat << \SHAR_EOF > 'probib.f' PROGRAM EXIB C C PARAMETER (NRRMAX = 1, MAXLR = 3, MAXL = 15) INTEGER MXLEV, NPD, NPTS, LENIWK, LENRWK, LENLWK PARAMETER (MXLEV=1, NPD=3, NPTS=190000) PARAMETER (LENIWK=NPTS*(9*MXLEV+24), + LENRWK=NPTS*NPD*(5*MXLEV+26+13), + 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, YF, ZD, XR, YB, ZU, DX, DY, DZ, + TOLS, TOLT, RINFO(2+3*NPD), RWK(LENRWK) C First call of VLUGR3 MNTR = 0 NPDE = 3 T = 0.0 TOUT = 1.0 DT = 0.001 XL = 0.0 YF = 0.0 ZD = 0.0 XR = 1.0 YB = 1.0 ZU = 1.0 DX = 0.1 DY = 0.1 DZ = 0.1 TOLS = 0.1 TOLT = 0.1 INFO(1) = 1 C MAXLEV INFO(2) = 4 C Domain a rectangular prism INFO(3) = 0 C Linear system solver print *, 'Lin.sys.solver BiCGStab or GCRO (0 / 10,11,12,13) ?' 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 solver 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 RINFO(5) = 1.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 CALL VLUGR3 (NPDE, T, TOUT, DT, XL,YF,ZD, XR,YB,ZU, DX, DY, DZ, + TOLS, TOLT, INFO, RINFO, RWK, LENRWK, IWK, LENIWK, LWK, LENLWK, + MNTR) PRINT *, 'VLUGR3 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 SUBROUTINE MONITR (T, DT, DTNEW, XL, YF, ZD, DXB, DYB, DZB, + LGRID, ISTRUC, LSOL, SOL) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER LGRID(0:*), ISTRUC(*), LSOL(*) DOUBLE PRECISION T, DT, DTNEW, XL, YF, ZD, DXB, DYB, DZB, 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 left/front/down point of (virtual) box C YF : IN. Y-coordinate of left/front/down point of (virtual) box C ZD : IN. Z-coordinate of left/front/down point of (virtual) box C DXB : IN. Cell width in X-direction of base grid C DYB : IN. Cell width in Y-direction of base grid C DZB : IN. Cell width in Z-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 C (LPLN,IPLN,LROW,IROW,ICOL) of refinement C level LEVEL for time T C ISTRUC : IN. (*) C ISTRUC(LGRID(LEVEL):.) contains (LPLN,IPLN,LROW,IROW,ICOL) C of grid level LEVEL, C LPLN : (0:LPLN(0)+1) C LPLN(0) = NPLNS: Actual # planes in LROW C LPLN(1:NPLNS): pointers to the start of a plane in LROW C LPLN(NPLNS+1) = NROWS+1: Total # rows in grid + 1 C IPLN : (NPLNS) C IPLN(IP): plane number of plane IP in virtual box C LROW : (NROWS+1) 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 box C ICOL : (NPTS) C ICOL(IPT): column number of grid point IPT in virtual box 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=170000, NPDE=3) DOUBLE PRECISION X(MAXPTS), Y(MAXPTS), Z(MAXPTS), UEX(MAXPTS*NPDE) C C----------------------------------------------------------------------- C INTEGER MAXLEV, LEVEL, LLPLN, LIPLN, LLROW, LIROW, LICOL, + NPLNS, NROWS, NPTS DOUBLE PRECISION DX, DY, DZ 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 DZ = DZB DO 10 LEVEL = 1, MAXLEV LLPLN = LGRID(LEVEL) NPLNS = ISTRUC(LLPLN) NROWS = ISTRUC(LLPLN+NPLNS+1)-1 LIPLN = LLPLN+NPLNS+2 LLROW = LIPLN+NPLNS NPTS = ISTRUC(LLROW+NROWS)-1 IF (NPTS .GT. MAXPTS) stop 'MONITR' LIROW = LLROW+NROWS+1 LICOL = LIROW+NROWS CALL SETXYZ (XL, YF, ZD, DX, DY, DZ, + ISTRUC(LLPLN), ISTRUC(LIPLN), + ISTRUC(LLROW), ISTRUC(LIROW), ISTRUC(LICOL), X, Y, Z) DX = DX/2 DY = DY/2 DZ = DZ/2 CALL PRERR (LEVEL, T, NPTS, NPDE, X, Y, Z, SOL(LSOL(LEVEL)+1), + UEX) 10 CONTINUE RETURN END SUBROUTINE PRERR (LEVEL, T, NPTS, NPDE, X, Y, Z, U, UEX) INTEGER LEVEL, NPTS, NPDE DOUBLE PRECISION T, X(NPTS), Y(NPTS), Z(NPTS), U(NPTS,NPDE), + UEX(NPTS,NPDE) INTEGER I,J DOUBLE PRECISION RMAX(3) CALL PDEIV (T, X, Y, Z, UEX, NPTS, NPDE) DO 1 J = 1, NPDE RMAX(J) = 0.0 DO 10 I = 1, NPTS RMAX(J) = MAX(RMAX(J),ABS(UEX(I,J)-U(I,J))) 10 CONTINUE 1 CONTINUE WRITE(28,'(''Error at T='',E9.3,'', level='',I1,'' :'', + I10,3E10.3)') T, LEVEL, NPTS, (RMAX(J), J=1, NPDE) RETURN END SUBROUTINE PDEIV (T, X, Y, Z, U, NPTS, NPDE) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER NPTS, NPDE DOUBLE PRECISION T, X(NPTS), Y(NPTS), Z(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,Y,Z : IN. Arrays of physical 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 INTEGER I DOUBLE PRECISION EPS PARAMETER (EPS = 2.0D-3) DO 10 I = 1, NPTS U(I,1) = 1-0.5/(1+EXP((-X(I)+Y(I)+Z(I)-0.75*T)/(4*EPS))) U(I,2) = 1.5-U(I,1) U(I,3) = 1.5-U(I,1) 10 CONTINUE RETURN END SUBROUTINE PDEF (T, X, Y, Z, U, + UT, UX, UY, UZ, UXX, UYY, UZZ, UXY, UXZ, UYZ, RES, NPTS, NPDE) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER NPTS, NPDE DOUBLE PRECISION T, X(NPTS), Y(NPTS), Z(NPTS), U(NPTS,NPDE), + UT(NPTS,NPDE), UX(NPTS,NPDE), UY(NPTS,NPDE), UZ(NPTS,NPDE), + UXX(NPTS,NPDE), UYY(NPTS,NPDE), UZZ(NPTS,NPDE), + UXY(NPTS,NPDE), UXZ(NPTS,NPDE), UYZ(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,Y,Z : IN. Physical coordinates of 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 UZ : IN. I C UXX : IN. I C UYY : IN. I Space derivatives of U on current grid C UZZ : IN. I C UXY : IN. I C UXZ : IN. I C UYZ : 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 INTEGER I DOUBLE PRECISION EPS PARAMETER (EPS = 2.0D-3) DO 10 I = 1, NPTS RES(I,1) = UT(I,1) + U(I,1)*UX(I,1) + + U(I,2)*UY(I,1) + U(I,3)*UZ(I,1) - + EPS*(UXX(I,1)+UYY(I,1)+UZZ(I,1)) RES(I,2) = UT(I,2) + U(I,1)*UX(I,2) + + U(I,2)*UY(I,2) + U(I,3)*UZ(I,2) - + EPS*(UXX(I,2)+UYY(I,2)+UZZ(I,2)) RES(I,3) = UT(I,3) + U(I,1)*UX(I,3) + + U(I,2)*UY(I,3) + U(I,3)*UZ(I,3) - + EPS*(UXX(I,3)+UYY(I,3)+UZZ(I,3)) 10 CONTINUE RETURN END SUBROUTINE PDEBC (T, X, Y, Z, U, UT, UX, UY, UZ, 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), Z(NPTS), U(NPTS,NPDE), + UT(NPTS,NPDE), UX(NPTS,NPDE), UY(NPTS,NPDE), UZ(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,Y,Z : IN. Physical coordinates of 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 UZ : IN. -I 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)+2) C LLBND(0) = NBNDS: total # physical planes in actual domain. C NB. edges and corners are stored for each plane they C belong to. C LLBND(1:NBNDS): pointers to a specific boundary in 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 : IN. (NBNDS) C ILBND(IB): type of boundary: C 1: Left plane -I C 2: Down plane I C 3: Right plane I max. first order derivative C 4: Up plane I C 5: Front plane I C 6: Back plane -I C LBND : IN. (NBIPTS) C LBND(IBPT): pointer to boundary point in actual grid C C----------------------------------------------------------------------- C INTEGER I, J, K, NBNDS DOUBLE PRECISION EPS, UI PARAMETER (EPS = 2.0D-3) NBNDS = LLBND(0) DO 10 J = 1, NBNDS DO 20 K = LLBND(J), LLBND(J+1)-1 I = LBND(K) UI = 1-0.5/(1+EXP((-X(I)+Y(I)+Z(I)-0.75*T)/(4*EPS))) RES(I,1) = U(I,1) - UI RES(I,2) = U(I,2) - (1.5-UI) RES(I,3) = U(I,3) - (1.5-UI) 20 CONTINUE 10 CONTINUE RETURN END SHAR_EOF fi # end of overwriting check if test -f 'probic.f' then echo shar: will not over-write existing file "'probic.f'" else cat << \SHAR_EOF > 'probic.f' PROGRAM EXIC C INTEGER MXLEV, NPD, NPTS, LENIWK, LENRWK, LENLWK PARAMETER (MXLEV=2, NPD=3, NPTS=40000) PARAMETER (LENIWK=NPTS*(7*MXLEV+25), + LENRWK=NPTS*NPD*(5*MXLEV+38*NPD+13), + 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, YF, ZD, XR, YB, ZU, DX, DY, DZ, + TOLS, TOLT, RINFO(2+3*NPD), RWK(LENRWK) C First call of VLUGR3 MNTR = 0 NPDE = 3 T = 0.0 TOUT = 1.0 DT = 0.001 XL = 0.0 YF = 0.0 ZD = 0.0 XR = 1.0 YB = 1.0 ZU = 1.0 DX = 0.1 DY = 0.1 DZ = 0.1 TOLS = 0.1 TOLT = 0.1 INFO(1) = 1 C MAXLEV INFO(2) = 3 C Domain a rectangular prism INFO(3) = 0 C Linear system solver print *, 'Lin.sys.solver BiCGStab or GCRO (0 / 10,11,12,13) ?' 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 solver 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 RINFO(5) = 1.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 CALL VLUGR3 (NPDE, T, TOUT, DT, XL,YF,ZD, XR,YB,ZU, DX, DY, DZ, + TOLS, TOLT, INFO, RINFO, RWK, LENRWK, IWK, LENIWK, LWK, LENLWK, + MNTR) PRINT *, 'VLUGR3 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 SUBROUTINE MONITR (T, DT, DTNEW, XL, YF, ZD, DXB, DYB, DZB, + LGRID, ISTRUC, LSOL, SOL) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER LGRID(0:*), ISTRUC(*), LSOL(*) DOUBLE PRECISION T, DT, DTNEW, XL, YF, ZD, DXB, DYB, DZB, 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 left/front/down point of (virtual) box C YF : IN. Y-coordinate of left/front/down point of (virtual) box C ZD : IN. Z-coordinate of left/front/down point of (virtual) box C DXB : IN. Cell width in X-direction of base grid C DYB : IN. Cell width in Y-direction of base grid C DZB : IN. Cell width in Z-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 C (LPLN,IPLN,LROW,IROW,ICOL) of refinement C level LEVEL for time T C ISTRUC : IN. (*) C ISTRUC(LGRID(LEVEL):.) contains (LPLN,IPLN,LROW,IROW,ICOL) C of grid level LEVEL, C LPLN : (0:LPLN(0)+1) C LPLN(0) = NPLNS: Actual # planes in LROW C LPLN(1:NPLNS): pointers to the start of a plane in LROW C LPLN(NPLNS+1) = NROWS+1: Total # rows in grid + 1 C IPLN : (NPLNS) C IPLN(IP): plane number of plane IP in virtual box C LROW : (NROWS+1) 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 box C ICOL : (NPTS) C ICOL(IPT): column number of grid point IPT in virtual box 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=250000, NPDE=1) DOUBLE PRECISION X(MAXPTS), Y(MAXPTS), Z(MAXPTS), UEX(MAXPTS*NPDE) C C----------------------------------------------------------------------- C INTEGER MAXLEV, LEVEL, LLPLN, LIPLN, LLROW, LIROW, LICOL, + NPLNS, NROWS, NPTS DOUBLE PRECISION DX, DY, DZ 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 DZ = DZB DO 10 LEVEL = 1, MAXLEV LLPLN = LGRID(LEVEL) NPLNS = ISTRUC(LLPLN) NROWS = ISTRUC(LLPLN+NPLNS+1)-1 LIPLN = LLPLN+NPLNS+2 LLROW = LIPLN+NPLNS NPTS = ISTRUC(LLROW+NROWS)-1 IF (NPTS .GT. MAXPTS) stop 'MONITR' LIROW = LLROW+NROWS+1 LICOL = LIROW+NROWS CALL SETXYZ (XL, YF, ZD, DX, DY, DZ, + ISTRUC(LLPLN), ISTRUC(LIPLN), + ISTRUC(LLROW), ISTRUC(LIROW), ISTRUC(LICOL), X, Y, Z) DX = DX/2 DY = DY/2 DZ = DZ/2 CALL PRERR (LEVEL, T, NPTS, NPDE, X, Y, Z, SOL(LSOL(LEVEL)+1), + UEX) 10 CONTINUE RETURN END SUBROUTINE PRERR (LEVEL, T, NPTS, NPDE, X, Y, Z, U, UEX) INTEGER LEVEL, NPTS, NPDE DOUBLE PRECISION T, X(NPTS), Y(NPTS), Z(NPTS), U(NPTS,NPDE), + UEX(NPTS,NPDE) INTEGER I,J DOUBLE PRECISION RMAX(3) CALL PDEIV (T, X, Y, Z, UEX, NPTS, NPDE) DO 1 J = 1, NPDE RMAX(J) = 0.0 DO 10 I = 1, NPTS RMAX(J) = MAX(RMAX(J),ABS(UEX(I,J)-U(I,J))) 10 CONTINUE 1 CONTINUE WRITE(28,'(''Error at T='',E9.3,'', level='',I1,'' :'', + I10,3E10.3)') T, LEVEL, NPTS, (RMAX(J), J=1, NPDE) RETURN END SUBROUTINE PDEIV (T, X, Y, Z, U, NPTS, NPDE) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER NPTS, NPDE DOUBLE PRECISION T, X(NPTS), Y(NPTS), Z(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,Y,Z : IN. Arrays of physical 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 INTEGER I DOUBLE PRECISION EPS PARAMETER (EPS = 5.0D-3) DO 10 I = 1, NPTS U(I,1) = 1-0.5/(1+EXP((-X(I)+Y(I)+Z(I)-0.75*T)/(4*EPS))) U(I,2) = 1.5-U(I,1) U(I,3) = 1.5-U(I,1) 10 CONTINUE RETURN END SUBROUTINE PDEF (T, X, Y, Z, U, + UT, UX, UY, UZ, UXX, UYY, UZZ, UXY, UXZ, UYZ, RES, NPTS, NPDE) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER NPTS, NPDE DOUBLE PRECISION T, X(NPTS), Y(NPTS), Z(NPTS), U(NPTS,NPDE), + UT(NPTS,NPDE), UX(NPTS,NPDE), UY(NPTS,NPDE), UZ(NPTS,NPDE), + UXX(NPTS,NPDE), UYY(NPTS,NPDE), UZZ(NPTS,NPDE), + UXY(NPTS,NPDE), UXZ(NPTS,NPDE), UYZ(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,Y,Z : IN. Physical coordinates of 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 UZ : IN. I C UXX : IN. I C UYY : IN. I Space derivatives of U on current grid C UZZ : IN. I C UXY : IN. I C UXZ : IN. I C UYZ : 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 INTEGER I DOUBLE PRECISION EPS PARAMETER (EPS = 5.0D-3) DO 10 I = 1, NPTS RES(I,1) = UT(I,1) + U(I,1)*UX(I,1) + + U(I,2)*UY(I,1) + U(I,3)*UZ(I,1) - + EPS*(UXX(I,1)+UYY(I,1)+UZZ(I,1)) RES(I,2) = UT(I,2) + U(I,1)*UX(I,2) + + U(I,2)*UY(I,2) + U(I,3)*UZ(I,2) - + EPS*(UXX(I,2)+UYY(I,2)+UZZ(I,2)) RES(I,3) = UT(I,3) + U(I,1)*UX(I,3) + + U(I,2)*UY(I,3) + U(I,3)*UZ(I,3) - + EPS*(UXX(I,3)+UYY(I,3)+UZZ(I,3)) 10 CONTINUE RETURN END SUBROUTINE PDEBC (T, X, Y, Z, U, UT, UX, UY, UZ, 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), Z(NPTS), U(NPTS,NPDE), + UT(NPTS,NPDE), UX(NPTS,NPDE), UY(NPTS,NPDE), UZ(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,Y,Z : IN. Physical coordinates of 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 UZ : IN. -I 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)+2) C LLBND(0) = NBNDS: total # physical planes in actual domain. C NB. edges and corners are stored for each plane they C belong to. C LLBND(1:NBNDS): pointers to a specific boundary in 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 : IN. (NBNDS) C ILBND(IB): type of boundary: C 1: Left plane -I C 2: Down plane I C 3: Right plane I max. first order derivative C 4: Up plane I C 5: Front plane I C 6: Back plane -I C LBND : IN. (NBIPTS) C LBND(IBPT): pointer to boundary point in actual grid C C----------------------------------------------------------------------- C INTEGER I, J, K, NBNDS DOUBLE PRECISION EPS, UI PARAMETER (EPS = 5.0D-3) NBNDS = LLBND(0) DO 10 J = 1, NBNDS DO 20 K = LLBND(J), LLBND(J+1)-1 I = LBND(K) UI = 1-0.5/(1+EXP((-X(I)+Y(I)+Z(I)-0.75*T)/(4*EPS))) RES(I,1) = U(I,1) - UI RES(I,2) = U(I,2) - (1.5-UI) RES(I,3) = U(I,3) - (1.5-UI) 20 CONTINUE 10 CONTINUE RETURN END SHAR_EOF fi # end of overwriting check if test -f 'probii.f' then echo shar: will not over-write existing file "'probii.f'" else cat << \SHAR_EOF > 'probii.f' PROGRAM EXII C INTEGER MXLEV, NPD, NPTS, LENIWK, LENRWK, LENLWK PARAMETER (MXLEV=2, NPD=2, NPTS=70000) C PARAMETER (LENIWK=NPTS*(7*MXLEV+43), C + LENRWK=NPTS*NPD*(5*MXLEV+38*NPD+13), PARAMETER (LENIWK=NPTS*(7*MXLEV+24), + LENRWK=NPTS*NPD*(5*MXLEV+37+21), + LENLWK=2*NPTS) C C----------------------------------------------------------------------- C INTEGER LUNDMP PARAMETER (LUNDMP = 89) C CHARACTER*2 NR INTEGER I, NPDE, INFO(7), IWK(LENIWK), MNTR LOGICAL LWK(LENLWK) DOUBLE PRECISION T, TOUT, TE, DT, XL, YF, ZD, XR, YB, ZU, DX, DY, + DZ, + TOLS, TOLT, RINFO(2+3*NPD), RWK(LENRWK) C First call of VLUGR3 MNTR = 0 NPDE = 2 T = 0.0 TOUT = 1.0 DT = 1D-5 XL = 0.0 YF = 0.0 ZD = 0.0 XR = 1.0 YB = 1.0 ZU = 1.0 DX = 0.1 DY = 0.1 DZ = 0.1 TOLS = 0.1 TOLT = 0.1 INFO(1) = 1 C MAXLEV INFO(2) = 4 C Domain a rectangular prism INFO(3) = 0 C Linear system solver print *, 'Lin.sys.solver BiCGStab or GCRO (0 / 10,11,12,13) ?' 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 solver info to unit # 61 INFO(7) = 61 C DTMIN = 1D-9 RINFO(1) = 1.0D-9 C DTMAX = TE - TS RINFO(2) = 0.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 DO 10 I = 1, 10 TE = T + TOUT/10 CALL VLUGR3 (NPDE, T, TE, DT, XL,YF,ZD, XR,YB,ZU, DX, DY, DZ, + TOLS, TOLT, INFO, RINFO, RWK, LENRWK, IWK, LENIWK, LWK, LENLWK, + MNTR) PRINT *, 'VLUGR3 returned with MNTR=', MNTR C C Save info on file WRITE(NR,'(I2.2)') I OPEN(UNIT=LUNDMP,FILE='DUMP'//NR,FORM='UNFORMATTED') CALL DUMP (LUNDMP, RWK, IWK) CLOSE(LUNDMP) 10 CONTINUE END SUBROUTINE MONITR (T, DT, DTNEW, XL, YF, ZD, DXB, DYB, DZB, + LGRID, ISTRUC, LSOL, SOL) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER LGRID(0:*), ISTRUC(*), LSOL(*) DOUBLE PRECISION T, DT, DTNEW, XL, YF, ZD, DXB, DYB, DZB, 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 left/front/down point of (virtual) box C YF : IN. Y-coordinate of left/front/down point of (virtual) box C ZD : IN. Z-coordinate of left/front/down point of (virtual) box C DXB : IN. Cell width in X-direction of base grid C DYB : IN. Cell width in Y-direction of base grid C DZB : IN. Cell width in Z-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 C (LPLN,IPLN,LROW,IROW,ICOL) of refinement C level LEVEL for time T C ISTRUC : IN. (*) C ISTRUC(LGRID(LEVEL):.) contains (LPLN,IPLN,LROW,IROW,ICOL) C of grid level LEVEL, C LPLN : (0:LPLN(0)+1) C LPLN(0) = NPLNS: Actual # planes in LROW C LPLN(1:NPLNS): pointers to the start of a plane in LROW C LPLN(NPLNS+1) = NROWS+1: Total # rows in grid + 1 C IPLN : (NPLNS) C IPLN(IP): plane number of plane IP in virtual box C LROW : (NROWS+1) 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 box C ICOL : (NPTS) C ICOL(IPT): column number of grid point IPT in virtual box 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=50000, NPDE=2) DOUBLE PRECISION X(MAXPTS), Y(MAXPTS), Z(MAXPTS), UEX(MAXPTS*NPDE) C C----------------------------------------------------------------------- C INTEGER MAXLEV, LEVEL, LLPLN, LIPLN, LLROW, LIROW, LICOL, + NPLNS, NROWS, NPTS DOUBLE PRECISION DX, DY, DZ C C Loop over the grid levels from coarse to fine. C Get physical coordinates of grid points C Compute ||err||_max and clip negative solution values MAXLEV = LGRID(0) DX = DXB DY = DYB DZ = DZB DO 10 LEVEL = 1, MAXLEV LLPLN = LGRID(LEVEL) NPLNS = ISTRUC(LLPLN) NROWS = ISTRUC(LLPLN+NPLNS+1)-1 LIPLN = LLPLN+NPLNS+2 LLROW = LIPLN+NPLNS NPTS = ISTRUC(LLROW+NROWS)-1 IF (NPTS .GT. MAXPTS) stop 'MONITR' LIROW = LLROW+NROWS+1 LICOL = LIROW+NROWS CALL SETXYZ (XL, YF, ZD, DX, DY, DZ, + ISTRUC(LLPLN), ISTRUC(LIPLN), + ISTRUC(LLROW), ISTRUC(LIROW), ISTRUC(LICOL), X, Y, Z) DX = DX/2 DY = DY/2 DZ = DZ/2 CALL PRERR (LEVEL, T, NPTS, NPDE, X, Y, Z, SOL(LSOL(LEVEL)+1), + UEX) 10 CONTINUE RETURN END SUBROUTINE PRERR (LEVEL, T, NPTS, NPDE, X, Y, Z, U, UEX) INTEGER LEVEL, NPTS, NPDE DOUBLE PRECISION T, X(NPTS), Y(NPTS), Z(NPTS), U(NPTS,NPDE), + UEX(NPTS,NPDE) INTEGER I DOUBLE PRECISION RMAX1, RMAX2 CALL PDEIV (T, X, Y, Z, UEX, NPTS, NPDE) RMAX1 = 0.0 RMAX2 = 0.0 DO 10 I = 1, NPTS RMAX1 = MAX(RMAX1,ABS(UEX(I,1)-U(I,1))) RMAX2 = MAX(RMAX2,ABS(UEX(I,2)-U(I,2))) IF (U(I,1) .LT. 0) U(I,1) = 0.0 IF (U(I,2) .LT. 0) U(I,2) = 0.0 10 CONTINUE WRITE(28,'("Error at T=",E9.3,", level=",I1," :",2E10.3,I10)') + T, LEVEL, RMAX1, RMAX2, NPTS RETURN END SUBROUTINE PDEIV (T, X, Y, Z, U, NPTS, NPDE) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER NPTS, NPDE DOUBLE PRECISION T, X(NPTS), Y(NPTS), Z(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,Y,Z : IN. Arrays of physical 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 DOUBLE PRECISION PI, S2, R, S, D, K1, K2, C2, EDKT PARAMETER (PI = 3.141592653589793, S2 = 1.414213562373095) PARAMETER (K1 = 1000.0, K2 = 1.0) PARAMETER (C2 = -K1/K2) INTEGER I R = (2+SIN(2*PI*T))/4 S = (2+S2/2*COS(2*PI*T))/4 DO 10 I = 1, NPTS D = EXP(-80*((X(I)-R)**2+(Y(I)-S)**2+(Z(I)-S)**2)) EDKT = EXP(-D*K2*T) U(I,1) = D/(K1+K2) * (K1*(1-C2)+(K1+K2)*C2*EDKT)/(1-C2+C2*EDKT) U(I,2) = D - U(I,1) 10 CONTINUE RETURN END SUBROUTINE PDEF (T, X, Y, Z, U, + UT, UX, UY, UZ, UXX, UYY, UZZ, UXY, UXZ, UYZ, RES, NPTS, NPDE) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER NPTS, NPDE DOUBLE PRECISION T, X(NPTS), Y(NPTS), Z(NPTS), U(NPTS,NPDE), + UT(NPTS,NPDE), UX(NPTS,NPDE), UY(NPTS,NPDE), UZ(NPTS,NPDE), + UXX(NPTS,NPDE), UYY(NPTS,NPDE), UZZ(NPTS,NPDE), + UXY(NPTS,NPDE), UXZ(NPTS,NPDE), UYZ(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,Y,Z : IN. Physical coordinates of 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 UZ : IN. I C UXX : IN. I C UYY : IN. I Space derivatives of U on current grid C UZZ : IN. I C UXY : IN. I C UXZ : IN. I C UYZ : 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 DOUBLE PRECISION PI, S2, K1, K2 PARAMETER (PI = 3.141592653589793, S2 = 1.414213562373095) PARAMETER (K1 = 1000.0, K2 = 1.0) INTEGER I DO 10 I = 1, NPTS RES(I,1) = UT(I,1) + + 2*PI*S2*((Y(I)+Z(I))/2-0.5)*UX(I,1) - + 2*PI/S2*(X(I)-0.5)*(UY(I,1)+UZ(I,1)) + - (-K2*U(I,1)*U(I,2) + K1*U(I,2)*U(I,2)) RES(I,2) = UT(I,2) + + 2*PI*S2*((Y(I)+Z(I))/2-0.5)*UX(I,2) - + 2*PI/S2*(X(I)-0.5)*(UY(I,2)+UZ(I,2)) + - (-K1*U(I,2)*U(I,2) + K2*U(I,1)*U(I,2)) 10 CONTINUE RETURN END SUBROUTINE PDEBC (T, X, Y, Z, U, UT, UX, UY, UZ, 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), Z(NPTS), U(NPTS,NPDE), + UT(NPTS,NPDE), UX(NPTS,NPDE), UY(NPTS,NPDE), UZ(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,Y,Z : IN. Physical coordinates of 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 UZ : IN. -I 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)+2) C LLBND(0) = NBNDS: total # physical planes in actual domain. C NB. edges and corners are stored for each plane they C belong to. C LLBND(1:NBNDS): pointers to a specific boundary in 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 : IN. (NBNDS) C ILBND(IB): type of boundary: C 1: Left plane -I C 2: Down plane I C 3: Right plane I max. first order derivative C 4: Up plane I C 5: Front plane I C 6: Back plane -I C LBND : IN. (NBIPTS) C LBND(IBPT): pointer to boundary point in actual grid C C----------------------------------------------------------------------- C DOUBLE PRECISION PI, S2, R, S, D, K1, K2, C2, EDKT, UI1 PARAMETER (PI = 3.141592653589793, S2 = 1.414213562373095) PARAMETER (K1 = 1000.0, K2 = 1.0) PARAMETER (C2 = -K1/K2) INTEGER I, J, K, NBNDS R = (2+SIN(2*PI*T))/4 S = (2+S2/2*COS(2*PI*T))/4 C NBNDS = LLBND(0) DO 10 J = 1, NBNDS C C Change the inflow boundaries IF (ILBND(J) .EQ. 1) THEN DO 20 K = LLBND(J), LLBND(J+1)-1 I = LBND(K) IF (Y(I)+Z(I) .GE. 1.0) THEN D = EXP(-80*((X(I)-R)**2+(Y(I)-S)**2+(Z(I)-S)**2)) EDKT = EXP(-D*K2*T) UI1 = D/(K1+K2) * + (K1*(1-C2)+(K1+K2)*C2*EDKT)/(1-C2+C2*EDKT) RES(I,1) = U(I,1) - UI1 RES(I,2) = U(I,2) - (D - UI1) ENDIF 20 CONTINUE ELSE IF (ILBND(J) .EQ. 2 .OR. ILBND(J) .EQ. 5) THEN DO 30 K = LLBND(J), LLBND(J+1)-1 I = LBND(K) IF (X(I) .LE. 0.5) THEN D = EXP(-80*((X(I)-R)**2+(Y(I)-S)**2+(Z(I)-S)**2)) EDKT = EXP(-D*K2*T) UI1 = D/(K1+K2) * + (K1*(1-C2)+(K1+K2)*C2*EDKT)/(1-C2+C2*EDKT) RES(I,1) = U(I,1) - UI1 RES(I,2) = U(I,2) - (D - UI1) ENDIF 30 CONTINUE ELSE IF (ILBND(J) .EQ. 3)THEN DO 40 K = LLBND(J), LLBND(J+1)-1 I = LBND(K) IF (Y(I)+Z(I) .LE. 1.0) THEN D = EXP(-80*((X(I)-R)**2+(Y(I)-S)**2+(Z(I)-S)**2)) EDKT = EXP(-D*K2*T) UI1 = D/(K1+K2) * + (K1*(1-C2)+(K1+K2)*C2*EDKT)/(1-C2+C2*EDKT) RES(I,1) = U(I,1) - UI1 RES(I,2) = U(I,2) - (D - UI1) ENDIF 40 CONTINUE ELSE IF (ILBND(J) .EQ. 4 .OR. ILBND(J) .EQ. 6) THEN DO 50 K = LLBND(J), LLBND(J+1)-1 I = LBND(K) IF (X(I) .GE. 0.5) THEN D = EXP(-80*((X(I)-R)**2+(Y(I)-S)**2+(Z(I)-S)**2)) EDKT = EXP(-D*K2*T) UI1 = D/(K1+K2) * + (K1*(1-C2)+(K1+K2)*C2*EDKT)/(1-C2+C2*EDKT) RES(I,1) = U(I,1) - UI1 RES(I,2) = U(I,2) - (D - UI1) ENDIF 50 CONTINUE ENDIF 10 CONTINUE RETURN END SHAR_EOF fi # end of overwriting check 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,YFW,ZDW, XRW,YBW,ZUW, DXB,DYB, + DZB, DTO COMMON /WRITIF/ MAXLVW, NPDEW, LRWKPS, LIWKPS, LRWKB, LIWKB COMMON /WRITLF/ FIRST, SECOND COMMON /WRITRF/ T0,TW,TEW,DTW, XLW,YFW,ZDW, XRW,YBW,ZUW, + DXB,DYB,DZB, DTO SAVE /WRITIF/, /WRITLF/, /WRITRF/ C C end INCLUDE 'CMNWRITEF' C C C----------------------------------------------------------------------- C INTEGER MXLEV, NPD, NPTS, LENIWK, LENRWK PARAMETER (MXLEV=3, NPD=1, NPTS=100000) PARAMETER (LENIWK=NPTS*(7*MXLEV+23), + 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, YFW, ZDW, DXB, DYB, DZB, + 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 VLUGR3 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 VLUGR3 C LENRWK : IN. Dimension of RWK. C IWK : OUT. Integer workstorage intended to pass to VLUGR3 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,YFW,ZDW, XRW,YBW,ZUW, DXB,DYB, + DZB, DTO COMMON /WRITIF/ MAXLVW, NPDEW, LRWKPS, LIWKPS, LRWKB, LIWKB COMMON /WRITLF/ FIRST, SECOND COMMON /WRITRF/ T0,TW,TEW,DTW, XLW,YFW,ZDW, XRW,YBW,ZUW, + DXB,DYB,DZB, 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,YFW,ZDW, XRW,YBW,ZUW, DXB,DYB,DZB, 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, YF, ZD, DXB, DYB, DZB, + LGRID, ISTRUC, LSOL, SOL) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER LUN, NPDE, LGRID(0:*), ISTRUC(*), LSOL(*) DOUBLE PRECISION T, XL, YF, ZD, DXB, DYB, DZB, 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 left/front/down point of virtual box C YF : IN. Y-coordinate of left/front/down point of virtual box C ZD : IN. Z-coordinate of left/front/down point of virtual box C DXB : IN. Cell width in X-direction of base grid C DYB : IN. Cell width in Y-direction of base grid C DZB : IN. Cell width in Z-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 C (LPLN,IPLN,LROW,IROW,ICOL) C of refinement level LEVEL for time T C ISTRUC : IN. (*) C ISTRUC(LGRID(LEVEL):.) contains (LPLN,IPLN,LROW,IROW,ICOL) C of grid level LEVEL, C LPLN : (0:LPLN(0)+1) C LPLN(0) = NPLNS: Actual # planes in grid C LPLN(1:NPLNS): pointers to the start of a plane in LROW C LPLN(NPLNS+1) = NROWS+1: Actual # rows in grid + 1 C IPLN : (NPLNS) C IPLN(IP): plane number of plane IP in virtual box C LROW : (NROWS+1) 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 box C ICOL : (NPTS) C ICOL(IPT): column number of grid point IPT in virtual box 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, LLPLN, LIPLN, LLROW, LIROW, LICOL, + NPLNS, NROWS, NPTS DOUBLE PRECISION DX, DY, DZ MAXLEV = LGRID(0) DX = DXB DY = DYB DZ = DZB DO 10 LEVEL = 1, MAXLEV LLPLN = LGRID(LEVEL) NPLNS = ISTRUC(LLPLN) NROWS = ISTRUC(LLPLN+NPLNS+1)-1 LIPLN = LLPLN+NPLNS+2 LLROW = LIPLN+NPLNS NPTS = ISTRUC(LLROW+NROWS)-1 LIROW = LLROW+NROWS+1 LICOL = LIROW+NROWS CALL PRSOLL (LUN, LEVEL, T, NPTS, NPDE, XL, YF, ZD, DX, DY, DZ, + ISTRUC(LLPLN), ISTRUC(LIPLN), ISTRUC(LLROW), ISTRUC(LIROW), + ISTRUC(LICOL), SOL(LSOL(LEVEL)+1)) DX = DX/2 DY = DY/2 DZ = DZ/2 10 CONTINUE RETURN END SUBROUTINE PRSOLL (LUN, LEVEL, T, NPTS, NPDE, XL, YF, ZD, + DX, DY, DZ, LPLN, IPLN, LROW, IROW, ICOL, U) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER LUN, LEVEL, NPTS, NPDE, LPLN(0:*), IPLN(*), + LROW(*), IROW(*), ICOL(*) DOUBLE PRECISION T, XL, YF, ZD, DX, DY, DZ, U(NPTS,NPDE) C Ccc PURPOSE: C Print solution and X-, Y- and Z-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 left/front/down point of virtual box C YF : IN. Y-coordinate of left/front/down point of virtual box C ZD : IN. Z-coordinate of left/front/down point of virtual box C DX : IN. Grid width in X-direction C DY : IN. Grid width in Y-direction C DZ : IN. Grid width in Z-direction C LPLN : IN. (0:LPLN(0)+1) C LPLN(0) = NPLNS: Actual # planes in LROW C LPLN(1:NPLNS): pointers to the start of a plane in LROW C LPLN(NPLNS+1) = NROWS+1: Total # rows in grid + 1 C IPLN : IN. (NPLNS) C IPLN(IP): plane number of plane IP in virtual box C LROW : IN. (NROWS+1) 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 box C ICOL : IN. (NPTS) C ICOL(IPT): column number of grid point IPT in virtual box C U : IN. Solution on this grid level C Ccc EXTERNALS USED: NONE C C----------------------------------------------------------------------- C INTEGER IC, IP, IPT, IR, NPLNS DOUBLE PRECISION X, Y, Z C NPLNS = LPLN(0) WRITE(LUN,'(//// A,T14,A,T30,A,T46,A,T62,A,T71,A //)') + 'Lev', 't', 'Z', 'Y', 'X', 'Solution' IP = 1 Z = ZD + IPLN(IP)*DZ IR = LPLN(IP) Y = YF + IROW(IR)*DY IPT = LROW(IR) X = XL + ICOL(IPT)*DX WRITE(LUN, + '(I3,T5,E12.5,T21,E12.5,T37,E12.5,T53,E12.5,T69,E12.5)') + LEVEL, T, Z, Y, X, U(IPT,1) DO 10 IC = 2, NPDE WRITE(LUN,'(T69,E12.5)') U(IPT,IC) 10 CONTINUE DO 14 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 15 IC = 2, NPDE WRITE(LUN,'(T69,E12.5)') U(IPT,IC) 15 CONTINUE 14 CONTINUE DO 20 IR = LPLN(IP)+1, LPLN(IP+1)-1 Y = YF + IROW(IR)*DY IPT = LROW(IR) X = XL + ICOL(IPT)*DX WRITE(LUN, + '(T37,E12.5,T53,E12.5,T69,E12.5)') + Y, X, U(IPT,1) DO 30 IC = 2, NPDE WRITE(LUN,'(T69,E12.5)') U(IPT,IC) 30 CONTINUE DO 40 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 50 IC = 2, NPDE WRITE(LUN,'(T69,E12.5)') U(IPT,IC) 50 CONTINUE 40 CONTINUE 20 CONTINUE DO 100 IP = 2, NPLNS Z = ZD + IPLN(IP)*DZ IR = LPLN(IP) Y = YF + 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)') + Z, Y, X, U(IPT,1) DO 110 IC = 2, NPDE WRITE(LUN,'(T69,E12.5)') U(IPT,IC) 110 CONTINUE DO 114 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 115 IC = 2, NPDE WRITE(LUN,'(T69,E12.5)') U(IPT,IC) 115 CONTINUE 114 CONTINUE DO 120 IR = LPLN(IP)+1, LPLN(IP+1)-1 Y = YF + IROW(IR)*DY IPT = LROW(IR) X = XL + ICOL(IPT)*DX WRITE(LUN, + '(T37,E12.5,T53,E12.5,T69,E12.5)') + Y, X, U(IPT,1) DO 130 IC = 2, NPDE WRITE(LUN,'(T69,E12.5)') U(IPT,IC) 130 CONTINUE DO 140 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 150 IC = 2, NPDE WRITE(LUN,'(T69,E12.5)') U(IPT,IC) 150 CONTINUE 140 CONTINUE 120 CONTINUE 100 CONTINUE RETURN END 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. 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,YFW,ZDW, XRW,YBW,ZUW, DXB,DYB, + DZB, DTO COMMON /WRITIF/ MAXLVW, NPDEW, LRWKPS, LIWKPS, LRWKB, LIWKB COMMON /WRITLF/ FIRST, SECOND COMMON /WRITRF/ T0,TW,TEW,DTW, XLW,YFW,ZDW, XRW,YBW,ZUW, + DXB,DYB,DZB, DTO SAVE /WRITIF/, /WRITLF/, /WRITRF/ C C end INCLUDE 'CMNWRITEF' C C C----------------------------------------------------------------------- C INTEGER MXLEV, NPD, NPTS, LENIWK, LENRWK PARAMETER (MXLEV=3, NPD=1, NPTS=100000) PARAMETER (LENIWK=NPTS*(7*MXLEV+23), + LENRWK=5*NPTS*NPD*MXLEV) C CHARACTER FILE*128 INTEGER IWK(LENIWK), + LSGNM1, LSGN, LSGNP1, LSUNM1, LSSN, LSUN, + LUNI, MAXLEV, NX, NXB, NY, NYB, NZ, NZB, 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((YBW - YFW)/DYB) NZB = NINT((ZUW - ZDW)/DZB) NX = NXB * 2**(UNILEV-1) NY = NYB * 2**(UNILEV-1) NZ = NZB * 2**(UNILEV-1) LUNI = LENRWK - (NX+1)*(NY+1)*(NZ+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, YF, ZD, DXB, DYB, DZB, NXB, NYB, NZB' PRINT *, TW, NPDEW, XLW, YFW, ZDW, DXB, DYB, DZB, NXB, NYB, NZB FILE = 'sol.dat' OPEN(UNIT=61,FILE=FILE) FILE = 'grid.dat' OPEN(UNIT=63,FILE=FILE) CALL WRUNI (61, 63, UNILEV, + TW, NPDEW, XLW, YFW, ZDW, DXB, DYB, DZB, NXB, NYB, NZB, + IWK(LSGN), IWK(LIWKPS), IWK(LSUN), RWK(LRWKPS), + RWK(LUNI), NX, NY, NZ) 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 VLUGR3 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 VLUGR3 C LENRWK : IN. Dimension of RWK. C IWK : OUT. Integer workstorage intended to pass to VLUGR3 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,YFW,ZDW, XRW,YBW,ZUW, DXB,DYB, + DZB, DTO COMMON /WRITIF/ MAXLVW, NPDEW, LRWKPS, LIWKPS, LRWKB, LIWKB COMMON /WRITLF/ FIRST, SECOND COMMON /WRITRF/ T0,TW,TEW,DTW, XLW,YFW,ZDW, XRW,YBW,ZUW, + DXB,DYB,DZB, 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,YFW,ZDW, XRW,YBW,ZUW, DXB,DYB,DZB, 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, YF, ZD, DXB, DYB, DZB, NXB, NYB, NZB, + LGRID, ISTRUC, LSOL, SOL, UNIFRM, NX, NY, NZ) C C----------------------------------------------------------------------- C Ccc PARAMETER SPECIFICATION: INTEGER LUNS, LUNG, UNILEV, + NPDE, NXB, NYB, NZB, LGRID(0:*), ISTRUC(*), LSOL(*), NX, NY, NZ DOUBLE PRECISION T, XL, YF, ZD, DXB, DYB, DZB, SOL(*), + UNIFRM(0:NX,0:NY,0:NZ,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. 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 left/front/down point of virtual box C YF : IN. Y-coordinate of left/front/down point of virtual box C ZD : IN. Z-coordinate of left/front/down point of virtual box C DXB : IN. Cell width in X-direction of base grid C DYB : IN. Cell width in Y-direction of base grid C DZB : IN. Cell width in Z-direction of base grid C NXB,NYB,NZB: IN. # gridcells in X-, Y- and Z-direction, resp., on grid C of base level 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 C (LPLN,IPLN,LROW,IROW,ICOL) C of refinement level LEVEL for time T C ISTRUC : IN. (*) C ISTRUC(LGRID(LEVEL):.) contains (LPLN,IPLN,LROW,IROW,ICOL) C of grid level LEVEL, C LPLN : (0:LPLN(0)+1) C LPLN(0) = NPLNS: Actual # planes in grid C LPLN(1:NPLNS): pointers to the start of a plane in LROW C LPLN(NPLNS+1) = NROWS+1: Actual # rows in grid + 1 C IPLN : (NPLNS) C IPLN(IP): plane number of plane IP in virtual box C LROW : (NROWS+1) 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 box C ICOL : (NPTS) C ICOL(IPT): column number of grid point IPT in virtual box 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,NZ: IN. # gridcells in X-, Y- and Z-direction, resp., on grid 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, IP, IPLN, IPT, IR, IROW, J, K, + LEVEL, LLPLN, LIPLN, LLROW, LIROW, LICOL, MAXLEV, + NPLNS, NROWS, NPTS DO 1 IC = 1, NPDE DO 1 IPLN = 0, NZ DO 1 IROW = 0, NY DO 1 ICOL = 0, NX UNIFRM(ICOL,IROW,IPLN,IC) = NONVAL 1 CONTINUE MAXLEV = LGRID(0) DO 10 LEVEL = 1, UNILEV IMUL = 2**(UNILEV-LEVEL) LLPLN = LGRID(LEVEL) NPLNS = ISTRUC(LLPLN) NROWS = ISTRUC(LLPLN+NPLNS+1)-1 LIPLN = LLPLN+NPLNS+2 LLROW = LIPLN+NPLNS NPTS = ISTRUC(LLROW+NROWS)-1 LIROW = LLROW+NROWS+1 LICOL = LIROW+NROWS DO 20 IP = 1, NPLNS IPLN = ISTRUC(LIPLN-1+IP)*IMUL DO 30 IR = ISTRUC(LLPLN+IP), ISTRUC(LLPLN+IP+1)-1 IROW = ISTRUC(LIROW-1+IR)*IMUL DO 40 IPT = ISTRUC(LLROW-1+IR), ISTRUC(LLROW+IR)-1 ICOL = ISTRUC(LICOL-1+IPT)*IMUL DO 50 IC = 1, NPDE UNIFRM(ICOL,IROW,IPLN,IC) = + SOL(LSOL(LEVEL)+(IC-1)*NPTS+IPT) 50 CONTINUE 40 CONTINUE 30 CONTINUE 20 CONTINUE 10 CONTINUE DO 100 LEVEL = 2, UNILEV IMUL = 2**(UNILEV-LEVEL) DO 110 K = IMUL, NZ, IMUL*2 DO 110 J = 0, NY, IMUL*2 DO 110 I = 0, NX, IMUL*2 IF (UNIFRM(I,J,K,1) .EQ. NONVAL) THEN DO 120 IC = 1, NPDE UNIFRM(I,J,K,IC) = + (UNIFRM(I,J,K-IMUL,IC)+UNIFRM(I,J,K+IMUL,IC))/2 120 CONTINUE ENDIF 110 CONTINUE DO 130 K = 0, NZ, IMUL DO 130 J = IMUL, NY, IMUL*2 DO 130 I = 0, NX, IMUL*2 IF (UNIFRM(I,J,K,1) .EQ. NONVAL) THEN DO 140 IC = 1, NPDE UNIFRM(I,J,K,IC) = + (UNIFRM(I,J-IMUL,K,IC)+UNIFRM(I,J+IMUL,K,IC))/2 140 CONTINUE ENDIF 130 CONTINUE DO 150 K = 0, NZ, IMUL DO 150 J = 0, NY, IMUL DO 150 I = IMUL, NX, IMUL*2 IF (UNIFRM(I,J,K,1) .EQ. NONVAL) THEN DO 160 IC = 1, NPDE UNIFRM(I,J,K,IC) = + (UNIFRM(I-IMUL,J,K,IC)+UNIFRM(I+IMUL,J,K,IC))/2 160 CONTINUE ENDIF 150 CONTINUE 100 CONTINUE DO 170 K = 0, NZ DO 170 J = 0, NY DO 170 I = 0, NX WRITE(LUNS,'(100E13.3)') (UNIFRM(I,J,K,IC), IC = 1, NPDE) 170 CONTINUE C C Grids DO 201 IPLN = 0, NZ DO 201 IROW = 0, NY DO 201 ICOL = 0, NX UNIFRM(ICOL,IROW,IPLN,1) = 0 201 CONTINUE DO 210 LEVEL = 1, UNILEV IMUL = 2**(UNILEV-LEVEL) LLPLN = LGRID(LEVEL) NPLNS = ISTRUC(LLPLN) NROWS = ISTRUC(LLPLN+NPLNS+1)-1 LIPLN = LLPLN+NPLNS+2 LLROW = LIPLN+NPLNS NPTS = ISTRUC(LLROW+NROWS)-1 LIROW = LLROW+NROWS+1 LICOL = LIROW+NROWS DO 220 IP = 1, NPLNS IPLN = ISTRUC(LIPLN-1+IP)*IMUL DO 230 IR = ISTRUC(LLPLN+IP), ISTRUC(LLPLN+IP+1)-1 IROW = ISTRUC(LIROW-1+IR)*IMUL DO 240 IPT = ISTRUC(LLROW-1+IR), ISTRUC(LLROW+IR)-1 ICOL = ISTRUC(LICOL-1+IPT)*IMUL UNIFRM(ICOL,IROW,IPLN,1) = LEVEL 240 CONTINUE 230 CONTINUE 220 CONTINUE 210 CONTINUE DO 300 LEVEL = 2, UNILEV IMUL = 2**(UNILEV-LEVEL) DO 310 K = IMUL, NZ, IMUL*2 DO 310 J = 0, NY, IMUL*2 DO 310 I = 0, NX, IMUL*2 IF (UNIFRM(I,J,K,1) .LT. LEVEL) THEN UNIFRM(I,J,K,1) = + MIN(UNIFRM(I,J,K-IMUL,1),UNIFRM(I,J,K+IMUL,1)) ENDIF 310 CONTINUE DO 320 K = 0, NZ, IMUL DO 320 J = IMUL, NY, IMUL*2 DO 320 I = 0, NX, IMUL*2 IF (UNIFRM(I,J,K,1) .LT. LEVEL) THEN UNIFRM(I,J,K,1) = + MIN(UNIFRM(I,J-IMUL,K,1),UNIFRM(I,J+IMUL,K,1)) ENDIF 320 CONTINUE DO 330 K = 0, NZ, IMUL DO 330 J = 0, NY, IMUL DO 330 I = IMUL, NX, IMUL*2 IF (UNIFRM(I,J,K,1) .LT. LEVEL) THEN UNIFRM(I,J,K,1) = + MIN(UNIFRM(I-IMUL,J,K,1),UNIFRM(I+IMUL,J,K,1)) ENDIF 330 CONTINUE 300 CONTINUE DO 350 K = 0, NZ DO 350 J = 0, NY DO 350 I = 0, NX WRITE(LUNG,'(I2)') NINT(UNIFRM(I,J,K,1)) 350 CONTINUE RETURN END SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Sp' then mkdir 'Sp' fi cd 'Sp' if test -f 'exmpl_28' then echo shar: will not over-write existing file "'exmpl_28'" else cat << \SHAR_EOF > 'exmpl_28' Error at T=0.100E-02, level=1 : 384 0.361E-03 0.361E-03 0.361E-03 Error at T=0.100E-02, level=2 : 2126 0.361E-03 0.361E-03 0.361E-03 Error at T=0.100E-02, level=3 : 9189 0.361E-03 0.361E-03 0.361E-03 Error at T=0.100E-02, level=4 : 43645 0.362E-03 0.362E-03 0.362E-03 Error at T=0.300E-02, level=1 : 384 0.102E-02 0.102E-02 0.102E-02 Error at T=0.300E-02, level=2 : 2126 0.102E-02 0.102E-02 0.102E-02 Error at T=0.300E-02, level=3 : 9697 0.103E-02 0.103E-02 0.103E-02 Error at T=0.300E-02, level=4 : 43645 0.104E-02 0.104E-02 0.104E-02 Error at T=0.699E-02, level=1 : 384 0.206E-02 0.206E-02 0.206E-02 Error at T=0.699E-02, level=2 : 2126 0.206E-02 0.206E-02 0.206E-02 Error at T=0.699E-02, level=3 : 9697 0.209E-02 0.209E-02 0.209E-02 Error at T=0.699E-02, level=4 : 43645 0.215E-02 0.215E-02 0.215E-02 Error at T=0.149E-01, level=1 : 384 0.293E-02 0.293E-02 0.293E-02 Error at T=0.149E-01, level=2 : 2126 0.293E-02 0.293E-02 0.293E-02 Error at T=0.149E-01, level=3 : 10139 0.314E-02 0.314E-02 0.314E-02 Error at T=0.149E-01, level=4 : 46261 0.478E-02 0.478E-02 0.478E-02 Error at T=0.308E-01, level=1 : 384 0.692E-03 0.692E-03 0.692E-03 Error at T=0.308E-01, level=2 : 2126 0.857E-03 0.857E-03 0.857E-03 Error at T=0.308E-01, level=3 : 10139 0.595E-02 0.595E-02 0.595E-02 Error at T=0.308E-01, level=4 : 45947 0.959E-02 0.959E-02 0.959E-02 Error at T=0.602E-01, level=1 : 384 0.882E-02 0.883E-02 0.883E-02 Error at T=0.602E-01, level=2 : 2126 0.884E-02 0.884E-02 0.884E-02 Error at T=0.602E-01, level=3 : 10197 0.182E-01 0.182E-01 0.182E-01 Error at T=0.602E-01, level=4 : 45921 0.192E-01 0.192E-01 0.192E-01 Error at T=0.896E-01, level=1 : 384 0.141E-01 0.141E-01 0.141E-01 Error at T=0.896E-01, level=2 : 2126 0.141E-01 0.141E-01 0.141E-01 Error at T=0.896E-01, level=3 : 10197 0.161E-01 0.161E-01 0.161E-01 Error at T=0.896E-01, level=4 : 48125 0.279E-01 0.279E-01 0.279E-01 Error at T=0.119E+00, level=1 : 384 0.117E-01 0.117E-01 0.117E-01 Error at T=0.119E+00, level=2 : 2198 0.302E-01 0.302E-01 0.302E-01 Error at T=0.119E+00, level=3 : 10825 0.302E-01 0.302E-01 0.302E-01 Error at T=0.119E+00, level=4 : 52205 0.339E-01 0.339E-01 0.339E-01 Error at T=0.149E+00, level=1 : 384 0.659E-02 0.659E-02 0.659E-02 Error at T=0.149E+00, level=2 : 2318 0.215E-01 0.215E-01 0.215E-01 Error at T=0.149E+00, level=3 : 11215 0.215E-01 0.215E-01 0.215E-01 Error at T=0.149E+00, level=4 : 49905 0.387E-01 0.387E-01 0.387E-01 Error at T=0.179E+00, level=1 : 384 0.278E-02 0.278E-02 0.278E-02 Error at T=0.179E+00, level=2 : 2318 0.119E-01 0.119E-01 0.119E-01 Error at T=0.179E+00, level=3 : 11047 0.362E-01 0.362E-01 0.362E-01 Error at T=0.179E+00, level=4 : 54293 0.416E-01 0.416E-01 0.416E-01 Error at T=0.208E+00, level=1 : 384 0.239E-01 0.239E-01 0.239E-01 Error at T=0.208E+00, level=2 : 2318 0.239E-01 0.239E-01 0.239E-01 Error at T=0.208E+00, level=3 : 11303 0.239E-01 0.239E-01 0.239E-01 Error at T=0.208E+00, level=4 : 51793 0.435E-01 0.435E-01 0.435E-01 Error at T=0.236E+00, level=1 : 384 0.384E-01 0.384E-01 0.384E-01 Error at T=0.236E+00, level=2 : 2318 0.384E-01 0.384E-01 0.384E-01 Error at T=0.236E+00, level=3 : 11635 0.388E-01 0.388E-01 0.388E-01 Error at T=0.236E+00, level=4 : 56249 0.446E-01 0.446E-01 0.446E-01 Error at T=0.266E+00, level=1 : 384 0.198E-01 0.198E-01 0.198E-01 Error at T=0.266E+00, level=2 : 2318 0.210E-01 0.210E-01 0.210E-01 Error at T=0.266E+00, level=3 : 11997 0.271E-01 0.271E-01 0.271E-01 Error at T=0.266E+00, level=4 : 53359 0.454E-01 0.454E-01 0.454E-01 Error at T=0.294E+00, level=1 : 384 0.143E-01 0.143E-01 0.143E-01 Error at T=0.294E+00, level=2 : 2310 0.172E-01 0.172E-01 0.172E-01 Error at T=0.294E+00, level=3 : 11789 0.406E-01 0.406E-01 0.406E-01 Error at T=0.294E+00, level=4 : 57831 0.458E-01 0.458E-01 0.458E-01 Error at T=0.323E+00, level=1 : 384 0.233E-01 0.233E-01 0.233E-01 Error at T=0.323E+00, level=2 : 2318 0.297E-01 0.297E-01 0.297E-01 Error at T=0.323E+00, level=3 : 12215 0.302E-01 0.302E-01 0.302E-01 Error at T=0.323E+00, level=4 : 54777 0.459E-01 0.459E-01 0.459E-01 Error at T=0.351E+00, level=1 : 384 0.170E-01 0.170E-01 0.170E-01 Error at T=0.351E+00, level=2 : 2334 0.396E-01 0.396E-01 0.396E-01 Error at T=0.351E+00, level=3 : 12515 0.414E-01 0.414E-01 0.414E-01 Error at T=0.351E+00, level=4 : 59145 0.459E-01 0.459E-01 0.459E-01 Error at T=0.381E+00, level=1 : 384 0.794E-02 0.794E-02 0.794E-02 Error at T=0.381E+00, level=2 : 2358 0.175E-01 0.175E-01 0.175E-01 Error at T=0.381E+00, level=3 : 12577 0.328E-01 0.328E-01 0.328E-01 Error at T=0.381E+00, level=4 : 56612 0.456E-01 0.456E-01 0.456E-01 Error at T=0.409E+00, level=1 : 384 0.897E-02 0.897E-02 0.897E-02 Error at T=0.409E+00, level=2 : 2358 0.209E-01 0.209E-01 0.209E-01 Error at T=0.409E+00, level=3 : 12101 0.416E-01 0.416E-01 0.416E-01 Error at T=0.409E+00, level=4 : 60640 0.453E-01 0.453E-01 0.453E-01 Error at T=0.439E+00, level=1 : 384 0.324E-01 0.324E-01 0.324E-01 Error at T=0.439E+00, level=2 : 2358 0.345E-01 0.345E-01 0.345E-01 Error at T=0.439E+00, level=3 : 12717 0.352E-01 0.352E-01 0.352E-01 Error at T=0.439E+00, level=4 : 58192 0.448E-01 0.448E-01 0.448E-01 Error at T=0.467E+00, level=1 : 384 0.386E-01 0.386E-01 0.386E-01 Error at T=0.467E+00, level=2 : 2358 0.390E-01 0.390E-01 0.390E-01 Error at T=0.467E+00, level=3 : 12953 0.413E-01 0.413E-01 0.413E-01 Error at T=0.467E+00, level=4 : 61382 0.457E-01 0.457E-01 0.457E-01 Error at T=0.496E+00, level=1 : 384 0.120E-01 0.120E-01 0.120E-01 Error at T=0.496E+00, level=2 : 2358 0.137E-01 0.137E-01 0.137E-01 Error at T=0.496E+00, level=3 : 12947 0.372E-01 0.372E-01 0.372E-01 Error at T=0.496E+00, level=4 : 58684 0.459E-01 0.459E-01 0.459E-01 Error at T=0.524E+00, level=1 : 384 0.215E-01 0.215E-01 0.215E-01 Error at T=0.524E+00, level=2 : 2342 0.229E-01 0.229E-01 0.229E-01 Error at T=0.524E+00, level=3 : 12607 0.407E-01 0.407E-01 0.407E-01 Error at T=0.524E+00, level=4 : 61690 0.464E-01 0.464E-01 0.464E-01 Error at T=0.554E+00, level=1 : 384 0.235E-01 0.235E-01 0.235E-01 Error at T=0.554E+00, level=2 : 2342 0.383E-01 0.383E-01 0.383E-01 Error at T=0.554E+00, level=3 : 12967 0.390E-01 0.390E-01 0.390E-01 Error at T=0.554E+00, level=4 : 58832 0.477E-01 0.477E-01 0.477E-01 Error at T=0.582E+00, level=1 : 384 0.145E-01 0.145E-01 0.145E-01 Error at T=0.582E+00, level=2 : 2374 0.369E-01 0.369E-01 0.369E-01 Error at T=0.582E+00, level=3 : 13185 0.396E-01 0.396E-01 0.396E-01 Error at T=0.582E+00, level=4 : 61722 0.479E-01 0.479E-01 0.479E-01 Error at T=0.612E+00, level=1 : 384 0.589E-02 0.589E-02 0.589E-02 Error at T=0.612E+00, level=2 : 2374 0.924E-02 0.924E-02 0.924E-02 Error at T=0.612E+00, level=3 : 13091 0.405E-01 0.405E-01 0.405E-01 Error at T=0.612E+00, level=4 : 58634 0.492E-01 0.492E-01 0.492E-01 Error at T=0.640E+00, level=1 : 384 0.159E-01 0.159E-01 0.159E-01 Error at T=0.640E+00, level=2 : 2374 0.232E-01 0.232E-01 0.232E-01 Error at T=0.640E+00, level=3 : 12699 0.379E-01