Forgot your password?
Document Actions
P. W. Hemker. A library of Euler multigrid routines.
P. W. Hemker. A library of Euler multigrid routines. Centrum voor Wiskunde en Informatica, Numerieka Wiskunde, 1986.
Size 167.4 kB - File type text/plain
Click here to get the file
File contents
#********************************************************************#
# #
# 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<U2),
'REF'[] 'FLUXAC' JACIJ = JACI[J,],
'REF' 'STATE' QGIJ = QGI[J]; QGIJ:=NUL ;
'FOR' K 'FROM' L3 'TO' U3
'DO' 'CASE' K+3
'IN' QGIJ+:= JACIJ[K]*QFM[J ] ,
(NOTL ! QGIJ+:= JACIJ[K]*QFI[J-1]),
QGIJ+:= JACIJ[K]*QFI[J ] ,
(NOTR ! QGIJ+:= JACIJ[K]*QFI[J+1]),
QGIJ+:= JACIJ[K]*QFP[J ]
'OUT' ERROR
'ESAC'
'OD'
'OD' ;
PET('FALSE',JAC[I])
'OD' ; G
'END'; #$E#
#E$#
# FIELD MANIPULATION # 'PR' EJECT 'PR'
#$H#
'PROC' INITIATE = ('FIELD' F, 'STATE' A, 'INT' TYPE )'FIELD':
#$B INI.FLD B$#
'BEGIN' 'REF'[,]'STATE' Q = (Q'OF'F)[@1,@1];
TYPE'OF'F:=TYPE ; 'INT' N1=1'UPB'Q,N2=2'UPB'Q;
'FOR' I 'TO' N1 'DO' 'REF'[]'STATE' FI=Q[I,];
'FOR' J 'TO' N2 'DO' FI[J] := A
'OD' 'OD'; F
'END'; #$E#
#E$#
'PROC'('REAL','REAL')'REAL' F0,F1,F2,F3,F4;
#$H#
'PROC' SETFIELD = ('FIELD' QQ, 'INT' TYPE, 'BOOL' STAR,
'PROC'('REAL','REAL')'REAL' F1,F2,F3,F4)'FIELD':
#$B SET.FLD B$#
'BEGIN'
'REAL' DDX = DDX'OF'QQ, DDY = DDY'OF'QQ ; TYPE'OF'QQ:=TYPE ;
'POINT' P; 'REF''REAL' X=X'OF'P, Y=Y'OF'P ; 'REAL' XX, YY ;
'PROC' HERE = ('REAL' XX,YY)'STATE':
( P:= MAPPING (XX,YY);
'STATE'(F1(X,Y),F2(X,Y),F3(X,Y),F4(X,Y)) );
'REF'[,]'STATE' Q = Q'OF'QQ;
'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;
(STAR
! 'IF' 'NOT' ( ILH*2 = IL+1 'AND' IUH*2 = IU 'AND'
JLH*2 = JL+1 'AND' JUH*2 = JU )
'THEN' ERROR; QQ
'ELSE' 'REAL' X23 = 2*DDX/3, Y23 = 2*DDY/3;
'FOR' IH 'FROM' ILH 'TO' IUH
'DO' 'INT'I = 2*IH-1; XX:= I*DDX ;
'REF'[]'STATE' QI=Q[I,],QIP=Q[I+1,];
'FOR' JH 'FROM' JLH 'TO' JUH
'DO' 'INT'J = 2*JH-1; YY:= J*DDY; 'INT'JP=J+1;
QI [J ] := HERE ( XX , YY-Y23 ) ;
QI [JP] := HERE ( XX-X23 , YY ) ;
QIP [JP] := HERE ( XX , YY+Y23 ) ;
QIP [J ] := HERE ( XX+X23 , YY )
'OD''OD'; QQ
'FI'
! 'FOR' I 'FROM' IL 'TO' IU
'DO' XX:= (I-0.5)*DDX; 'REF'[]'STATE'QI=Q[I,] ;
'FOR' J 'FROM' JL 'TO' JU
'DO' QI[J]:= HERE( XX,(J-0.5)*DDY )
'OD' 'OD'; QQ
) 'END' ; #$E#
#E$#
#$H#
'OP' 'RESTRICT' = ('FIELD' A)'FIELD':
#$B RESTR B$#
'BEGIN' 'REAL' W = ( TYPE'OF'A ! 0.25,0.25 ! 1.0);
'REF'[,]'STATE' Q = Q'OF'A;
'INT' L2= 1'LWB'Q, U2= 1'UPB'Q, L3= 2'LWB'Q, U3= 2'UPB'Q;
'INT' LL2:= (L2+1)'OVER'2, UU2:= U2'OVER'2, TI, TJ,
LL3:= (L3+1)'OVER'2, UU3:= U3'OVER'2;
'HEAP'[LL2:UU2,LL3:UU3]'STATE' QN;
'FOR' I 'FROM' LL2 'TO' UU2 'DO' TI:= I+I;
'REF'[]'STATE' QNI=QN[I,],QTIM=Q[TI-1,],QTI=Q[TI,];
'FOR' J 'FROM' LL3 'TO' UU3 'DO' TJ:= J+J;
QNI[J]:= W* (QTIM[TJ-1]+QTIM[TJ]+QTI[TJ-1]+QTI[TJ])
'OD' 'OD';
'FIELD'(2*(DDX'OF'A),2*(DDY'OF'A),'HEAP''INT':=TYPE'OF'A,QN)
'END'; #$E#
#E$#
#$H#
'OP' 'PROLON' = ('FIELD'A)'FIELD':
#$B PROL.1 B$#
'BEGIN' 'REF'[,]'STATE' Q = Q'OF'A;
'INT' L2= 1'LWB'Q, U2= 1'UPB'Q, L3= 2'LWB'Q, U3= 2'UPB'Q;
'INT' TI, TJ;
'HEAP'[2*L2-1:2*U2,2*L3-1:2*U3]'STATE' QN;
'FOR' I 'FROM' L2 'TO' U2 'DO' TI:= I+I;
'REF'[]'STATE' QNTIM=QN[TI-1,],QNTI=QN[TI,],QI=Q[I,];
'FOR' J 'FROM' L3 'TO' U3 'DO' TJ:= J+J;
QNTIM[TJ-1]:=QNTIM[TJ]:=QNTI[TJ-1]:=QNTI[TJ]:=QI[J]
'OD' 'OD';
'FIELD'((DDX'OF'A)/2,(DDY'OF'A)/2,'HEAP''INT':=TYPE'OF'A,QN)
'END'; #$E#
#E$#
'PR' EJECT 'PR'
#$H#
'OP' 'PROLON2' = ('FIELD' F)'FIELD':
#$B PROLO2 B$#
'BEGIN' 'REF'[,]'STATE' QQ = (Q'OF' F);
'INT' L2= 1'LWB'QQ, U2= 1'UPB'QQ,
L3= 2'LWB'QQ, U3= 2'UPB'QQ;
'IF' 'ODD'(U2-L2+1) 'OR' 'ODD'(U3-L3+1)
'THEN' 'PROLON' F
'ELSE' 'HEAP'[2*L2-1:2*U2,2*L3-1:2*U3]'STATE'QQQ;
'INT' TI, TJ, PI, PJ;
'REAL' X,Y;
'STATE' A,B,C,D,P,Q,R,S;
'FOR' I 'FROM' L2+1 'BY' 2 'TO' U2 'DO' TI:= I+I;
'REF'[]'STATE' QQIM=QQ[I-1,],QQI=QQ[I,] ;
'FOR' J 'FROM' L3+1 'BY' 2 'TO' U3 'DO' TJ:= J+J;
A:= QQIM[J-1] ; B := QQI[J-1]; C := QQIM[J];
D:= QQI [J]-B-C+A; B-:= A ; C-:= A ; X:=-0.75 ;
'FOR' K 'TO' 4 'DO' P:=C+(X+:=0.5)*D ; Q:=A+X*B ;
Y:=-0.75;'REF'[]'STATE'QQQK=QQQ[TI-4+K,];
'FOR' L 'TO' 4 'DO' QQQK[TJ-4+L] := Q+(Y+:=0.5)*P
'OD' 'OD'
'OD' 'OD';
'FIELD'((DDX'OF'F)/2,(DDY'OF'F)/2,'HEAP''INT':=TYPE'OF'F,
QQQ)
'FI'
'END'; #$E#
#E$#
#$H#
'OP' 'PROLIN2' = ('FIELD' F)'FIELD':
#$B PROLI2 B$#
'BEGIN' 'REF'[,]'STATE' QQ = (Q'OF' F);
'INT' L1= 1'LWB'QQ, U1= 1'UPB'QQ,
L2= 2'LWB'QQ, U2= 2'UPB'QQ;
'IF' (L1=U1)'OR'(L2=U2)
'THEN' 'PROLON'F
'ELSE' 'HEAP'[2*L1-1:2*U1,2*L2-1:2*U2]'STATE'QQQ;
'INT' TI, TJ, IP, JP,
'STATE' P,Q,R,S,
'BOOL' LWBI, UPBI;
# 25 15 5 -5 #
# 15 9 3 -3 # # /16 #
# 5 3 1 -1 #
# -5 -3 -1 1 #
'REAL' Z25 = 25/16, Z15 = 15/16, Z9 = 9/16,
Z5 = 5/16, Z3 = 3/16, Z1 = 1/16;
'FOR' I 'FROM' L1+1 'TO' U1
'DO' IP:=2*I-2;
LWBI:=(I=L1+1);
UPBI:=(I=U1 );
'FOR' J 'FROM' L2+1 'TO' U2
'DO' JP:= 2*J-2;
P:= QQ[I-1,J-1]; Q:= QQ[I ,J-1];
R:= QQ[I-1,J ]; S:= QQ[I ,J ];
QQQ[IP ,JP ]:=Z9*P+Z3*(Q+R)+Z1*S;
QQQ[IP+1,JP ]:=Z9*Q+Z3*(P+S)+Z1*R;
QQQ[IP ,JP+1]:=Z9*R+Z3*(P+S)+Z1*Q;
QQQ[IP+1,JP+1]:=Z9*S+Z3*(Q+R)+Z1*P;
(J=L2+1!(LWBI ! QQQ[IP-1,JP-1]:=
Z25*P-Z5*(Q+R)+Z1*S);
(UPBI ! QQQ[IP+2,JP-1]:=
Z25*Q-Z5*(P+S)+Z1*R);
QQQ[IP ,JP-1]:=
Z15*P+Z5*Q-Z3*R-Z1*S;
QQQ[IP+1,JP-1]:=
Z15*Q+Z5*P-Z3*S-Z1*R );
(J=U2 ! (LWBI ! QQQ[IP-1,JP+2]:=
Z25*R-Z5*(P+S)+Z1*Q);
(UPBI ! QQQ[IP+2,JP+2]:=
Z25*S-Z5*(Q+R)+Z1*P);
QQQ[IP ,JP+2]:=
Z15*R+Z5*S-Z3*P-Z1*Q;
QQQ[IP+1,JP+2]:=
Z15*S+Z5*R-Z3*Q-Z1*P );
(LWBI ! QQQ[IP-1,JP ]:=
Z15*P+Z5*R-Z3*Q-Z1*S;
QQQ[IP-1,JP+1]:=
Z15*R+Z5*P-Z3*S-Z1*Q );
(UPBI ! QQQ[IP+2,JP ]:=
Z15*Q+Z5*S-Z3*P-Z1*R;
QQQ[IP+2,JP+1]:=
Z15*S+Z5*Q-Z3*R-Z1*P )
'OD' 'OD';
'FIELD'((DDX'OF'F)/2,(DDY'OF'F)/2,'HEAP''INT':=TYPE'OF'F,
QQQ)
'FI'
'END'; #$E#
#E$#
#$H#
'OP' 'CLLPNT' = ('FIELD' F)'FIELD':
#$B CLLPNT B$#
'BEGIN' 'FIELD' WF:= ADD WALL('TRUE',F);
'REF'[,]'STATE' QF = (Q'OF' WF);
'INT' L1= 1'LWB'QF, U1= 1'UPB'QF,
L2= 2'LWB'QF, U2= 2'UPB'QF;
'INT' U1M = U1-1, U2M = U2-1;
'HEAP'[L1:U1M,L2:U2M]'STATE' QPF;
'FOR' I 'FROM' L1 'TO' U1M
'DO' 'FOR' J 'FROM' L2 'TO' U2M
'DO' QPF[I,J]:= 0.25*( QF[I ,J ]+QF[I ,J+1]+
QF[I+1,J ]+QF[I+1,J+1] )
'OD'
'OD' ;
Q'OF'WF:='NIL';
'FIELD'( DDX'OF'F, DDY'OF'F, 'HEAP''INT':=TYPE'OF'F, QPF)
'END' ; #$E#
#E$#
#$H#
'OP' 'PROSCND' = ('FIELD' C) 'FIELD':
#$B PROSCND B$#
'IF' 'REF'[,]'STATE' QC = (Q'OF' C);
'ODD'(1'UPB'QC-1'LWB'QC+1)'OR''ODD'(2'UPB'QC-2'LWB'QC+1)
'THEN' 'PROLON'C
'ELSE' 'FIELD' PC:= 'CLLPNT'C;
'REF'[,]'STATE' QP = (Q'OF'PC);
'INT' L1= 1'LWB'QP, U1= 1'UPB'QP,
L2= 2'LWB'QP, U2= 2'UPB'QP;
'HEAP'[2*(L1+1)-1:2*U1,2*(L2+1)-1:2*U2]'STATE'QF;
'STATE' S11, S12, S13, S21, S22, S23, S31, S32, S33;
'REAL' Z1:=1/16, Z2:=1/8, Z3 :=3/16, Z4 := 1/4, Z5 := 5/16,
Z6:=3/8, Z9:=9/16, Z10:=5/8, Z15:=15/16, Z25:=25/16;
( TYPE'OF'C = 0 ! Z1 *:=0.25; Z2 *:=0.25; Z3 *:=0.25;
Z4 *:=0.25; Z5 *:=0.25; Z6 *:=0.25;
Z9 *:=0.25; Z10*:=0.25; Z15*:=0.25;
Z25*:=0.25 );
'INT' I1:= 2*(L1+1)-1;
'FOR' I 'FROM' L1 'BY' 2 'TO' (U1-2)
'DO' 'INT' I2 = I1+1, I3 = I1+2, I4 = I1+3;
S13:= QP[I ,L2]; S23:= QP[I+1,L2]; S33:= QP[I+2,L2];
'INT' J1:= 2*(L2+1)-1;
'FOR' J 'FROM' L2 'BY' 2 'TO' (U2-2)
'DO' 'INT' J2 = J1+1, J3 = J1+2, J4 = J1+3;
S11:= S13; S12:= QP[I ,J+1]; S13:= QP[I ,J+2];
S21:= S23; S22:= QP[I+1,J+1]; S23:= QP[I+1,J+2];
S31:= S33; S32:= QP[I+2,J+1]; S33:= QP[I+2,J+2];
QF[I1,J1]:= Z4 *S11 + Z6 *S12 - Z2 *S13
+ Z6 *S21 + Z9 *S22 - Z3 *S23
- Z2 *S31 - Z3 *S32 + Z1 *S33;
QF[I1,J2]:= Z10*S12 - Z2 *S13
+ Z15*S22 - Z3 *S23
- Z5 *S32 + Z1 *S33;
QF[I1,J3]:= - Z2 *S11 + Z10*S12
- Z3 *S21 + Z15*S22
+ Z1 *S31 - Z5 *S32;
QF[I1,J4]:= - Z2 *S11 + Z6 *S12 + Z4 *S13
- Z3 *S21 + Z9 *S22 + Z6 *S23
+ Z1 *S31 - Z3 *S32 - Z2 *S33;
QF[I2,J1]:= Z10*S21 + Z15*S22 - Z5 *S23
- Z2 *S31 - Z3 *S32 + Z1 *S33;
QF[I2,J2]:= Z25*S22 - Z5 *S23
- Z5 *S32 + Z1 *S33;
QF[I2,J3]:= - Z5 *S21 + Z25*S22
+ Z1 *S31 - Z5 *S32;
QF[I2,J4]:= - Z5 *S21 + Z15*S22 + Z10*S23
+ Z1 *S31 - Z3 *S32 - Z2 *S33;
QF[I3,J1]:= - Z2 *S11 - Z3 *S12 + Z1 *S13
+ Z10*S21 + Z15*S22 - Z5 *S23;
QF[I3,J2]:= - Z5 *S12 + Z1 *S13
+ Z25*S22 - Z5 *S23;
QF[I3,J3]:= Z1 *S11 - Z5 *S12
- Z5 *S21 + Z25*S22;
QF[I3,J4]:= Z1 *S11 - Z3 *S12 - Z2 *S13
- Z5 *S21 + Z15*S22 + Z10*S23;
QF[I4,J1]:= - Z2 *S11 - Z3 *S12 + Z1 *S13
+ Z6 *S21 + Z9 *S22 - Z3 *S23
+ Z4 *S31 + Z6 *S32 - Z2 *S33;
QF[I4,J2]:= - Z5 *S12 + Z1 *S13
+ Z15*S22 - Z3 *S23
+ Z10*S32 - Z2 *S33;
QF[I4,J3]:= Z1 *S11 - Z5 *S12
- Z3 *S21 + Z15*S22
- Z2 *S31 + Z10*S32;
QF[I4,J4]:= Z1 *S11 - Z3 *S12 - Z2 *S13
- Z3 *S21 + Z9 *S22 + Z6 *S23
- Z2 *S31 + Z6 *S32 + Z4 *S33;
J1+:= 4
'OD' ;
I1+:= 4
'OD' ;
Q'OF'PC:= 'NIL';
'FIELD'((DDX'OF'C)/2,(DDY'OF'C)/2,'HEAP''INT':=TYPE'OF'C,QF)
'FI' ; #$E#
#E$#
# FIRST-ORDER STAR-INTERPOLATION # 'PR' EJECT 'PR'
#$H#
'OP' 'STARPROL' = ('FIELD' F)'FIELD':
#$B STRPRL B$#
'IF' 'REF'[,]'STATE' QC = Q'OF'F;
'INT' IL= 1'LWB'QC, IU= 1'UPB'QC,
JL= 2'LWB'QC, JU= 2'UPB'QC;
'ODD'(IU-IL+1) 'OR' 'ODD'(JU-JL+1)
'THEN' 'PROLON' F
'ELSE' 'HEAP'[2*IL-1:2*IU,2*JL-1:2*JU]'STATE'QF;
'INT' TI, TJ;
'FOR' I 'FROM' IL+1 'BY' 2 'TO' IU 'DO' TI:= I+I;
'FOR' J 'FROM' JL+1 'BY' 2 'TO' JU 'DO' TJ:= J+J;
QF[TI-3,TJ-2]:= QF[TI-2,TJ-2]:=
QF[TI-3,TJ-1]:= QF[TI-3,TJ ]:= QC[I-1,J ];
QF[TI-2,TJ ]:= QF[TI-2,TJ-1]:=
QF[TI-1,TJ ]:= QF[TI ,TJ ]:= QC[I ,J ];
QF[TI ,TJ-1]:= QF[TI-1,TJ-1]:=
QF[TI ,TJ-2]:= QF[TI ,TJ-3]:= QC[I ,J-1];
QF[TI-1,TJ-3]:= QF[TI-1,TJ-2]:=
QF[TI-2,TJ-3]:= QF[TI-3,TJ-3]:= QC[I-1,J-1]
'OD' 'OD';
'FIELD'((DDX'OF'F)/2,(DDY'OF'F)/2,'HEAP''INT':=TYPE'OF'F,
QF)
'FI'; #$E#
#E$#
# STATE FUNCTIONS # 'PR' EJECT 'PR'
# THESE FUNCTIONS ASSUME STATE (1) #
'PROC' DENSITY = ('STATE' Q)'REAL' : C1'OF'Q ;
'PROC' TOTENERGY = ('STATE' Q)'REAL' : C4'OF'Q ;
'PROC' XVELOCITY = ('STATE' Q)'REAL' : C2'OF'Q/C1'OF'Q ;
'PROC' YVELOCITY = ('STATE' Q)'REAL' : C3'OF'Q/C1'OF'Q ;
'PROC' XMOMENT = ('STATE' Q)'REAL' : C2'OF'Q ;
'PROC' YMOMENT = ('STATE' Q)'REAL' : C3'OF'Q ;
'PROC' LOCMACH = ('STATE' Q)'REAL': SPEED(Q)/SOUND(Q);
'PROC' PRESSURE = ('STATE' Q)'REAL':
('REAL' R=C1'OF'Q; 'REAL' RU=C2'OF'Q, RV=C3'OF'Q;
((RU*RU+RV*RV)/(-2*R) + C4'OF'Q)*GAMM1 );
'PROC' DIREC = ('STATE' Q)'REAL':
('REAL' U=C2'OF'Q, V= C3'OF'Q;
'REAL' D=(U=0!'SIGN'V!2*ARCTAN(V/U)/PI);
U<0! ( V<0 ! D-2 ! D+2) ! D );
'PROC' ENTROPY = ('STATE' Q)'REAL':
('REAL' R=C1'OF'Q, RU=C2'OF'Q, RV=C3'OF'Q;
((RU*RU+RV*RV)/(-2*R)+ C4'OF'Q)*GAMM1* R**-GAM);
'PROC' TOTENTHALPY = ('STATE' Q)'REAL':
('REAL' R=C1'OF'Q, RU=C2'OF'Q, RV=C3'OF'Q;
((RU*RU+RV*RV)*GAMM1/(-2*R)+ GAM*C4'OF'Q) /R );
'PROC' SOUND = ('STATE' Q)'REAL': SQRT(GAM*GAMM1*TEMPERAT(Q));
'PROC' SPEED = ('STATE' Q)'REAL':
SQRT(C2'OF'Q*C2'OF'Q+C3'OF'Q*C3'OF'Q)/C1'OF'Q;
'PROC' TEMPERAT= ('STATE' Q)'REAL':
('REAL' R=C1'OF'Q, RU=C2'OF'Q, RV=C3'OF'Q;
# I.E. P/(R*GAMM1) # ((RU*RU+RV*RV)/(-2*R) + C4'OF'Q)/R );
[]'PROC'('STATE')'REAL' VAR =
(DENSITY, XMOMENT, YMOMENT, TOTENERGY,
XVELOCITY, YVELOCITY, SOUND, ENTROPY,
PRESSURE, SPEED, LOCMACH, DIREC,
TOTENTHALPY, TEMPERAT);
'MODE' 'VAR' = 'STRUCT'('BOOL' INTERN, 'REAL' DDX,DDY,
'INT' VAR, 'REF'[,]'REAL'Q);
#$H#
'PROC' VARIABLE = ('STRING' NAME, 'BOOL' INTERN, 'FIELD' F) 'VAR':
#$B VAR B$#
'BEGIN' 'INT' VARNO:= 0;
( NAME > "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) )) );
RESTOL<RNEW 'AND' RNEW<=ROLD 'AND' ACCOUNT OK
'DO' 'FIELD' W:= 'NUL''COPY'QF, RES; ROLD:= RNEW;
'REAL' OLDNORM:= RNEW, ONORM:= RNEW, NORM:= RNEW;
(MONINWT ! PRINT((NEWLINE,"LINSOLVE"+21*" "+
FLOAT(OLDNORM,10,3,3)+" "+
FIXED(-4*LN(OLDNORM)/LN(10.0),6,1) )));
'INT' CM;
'FOR' C 'TO' MLINSOLVE
'WHILE' LTOL<NORM 'AND' NORM<=ONORM
'DO' CM:=C; ONORM:= NORM;
RES:= 'COPY' RHS;
( C>1 ! 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<U2),
'REF'[]'FLUXAC' JACIJ = JACI[J,],
'REF''STATE' QGIJ = QGI[J];
'FOR' N 'FROM' L3 'TO' U3
'DO' 'CASE' N+3
'IN' QGIJ-:= JACIJ[N]*QFM[J ] ,
(NOTL ! QGIJ-:= JACIJ[N]*QFI[J-1]),
'SKIP',
(NOTR ! QGIJ-:= JACIJ[N]*QFI[J+1]),
QGIJ-:= JACIJ[N]*QFP[J ]
'OUT' ERROR
'ESAC'
'OD' ;
QGIJ*:= D1;
QGIJ := QGIJ/JACIJ[0]+D2*QFI[J]
'OD' ;
PET('FALSE',JAC[I])
'OD' ;
F:= G
'OD'
'END'; #$E#
#E$#
'PR' EJECT 'PR'
#$H#
'PROC' LINGAUS = ('INT' K, L, 'REF''FIELD' F, RHS) 'VOID':
#$B LINGS B$#
'BEGIN'# BY PDZ; INSERTED 850102 #
'REF''NETMAT' JAC = JACOBIAN[L],
'REF'[,]'STATE' QF = Q'OF'F, QR = Q'OF'RHS;
'INT' L1 = 1'LWB'QF, U1 = 1'UPB'QF,
L2 = 2'LWB'QF, U2 = 2'UPB'QF;
( L1/= 'LWB'JAC 'OR' U1/= 'UPB'JAC 'OR'
L1/=1'LWB'QR 'OR' U1/=1'UPB'QR 'OR'
L2/=2'LWB'QR 'OR' U2/=2'UPB'QR ! ERROR);
'INT' ISL, IS, ISU, JSL, JS, JSU;
(LINBACKW ! ISL:=U1; ISU:=L1; IS:=-1 ! ISL:=L1; ISU:=U1; IS:= 1);
(LINREVER ! JSL:=U2; JSU:=L2; JS:=-1 ! JSL:=L2; JSU:=U2; JS:= 1);
'TO' ( LINSYM ! 2*K ! K )
'DO' 'FOR' I 'FROM' ISL 'BY' IS 'TO' ISU
'DO' 'BOOL' LWB = (I=L1), UPB = (I=U1);
'REF'[,]'FLUXAC' JACI = GET(JAC[I]),
'REF'[]'STATE' QRI = QR[ 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' JSL 'BY' JS 'TO' JSU
'DO' 'BOOL' NOTL = (J>L2), NOTR = (J<U2),
'REF'[]'FLUXAC' JACIJ = JACI[J,],
'STATE' RIJ:= QRI[J] ;
'FOR' N 'FROM' L3 'TO' U3
'DO' 'CASE' N+3
'IN' RIJ-:= JACIJ[N]*QFM[J ] ,
(NOTL ! RIJ-:= JACIJ[N]*QFI[J-1]),
'SKIP',
(NOTR ! RIJ-:= JACIJ[N]*QFI[J+1]),
RIJ-:= JACIJ[N]*QFP[J ]
'OUT' ERROR
'ESAC'
'OD' ;
QFI[J]:= RIJ/JACIJ[0]
'OD' ;
PET('FALSE',JAC[I])
'OD';
(LINSYM ! 'INT' IJ:= ISL; ISL:=ISU; ISU:=IJ; IS:=-IS;
IJ:= JSL; JSL:=JSU; JSU:=IJ; JS:=-JS )
'OD'
'END'; #$E#
#E$#
# MDCP ROUTINES # 'PR' EJECT 'PR'
#$H#
'PROC' TAU CORRECT = ('BOOL' SHOW, EXTRA, 'REF''FIELD' QL,RL)'VOID':
#$B TAU.COR B$#
'BEGIN' 'FIELD' TAU; 'BOOL' MORE = SHOW 'OR' EXTRA;
Q'OF'RL:= 'NIL'; RL:= 'NULRHS' QL;
RESIDU2('FALSE',QL,RL);
( SHOW
! PRINT((NEWLINE, "TAU CORRECT",
NEWLINE, "NORM VAN N2(QH)"));
'PRINTNORM' RL );
(MORE! TAU:= 'RESTRICT' RL );
RESIDU1( 'TRUE',QL,RL);
( SHOW ! PRINT((NEWLINE,
"NORM VAN NH1(QH)-N2(QH)"));
'PRINTNORM' RL );
(MORE! RESIDU2('TRUE','RESTRICT'QL,TAU);
( SHOW ! PRINT((NEWLINE,
"NORM VAN TAU[2H,H]"));
'PRINTNORM' TAU);
(EXTRA! TAU*:= (1/3);
TAU := 'PROSCND' TAU;
RL +:= TAU
))
'END'; #$E#
#E$#
#$H#
'PROC' MDCP SMOOTH = ('FIELD' QF )'VOID':
#$B MDCP.SM B$#
'BEGIN' 'REF'[,]'STATE' Q = Q'OF'QF; STATE Z ( 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;
[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,QIO[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 );
QIJ +:= -0.5 * F /FAC
'OD'; QIM:= QIO; HORIN:= HORIZ
'OD'
'END'; #$E#
#E$#
# MULTIGRID (I) # 'PR' EJECT 'PR'
'PROC' LIN RELAX:= ('INT' K,L, 'REF''FIELD' Q, R)'VOID': 'SKIP';
'INT' PCS:=1, SIGMACS:=1, QCS:=1;
#$H#
'PROC' MGCS = ('INT' K,L, 'REF''FIELD' Q, RHS)'VOID':
#$B MGCS B$#
'BEGIN' 'REF''NETMAT' JAC = JACOBIAN[L];
STATE Z(Q);
LIN RELAX(PCS,L, Q, RHS);
(MONICS ! MONITOR("MGCSPRE ",RHS-JAC*Q));
'IF' L>1 '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);
(L<LEVELS ! (STATE2! STATE E(SOL[L]) );
SOL[L+1]:= 'PROLON2' SOL[L] );
( MONIFAS ! MONITOR("FAS ",RHS[L]'MINRESIDU'SOL[L]) )
'OD';
SOL[LEVELS]
'END'; #$E#
#E$#
# MULTIGRID (III) # 'PR' EJECT 'PR'
#$H#
'PROC' STARFAS = ('INT' L, 'REF'[]'FIELD'Q, RHS) 'VOID':
#$B STARFAS B$#
'BEGIN' 'PROC'('BOOL','FIELD','FIELD')'FIELD' CGRESIDU= RESIDU;
RESIDU:= STARRES1;
(MONIREL! MONITOR("RELINA* ",RHS[L]'MINRESIDU'Q[L]));
STARRELAX(PMG,L, Q[L], RHS[L]);
'IF' L>1 '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);
(L<LEVELS ! SOL[L+1]:= 'STARPROL' SOL[L] );
( MONIFAS ! RESIDU:= STARRES1;
MONITOR("FAS ",RHS[L]'MINRESIDU'SOL[L]);
RESIDU:= CGRESIDU )
'OD';
SOL[LEVELS]
'END'; #$E#
#E$#
# STANDARD PROBLEMS. # 'PR' EJECT 'PR'
'REAL' RR,UU,VV,PP,CC,ZZ, UU1,UU2,VV1,VV2,CC1,CC2,ZZ1,ZZ2, DIST;
#$H#
'PROC' PROBLEM1 = 'FIELD':
#$B PROB1 B$#
'BEGIN' PRINT(("STANDAARDPROBLEEM 1"+10*" "+"VSN. 850515", NEWLINE,
"BUMP, TRANSSOON", NEWLINE,
"EQUIDISTANT"));
'FIELD' Q1 = COMPFIELD (1.0,1.0, -2,2, 0,2);
MAPPING:= ('REAL' XX,YY)'POINT': (XX+0.5,YY);
WALLS:= WALLSR; FLUXT:= FLUXTR;
OOST :=('WALL' W,'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE':
WALL (W, 'FALSE', 0, QIN,QAC) ;
WEST :=('WALL'W, 'REAL'XX,YY,'STATE'QIN,'REF''STATAC'QAC)'STATE':
('REAL' X = X'OF'MAPPING(XX,YY);
WALL (W, 'TRUE' , ('ABS'X>0.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]<NRMV[C,0]!NRMH[C,0]!NRMV[C,0])
'OD'; NRM
'END'; #$E#
#E$#
#$H#
'OP' 'PRINTTV' = ('FIELD' F)'VOID':
#$B PRINTTV B$#
'BEGIN' [,]'REAL' N = TOTVAR( F );
PRINT((NEWLINE,"TOT.VAR.: ",SF(N[,1]),
NEWLINE,"MAX.VAR.: ",SF(N[,0]),NEWLINE ))
'END' ; #$E#
#E$#
# OUTPUT ROUTINES # 'PR' EJECT 'PR'
#$H#
'OP' 'REPORT' = ('FIELD' F)'FIELD':
#$B REPORT B$#
'BEGIN' 'REAL' DDX = DDX'OF'F, DDY = DDY'OF'F;
'REF'[,]'STATE' Q = Q'OF'F;
'INT' L2= 1'LWB'Q, U2= 1'UPB'Q, L3= 2'LWB'Q, U3= 2'UPB'Q;
'STRING' LS2= WHOLE(L2,0), US2= WHOLE(U2,0),
LS3= WHOLE(L3,0), US3= WHOLE(U3,0);
'POINT' NW= MAPPING(XXN,YYW), NO= MAPPING(XXN,YYO),
ZW= MAPPING(XXZ,YYW), ZO= MAPPING(XXZ,YYO);
PRINT((NEWLINE,"VELD: [ " + LS2 + " : " + US2 +" ] * [ "+
LS3 + " : " + US3 +" ] ",
NEWLINE,"X IN ("+ ST((L2-1)*DDX)+ST(U2*DDX) + " ) "+
" IN ("+ ST( 'POINT'((XXN,XXZ)) )+ " ) "+
" NW ("+ ST( NW )+ " ) "+
" NO ("+ ST( NO )+ " ) ",
NEWLINE,"Y IN ("+ ST((L3-1)*DDY)+ST(U3*DDY) + " ) "+
" IN ("+ ST( 'POINT'((YYW,YYO)) )+ " ) "+
" ZW ("+ ST( ZW )+ " ) "+
" ZO ("+ ST( ZO )+ " ) "
)); F
'END';#$E#
#E$#
#$H#
'OP' 'PRINT' = ('FIELD' F)'VOID':
#$B PRNT.1 B$#
'FOR' K 'TO' 4
'DO' K'PRINT'F
'OD' ; #$E#
#E$#
#$H#
'OP' 'PRINT' = ('INT' K,'FIELD' F)'VOID':
#$B PRNT.2 B$#
'BEGIN' 'REF'[,]'STATE' Q = (Q'OF'F)[@1,@1];
'STRING' S:= ( TYPE'OF'F <= 1
! (K!"MASS","X-MOMT","Y-MOMT","ENERGY")
! (K!"X-SPD","Y-SPD","SND-SPD","LN(ENTROPY)") );
( TYPE'OF'F = 0 ! S+:= "-FLUX");
( TYPE'OF'F = 1 ! S+:= "-DENSITY");
PRINT((NEWLINE,S));
'FOR' I 'TO' 1'UPB'Q
'DO' PRINT((NEWLINE, ST( []'REAL' ((
(K!C1'OF'Q[I,],C2'OF'Q[I,],C3'OF'Q[I,],C4'OF'Q[I,])
)) )
))
'OD'
'END'; #$E#
#E$#
#$H#
'OP' 'PRINT' = ('PROC'('STATE')'REAL' P,'FIELD' F)'VOID':
#$B OPPRINT B$#
'BEGIN' 'REF'[,]'STATE' Q = (Q'OF'F)[@1,@1];
STATE E (F);
'FOR' I 'TO' 1'UPB'Q
'DO' 'STRING' S:= "";
'FOR' J 'TO' 2'UPB'Q 'DO' S+:=ST(P(Q[I,J])) 'OD';
PRINT((NEWLINE,S))
'OD'
'END'; #$E#
#E$#
'OP' 'FIX' = ('REAL'X)'STRING' : FIXED(X,7,3 ) ;
'OP' 'FLT' = ('REAL'X)'STRING' : FLOAT(X,7,3,1) ;
'MODE' 'UNI' = 'UNION'('REAL','POINT','WALL',[]'REAL','STATE' ) ;
#$H#
'PROC' ST = ('UNI' U)'STRING':
#$B ST.UNI B$#
'CASE' U 'IN'
('REAL' X): " "+'FIX' X ,
('POINT' P): " "+'FIX'X'OF'P+" "+'FIX'Y'OF'P ,
('WALL' W): " " + 'FIX' C'OF'W + " " + 'FIX' S'OF'W +
" " + 'FIX'DS'OF'W ,
([]'REAL' A): ('STRING'S:="";'FOR'I'FROM''LWB'A'TO''UPB'A
'DO' S+:=" "+'FIX' A[I] 'OD' ; S ) ,
('STATE' S): " " + 'FIX'C1'OF'S + " " + 'FIX'C2'OF'S +
" " + 'FIX'C3'OF'S + " " + 'FIX'C4'OF'S
'ESAC' ; #$E#
#E$#
#$H#
'PROC' SF = ('UNI' U)'STRING':
#$B SF.UNI B$#
'CASE' U 'IN'
('REAL' X): " "+'FLT' X ,
('POINT' P): " "+'FLT'X'OF'P+" "+'FLT'Y'OF'P ,
('WALL' W): " " + 'FLT' C'OF'W + " " + 'FLT' S'OF'W +
" " + 'FLT'DS'OF'W ,
([]'REAL' A): ('STRING'S:="";'FOR'I'FROM''LWB'A'TO''UPB'A
'DO' S+:=" "+'FLT' A[I] 'OD' ; S ) ,
('STATE' S): " " + 'FLT'C1'OF'S + " " + 'FLT'C2'OF'S +
" " + 'FLT'C3'OF'S + " " + 'FLT'C4'OF'S
'ESAC' ; #$E#
#E$#
#$H#
'PROC' FLUXNORM = ('FIELD'F,'REF''FLUX' NORM,RATIO)'REAL':
#$B FLUXXN B$#
'BEGIN' 'REF'[,]'STATE' Q = (Q'OF'F)[@1,@1];
'INT' N1=1'UPB'Q, N2=2'UPB'Q;
'FOR' K 'TO' 4
'DO' 'REF'[,]'REAL' QK=(K!C1'OF'Q,C2'OF'Q,C3'OF'Q,C4'OF'Q);
'REAL' T,S:= 0, S1:= 0;
'FOR' I 'TO' N1 'DO' 'FOR'J 'TO' N2 'DO'
T:= 'ABS' QK[I,J]; S1+:=T; (T>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]<MIN! MIN:=X[I]) 'OD';
FACTOR:= W/(MAX-MIN);
PRINT((NEWLINE,NEWLINE,"RANGE = ",MIN," TOT ",MAX));
'FOR'I'FROM'L'TO'U
'DO' K:= 'ROUND'( (X[I]-MIN)*FACTOR );
PRINT((NEWLINE,K*" ","*"))
'OD'
'END'; #$E#
#E$#
'PROC' MONITOR := ('STRING' T,'FIELD' F)'VOID':
'BEGIN' PRINT((NEWLINE,">>> 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<LEVELS
!PRINT((NEWPAGE));
Q[L+1]:= 'PROLON2' STATE E(Q[L])
)
'OD'
'END'