#********************************************************************# # # # PWH/860414 # # # # # # A LIBRARY OF EULER MULTIGRID ROUTINES # # # # A PROGRAMME TO STUDY THE NUMERICAL BEHAVIOUR # # # # OF A MULTIGRID TECHNIQUE FOR THE SOLUTION OF # # THE EULER EQUATIONS FOR COMPRESSIBLE FLUID FLOW. # # # # P.W. Hemker # # # # # # REFERENCES: # # # # -- Hemker, P.W. & Spekreijse, S.P., # # Multigrid solution of the steady Euler equations. In: # # Advances in Multi-Grid methods, (D.Braess, W.Hackbusch & # # U.Trottenberg eds) pp. 33-44, Proceedings of the # # conference held in Oberwolfach, Dec. 1984, Vieweg Publ. # # Comp., Braunschweig, 1985 # # # # -- Hemker, P.W., # # Multigrid and improved accuracy for the Euler equations. # # In: Discretization in differential equations and # # enclosures, (E.Adams, R.Ansorge, Ch.Groszmann and # # H.-G.Roos eds) pp. 123-136, Akademie-Verlag Berlin, 1987 # # # # -- Hemker, P.W. & Spekreijse, S., # # Multiple Grid and Osher's Scheme for the Efficient # # Solution of the Steady Euler Equations. J. Appl. Num. # # Math. 2 (1986) 475-493. # # # # -- Hemker, P.W., # # Defect correction and higher order schemes for the # # multigrid solution of the steady Euler equations. In: # # Multi-Grid Methods II (W. Hackbusch and U. Trottenberg # # eds) pp. 149-165, Lecture Notes in Mathematics 1228, # # Springer Verlag, 1986. # # # # -- Hemker, P.W. & Johnson, G.M., # # Multigrid approaches to the Euler equations. Chapter 3 in: # Multigrid Methods, (S.F. McCormick ed.) pp. 57-72, SIAM # # Frontiers in Appl. Math. Vol.3, 1987. # # # # # # # # HISTORY: # # # # PF=EUMG001 VSN.840222 # # PF=EUMG002 VSN.840306 # # PF=EUMG003 VSN.840308 # # PF=EUMG004 VSN.840316 # # PF=EUMG005 (BSTATE) VSN.840327 # # PF=EUMG006 (2ND ORDER) VSN.840517 # # PF=EUMG007 (SEIDEL)(S) VSN.840604 # # PF=EUMG008 (ILLU ) VSN.840608 # # PF=EUMG009 (LIBRARY) VSN.840613 # # PF=EUMG010 (STORAGE) VSN.840615 # # PF=EUMG011 (U,V,C,Z) VSN.840702 # # PF=EUMG012 (LINSYS)(FTN,P) VSN.840709 # # PF=EUMG013 (MAPPING) VSN.840724 # # PF=EUMG014 (NEW FIELDS)(S) VSN.840823 # # (WALLS,FLUXT) # # PF=EUMG015 (OSH-NAT) VSN.840927 # # (MGCS) # # PF=EUMG016 (RB4/2 RELAX) VSN.841114 # # (FMG,PROLON2) # # PF=EUMG017 (RAP)(P) VSN.850102 # # (NEWTON)(P) # # (SEIDEL:LEX/RB) # # PF=EUMG018 (NEWTON IT) VSN.850322 # # (RESIDU2ND) # # PF=EUMG019 (ST.PROBLEMS) VSN.850515 # # (SMOOTH) # # PF=EUMG020 (VL 2ND ORDER) VSN.850731 # # (JACOBI) # # PF=EUMG021 (STANDARD PROBLS) VSN.850930 # # PF=EUMG022 (STRUCTS) (JR) VSN.860110 # # (ST/SF EXTENDED)(JR) # # PF=EUMG023 (TOTVAR/SETFIELD) VSN.860224 # # (PROBLEM 3B) # # (STAR ROUTINES ) # # PF=EUMG024 (MDCP) VSN.860414 # # (2ND INTERP)(P) # # # #********************************************************************# EULIB: 'PR' EJECT 'PR' 'BEGIN' # TWO-DIMENSIONAL NON-ISENTROPIC EULERIAN FLOW # 'PRIO' 'MINRESIDU' = 1, 'PLUSRESIDU' = 1, 'CORRECT' = 1, 'PRINT' = 1, -* = 4; # TYPE # 'MODE' 'FLUX' = 'STRUCT'('REAL' C1,C2,C3,C4); # 0 : # 'MODE' 'STATE' = 'STRUCT'('REAL' C1,C2,C3,C4); # 1 :(R,RU,RV,RE) # # 2 :(U, V, C, Z) # 'MODE' 'FIELD' = 'STRUCT'('REAL' DDX,DDY, 'REF''INT' TYPE, 'REF'[,]'STATE' Q); 'MODE' 'FLUXAC'= 'STRUCT'('REAL'C11,C12,C13,C14,C21,C22,C23,C24, C31,C32,C33,C34,C41,C42,C43,C44); 'MODE' 'STATAC'= 'STRUCT'('REAL'C11,C12,C13,C14,C21,C22,C23,C24, C31,C32,C33,C34,C41,C42,C43,C44); 'MODE' 'LINMAT'= 'STRUCT'('REF''REF'[,]'FLUXAC' ARRAY, 'INT' JL,JU,KL,KU,P); 'MODE' 'NETMAT'= 'REF'[]'LINMAT'; 'MODE' 'POINT' = 'STRUCT'('REAL' X, Y); 'MODE' 'WALL' = 'STRUCT'('REAL' C, S,DS); # .--------------------------------> YY # # ! -----> # # ! J # # ! ! JL . . . . . JU # # ! ! NW NO # # ! ! IL * * * * * * * # # ! ! . * * * * * * * MAPPING # # ! ! . * * * * * * * ----------> (X,Y) # # ! ! . * * * * * * * # # ! ! IU * * * * * * * # # ! ! ZW ZO # # XX ! V # # V I # # # # COMPUTATIONAL DOMAIN PHYSICAL DOMAIN # #$H# 'PROC' COMP FIELD = ('REAL' DDX,DDY, 'INT' XXL,XXU,YYL,YYU)'FIELD': #$B CMP.FLD B$# 'BEGIN' XXN:= XXL*DDX; XXZ:= XXU*DDX; YYW:= YYL*DDY; YYO:= YYU*DDY; (DDX,DDY,'HEAP''INT', 'HEAP'[(XXL+1):XXU,(YYL+1):YYU]'STATE') 'END'; #$E# #E$# 'REAL' XXN,XXZ,YYW,YYO; # GRENZEN VAN HET REKENGEBIED # 'PROC' MAPPING := ('REAL' XX,YY)'POINT': (XX,YY); # BACK-STORAGE SYSTEM # 'PR' EJECT 'PR' 'PROC' BACKSTORE = 'VOID': # PREPARES THE BACK-FILE # 'PR' XREF A68FTN,START'PR' 'SKIP'; 'PROC' PPUT = ('REF''REAL' A, 'REF''INT' SIZE,POS)'VOID': 'PR' XREF A68FTN,PUT 'PR' 'SKIP'; 'PROC' GGET = ('REF''REAL' A, 'REF''INT' SIZE,POS)'VOID': 'PR' XREF A68FTN,GET 'PR' 'SKIP'; 'PROC' GGEN = ('REF''REAL' A, 'REF''INT' SIZE,POS)'VOID': 'PR' XREF A68FTN,GEN 'PR' 'SKIP'; 'INT' LONGLINE:= 5; #$H# 'PROC' PET = ('BOOL' CHANGE, 'LINMAT' L)'VOID': #$B PET.68 B$# 'IF' 'INT'P:= P'OF'L; P>0 'THEN' 'IF' CHANGE 'THEN' 'INT' JL = JL'OF'L, KL = KL'OF'L; 'INT'S:=(KU'OF'L-KL+1)*(JU'OF'L-JL+1)*16; 'REF'[,]'FLUXAC' ARRAY = ARRAY'OF'L; PPUT(C11'OF'ARRAY[JL,KL],S,P) 'FI' ; ARRAY'OF'L:= 'REF'[,]'FLUXAC'('NIL') 'FI' ; #$E# #E$# #$H# 'PROC' GET = ('LINMAT' L)'REF'[,]'FLUXAC': #$B GET.68 B$# 'IF' 'REF'[,]'FLUXAC'(ARRAY'OF'L) 'ISNT' 'REF'[,]'FLUXAC'('NIL') 'THEN' ARRAY'OF'L 'ELSE' 'INT' JL=JL'OF'L, JU=JU'OF'L, KL=KL'OF'L, KU=KU'OF'L; 'INT' P:= P'OF'L, S:= (KU-KL+1)*(JU-JL+1)*16; 'REF'[,]'FLUXAC' ARRAY = 'HEAP'[JL:JU,KL:KU]'FLUXAC'; ARRAY'OF'L:= ARRAY; GGET(C11'OF'ARRAY[JL,KL],S,P); ARRAY 'FI' ; #$E# #E$# #$H# 'PROC' GEN = ('REF''LINMAT' L, 'INT' JL,JU,KL,KU)'REF'[,]'FLUXAC': #$B GEN.68 B$# 'BEGIN' 'INT' S,P:= 0, 'INT' LJ = JU-JL+1, 'REF'[,]'FLUXAC' ARRAY = 'HEAP'[JL:JU,KL:KU]'FLUXAC'; 'IF' LJ>LONGLINE 'THEN' S:= (KU-KL+1)*LJ*16; GGEN(C11'OF'ARRAY[JL,KL],S,P) 'FI' ; 'REF''REF'[,]'FLUXAC' ARR = 'HEAP''REF'[,]'FLUXAC':= ARRAY; L:= (ARR,JL,JU,KL,KU,P); ARRAY 'END' ; #$E# #E$# #$H# 'PROC' ATTACH = ('STRING' LFN, 'INT' NR FIELDS )'VOID': #$B ATTACH B$# 'BEGIN' OUTCHECK:= 0; OUTTOP:= NR FIELDS; WRITE:= NR FIELDS=0; ESTABLISH(FILE,LFN,BINARYSEQUENTIALCHANNEL,1,5000,3500) 'END'; #$E# #E$# 'FILE' FILE; 'BOOL' WRITE; 'INT' OUTCHECK, OUTTOP; #$H# 'PROC' OUT = ('INT' NO, 'FIELD' F)'VOID': #$B FLD.OUT B$# 'BEGIN' PRINT((NEWLINE, "OUT: " ,WHOLE(NO,0) )); 'REF' [,]'STATE' Q = Q'OF'F; 'INT' L1=1'LWB'Q, U1=1'UPB'Q, L2=2'LWB'Q, U2=2'UPB'Q, TYPE = TYPE'OF'F; 'REAL' DDX = DDX'OF'F, DDY = DDY'OF'F; PUTBIN(FILE, (NO,L1,U1,L2,U2,TYPE,DDX,DDY,NEWLINE) ); PUTBIN(FILE, (Q, 221141,NEWLINE) ) 'END'; #$E# #E$# #$H# 'PROC' IN = ('INT' NO, 'REF''FIELD' F)'VOID': #$B FLD.IN B$# 'BEGIN' PRINT((" IN : " ,WHOLE(NO,0) )); 'INT' L1,U1,L2,U2,TYPE,NOC; 'REAL' DDX, DDY; GETBIN(FILE, (NOC,L1,U1,L2,U2,TYPE,DDX,DDY,NEWLINE) ); ( NOC /= NO ! PRINT((NEWLINE,"IN-FOUT: ", NO, NOC)) ); 'HEAP' [L1:U1,L2:U2] 'STATE' Q; GETBIN(FILE, (Q, NOC,NEWLINE) ); ( NOC /= 221141 ! PRINT((NEWLINE,"IN-FOUT SLOT:", NO)) ); F:= 'FIELD'(DDX,DDY,'HEAP''INT':=TYPE,Q) 'END'; #$E# #E$# #$H# 'PROC' REWIND = 'VOID': #$B REWIND B$# 'BEGIN' ( WRITE ! OUTTOP:= OUTCHECK; PUTBIN(FILE, (NEWLINE)) ); WRITE:='FALSE'; RESET(FILE); OUTCHECK:=0 'END'; #$E# #E$# #$H# 'PROC' SKIPF = ('INT' N)'VOID': #$B SKIPF B$# 'BEGIN' 'FIELD' SF; 'IF' OUTCHECK+N <= OUTTOP 'THEN' 'TO' N 'DO' IN(OUTCHECK+:=1,SF) 'OD' 'FI' 'END' ; #$E# #E$# #$H# 'PROC' IN FRONT = ('INT' N)'VOID': #$B IN.FRNT B$# 'BEGIN' PRINT((NEWLINE, "FRONT : " ,WHOLE(N,0) )); ( N > OUTTOP ! ERROR ); ( OUTCHECK>=N ! REWIND ); SKIPF(N-1-OUTCHECK) 'END'; #$E# #E$# #$H# 'PROC' TAKE = ('INT' N, 'REF''FIELD' SF)'VOID': #$B FL.TAKE B$# 'BEGIN' IN FRONT (N); IN(OUTCHECK+:=1,SF); PRINT((" TAKEN : " ,WHOLE(OUTCHECK,0) )) 'END'; #$E# #E$# # PICTURE PREPARATION # 'PR' EJECT 'PR' 'FILE' TRANSPORT; 'PROC' TRANSA = 'VOID': ESTABLISH(TRANSPORT,"TAPE10",ZTYPE CHANNEL,1,100000,40); #$H# 'PROC' TRANSB = 'VOID': #$B TRANSB B$# 'BEGIN' 'REAL' ENDOFFILE = 77777; PUT ( TRANSPORT,(FLOAT(ENDOFFILE,12,6,2),NEWLINE) ) ; CLOSE(T RANSPORT) 'END' ; #$E# #E$# # INFOPL MEANING # # IF X-DIM AND Y-DIM > 1 # # +/-3 3D-PLOT # # +/-2 CONTOUR-PLOT # # IF X-DIM OR Y-DIM = 1 # # -1 CURVE -PLOT, NO LEGEND # # 0 CURVE -PLOT, 1ST ORDER # # 1 CURVE -PLOT, 1 DCP # # 2 CURVE -PLOT, 2 DCP # # 3 CURVE -PLOT, 3 DCP # # 4 CURVE -PLOT, 4 DCP # # 5 CURVE -PLOT, EXACT # #$H# 'PROC' TRANSV = ('INT' INFOPL, 'VAR' VARIA) 'VOID': #$B TRANSV B$# ([1:20]'REAL' DUMP, 'REF'[,]'REAL' Q = Q'OF'VARIA, 'INT' INT = (INTERN'OF'VARIA!0!1); 'REAL' R, 'INT' KEY, 'INT' L1 = 1'LWB'Q, U1 = 1'UPB'Q, L2 = 2'LWB'Q, U2 = 2'UPB'Q; 'FOR' K 'FROM' 1 'TO' 20 'DO' DUMP[K]:= 0 'OD'; 'IF' (U1-L1+1)>1 'AND' (U2-L2+1)>1 'THEN' KEY:= 54321; DUMP[1:12]:= (KEY,(U1-L1+1),(U2-L2+1),L1,U1,L2,U2, DDX'OF'VARIA,DDY'OF'VARIA,INFOPL,VAR'OF'VARIA,INT) 'ELSE' KEY:= 76543; DUMP[1:8]:= ((U1-L1+1)>1 !(KEY,(U1-L1+1),L1,U1,DDX'OF'VARIA,INFOPL, VAR'OF'VARIA,INT) !(KEY,(U2-L2+1),L2,U2,DDY'OF'VARIA,INFOPL, VAR'OF'VARIA,INT) ) 'FI' ; 'FOR' K 'TO' (KEY=54321!20!10) 'DO' ( 'ABS'(R:= DUMP[K]) >= 1.0E10 ! ERROR ! PUT(TRANSPORT,(FLOAT(R,12,6,2),NEWLINE)) ) 'OD'; 'FOR' I 'FROM' L1 'TO' U1 'DO' 'FOR' J 'FROM' L2 'TO' U2 'DO' R:= Q[I,J]; ('ABS'R <= 1.0E-10 ! R:=0 ); ('ABS'R >= 1.0E+10 ! ERROR); PUT(TRANSPORT,(FLOAT(R,12,6,2),NEWLINE)) 'OD' 'OD'); #$E# #E$# # OPERATORS # 'PR' EJECT 'PR' 'OP' ** = ('REAL' A,B)'REAL': EXP( B*LN(A) ); 'OP' - = ('FLUX' F)'FLUX':(-C1'OF'F,-C2'OF'F,-C3'OF'F,-C4'OF'F); 'OP' + = ('FLUX' A,B)'FLUX': (C1'OF'A+C1'OF'B,C2'OF'A+C2'OF'B,C3'OF'A+C3'OF'B,C4'OF'A+C4'OF'B); 'OP' - = ('FLUX' A,B)'FLUX': (C1'OF'A-C1'OF'B,C2'OF'A-C2'OF'B,C3'OF'A-C3'OF'B,C4'OF'A-C4'OF'B); 'OP' * = ('FLUX' A,B)'REAL': C1'OF'A*C1'OF'B+C2'OF'A*C2'OF'B+C3'OF'A*C3'OF'B+C4'OF'A*C4'OF'B ; 'OP' * = ('REAL' B, 'FLUX' A)'FLUX': ( B*C1'OF'A , B*C2'OF'A , B*C3'OF'A , B*C4'OF'A ); 'OP' / = ('FLUX' A, 'REAL' B)'FLUX': ( C1'OF'A/B , C2'OF'A/B , C3'OF'A/B , C4'OF'A/B ) ; 'OP' +:= = ('REF''FLUX' A, 'FLUX' B)'REF''FLUX': ((C1'OF'A +:= C1'OF'B , C2'OF'A +:= C2'OF'B , C3'OF'A +:= C3'OF'B , C4'OF'A +:= C4'OF'B );A ) ; 'OP' -:= = ('REF''FLUX' A, 'FLUX' B)'REF''FLUX': ((C1'OF'A -:= C1'OF'B , C2'OF'A -:= C2'OF'B , C3'OF'A -:= C3'OF'B , C4'OF'A -:= C4'OF'B );A ) ; 'OP' *:= = ('REF''FLUX' A, 'REAL' B)'REF''FLUX': ((C1'OF'A*:=B , C2'OF'A*:=B , C3'OF'A*:=B , C4'OF'A*:=B);A); 'OP' /:= = ('REF''FLUX' A, 'REAL' B)'REF''FLUX': ((C1'OF'A/:=B , C2'OF'A/:=B , C3'OF'A/:=B , C4'OF'A/:=B);A); 'OP' 'ABS' = ('FLUX' A)'REAL': 'ABS'C1'OF'A+ 'ABS'C2'OF'A+ 'ABS'C3'OF'A+ 'ABS'C4'OF'A ; #$H# 'OP' - = ('FLUXAC' A)'FLUXAC' : #$B MINFA B$# 'FLUXAC'( -C11'OF'A,-C12'OF'A,-C13'OF'A,-C14'OF'A , -C21'OF'A,-C22'OF'A,-C23'OF'A,-C24'OF'A , -C31'OF'A,-C32'OF'A,-C33'OF'A,-C34'OF'A , -C41'OF'A,-C42'OF'A,-C43'OF'A,-C44'OF'A ) ; #$E# #E$# #$H# 'OP' + = ('FLUXAC' A,B)'FLUXAC': #$B FAPLFA B$# ( C11'OF'A + C11'OF'B , C12'OF'A + C12'OF'B , C13'OF'A + C13'OF'B , C14'OF'A + C14'OF'B , C21'OF'A + C21'OF'B , C22'OF'A + C22'OF'B , C23'OF'A + C23'OF'B , C24'OF'A + C24'OF'B , C31'OF'A + C31'OF'B , C32'OF'A + C32'OF'B , C33'OF'A + C33'OF'B , C34'OF'A + C34'OF'B , C41'OF'A + C41'OF'B , C42'OF'A + C42'OF'B , C43'OF'A + C43'OF'B , C44'OF'A + C44'OF'B ) ; #$E# #E$# #$H# 'OP' - = ('FLUXAC' A,B)'FLUXAC': #$B FAMNFA B$# ( C11'OF'A - C11'OF'B , C12'OF'A - C12'OF'B , C13'OF'A - C13'OF'B , C14'OF'A - C14'OF'B , C21'OF'A - C21'OF'B , C22'OF'A - C22'OF'B , C23'OF'A - C23'OF'B , C24'OF'A - C24'OF'B , C31'OF'A - C31'OF'B , C32'OF'A - C32'OF'B , C33'OF'A - C33'OF'B , C34'OF'A - C34'OF'B , C41'OF'A - C41'OF'B , C42'OF'A - C42'OF'B , C43'OF'A - C43'OF'B , C44'OF'A - C44'OF'B ) ; #$E# #E$# #$H# 'OP' * = ('REAL' B, 'FLUXAC' A)'FLUXAC': #$B RTIMFA B$# ( B*C11'OF'A,B*C12'OF'A,B*C13'OF'A,B*C14'OF'A , B*C21'OF'A,B*C22'OF'A,B*C23'OF'A,B*C24'OF'A , B*C31'OF'A,B*C32'OF'A,B*C33'OF'A,B*C34'OF'A , B*C41'OF'A,B*C42'OF'A,B*C43'OF'A,B*C44'OF'A ) ; #$E# #E$# #$H# 'OP' / = ('FLUXAC' A,'REAL' B )'FLUXAC': #$B FADIVR B$# ( C11'OF'A/B,C12'OF'A/B,C13'OF'A/B,C14'OF'A/B , C21'OF'A/B,C22'OF'A/B,C23'OF'A/B,C24'OF'A/B , C31'OF'A/B,C32'OF'A/B,C33'OF'A/B,C34'OF'A/B , C41'OF'A/B,C42'OF'A/B,C43'OF'A/B,C44'OF'A/B ) ; #$E# #E$# #$H# 'OP' +:= = ('REF''FLUXAC' A, 'FLUXAC' B)'REF''FLUXAC': #$B FAPBFA B$# 'BEGIN' (C11'OF'A +:= C11'OF'B , C12'OF'A +:= C12'OF'B , C13'OF'A +:= C13'OF'B , C14'OF'A +:= C14'OF'B , C21'OF'A +:= C21'OF'B , C22'OF'A +:= C22'OF'B , C23'OF'A +:= C23'OF'B , C24'OF'A +:= C24'OF'B , C31'OF'A +:= C31'OF'B , C32'OF'A +:= C32'OF'B , C33'OF'A +:= C33'OF'B , C34'OF'A +:= C34'OF'B , C41'OF'A +:= C41'OF'B , C42'OF'A +:= C42'OF'B , C43'OF'A +:= C43'OF'B , C44'OF'A +:= C44'OF'B );A 'END' ; #$E# #E$# #$H# 'OP' -:= = ('REF''FLUXAC' A, 'FLUXAC' B)'REF''FLUXAC': #$B FAMBFA B$# 'BEGIN' (C11'OF'A -:= C11'OF'B , C12'OF'A -:= C12'OF'B , C13'OF'A -:= C13'OF'B , C14'OF'A -:= C14'OF'B , C21'OF'A -:= C21'OF'B , C22'OF'A -:= C22'OF'B , C23'OF'A -:= C23'OF'B , C24'OF'A -:= C24'OF'B , C31'OF'A -:= C31'OF'B , C32'OF'A -:= C32'OF'B , C33'OF'A -:= C33'OF'B , C34'OF'A -:= C34'OF'B , C41'OF'A -:= C41'OF'B , C42'OF'A -:= C42'OF'B , C43'OF'A -:= C43'OF'B , C44'OF'A -:= C44'OF'B );A 'END' ; #$E# #E$# #$H# 'OP' *:= = ('REF''FLUXAC' A, 'REAL' B)'REF''FLUXAC': #$B FATBRL B$# 'BEGIN' (C11'OF'A*:=B,C12'OF'A*:=B,C13'OF'A*:=B,C14'OF'A*:=B, C21'OF'A*:=B,C22'OF'A*:=B,C23'OF'A*:=B,C24'OF'A*:=B, C31'OF'A*:=B,C32'OF'A*:=B,C33'OF'A*:=B,C34'OF'A*:=B, C41'OF'A*:=B,C42'OF'A*:=B,C43'OF'A*:=B,C44'OF'A*:=B);A 'END' ; #$E# #E$# #$H# 'OP' /:= = ('REF''FLUXAC' A, 'REAL' B)'REF''FLUXAC': #$B FADBRL B$# 'BEGIN' (C11'OF'A/:=B,C12'OF'A/:=B,C13'OF'A/:=B,C14'OF'A/:=B, C21'OF'A/:=B,C22'OF'A/:=B,C23'OF'A/:=B,C24'OF'A/:=B, C31'OF'A/:=B,C32'OF'A/:=B,C33'OF'A/:=B,C34'OF'A/:=B, C41'OF'A/:=B,C42'OF'A/:=B,C43'OF'A/:=B,C44'OF'A/:=B);A 'END' ; #$E# #E$# #$H# 'OP' * = ('FLUXAC' A, 'FLUX' B)'FLUX': #$B FATIMF B$# 'BEGIN' 'REAL' B1=C1'OF'B,B2=C2'OF'B,B3=C3'OF'B,B4=C4'OF'B ; 'FLUX'(C11'OF'A*B1+C12'OF'A*B2+C13'OF'A*B3+C14'OF'A*B4, C21'OF'A*B1+C22'OF'A*B2+C23'OF'A*B3+C24'OF'A*B4, C31'OF'A*B1+C32'OF'A*B2+C33'OF'A*B3+C34'OF'A*B4, C41'OF'A*B1+C42'OF'A*B2+C43'OF'A*B3+C44'OF'A*B4) 'END' ; #$E# #E$# #$H# 'OP' * = ('FLUX' B, 'FLUXAC' A)'FLUX': #$B FTIMFA B$# 'BEGIN' 'REAL' B1=C1'OF'B,B2=C2'OF'B,B3=C3'OF'B,B4=C4'OF'B ; 'FLUX'(B1*C11'OF'A+B2*C21'OF'A+B3*C31'OF'A+B4*C41'OF'A, B1*C12'OF'A+B2*C22'OF'A+B3*C32'OF'A+B4*C42'OF'A, B1*C13'OF'A+B2*C23'OF'A+B3*C33'OF'A+B4*C43'OF'A, B1*C14'OF'A+B2*C24'OF'A+B3*C34'OF'A+B4*C44'OF'A) 'END' ; #$E# #E$# #$H# 'OP' * = ('FLUXAC' A,B)'FLUXAC' : #$B FATIFA B$# 'FLUXAC' ( C11'OF'A*C11'OF'B+C12'OF'A*C21'OF'B+C13'OF'A*C31'OF'B+C14'OF'A*C41'OF'B, C11'OF'A*C12'OF'B+C12'OF'A*C22'OF'B+C13'OF'A*C32'OF'B+C14'OF'A*C42'OF'B, C11'OF'A*C13'OF'B+C12'OF'A*C23'OF'B+C13'OF'A*C33'OF'B+C14'OF'A*C43'OF'B, C11'OF'A*C14'OF'B+C12'OF'A*C24'OF'B+C13'OF'A*C34'OF'B+C14'OF'A*C44'OF'B, C21'OF'A*C11'OF'B+C22'OF'A*C21'OF'B+C23'OF'A*C31'OF'B+C24'OF'A*C41'OF'B, C21'OF'A*C12'OF'B+C22'OF'A*C22'OF'B+C23'OF'A*C32'OF'B+C24'OF'A*C42'OF'B, C21'OF'A*C13'OF'B+C22'OF'A*C23'OF'B+C23'OF'A*C33'OF'B+C24'OF'A*C43'OF'B, C21'OF'A*C14'OF'B+C22'OF'A*C24'OF'B+C23'OF'A*C34'OF'B+C24'OF'A*C44'OF'B, C31'OF'A*C11'OF'B+C32'OF'A*C21'OF'B+C33'OF'A*C31'OF'B+C34'OF'A*C41'OF'B, C31'OF'A*C12'OF'B+C32'OF'A*C22'OF'B+C33'OF'A*C32'OF'B+C34'OF'A*C42'OF'B, C31'OF'A*C13'OF'B+C32'OF'A*C23'OF'B+C33'OF'A*C33'OF'B+C34'OF'A*C43'OF'B, C31'OF'A*C14'OF'B+C32'OF'A*C24'OF'B+C33'OF'A*C34'OF'B+C34'OF'A*C44'OF'B, C41'OF'A*C11'OF'B+C42'OF'A*C21'OF'B+C43'OF'A*C31'OF'B+C44'OF'A*C41'OF'B, C41'OF'A*C12'OF'B+C42'OF'A*C22'OF'B+C43'OF'A*C32'OF'B+C44'OF'A*C42'OF'B, C41'OF'A*C13'OF'B+C42'OF'A*C23'OF'B+C43'OF'A*C33'OF'B+C44'OF'A*C43'OF'B, C41'OF'A*C14'OF'B+C42'OF'A*C24'OF'B+C43'OF'A*C34'OF'B+C44'OF'A*C44'OF'B ) ; #$E# #E$# #$H# 'OP' -* = ('FLUXAC' A,B)'FLUXAC': #$B FAMTFA B$# 'BEGIN' -(A*B) 'END' ; #$E# #E$# #$H# 'OP' / = ('INT' ONE, 'FLUXAC' A)'FLUXAC': #$B GAUSIN B$# 'BEGIN' # 4*4 MATRIX INVERSION # (ONE /= 1 ! ERROR ); 'PROC' INV = ('REF''REAL' A, B, PIVOT)'VOID': 'PR' XREF A68FTN,INVERT 'PR' 'SKIP'; [1:4,1:4]'REAL' AA,B ;'REAL' PIVOT ; AA := ( (C11'OF'A,C21'OF'A,C31'OF'A,C41'OF'A) , (C12'OF'A,C22'OF'A,C32'OF'A,C42'OF'A) , (C13'OF'A,C23'OF'A,C33'OF'A,C43'OF'A) , (C14'OF'A,C24'OF'A,C34'OF'A,C44'OF'A) ); INV(AA[1,1],B[1,1],PIVOT); (PIVOT<1.0E-5 ! PRINT((NEWLINE," MAYBE SINGULAR MATRIX,", " PIVOT = ",PIVOT)); (PIVOT=0!ERROR) ); 'FLUXAC' ( B[1,1],B[2,1],B[3,1],B[4,1], B[1,2],B[2,2],B[3,2],B[4,2], B[1,3],B[2,3],B[3,3],B[4,3], B[1,4],B[2,4],B[3,4],B[4,4] ) 'END';#$E# #E$# #$H# 'OP' / = ('FLUX' Y, 'FLUXAC' A)'STATE': #$B GAUSEL B$# 'BEGIN' 'PROC' GAUSEL = ('REF''REAL' C,X, PIVOT)'VOID': 'PR' XREF A68FTN,GAUSS 'PR' 'SKIP'; [1:4] 'REAL' X, [1:5,1:4] 'REAL' C, 'REAL' PIVOT; # A68 <-> FORTRAN TRANPOSITION OF ROWS AND COLUMNS. # C := ( (C11'OF'A,C21'OF'A,C31'OF'A,C41'OF'A), (C12'OF'A,C22'OF'A,C32'OF'A,C42'OF'A), (C13'OF'A,C23'OF'A,C33'OF'A,C43'OF'A), (C14'OF'A,C24'OF'A,C34'OF'A,C44'OF'A), (C1 'OF'Y,C2 'OF'Y,C3 'OF'Y,C4 'OF'Y) ); GAUSEL(C[1,1],X[1],PIVOT); (PIVOT<1.0E-5 ! PRINT((NEWLINE," MAYBE SINGULAR MATRIX,", " PIVOT = ",PIVOT)); (PIVOT=0!ERROR)); 'STATE' ( X[1],X[2],X[3],X[4] ) 'END'; #$E# #E$# #$H# 'OP' / = ('REF'[]'FLUX' B, 'REF'[,]'FLUXAC' A)'REF'[]'FLUX': #$B DIV4 B$# 'BEGIN' [1:4,1:4]'FLUXAC' MAT; 'HEAP'[1:4]'FLUX' VEC; 'PROC' SOL4 = ('REF'[,]'FLUXAC' A, 'REF'[]'FLUX' R)'VOID': # NB: OVERSCHRIJFT A EN R !!!!!!!!!!!!!!! # 'BEGIN' 'FOR' I 'TO' 4 'DO' A[I,I]:= 1/A[I,I]; 'FOR' J 'FROM' I+1 'TO' 4 'DO' A[I,J] := A[I,I]*A[I,J]; 'FOR' K 'FROM' I+1 'TO' 4 'DO' A[K,J]-:= A[K,I]*A[I,J] 'OD' 'OD'; R[I ] := A[I,I]*R[I ]; 'FOR' K 'FROM' I+1 'TO' 4 'DO' R[K ]-:= A[K,I]*R[I ] 'OD' 'OD'; 'FOR' I 'FROM' 3 'BY' -1 'TO' 1 'DO''FOR' J 'FROM' I+1 'TO' 4 'DO' R[I ]-:= A[I,J]*R[J ] 'OD''OD' 'END'; 'FOR' I 'TO' 4 'DO' VEC[I]:= B[I]; 'FOR' J 'TO' 4 'DO' MAT[I,J]:= A[I,J] 'OD' 'OD'; SOL4(MAT,VEC); VEC 'END'; #$E# #E$# # COLUMN HANDLING PROCEDURES # 'PR' EJECT 'PR' 'PRIO' 'FC'=9 ; #$H# 'OP' 'FC' = ('INT' I,'FLUXAC'A)'FLUX': # FREE A COLUMN # #$B FRE.COL B$# (I! 'FLUX' (C11'OF'A,C21'OF'A,C31'OF'A,C41'OF'A), 'FLUX' (C12'OF'A,C22'OF'A,C32'OF'A,C42'OF'A), 'FLUX' (C13'OF'A,C23'OF'A,C33'OF'A,C43'OF'A), 'FLUX' (C14'OF'A,C24'OF'A,C34'OF'A,C44'OF'A) ! NUL ) ; #$E# #E$# # ASSIGNMENT OF FLUX B TO COLUMN I OF A # #$H# 'PROC' CASS = ('REF''FLUXAC'A,'INT'I,'FLUX'B)'REF''FLUXAC': #$B COL.ASS B$# ('IF' A'ISNT''NIL' 'THEN''REAL'B1=C1'OF'B,B2=C2'OF'B,B3=C3'OF'B,B4=C4'OF'B; (I!((C11'OF'A:=B1,C21'OF'A:=B2,C31'OF'A:=B3,C41'OF'A:=B4)), ((C12'OF'A:=B1,C22'OF'A:=B2,C32'OF'A:=B3,C42'OF'A:=B4)), ((C13'OF'A:=B1,C23'OF'A:=B2,C33'OF'A:=B3,C43'OF'A:=B4)), ((C14'OF'A:=B1,C24'OF'A:=B2,C34'OF'A:=B3,C44'OF'A:=B4)) ! PRINT((NEWLINE,"WRONG INDEX IN -CASS- CALL"));ERROR ) 'FI';A) ; #$E# #E$# # ADDITION OF FLUX B TO COLUMN I OF A # #$H# 'PROC' CADD = ('REF''FLUXAC'A,'INT'I,'FLUX'B)'REF''FLUXAC': #$B COL.ADD B$# ('IF' A'ISNT''NIL' 'THEN''REAL'B1=C1'OF'B,B2=C2'OF'B,B3=C3'OF'B,B4=C4'OF'B; (I!((C11'OF'A+:=B1,C21'OF'A+:=B2,C31'OF'A+:=B3,C41'OF'A+:=B4)), ((C12'OF'A+:=B1,C22'OF'A+:=B2,C32'OF'A+:=B3,C42'OF'A+:=B4)), ((C13'OF'A+:=B1,C23'OF'A+:=B2,C33'OF'A+:=B3,C43'OF'A+:=B4)), ((C14'OF'A+:=B1,C24'OF'A+:=B2,C34'OF'A+:=B3,C44'OF'A+:=B4)) ! PRINT((NEWLINE,"WRONG INDEX IN -CADD- CALL"));ERROR ) 'FI';A) ; #$E# #E$# # SUBSTRACTION OF FLUX B TO COLUMN I OF A # #$H# 'PROC' CSUB = ('REF''FLUXAC'A,'INT'I,'FLUX'B)'REF''FLUXAC': #$B COL.SUB B$# ('IF' A'ISNT''NIL' 'THEN''REAL'B1=C1'OF'B,B2=C2'OF'B,B3=C3'OF'B,B4=C4'OF'B; (I!((C11'OF'A-:=B1,C21'OF'A-:=B2,C31'OF'A-:=B3,C41'OF'A-:=B4)), ((C12'OF'A-:=B1,C22'OF'A-:=B2,C32'OF'A-:=B3,C42'OF'A-:=B4)), ((C13'OF'A-:=B1,C23'OF'A-:=B2,C33'OF'A-:=B3,C43'OF'A-:=B4)), ((C14'OF'A-:=B1,C24'OF'A-:=B2,C34'OF'A-:=B3,C44'OF'A-:=B4)) ! PRINT((NEWLINE,"WRONG INDEX IN -CSUB- CALL"));ERROR ) 'FI';A) ; #$E# #E$# # MULTIPLICATION OF COLUMN I OF A BY A FACTOR B # #$H# 'PROC' CMUL = ('REF''FLUXAC'A,'INT'I,'REAL'B)'REF''FLUXAC': #$B COL.MUL B$# ( 'IF' (A'ISNT''NIL')'AND'(B/=1.0) 'THEN' (I!((C11'OF'A*:=B,C21'OF'A*:=B,C31'OF'A*:=B,C41'OF'A*:=B)), ((C12'OF'A*:=B,C22'OF'A*:=B,C32'OF'A*:=B,C42'OF'A*:=B)), ((C13'OF'A*:=B,C23'OF'A*:=B,C33'OF'A*:=B,C43'OF'A*:=B)), ((C14'OF'A*:=B,C24'OF'A*:=B,C34'OF'A*:=B,C44'OF'A*:=B)) ! PRINT((NEWLINE,"WRONG INDEX IN -CMUL- CALL"));ERROR ) 'FI';A) ; #$E# #E$# # DIVISION OF COLUMN I OF A BY A FACTOR B # #$H# 'PROC' CDIV = ('REF''FLUXAC'A,'INT'I,'REAL'B)'REF''FLUXAC': #$B COL.DIV B$# ('IF' (A'ISNT''NIL')'AND'(B/=0.0)'AND'(B/=1.0) 'THEN' (I!((C11'OF'A/:=B,C21'OF'A/:=B,C31'OF'A/:=B,C41'OF'A/:=B)), ((C12'OF'A/:=B,C22'OF'A/:=B,C32'OF'A/:=B,C42'OF'A/:=B)), ((C13'OF'A/:=B,C23'OF'A/:=B,C33'OF'A/:=B,C43'OF'A/:=B)), ((C14'OF'A/:=B,C24'OF'A/:=B,C34'OF'A/:=B,C44'OF'A/:=B)) ! PRINT((NEWLINE,"WRONG INDEX IN -CDIV- CALL"));ERROR ) 'FI';A) ; #$E# #E$# # ROW HANDLING PROCEDURES # 'PR' EJECT 'PR' 'PRIO' 'FR'=9 ; #$H# 'OP' 'FR' = ('INT' I,'FLUXAC'A)'FLUX': # FREE A ROW # #$B FRE.ROW B$# (I! 'FLUX' (C11'OF'A,C12'OF'A,C13'OF'A,C14'OF'A), 'FLUX' (C21'OF'A,C22'OF'A,C23'OF'A,C24'OF'A), 'FLUX' (C31'OF'A,C32'OF'A,C33'OF'A,C34'OF'A), 'FLUX' (C41'OF'A,C42'OF'A,C43'OF'A,C44'OF'A) ! NUL ) ; #$E# #E$# # ASSIGNMENT OF FLUX B TO ROW I OF A # #$H# 'PROC' RASS = ('REF''FLUXAC'A,'INT'I,'FLUX'B)'REF''FLUXAC': #$B ROW.ASS B$# ('IF' A'ISNT''NIL' 'THEN''REAL'B1=C1'OF'B,B2=C2'OF'B,B3=C3'OF'B,B4=C4'OF'B; (I!((C11'OF'A:=B1,C12'OF'A:=B2,C13'OF'A:=B3,C14'OF'A:=B4)), ((C21'OF'A:=B1,C22'OF'A:=B2,C23'OF'A:=B3,C24'OF'A:=B4)), ((C31'OF'A:=B1,C32'OF'A:=B2,C33'OF'A:=B3,C34'OF'A:=B4)), ((C41'OF'A:=B1,C42'OF'A:=B2,C43'OF'A:=B3,C44'OF'A:=B4)) ! PRINT((NEWLINE,"WRONG INDEX IN -RASS- CALL"));ERROR ) 'FI';A) ; #$E# #E$# # ADDITION OF FLUX B TO ROW I OF A # #$H# 'PROC' RADD = ('REF''FLUXAC'A,'INT'I,'FLUX'B)'REF''FLUXAC': #$B ROW.ADD B$# ('IF' A'ISNT''NIL' 'THEN''REAL'B1=C1'OF'B,B2=C2'OF'B,B3=C3'OF'B,B4=C4'OF'B; (I!((C11'OF'A+:=B1,C12'OF'A+:=B2,C13'OF'A+:=B3,C14'OF'A+:=B4)), ((C21'OF'A+:=B1,C22'OF'A+:=B2,C23'OF'A+:=B3,C24'OF'A+:=B4)), ((C31'OF'A+:=B1,C32'OF'A+:=B2,C33'OF'A+:=B3,C34'OF'A+:=B4)), ((C41'OF'A+:=B1,C42'OF'A+:=B2,C43'OF'A+:=B3,C44'OF'A+:=B4)) ! PRINT((NEWLINE,"WRONG INDEX IN -RADD- CALL"));ERROR ) 'FI';A) ; #$E# #E$# # SUBSTRACTION OF FLUX B TO ROW I OF A # #$H# 'PROC' RSUB = ('REF''FLUXAC'A,'INT'I,'FLUX'B)'REF''FLUXAC': #$B ROW.SUB B$# ('IF' A'ISNT''NIL' 'THEN''REAL'B1=C1'OF'B,B2=C2'OF'B,B3=C3'OF'B,B4=C4'OF'B; (I!((C11'OF'A-:=B1,C12'OF'A-:=B2,C13'OF'A-:=B3,C14'OF'A-:=B4)), ((C21'OF'A-:=B1,C22'OF'A-:=B2,C23'OF'A-:=B3,C24'OF'A-:=B4)), ((C31'OF'A-:=B1,C32'OF'A-:=B2,C33'OF'A-:=B3,C34'OF'A-:=B4)), ((C41'OF'A-:=B1,C42'OF'A-:=B2,C43'OF'A-:=B3,C44'OF'A-:=B4)) ! PRINT((NEWLINE,"WRONG INDEX IN -RSUB- CALL"));ERROR ) 'FI';A) ; #$E# #E$# # MULTIPLICATION OF ROW I OF A BY A FACTOR B # #$H# 'PROC' RMUL = ('REF''FLUXAC'A,'INT'I,'REAL'B)'REF''FLUXAC': #$B ROW.MUL B$# ( 'IF' (A'ISNT''NIL')'AND'(B/=1.0) 'THEN' (I!((C11'OF'A*:=B,C12'OF'A*:=B,C13'OF'A*:=B,C14'OF'A*:=B)), ((C21'OF'A*:=B,C22'OF'A*:=B,C23'OF'A*:=B,C24'OF'A*:=B)), ((C31'OF'A*:=B,C32'OF'A*:=B,C33'OF'A*:=B,C34'OF'A*:=B)), ((C41'OF'A*:=B,C42'OF'A*:=B,C43'OF'A*:=B,C44'OF'A*:=B)) ! PRINT((NEWLINE,"WRONG INDEX IN -RMUL- CALL"));ERROR ) 'FI';A) ; #$E# #E$# # DIVISION OF ROW I OF A BY A FACTOR B # #$H# 'PROC' RDIV = ('REF''FLUXAC'A,'INT'I,'REAL'B)'REF''FLUXAC': #$B ROW.DIV B$# ( 'IF' (A'ISNT''NIL')'AND'(B/=0.0)'AND'(B/=1.0) 'THEN' (I!((C11'OF'A/:=B,C12'OF'A/:=B,C13'OF'A/:=B,C14'OF'A/:=B)), ((C21'OF'A/:=B,C22'OF'A/:=B,C23'OF'A/:=B,C24'OF'A/:=B)), ((C31'OF'A/:=B,C32'OF'A/:=B,C33'OF'A/:=B,C34'OF'A/:=B)), ((C41'OF'A/:=B,C42'OF'A/:=B,C43'OF'A/:=B,C44'OF'A/:=B)) ! PRINT((NEWLINE,"WRONG INDEX IN -RDIV- CALL"));ERROR ) 'FI';A) ; #$E# #E$# # ILLU RELAX # 'PR' EJECT 'PR' #$H# 'PROC' ILLU RELAX = ('NETMAT' DEC,JAC, 'FIELD' QF,RF )'VOID': #$B ILLUR B$# 'BEGIN' 'REF'[,]'STATE'QQ = Q'OF'QF, RH=Q'OF'RF; 'INT' L1= 1'LWB'QQ, U1= 1'UPB'QQ, L2= 2'LWB'QQ, U2= 2'UPB'QQ; 'LOC'[L2:U2]'STATE' DUI; 'PROC' SOLL = ('REF'[]'FLUXAC' L,D,U, 'REF'[]'STATE' Z )'VOID': ('FOR'J'FROM' L2+1 'TO' U2 'DO' Z[J]+:= L[J]*Z[J-1] 'OD'; 'FOR'J'FROM' L2 'TO' U2 'DO' Z[J] := D[J]*Z[J ] 'OD'; 'FOR'J'FROM' U2-1 'BY' -1 'TO' L2 'DO' Z[J]+:= U[J]*Z[J+1] 'OD' ); ( 'REF'[,]'FLUXAC' DECL1 = GET(DEC[L1]); SOLL(DECL1[,-1],DECL1[,0],DECL1[,1],RH[L1,]); PET('FALSE',DEC[L1]) ); 'FOR' I 'FROM' L1+1 'TO' U1 'DO' 'REF'[,]'FLUXAC' DECI = GET(DEC[I]), JACI = GET(JAC[I]); 'FOR' J 'FROM' L2 'TO' U2 'DO' RH[I,J]-:= JACI[J,-2]*RH[I-1,J ] 'OD'; SOLL(DECI[,-1],DECI[,0],DECI[,1],RH[I,]); PET('FALSE',DEC[I]); PET('FALSE',JAC[I]) 'OD'; 'FOR' J 'FROM' L2 'TO' U2 'DO' QQ[U1,J]+:= ( DUI[J] := RH[U1,J] ) 'OD'; 'FOR' I 'FROM' U1-1 'BY' -1 'TO' L1 'DO' 'REF'[,]'FLUXAC' DECI = GET(DEC[I]), JACI = GET(JAC[I]); 'FOR' J 'FROM' L2 'TO' U2 'DO' DUI[J] := JACI[J, 2]*DUI[J] 'OD'; SOLL(DECI[,-1],DECI[,0],DECI[,1],DUI); PET('FALSE',DEC[I]); PET('FALSE',JAC[I]); 'FOR' J 'FROM' L2 'TO' U2 'DO' QQ[I,J]+:= ( DUI[J] := RH[I,J]-DUI[J] )'OD' 'OD' 'END'; #$E# #E$# #$H# 'PROC' ILLUDEC = ('NETMAT' JAC, 'REF''NETMAT' DECN )'VOID': # [L1:U1,L2:U2,-1:+1]'FLUXAC' DEC # # [L1:U1]'LINMAT' DEC # #$B ILLUD B$# 'BEGIN' 'INT' L1 = 'LWB'JAC, U1 = 'UPB'JAC; 'INT' L2 = JL'OF'JAC[L1], U2 = JU'OF'JAC[L1]; DECN:= 'HEAP'[L1:U1]'LINMAT'; [L2:U2,-1:+1]'FLUXAC' D, DINV; 'NETMAT' DEC = DECN; (L1/='LWB'DEC 'OR' U1/='UPB'DEC ! ERROR); 'PROC' DECC = ('REF'[,]'FLUXAC' DEC, D)'VOID': ( DEC[L2, 0]:= 1/ D[L2, 0]; 'FOR' J 'FROM' L2+1 'TO' U2 'DO' DEC[J ,-1]:= D[J ,-1] -*DEC[J-1, 0]; DEC[J , 0]:= 1/(D[J , 0] + DEC[J ,-1]*D[J-1,1]) 'OD' ; 'FOR' J 'FROM' L2 'TO' U2-1 'DO' DEC[J ,+1]:= DEC[J,0] -*D[J,1] 'OD' ); ('REF'[,]'FLUXAC' JACL1 = GET(JAC[L1]), DECL1 = GEN(DEC[L1],L2,U2,-1,1); DECC(DECL1,JACL1) ); 'FOR' I 'FROM' L1 'TO' U1-1 'DO' 'REF'[,]'FLUXAC' JACI = GET(JAC[I ]), JACIP = GET(JAC[I+1]), DECI = GET(DEC[I ]), DECIP = GEN(DEC[I+1],L2,U2,-1,1); DINV[U2,0] := DECI[U2,0]; 'FOR' J 'FROM' U2-1 'BY' -1 'TO' L2 'DO' DINV[ J,0]:= DECI[J,0] +DECI[J,1]*DINV[J+1,0]*DECI[J+1,-1] 'OD'; 'FOR' J 'FROM' L2+1 'TO' U2 'DO' DINV[J ,-1]:= DINV[J ,0]*DECI[J,-1]; DINV[J-1, 1]:= DECI[J-1,1]*DINV[J, 0] 'OD'; 'FOR' K 'FROM' -1 'TO' 1 'DO' 'FOR' J 'FROM' (K=-1!L2+1!L2) 'TO' (K=1!U2-1!U2) 'DO' D[J,K] := JACIP[J, K] -JACIP[J,-2]*DINV[J,K]*JACI[J+K,2] 'OD' 'OD'; DECC(DECIP,D); PET('FALSE',JAC[ I]); PET('TRUE',DEC[ I]) 'OD' ; PET('FALSE',JAC[U1]); PET('TRUE',DEC[U1]) 'END'; #$E# #E$# # CONSTANTS AND 'NUL' # 'PR' EJECT 'PR' 'REAL' GAM = 1.4 ; 'REAL' OGAM = 1/GAM , OTGAM = 0.5/ GAM , GAMM1 = GAM-1 , OGAMM1 = 1/(GAM-1) , HGAMM1 = (GAM-1)/2 , TOGAMM1 = 2/(GAM-1) , TOGAMP1 = 2/(GAM+1) , OGAMGAMM1 = 1/(GAM*(GAM-1)) , GAMM1OGAMP1= (GAM-1)/(GAM+1) ; 'FLUX' NUL = (0,0,0,0); 'FLUXAC' NULFLUXAC = (0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0 ); 'OP' 'NUL' = ('REF''FLUXAC' A)'VOID' : A:=NULFLUXAC ; #$H# 'OP' 'NUL' = ('REF'[,]'STATE' A )'VOID': #$B NL.RRS B$# 'BEGIN' 'INT' IL=1'LWB'A,IU=1'UPB'A,JL=2'LWB'A,JU=2'UPB'A ; 'FOR' I 'FROM' IL 'TO' IU 'DO' 'REF'[]'STATE'AI=A[I,]; 'FOR' J 'FROM' JL 'TO' JU 'DO' AI[J] := NUL 'OD' 'OD' 'END' ; #$E# #E$# #$H# 'OP' 'NUL' = ('REF'[]'FLUXAC' A )'VOID': #$B NL.RFA B$# 'FOR' I 'FROM' 'LWB'A 'TO' 'UPB'A 'DO' A[I] := NULFLUXAC 'OD' ; #$E# #E$# #$H# 'OP' 'NUL' = ('REF'[,]'FLUXAC' A )'VOID': #$B NL.RRFA B$# 'BEGIN' 'INT' IL=1'LWB'A,IU=1'UPB'A,JL=2'LWB'A,JU=2'UPB'A ; 'FOR' I 'FROM' IL 'TO' IU 'DO' 'REF'[]'FLUXAC'AI=A[I,]; 'FOR' J 'FROM' JL 'TO' JU 'DO' AI[J] := NULFLUXAC 'OD' 'OD' 'END' ; #$E# #E$# #$H# 'OP' 'NUL' = ('FIELD' A)'FIELD': #$B NL.FLD B$# 'BEGIN' 'NUL'(Q'OF'A) ; A 'END' ; #$E# #E$# #$H# 'OP' 'NULRHS' = ('FIELD' Q)'FIELD': #$B NL.RHS B$# 'BEGIN' 'FIELD' F = 'NUL''COPY'Q ; TYPE'OF'F := 0 ; F 'END' ; #$E# #E$# 'PR' EJECT 'PR' #$H# 'OP' 'COPY' = ('FIELD' F)'FIELD': #$B COPY.FL B$# 'BEGIN' 'REF'[,]'STATE' Q = Q'OF'F; (DDX'OF'F, DDY'OF'F, 'HEAP''INT':= TYPE'OF'F, 'HEAP'[1'LWB'Q:1'UPB'Q,2'LWB'Q:2'UPB'Q]'STATE':= Q ) 'END'; #$E# #E$# #$H# 'OP' 'NORM' = ('FIELD'F)[,]'REAL': #$B NORM B$# 'BEGIN' 'REF'[,]'STATE' Q = (Q'OF'F)[@1,@1]; 'INT' N1=1'UPB'Q, N2=2'UPB'Q; [1:4,-1:2]'REAL' NRM; 'FOR' K 'TO' 4 'DO' 'REF'[,]'REAL' QK=(K!C1'OF'Q,C2'OF'Q,C3'OF'Q,C4'OF'Q); 'REAL' T, S1:= 0, S2:= 0, MAX:= 0; 'FOR' I 'TO' N1 'DO' 'FOR'J 'TO' N2 'DO' T:= 'ABS' QK[I,J]; S1+:=T; S2+:=T*T; ( T>MAX ! MAX:=T ) 'OD' 'OD'; NRM[K,'AT'1]:= (S1,MAX, S1/(N1*N2), SQRT(S2/(N1*N2)) ) 'OD';NRM 'END'; #$E# #E$# #$H# 'OP' 'NORM1' = ('FIELD'F)[]'REAL': #$B NORM1 B$# 'BEGIN' [,]'REAL'N= 'NORM'F; N [ , (TYPE'OF'F=0!-1!1) ] 'END' ; #$E# #E$# #$H# 'OP' 'PRINTNORM' = ('FIELD'F)'VOID': #$B NRM.PRT B$# 'BEGIN' 'INT' TYPE=TYPE'OF'F; [,]'REAL' N = 'NORM' F; PRINT((NEWLINE,"INF-NORM: ",N[, 0], NEWLINE," 2-NORM: ",N[, 2], NEWLINE," 1-NORM: ",N[,(TYPE=0!-1!1)],NEWLINE)) 'END' ; #$E# #E$# #$H# 'OP' + = ('FIELD' A,B)'FIELD': #$B FLPLFL B$# 'BEGIN' 'COPY' A +:= B 'END' ; #$E# #E$# #$H# 'OP' - = ('FIELD' A,B)'FIELD': #$B FLMNFL B$# 'BEGIN' 'COPY' A -:= B 'END' ; #$E# #E$# #$H# 'OP' +:= = ('FIELD' A, 'FIELD' B)'FIELD': #$B PL.AFL B$# 'BEGIN' (TYPE'OF'A /= TYPE'OF'B ! ERROR ); 'REF'[,]'STATE' QA = Q'OF'A, QB = Q'OF'B; 'INT' L1=1'LWB'QB, U1=1'UPB'QB, L2=2'LWB'QB, U2=2'UPB'QB; 'FOR' I 'FROM' L1 'TO' U1 'DO' 'REF'[]'STATE'AI=QA[I,],BI=QB[I,]; 'FOR' J 'FROM' L2 'TO' U2 'DO' 'REF''STATE' AIJ=AI[J],BIJ=BI[J]; ( C1'OF'AIJ +:= C1'OF'BIJ , C2'OF'AIJ +:= C2'OF'BIJ , C3'OF'AIJ +:= C3'OF'BIJ , C4'OF'AIJ +:= C4'OF'BIJ ) 'OD' 'OD'; A 'END'; #$E# #E$# #$H# 'OP' -:= = ('FIELD' A, 'FIELD' B)'FIELD': #$B MN.AFL B$# 'BEGIN' (TYPE'OF'A /= TYPE'OF'B ! ERROR ); 'REF'[,]'STATE' QA = Q'OF'A, QB = Q'OF'B; 'INT' L1=1'LWB'QB, U1=1'UPB'QB, L2=2'LWB'QB, U2=2'UPB'QB; 'FOR' I 'FROM' L1 'TO' U1 'DO' 'REF'[]'STATE'AI=QA[I,],BI=QB[I,]; 'FOR' J 'FROM' L2 'TO' U2 'DO' 'REF''STATE' AIJ=AI[J],BIJ=BI[J]; ( C1'OF'AIJ -:= C1'OF'BIJ , C2'OF'AIJ -:= C2'OF'BIJ , C3'OF'AIJ -:= C3'OF'BIJ , C4'OF'AIJ -:= C4'OF'BIJ ) 'OD' 'OD'; A 'END'; #$E# #E$# #$H# 'OP' *:= = ('FIELD' A, 'REAL' B)'FIELD': #$B TM.AFL B$# 'IF' B = 1.0 'THEN' A 'ELSE' 'REF'[,]'STATE' QA = (Q'OF'A)[@1,@1]; 'INT' N1 = 1'UPB'QA, N2 = 2'UPB'QA; 'FOR' I 'TO' N1 'DO' 'REF'[]'STATE' AI=QA[I,]; 'FOR' J 'TO' N2 'DO' 'REF''STATE' AIJ =AI[J] ; ( C1'OF'AIJ*:=B,C2'OF'AIJ*:=B,C3'OF'AIJ*:=B,C4'OF'AIJ*:=B ) 'OD' 'OD' ; A 'FI' ; #$E# #E$# #$H# 'OP' * = ('NETMAT' JAC, 'FIELD' F) 'FIELD': #$B TIM B$# 'BEGIN'# PDZ # 'REF'[,]'STATE' QF = Q'OF'F; 'INT' L1 = 1'LWB'QF, U1 = 1'UPB'QF, L2 = 2'LWB'QF, U2 = 2'UPB'QF; (L1/='LWB'JAC 'OR' U1/='UPB'JAC ! ERROR); 'HEAP'[L1:U1,L2:U2]'STATE' QG; 'FIELD' G = (DDX'OF'F,DDY'OF'F, 'HEAP''INT':= 0, QG); 'FOR' I 'FROM' L1 'TO' U1 'DO' 'BOOL' LWB = (I=L1), UPB = (I=U1); 'REF'[,]'FLUXAC' JACI = GET(JAC[I]), 'REF' []'STATE' QFM = QF[(LWB!L1!I-1),], QFI = QF[ I ,], QFP = QF[(UPB!U1!I+1),], QGI = QG[ I ,]; 'INT' L3 = ('INT'L=2'LWB'JACI; (LWB!(L<-1!-1!L)!L)), U3 = ('INT'U=2'UPB'JACI; (UPB!(U> 1! 1!U)!U)); 'FOR' J 'FROM' L2 'TO' U2 'DO' 'BOOL' NOTL = (J>L2), NOTR = (J "YVELOCITY" ! ERROR; VARNO:= 0 ); ( NAME <= "YVELOCITY" ! VARNO:= 6); ( NAME <= "YMOMENT" ! VARNO:= 3); ( NAME <= "XVELOCITY" ! VARNO:= 5); ( NAME <= "XMOMENT" ! VARNO:= 2); ( NAME <= "TOTENTHALPY" ! VARNO:= 13); ( NAME <= "TOTENERGY" ! VARNO:= 4); ( NAME <= "TEMPERATURE" ! VARNO:= 14); ( NAME <= "SPEED" ! VARNO:= 10); ( NAME <= "SOUND" ! VARNO:= 7); ( NAME <= "PRESSURE" ! VARNO:= 9); ( NAME <= "LOCMACH" ! VARNO:= 11); ( NAME <= "ENTROPY" ! VARNO:= 8); ( NAME <= "DIRECTION" ! VARNO:= 12); ( NAME <= "DENSITY" ! VARNO:= 1); ( NAME = "" ! ERROR; VARNO:= 0); ( TYPE'OF'F /=0 ! STATE E(F) ); [,]'STATE' Q = (Q'OF'F)[@1,@1]; 'INT' L1Q = 1'LWB'Q, U1Q= 1'UPB'Q, L2Q = 2'LWB'Q, U2Q= 2'UPB'Q; 'HEAP'[L1Q:U1Q,L2Q:U2Q]'REAL' V; 'PROC'('STATE')'REAL' COMP = VAR[VARNO]; 'FOR' I 'FROM' L1Q 'TO' U1Q 'DO' 'FOR' J 'FROM' L2Q 'TO' U2Q 'DO' V[I,J]:= COMP(Q[I,J]) 'OD' 'OD'; 'VAR'(INTERN,DDX'OF'F, DDY'OF'F, VARNO, V) 'END'; #$E# #E$# # FLUX FUNCTION # 'PR' EJECT 'PR' #$H# 'PROC' ZFLUX = ('STATE' Q, 'REF''FLUXAC' FAC )'FLUX': #$B ZFLUX B$# 'BEGIN' 'REAL' U:=C1'OF'Q,V:=C2'OF'Q,C:=C3'OF'Q,Z:=C4'OF'Q, C2,U2,POR,R,RU,RUV,RUUP,W2,EPR,RUEPR,OC,OCU,OS; C2 := C*C; U2 := U*U; POR := C2*OGAM; R := EXP((LN(POR)-Z)*OGAMM1); RU := R*U; RUV := RU*V; RUUP := R*(U2+POR); W2 := (V*V+U2)*0.5; EPR := C2*OGAMM1 + W2; RUEPR:= RU*EPR; ( FAC 'ISNT''NIL' ! OS:= -OGAMM1; OC:= TOGAMM1*R/C; OCU:= OC*U; FAC := 'FLUXAC' ( R, 0, OCU , OS*RU , 2*RU, 0, OC * (U2+C2), OS*RUUP , R*V, RU, OCU* V , OS*RUV , (U2+EPR)*R,RUV, OCU*(EPR+C2), OS*RUEPR ) ); 'FLUX'( (RU, RUUP, RUV, RUEPR) ) 'END';#$E# #E$# # SIMULTANEOUS STATE FUNCTIONS # 'PR' EJECT 'PR' 'BOOL' WARNING:= 'TRUE'; 'PROC' REVEAL = ('REAL' C,Z, 'REF''REAL' P,R)'VOID': ('REAL'POR:= C*C/GAM; R:= EXP((LN(POR)-Z)*OGAMM1); P:= POR*R); #$H# 'OP' 'EZ' = ('STATE' A)'STATE': #$B EZ B$# ('REAL' R ,U:=C1'OF'A, V:= C2'OF'A, P, POR:= C3'OF'A*C3'OF'A/GAM; R:= EXP((LN(POR)-C4'OF'A)*OGAMM1); P:= POR*R; 'STATE' (R,R*U,R*V,(U*U+V*V)*R/2 + P/GAMM1) ); #$E# #E$# #$H# 'OP' 'ZE' = ('STATE' A)'STATE': #$B ZE B$# ('REAL' R:= C1'OF'A, M:= C2'OF'A, N:= C3'OF'A, P; P:= ((M*M+N*N)*-0.5/R + C4'OF'A)*GAMM1; 'STATE' (M/R,N/R,SQRT(GAM*P/R),LN(P)-GAM*LN(R)) ); #$E# #E$# #$H# 'OP' 'EP' = ([]'REAL' A)'STATE': #$B EP B$# ('REAL' R = A[1] ; 'REAL' RU = R*A[2] , RV = R*A[3] ; 'STATE'(R,RU,RV, (RU*RU+RV*RV)*0.5/R + A[4]*OGAMM1) ); #$E# #E$# #$H# 'OP' 'PE' = ('STATE' A)[]'REAL': #$B PE B$# ('REAL' R:= C1'OF'A; 'REAL' U = C2'OF'A/R, V = C3'OF'A/R; 'REAL' P:= ((U*U+V*V)*-0.5*R + C4'OF'A)*GAMM1; ( P<=0 'OR' R<=0 ! R:= P:= 1.0E-10; ( WARNING ! PRINT((NEWLINE,"NEGATIVE RHO OR P "+50*"!")); WARNING := 'FALSE' )); []'REAL'(R,U,V,P, SQRT(GAM*P/R), LN(P)-GAM*LN(R) ) ); #$E# #E$# #$H# 'OP' 'PZ' = ('STATE' A)[]'REAL': #$B PZ B$# ( 'REAL' P,R; REVEAL(C3'OF'A,C4'OF'A,P,R); []'REAL'(R,C1'OF'A,C2'OF'A,P,C3'OF'A,C4'OF'A) ); #$E# #E$# #$H# 'OP' 'ZP' = ([]'REAL' A)'STATE': #$B ZP B$# ( 'ZE''EP' A ); #$E# #E$# 'PR' EJECT 'PR' #$H# 'PROC' STATE Z = ('FIELD' F)'FIELD': #$B STT.Z B$# 'IF' TYPE'OF'F = 0 'THEN' ERROR; 'SKIP' 'ELIF' TYPE'OF'F = 1 'THEN' TYPE'OF'F:= 2; 'REF'[,]'STATE'Q = Q'OF'F; 'FOR' I 'FROM' 1'LWB'Q 'TO' 1'UPB'Q 'DO' 'REF'[]'STATE' QI=Q[I,]; 'FOR' J 'FROM' 2'LWB'Q 'TO' 2'UPB'Q 'DO' QI[J]:='ZE'QI[J] 'OD''OD'; F 'ELSE' F 'FI';#$E# #E$# #$H# 'PROC' STATE E = ('FIELD' F)'FIELD': #$B STT.E B$# 'IF' TYPE'OF'F = 0 'THEN' ERROR; 'SKIP' 'ELIF' TYPE'OF'F = 2 'THEN' TYPE'OF'F:= 1; 'REF'[,]'STATE'Q = Q'OF'F; 'FOR' I 'FROM' 1'LWB'Q 'TO' 1'UPB'Q 'DO' 'REF'[]'STATE'QI=Q[I,]; 'FOR' J 'FROM' 2'LWB'Q 'TO' 2'UPB'Q 'DO' QI[J]:='EZ'QI[J] 'OD''OD'; F 'ELSE' F 'FI';#$E# #E$# # NUMERICAL FLUX # 'PR' EJECT 'PR' 'PROC' NFLUX:= ('STATE' QO,Q1, 'REF''FLUXAC' FAC0,FAC1)'FLUX': ( ERROR; NUL ); 'INT' SOS:= 1; 'BOOL' CAVITATION:= 'FALSE'; #$H# 'PROC' OSHER = ('STATE' Q0,Q1,'REF' 'FLUXAC' FAC0,FAC1)'FLUX': #$B OSHER B$# 'BEGIN' 'REAL' U0=C1'OF'Q0,V0=C2'OF'Q0,C0=C3'OF'Q0,Z0=C4'OF'Q0, U1=C1'OF'Q1,V1=C2'OF'Q1,C1=C3'OF'Q1,Z1=C4'OF'Q1, SGN = SOS, SHGAMM1 = SOS*HGAMM1; 'BOOL' FACP = FAC1 'ISNT''NIL', FACM = FAC0 'ISNT''NIL'; 'REAL' PSI0 = U0-C0/SHGAMM1, PSI1 = U1+C1/SHGAMM1, ALFA = EXP((Z1-Z0)*OTGAM); 'REAL' OALFAP1=1/(1+ALFA); 'REAL' UH = (PSI1+ALFA*PSI0)*OALFAP1 +1.0E-200, CL = SHGAMM1*(PSI1-PSI0)*OALFAP1; (CL<0! # CAVITATION # ERROR ); 'REAL' CR = ALFA*CL, LAM0 = U0+SGN*C0, LAML = UH+SGN*CL, LAM1 = U1-SGN*C1, LAMR = UH-SGN*CR; 'FLUX' FLUC , F := NUL; 'REF''FLUXAC' FLAC = (FACM 'OR' FACP ! (FACP !FAC1:=NULFLUXAC); (FACM !FAC0:=NULFLUXAC); 'LOC''FLUXAC' ! 'NIL' ); (UH*LAM0*LAML*LAMR*LAM1=0 ! ERROR); (CL<0 ! CAVITATION := 'TRUE' ; EXCEPTION ) ; 'IF' LAM0>0 'THEN' (FACM ! F+:=ZFLUX(Q0, FAC0) ! F+:=ZFLUX(Q0,'NIL') ) 'FI'; 'IF' LAM1<0 'THEN' (FACP ! F+:=ZFLUX(Q1, FAC1) ! F+:=ZFLUX(Q1,'NIL') ) 'FI'; 'IF' LAM0*LAML < 0 'THEN' 'REAL' U= GAMM1OGAMP1*PSI0; FLUC:= ZFLUX((U,V0,-SGN*U,Z0),(FACM ! FLAC ! 'NIL')); (FACM!'FLUX' K = (1'FC'FLAC)-SGN*(3'FC'FLAC) ; CASS(FLAC,1,GAMM1OGAMP1 *K ) ; CASS(FLAC,3,(-SGN*TOGAMP1)*K ) ; (LAM0<0 ! FAC0+:=FLAC ! FAC0-:=FLAC )); (LAM0<0 ! F+:=FLUC ! F-:=FLUC ) 'FI'; 'IF' LAMR*LAM1 < 0 'THEN' 'REAL' U= GAMM1OGAMP1*PSI1; FLUC:= ZFLUX((U,V1,SGN*U,Z1),(FACP ! FLAC ! 'NIL' )); (FACP!'FLUX' K= 1'FC'FLAC+SGN*3'FC'FLAC; CASS(FLAC,1,GAMM1OGAMP1*K); CASS(FLAC,3,(SGN*TOGAMP1)*K); (LAM1>0 ! FAC1+:=FLAC ! FAC1-:=FLAC )); (LAM1>0 ! F+:=FLUC ! F-:=FLUC ) 'FI'; 'IF' LAML*UH<0 'THEN' FLUC:= ZFLUX((UH,V0,CL,Z0),FLAC); (FACM!'REAL' DUM= CR*OALFAP1*OGAMGAMM1; 'FLUX' K1= (ALFA*OALFAP1)*(1'FC'FLAC) -(SHGAMM1*OALFAP1)*(3'FC'FLAC), K2= DUM*(1'FC'FLAC+SHGAMM1*(3'FC'FLAC)); CSUB(FAC0,1,SGN*K1 ); CSUB(FAC0,2,SGN*(2'FC'FLAC) ); CSUB(FAC0,3,-TOGAMM1*K1 ); CSUB(FAC0,4,SGN*(4'FC'FLAC)+K2) ); (FACP!'REAL' DUM=-CR*OGAMGAMM1; 'FLUX' K = OALFAP1* (1'FC'FLAC+ SHGAMM1 *(3'FC'FLAC)); CSUB(FAC1,1, SGN *K); CSUB(FAC1,3, TOGAMM1*K); CSUB(FAC1,4, DUM *K) ); (SGN>0! F-:= FLUC ! F+:= FLUC) 'ELIF' UH*LAMR<0 'THEN' FLUC:= ZFLUX((UH,V1,CR,Z1),FLAC); (FACM!'REAL' DUM= CL*OGAMGAMM1; 'FLUX'K= ALFA*OALFAP1* (1'FC'FLAC- SHGAMM1 *(3'FC'FLAC)); CSUB(FAC0,1, SGN*K); CSUB(FAC0,3, -TOGAMM1*K); CSUB(FAC0,4, DUM*K) ); (FACP!'REAL' DUM=-CR*OALFAP1*OGAMGAMM1; 'FLUX' K1= OALFAP1* (1'FC'FLAC+(SHGAMM1*ALFA)*3'FC'FLAC), K2= DUM*(1'FC'FLAC- SHGAMM1 *3'FC'FLAC); CSUB(FAC1,1, SGN*K1); CSUB(FAC1,2, SGN*2'FC'FLAC); CSUB(FAC1,3, TOGAMM1*K1); CSUB(FAC1,4, SGN*4'FC'FLAC+K2) ); (SGN>0! F-:= FLUC ! F+:= FLUC) 'FI'; EXCEPTION : 'SKIP' ; 'FLUX'(F) 'END';#$E# #E$# 'PR' EJECT 'PR' #$H# 'PROC' POSINT = ('STATE' Q0,Q1, 'REF''FLUXAC' FAC0,FAC1)'FLUX': #$B POSINT B$# 'BEGIN' 'FLUX' F:= OSHER(Q0,Q1,FAC0,FAC1); (FAC1 'IS' 'NIL' ! F:=ZFLUX(Q1,'NIL')-F ! 'FLUXAC' FLAC; F:=ZFLUX(Q1,FLAC)-F; FAC1:=FLAC-FAC1 ); (FAC0 'ISNT''NIL' ! FAC0:= -FAC0); 'FLUX'(F) 'END'; #$E# #E$# #$H# 'PROC' NEGINT = ('STATE' Q0,Q1, 'REF' 'FLUXAC' FAC0,FAC1)'FLUX': #$B NEGINT B$# 'BEGIN' 'FLUX' F:= OSHER(Q0,Q1,FAC0,FAC1); (FAC0 'IS' 'NIL' ! F-:= ZFLUX(Q0,'NIL') ! 'FLUXAC' FLAC; F-:= ZFLUX(Q0, FLAC); FAC0-:=FLAC); 'FLUX'(F) 'END'; #$E# #E$# # BOUNDARY STATES # 'PR' EJECT 'PR' 'PROC'('WALL','REAL','REAL','STATE','REF''STATAC')'STATE' NOORD, ZUID, OOST, WEST; #$H# 'PROC' FIXED UVCZ = ('WALL' W, 'REAL' U,V,C,Z,'STATE' QIN, 'REF''STATAC' QAC)'STATE': #$B BC.FIX B$# ( (QAC 'ISNT''NIL' ! QAC:=NULFLUXAC); 'STATE'(U,V,C,Z) ); #$E# #E$# #$H# 'PROC' SUB OUT PRESS = ('WALL' W, 'BOOL' LEFT, 'REAL' P1, 'STATE' QIN, 'REF''STATAC' QAC)'STATE': #$B BC.RSO B$# 'BEGIN' # NOT CHECKED: LEFT='TRUE' # 'REAL' C =C'OF'W, S =S'OF'W, C0=C3'OF'QIN, Z=C4'OF'QIN; 'REAL' U0=C*C1'OF'QIN+S*C2'OF'QIN,V=-S*C1'OF'QIN+C*C2'OF'QIN; 'REAL' R1,C1,U1, STOGM:= ( LEFT ! -TOGAMM1 ! TOGAMM1 ); R1:=EXP((LN(P1)-Z)*OGAM); C1:=SQRT(GAM*P1/R1); U1:=U0+STOGM*(C0-C1); 'IF' QAC 'ISNT''NIL' 'THEN' 'REAL' DUM=C1*OGAM; QAC:=( 1 ,0 ,C*STOGM ,-C*OGAMM1*DUM , 0 ,1 ,S*STOGM ,-S*OGAMM1*DUM , 0 ,0 ,0 , 0.5*DUM , 0 ,0 ,0 , 1 ) 'FI'; 'STATE'(C*U1-S*V,S*U1+C*V,C1,Z) 'END'; #$E# #E$# #$H# 'PROC' WALL = ('WALL' W, 'BOOL' LEFT, 'REAL' NOR,'STATE' QIN, 'REF''STATAC' QAC)'STATE': #$B BC.LW B$# 'BEGIN' 'REAL' C = C'OF'W, S= S'OF'W, C1= C3'OF'QIN, Z = C4'OF'QIN; 'REAL' U1 =C*C1'OF'QIN+S*C2'OF'QIN,V= -S*C1'OF'QIN+C*C2'OF'QIN, SHG= (LEFT!-HGAMM1!+HGAMM1); 'REAL' U0:=NOR*V,C0; C0:=C1-SHG*(U0-U1); 'IF' QAC 'ISNT''NIL' 'THEN' QAC:=( S*( -C*NOR+S) , C*( C*NOR-S ) ,0 ,0 , S*( -S*NOR-C) , C*( S*NOR+C ) ,0 ,0 , (C+S*NOR)*SHG , (S-C*NOR)*SHG ,1 ,0 , 0 , 0 ,0 ,1 ) 'FI'; 'STATE'(C*U0-S*V,S*U0+C*V,C0,Z) 'END'; #$E# #E$# 'PR' EJECT 'PR' #$H# 'PROC' SUB IN UVZ = ('WALL' W,'BOOL' LEFT,'REAL' UX,VY,Z0,'STATE' QIN, 'REF''STATAC' QAC)'STATE': #$B BC.LSI B$# 'BEGIN' # NOT CHECKED: LEFT='FALSE' # 'REAL' C =C'OF'W, S =S'OF'W, C1=C3'OF'QIN, Z1=C4'OF'QIN; 'REAL' U0=C*UX +S*VY , V0=-S*UX +C*VY , U1=C*C1'OF'QIN+S*C2'OF'QIN,V=-S*C1'OF'QIN+C*C2'OF'QIN, SHG= (LEFT !HGAMM1 !-HGAMM1); 'REAL' C0, DD, D1, D2; ( SOS>0 #I.E. OSHERS SCHEME # !'REAL' P1,R1,RH,CH; REVEAL(C1,Z1,P1,R1); RH:=EXP(OGAM*(LN(P1)-Z0)); CH:=SQRT(GAM*P1/RH); C0:=CH+SHG*(U0-U1); DD:= R1*C1/(RH*CH); D1:= -SHG; D2:= -CH*OTGAM !'REAL' PH,RH,CH,R0; CH:= C1+SHG*(U0-U1); REVEAL(CH,Z1,PH,RH); R0:= EXP(OGAM*(LN(PH)-Z0)); C0:= SQRT(GAM*PH/R0); DD:= RH*CH/(R0*C0); D1:= -SHG*DD; D2:= -C0*OTGAM ); ( QAC 'ISNT''NIL' ! QAC:= ( 0 , 0 , 0 ,0 , 0 , 0 , 0 ,0 , D1*C, D1*S, DD ,D2 , 0 , 0 , 0 ,0 ) );'STATE'(UX,VY,C0,Z0) 'END'; #$E# #E$# # EXTEND FIELD WITH WALL STATES # 'PR' EJECT 'PR' #$H# 'PROC' ADD WALL STATES = ('BOOL' SECOND, 'FIELD' QF)'FIELD': #$B ADDW B$# 'BEGIN' 'REF'[,]'STATE' Q = Q'OF'QF; 'REAL' DDX = DDX'OF'QF, DDY = DDY'OF'QF; 'INT' IL=1'LWB'Q,IU=1'UPB'Q,JL=2'LWB'Q,JU=2'UPB'Q; STATE Z ( QF ); [IL-1:IU+1,JL-1:JU+1]'STATE' S; S[IL:IU,JL:JU]:= Q[IL:IU,JL:JU]; 'STATE' QQ; 'FIELD' NEW = 'FIELD'((DDX'OF'QF, DDY'OF'QF, 'HEAP''INT':=TYPE'OF'QF, S)); 'BOOL' PLUS = 'TRUE'; 'REAL' XX:= (IL-1.5)*DDX; [JL-1:JU] 'WALL' HORI,VERI ; [JL-1:JU] 'POINT' COORD ; ( 'REF'[]'STATE' QI = Q[IL,], QIP = Q[IL+1,]; WALLS('TRUE',PLUS,(IL-1)*DDX,DDY,COORD,'NIL',HORI); 'REAL' YY:= (JL-1.5)*DDY; 'FOR' J 'FROM' JL 'TO' JU 'DO' QQ:= QI[J]; (SECOND! QQ+:= (QQ-QIP[J])/2.0); S[IL-1,J]:=NOORD(HORI[J],XXN,YY+:=DDY,QQ,'NIL') 'OD' ); 'FOR' I 'FROM' IL 'TO' IU 'DO' 'REF'[]'STATE' QI = Q[I,]; WALLS('TRUE',PLUS,I*DDX,DDY,COORD, VERI,HORI); QQ:= QI[JL]; (SECOND! QQ+:= (QQ-QI[JL+1])/2.0); S[I,JL-1]:=WEST(VERI[JL-1],XX+:=DDX,YYW,QQ,'NIL'); QQ:= QI[JU]; (SECOND! QQ+:= (QQ-QI[JU-1])/2.0); S[I,JU+1]:=OOST(VERI[JU ],XX ,YYO,QQ,'NIL') 'OD'; ( 'REF'[]'STATE' QI = Q[IU,], QIM = Q[IU-1,]; 'REAL' YY:= (JL-1.5)*DDY; 'FOR' J 'FROM' JL 'TO' JU 'DO' QQ:= QI[J]; (SECOND! QQ+:= (QQ-QIM[J])/2.0); S[IU+1,J]:= ZUID(HORI[J],XXZ,YY+:=DDY,QQ,'NIL') 'OD' ); S[IL-1,JL-1] := S[IL-1,JL] + S[IL,JL-1] - S[IL,JL]; S[IL-1,JU+1] := S[IL-1,JU] + S[IL,JU+1] - S[IL,JU]; S[IU+1,JL-1] := S[IU+1,JL] + S[IU,JL-1] - S[IU,JL]; S[IU+1,JU+1] := S[IU+1,JU] + S[IU,JU+1] - S[IU,JU]; NEW 'END'; #$E# #E$# 'PR' EJECT 'PR' #$H# 'PROC' ADD WALL = ('BOOL' SECOND, 'FIELD' F)'FIELD': #$B ADDWSH B$# 'BEGIN' 'REF'[,]'STATE' Q = Q'OF'F; 'INT' L1 = 1'LWB' Q, U1 = 1'UPB'Q, L2 = 2'LWB'Q, U2 = 2'UPB'Q; 'INT' L1S = L1-1, U1S = U1+1, L2S = L2-1, U2S = U2+1; 'HEAP'[L1S:U1S,L2S:U2S]'STATE' S; S[L1:U1,L2:U2]:= Q[L1:U1,L2:U2]; 'IF' SECOND 'AND' (U1-L1+1>=3) 'AND' (U2-L2+1>=3) 'THEN' 'FOR' J 'FROM' L2S+1 'TO' U2S-1 'DO' S[L1S,J]:= 3.0*Q[L1,J]-3.0*Q[L1+1,J]+Q[L1+2,J]; S[U1S,J]:= 3.0*Q[U1,J]-3.0*Q[U1-1,J]+Q[U1-2,J] 'OD' ; 'FOR' I 'FROM' L1S+1 'TO' U1S-1 'DO' S[I,L2S]:= 3.0*Q[I,L2]-3.0*Q[I,L2+1]+Q[I,L2+2]; S[I,U2S]:= 3.0*Q[I,U2]-3.0*Q[I,U2-1]+Q[I,U2-2] 'OD' 'ELSE' 'FOR' J 'FROM' L2S+1 'TO' U2S-1 'DO' S[L1S,J]:= 2.0*Q[L1,J]-Q[L1+1,J]; S[U1S,J]:= 2.0*Q[U1,J]-Q[U1-1,J] 'OD' ; 'FOR' I 'FROM' L1S+1 'TO' U1S-1 'DO' S[I,L2S]:= 2.0*Q[I,L2]-Q[I,L2+1]; S[I,U2S]:= 2.0*Q[I,U2]-Q[I,U2-1] 'OD' 'FI' ; S[L1S,L2S]:= S[L1S ,L2S+1] + S[L1S+1,L2S ] - S[L1S+1,L2S+1]; S[L1S,U2S]:= S[L1S+1,U2S ] + S[L1S ,U2S-1] - S[L1S+1,U2S-1]; S[U1S,L2S]:= S[U1S ,L2S+1] + S[U1S-1,L2S ] - S[U1S-1,L2S+1]; S[U1S,U2S]:= S[U1S-1,U2S ] + S[U1S ,U2S-1] - S[U1S-1,U2S-1]; 'FIELD'(DDX'OF'F,DDY'OF'F,'HEAP''INT':=TYPE'OF'F,S) 'END' ; #$E# #E$# # WALLS # 'PR' EJECT 'PR' 'PROC' WALLS := ('BOOL' FORWARD,PLUS, 'REAL' XX, DDY, 'REF'[]'POINT' COORD, 'REF'[]'WALL' VERI,HORI)'VOID': (ERROR;'SKIP'); #$H# 'PROC' WALLSG = ('BOOL' FORWARD,PLUS, 'REAL' XX, DDY, 'REF'[]'POINT' COORD, 'REF'[]'WALL' VERI,HORI)'VOID': #$B WALLSG B$# 'BEGIN' 'BOOL' VER = VERI 'ISNT' 'NIL'; 'INT' JLM = 'LWB' HORI, JU = 'UPB' HORI; 'REAL' DX,DY,DS; 'POINT' PT,PB,PO; 'REF''REAL' XT= X'OF'PT, XB= X'OF'PB, XO= X'OF'PO, YT= Y'OF'PT, YB= Y'OF'PB, YO= Y'OF'PO; # ALLEEN AFBEELDINGEN MET POSITIEVE ORIENTATIE ZYN TOEGESTAAN # 'FOR' J 'FROM' JLM 'TO' JU 'DO' 'IF' VER 'THEN' PT:= COORD[J]; COORD[J]:= PB:= MAPPING(XX,J*DDY); (FORWARD ! DX:= XT-XB; DY:= YT-YB ! DX:= XB-XT; DY:= YB-YT ); DS:= SQRT(DX*DX+DY*DY); VERI[J]:= 'WALL'(DY/DS,-DX/DS,(PLUS!DS!-DS)) 'ELSE' COORD[J]:= PB:= MAPPING(XX,J*DDY) 'FI'; 'IF' J > JLM 'THEN' DX:= XB-XO; DY:= YB-YO; DS:= SQRT(DX*DX+DY*DY); HORI[J]:= 'WALL'( DY/DS,-DX/DS,(PLUS!DS!-DS)) 'FI'; PO:= PB 'OD' 'END'; #$E# #E$# 'PR' EJECT 'PR' #$H# 'PROC' WALLSR = ('BOOL' FORWARD,PLUS, 'REAL' XX, DDY, 'REF'[]'POINT' COORD, 'REF'[]'WALL' VERT,HORI)'VOID': #$B WALLSR B$# 'IF' HORI 'IS' 'NIL' 'THEN' 'FOR' J 'FROM' 'LWB' COORD 'TO' 'UPB' COORD 'DO' COORD[J]:= MAPPING(XX,J*DDY) 'OD' 'ELSE' 'BOOL' VER = VERT 'ISNT' 'NIL'; 'INT' JLM = 'LWB' HORI, JU = 'UPB' HORI; 'POINT' PB; 'REF''REAL' XB = X'OF'PB, YB = Y'OF'PB; 'REAL' XT,YO, DX,DY,DS; 'WALL' VERTI; XT:= X'OF'COORD[JLM]; COORD[JLM]:= PB:= MAPPING(XX, JLM*DDY); 'IF' VER 'THEN' DX:= (FORWARD! XT-XB ! XB-XT); DS:= 'ABS'DX; VERT[JLM]:= VERTI:= 'WALL'(0,-'SIGN'DX,(PLUS!DS!-DS)) 'FI'; 'FOR' J 'FROM' JLM+1 'TO' JU 'DO' YO:= YB; ( VER ! VERT[J]:= VERTI); COORD[J]:= PB:= MAPPING(XX,J*DDY); DY:= YB-YO; DS:='ABS'DY; HORI [J]:= 'WALL'('SIGN'DY, 0,(PLUS!DS!-DS)) 'OD' 'FI'; #$E# #E$# # FLUXT # 'PR' EJECT 'PR' 'PROC' FLUXT:= ('WALL' W, 'STATE' Q0,Q1,'REF''FLUXAC' NF0,NF1) 'FLUX': (ERROR;'SKIP'); #$H# 'PROC' FLUXTR = ('WALL' W, 'STATE' Q0,Q1,'REF''FLUXAC' NF0,NF1) 'FLUX': #$B FLUXTR B$# 'IF' 'REAL' DS = DS'OF'W, ADS = 'ABS'DS'OF'W; C'OF'W > 0.999 'THEN' 'FLUX' N := NFLUX(Q0,Q1,NF0,NF1); (NF0'ISNT''NIL'! NF0 *:= ADS ); (NF1'ISNT''NIL'! NF1 *:= ADS ); N*:=DS 'ELIF' S'OF'W > 0.999 'THEN' 'REAL' U0=C1'OF'Q0, V0=C2'OF'Q0, U1=C1'OF'Q1, V1=C2'OF'Q1; 'FLUX' NF = NFLUX ( (V0, -U0, C3'OF'Q0, C4'OF'Q0), (V1, -U1, C3'OF'Q1, C4'OF'Q1), NF0, NF1); 'FLUX' NT; ( NF0 'ISNT' 'NIL' ! NT:= -(2'FC'NF0) ; CASS(NF0,2,1'FC'NF0) ; CASS(NF0,1,NT); NT:= -(3'FR'NF0) ; RASS(NF0,3,2'FR'NF0) ; RASS(NF0,2,NT); NF0 *:= ADS ); ( NF1 'ISNT' 'NIL' ! NT:= -(2'FC'NF1) ; CASS(NF1,2,1'FC'NF1) ; CASS(NF1,1,NT); NT:= -(3'FR'NF1) ; RASS(NF1,3,2'FR'NF1) ; RASS(NF1,2,NT); NF1 *:= ADS ); 'FLUX' ( DS*C1'OF'NF,-DS*C3'OF'NF, DS*C2'OF'NF, DS*C4'OF'NF ) 'ELSE' ERROR; 'SKIP' 'FI' ; #$E# #E$# 'PR' EJECT 'PR' #$H# 'PROC' FLUXTG = ('WALL' W, 'STATE' Q0,Q1,'REF''FLUXAC' NF0,NF1) 'FLUX': #$B FLUXTG B$# 'BEGIN' 'REAL' C = C'OF'W, S = S'OF'W, DS = DS'OF'W, ADS = 'ABS'DS'OF'W, U0=C1'OF'Q0, V0=C2'OF'Q0, U1=C1'OF'Q1, V1=C2'OF'Q1; 'FLUX' NF = NFLUX ( (C*U0+S*V0, -S*U0+C*V0, C3'OF'Q0, C4'OF'Q0), (C*U1+S*V1, -S*U1+C*V1, C3'OF'Q1, C4'OF'Q1),NF0,NF1); 'STATE' MT,NT; ( NF0 'ISNT' 'NIL' ! MT := 1'FC'NF0; NT := 2'FC'NF0; CASS(NF0,1,C*MT-S*NT); CASS(NF0,2,S*MT+C*NT); MT := 2'FR'NF0; NT := 3'FR'NF0; RASS(NF0,2,C*MT-S*NT); RASS(NF0,3,S*MT+C*NT); NF0 *:= ADS ); ( NF1 'ISNT' 'NIL' ! MT := 1'FC'NF1; NT := 2'FC'NF1; CASS(NF1,1,C*MT-S*NT);CASS(NF1,2,S*MT+C*NT); MT := 2'FR'NF1; NT := 3'FR'NF1; RASS(NF1,2,C*MT-S*NT);RASS(NF1,3,S*MT+C*NT); NF1 *:= ADS ); 'REAL' MS = DS*C2'OF'NF, NS = DS*C3'OF'NF; 'FLUX' ( DS*C1'OF'NF, C*MS-S*NS, S*MS+C*NS, DS*C4'OF'NF ) 'END'; #$E# #E$# # RESIDUAL # 'PR' EJECT 'PR' 'BOOL' CENTRAL2:= 'FALSE', ORDER2:= 'FALSE'; #$H# 'OP' 'MINRESIDU' = ('FIELD' F,Q)'FIELD': #$B MINRES B$# 'BEGIN' 'FIELD' FF = 'COPY'F ; RESIDU('FALSE',Q,FF); FF 'END' ; #$E# #E$# #$H# 'OP' 'PLUSRESIDU' = ('FIELD' F,Q)'FIELD': #$B PLUSRES B$# 'BEGIN' 'FIELD' FF = 'COPY'F ; RESIDU('TRUE',Q,FF); FF 'END' ; #$E# #E$# #$H# 'OP' 'RESIDU' = ('FIELD' Q)'FIELD': #$B RESIDU B$# 'BEGIN' 'FIELD' F = 'NULRHS'Q ; RESIDU('TRUE',Q,F); F 'END' ; #$E# #E$# 'PROC' RESIDU :=('BOOL'PLUS,'FIELD'QF,FF)'FIELD': 'SKIP'; 'PROC' RESIDUHO:=('BOOL'PLUS,'FIELD'QF,FF)'FIELD':RESIDU(PLUS,QF,FF); #$H# 'PROC' RESIDU1 = ('BOOL'PLUS,'FIELD' QF,RF)'FIELD': #$B RESID1 B$# 'BEGIN' 'REF'[,]'STATE' Q = Q'OF'QF, RHS = Q'OF'RF; 'REAL' DDX = DDX'OF'QF, DDY = DDY'OF'QF; 'INT' IL=1'LWB'Q,IU=1'UPB'Q,JL=2'LWB'Q,JU=2'UPB'Q; 'REAL' XX:= (IL-1.5)*DDX; [JL-1:JU] 'WALL' HORI,VERI; [JL-1:JU]'POINT' COORD; 'STATE' QOUT; STATE Z ( QF ); 'FLUX' NF; ( 'REF'[] 'STATE' QI = Q[IL,], RHSI = RHS[IL,]; WALLS('TRUE',PLUS,(IL-1)*DDX,DDY,COORD,'NIL',HORI); 'REAL' YY:= (JL-1.5)*DDY; 'FOR' J 'FROM' JL 'TO' JU 'DO' QOUT := NOORD(HORI[J],XXN,YY+:=DDY,QI[J],'NIL'); NF := FLUXT(HORI[J],QOUT,QI[J],'NIL','NIL'); RHSI[J]-:= NF 'OD' ); 'FOR' I 'FROM' IL 'TO' IU 'DO' 'REF'[]'STATE' QI = Q[I,], RHSI = RHS[I,]; WALLS('TRUE',PLUS,I*DDX,DDY,COORD, VERI,HORI); QOUT := WEST(VERI[JL-1],XX+:=DDX,YYW,QI[JL],'NIL'); NF := FLUXT(VERI[JL-1],QOUT,QI[JL],'NIL','NIL'); RHSI[JL]-:= NF; 'FOR' J 'FROM' JL 'TO' JU-1 'DO' # VERT. WALLS # NF:= FLUXT(VERI[J],QI[J],QI[J+1],'NIL','NIL'); RHSI[J ]+:= NF; RHSI[J+1]-:= NF 'OD'; QOUT := OOST(VERI[JU], XX,YYO,QI[JU],'NIL'); NF := FLUXT(VERI[JU],QI[JU],QOUT,'NIL','NIL'); RHSI [JU]+:= NF; 'IF' I < IU 'THEN' 'REF'[]'STATE' QIP = Q[I+1,], RHSIP = RHS[I+1,]; 'FOR' J 'FROM' JL 'TO' JU 'DO' #HORIZ. WALLS # NF:= FLUXT(HORI[J],QI[J],QIP[J],'NIL','NIL'); RHSI [J ]+:= NF; RHSIP[J ]-:= NF 'OD' 'FI' 'OD'; ( 'REF'[] 'STATE' QI = Q[IU,], RHSI = RHS[IU,]; 'REAL' YY:= (JL-1.5)*DDY; 'FOR' J 'FROM' JL 'TO' JU 'DO' QOUT := ZUID(HORI[J],XXZ,YY+:=DDY,QI[J],'NIL'); NF := FLUXT(HORI[J],QI[J],QOUT,'NIL','NIL'); RHSI[J]+:= NF 'OD' );'SKIP' 'END'; #$E# #E$# # SECOND ORDER RESIDUAL (I) # 'PR' EJECT 'PR' 'PROC' RESIDU2:= ('BOOL'PLUS,'FIELD' QF,RF)'FIELD': ( ERROR; 'SKIP'); #$H# 'PROC' CENTRALX= ('BOOL' P, 'FIELD' Q,F) 'FIELD': #$B CNTR.X B$# 'BEGIN' NFLUX:= ('STATE' Q0,Q1, 'REF''FLUXAC' D0,D1)'FLUX': ZFLUX(0.5*(Q0+Q1), 'NIL'); RESIDU1 (P,Q,F); NFLUX:= OSHER; 'SKIP' 'END'; #$E# #E$# #$H# 'PROC' CENTRALY= ('BOOL' P, 'FIELD' Q,F) 'FIELD': #$B CNTR.Y B$# 'BEGIN' NFLUX:= ('STATE' Q0,Q1, 'REF''FLUXAC' D0,D1)'FLUX': 0.5*(ZFLUX(Q0,'NIL') + ZFLUX(Q1,'NIL')); RESIDU1 (P,Q,F); NFLUX:= OSHER; 'SKIP' 'END'; #$E# #E$# # RESIDUAL SMOOTHERS # 'PR' EJECT 'PR' #$H# 'PROC' SMOOTHBOX = ('FIELD' A)'FIELD': #$B SMOBX B$# 'BEGIN' 'INT' TYPE = TYPE'OF'A; TYPE'OF'A:= 1; 'FIELD' B = 'PROLON' 'RESTRICT'A; TYPE'OF'A:= TYPE; Q'OF'A := Q'OF'B; A 'END'; #$E# #E$# #$H# 'PROC' SMOOTH5 = ('FIELD' A)'FIELD': #$B SMO5 B$# 'BEGIN' 'REF'[,]'STATE' Q=Q'OF'A; 'INT' IL=1'LWB'Q,IU=1'UPB'Q,JL=2'LWB'Q,JU=2'UPB'Q; [JL:JU]'STATE'QIM:=Q[IL,JL:JU] ; 'STATE' QIJM,QIJOLD ; 'FOR' I 'FROM' IL 'TO' IU-1 'DO' 'REF'[]'STATE'QI=Q[I,],QIP=Q[I+1,] ; QIJM:=QI[JL] ; 'FOR' J 'FROM' JL 'TO' JU-1 'DO' 'REF''STATE' QIJ=QI[J]; QIJOLD:=QIJ ; (((((QIJ*:=4.0)+:=QIJM)+:=QIM[J])+:=QI[J+1])+:=QIP[J])*:=0.125; QIJM := QIJOLD ; QIM[J]:=QIJOLD 'OD' ; QIJOLD:=QI[JU] ; ((((QI[JU]*:=5.0)+:=QIJM)+:=QIM[JU])+:=QIP[JU])*:=0.125 ; QIM[JU]:=QIJOLD 'OD' ; 'REF'[]'STATE' QI=Q[IU,] ; QIJM:=QI[JL] ; 'FOR' J 'FROM' JL 'TO' JU-1 'DO' 'REF''STATE' QIJ=QI[J] ; QIJOLD := QIJ ; ((((QIJ*:=5.0)+:=QIJM)+:=QIM[J])+:= QI[J+1])*:=0.125 ; QIJM := QIJOLD 'OD' ; ((( QI[JU]*:=6.0)+:=QIJM)+:=QIM[JU])*:=0.125 ; A 'END'; #$E# #E$# #$H# 'PROC' JACSMOOTH = ('REAL' DAMPING, 'FIELD' QF, RHS)'FIELD': #$B JACSMO B$# 'BEGIN' 'FIELD' JJJ = 'COPY' QF; 'REF'[,]'STATE' Q = Q'OF'QF, FF = Q'OF'RHS, JJ= Q'OF'JJJ; STATE Z ( QF ); 'REAL' DDX = DDX'OF'QF, DDY = DDY'OF'QF; 'INT' IL=1'LWB'FF,IU=1'UPB'FF,JL=2'LWB'FF,JU=2'UPB'FF; [JL-1:JU] 'WALL' HORIN,HORIZ,VERI ; [JL-1:JU] 'POINT' COORD ; [JL :JU] 'STATE' QIM,QIO ; 'STATE' UC, QOUT; 'STATAC' QAC ; 'FLUX' NF,F ; 'FLUXAC' NFAC0,NFAC1,FAC ; WALLS( 'TRUE','TRUE',(IL-1)*DDX,DDY,COORD,'NIL',HORIN); 'FOR' I 'FROM' IL 'TO' IU 'DO''REAL'XX=(I-0.5)*DDX; WALLS( 'TRUE','TRUE', I *DDX,DDY,COORD, VERI,HORIZ); QIO:= Q[I,]; 'FOR' J 'FROM' JL 'TO' JU 'DO''REAL'YY=(J-0.5)*DDY; 'REF''STATE' QIJ=Q[I,J]; UC:=QIJ; F :=(I=IU ! QOUT := ZUID(HORIZ[J],XXZ,YY,UC,QAC); FLUXT(HORIZ[J],UC,QOUT,NFAC0,NFAC1) ! FLUXT(HORIZ[J],UC,Q[I+1,J],NFAC0,'NIL')); FAC :=(I=IU ! NFAC0+:= NFAC1*QAC ! NFAC0 ); F -:=(I=IL ! QOUT := NOORD(HORIN[J],XXN,YY,UC,QAC); FLUXT(HORIN[J],QOUT,UC,NFAC0,NFAC1) ! FLUXT(HORIN[J],QIM[J],UC,'NIL',NFAC1)); FAC-:=(I=IL ! NFAC1+:= NFAC0*QAC ! NFAC1 ); F +:=(J=JU ! QOUT := OOST(VERI[J],XX,YYO,UC,QAC); FLUXT(VERI[J],UC,QOUT,NFAC0,NFAC1) ! FLUXT(VERI[J],UC,Q[I,J+1],NFAC0,'NIL')); FAC+:=(J=JU ! NFAC0+:= NFAC1*QAC ! NFAC0 ); F -:=(J=JL ! QOUT := WEST(VERI[J-1],XX,YYW,UC,QAC); FLUXT(VERI[J-1],QOUT,UC,NFAC0,NFAC1) ! FLUXT(VERI[J-1],QIO[J-1],UC,'NIL',NFAC1)); FAC-:=(J=JL ! NFAC1+:=NFAC0*QAC ! NFAC1 ); JJ[I,J] := (FF[I,J]-F) /FAC; JJ[I,J]*:= DAMPING 'OD'; QIM:= QIO; HORIN:= HORIZ 'OD'; JJJ 'END'; #$E# #E$# # SECOND ORDER RESIDUAL (II) # 'PR' EJECT 'PR' #$H# 'PROC' SUPERBOX = ('BOOL'PLUS,'FIELD' QF,RF)'FIELD': #$B SUPERB B$# 'IF' 'REF'[,]'STATE' QQ = Q'OF'QF, RHS = Q'OF'RF; 'REAL' DDX = DDX'OF'QF, DDY = DDY'OF'QF; 'INT' IL=1'LWB'QQ,IU=1'UPB'QQ,JL=2'LWB'QQ,JU=2'UPB'QQ; 'NOT'( 'ODD'(IU-IL) 'AND' 'ODD'(JU-JL) ) 'THEN' RESIDU1(PLUS,QF,RF) 'ELSE' [0:2,JL-1:JU] 'WALL' HOR; [1:2,JL-1:JU] 'WALL' VER; [JL-1:JU]'POINT' COORD ; [JL :JU]'STATE' OUTW ; [1:2] 'STATE' LEFT ; [1:2,1:2]'STATE' NS,WE ; 'STATE' QOUT ; STATE Z ( QF ); 'FLUX' NF; 'INT' II,JJ,JM; 'STATE' A,B,C,D,P,Q,R,S; 'REAL' XX:= (IL-1.5)*DDX, YY:= (JL-1.5)*DDY; WALLS('TRUE',PLUS,(IL-1)*DDX,DDY,COORD,'NIL',HOR[0,]); 'FOR' I 'FROM' IL 'BY' 2 'TO' IU 'DO' WALLS('TRUE',PLUS, I *DDX,DDY,COORD,VER[1,],HOR[1,]); WALLS('TRUE',PLUS,(I+1)*DDX,DDY,COORD,VER[2,],HOR[2,]); 'FOR' J 'FROM' JL 'BY' 2 'TO' JU 'DO' JM:= J-1; # INTERPOLATIE GEDEELTE # P:= QQ[ I ,J ] ; R:= QQ[ I ,J+1] ; Q:= QQ[ I+1,J ] ; S:= QQ[ I+1,J+1] ; A:= 0.25*(P+Q+R+S); C:= 0.25*(R+S-P-Q); B:= 0.25*(Q+S-P-R); D:= 0.25*(P+S-R-Q); # Q = A + B.X + C.Y + D.XY # 'FOR' K 'TO' 2 # K=1: N,W; K=2: S,E # 'DO' 'REAL' Z = ( K ! -2.0 , 2.0 ); P:=A+Z*B ; Q:=C+Z*D ; NS[K,1]:=P-Q ; NS[K,2]:=P+Q ; P:=A+Z*C ; Q:=B+Z*D ; WE[K,1]:=P-Q ; WE[K,2]:=P+Q 'OD'; # NAGEGAAN 850315 # 'FOR' L 'TO' 2 'DO' JJ:= J+L-1; ( I=IL !OUTW[JJ] :=NOORD(HOR[0,JJ],XXN,YY+:=DDY,NS[1,L] ,'NIL'); NF :=FLUXT(HOR[0,JJ],OUTW[JJ],NS[1,L],'NIL','NIL'); RHS[I ,JJ] -:= NF !NF :=FLUXT(HOR[0,JJ],OUTW[JJ],NS[1,L],'NIL','NIL'); RHS[I-1,JJ] +:= NF; RHS[I ,JJ] -:= NF ); OUTW[JJ]:= NS[2,L] 'OD'; 'FOR' K 'TO' 2 'DO' II:= I+K-1; (J=JL !LEFT[K] :=WEST (VER[K,JM],XX+:=DDX,YYW,WE[1,K] ,'NIL'); NF :=FLUXT(VER[K,JM],LEFT[K],WE[1,K],'NIL','NIL'); RHS[II,J ] -:= NF !NF :=FLUXT(VER[K,JM],LEFT[K],WE[1,K],'NIL','NIL'); RHS[II,JM] +:= NF; RHS[II,J ] -:= NF ); LEFT[K]:= WE[2,K] 'OD'; NF:= FLUXT( VER[1,J ],A-B,A-B,'NIL','NIL' ); RHS[I ,J ]+:= NF; RHS[I ,J+1] -:= NF; NF:= FLUXT( VER[2,J ],A+B,A+B,'NIL','NIL' ); RHS[I+1,J ]+:= NF; RHS[I+1,J+1] -:= NF; NF:= FLUXT( HOR[1,J ],A-C,A-C,'NIL','NIL' ); RHS[I ,J ]+:= NF; RHS[I+1,J ] -:= NF; NF:= FLUXT( HOR[1,J+1],A+C,A+C,'NIL','NIL'); RHS[I ,J+1]+:= NF; RHS[I+1,J+1] -:= NF 'OD'; HOR[0,]:= HOR[2,]; 'FOR' K 'TO' 2 'DO' II:= I+K-1; WE[1,K]:= OOST (VER[K,JU],XX ,YYO,WE[2,K],'NIL'); NF := FLUXT(VER[K,JU],LEFT[K],WE[1,K],'NIL','NIL'); RHS[II,JU]+:= NF 'OD' 'OD'; YY:= (JL-1.5)*DDY; 'FOR' J 'FROM' JL 'TO' JU 'DO' QOUT := ZUID (HOR[2,J],XXZ,YY+:=DDY,OUTW[J],'NIL'); NF := FLUXT(HOR[2,J],OUTW[J],QOUT,'NIL','NIL'); RHS[IU,J]+:= NF 'OD';'SKIP' 'FI' ; #$E# #E$# # SECOND ORDER RESIDUAL (III) # 'PR' EJECT 'PR' 'REAL' KAPPA:= 1/3; #$H# 'PROC' RESIDUVL = ('BOOL'PLUS,'FIELD' QF,RF)'FIELD': #$B RESIVL B$# 'BEGIN' 'REF'[,]'STATE' Q = Q'OF'QF, RHS = Q'OF'RF; 'REAL' DDX = DDX'OF'QF, DDY = DDY'OF'QF; 'INT' IL=1'LWB'Q,IU=1'UPB'Q,JL=2'LWB'Q,JU=2'UPB'Q; 'REAL' XX:= (IL-1.5)*DDX, YY:= (JL-1.5)*DDY; [JL-1:JU]'WALL' HORIN,HORIZ,VERI; [JL-1:JU]'POINT' COORD; 'REAL' KP = (1+KAPPA)/4, KM = (1-KAPPA)/4; [JL:JU]'STATE' QNO,DIM; STATE Z ( QF ); 'STATE' DJM,DIP,DJP, QNI,QZI, QWO,QWI,QOI, QIJ; 'FLUX' NF; # KAPPA = -1 : ZUIVER UPWIND, KAPPA = +1 : CENTRALE DIFFERENTIES (MET DISCONTINUE AANSLUITING OP DE RANDWAARDEN, VOLGENS DE RIEMANN INVARIANTEN) BETEKENIS VAN DE VARIABELEN: ---------> J ^ ! NO ! ! .----------------- ! DIM ! ! NI ! ! ! ! ! ! ! ! WO ! WI QIJ OI V I ! ! ^ V ! ! ! ZI ! ! ! DIP <--------------> <-----------.--------> DJM DJP ! ! ! V # WALLS('TRUE',PLUS,(IL-1)*DDX,DDY,COORD,'NIL',HORIN); 'FOR' I 'FROM' IL 'TO' IU 'DO' 'REF'[]'STATE' QI = Q[I,] , RHSI = RHS[I,] ; WALLS('TRUE',PLUS,I*DDX,DDY,COORD, VERI,HORIZ); 'FOR' J 'FROM' JL 'TO' JU 'DO' 'FLUX' QIJ = QI[J]; # VERT. WALLS # DJP:= ( J=JU ! DJM ! QI[J+1] - QIJ); ( J=JL ! QWI:= QIJ -0.5*DJP; QOI:= QIJ +0.5*DJP; QWO:= WEST(VERI[JL-1],XX+:=DDX,YYW,QWI,'NIL'); NF := FLUXT(VERI[J-1],QWO,QWI,'NIL','NIL') ! QWI:= QIJ - SCHEME(DJP,DJM); QOI:= QIJ + SCHEME(DJM,DJP); NF := FLUXT(VERI[J-1],QWO,QWI,'NIL','NIL'); RHSI[J-1]+:= NF );RHSI[J ]-:= NF; DJM:= DJP; QWO:= QOI; # HORIZ. WALLS # DIP:= ( I=IU ! DIM[J] ! Q[I+1,J] - QIJ); ( I=IL ! QNI:= QIJ -0.5*DIP; QZI:= QIJ +0.5*DIP; QNO[J]:= NOORD(HORIN[J],XXN,YY+:=DDY,QNI,'NIL'); NF := FLUXT(HORIN[J],QNO[J],QNI,'NIL','NIL') ! QNI:= QIJ - SCHEME(DIP ,DIM[J]); QZI:= QIJ + SCHEME(DIM[J],DIP ); NF := FLUXT(HORIN[J],QNO[J],QNI,'NIL','NIL'); RHS[I-1,J]+:= NF );RHSI [J]-:= NF; DIM[J]:= DIP; QNO[J]:= QZI 'OD'; QWI := OOST (VERI[JU], XX,YYO, QWO ,'NIL'); NF := FLUXT(VERI[JU],QWO,QWI,'NIL','NIL'); RHSI[JU]+:= NF; HORIN:= HORIZ 'OD'; YY:= (JL-1.5)*DDY; 'FOR' J 'FROM' JL 'TO' JU 'DO' QNI := ZUID (HORIZ[J],XXZ,YY+:=DDY,QNO[J],'NIL'); NF := FLUXT(HORIZ[J],QNO[J],QNI,'NIL','NIL'); RHS[IU,J]+:= NF 'OD'; 'SKIP' 'END'; #$E# #E$# # DIFFERENT VAN LEER TYPE SCHEMES # 'PR' EJECT 'PR' 'PROC' SCHEME := ('STATE' DM,DP)'STATE': (0.25*(1.0-KAPPA))*DM + (0.25*(1.0+KAPPA))*DP; #$H# 'PROC' ONESIDED = ('STATE' DM,DP)'STATE': #$B SHM.ONE B$# 'BEGIN' 0.5*DM 'END' ; #$E# #E$# #$H# 'PROC' FROMM = ('STATE' DM,DP)'STATE': #$B SHM.FRM B$# 'BEGIN' 0.25*(DM+DP) 'END' ; #$E# #E$# #$H# 'PROC' UPWIND BIASED = ('STATE' DM,DP)'STATE': #$B SHM.BIA B$# 'BEGIN' (1/3)*DM + (2/3)*DP 'END' ; #$E# #E$# #$H# 'PROC' CENTERED = ('STATE' DM,DP)'STATE': #$B SHM.CEN B$# 'BEGIN' 0.5*DP 'END' ; #$E# #E$# #$H# 'PROC' ALBADA = ('STATE' DM,DP)'STATE': #$B SHM.ALB B$# 'BEGIN' [1:4]'REAL' U, 'REAL' A,B,R,F, []'REAL' M = (C1'OF'DM, C2'OF'DM, C3'OF'DM, C4'OF'DM), P = (C1'OF'DP, C2'OF'DP, C3'OF'DP, C4'OF'DP); 'FOR' K 'TO' 4 'DO' A:= P[K]; B:= M[K]; ( 'ABS'A > 'ABS'B ! R:=B/A; R:=R*R; F:= B+A*R ! R:=A/B; R:=R*R; F:= A+B*R ); U[K]:= F/(2.0 + 2.0*R) 'OD'; 'STATE' (U[1],U[2],U[3],U[4]) 'END'; #$E# #E$# # FIRST ORDER STAR-RESIDUAL # 'PR' EJECT 'PR' #$H# 'PROC' STARRES1 = ('BOOL'PLUS,'FIELD' QF,RF)'FIELD': #$B STAR1 B$# 'IF' 'REF'[,]'STATE' Q = Q'OF'QF, R = Q'OF'RF; 'REAL' DDX = DDX'OF'QF, DDY = DDY'OF'QF; 'INT' IL=1'LWB'Q,IU=1'UPB'Q,JL=2'LWB'Q,JU=2'UPB'Q; 'INT' ILH = (IL+1)'OVER'2, IUH = IU'OVER'2, JLH = (JL+1)'OVER'2, JUH = JU'OVER'2; 'NOT'(ILH*2=IL+1'AND' IUH*2=IU'AND' JLH*2=JL+1'AND' JUH*2=JU ) 'THEN' RESIDU1(PLUS,QF,RF) 'ELSE' 'REAL' DX, DY, DS, XX,YY, DDXH:= 2*DDX, DDYH:= 2*DDY; 'WALL' CW; [JLH-1:JUH]'POINT' COORDN, COORDZ; [JLH-1:JUH] 'WALL' HORIN, HORIZ, VERT; 'FLUX' NF; 'STATE' QOUT; STATE Z ( QF ); WALLS('TRUE',PLUS,(ILH-1)*DDXH,DDYH,COORDZ,'NIL',HORIN); 'FOR' IH 'FROM' ILH 'TO' IUH 'DO' 'INT' I = 2*IH-1; XX:= I*DDX; COORDN := COORDZ; WALLS('TRUE',PLUS,IH*DDXH,DDYH,COORDZ,VERT,HORIZ); 'FOR' JH 'FROM' JLH 'TO' JUH 'DO' 'INT' J = 2*JH-1, JMH = JH-1; 'REF''STATE' QA=Q[I ,J ], QB=Q[I ,J+1], QC=Q[I+1,J+1], QD=Q[I+1,J ]; 'REF''FLUX' RA=R[I ,J ], RB=R[I ,J+1], RC=R[I+1,J+1], RD=R[I+1,J ]; 'POINT' MID = MAPPING(XX,YY:=J*DDY); 'PROC' DIAGH = ('POINT' CP)'WALL': 'BEGIN' DX := X'OF'CP - X'OF'MID; DY := Y'OF'CP - Y'OF'MID; DS := SQRT(DX*DX+DY*DY); 'WALL'(DY/DS,-DX/DS,(PLUS!DS!-DS) ) 'END'; ( IH=ILH ! QOUT :=NOORD(HORIN[JH],XXN,YY,QB,'NIL'); RB -:=FLUXT(HORIN[JH],QOUT ,QB,'NIL','NIL') ! NF :=FLUXT(HORIN[JH],Q[I-1,J],QB,'NIL','NIL'); RB -:= NF; R[I-1,J] +:= NF ); (JH=JLH ! QOUT :=WEST (VERT[JMH],XX,YYW,QA,'NIL'); RA -:=FLUXT(VERT[JMH],QOUT,QA,'NIL','NIL') ! NF :=FLUXT(VERT[JMH],Q[I+1,J-1],QA,'NIL','NIL'); RA -:= NF; R[I+1,J-1] +:= NF ); CW:= DIAGH(COORDN[JMH]); NF:= FLUXTG(CW,QA,QB,'NIL','NIL'); RA+:= NF; RB-:= NF; CW:= DIAGH(COORDN[JH ]); NF:= FLUXTG(CW,QB,QC,'NIL','NIL'); RB+:= NF; RC-:= NF; CW:= DIAGH(COORDZ[JH ]); NF:= FLUXTG(CW,QC,QD,'NIL','NIL'); RC+:= NF; RD-:= NF; CW:= DIAGH(COORDZ[JMH]); NF:= FLUXTG(CW,QD,QA,'NIL','NIL'); RD+:= NF; RA-:= NF 'OD';HORIN := HORIZ; QOUT := OOST (VERT[JUH],XX,YYO,Q[2*IH,JU],'NIL'); R[I+1,JU]+:= FLUXT(VERT[JUH],Q[I+1,JU],QOUT,'NIL','NIL') 'OD'; 'FOR' JH 'FROM' JLH 'TO' JUH 'DO' 'INT' J= 2*JH-1; YY:= J*DDY; QOUT := ZUID (HORIZ[JH],XXZ,YY,Q[IU,J],'NIL'); R[IU,J] +:= FLUXT(HORIZ[JH],Q[IU,J],QOUT,'NIL','NIL') 'OD';'SKIP' 'FI' ; #$E# #E$# # SECOND ORDER STAR-RESIDUAL # 'PR' EJECT 'PR' #$H# 'PROC' STARRES2 = ('BOOL'PLUS,'FIELD' QF,RF)'FIELD': #$B STAR2 B$# 'IF' 'REF'[,]'STATE' Q = Q'OF'QF, R = Q'OF'RF; 'REAL' DDX = DDX'OF'QF, DDY = DDY'OF'QF; 'INT' IL=1'LWB'Q,IU=1'UPB'Q,JL=2'LWB'Q,JU=2'UPB'Q; 'INT' ILH = (IL+1)'OVER'2, IUH = IU'OVER'2, JLH = (JL+1)'OVER'2, JUH = JU'OVER'2; 'NOT'(ILH*2=IL+1'AND' IUH*2=IU'AND' JLH*2=JL+1'AND' JUH*2=JU ) 'THEN' RESIDU1(PLUS,QF,RF) 'ELSE' 'REAL' DX, DY, DS, XX,YY, DDXH:= 2*DDX, DDYH:= 2*DDY; 'WALL' CW; [JLH-1:JUH]'POINT' COORDN, COORDZ; [JLH-1:JUH] 'WALL' HORIN, HORIZ, VERT; 'FLUX' NF; 'STATE' QIA,QIB,QIC,QID,QUO, QDX,QDY, QVER, QOUT; [JLH:JUH]'STATE' QHOR; STATE Z ( QF ); WALLS('TRUE',PLUS,(ILH-1)*DDXH,DDYH,COORDZ,'NIL',HORIN); 'FOR' IH 'FROM' ILH 'TO' IUH 'DO' 'INT' I = 2*IH-1; XX:= I*DDX; COORDN := COORDZ; WALLS('TRUE',PLUS,IH*DDXH,DDYH,COORDZ,VERT,HORIZ); 'FOR' JH 'FROM' JLH 'TO' JUH 'DO' 'INT' J = 2*JH-1, JMH = JH-1; YY:= J*DDY; 'REF''STATE' QA=Q[I ,J ], QB=Q[I ,J+1], QC=Q[I+1,J+1], QD=Q[I+1,J ]; 'REF''FLUX' RA=R[I ,J ], RB=R[I ,J+1], RC=R[I+1,J+1], RD=R[I+1,J ]; 'POINT' MID = MAPPING(XX,YY); 'PROC' DIAGH = ('POINT' CP)'WALL': 'BEGIN' DX := X'OF'CP - X'OF'MID; DY := Y'OF'CP - Y'OF'MID; DS := SQRT(DX*DX+DY*DY); 'WALL'(DY/DS,-DX/DS,(PLUS!DS!-DS) ) 'END'; QIA:= QIB:= QIC:= QID:= 0.25*(QA+QB+QC+QD); QDX:= 0.75*(QD-QB); QDY:= 0.75*(QC-QA); QIA-:= QDY; QIB-:= QDX; QIC+:= QDY; QID+:= QDX; ( IH=ILH ! QOUT :=NOORD(HORIN[JH],XXN,YY,QIB,'NIL'); RB -:=FLUXT(HORIN[JH],QOUT ,QIB,'NIL','NIL') ! NF :=FLUXT(HORIN[JH],QHOR[JH],QIB,'NIL','NIL'); RB -:= NF; R[I-1,J] +:= NF ); QHOR[JH]:= QID; (JH=JLH ! QOUT :=WEST (VERT[JMH],XX,YYW,QIA,'NIL'); RA -:=FLUXT(VERT[JMH],QOUT,QIA,'NIL','NIL') ! NF :=FLUXT(VERT[JMH],QVER,QIA,'NIL','NIL'); RA -:= NF; R[I+1,J-1] +:= NF ); QVER:=QIC; NFLUX:= ('STATE' Q0,Q1, 'REF''FLUXAC' F0,F1)'FLUX': ZFLUX(0.5*(Q0+Q1),'NIL'); CW:= DIAGH(COORDN[JMH]); NF:= FLUXTG(CW,QIA,QIB,'NIL','NIL'); RA+:= NF; RB-:= NF; CW:= DIAGH(COORDN[JH ]); NF:= FLUXTG(CW,QIB,QIC,'NIL','NIL'); RB+:= NF; RC-:= NF; CW:= DIAGH(COORDZ[JH ]); NF:= FLUXTG(CW,QIC,QID,'NIL','NIL'); RC+:= NF; RD-:= NF; CW:= DIAGH(COORDZ[JMH]); NF:= FLUXTG(CW,QID,QIA,'NIL','NIL'); RD+:= NF; RA-:= NF; NFLUX:= OSHER 'OD';HORIN := HORIZ; QOUT := OOST (VERT[JUH],XX,YYO,QIC,'NIL'); R[I+1,JU]+:= FLUXT(VERT[JUH],QVER,QOUT,'NIL','NIL') 'OD'; 'FOR' JH 'FROM' JLH 'TO' JUH 'DO' 'INT' J= 2*JH-1; YY:= J*DDY; QOUT := ZUID (HORIZ[JH],XXZ,YY,QID,'NIL'); R[IU,J] +:= FLUXT(HORIZ[JH],QHOR[JH],QOUT,'NIL','NIL') 'OD';'SKIP' 'FI' ; #$E# #E$# # SECOND ORDER OPERATORS # 'PR' EJECT 'PR' #$H# 'OP' 'RESIDU2' = ('FIELD' Q)'FIELD': #$B RESD2 B$# 'BEGIN' 'FIELD' F = 'NULRHS'Q; RESIDU2('TRUE',Q,F); F 'END'; #$E# #E$# #$H# 'OP' 'CORRECT' = ('FIELD' QL,RL)'FIELD': #$B CORREC B$# 'BEGIN' RESIDU ( 'TRUE',QL,RL); RESIDU2('FALSE',QL,RL); RL 'END'; #$E# #E$# #$H# 'OP' 'TAU' = ('FIELD' Q)'FIELD': #$B TAU.FL B$# 'BEGIN' 'RESTRICT''RESIDU'Q'MINRESIDU''RESTRICT'Q 'END' ; #$E# #E$# # LINEAR SYSTEM CONCTRUCTION # 'PR' EJECT 'PR' #$H# 'PROC' LINSYS = ('BOOL'PLUS,'FIELD' QF,RF, 'REF''NETMAT' JACN)'VOID': #$B LINSYS B$# 'BEGIN' 'REF'[,]'STATE' Q = Q'OF'QF, RHS = Q'OF'RF; 'INT' IL=1'LWB'Q,IU=1'UPB'Q,JL=2'LWB'Q,JU=2'UPB'Q; 'REAL' DDX = DDX'OF'QF, DDY = DDY'OF'QF; JACN:= 'HEAP'[IL:IU]'LINMAT'; 'NETMAT' JAC = JACN; 'REAL' XX:= (IL-1.5)*DDX; [JL-1:JU] 'WALL' HORI,VERI; [JL-1:JU]'POINT' COORD; 'STATE' QOUT; STATE Z ( QF ); 'STATAC' AC; 'FLUX' NF; 'FLUXAC' NFAC0,NFAC1; ( 'REF'[] 'STATE' QI = Q[IL,], RHSI = RHS[IL,]; 'REF'[,]'FLUXAC' JACI = GEN(JAC[IL],JL,JU,-2,+2); 'NUL' JACI; WALLS('TRUE',PLUS,(IL-1)*DDX,DDY,COORD,'NIL',HORI); 'REAL' YY:= (JL-1.5)*DDY; 'FOR' J 'FROM' JL 'TO' JU 'DO' QOUT := NOORD(HORI[J],XXN,YY+:=DDY,QI[J],AC); NF := FLUXT(HORI[J],QOUT,QI[J],NFAC0,NFAC1); NFAC1+:= NFAC0*AC; RHSI[ J]-:= NF; JACI[J,0]-:= NFAC1 'OD' ); 'FOR' I 'FROM' IL 'TO' IU 'DO' 'REF'[]'STATE' QI = Q[I,], RHSI = RHS[I,]; 'REF'[,]'FLUXAC' JACI = GET(JAC[I]); WALLS('TRUE',PLUS,I*DDX,DDY,COORD, VERI,HORI); QOUT := WEST(VERI[JL-1],XX+:=DDX,YYW,QI[JL],AC); NF := FLUXT(VERI[JL-1],QOUT,QI[JL],NFAC0,NFAC1); NFAC1+:= NFAC0*AC; RHSI[ JL]-:= NF; JACI[JL,0]-:= NFAC1; 'FOR' J 'FROM' JL 'TO' JU-1 'DO' # VERT. WALLS # NF:= FLUXT(VERI[J],QI[J],QI[J+1],NFAC0,NFAC1); RHSI[ J ]+:= NF; RHSI[ J+1]-:= NF; JACI[J , 0]+:= NFAC0; JACI[J , 1]+:= NFAC1; JACI[J+1, 0]-:= NFAC1; JACI[J+1,-1]-:= NFAC0 'OD'; QOUT := OOST(VERI[JU], XX,YYO,QI[JU],AC); NF := FLUXT(VERI[JU],QI[JU],QOUT,NFAC0,NFAC1); NFAC0+:= NFAC1*AC; RHSI [JU ]+:= NF; JACI [JU,0]+:= NFAC0; 'IF' I < IU 'THEN' 'REF'[]'STATE' QIP = Q[I+1,], RHSIP = RHS[I+1,]; 'REF'[,]'FLUXAC' JACIP = GEN(JAC[I+1],JL,JU,-2,2); 'NUL' JACIP; 'FOR' J 'FROM' JL 'TO' JU 'DO' #HORIZ. WALLS # NF:= FLUXT(HORI[J],QI[J],QIP[J],NFAC0,NFAC1); RHSI [J ]+:= NF; RHSIP[J ]-:= NF; JACI [J, 0]+:= NFAC0 + DIAG(QI[J]); JACI [J, 2]+:= NFAC1; JACIP[J, 0]-:= NFAC1; JACIP[J,-2]-:= NFAC0 'OD'; PET('TRUE',JAC[I]) 'FI' 'OD'; ( 'REF'[] 'STATE' QI = Q[IU,], RHSI = RHS[IU,]; 'REF'[,]'FLUXAC' JACI = GET(JAC[IU]); 'REAL' YY:= (JL-1.5)*DDY; 'FOR' J 'FROM' JL 'TO' JU 'DO' QOUT := ZUID(HORI[J],XXZ,YY+:=DDY,QI[J],AC); NF := FLUXT(HORI[J],QI[J],QOUT,NFAC0,NFAC1); NFAC0+:= NFAC1*AC; RHSI[J ]+:= NF; JACI[J,0]+:= NFAC0 + DIAG(QI[J]) 'OD';PET('TRUE',JAC[IU]) ) 'END'; #$E# #E$# 'PROC' DIAG := ('STATE' Q) 'FLUXAC': 'BEGIN' PRINT((NEWLINE,">>> DIAG STILL UNDEFINED <<<", NEWLINE,">>> ADVISE : DO DIAG:=DIAG1 <<<", NEWLINE)); ERROR ; NULFLUXAC 'END' ; #$H# 'PROC' DIAG1= ('STATE' Q) 'FLUXAC': #$B DIAG1 B$# 'BEGIN' 'REAL' U:=C1'OF'Q, V:=C2'OF'Q, C:=C3'OF'Q, Z:=C4'OF'Q, P,R,ROC,ROG,W2,CG; REVEAL(C,Z,P,R); R *:= ODIAGSTEP; ROC:= TOGAMM1*R/C; ROG:= -R*OGAMM1; W2 := 0.5*(U*U+V*V); CG := C*C*OGAMM1; 'FLUXAC'( 0 ,0 ,ROC ,ROG , R ,0 ,ROC*U ,ROG*U , 0 ,R ,ROC*V ,ROG*V , R*U ,R*V ,ROC*(W2+CG),ROG*(W2+CG*OGAM) ) 'END'; #$E# #E$# #$H# 'PROC' PLUS DIAG = ('FIELD' QF,DELTAQ,RF) 'VOID': #$B PLUSD B$# 'BEGIN' 'REF'[,]'STATE' Q =Q'OF'QF, DEL = Q'OF'DELTAQ, RHS =Q'OF'RF; 'INT' IL=1'LWB'Q,IU=1'UPB'Q,JL=2'LWB'Q,JU=2'UPB'Q; 'FOR' I 'FROM' IL 'TO' IU 'DO''REF'[]'STATE' QI=Q[I,], DELI=DEL[I,], RHSI=RHS[I,]; 'FOR' J 'FROM' JL 'TO' JU 'DO' RHSI[J] +:= DIAG(QI[J])*DELI[J] 'OD' 'OD' 'END'; #$E# #E$# # GALERKIN OPERATOR CONSTRUCTION # 'PR' EJECT 'PR' #$H# 'OP' 'RAP' = ('NETMAT' F) 'NETMAT': #$B RAP B$# 'BEGIN'# BY PDZ; INSERTED 850102 # 'INT' ILF = 'LWB'F, IUF = 'UPB'F; 'INT' ILC = (ILF+1)'OVER'2, IUC = IUF'OVER'2; 'NETMAT' C = 'HEAP'[ILC:IUC]'LINMAT'; 'INT' I:= ILF; 'FOR' IC 'FROM' ILC 'TO' IUC 'DO' 'REF'[,]'FLUXAC' FI = GET(F[I]), FIP = GET(F[I+1]); 'INT' JLF = 1'LWB'FI, JUF = 1'UPB'FI; 'INT' JLC = (JLF+1)'OVER'2, JUC = JUF'OVER'2; 'REF'[,]'FLUXAC' CIC = GEN(C[IC],JLC,JUC,-2,+2); 'NUL'CIC; 'INT' J:= JLF; 'FOR' JC 'FROM' JLC 'TO' JUC 'DO' 'REF'[]'FLUXAC' FIJ = FI [J ,], FIJP = FI [J+1,], FIPJ = FIP[J ,], FIPJP = FIP[J+1,], CICJC = CIC[JC,]; (CICJC[-2]+:=FIJ [-2])+:=FIJP [-2]; (CICJC[-1]+:=FIJ [-1])+:=FIPJ [-1]; ((CICJC[ 0]+:=FIJ [ 0])+:=FIJ [ 1])+:=FIJ [ 2]; ((CICJC[ 0]+:=FIJP [ 0])+:=FIJP [-1])+:=FIJP [ 2]; ((CICJC[ 0]+:=FIPJ [ 0])+:=FIPJ [ 1])+:=FIPJ [-2]; ((CICJC[ 0]+:=FIPJP[ 0])+:=FIPJP[-1])+:=FIPJP[-2]; (CICJC[ 1]+:=FIJP [ 1])+:=FIPJP[ 1]; (CICJC[ 2]+:=FIPJ [ 2])+:=FIPJP[ 2]; J+:= 2 'OD' ; PET('FALSE',F[I ]); PET('FALSE',F[I+1]); PET('TRUE' ,C[IC]); I+:= 2 'OD'; C 'END'; #$E# #E$# # RELAX # 'PR' EJECT 'PR' 'PROC' RELAX := ('INT' K,L, 'REF''FIELD' Q, RHS)'VOID': ERROR; 'BOOL' BACKWARD:= 'FALSE', REVERSE := 'FALSE', SYMMETRIC:= 'FALSE', RED := 'FALSE', REDBLACK:= 'FALSE'; 'REAL' DAMPING := 1, NLGSTOL := 0.99, TIMESTEP:= 1.0E10, ODIAGSTEP:= 0.00; 'REF'[]'REF'[,]'REAL' VOLUMES; # VOOR ELK LEVEL EEN 'REF'[,]'REAL'# #$H# 'PROC' STEPPING = ('INT' K,L, 'REF''FIELD' Q, RHS)'VOID': #$B STEPP B$# 'TO' K 'WHILE' WARNING 'DO' 'FIELD' R := RHS'MINRESIDU'Q; 'REF'[,]'STATE' QQ = Q'OF'Q, FF = Q'OF'R; ( VOLUMES[L] 'IS' 'NIL' ! VOLUMES[L]:= VOLUME(Q) ); 'REF'[,]'REAL' VOL = VOLUMES[L]; (MONI ! MONITOR("STEPPNG ",R) ); 'INT' IL=1'LWB'QQ,IU=1'UPB'QQ,JL=2'LWB'QQ,JU=2'UPB'QQ; 'FOR' I 'FROM' IL 'TO' IU 'DO' 'REF'[]'REAL' VOLI=VOL[I,]; 'REF'[]'STATE'QQI=QQ[I,],FFI=FF [I,]; 'FOR' J 'FROM' JL 'TO' JU 'DO' QQI[J]+:= (TIMESTEP/VOLI[J]) * FFI[J] 'OD' 'OD' 'OD' ; #$E# #E$# 'PROC' VOLM = ('POINT' A,B,C,D )'REAL' : 'ABS' (0.5 * ( (X'OF'C-X'OF'A) * (Y'OF'D-Y'OF'B) - (Y'OF'C-Y'OF'A) * (X'OF'D-X'OF'B) ) ) ; #$H# 'PROC' VOLUME = ('FIELD' A) 'REF'[,]'REAL': #$B VOLUME B$# 'BEGIN' 'REF'[,]'STATE' Q = Q'OF'A; 'INT' IL=1'LWB'Q,IU=1'UPB'Q,JL=2'LWB'Q,JU=2'UPB'Q; 'REAL' DDX = DDX'OF'A, DDY = DDY'OF'A; 'HEAP' [IL:IU,JL:JU]'REAL' VOL; [JL-1:JU]'POINT' COORD ; 'POINT' P1,P2 ; 'REAL' XX:= (IL-1)*DDX,YY:=(JL-2)*DDY ; 'FOR' J 'FROM' JL-1 'TO' JU 'DO' COORD[J]:= MAPPING(XX,YY+:=DDY) 'OD' ; 'FOR' I 'FROM' IL 'TO' IU 'DO' YY:=(JL-1)*DDY ; 'REF'[]'REAL' VOLI=VOL[I,] ; P1:=MAPPING(XX+:=DDX,YY) ; 'FOR' J 'FROM' JL 'TO' JU 'DO' P2 := MAPPING(XX,YY+:=DDY) ; VOLI[J]:= VOLM ( COORD[J-1],P1,P2,COORD[J] ) ; COORD[J-1]:=P1 ; P1:=P2 'OD';COORD[JU]:=P2 'OD'; VOL 'END'; #$E# #E$# 'PR' EJECT 'PR' #$H# 'PROC' SOL4 = ([,]'REF''FLUXAC' A, []'REF''FLUX' R)'VOID': # NB: OVERSCHRIJFT A EN R !!!!!!!!!!!!!!! # #$B SOL4 B$# 'BEGIN' 'FOR' I 'TO' 4 'DO' A[I,I]:= 1/A[I,I]; 'FOR' J 'FROM' I+1 'TO' 4 'DO' A[I,J] := A[I,I]*A[I,J]; 'FOR' K 'FROM' I+1 'TO' 4 'DO' A[K,J]-:= A[K,I]*A[I,J] 'OD' 'OD'; R[I ] := A[I,I]*R[I ]; 'FOR' K 'FROM' I+1 'TO' 4 'DO' R[K ]-:= A[K,I]*R[I ] 'OD' 'OD'; 'FOR' I 'FROM' 3 'BY' -1 'TO' 1 'DO''FOR' J 'FROM' I+1 'TO' 4 'DO' R[I ]-:= A[I,J]*R[J ] 'OD''OD' 'END'; #$E# #E$# #$H# 'PROC' RB4 RELAX = ('INT' K,L, 'REF''FIELD' QF, RF)'VOID': #$B RB4REL B$# 'BEGIN' 'REF'[,]'STATE' Q = Q'OF'QF; 'INT' L1= 1'LWB'Q, U1= 1'UPB'Q, L2= 2'LWB'Q, U2= 2'UPB'Q; 'BOOL' BLACK:= 'FALSE'; 'REAL' DAMP = DAMPING; 'REF''NETMAT' JAC= JACOBIAN[L]; 'FIELD' RHS; 'INT' IP,JP; [-2:2]'FLUXAC' JIJ,JIJP,JIPJ,JIPJP; 'FLUXAC' NUL1,NUL2,NUL3,NUL4; 'FLUX' RIJ,RIJP,RIPJ,RIPJP; [,]'REF''FLUXAC'JS=((JI J [ 0],JI J [ 1],JI J [ 2], NUL1 ), (JI JP[-1],JI JP[ 0], NUL2 ,JI JP[ 2]), (JIPJ [-2], NUL3 ,JIPJ [ 0],JIPJ [ 1]), ( NUL4 ,JIPJP[-2],JIPJP[-1],JIPJP[ 0])); []'REF''FLUX' RS=(RIJ,RIJP,RIPJ,RIPJP); ('ODD'(1+U1-L1) 'OR' 'ODD'(1+U2-L2) ! ERROR ); 'FOR' KK 'TO' ( REDBLACK ! 2*K ! K ) 'DO' # HET VOLGENDE KAN VEEL EFFICIENTER !! # 'FIELD' QFOLD = ( MONILIN ! 'COPY'QF! 'SKIP'); ( NEWLINSYS[L] ! LINSYS('FALSE',QF,RHS:= 'COPY'RF, JAC) ! RHS:= RF'MINRESIDU'QF ); ( MONI ! MONITOR("RB4-RELAX ",RHS) ); 'REF'[,]'STATE' R=Q'OF'RHS; ( REDBLACK ! BLACK:= RED:= 'NOT'RED ); 'FOR' I 'FROM' L1 'BY' 2 'TO' U1 'DO' 'REF'[,]'FLUXAC' JACI = GET(JAC[I ]), JACIP= GET(JAC[I+1]); 'FOR' J 'FROM' ( REDBLACK 'AND' BLACK ! L2+2 ! L2 ) 'BY' ( REDBLACK ! 4 ! 2 ) 'TO' U2 'DO' IP:=I+1; JP:= J+1; JI J := JACI [J ,]; RI J := R[I ,J ]; JI JP:= JACI [JP,]; RI JP:= R[I ,JP]; JIPJ := JACIP[J ,]; RIPJ := R[IP,J ]; JIPJP:= JACIP[JP,]; RIPJP:= R[IP,JP]; NUL1:= NUL2:= NUL3:= NUL4:= NULFLUXAC; SOL4 ( JS, RS); 'REAL' FACTOR:= DAMP; 'WHILE' FACTOR*C3'OF'RI J /C3'OF'Q[I ,J ]<-0.99 'OR' FACTOR*C3'OF'RIPJ /C3'OF'Q[IP,J ]<-0.99 'OR' FACTOR*C3'OF'RI JP/C3'OF'Q[I ,JP]<-0.99 'OR' FACTOR*C3'OF'RIPJP/C3'OF'Q[IP,JP]<-0.99 'DO' FACTOR/:= 2 'OD'; (FACTOR/=DAMP !PRINT((NEWLINE,"RB-WARNING",I,J,FACTOR))); Q[I ,J ]+:= FACTOR*RIJ ; Q[I ,JP]+:= FACTOR*RIJP; Q[IP,J ]+:= FACTOR*RIPJ; Q[IP,JP]+:= FACTOR*RIPJP 'OD';BLACK:= 'NOT'BLACK; PET('FALSE',JAC[I]); PET('FALSE',JAC[I+1]) 'OD'; ( MONILIN ! PRINT((NEWLINE,"LINEAR CHECK RB-RELAX",NEWLINE)); RHS-:= JAC*(QF-QFOLD); 'PRINT' RHS; MONITOR("RB-LINEAR ",RHS) ) 'OD' 'END'; #$E# #E$# 'PR' EJECT 'PR' #$H# 'PROC' RB2 RELAX = ('INT' K,L, 'REF''FIELD' QF, RF)'VOID': #$B RB2REL B$# 'BEGIN' 'REF'[,]'STATE' Q = Q'OF'QF; 'INT' L1= 1'LWB'Q, U1= 1'UPB'Q, L2= 2'LWB'Q, U2= 2'UPB'Q; 'BOOL' BLACK:= 'FALSE'; 'REAL' DAMP = DAMPING; 'REF''NETMAT' JAC= JACOBIAN[L]; 'FIELD' RHS; 'STATE' X; 'FOR' KK 'TO' ( REDBLACK ! 2*K ! K ) 'DO' # KAN EFFICIENTER # 'FIELD' QFOLD = ( MONILIN ! 'COPY'QF! 'SKIP'); ( NEWLINSYS[L] ! LINSYS('FALSE',QF,RHS:= 'COPY'RF, JAC) ! RHS:= RF'MINRESIDU'QF ); ( MONI ! MONITOR("RB2-RELAX ",RHS) ); 'REF'[,]'STATE' R=Q'OF'RHS; ( REDBLACK ! BLACK:= RED:= 'NOT' RED ); 'FOR' I 'FROM' L1 'TO' U1 'DO' 'REF'[,]'FLUXAC' JACI = GET(JAC[I ]); 'FOR' J 'FROM' ( REDBLACK 'AND' BLACK ! L2+1 ! L2 ) 'BY' ( REDBLACK ! 2 ! 1 ) 'TO' U2 'DO' X:= R[I,J]/JACI[J,0]; 'REAL' FACTOR:= DAMP; 'WHILE' FACTOR*C3'OF'X/C3'OF'Q[I,J]<-0.99 'DO' FACTOR/:= 2 'OD'; (FACTOR/=DAMP !PRINT((NEWLINE,"RB-WARNING",I,J,FACTOR))); Q[I,J]+:= ( FACTOR=1! X ! FACTOR*X ) 'OD';BLACK:= 'NOT'BLACK; PET('FALSE',JAC[I]) 'OD'; ( MONILIN ! PRINT((NEWLINE,"LINEAR CHECK RB-RELAX",NEWLINE)); RHS-:= JAC*(QF-QFOLD); 'PRINT' RHS; MONITOR("RB-LINEAR ",RHS) ) 'OD' 'END'; #$E# #E$# 'PR' EJECT 'PR' #$H# 'PROC' JACOBI = ('INT' K,L, 'REF''FIELD' QF, RHS)'VOID': #$B JACOBI B$# 'BEGIN' 'REF'[,]'STATE' Q = Q'OF'QF, FF = Q'OF'RHS; STATE Z ( QF ); 'REAL' DDX = DDX'OF'QF, DDY = DDY'OF'QF; 'INT' IL=1'LWB'FF,IU=1'UPB'FF,JL=2'LWB'FF,JU=2'UPB'FF; [JL-1:JU] 'WALL' HORIN,HORIZ,VERI ; [JL-1:JU] 'POINT' COORD ; [JL :JU] 'STATE' QIM,QIO ; 'STATE' UC, QOUT; 'STATAC' QAC ; 'FLUX' NF,F,DELTA ; 'FLUXAC' NFAC0,NFAC1,FAC ; 'REAL' ERROR,FACTOR,DAMP:=DAMPING ; 'TO' K 'DO' WALLS( 'TRUE','TRUE',(IL-1)*DDX,DDY,COORD,'NIL',HORIN); 'FOR' I 'FROM' IL 'TO' IU 'DO''REAL'XX=(I-0.5)*DDX; WALLS( 'TRUE','TRUE', I *DDX,DDY,COORD, VERI,HORIZ); QIO:= Q[I,]; 'FOR' J 'FROM' JL 'TO' JU 'DO''REAL'YY=(J-0.5)*DDY; 'REF''STATE' QIJ=Q[I,J]; UC:=QIJ; 'TO' 13 'WHILE' F :=(I=IU ! QOUT := ZUID(HORIZ[J],XXZ,YY,UC,QAC); FLUXT(HORIZ[J],UC,QOUT,NFAC0,NFAC1) ! FLUXT(HORIZ[J],UC,Q[I+1,J],NFAC0,'NIL')); FAC :=(I=IU ! NFAC0+:= NFAC1*QAC ! NFAC0 ); F -:=(I=IL ! QOUT := NOORD(HORIN[J],XXN,YY,UC,QAC); FLUXT(HORIN[J],QOUT,UC,NFAC0,NFAC1) ! FLUXT(HORIN[J],QIM[J],UC,'NIL',NFAC1)); FAC-:=(I=IL ! NFAC1+:= NFAC0*QAC ! NFAC1 ); F +:=(J=JU ! QOUT := OOST(VERI[J],XX,YYO,UC,QAC); FLUXT(VERI[J],UC,QOUT,NFAC0,NFAC1) ! FLUXT(VERI[J],UC,Q[I,J+1],NFAC0,'NIL')); FAC+:=(J=JU ! NFAC0+:= NFAC1*QAC ! NFAC0 ); F -:=(J=JL ! QOUT := WEST(VERI[J-1],XX,YYW,UC,QAC); FLUXT(VERI[J-1],QOUT,UC,NFAC0,NFAC1) ! FLUXT(VERI[J-1],QIO[J-1],UC,'NIL',NFAC1)); FAC-:=(J=JL ! NFAC1+:=NFAC0*QAC ! NFAC1 ); DELTA := FF[I,J]-F; ERROR := 'ABS'DELTA; DELTA := DELTA/FAC; FACTOR:= 1.0; 'WHILE' FACTOR * 'ABS'(C3'OF'DELTA/C3'OF'UC) > 0.99 'DO' FACTOR/:=2 'OD'; ( FACTOR/=1.0 ! DELTA*:= FACTOR ); UC+:= DELTA; (MONIGS! PRINT((NEWLINE, "NEWT."+SF(ERROR)+WHOLE(I,5)+WHOLE(J,5) )); (FACTOR/=DAMP! PRINT((" "+SF(FACTOR) )) )); ERROR > NLGSTOL 'DO' 'SKIP' 'OD'; QIJ+:= DAMP*(UC-QIJ) 'OD'; QIM:= QIO; HORIN:= HORIZ 'OD'; (CAVITATION! PRINT((NEWLINE,"CAVITATION OCCURRED",NEWLINE)); CAVITATION:= 'FALSE') 'OD'; (MONI ! MONITOR("JACOBI ",RHS'MINRESIDU'QF)) 'END'; #$E# #E$# 'PR' EJECT 'PR' #$H# 'PROC' SEIDEL = ('INT' K,L, 'REF''FIELD' QF, RHS)'VOID': #$B SEIDEL B$# 'BEGIN' 'REF'[,]'STATE' Q = Q'OF'QF, FF = Q'OF'RHS; STATE Z ( QF ); 'REAL' DDX = DDX'OF'QF, DDY = DDY'OF'QF; 'INT' IL=1'LWB'FF,IU=1'UPB'FF,JL=2'LWB'FF,JU=2'UPB'FF; 'INT' ISL,IS,ISU,JSL,JS,JSU; (BACKWARD! ISL:=IU;ISU:=IL;IS:=-1 ! ISL:=IL;ISU:=IU;IS:=1); (REVERSE ! JSL:=JU;JSU:=JL;JS:=-1 ! JSL:=JL;JSU:=JU;JS:=1); [JL-1:JU] 'WALL' HORIN,HORIZ,VERI; [JL-1:JU]'POINT' COORD; 'STATE' QOUT; 'STATAC' QAC; 'FLUX' NF,F,DELTA; 'FLUXAC' NFAC0,NFAC1,FAC; 'REAL' ERROR,FACTOR,DAMP:=DAMPING; 'FOR' KK 'TO' ( SYMMETRIC'OR'REDBLACK ! 2*K ! K) 'DO' 'BOOL' SHIFT:= RED = 'ODD'KK; ( IS=1 ! WALLS( 'TRUE','TRUE',(IL-1)*DDX,DDY,COORD,'NIL',HORIZ) ! WALLS('FALSE','TRUE', IU *DDX,DDY,COORD,'NIL',HORIN) ); 'FOR' I 'FROM' ISL 'BY' IS 'TO' ISU 'DO''REAL'XX=(I-0.5)*DDX; ( IS=1 ! HORIN:= HORIZ; WALLS( 'TRUE','TRUE', I *DDX,DDY,COORD, VERI,HORIZ) ! HORIZ:= HORIN; WALLS('FALSE','TRUE',( I-1)*DDX,DDY,COORD, VERI,HORIN) ); 'INT' JJSL:= JSL, JJS:= JS; (REDBLACK! (SHIFT:='NOT'SHIFT! JJSL+:=JS ); JJS+:=JS ); 'FOR' J 'FROM' JJSL 'BY' JJS 'TO' JSU 'DO''REAL'YY=(J-0.5)*DDY; 'TO' 13 'WHILE' 'REF''STATE' UC=Q[I,J]; F :=(I=IU ! QOUT := ZUID(HORIZ[J],XXZ,YY,UC,QAC); FLUXT(HORIZ[J],UC,QOUT,NFAC0,NFAC1) ! FLUXT(HORIZ[J],UC,Q[I+1,J],NFAC0,'NIL')); FAC :=(I=IU ! NFAC0+:= NFAC1*QAC ! NFAC0 ); F -:=(I=IL ! QOUT := NOORD(HORIN[J],XXN,YY,UC,QAC); FLUXT(HORIN[J],QOUT,UC,NFAC0,NFAC1) ! FLUXT(HORIN[J],Q[I-1,J],UC,'NIL',NFAC1)); FAC-:=(I=IL ! NFAC1+:= NFAC0*QAC ! NFAC1 ); F +:=(J=JU ! QOUT := OOST(VERI[J],XX,YYO,UC,QAC); FLUXT(VERI[J],UC,QOUT,NFAC0,NFAC1) ! FLUXT(VERI[J],UC,Q[I,J+1],NFAC0,'NIL')); FAC+:=(J=JU ! NFAC0+:= NFAC1*QAC ! NFAC0 ); F -:=(J=JL ! QOUT := WEST(VERI[J-1],XX,YYW,UC,QAC); FLUXT(VERI[J-1],QOUT,UC,NFAC0,NFAC1) ! FLUXT(VERI[J-1],Q[I,J-1],UC,'NIL',NFAC1)); FAC-:=(J=JL ! NFAC1+:=NFAC0*QAC ! NFAC1 ); DELTA := FF[I,J]-F; ERROR := 'ABS'DELTA; DELTA := DELTA/FAC; FACTOR:= DAMP; 'WHILE' FACTOR * 'ABS'(C3'OF'DELTA/C3'OF'UC) > 0.99 'DO' FACTOR/:=2 'OD'; ( FACTOR/=1 ! DELTA*:= FACTOR ); UC+:= DELTA; 'IF' MONIGS 'THEN' PRINT(( NEWLINE,"NEWT." + SF( ERROR ) + WHOLE(I,5) + WHOLE(J,5) )); (FACTOR/=DAMP!PRINT((" "+SF( FACTOR ) )) ) 'FI'; ERROR > NLGSTOL 'DO' 'SKIP' 'OD' 'OD' 'OD'; (CAVITATION! PRINT((NEWLINE,"CAVITATION OCCURRED",NEWLINE)); CAVITATION:= 'FALSE') ; 'IF' SYMMETRIC 'THEN' 'INT' IJ:= ISL; ISL:=ISU; ISU:=IJ; IS:=-IS; IJ:= JSL; JSL:=JSU; JSU:=IJ; JS:=-JS 'FI' 'OD'; (MONI ! MONITOR("SEIDEL ",RHS'MINRESIDU'QF)) 'END'; #$E# #E$# # STAR-RELAXATIE # 'BOOL' BLOCK ; #$H# 'PROC' STARRELAX = ('INT' K,L, 'REF''FIELD' QF, RHS)'VOID': #$B STRREL B$# 'BEGIN' 'REF'[,]'STATE' QQ = Q'OF'QF, FF = Q'OF'RHS; STATE Z ( QF ); 'REAL' DDXH = 2*DDX'OF'QF, DDYH = 2*DDY'OF'QF; 'INT' ILH=(1'LWB'FF+1)'OVER'2,IUH=(1'UPB'FF)'OVER'2, JLH=(2'LWB'FF+1)'OVER'2,JUH=(2'UPB'FF)'OVER'2; 'INT' ISL,IS,ISU,JSL,JS,JSU; (BACKWARD!ISL:=IUH;ISU:=ILH;IS:=-1 !ISL:=ILH;ISU:=IUH;IS:=1); (REVERSE !JSL:=JUH;JSU:=JLH;JS:=-1 !JSL:=JLH;JSU:=JUH;JS:=1); [JLH-1:JUH] 'WALL' HORIN,HORIZ,VERT; [JLH-1:JUH]'POINT' COORDN,COORDZ; 'REAL' ERROR,FACTOR, DAMP:=DAMPING*0.5; [1:4]'STATE' QS,DELTA; [1:4]'FLUX' RS; [1:4,1:4]'FLUXAC' JAC ; 'REF''STATE' QA= QS[ 1],QB= QS[ 2],QC= QS[ 3],QD= QS[ 4]; 'REF''FLUX' RA= RS[ 1],RB= RS[ 2],RC= RS[ 3],RD= RS[ 4]; 'REF''FLUXAC' AA=JAC[1,1],AB=JAC[1,2],AC=JAC[1,3],AD=JAC[1,4], BA=JAC[2,1],BB=JAC[2,2],BC=JAC[2,3],BD=JAC[2,4], CA=JAC[3,1],CB=JAC[3,2],CC=JAC[3,3],CD=JAC[3,4], DA=JAC[4,1],DB=JAC[4,2],DC=JAC[4,3],DD=JAC[4,4]; # NOORD (I-1,J) O---------------------O ! \ (I,J+1) / ! ! \ / ! B ! \ / ! ! ! \ / (I+1, ! ! (I+1,J-1) ! (I,J) O J+1) ! (I,J+2) A ---.--- C WEST ! / \ ! OOST ! ! / \ ! ! ! / \ ! D O / (I+1,J) \ ! O---------------------O (I+2,J+1) ZUID # #!# 'FOR' KK 'TO' ( SYMMETRIC'OR'REDBLACK ! 2*K ! K) 'DO' 'BOOL' SHIFT:= RED = 'ODD'KK; ( IS=1 !WALLS( 'TRUE','TRUE',(ILH-1)*DDXH,DDYH,COORDZ,'NIL',HORIZ) !WALLS('FALSE','TRUE', IUH *DDXH,DDYH,COORDN,'NIL',HORIN) ); #!# 'FOR' IH 'FROM' ISL 'BY' IS 'TO' ISU 'DO''REAL'XX=(IH-0.5)*DDXH; ( IS=1 ! HORIN:= HORIZ; COORDN:= COORDZ; WALLS( 'TRUE','TRUE', IH *DDXH,DDYH,COORDZ, VERT,HORIZ) ! HORIZ:= HORIN; COORDZ:= COORDN; WALLS('FALSE','TRUE',( IH-1)*DDXH,DDYH,COORDN, VERT,HORIN) ); 'INT' JJSL:= JSL, JJS:= JS; (REDBLACK! (SHIFT:='NOT'SHIFT! JJSL+:=JS ); JJS+:=JS ); #!# 'FOR' JH 'FROM' JJSL 'BY' JJS 'TO' JSU 'DO''REAL'YY=(JH-0.5)*DDYH; 'BOOL' OK:= BLOCK; 'INT' I = 2*IH-1, J = 2*JH-1, JMH = JH-1; 'POINT' MID = MAPPING(XX,YY); QS:= (QQ[I ,J ], QQ[I ,J+1], QQ[I+1,J+1], QQ[I+1,J ]); 'PROC' DIAGH = ('POINT' CP)'WALL': 'BEGIN' 'REAL' DX := X'OF'CP - X'OF'MID, DY := Y'OF'CP - Y'OF'MID, DS; DS := SQRT(DX*DX+DY*DY); 'WALL'(DY/DS,-DX/DS,DS) 'END'; #$!# 'TO' 8 'WHILE' 'PROC' RAJAC = 'VOID': 'BEGIN' # BEREKENT LOKAAL: RS:= "RHS-N(Q)" EN JAC:= "DN/DQ" # 'WALL' CW; 'STATE' QOUT; 'STATAC'QAC; 'FLUX' NF; 'FLUXAC' NF0, NF1; RA:= FF[I ,J ]; RB:= FF[I ,J+1]; RC:= FF[I+1,J+1]; RD:= FF[I+1,J ]; 'FOR' K 'TO' 4 'DO' 'FOR' L 'TO' 4 'DO' JAC[K,L]:= NULFLUXAC 'OD' 'OD'; ( IH=ILH ! QOUT :=NOORD(HORIN[JH],XXN,YY,QB,QAC); RB +:=FLUXT(HORIN[JH], QOUT ,QB,NF0,NF1); BB-:= (NF1+:=NF0*QAC) ! RB +:=FLUXT(HORIN[JH],QQ[I-1,J],QB,'NIL',NF1); BB-:= NF1 ); ( IH=IUH ! QOUT :=ZUID (HORIZ[JH],XXZ,YY,QD,QAC); RD -:=FLUXT(HORIZ[JH],QD,QOUT,NF0,NF1); DD+:= (NF0+:=NF1*QAC) ! RD -:=FLUXT(HORIZ[JH],QD,QQ[I+2,J+1],NF0,'NIL'); DD+:= NF0 ); (JH=JUH ! QOUT :=OOST (VERT[JH],XX,YYO,QC,QAC); RC -:=FLUXT(VERT[JH],QC,QOUT,NF0,NF1); CC+:= (NF0+:= NF1*QAC) ! RC -:=FLUXT(VERT[JH],QC,QQ[I,J+2],NF0,'NIL'); CC+:= NF0 ); (JH=JLH ! QOUT :=WEST (VERT[JMH],XX,YYW,QA,QAC); RA +:=FLUXT(VERT[JMH],QOUT,QA,NF0,NF1); AA-:= (NF1+:= NF0*QAC) ! RA +:=FLUXT(VERT[JMH],QQ[I+1,J-1],QA,'NIL',NF1); AA-:= NF1 ); CW:= DIAGH(COORDN[JMH]); NF:= FLUXTG(CW,QA,QB,NF0,NF1); RA-:=NF; RB+:=NF; AA+:=NF0; AB+:=NF1; BA-:=NF0; BB-:= NF1; CW:= DIAGH(COORDN[JH] ); NF:= FLUXTG(CW,QB,QC,NF0,NF1); RB-:=NF; RC+:=NF; BB+:=NF0; BC+:=NF1; CB-:=NF0; CC-:= NF1; CW:= DIAGH(COORDZ[JH] ); NF:= FLUXTG(CW,QC,QD,NF0,NF1); RC-:=NF; RD+:=NF; CC+:=NF0; CD+:=NF1; DC-:=NF0; DD-:= NF1; CW:= DIAGH(COORDZ[JMH]); NF:= FLUXTG(CW,QD,QA,NF0,NF1); RD-:=NF; RA+:=NF; DD+:=NF0; DA+:=NF1; AD-:=NF0; AA-:= NF1 'END' # RAJAC# ; RAJAC; ERROR := 'ABS'RA + 'ABS'RB + 'ABS'RC + 'ABS'RD; ( MONIGS ! PRINT((NEWLINE,"GS ERROR VOOR "+WHOLE(JUH-JLH+1,5)+ WHOLE(IH,5)+WHOLE(JH,5), ERROR)) ); 'IF' ( OK # NEWTON OR DAMPED JACOBI ? # ! DELTA := RS/JAC; 'FOR' C 'TO' 4 'WHILE' OK 'DO' OK:= 'ABS'(C3'OF'DELTA[C]/C3'OF'QS[C]) < 0.99 'OD' ); OK 'THEN' 'FOR' C 'TO' 4 # NEWTON GS # 'DO' QS[C]+:= DELTA[C] 'OD' 'ELSE' 'FOR' C 'TO' 4 # DAMPED JACOBI # 'DO' FACTOR:= DAMP; 'WHILE' 'ABS'(C3'OF'DELTA[C]/C3'OF'QS[C]) > 0.9/FACTOR 'DO' FACTOR/:=2 'OD'; ( FACTOR<=0.1 'OR' ( BLOCK'AND'MONIGS) ! PRINT((NEWLINE,"DAMPED JACOBI FACTOR",FACTOR)) ); ( FACTOR >0.1 ! QS[C]+:= ( DELTA[C]*:= FACTOR ) ) 'OD' 'FI'; # 'IF' ( 'FALSE'; MONIGS 'AND' (BLOCK'AND''NOT'OK) ); 'FALSE' 'THEN' [1:4]'FLUX' ROLD:= RS, RNEW; [1:4]'STATE' QOLD:= (QQ[I ,J ], QQ[I ,J+1], QQ[I+1,J+1], QQ[I+1,J ]); [1:4,1:4]'FLUXAC' JACOLD:= JAC; PRINT((NEWLINE,NEWLINE,"IH,JH:"+WHOLE(IH,0)+" "+WHOLE(JH,0), NEWLINE,"OPLOSSING VOOR CORRECTIE ", NEWLINE,QOLD, NEWLINE,"RESIDU VOOR CORRECTIE ",'ABS'ROLD, NEWLINE,ROLD)); ('REAL' S ='ABS'(ROLD-JAC*DELTA); S > 1.0E-10 ! PRINT((NEWLINE,"DAMPING OR BAD LIN. SOLVING",S)) ); RAJAC; RNEW:= RS; PRINT((NEWLINE,"RESIDU NA CORRECTIE ",'ABS'RNEW, NEWLINE,"CONTROLE NA DE SLAG A ",'ABS'(RNEW-ROLD), NEWLINE,RNEW-ROLD, NEWLINE,"CONTROLE NA DE SLAG B ",'ABS'(JACOLD*DELTA), NEWLINE,JACOLD*DELTA, NEWLINE,(RNEW-ROLD)/(JACOLD*DELTA), NEWLINE,"RELAT VERBETERING", 'ABS'(RNEW-ROLD+JACOLD*DELTA)/'ABS'(RNEW-ROLD) )) 'FI'; # ERROR > NLGSTOL 'AND' OK 'DO' 'SKIP' 'OD'; QQ[I ,J ]:= QA; QQ[I ,J+1]:= QB; QQ[I+1,J+1]:= QC; QQ[I+1,J ]:= QD 'OD' 'OD'; (CAVITATION! PRINT((NEWLINE,"CAVITATION OCCURRED",NEWLINE)); CAVITATION:= 'FALSE') ; 'IF' SYMMETRIC 'THEN' 'INT' IJ:= ISL; ISL:=ISU; ISU:=IJ; IS:=-IS; IJ:= JSL; JSL:=JSU; JSU:=IJ; JS:=-JS 'FI' 'OD'; (MONI ! MONITOR("STARRELAX ",RHS'MINRESIDU'QF)) 'END'; #$E# #E$# # NEWTON ITERATION # 'PR' EJECT 'PR' #$H# 'PROC' NEWTON = ('INT' K,L, 'REF''FIELD' QF,RF) 'VOID': #$B NEWTON B$# 'BEGIN''REF''NETMAT' JAC= JACOBIAN[L], DEC= DECOMP[L], 'FIELD' RHS, 'INT' IL= 1'LWB'(Q'OF'QF), IU= 1'UPB'(Q'OF'QF); 'REAL' ROLD:= MAXREAL, RNEW,LTOL; 'FOR' KK 'TO' K 'WHILE'( KK=1 'OR' NEW LIN SYS[L] ! (MONINWT !PRINT((NEWLINE,"UPDATE JACOBIANS ",NEWLINSYS))); LINSYS ('FALSE',QF,RHS:='COPY'RF,JAC); RNEW:= FLUXNORM(RHS,'LOC''FLUX','LOC''FLUX'); NEWLINSYS[L]:= 'FALSE'; ( ILLU ! ILLUDEC(JAC,DEC) ); 'FOR' LEV 'FROM' L-1 'BY'-1 'TO' 1 'WHILE' NEWLINSYS[LEV] 'DO' JACOBIAN [LEV]:= 'RAP'JACOBIAN[LEV+1]; NEWLINSYS[LEV]:= 'FALSE'; ( ILLU ! ILLUDEC(JACOBIAN[LEV],DECOMP[LEV]) ) 'OD' ); LTOL:= ('REAL' Z = (RNEW<1.0 ! RNEW*RNEW ! RNEW); (Z>LINRESTOL!Z!LINRESTOL) ); ( MONINWT ! PRINT((NEWLINE,"NEWTON: "+FLOAT(RNEW,10,3,3)+ "DESIRED RESIDUAL = "+FLOAT(LTOL,10,3,3) )) ); RESTOL1 ! PLUS DIAG (QF,W,RES) ); LINSOLVE(1,L,W,RES); NORM:= FLUXNORM(RES-:=JAC*W,'LOC''FLUX','LOC''FLUX'); (MONINWT! PRINT((NEWLINE,"LINSOLVE "+20*" "+ FLOAT(NORM,10,3,3)+" "+ FIXED(-4*LN(NORM)/LN(10.0),6,1) ))); ( MONI 'OR' NORM>=OLDNORM ! MONITOR("LINSOLV-NL",RF'MINRESIDU'(QF+W)) ) 'OD'; (MONINWT ! PRINT((NEWLINE,"LINSOLVE CONV = ", FLOAT( EXP((LN(NORM/OLDNORM))/CM),10,3,3) )) ); ( NORM > OLDNORM # DIVERGENCE # ! ODIAGSTEP*:= 10; LINEARIZE; PRINT((NEWLINE, "NEWTON STEP NOT ACCEPTED! ")) ! QF+:= W ); ( ONORM < 1.5*NORM # SLOW CONVERGENCE # ! ODIAGSTEP*:= 2; LINEARIZE; PRINT((NEWLINE,"STEP CHANGED:"+ SF( ODIAGSTEP )+50*"*",NEWLINE)) ); RHS := RF'MINRESIDU'QF; RNEW:= FLUXNORM(RHS,'LOC''FLUX','LOC''FLUX'); ( # VGL. LIN-NONLIN RESIDU: # RNEW/NORM > 1.2 ! LINEARIZE ) 'OD' 'END'; #$E# #E$# # APPROXIMATE LINEAR SOLVERS # 'PR' EJECT 'PR' 'REF'[]'NETMAT' JACOBIAN,DECOMP; 'REF'[]'BOOL' NEW LIN SYS; 'BOOL' ILLU:= 'TRUE'; 'INT' MLINSOLVE:= 13; 'REAL' LINRESTOL:= 1.0E-11, RESTOL:= 1.0E-11; 'PROC' LINSOLVE := ('INT' K,L, 'REF''FIELD' Q, R)'VOID':'SKIP'; #$H# 'PROC' LINEARIZE = 'VOID': #$B LNRIZE B$# 'FOR' I 'TO' 'UPB'NEWLINSYS 'DO' NEWLINSYS[I]:='TRUE' 'OD'; #$E# #E$# #$H# 'PROC' INIZ LIN = ('INT' LEVELS)'VOID': #$B INIZLIN B$# 'BEGIN' JACOBIAN := 'HEAP'[1:LEVELS]'NETMAT'; DECOMP := 'HEAP'[1:LEVELS]'NETMAT'; NEWLINSYS:= 'HEAP'[1:LEVELS]'BOOL'; LINEARIZE 'END'; #$E# #E$# # LINEAR RELAXATION # 'BOOL' LINSYM:= 'TRUE', LINREVER:= 'FALSE', LINBACKW:= 'FALSE'; 'REAL' DAMP:= 1.0; #$H# 'PROC' LU RELAX = ('INT' K, L, 'REF''FIELD' QF, RHS) 'VOID': #$B LUREL B$# 'BEGIN'# BY PDZ; INSERTED 850102 # 'TO' K 'DO' 'FIELD' RES:= RHS-JACOBIAN[L]*QF; ILLURELAX(DECOMP[L],JACOBIAN[L],QF,RES) 'OD' 'END'; #$E# #E$# 'PR' EJECT 'PR' #$H# 'PROC' LINJACO = ('INT' K, L, 'REF''FIELD' F, RHS) 'VOID': #$B LINJC B$# 'BEGIN'# BY PDZ; INSERTED 850102 # 'REF''NETMAT' JAC = JACOBIAN[L], 'REF'[,]'STATE' QR = Q'OF'RHS; 'INT' L1 = 1'LWB'QR, U1 = 1'UPB'QR, L2 = 2'LWB'QR, U2 = 2'UPB'QR; ( L1/= 'LWB'JAC 'OR' U1/= 'UPB'JAC ! ERROR ); 'REAL' D1 = 1/DAMP, D2 = (DAMP-1)/DAMP; 'TO' K 'DO' 'FIELD' G = 'COPY'RHS; TYPE'OF'G:= TYPE'OF'F; 'REF'[,]'STATE' QG = Q'OF'G, QF = Q'OF'F; ( L1/=1'LWB'QF 'OR' U1/=1'UPB'QF 'OR' L2/=2'LWB'QF 'OR' U2/=2'UPB'QF ! ERROR ); 'FOR' I 'FROM' L1 'TO' U1 'DO' 'BOOL' LWB = (I=L1), UPB = (I=U1); 'REF'[,]'FLUXAC' JACI = GET(JAC[I]), 'REF'[]'STATE' QRI = QR[ I ,], QGI = QG[ I ,], QFM = QF[(LWB!L1!I-1),], QFI = QF[ I ,], QFP = QF[(UPB!U1!I+1),]; 'INT' L3 = ('INT'L=2'LWB'JACI; (LWB!(L<-1!-1!L)!L)), U3 = ('INT'U=2'UPB'JACI; (UPB!(U> 1! 1!U)!U)); 'FOR' J 'FROM' L2 'TO' U2 'DO' 'BOOL' NOTL = (J>L2), NOTR = (J 1! 1!U)!U)); 'FOR' J 'FROM' JSL 'BY' JS 'TO' JSU 'DO' 'BOOL' NOTL = (J>L2), NOTR = (J1 'AND' SIGMACS>0 'THEN' 'FIELD' QC,RC; QC:= 'NUL''RESTRICT' Q; RC:= 'RESTRICT' (RHS-JAC*Q); 'TO' SIGMACS 'DO' MGCS (1,L-1,QC,RC) 'OD'; Q+:= 'PROLON' QC; (MONICS ! 'FIELD' RR = RHS-JAC*Q; MONITOR("MGCSCGC1 ",RR); MONITOR("MGCSCGC2 ",'RESTRICT'RR) ) 'FI'; LIN RELAX(QCS,L, Q, RHS); (MONICS ! MONITOR("MGCSPOST ",RHS-JAC*Q)) 'END'; #$E# #E$# # MULTIGRID (II) # 'PR' EJECT 'PR' 'INT' PMG, SIGMA, QMG; 'BOOL' FASAB:= 'FALSE', STATE2:= 'TRUE', GALERKIN:= 'FALSE'; #$H# 'PROC' FAS = ('INT' L, 'REF'[]'FIELD'Q, RHS) 'VOID': #$B FAS B$# 'BEGIN''PROC' ADDPROL = ('FIELD' QF,QC,WC)'VOID': 'BEGIN' 'REF'[,]'STATE' QQ = (Q'OF'QF)[@1,@1], QQC = (Q'OF'QC)[@1,@1], QWC = (Q'OF'WC)[@1,@1]; 'INT' N1:= 1'UPB'QQC, N2:= 2'UPB'QQC, TI,TJ,PI,PJ; 'REAL' X,Y; 'STATE' A,B,C,D,P,Q,R,S; 'IF' GALERKIN 'OR' 'ODD'N1 'OR' 'ODD'N2 'THEN''FOR' I 'TO' N1 'DO' TI:= I+I; 'FOR' J 'TO' N2 'DO' TJ:= J+J; C := QQC[I,J] - QWC[I,J]; QQ[TI-1,TJ-1]+:= C; QQ[TI-1,TJ ]+:= C; QQ[TI ,TJ-1]+:= C; QQ[TI ,TJ ]+:= C 'OD' 'OD' 'ELSE''FOR' I 'FROM' 2 'BY' 2 'TO' N1 'DO' TI:= I+I; 'FOR' J 'FROM' 2 'BY' 2 'TO' N2 'DO' TJ:= J+J; P:= QQC[I-1,J-1] - QWC[I-1,J-1]; Q:= QQC[I ,J-1] - QWC[I ,J-1]; R:= QQC[I-1,J ] - QWC[I-1,J ]; S:= QQC[I ,J ] - QWC[I ,J ]; A:= 0.25*(P+Q+R+S); B:= 0.25*(Q+S-P-R); C:= 0.25*(R+S-P-Q); D:= 0.25*(P+S-R-Q); 'FOR' K 'TO' 4 'DO' X:= K-2.5; PI:= TI-4+K; P:= X*B+A; Q := X*D+C; 'FOR' L 'TO' 4 'DO' Y:= L-2.5; PJ:= TJ-4+L; QQ[PI,PJ]+:= Y*Q+P 'OD' 'OD' 'OD' 'OD' 'FI' 'END'; (MONIREL! MONITOR("RELINA ",RHS[L]'MINRESIDU'Q[L])); RELAX(PMG,L, Q[L], RHS[L]); 'IF' L>1 'AND' SIGMA>0 'THEN' (MONIREL! MONITOR("RELOUTA ",RHS[L]'MINRESIDU'Q[L])); ( FASAB ! STATE E(Q[L]); Q[L-1]:='RESTRICT'Q[L] ); 'FIELD' W = 'COPY' Q[L-1]; RESIDUHO('FALSE',Q[L],RHS[L-1]:= 'COPY'RHS[L]); RESIDU( 'TRUE',Q[L-1],RHS[L-1]:= 'RESTRICT'RHS[L-1]); # RHS[L-1]:= 'RESTRICT'RHS[L]; 'FIELD' TAU:= 'RESIDU'Q[L-1] - 'RESTRICT''RESIDU'Q[L]; RHS[L-1]+:= TAU; # 'TO' SIGMA 'DO' FAS (L-1,Q,RHS) 'OD'; ADDPROL(STATE E(Q[L]), STATE E(Q[L-1]), STATE E(W)); # Q[L] +:= 'PROLON2' STATE E (Q[L-1]-W); # (MONIREL! MONITOR("RELINB ",RHS[L]'MINRESIDU'Q[L])) 'FI'; RELAX(QMG,L, Q[L], RHS[L]); (MONIREL! MONITOR("RELOUTB ",RHS[L]'MINRESIDU'Q[L])) 'END'; #$E# #E$# #$H# 'PROC' FMG = ('REF'[]'FIELD' SOL,RHS)'FIELD': #$B FMG B$# 'BEGIN' 'INT' LEVELS = 'UPB'SOL; ('LWB'SOL /=1 ! ERROR); ('LWB'RHS /=1 'OR' 'UPB'RHS /= LEVELS ! ERROR); 'FOR' L 'TO' LEVELS 'DO' RHS[L]:='NULRHS'SOL[L]; FAS (L, SOL, RHS); (L1 'AND' SIGMA>0 'THEN' (MONIREL! MONITOR("RELOUTA* ",RHS[L]'MINRESIDU'Q[L])); 'FIELD' W = 'COPY' ( Q[L-1]:= 'RESTRICT'Q[L] ); RESIDU('FALSE',Q[L],RHS[L-1]:= 'COPY'RHS[L]); RESIDU:= CGRESIDU; RESIDU( 'TRUE',Q[L-1],RHS[L-1]:= 'RESTRICT'RHS[L-1]); 'TO' SIGMA 'DO' FAS (L-1,Q,RHS) 'OD'; RESIDU:= STARRES1; Q[L] +:= 'PROLON' (Q[L-1]-W); (MONIREL! MONITOR("RELINB* ",RHS[L]'MINRESIDU'Q[L])) 'FI'; STARRELAX(QMG,L, Q[L], RHS[L]); (MONIREL! MONITOR("RELOUTB* ",RHS[L]'MINRESIDU'Q[L])); RESIDU:= CGRESIDU 'END'; #$E# #E$# #$H# 'PROC' STARFMG = ('REF'[]'FIELD' SOL,RHS)'FIELD': #$B STARFMG B$# 'BEGIN' 'PROC'('BOOL','FIELD','FIELD')'FIELD' CGRESIDU= RESIDU; 'INT' LEVELS = 'UPB'SOL; ('LWB'SOL /=1 ! ERROR); ('LWB'RHS /=1 'OR' 'UPB'RHS /= LEVELS ! ERROR); 'FOR' L 'TO' LEVELS 'DO' RHS[L]:='NULRHS'SOL[L]; STARFAS (L, SOL, RHS); (L0.5!0!8.0*0.042*X), QIN,QAC) ); ZUID :=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE': SUB OUT PRESS(W,'FALSE',PP,QIN,QAC); NOORD:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE': SUB IN UVZ (W,'TRUE',UU,VV,ZZ,QIN,QAC); RR := 1.0; UU:= 0.85; VV:= 0; PP:= 1/GAM; CC := SQRT(GAM*PP/RR); ZZ := LN(PP)-GAM*LN(RR); INITIATE('REPORT'Q1,'EP'([]'REAL'(RR,UU,VV,PP)),1); PRINT((NEWLINE,NEWLINE)); Q1 'END'; #$E# #E$# #$H# 'PROC' PROBLEM2 = 'FIELD': #$B PROB2 B$# 'BEGIN' PRINT(("STANDAARDPROBLEEM 2"+10*" "+"VSN. 850515", NEWLINE, "BUMP, TRANSSOON", NEWLINE, "NON-EQUIDISTANT")); 'FIELD' Q1 = COMPFIELD (1.0,1.0, -2,3, 0,2); OOST:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE': WALL (W, 'FALSE', 0, QIN,QAC) ; ZUID:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE': SUB OUT PRESS(W,'FALSE',PP,QIN,QAC); NOORD:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE': SUB IN UVZ (W,'TRUE',UU,VV,ZZ,QIN,QAC); MAPPING:= ('REAL' XX,YY) 'POINT': 'BEGIN''REAL' X,Y; ( XX<=-1.37 ! X:= 2.27*XX+2.54; X:=-0.57-0.15*(EXP(-1.650*(X+0.57))-1.0) !:XX>= 2.15 ! X:= 2.88*XX-5.64; X:= 0.55+0.11*(EXP( 1.286*(X-0.55))-1.0) ! X:= 0.32*XX-0.14 ); Y:= (YY>1.99 !2.0 !:YY<0.01 !0! 0.19*(EXP(1.222*YY)-1.0) ); 'POINT'(X, ('ABS'X<=0.5! Y+:= 0.084*(0.25-X*X)*(2.0-Y) ! Y)) 'END'; WEST:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE': WALL (W, 'TRUE' , 0, QIN,QAC) ; WALLS:= WALLSG; FLUXT:= FLUXTG; RR := GAM; UU:= 0.85; VV:= 0; PP:= 1.05; CC := SQRT(GAM*PP/RR); ZZ := LN(PP)-GAM*LN(RR); INITIATE('REPORT'Q1,'EP'([]'REAL'(RR,UU,VV,PP)),1); PRINT((NEWLINE,NEWLINE)); Q1 'END'; #$E# #E$# #$H# 'PROC' PROBLEM3 = 'FIELD': #$B PROB3 B$# 'BEGIN' PRINT(("STANDAARDPROBLEEM 3"+10*" "+"VSN. 850515", NEWLINE, "CONTACT DISCONTINUITY ", NEWLINE, "EQUIDISTANT")); 'FIELD' Q1 = COMPFIELD (1.0,1.0, 0,2, 0,2); MAPPING:= ('REAL' XX,YY) 'POINT': (XX,YY); PP:=1.0; UU1:=0.6; VV1:=-0.6; CC1:=1.0; UU2:=0.3; VV2:=-0.3; CC2:=0.6; 'REAL' RR1:=GAM*PP/(CC1*CC1), RR2:=GAM*PP/(CC2*CC2); ZZ1:=LN(PP*(RR1**(-GAM))); ZZ2:=LN(PP*(RR2**(-GAM))); RR :=(RR1+RR2)/2.0; UU:=(UU1+UU2)/2.0; VV:=(VV1+VV2)/2.0; WALLS:= WALLSR; FLUXT:= FLUXTR; OOST:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE': SUB IN UVZ (W,'FALSE',UU1,VV1,ZZ1,QIN,QAC); WEST:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE': SUB OUT PRESS(W,'TRUE',PP,QIN,QAC) ; ZUID:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE': SUB OUT PRESS(W,'FALSE',PP,QIN,QAC); NOORD:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE': SUB IN UVZ (W,'TRUE',UU2,VV2,ZZ2,QIN,QAC); INITIATE('REPORT'Q1,'EP'([]'REAL'(RR,UU,VV,PP)),1); PRINT((NEWLINE,NEWLINE)); Q1 'END'; #$E# #E$# 'REAL' PHI0; #$H# 'PROC' PROBLEM3B = 'FIELD': #$B PROB3B B$# 'BEGIN' PRINT(("STANDAARDPROBLEEM 3B"+10*" "+"VSN. 860217", NEWLINE, "CONTACT DISCONTINUITY,PHI=PI*"+FIXED(PHI0/PI,6,2), NEWLINE,"EQUIDISTANT")); 'FIELD' Q1 = COMPFIELD (1.0,1.0, 0,2, 0,2); MAPPING:= ('REAL' XX,YY) 'POINT': (XX,YY); PP:=1.0; CC1:=0.6; UU1:=0.3*SQRT(2.0)*SIN(PHI0); VV1:=0.3*SQRT(2.0)*COS(PHI0); PP:=1.0; CC2:=1.0; UU2:=0.6*SQRT(2.0)*SIN(PHI0); VV2:=0.6*SQRT(2.0)*COS(PHI0); F0:= (('REAL' X,Y)'REAL': 0.0 ); F1:= ('REAL'X,Y)'REAL':(Y>X*TAN(PHI0)!UU1!UU2)*(1+F0(X,Y)); F2:= ('REAL'X,Y)'REAL':(Y>X*TAN(PHI0)!VV1!VV2)*(1+F0(X,Y)); F3:= ('REAL'X,Y)'REAL':(Y>X*TAN(PHI0)!CC1!CC2)*(1+F0(X,Y)); F4: =('REAL'X,Y)'REAL':(Y>X*TAN(PHI0)!ZZ1!ZZ2)*(1+F0(X,Y)); 'REAL' RR1:= GAM*PP/(CC1*CC1); ZZ1:= LN(PP)-GAM*LN(RR1); 'REAL' RR2:= GAM*PP/(CC2*CC2); ZZ2:= LN(PP)-GAM*LN(RR2); RR :=(RR1+RR2)/2.0; UU:=(UU1+UU2)/2.0; VV:=(VV1+VV2)/2.0; WALLS:= WALLSR; FLUXT:= FLUXTR; OOST:= WEST:= ZUID:= NOORD:= ('WALL' W,'REAL'X,Y,'STATE'QIN,'REF''STATAC'QAC)'STATE': FIXED UVCZ (W, F1(X,Y),F2(X,Y),F3(X,Y),F4(X,Y) ,QIN,QAC); INITIATE('REPORT'Q1,'EP'([]'REAL'(RR,UU,VV,PP)),1); PRINT((NEWLINE,NEWLINE)); Q1 'END'; #$E# #E$# #$H# 'PROC' PROBLEM4 = 'FIELD': #$B PROB4 B$# 'BEGIN' PRINT(("STANDAARDPROBLEEM 4"+10*" "+"VSN. 850515", NEWLINE, "GLADDE SUBSONE STROMING", NEWLINE, "EQUIDISTANT")); 'FIELD' Q1 = COMPFIELD (0.5,0.5, -2,2, 0,2); MAPPING:= ('REAL' XX,YY) 'POINT': (XX,YY); OOST:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE': WALL (W, 'FALSE', 0, QIN,QAC) ; ZUID:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE': SUB OUT PRESS(W,'FALSE',PP,QIN,QAC); NOORD:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE': SUB IN UVZ (W,'TRUE',UU,VV,ZZ,QIN,QAC); WEST:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE': WALL (W, 'TRUE' , DIST*SIN(PI*XX), QIN,QAC); WALLS:= WALLSR; FLUXT:= FLUXTR; DIST:=0.020; RR := 1.0; UU:= 0.75; VV:= 0; PP:= 1/GAM; CC := SQRT(GAM*PP/RR); ZZ := LN(PP)-GAM*LN(RR); INITIATE('REPORT'Q1,'EP'([]'REAL'(RR,UU,VV,PP)),1); PRINT((NEWLINE,NEWLINE)); Q1 'END'; #$E# #E$# #$H# 'PROC' PROBLEM5 = 'FIELD': #$B PROB5 B$# 'BEGIN' PRINT(("STANDAARDPROBLEEM 5"+10*" "+"VSN. 850731", NEWLINE, "REFLECTERENDE SCHEVE SCHOKGOLF IN KANAAL", NEWLINE, "EQUIDISTANT ALS BIJ OSHER")); 'FIELD' Q1 = COMPFIELD (4/6,1/2, 0,6, 0,2); MAPPING:= ('REAL' XX,YY) 'POINT': (XX,YY); NOORD:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE': FIXED UVCZ (W,UU1,VV1,CC1,ZZ1,QIN,QAC); WEST:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE': WALL (W, 'TRUE' , 0, QIN,QAC) ; OOST:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE': FIXED UVCZ (W,UU2,VV2,CC2,ZZ2,QIN,QAC); ZUID:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE': SUB OUT PRESS (W,'FALSE',PP,QIN,QAC); WALLS:= WALLSR; FLUXT:= FLUXTR; # TESTPROBLEEM: OBLIQUE SHOCK, EXACTE OPLOSSING: 'REAL' U1:=2.9,V1:=0.0,C1:=1.0,Z1:=-0.47,P1:=1.0,R1:=1.4,S1:=0.62; 'REAL' U2:=2.62,V2:=-0.51,C2:=1.12,Z2:=-0.45,P2:=2.14,R2:=2.38, S2:=0.64; 'REAL' U3:=2.41,V3:=0.0,C3:=1.23,Z3:=-0.45,P3:=4.0,R3:=3.70,S3:=0.64; # UU1:=2.9 ; VV1:= 0.0 ; CC1:=1.0 ; ZZ1:=-0.47; UU2:=2.62; VV2:=-0.51; CC2:=1.12; ZZ2:=-0.45; PP :=4.0 ; UU:=(UU1+UU2)/2.0; VV:=(VV1+VV2)/2.0; CC:=(CC1+CC2)/2.0; ZZ:=(ZZ1+ZZ2)/2.0; INITIATE('REPORT'Q1,'EZ''STATE'(UU,VV,CC,ZZ),1); PRINT((NEWLINE,NEWLINE)); Q1 'END'; #$E# #E$# #?$H? 'PROC' PROBLEM6 = 'FIELD': ?$B PROB5 B$? ? NOG NIET GEBRUIKT TESTPROBLEEM, 850513 ? 'BEGIN' PRINT(("STANDAARDPROBLEEM 6"+10*" "+"VSN. 850515", NEWLINE, "REFLECTERENDE SCHEVE SCHOKGOLF IN KANAAL", NEWLINE, "NIET-EQUIDISTANTE INLAAT, EQUIDISTANTE UITLAAT")); 'FIELD' Q1 = COMPFIELD (1.0,1.0, 0,8, 0,1); MAPPING:= ('REAL' XX,YY) 'POINT': ( XX, ( XX<1.0 ! YY*(1.0+0.525*(1.0-XX)) ! YY ) ); NOORD:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE': FIXED UVCZ (W,UU,VV,CC,ZZ,QIN,QAC); WEST:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE': WALL (W, 'TRUE' , 0, QIN,QAC) ; OOST:=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE': WALL (W, 'FALSE', 0, QIN,QAC) ; WALLS:= WALLSG; FLUXT:= FLUXTG; RR:= 1.0; UU:= 2.0; VV:= 0.0; PP:= 1.0/GAM; CC:= SQRT(GAM*PP/RR); ZZ:= LN(PP)-GAM*LN(RR); INITIATE('REPORT'Q1,'EP'([]'REAL'(RR,UU,VV,PP)),1); PRINT((NEWLINE,NEWLINE)); Q1 'END';?$E? ?E$? # # RESULTAAT ANALYSE # 'PR' EJECT 'PR' #$H# 'PROC' PRINTORDER = ([]'INT' NOS )'VOID': #$B PRNT.OR B$# 'BEGIN' ('UPB' NOS /=3 ! ERROR); [1:3]'FIELD' ALLF; PRINT((NEWLINE,"OBSERVED ORDER OF ACCURACY")); 'FOR' L 'TO' 3 'DO' TAKE( NOS[L], ALLF[L] ) 'OD'; ALLF[3]:= 'RESTRICT' ALLF[3]; ALLF[2]-:= ALLF[3]; ALLF[1]-:= 'RESTRICT' ALLF[3]; ALLF[1]-:= 'RESTRICT' ALLF[2]; [1:4 ]'REAL' DEN:= 'NORM1'ALLF[2], NUM:= 'NORM1'ALLF[1]; PRINT((NEWLINE,"ORDER=")); 'FOR' I 'TO'4 'DO' 'REAL' FACTOR := NUM[I]/DEN[I]; PRINT((" "+WHOLE(I,0)+" ("+FIXED(FACTOR,6,2)+ ") "+(FACTOR>0!FIXED(LN(FACTOR)/LN(2.0),6,2)!"******") )) 'OD' 'END'; #$E# #E$# #$H# 'PROC' TOTVAR = ('FIELD' QF )[,]'REAL': #$B TOTVAR B$# 'BEGIN' 'REF'[,]'STATE' Q = Q'OF'QF; 'REAL' DDX = DDX'OF'QF, DDY = DDY'OF'QF; 'INT' IL=1'LWB'Q,IU=1'UPB'Q,JL=2'LWB'Q,JU=2'UPB'Q; [1:4,0:1]'REAL' NRM; 'FIELD' QR = (DDX,DDY,'LOC''INT':=1,Q[ IL+1:IU'AT'1,]), QL = (DDX,DDY,'LOC''INT':=1,Q[ IL:IU-1'AT'1,]), QU = (DDX,DDY,'LOC''INT':=1,Q[,JL+1:JU'AT'1 ]), QD = (DDX,DDY,'LOC''INT':=1,Q[,JL:JU-1'AT'1 ]); [,]'REAL' NRMH = 'NORM' (QR-QL), NRMV = 'NORM' (QU-QD); 'FOR' C 'TO' 4 'DO' NRM[C,1]:= (NRMH[C,1]+NRMV[C,1])* 0.5*SQRT((IU-IL+1)*(JU-JL+1)) ; NRM[C,0]:= (NRMH[C,0]S ! S:=T ) 'OD' 'OD'; (K!C1'OF'NORM ,C2'OF'NORM ,C3'OF'NORM ,C4'OF'NORM ) := S1; (K!C1'OF'RATIO,C2'OF'RATIO,C3'OF'RATIO,C4'OF'RATIO) := S*N1*N2/S1 'OD';'ABS' NORM 'END';#$E# #E$# #$H# 'PROC' ACCOUNT OK = 'BOOL': #$B ACCOU B$# 'BEGIN' 'REAL' E,D,C:= CLOCK; D:=C-TIME; E:=TIMELIMIT-C; 'BOOL' THROUGH = E > 1.2*D; TIME:= C; THROUGH 'END'; #$E# #E$# #$H# 'PROC' HISTO = ('PROC'('STATE')'REAL' P, 'REF'[]'STATE' S)'VOID': #$B HISTO B$# 'BEGIN' 'INT' L= 'LWB'S, U='UPB'S; [L:U]'REAL'X; 'FOR'I'FROM'L'TO'U 'DO'X[I]:=P(S[I])'OD'; 'REAL' FACTOR, MAX:= 0, MIN:= 0; 'INT' K, W:= (U-L<30!40!:U-L<60!80!120); 'FOR'I'FROM'L'TO'U 'DO'(X[I]>MAX!MAX:=X[I]!: X[I]>> MONITOR STILL UNDEFINED <<<", NEWLINE,">>> ADVISE : DO MONITOR := MONIT1 <<<", NEWLINE)); ERROR ; 'SKIP' 'END' ; #$H# 'PROC' MONIT1 = ('STRING' TEXT, 'FIELD' R)'VOID': #$B MONIT1 B$# 'BEGIN' 'FLUX' NORM,RATIO; FLUXNORM(R,NORM,RATIO); 'REAL' S = ( C1'OF'NORM + C2'OF'NORM + C3'OF'NORM + C4'OF'NORM ) /4; PRINT((NEWLINE, TEXT+" ", ST([]'REAL'((DDX'OF'R,DDY'OF'R))),SPACE, FLOAT(S,10,3,3)+" < "+ SF(NORM)+ST(RATIO) )) 'END'; #$E# #E$# 'BOOL' MONIREL:= 'FALSE', MONICS := 'FALSE', MONINWT:= 'FALSE', MONILIN:= 'FALSE', MONIFAS:= 'FALSE', MONI := 'FALSE', MONIGS := 'FALSE'; 'REAL' TIME:= CLOCK; PRINT((">> EULER A68-LIBRARY VSN.860414 ",NEWLINE, ">> FLUX ETC. ARE STRUCTS (JR) ",NEWLINE, NEWLINE,NEWLINE)) ; 'PR' PROG 'PR' 'SKIP' ; (CAVITATION! PRINT((NEWLINE,"CAVITATION OCCURRED",NEWLINE)) ) 'END' LIBRARY(EULIB,NEW) ADD(*,LGO) FINISH. ENDRUN. LIBRARY(EULIB,OLD) ADD(*,LGO) FINISH. ENDRUN. 'BEGIN' # EUTEST VSN.850724 PWH/840305 # # REFERENTIE OPLOSSING, EERSTE ORDE, PROB. 4 # # GLAD TESTPROBLEEM # RESIDU := RESIDU1 ; RESIDU2 := CENTRALX ; NFLUX := OSHER ; SOS:=-1 ; GALERKIN := 'TRUE' ; SYMMETRIC := 'TRUE' ; FASAB := 'FALSE' ; PMG:=1;QMG:=1;SIGMA:=1 ; RELAX := SEIDEL ; NLGSTOL := 0.1 ; 'PROC' SHOW = ('INT' L)'VOID': 'IF' L=LEVELS 'THEN' 'FIELD' RR:= Q[L]; 'FLUX' NORM,RATIO; 'REAL' N; PRINT((NEWLINE,NEWLINE, "OPLOSSING OP NIVEAU ",WHOLE(L,0) )); 'PRINT' RR; RR:= ADD WALL STATES('FALSE',Q[L]); STATE E (RR); PRINT((NEWLINE,NEWLINE,"PRESSURE",NEWLINE)); PRESSURE'PRINT'RR; Q'OF'RR:= 'NIL'; RR := 'RESIDU' Q[L]; FLUXNORM(RR,NORM,RATIO); PRINT((NEWLINE,NEWLINE, "RESIDU1 OP NIVEAU ",WHOLE(L,0), " ", N:= 'ABS'NORM)); PRINT((" FACTOR=10**-"+WHOLE(N:= 'ROUND'(LN(N)/LN(10.0)),0))); 'PRINT' (RR*:= 10.0**(1-N)); RR := 'RESIDU2' Q[L]; FLUXNORM(RR,NORM,RATIO); PRINT((NEWLINE,NEWLINE, "RESIDU2 OP NIVEAU ",WHOLE(L,0), " ", N:= 'ABS'NORM)); PRINT((" FACTOR=10**-"+WHOLE(N:= 'ROUND'(LN(N)/LN(10.0)),0))); 'PRINT' (RR*:= 10.0**(1-N)); RR:= Q[L] 'CORRECT''NUL'R[L]; FLUXNORM(RR,NORM,RATIO); PRINT((NEWLINE,NEWLINE, "CORRECTIE OP NIVEAU ",WHOLE(L,0), " ", N:= 'ABS'NORM)); PRINT((" FACTOR=10**-"+WHOLE(N:= 'ROUND'(LN(N)/LN(10.0)),0))); 'PRINT' (RR*:= 10.0**(1-N)); RR:= 'RESTRICT' RR; FLUXNORM(RR,NORM,RATIO); PRINT((NEWLINE,NEWLINE, "RESTRICTBAR VAN CORRECTIE OP NIVEAU ",WHOLE(L,0), " ", N:= 'ABS'NORM)); PRINT((" FACTOR=10**-"+WHOLE(N:= 'ROUND'(LN(N)/LN(10.0)),0))); 'PRINT' (RR*:= 10.0**(1-N)); (L>1 ! RR := 'TAU1' Q[L]; FLUXNORM(RR,NORM,RATIO); PRINT((NEWLINE,NEWLINE, "TAU1 OP NIVEAU ",WHOLE(L,0), " ", N:= 'ABS'NORM)); PRINT((" FACTOR=10**-"+WHOLE(N:= 'ROUND'(LN(N)/LN(10.0)),0))); 'PRINT' (RR*:= 10.0**(1-N)); RR:= 'TAU2' Q[L]; FLUXNORM(RR,NORM,RATIO); PRINT((NEWLINE,NEWLINE, "TAU2 OP NIVEAU ",WHOLE(L,0), " ", N:= 'ABS'NORM)); PRINT((" FACTOR=10**-"+WHOLE(N:= 'ROUND'(LN(N)/LN(10.0)),0))); 'PRINT' (RR*:= 10.0**(1-N)); RR := R[L-1]; FLUXNORM(RR,NORM,RATIO); PRINT((NEWLINE,NEWLINE, "RHS OP NIVEAU ",WHOLE(L-1,0), " ", N:= 'ABS'NORM)); PRINT((" FACTOR=10**-"+WHOLE(N:= 'ROUND'(LN(N)/LN(10.0)),0))); 'PRINT' (RR*:= 10.0**(1-N)) ) 'FI'; 'INT' LEVELS = 3; [1:LEVELS]'FIELD'Q,R; Q[1]:= PROBLEM4; 'FOR' L 'TO'LEVELS 'DO' R[L]:= 'NULRHS'Q[L]; 'TO' 4 'WHILE' MONITOR("FAS ",'RESIDU'(Q[L])); ACCOUNT OK 'DO' FAS(L, Q, R) 'OD'; SHOW(L); (L