Forgot your password?
Document Actions
P. W. Hemker and P. M. de Zeeuw. A library of multigrid routines.
P. W. Hemker and P. M. de Zeeuw. A library of multigrid routines. Centrum voor Wiskunde en Informatica, Numerieka Wiskunde, 1982.
Size 103.5 kB - File type text/plain
Click here to get the file
File contents
# *****************************************************************#
# #
# #
# A LIBRARY OF MULTIGRID ROUTINES #
# #
# #
# #
# #
# THIS ALGOL 68 PRELUDE IS A COLLECTION OF PROGRAM MODULES, WHICH #
# FORM A LIBRARY FOR NUMERICAL EXPERIMENTS WITH MULTIGRID METHODS.#
# THE PRELUDE WAS DEVELOPED FROM 800918 ON AND IS CONTINUOUSLY #
# CHANGING. HOWEVER, AT THIS STAGE (820224), IT HAS TEMPORARILY #
# REACHED A STATIONARY STATE AND IS USED FOR TESTING DIFFERENT #
# MULTIGRID STRATEGIES. #
# #
# SINCE THE PRESENT VERSION WORKS ONLY FOR RECTANGULAR REGIONS #
# AND WITHOUT BACK-STORAGE FACILITIES, IT IS STILL SURVEYABLE AND #
# WRITTEN IN A RATHER READABLE KIND OF ALGOL 68. THEREFORE IT IS #
# A CONCISE AND PRECISE DESCRIPTION OF MANY ESSENTIAL PARTS OF #
# THE ALGORITHMS. #
# OTHER VERSIONS OF THE ROUTINES, FOR MORE GENERAL DOMAINS AND #
# WITH THE USE OF DISC STORAGE FACILITIES, ARE ALSO AVAILABLE #
# BUT NOT INCLUDED IN THIS LIBRARY. #
# WE KNOW THAT, AT A NUMBER OF SPOTS, CERTAINLY NOT THE OPTIMAL #
# ALGORITHM HAS BEEN USED AND THE CODE CAN BE IMPROVED. #
# FURTHER WE BELIEVE THAT ANY PIECE OF SOFTWARE OF SOME LENGTH #
# CONTAINS AT LEAST ONE BUG. #
# WITH THIS RESERVE WE GIVE HERE THE PRELIMINARY VERSION OF OUR #
# MULTIGRID LIBRARY. #
# #
# #
# P.W.HEMKER #
# P.M.DE ZEEUW #
# #
# #
# #
# #
# *****************************************************************#
'PR' EJECT 'PR'
# *****************************************************************#
# #
# #
# CONTENTS: #
# #
# MODE DECLARATIONS 3 #
# ZERO AND COPY OPERATI0NS 4 #
# ELEMENTARY NET OPERATI0NS 5 #
# ELEMENTARY GRID OPERATI0NS 6 #
# NORM OPERATORS 7 #
# PROLONGATIONS 8 #
# RESTRICTIONS 11 #
# FINITE ELEMENT METHOD (FEM) 14 #
# FINITE ELEMENT METHOD (S-U P-G FEM) 17 #
# GALERKIN OPERATOR CONSTRUCTION 23 #
# SOLVE BY GAUSSIAN ELIMINATION 25 #
# OPERATOR EVALUATION 27 #
# RELAXATION METHODS 29 #
# POINT RELAXATION PROCEDURES 30 #
# LINE RELAXATION PROCEDURES 31 #
# DAMPED JACOBI RELAXATION 32 #
# ILU RELAXATION 33 #
# ILLU RELAXATION 34 #
# OTHER POINTWISE RELAXATION 38 #
# GRID STYLE PROCEDURES 40 #
# CENTRAL PROCEDURES GRID STYLE 41 #
# CENTRAL PROCEDURES NEW STYLE 43 #
# NAG CONTRIBUTION : THE MORE EFFICIENT VERSION 45 #
# PRELIMINARY NAG MGM ROUTINE 50 #
# BLACK-BOX SOLUTION PROCEDURE 52 #
# ASYMMETRIC TRANSFERS 53 #
# COLLECTION OF PRIMARY TESTPROBLEMS 56 #
# OUTPUT PROCEDURES 57 #
# GLOBAL PARAMETERS 58 #
# #
# *****************************************************************#
'PR' EJECT 'PR'
MGLIB: # A MODULAR 2-D MULTIGRID ALGOL 68 PROGRAM #
# FOR EXPERIMENTAL PURPOSES. PWH/800918 #
# VSN 830421 PF: MGTEXT #
'BEGIN'
# MODE AND PRIORITY DECLARATIONS #
'MODE' 'VEC' = 'REF'[ ]'REAL';
'MODE' 'NET' = 'REF'[, ]'REAL';
'MODE' 'NETMAT' = 'REF'[,,]'REAL';
'MODE' 'GRID' = 'STRUCT'('INT' LEVEL, 'NET' NET );
'REAL' # GLOBAL FOR MODE GRID #
X0:=0.0, Y0:=0.0, XH0:=1.0, YH0:=1.0;
'MODE' 'PROBLEM' = 'STRUCT'('BOOL'LINEAR,
'REF'[,]'REAL' OMEGA,
'PROC'('REF''REAL','REF''REAL','REF''REAL','REF''REAL',
'REAL','REAL' )'VOID' PRINCIPLE PART,
'PROC'('REF''REAL','REF''REAL','REF''REAL','REF''REAL',
'REAL','REAL','REAL')'VOID' LOWER ORDER,
'PROC'('REF''REAL',
'REAL','REAL' )'VOID' BOUNDARY PP,
'PROC'('REF''REAL','REF''REAL',
'REAL','REAL','REAL')'VOID' BOUNDARY LO);
'PRIO' 'NORM' = 8, 'RESIDUALNORM' = 8;
# ZERO AND COPY OPERATI0NS # 'PR' EJECT 'PR'
'OP' 'ZERO' = ('VEC' A ) 'VEC':
('FOR' I 'FROM' 'LWB'A 'TO' 'UPB'A
'DO' A[I]:= 0.0 'OD'; A );
'OP' 'ZERO' = ('NET' A ) 'NET':
('FOR' I 'FROM' 1'LWB'A 'TO' 1'UPB'A
'DO' 'ZERO' A[I,] 'OD'; A );
'OP' 'ZERO' = ('NETMAT' A ) 'NETMAT':
('FOR' I 'FROM' 1'LWB'A 'TO' 1'UPB'A
'DO' 'ZERO' A[I,,] 'OD'; A );
'OP' 'ZERO' = ('GRID' A )'GRID': ('ZERO' (NET'OF'A); A);
'OP' 'COPY' = ('VEC' N ) 'VEC':
('HEAP' ['LWB'N:'UPB'N]'REAL' :=N);
'OP' 'COPY' = ('NET' N ) 'NET':
('HEAP' [1'LWB'N:1'UPB'N,2'LWB'N:2'UPB'N]'REAL' :=N);
'OP' 'COPY' = ('NETMAT' N ) 'NETMAT':
('HEAP' [1'LWB'N:1'UPB'N,2'LWB'N:2'UPB'N,
3'LWB'N:3'UPB'N]'REAL' :=N);
'OP' 'COPY' = ('GRID' A )'GRID':
('NET'N= NET'OF'A; (LEVEL'OF'A,
'HEAP' [1'LWB'N:1'UPB'N,2'LWB'N:2'UPB'N]'REAL' :=N));
# ELEMENTARY NET OPERATI0NS # 'PR' EJECT 'PR'
'OP' + = ('NET' AA,BB )'NET':
#$B NET+ B$#
( 'INT' L1 = 1'LWB'AA, L2 = 2'LWB'AA,
U1 = 1'UPB'AA, U2 = 2'UPB'AA;
'HEAP'[L1:U1,L2:U2]'REAL' CC;
'FOR' I 'FROM' L1 'TO' U1 'DO'
'FOR' J 'FROM' L2 'TO' U2 'DO'
CC[I,J]:= AA[I,J] + BB[I,J] 'OD''OD';
CC )
#$E#;
#E$#
'OP' - = ('NET' AA,BB )'NET':
#$B NET- B$#
( 'INT' L1 = 1'LWB'AA, L2 = 2'LWB'AA,
U1 = 1'UPB'AA, U2 = 2'UPB'AA;
'HEAP'[L1:U1,L2:U2]'REAL' CC;
'FOR' I 'FROM' L1 'TO' U1 'DO'
'FOR' J 'FROM' L2 'TO' U2 'DO'
CC[I,J]:= AA[I,J] - BB[I,J] 'OD''OD';
CC )
#$E#;
#E$#
'OP' +:= = ('NET' AA,BB )'NET':
#$B NET= B$#
( 'INT' L1 = 1'LWB'AA, L2 = 2'LWB'AA,
U1 = 1'UPB'AA, U2 = 2'UPB'AA;
'FOR' I 'FROM' L1 'TO' U1 'DO'
'FOR' J 'FROM' L2 'TO' U2 'DO'
AA[I,J]+:= BB[I,J] 'OD''OD';
AA )
#$E#;
#E$#
'OP' * = ('REAL' S,'NET' N ) 'NET':
#$B NET* B$#
( 'INT' L1 = 1'LWB' N, L2 = 2'LWB'N,
U1 = 1'UPB' N, U2 = 2'UPB'N;
'HEAP'[L1:U1,L2:U2]'REAL' NN;
'FOR' I 'FROM' L1 'TO' U1 'DO'
'FOR' J 'FROM' L2 'TO' U2 'DO'
NN[I,J]:= S * N[I,J] 'OD' 'OD';
NN )
#$E#;
#E$#
'OP' / = ('NET' G,'REAL' S) 'NET':
#$B NET/ B$#
( 'REAL' T = 1.0/S; T * G )
#$E#;
#E$#
# ELEMENTARY GRID OPERATI0NS # 'PR' EJECT 'PR'
'OP' + = ('GRID' A,B )'GRID':
#$B GRID+ B$#
( (LEVEL'OF'A /= LEVEL'OF'B ! ERROR );
(LEVEL'OF'A, NET'OF'A + NET'OF'B ) )
#$E#;
#E$#
'OP' - = ('GRID' A,B )'GRID':
#$B GRID- B$#
( (LEVEL'OF'A /= LEVEL'OF'B ! ERROR );
(LEVEL'OF'A, NET'OF'A - NET'OF'B ) )
#$E#;
#E$#
'OP' - = ('GRID' A,'PROC'('REAL','REAL')'REAL'SOL )'GRID':
#$B GRID-- B$#
( 'INT' L:= LEVEL'OF'A; L:= 2**L;
'NET' AA= NET'OF'A;
'INT' L1 = 1'LWB'AA, L2 = 2'LWB'AA,
U1 = 1'UPB'AA, U2 = 2'UPB'AA;
'REAL' H1 = XH0/L, H2 = YH0/L;
'HEAP'[L1:U1,L2:U2]'REAL' CC;
'REAL' X:= X0+L1*H1, Y;
'FOR' I 'FROM' L1 'TO' U1 'DO' Y:= Y0+L2*H2;
'FOR' J 'FROM' L2 'TO' U2 'DO'
CC[I,J]:= AA[I,J] - SOL(X,Y); Y+:= H1 'OD';
X+:= H2 'OD';
(LEVEL'OF'A, CC ) )
#$E#;
#E$#
'OP' +:= = ('REF''GRID' A,'GRID' B )'REF' 'GRID':
#$B GRID= B$#
( (LEVEL'OF'A /= LEVEL'OF'B ! ERROR );
NET'OF'A +:= NET'OF'B; A)
#$E#;
#E$#
'OP' * = ('REAL' S,'GRID' G ) 'GRID':
#$B GRID* B$#
(LEVEL'OF'G, S*NET'OF'G)
#$E#;
#E$#
'OP' / = ('GRID' G,'REAL' S) 'GRID':
#$B GRID/ B$#
( 'REAL' T = 1.0/S; T * G )
#$E#;
#E$#
# NORM OPERATORS # 'PR' EJECT 'PR'
'OP' 'NORM' = ('BOOL' TRUE, 'NET' N) []'REAL':
#$B NRMNT B$#
('REAL' TEST = ( TRUE ! MAXREAL ! 1.0E+15);
'REAL'S:= 0,S1:=0,S2:=0,T;
'INT' L1:= 1'LWB'N, B1:= 1'UPB'N,
L2:= 2'LWB'N, B2:= 2'UPB'N,NN;
NN:= (B1-L1+1)*(B2-L2+1);
'FOR' I 'FROM' L1 'TO' B1 'DO'
'FOR' J 'FROM' L2 'TO' B2 'DO'
( T:= 'ABS'N[I,J]; T < TEST
! S1+:=T; S2+:=T*T; (T>S ! S:=T )
! NN-:= 1 )
'OD' 'OD'; (S1/NN,SQRT(S2/NN),S) )
#$E#;
#E$#
'OP' 'NORM' = ('NET' N) []'REAL':
#$B NRM B$#
'TRUE''NORM'N
#$E#;
#E$#
'PROC' RESIDUALNORM = ('NETMAT' M,'NET' U,F)[]'REAL':
# U = ITERAND, F = RIGHT HAND SIDE #
#$B NRMR B$#
'BEGIN' 'INT' L1= 1'LWB'U, L2= 2'LWB'U, L3=-3,
U1= 1'UPB'U, U2= 2'UPB'U;
'INT' NN:= (U1-L1+1)*(U2-L2+1);
'REAL' TEST = 1.0E+15;
'REAL' V, S1:= 0, S2:= 0, S3:=0, T;
'REF'[ ]'REAL' UIM:= U[L1,@L2], UI, UIP:= U[L1,@L2];
'FOR' I 'FROM' L1 'TO' U1
'DO' ( UI:= UIP; I = U1 ! 'SKIP' ! UIP:= U[I+1,@L2] );
'REF' []'REAL' FI = F[I,@L2];
'REF'[,]'REAL' MI = M[I,@L2,@L3];
'INT' JM:= L2, JJ, JP:= L2;
'FOR' J 'FROM' L2 'TO' U2
'DO' ( JJ:= JP; J=U2 ! 'SKIP' ! JP+:= 1 );
'REF'[]'REAL' MIJ = MI[JJ,@L3];
( MIJ[0] > TEST ! NN-:= 1 !
T := MIJ[-3]*UIM[JJ] + MIJ[-2]*UIM[JP] +
MIJ[-1]*UI [JM] + MIJ[ 0]*UI [JJ] + MIJ[ 1]*UI [JP] +
MIJ[ 2]*UIP[JM] + MIJ[ 3]*UIP[JJ] ;
T := 'ABS' ( FI[JJ] - T );
S1+:= T; S2+:= T*T; ( T>S3 ! S3:= T ));
JM:= JJ
'OD';
UIM:= UI
'OD'; (S1/NN,SQRT(S2/NN),S3)
'END'
#$E#;
#E$#
# PROLONGATIONS # 'PR' EJECT 'PR'
'PROC' LIN INT POL = ('NET' NET )'NET':
#$B PROL7 B$#
'BEGIN' 'INT' L1 = 1'LWB'NET, L2 = 2'LWB'NET,
B1 = 1'UPB'NET, B2 = 2'UPB'NET;
'HEAP'[2*L1:2*B1,2*L2:2*B2]'REAL' FINE;
'INT' JJ; 'REAL' U2,U3,U4;
'REF'[]'REAL' UIP= NET[L1,@L2],
UPP= FINE[2*L1,@2*L2];
JJ:= 2*L2; UPP[JJ]:= U4:= UIP[L2];
'FOR' JP 'FROM' L2+1 'TO' B2
'DO' U3:= U4; U4:= UIP[JP];
UPP[JJ+:=1]:= (U3+U4)/2;
UPP[JJ+:=1]:= U4
'OD';
'FOR' IP 'FROM' L1+1 'TO' B1
'DO' 'REF'[]'REAL' UI = NET [IP-1 ,@ L2],
UIP = NET [IP ,@ L2],
UMM = FINE[2*IP-1,@2*L2],
UPP = FINE[2*IP ,@2*L2];
JJ:= 2*L2; U2:= UI[L2]; U4:= UIP[L2];
UMM[JJ]:= (U2+U4)/2; UPP[JJ]:= U4;
'FOR' JP 'FROM' L2+1 'TO' B2
'DO' JJ+:= 1; U2:= UI [JP];
U3:= U4; U4:= UIP[JP];
UMM[JJ] := (U2+U3)/2;
UPP[JJ] := (U3+U4)/2;
JJ+:= 1;
UMM[JJ] := (U2+U4)/2;
UPP[JJ] := U4
'OD'
'OD'; FINE
'END'
#$E#;
#E$#
'PROC' BIL INT POL = ('NET' NET) 'NET':
#$B PROL9 B$#
'BEGIN' 'INT' L1 = 1'LWB'NET, L2 = 2'LWB'NET,
B1 = 1'UPB'NET, B2 = 2'UPB'NET;
'HEAP'[2*L1:2*B1,2*L2:2*B2]'REAL' FINE;
'INT' JJ; 'REAL' U1,U2,U3,U4;
'REF'[]'REAL' UIP= NET[L1,@L2],
UPP= FINE[2*L1,@2*L2];
JJ:= 2*L2; UPP[JJ]:= U4:= UIP[L2];
'FOR' JP 'FROM' L2+1 'TO' B2
'DO' U3:= U4; U4:= UIP[JP];
UPP[JJ+:=1]:= (U3+U4)/2;
UPP[JJ+:=1]:= U4
'OD';
'FOR' IP 'FROM' L1+1 'TO' B1
'DO' 'REF'[]'REAL' UI = NET [IP-1 ,@ L2],
UIP = NET [IP ,@ L2],
UMM = FINE[2*IP-1,@2*L2],
UPP = FINE[2*IP ,@2*L2];
JJ:= 2*L2; U2:= UI[L2]; U4:= UIP[L2];
UMM[JJ]:= (U2+U4)/2; UPP[JJ]:= U4;
'FOR' JP 'FROM' L2+1 'TO' B2
'DO' JJ+:= 1;
U1:= U2; U2:= UI [JP];
U3:= U4; U4:= UIP[JP];
UMM[JJ] := (U1+U2+U3+U4)/4;
UPP[JJ] := (U3+U4)/2;
JJ+:= 1;
UMM[JJ] := (U2+U4)/2;
UPP[JJ] := U4
'OD'
'OD'; FINE
'END'
#$E#;
#E$#
# ORIENTATION:
L2 B2
--------------------------> Y
! L1 JM JJ JP
!
! UIM O -- O -- O
! ! / ! / !
! UII O -- U1 -- U2
! ! / ! / !
X ! UIP O -- U3 -- U4
V B1
----------------------> Y
! L2 B2
! L1 STAR ORIENTATION:
! U1 ------ U2 -----------------> Y
! ! / ! ! J
! ! / ! ! -3 -2
! U3 ------ U4 ! -1 0 1
! B1 ! 2 3
V X ! I
V
#
'PR' EJECT 'PR'
'PROC' SQR INT POL = ('NET' NET )'NET':
#$B PROLS B$#
'IF' 'INT' L1 = 1'LWB'NET, L2 = 2'LWB'NET,
B1 = 1'UPB'NET, B2 = 2'UPB'NET;
'ODD'(B1-L1) 'OR' 'ODD'(B2-L2)
'THEN' LIN INT POL (NET)
'ELSE' 'INT' LL1 = 2*L1, LL2 = 2*L2;
'HEAP'[LL1:2*B1,LL2:2*B2]'REAL' FINE;
'INT' JJ, JP;
'REAL' X1, X2, X3, Y1, Y2, Y3, Z1, Z2, Z3, YY2, YY3, ZZ2, ZZ3;
'REF'[]'REAL' UI= NET[ L1,@L2], FI= FINE[LL1,];
FI[LL2]:= X1:= UI[L2]; JJ:= LL2+1;
'FOR' J 'FROM' L2+1 'BY' 2 'TO' B2-1
'DO' X2:= UI[J]; X3:= UI[J+1];
FI[JJ:JJ+3] :=( ( 3*(X1 + 2*X2) - X3 )/8, X2,
( -X1 + 3*(2*X2 + X3 ))/8, X3 );
JJ +:= 4; X1:= X3
'OD';
'FOR' II 'FROM' L1+1 'BY' 2 'TO' B1-1
'DO' 'REF'[ ]'REAL' UIM= NET[II-1,@L2], UII= NET[II ,@L2],
UIP= NET[II+1,@L2];
'REF'[,]'REAL' FINEI = FINE[2*II-1:2*II+2,@LL2];
X3:= UIM[L2] /8;
Y3:= ( YY3:= UII[L2] )/4;
Z3:= ( ZZ3:= UIP[L2] )/8;
FINEI[,LL2]:= ( 3*(X3+Y3) - Z3, YY3, 3*(Y3+Z3) - X3, ZZ3 );
'FOR' JJ 'FROM' L2+1 'BY' 2 'TO' B2-1
'DO' JP:= JJ+1; X1:= X3; Y1:= Y3; Z1:= Z3;
X2:= UIM[JJ] /4; X3:= UIM[JP] /8;
Y2:= ( YY2:= UII[JJ] )/4; Y3:= ( YY3:= UII[JP] )/4;
Z2:= ( ZZ2:= UIP[JJ] )/4; Z3:= ( ZZ3:= UIP[JP] )/8;
FINEI[,2*JJ-1:2*JJ+2]:=
((2*(X2+Y1)-Z1+Y2-X3,
2*(X2+Y2)-X1+Y1-Z1,
3*(X3+Y2)-Z1,
3*(X3+Y3)-Z3 ),
(2*(Y1+Y2)-X1+X2-X3, YY2,
2*(Y2+Y3)-Z1+Z2-Z3, YY3 ),
(3*(Z1+Y2)-X3,
2*(Y2+Z2)-X3+Y3-Z3,
2*(Z2+Y3)-X3+Y2-Z1,
3*(Z3+Y3)-X3 ),
(3*(Z1+Z2)-Z3, ZZ2, 3*(Z3+Z2)-Z1, ZZ3 ))
'OD' 'OD';
FINE
'FI'
#$E#;
#E$#
# RESTRICTIONS # 'PR' EJECT 'PR'
'PROC' INJECTION = ('NET' NET)'NET':
#$B RESTR0 B$#
'BEGIN' 'INT' NL1 = (1'LWB'NET)'OVER'2, NL2 = (2'LWB'NET)'OVER'2,
NU1 = (1'UPB'NET)'OVER'2, NU2 = (2'UPB'NET)'OVER'2;
'HEAP'[NL1:NU1,NL2:NU2]'REAL' COARSE;
'FOR' I 'FROM' NL1 'TO' NU1 'DO'
'FOR' J 'FROM' NL2 'TO' NU2 'DO'
COARSE[I,J]:= NET[2*I,2*J] 'OD' 'OD';
COARSE
'END'
#$E#;
#E$#
'PROC' INJ WEIGHT = ('NET' NET)'NET':
#$B RESTR1 B$#
'BEGIN' 'INT' NL1 = (1'LWB'NET)'OVER'2, NL2 = (2'LWB'NET)'OVER'2,
NU1 = (1'UPB'NET)'OVER'2, NU2 = (2'UPB'NET)'OVER'2;
'HEAP'[NL1:NU1,NL2:NU2]'REAL' COARSE;
COARSE[NL1,NL2]:= 2.0*NET[2*NL1,2*NL2];
COARSE[NU1,NL2]:= 2.5*NET[2*NU1,2*NL2];
COARSE[NL1,NU2]:= 2.5*NET[2*NL1,2*NU2];
COARSE[NU1,NU2]:= 2.0*NET[2*NU1,2*NU2];
'FOR' J 'FROM' NL2+1 'TO' NU2-1
'DO' COARSE[NL1,J]:= 3.0*NET[2*NL1,2*J];
COARSE[NU1,J]:= 3.0*NET[2*NU1,2*J]
'OD';
'FOR' I 'FROM' NL1+1 'TO' NU1-1
'DO' COARSE[I,NL2]:= 3*NET[2*I,2*NL2];
COARSE[I,NU2]:= 3*NET[2*I,2*NU2];
'FOR' J 'FROM' NL2+1 'TO' NU2-1
'DO' COARSE[I,J] := 4*NET[2*I,2*J] 'OD'
'OD'; COARSE
'END'
#$E#;
#E$#
'PR' EJECT 'PR'
'PROC' HALF WEIGHT = ('NET' NET)'NET':
#$B RESTR5 B$#
'BEGIN' 'INT' L1 = (1'LWB'NET), L2 = (2'LWB'NET),
U1 = (1'UPB'NET), U2 = (2'UPB'NET);
'INT' NL1 = L1'OVER'2, NL2 = L2'OVER'2,
NU1 = U1'OVER'2, NU2 = U2'OVER'2;
'HEAP'[NL1:NU1,NL2:NU2]'REAL' COARSE;
'INT' II,IM,IP, JJ,JM,JP;
[L2:U2]'REAL' ZERO; 'ZERO' ZERO;
'FOR' I 'FROM' NL1 'TO' NU1
'DO' II:= 2*I; IP:= II+1; IM:= II-1;
'REF'[]'REAL' CI = COARSE[I,@NL2],
UIM= ( IM<L1 ! ZERO ! NET[IM,@L2]),
UII= NET[II,@L2] ,
UIP= ( IP>U1 ! ZERO ! NET[IP,@L2]);
JP:= L2+1;
CI[NL2]:= 0.5 *( UIM[L2]+UIM[JP]
+2*UII[L2]+UII[JP]
+UIP[L2]);
'FOR' J 'FROM' NL2+1 'TO' NU2-1
'DO' JJ:= 2*J; JP:= JJ+1; JM:= JJ-1;
CI[J]:= 0.5 *( UIM[JJ]
+UII[JM]+4*UII[JJ]+UII[JP]
+UIP[JJ] )
'OD';
JM:=U2-1;
CI[NU2]:= 0.5 *( UIM[U2]
+UII[JM]+2*UII[U2]
+UIP[JM]+ UIP[U2])
'OD' ;
COARSE
'END'
#$E#;
#E$#
'PROC' LIN WEIGHT = ('NET' NET)'NET':
#$B RESTR7 B$#
'BEGIN' 'INT' L1 = (1'LWB'NET), L2 = (2'LWB'NET),
U1 = (1'UPB'NET), U2 = (2'UPB'NET);
'INT' NL1 = L1'OVER'2, NL2 = L2'OVER'2,
NU1 = U1'OVER'2, NU2 = U2'OVER'2;
'HEAP'[NL1:NU1,NL2:NU2]'REAL' COARSE;
'INT' II,IM,IP, JJ,JM,JP;
[L2:U2]'REAL' ZERO; 'ZERO' ZERO;
'FOR' I 'FROM' NL1 'TO' NU1
'DO' II:= 2*I; IP:= II+1; IM:= II-1;
'REF'[]'REAL' CI = COARSE[I,@NL2],
UIM= ( IM<L1 ! ZERO ! NET[IM,@L2]),
UII= NET[II,@L2] ,
UIP= ( IP>U1 ! ZERO ! NET[IP,@L2]);
JP:= L2+1;
CI[NL2]:= 0.5 *( UIM[L2]+UIM[JP]
+2*UII[L2]+UII[JP]
+UIP[L2]);
'FOR' J 'FROM' NL2+1 'TO' NU2-1
'DO' JJ:= 2*J; JP:= JJ+1; JM:= JJ-1;
CI[J]:= 0.5 *( UIM[JJ]+UIM[JP]
+UII[JM]+2*UII[JJ]+UII[JP]
+UIP[JM]+ UIP[JJ] )
'OD';
JM:=U2-1;
CI[NU2]:= 0.5 *( UIM[U2]
+UII[JM]+2*UII[U2]
+UIP[JM]+ UIP[U2])
'OD'; COARSE
'END'
#$E#;
#E$#
'PROC' FULL WEIGHT = ('NET' NET)'NET':
#$B RESTR9 B$#
'BEGIN' 'INT' L1 = (1'LWB'NET), L2 = (2'LWB'NET),
U1 = (1'UPB'NET), U2 = (2'UPB'NET);
'INT' NL1 = L1'OVER'2, NL2 = L2'OVER'2,
NU1 = U1'OVER'2, NU2 = U2'OVER'2;
'HEAP'[NL1:NU1,NL2:NU2]'REAL' COARSE;
'INT' II,IM,IP, JJ,JM,JP;
[L2:U2]'REAL' ZERO; 'ZERO' ZERO;
'FOR' I 'FROM' NL1 'TO' NU1
'DO' II:= 2*I; IP:= II+1; IM:= II-1;
'REF'[]'REAL' CI = COARSE[I,@NL2],
UIM= ( IM<L1 ! ZERO ! NET[IM,@L2]),
UII= NET[II,@L2] ,
UIP= ( IP>U1 ! ZERO ! NET[IP,@L2]);
JP:= L2+1;
CI[NL2]:= 0.5 *( UIM[L2]+UIM[JP]
+2*UII[L2]+UII[JP]
+UIP[L2]);
'FOR' J 'FROM' NL2+1 'TO' NU2-1
'DO' JJ:= 2*J; JP:= JJ+1; JM:= JJ-1;
CI[J]:= 0.25*( UIM[JM]+2*UIM[JJ]+ UIM[JP]
+2*UII[JM]+4*UII[JJ]+2*UII[JP]
+ UIP[JM]+2*UIP[JJ]+ UIP[JP])
'OD';
JM:= U2-1;
CI[NU2]:= 0.5 *( UIM[U2]
+UII[JM]+2*UII[U2]
+UIP[JM]+ UIP[U2])
'OD' ;
COARSE
'END'
#$E#;
#E$#
# FINITE ELEMENT METHOD (FEM) # 'PR' EJECT 'PR'
# FEM DISCRETIZATION FOR MG /PWH800925 / VSN.810311 #
'PROC' FEM = ('PROBLEM' PROBLEM, 'GRID' U,
'REF''NET' RHS, 'REF''NETMAT' JAC )'VOID':
#$B FEM B$#
'BEGIN' 'REF'[,]'REAL' NET = NET'OF' U;
'INT' L1 = 1'LWB' NET, L2 = 2'LWB'NET,
B1 = 1'UPB' NET, B2 = 2'UPB'NET;
'HEAP' [L1:B1,L2:B2,-3:3 ]'REAL' STIFF; 'ZERO' STIFF;
'HEAP' [L1:B1,L2:B2]'REAL' FORCE; 'ZERO' FORCE;
'REAL' H = XH0/(2**(LEVEL'OF'U)),
K = YH0/(2**(LEVEL'OF'U));
'PROC'('REF''REAL','REF''REAL','REF''REAL','REF''REAL',
'REAL','REAL' )'VOID' PRINCIPLE PART
= PRINCIPLE PART 'OF' PROBLEM;
'PROC'('REF''REAL','REF''REAL','REF''REAL','REF''REAL',
'REAL','REAL','REAL')'VOID' LOWER ORDER
= LOWER ORDER 'OF' PROBLEM;
[]'INT' ORDER = (0,1,3,-1,0,2,-3,-2,0,
0,2,3,-2,0,1,-3,-1,0);
'REAL' H3 = H/3, K3 = K/3;
[]'REAL' PXT1 = (-1/H, 0 ,+1/H, 0), PXT2 = ( 0,-1/H, 0 ,+1/H),
PYT1 = (-1/K,+1/K, 0 , 0), PYT2 = ( 0, 0 ,-1/K,+1/K);
'INT' LB,UB;
'REAL' XN,XS,YW,YE, AXX,AXY,AYX,AYY,
BBX,BBY, KK,PE,FF, WEIGHT:= H*K/6;
[1:4]'REAL' PX,PY,BX,BY,C,F;
[L2:B2]'REAL' BXI,BXIP,BYI,BYIP,CI,CIP,FI,FIP;
XS:= X0+L1*H; 'REF' []'REAL' UIP= NET[L1,@L2];
'FOR'J 'FROM'L2 'TO' B2
'DO' LOWER ORDER
(BXIP[J],BYIP[J],CIP[J],FIP[J],XS,Y0+J*K,UIP[J])
'OD';
'FOR' IP 'FROM' L1+1 'TO' B1
'DO' 'REF' []'REAL' UI = NET[IP-1,@L2],
UIP= NET[IP ,@L2];
XN:= XS; XS := X0 + IP*H;
BXI:= BXIP; BYI:= BYIP; CI:=CIP; FI:= FIP;
'INT' JJ:= L2;
'REAL' U1,U2:= UI [L2],
U3,U4:= UIP[L2];
YE:= Y0+L2*K;
LOWER ORDER (BXIP[L2],BYIP[L2],CIP[L2],FIP[L2],XS,YE,U4);
'FOR' JP 'FROM' L2+1 'TO' B2
'DO' YW:= YE; YE := Y0 + JP*K;
U1:= U2; U2 := UI [JP];
U3:= U4; U4 := UIP[JP];
LOWER ORDER (BXIP[JP],BYIP[JP],CIP[JP],FIP[JP],XS,YE,U4);
BX:= (BXI[JJ],BXI[JP],BXIP[JJ],BXIP[JP]);
BY:= (BYI[JJ],BYI[JP],BYIP[JJ],BYIP[JP]);
C := (CI [JJ],CI [JP],CIP [JJ],CIP [JP]);
F := (FI [JJ],FI [JP],FIP [JJ],FIP [JP]);
'INT' LL:= 0;
'FOR' TRIANGLE 'TO' 2
'DO' 'CASE' TRIANGLE
'IN' (LB:= 1; UB:= 3; PX:= PXT1; PY:= PYT1;
PRINCIPLE PART(AXX,AXY,AYX,AYY,XN+H3,YW+K3);
BBX:= (BX[1]+BX[2]+BX[3])/3;
BBY:= (BY[1]+BY[2]+BY[3])/3 ),
(LB:= 2; UB:= 4; PX:= PXT2; PY:= PYT2;
PRINCIPLE PART(AXX,AXY,AYX,AYY,XS-H3,YE-K3);
BBX:= (BX[2]+BX[3]+BX[4])/3;
BBY:= (BY[2]+BY[3]+BY[4])/3 )
'ESAC';
KK:= HB(AXX,AXY,AYX,AYY,BBX,BBY,H,K);
AXX*:=3; AXY*:=3; AYX*:=3; AYY*:=3;
'INT' EI:= 0;
'FOR' II 'FROM' IP-1 'TO' IP 'DO'
'FOR' JJ 'FROM' JP-1 'TO' JP 'DO'
'IF' EI+:= 1; EI = (TRIANGLE!4,1)
'THEN' 'SKIP'
'ELSE' FF:= Q2 * F[EI];
'FOR' EJ 'FROM' LB 'TO' UB
'DO' LL+:= 1;
PE := KK*(BX[EJ]*PX[EI]+BY[EJ]*PY[EI]);
FF+:= (Q0+PE)*F[EJ];
STIFF[II,JJ,ORDER[LL]] +:=
(PX[EJ] *(AXX*PX[EI]+AXY*PY[EI] + BX[EI])+
PY[EJ] *(AYX*PX[EI]+AYY*PY[EI] + BY[EI])+
C[EJ]*( (EI=EJ!Q1!Q0) + PE ) )*WEIGHT
'OD';
FORCE[II,JJ] +:= FF*WEIGHT
'FI'
'OD' 'OD'
'OD';
JJ:= JP
'OD' 'OD';
BOUNDARY CONDITIONS ( PROBLEM,U,FORCE,STIFF );
RHS:= FORCE; JAC:= STIFF
'END'
#$E#;
#E$#
# DIFFERENT OPTIONS FOR DISCRETIZATION # 'PR' EJECT 'PR'
# BROOKS AND HUGHES ANISOTROPIC DIFFUSION #
#
'REAL' Q0:= 0, Q1:= 1, Q2:= 1;
'PROC' LUMP = ('BOOL' B)'VOID':
( B ! Q0:=0; Q1:=1; Q2:=1 ! Q0:=0.25; Q1:=0.50; Q2:=0.25 );
'REAL' HB FACTOR:= 1/3;
'PROC' HB:= ('REF''REAL' AXX,AXY,AYX,AYY,'REAL' BX,BY,H,K)'REAL':0.0;
#
'PROC' KAPPA = ('REF''REAL' AXX,AXY,AYX,AYY, 'REAL' BX,BY,H,K)'REAL':
#$B HBW B$#
'BEGIN' # ## ADDS ARTIFICIAL DIFFUSION AND COMPUTES K # ##
'REAL' BEB = (BX*AXX+BY*AYX)*BX + 1.0E-100 +
(BX*AXY+BY*AYY)*BY ,
BB = BX*BX + BY*BY ,
BH = SQRT(BX*BX*H*H + BY*BY*K*K);
'REAL' KK = ( BB=0 ! 0.0 ! BH*HBFACTOR*M(BH*BB/BEB)/BB );
AXX+:= KK*BX*BX; AXY+:= KK*BX*BY;
AYX+:= KK*BX*BY; AYY+:= KK*BY*BY; KK
'END'
#$E#;
#E$#
'PROC' M = ('REAL' A)'REAL':
#$B ILINM B$#
'IF' 'REAL' X, W:= 'ABS'A; W< 0.2
'THEN' W*:= W; (((( W - 9.9 ) * W + 99.0 ) * W
- 1039.5 ) * W + 15592.5) * A /46777.5
'ELSE' X:= ( W-1.0)/W;
( W< 18.0 ! X +:= 2.0/(EXP(W+W) - 1.0));
( A > 0.0 ! X ! -X )
'FI'
#$E#;
#E$#
#
---------------------->
! J Y
! I STAR ORIENTATION:
! 1 -- N -- 2 -----------------> Y
! ! / ! ! J
! W / E !
! ! / ! ! -3 -2
! 3 -- S -- 4 ! -1 0 1
! X ! 2 3
V X ! I
V
#
# FINITE ELEMENT METHOD (S-U P-G FEM) # 'PR' EJECT 'PR'
# FEM DISCRETIZATION FOR MG /PWH800925 / VSN.810311 #
# PIECEWISE CONSTANT COEFFICIENTS VSN.830107 #
# INCL. STREAMLINE-UPWIND PETROV-GALERKIN #
'PROC' FEMD = ('PROBLEM' PROBLEM, 'GRID' U,
'REF''NET' RHS, 'REF''NETMAT' JAC )'VOID':
#$B FEMW B$#
'BEGIN' 'REF'[,]'REAL' NET = NET'OF' U;
'INT' L1 = 1'LWB' NET, L2 = 2'LWB'NET,
B1 = 1'UPB' NET, B2 = 2'UPB'NET;
'HEAP' [L1:B1,L2:B2,-3:3 ]'REAL' STIFF; 'ZERO' STIFF;
'HEAP' [L1:B1,L2:B2]'REAL' FORCE; 'ZERO' FORCE;
'REAL' H = XH0/(2**(LEVEL'OF'U)),
K = YH0/(2**(LEVEL'OF'U));
'PROC'('REF''REAL','REF''REAL','REF''REAL','REF''REAL',
'REAL','REAL' )'VOID' PRINCIPLE PART
= PRINCIPLE PART 'OF' PROBLEM;
'PROC'('REF''REAL','REF''REAL','REF''REAL','REF''REAL',
'REAL','REAL','REAL')'VOID' LOWER ORDER
= LOWER ORDER 'OF' PROBLEM;
'BOOL' LIN = LINEAR'OF'PROBLEM;
[]'INT' ORDER = (0,1,3,-1,0,2,-3,-2,0,
0,2,3,-2,0,1,-3,-1,0);
'REAL' H3 = H/3, K3 = K/3;
[]'REAL' PXT1 = (-1/H, 0 ,+1/H, 0), PXT2 = ( 0,-1/H, 0 ,+1/H),
PYT1 = (-1/K,+1/K, 0 , 0), PYT2 = ( 0, 0 ,-1/K,+1/K);
'INT' LB,UB;
'REAL' XN,YW, AXX,AXY,AYX,AYY, WEIGHT:= H*K/6;
'REAL' XX,YY,UU, BX ,BY ,C ,F , KK,PE;
[1:4]'REAL' PX,PY;
'FOR' I 'FROM' L1 'TO' B1-1
'DO' XN:= X0 + I*H;
'REF' []'REAL' UI = NET[I ,@L2],
UIP= NET[I+1,@L2];
'REAL' U1,U2:= UI [L2],
U3,U4:= UIP[L2];
'FOR' J 'FROM' L2 'TO' B2-1
'DO' YW:= Y0 + J*K;
U1:= U2; U2 := UI [J+1];
U3:= U4; U4 := UIP[J+1];
'INT' LL:= 0;
'FOR' TRIANGLE 'TO' 2
'DO' 'CASE' TRIANGLE
'IN' (LB:= 1; UB:= 3;
PX:= PXT1; PY:= PYT1;
XX:= XN+H3; YY:= YW+K3;
(LIN!'SKIP'!UU:= (U1+U2+U3)/3) ),
(LB:= 2; UB:= 4;
PX:= PXT2; PY:= PYT2;
XX:= XX+H3; YY:= YY+K3;
(LIN!'SKIP'!UU:= (U2+U3+U4)/3) )
'ESAC';
PRINCIPLE PART(AXX,AXY,AYX,AYY,XX,YY);
LOWER ORDER (BX,BY,C,F,XX,YY,UU);
KK:= HB(AXX,AXY,AYX,AYY,BX,BY,H,K);
AXX*:=3; AXY*:=3; AYX*:=3; AYY*:=3;
'INT' EI:= 0;
'FOR' II 'FROM' I 'TO' I+1 'DO'
'FOR' JJ 'FROM' J 'TO' J+1 'DO'
'IF' EI+:= 1; EI = (TRIANGLE!4,1)
'THEN' 'SKIP'
'ELSE' PE := KK*(BX*PX[EI]+BY*PY[EI]);
'FOR' EJ 'FROM' LB 'TO' UB
'DO' LL+:= 1;
STIFF[II,JJ,ORDER[LL]] +:=
(PX[EJ] *(AXX*PX[EI]+AXY*PY[EI] + BX)+
PY[EJ] *(AYX*PX[EI]+AYY*PY[EI] + BY)+
C*( (EI=EJ!0.5!0.25) + PE ) )*WEIGHT
'OD';
FORCE[II,JJ] +:= (1.0+3.0*PE) *F*WEIGHT
'FI'
'OD' 'OD'
'OD' 'OD' 'OD';
BOUNDARY CONDITIONS ( PROBLEM,U,FORCE,STIFF );
RHS:= FORCE; JAC:= STIFF
'END'
#$E#;
#E$#
# BOUNDARY CONDITION OPTIONS # 'PR' EJECT 'PR'
'PROC' BC PW LIN = ('PROBLEM' PROBLEM, 'GRID' U,
'NET' FORCE, 'NETMAT' STIFF )'VOID':
#$B BCPL B$#
'BEGIN' 'REF'[,]'REAL' NET = NET'OF' U;
'INT' L1 = 1'LWB' NET, L2 = 2'LWB'NET,
B1 = 1'UPB' NET, B2 = 2'UPB'NET;
'REAL' H = XH0/(2**(LEVEL'OF'U)), SHIFT = 0.001,
K = YH0/(2**(LEVEL'OF'U));
'REF'[,]'REAL' OMEGA = OMEGA 'OF' PROBLEM ;
'PROC'('REF''REAL','REF''REAL','REF''REAL','REF''REAL',
'REAL','REAL' )'VOID' PRINCIPLE PART
= PRINCIPLE PART 'OF' PROBLEM;
'PROC'('REF''REAL',
'REAL','REAL' )'VOID' BOUNDARY PP
= BOUNDARY PP 'OF' PROBLEM;
'PROC'('REF''REAL','REF''REAL',
'REAL','REAL','REAL')'VOID' BOUNDARY LO
= BOUNDARY LO 'OF' PROBLEM;
'PROC' BOUNDARY SEGMENT = ('REF' 'INT' I,J, 'INT' DIRECTION,
'REF' 'REAL' X,Y,BZ,CZ)'VOID':
'BEGIN' 'INT' IN,JN;
'REAL' XN,YN,AZN,BZN,CZN,HSTEP,XHALF,YHALF,P,Q,
A11,A12,A21,A22,PPP,PLO;
'INT' AD = 'ABS' DIRECTION, ND= 'SIGN' DIRECTION;
( AD = 1 ! IN := I ; XN:= X
! IN := I + ND; XN:= X0 + IN*H; HSTEP:= H/2 );
( AD = 3 ! JN := J ; YN:= Y
! JN := J + ( AD ! ND, -ND ! ERROR;0 );
YN:= Y0 + JN*K; HSTEP:= K/2 );
( AD = 2 ! HSTEP:= 0.5*SQRT(H*H+K*K) );
XHALF:= (X+XN)*0.5 - (YN-Y)*SHIFT;
YHALF:= (Y+YN)*0.5 + (XN-X)*SHIFT;
PRINCIPLE PART (A11,A12,A21,A22,XHALF,YHALF);
BOUNDARY PP (AZN, XHALF,YHALF);
BOUNDARY LO (BZN,CZN, XN, YN, NET[IN,JN]);
'CASE' AD
'IN' ( P:= A11; Q:= A12 ),
( P:= (A11+A12+A21+A22)/2; Q:=(-A11+A12-A21+A22)/2 ),
( P:= A22; Q:= -A21 )
'ESAC';
PPP:= ( P*AZN - Q)/2; PLO:= P*HSTEP;
STIFF[ I, J, 0 ] +:= -PPP + PLO*BZ ;
STIFF[ I, J, DIRECTION] +:= PPP ;
FORCE[ I, J ] +:= PLO*CZ ;
STIFF[IN,JN, 0 ] +:= +PPP + PLO*BZN;
STIFF[IN,JN,-DIRECTION] +:= -PPP ;
FORCE[IN,JN ] +:= PLO*CZN;
I:= IN; J:= JN; X:= XN; Y:= YN; BZ:= BZN; CZ:= CZN
'END';
'IF' OMEGA :=: 'NIL'
'THEN' 'INT' I:= L1, J:= L2; 'REAL' X:= X0+I*H, Y:= Y0+J*K, BZ,CZ;
BOUNDARY LO (BZ,CZ,X,Y,NET[I,J]);
'WHILE' I<B1 'DO' BOUNDARY SEGMENT(I,J, 3,X,Y,BZ,CZ) 'OD';
'WHILE' J<B2 'DO' BOUNDARY SEGMENT(I,J, 1,X,Y,BZ,CZ) 'OD';
'WHILE' I>L1 'DO' BOUNDARY SEGMENT(I,J,-3,X,Y,BZ,CZ) 'OD';
'WHILE' J>L2 'DO' BOUNDARY SEGMENT(I,J,-1,X,Y,BZ,CZ) 'OD'
'ELIF' 1'UPB'OMEGA >1
'THEN' 'INT' I:= 'ROUND'((OMEGA[1,1]-X0)/H ),
J:= 'ROUND'((OMEGA[1,2]-Y0)/K ),IN,JN,DIR;
'REAL' X:= X0+I*H, Y:= Y0+J*K, BZ,CZ;
BOUNDARY LO (BZ,CZ,X,Y,NET[I,J]);
'FOR' S 'FROM' 2 'TO' 1'UPB'OMEGA
'DO' IN:= 'ROUND'((OMEGA[S,1]-X0)/H );
JN:= 'ROUND'((OMEGA[S,2]-Y0)/K );
DIR:= 'SIGN'(JN-J) + 3*'SIGN'(IN-I);
'WHILE' BOUNDARY SEGMENT(I,J,DIR,X,Y,BZ,CZ);
I /= IN 'OR' J /= JN
'DO' 'SKIP' 'OD'
'OD'
'FI'
'END'
#$E#;
#E$#
'PR' EJECT 'PR'
'PROC' BC PW CON = ('PROBLEM' PROBLEM, 'GRID' UG,
'NET' FORCE, 'NETMAT' STIFF )'VOID':
#$B BCPC B$#
'BEGIN' 'REF'[,]'REAL' NET = NET'OF' UG;
'INT' L1 = 1'LWB' NET, L2 = 2'LWB'NET,
B1 = 1'UPB' NET, B2 = 2'UPB'NET;
'REAL' H = XH0/(2**(LEVEL'OF'UG)),
K = YH0/(2**(LEVEL'OF'UG));
'PROC'('REF''REAL','REF''REAL',
'REAL','REAL','REAL')'VOID' BOUNDARY LO
= BOUNDARY LO 'OF' PROBLEM;
'PROC' BOUNDARY SEGMENT = ('INT' DIRECTION )'VOID':
'BEGIN' 'REAL' XO:= X, YO:= Y, UO:= U, B, C, W;
'INT' IO:= I, JO:= J, ND:= 'SIGN' DIRECTION;
W:= ( 'ABS' DIRECTION = 1
! J+:= ND; Y:= Y0+J*K; K/2
! I+:= ND; X:= X0+I*H; H/2
);
BOUNDARY LO (B,C,(X+XO)/2,(Y+YO)/2,
(LIN! 'SKIP' ! ((U:=NET[I,J]) + UO)/2 ));
# ## DIRICHLET: B=LARGE, C=LARGE*BC # ##
# ## NEUMANN : B= 0.0, C= H # ##
STIFF[IO,JO, 0] +:= W*B;
FORCE[IO,JO ] +:= W*C;
STIFF[I ,J , 0] +:= W*B;
FORCE[I ,J ] +:= W*C
'END';
'BOOL' LIN = LINEAR'OF'PROBLEM;
'INT' I:= L1, J:= L2;
'REAL' X:= X0+I*H, Y:= Y0+J*K, U:= NET[I,J];
'WHILE' I<B1 'DO' BOUNDARY SEGMENT( 3) 'OD';
'WHILE' J<B2 'DO' BOUNDARY SEGMENT( 1) 'OD';
'WHILE' I>L1 'DO' BOUNDARY SEGMENT(-3) 'OD';
'WHILE' J>L2 'DO' BOUNDARY SEGMENT(-1) 'OD'
'END'
#$E#;
#E$#
'PR' EJECT 'PR'
'PROC' BC DIRICH = ('PROBLEM' PROBLEM, 'GRID' U,
'NET' FORCE, 'NETMAT' STIFF )'VOID':
#$B BCDIR B$#
'BEGIN' 'REF'[,]'REAL' NET = NET'OF' U;
'INT' L1 = 1'LWB' NET, L2 = 2'LWB'NET,
B1 = 1'UPB' NET, B2 = 2'UPB'NET;
'REAL' H = XH0/(2**(LEVEL'OF'U)),
K = YH0/(2**(LEVEL'OF'U));
'PROC'('REF''REAL','REF''REAL',
'REAL','REAL','REAL')'VOID' BOUNDARY LO
= BOUNDARY LO 'OF' PROBLEM;
'PROC' BOUNDARY SEGMENT = ('INT' DIRECTION )'VOID':
'BEGIN' 'INT' ND = 'SIGN' DIRECTION;
( 'ABS' DIRECTION = 1
! J+:= ND; Y:= Y0+J*K
! I+:= ND; X:= X0+I*H
);
# ## DIRICHLET: B=LARGE, C=LARGE*BC # ##
BOUNDARY LO ( STIFF[I,J,0], FORCE[I,J], X, Y,
( LIN ! 'SKIP' ! NET[I,J]) )
'END';
'BOOL' LIN = LINEAR'OF'PROBLEM;
'INT' I:= L1, J:= L2;
'REAL' X:= X0+I*H, Y:= Y0+J*K;
BOUNDARY LO ( STIFF[I,J,0], FORCE[I,J], X, Y,
( LIN ! 'SKIP' ! NET[I,J]) );
'WHILE' I<B1 'DO' BOUNDARY SEGMENT( 3) 'OD';
'WHILE' J<B2 'DO' BOUNDARY SEGMENT( 1) 'OD';
'WHILE' I>L1 'DO' BOUNDARY SEGMENT(-3) 'OD';
'WHILE' J>L2 'DO' BOUNDARY SEGMENT(-1) 'OD'
'END'
#$E#;
#E$#
# GALERKIN OPERATOR CONSTRUCTION # 'PR' EJECT 'PR'
'OP' 'RAP' = ('NETMAT' AFI )'NETMAT':
#$B RAP B$#
'BEGIN' 'INT' L1 = (1'LWB'AFI)'OVER'2, U1 = (1'UPB'AFI)'OVER'2,
L2 = (2'LWB'AFI)'OVER'2, U2 = (2'UPB'AFI)'OVER'2;
'HEAP' [L1:U1,L2:U2,-3:3]'REAL' ACO; 'ZERO' ACO;
'INT' TI,TK,TIM,TKM;
[1:3,1:3,-3:3]'REAL' FINE;
'REF'[]'REAL'
A = FINE[1,1,@-3], B = FINE[1,2,@-3], C = FINE[1,3,@-3],
D = FINE[2,1,@-3], E = FINE[2,2,@-3], F = FINE[2,3,@-3],
G = FINE[3,1,@-3], H = FINE[3,2,@-3], J = FINE[3,3,@-3];
'REF''REAL'
AA =A[ 0], AB =A[ 1], AD =A[ 3],
BA =B[-1], BB =B[ 0], BC =B[ 1], BD =B[ 2], BE =B[ 3],
CB =C[-1], CC =C[ 0], CE =C[ 2], CF =C[ 3],
DA =D[-3], DB =D[-2], DD =D[ 0], DE =D[ 1], DG =D[ 3],
EB =E[-3], EC =E[-2], ED =E[-1], EE =E[ 0],
EF =E[ 1], EG =E[ 2], EH =E[ 3],
FC =F[-3], FE =F[-1], FF =F[ 0], FH =F[ 2], FJ =F[ 3],
GD =G[-3], GE =G[-2], GG =G[ 0], GH =G[ 1],
HE =H[-3], HF =H[-2], HG =H[-1], HH =H[ 0], HJ =H[ 1],
JF =J[-3], JH =J[-1], JJ =J[ 0];
TI:= 2*L1; TK:= 2*L2;
FINE[3,3,]:= AFI[TI,TK,];
# ##JJ# ##ACO[L1,L2,0]:= 4*JJ;
'FOR' K 'FROM' L2+1 'TO' U2
'DO' TKM:= TK; TK:= 2*K;
FINE[3,1:3,]:= AFI[TI,TKM:TK,];
# ##GG# ##ACO[L1,K-1,0] +:= (GH+HG)*2 + HH;
# ##JJ# ##ACO[L1,K ,0] +:= (JH+HJ)*2 + HH +JJ*4;
# ##GJ# ##ACO[L1,K-1,1] +:= (GH+HJ)*2 + HH;
# ##JG# ##ACO[L1,K ,-1] +:= (HG+JH)*2 + HH
'OD';
'FOR' I 'FROM' L1+1 'TO' U1
'DO' TIM:= TI; TI:= 2*I; TK:= 2*L2;
FINE[1:3,3,]:= AFI[TIM:TI,TK,];
# ##CC# ##ACO[I-1,L2,0] +:= (CF+FC)*2 + FF;
# ##JJ# ##ACO[I ,L2,0] +:= (JF+FJ)*2 + FF + JJ*4;
# ##CJ# ##ACO[I-1,L2,3] +:= (CF+FJ)*2 + FF;
# ##JC# ##ACO[I ,L2,-3] +:= (FC+JF)*2 + FF;
'FOR' K 'FROM' L2+1 'TO' U2
'DO' TKM:= TK; TK:= 2*K;
FINE[1:3,1:3,]:= AFI[TIM:TI,TKM:TK,];
'REF'[]'REAL' A = ACO[I-1,K-1,@-3],
C = ACO[I-1,K ,@-3],
G = ACO[I ,K-1,@-3],
J = ACO[I ,K ,@-3];
# ##AA# ## A[ 0]+:= BD + DB;
# ##CC# ## C[ 0]+:= (CE+EC+CF+FC)*2 + EE+FF+BE+EB+EF+FE;
# ##GG# ## G[ 0]+:= (GE+EG+GH+HG)*2 + EE+HH+DE+ED+EH+HE;
# ##JJ# ## J[ 0]+:= JJ*4+ (JF+FJ+JH+HJ)*2 + FF+HH+FH+HF;
# ##AC# ## A[ 1]+:= BE+DB+DE;
# ##CA# ## C[-1]+:= EB+BD+ED;
# ##AG# ## A[ 3]+:= BD+DE+BE;
# ##GA# ## G[-3]+:= DB+ED+EB;
# ##GC# ## G[-2]+:= (GE+EC)*2 + EE+HE+DE+HF+DB+EF+EB;
# ##CG# ## C[ 2]+:= (EG+CE)*2 + EE+EH+ED+FH+BD+FE+BE;
# ##GJ# ## G[ 1]+:= (GH+HJ)*2 + HH+EH+HF+EF;
# ##JG# ## J[-1]+:= (HG+JH)*2 + HH+HE+FH+FE;
# ##CJ# ## C[ 3]+:= (CF+FJ)*2 + FF+EH+EF+FH;
# ##JC# ## J[-3]+:= (FC+JF)*2 + FF+HE+FE+HF
'OD' 'OD';
'FOR' I 'FROM' L1 'TO' U1 'DO'
'FOR' K 'FROM' L2 'TO' U2 'DO'
'FOR' L 'FROM' -3 'TO' 3 'DO' ACO[I,K,L] *:= 0.25
'OD' 'OD' 'OD'; ACO
'END'
#$E#;
#E$#
# ORIENTATION:
ACO = COARSE K-1 K
--------------------------> Y
!
! FINE 1 2 3
! -3 -2
! I-1 1 A -- B -- C ! /
! ! / ! / ! -1 - 0 - 1
! 2 D -- E -- F / !
! ! / ! / ! 2 3
! I 3 G -- H -- J
X V
THIS ALGORITHM CAN BE EFFECTUATED BY 87 ADDITIONS AND
7 MULTIPLICATIONS BY 0.25
PER COARSE GRID POINT
#
# SOLVE BY GAUSSIAN ELIMINATION # 'PR' EJECT 'PR'
'PROC' SOLVE = ('NETMAT' C, 'NET' U,F)'VOID':
#$B SOLVE B$#
'BEGIN' 'NETMAT' B = C [@1,@1,@3'LWB'C];
'NET' RHS = F [@1,@1];
'INT' N1 = 1'UPB'RHS, N2= 2'UPB'RHS;
'INT' N = N1*N2, S = 3'UPB'B, N3 = N2-3;
# ## S=4 : NINE-STAR ; S=3 : SEVEN-STAR # ##
'PRIO' <> = 4;
'OP' <> = ('VEC' A,B) 'VOID':
'FOR' I 'FROM' 'LWB' A 'TO' 'UPB' A
'DO' 'REAL' S = A[I]; A[I]:= B[I]; B[I]:= S 'OD';
'REF'[,]'REAL' SOLUT = U[@1,@1];
[1:N, 1:N+1]'REAL' A;
'VEC' V = A[,N+1];
'FOR' I 'TO' N1 'DO'
'FOR' J 'TO' N2 'DO'
'BEGIN' 'INT' K = (I-1)*N2 + J;
'REF'[]'REAL' AK= A[K,], BIJ= B[I,J,@-S];
'FOR' L 'TO' N 'DO' AK[L]:= 0.0 'OD';
'FOR' Z 'FROM' -S 'TO' S
'DO' 'IF' I= 1 'AND' Z<-1 'THEN' 'SKIP'
'ELIF' I=N1 'AND' Z> 1 'THEN' 'SKIP'
'ELIF' J= 1 'AND' (Z=-4 'OR' Z=-1 'OR' Z= 2)
'THEN' 'SKIP'
'ELIF' J=N2 'AND' (Z= 4 'OR' Z= 1 'OR' Z=-2)
'THEN' 'SKIP'
'ELSE'
AK[(Z<-1 ! K+Z-N3 !: Z>1 ! K+Z+N3 ! K+Z)]
:= BIJ[Z]
'FI'
'OD';
AK[N+1] := RHS[I,J]
'END' 'OD' 'OD';
# ## 'FOR' K 'TO' N
'DO' PRINT(NEWLINE);
'FOR' L 'TO' N+1
'DO' PRINT((FLOAT(A[K,L],7,1,2)," "))
'OD' 'OD';PRINT((NEWLINE,NEWLINE));
# ##
'PR' EJECT 'PR'
'FOR' J 'TO' N
'DO' 'INT' JP1= J+1; 'INT' PJ:= J;
'REAL' SI,S:= 'ABS' A[J,J];
'FOR' I 'FROM' JP1 'TO' N
'DO' ((SI:='ABS'A[I,J]) >S ! S:=SI; PJ:=I ) 'OD';
(J /= PJ ! A[PJ,] <> A[J,] ); S := A[J,J];
'FOR' I 'FROM' JP1 'TO' N
'DO' SI:= A[I,J]/S;
'FOR' K 'FROM' J 'TO' N+1
'DO' A[I,K] -:= A[J,K]*SI 'OD'
'OD'
'OD';
'FOR' J 'FROM' N 'BY' -1 'TO' 1
'DO' V[J] /:= A[J,J];
'FOR' I 'FROM' J-1 'BY' -1 'TO' 1
'DO' V[I] -:= A[I,J]*V[J] 'OD'
'OD';
'FOR' I 'TO' N1
'DO' 'FOR' J 'TO' N2
'DO' SOLUT[I,J]:=V[(I-1)*N2 + J] 'OD'
'OD'
'END'
#$E#;
#E$#
# OPERATOR EVALUATION # 'PR' EJECT 'PR'
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #
# EVERYWHERE THE MATRIX DOES NOT DEFINE THE NETMAT M, #
# !!!!!!!!!!! M SHOULD CONTAIN ZEROES !!!!!!!!!!!!!!! #
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #
'PROC' RESIDUAL = ('NETMAT' M, 'NET' U,F )'NET':
#$B RESID B$#
'BEGIN''INT' L1= 1'LWB'U, L2= 2'LWB'U,
U1= 1'UPB'U, U2= 2'UPB'U;
'HEAP'[L1:U1,L2:U2]'REAL' S;
'REF'[ ]'REAL' UIM:= U[L1,@L2], UI, UIP:= U[L1,@L2];
'FOR' I 'FROM' L1 'TO' U1
'DO' ( UI:= UIP; I = U1 ! 'SKIP' ! UIP:= U[I+1,@L2] );
'REF' []'REAL' SI = S[I,@L2], FI = F[I,@L2];
'REF'[,]'REAL' MI = M[I,@L2,@-3];
'INT' JM:= L2, JJ, JP:= L2;
'FOR' J 'FROM' L2 'TO' U2
'DO' ( JJ:= JP; J=U2 ! 'SKIP' ! JP+:= 1 );
'REF'[]'REAL' MIJ = MI[JJ,@-3];
SI[JJ]:= FI[JJ] -
( MIJ[-3]*UIM[JJ] + MIJ[-2]*UIM[JP] +
MIJ[-1]*UI [JM] + MIJ[ 0]*UI [JJ] + MIJ[ 1]*UI [JP] +
MIJ[ 2]*UIP[JM] + MIJ[ 3]*UIP[JJ] );
JM:= JJ
'OD';
UIM:= UI
'OD';S
'END'
#$E#;
#E$#
'PR' EJECT 'PR'
'OP' * = ('NETMAT' M, 'NET' U)'NET':
#$B NETMAT* B$#
'BEGIN' 'INT' L1= 1'LWB'U, L2= 2'LWB'U, L3= 3'LWB'M,
U1= 1'UPB'U, U2= 2'UPB'U;
'HEAP'[L1:U1,L2:U2]'REAL' S;
'REF'[ ]'REAL' UIM:= U[L1,@L2], UI, UIP:= U[L1,@L2];
'FOR' I 'FROM' L1 'TO' U1
'DO' ( UI:= UIP; I = U1 ! 'SKIP' ! UIP:= U[I+1,@L2] );
'REF' []'REAL' SI = S[I,@L2];
'REF'[,]'REAL' MI = M[I,@L2,@L3];
'INT' JM:= L2, JJ, JP:= L2;
'FOR' J 'FROM' L2 'TO' U2
'DO' ( JJ:= JP; J=U2 ! 'SKIP' ! JP+:= 1 );
'REF'[]'REAL' MIJ = MI[JJ,@L3];
'REF' 'REAL' SIJ = SI[JJ]:=
+ MIJ[-3]*UIM[JJ] + MIJ[-2]*UIM[JP] +
MIJ[-1]*UI [JM] + MIJ[ 0]*UI [JJ] + MIJ[ 1]*UI [JP] +
MIJ[ 2]*UIP[JM] + MIJ[ 3]*UIP[JJ] ;
JM:= JJ
'OD';
UIM:= UI
'OD';S
'END'
#$E#;
#E$#
# RELAXATION METHODS # 'PR' EJECT 'PR'
'PROC' DIRICHLET RELAX = ('PROBLEM' P, 'REF''GRID' UG,
'NETMAT' M, 'NET' F)'VOID':
# SUBSTITUTES BOUNDARY VALUES,
USEFULL IN CASE OF DIRICHLET BOUNDARY CONDITIONS #
#$B RELAXD B$#
'IF' 'REF'[, ]'REAL' OMEGA = OMEGA'OF'P;
'REF'[, ]'REAL' U = NET'OF'UG;
OMEGA :=: 'NIL'
'THEN' 'INT' L1=1'LWB'U, U1=1'UPB'U, L2=2'LWB'U, U2=2'UPB'U;
'FOR' I 'FROM' L1 'BY' U1-L1 'TO' U1
'DO' 'FOR' J 'FROM' L2 'TO' U2
'DO' U[I,J]:= F[I,J]/M[I,J,0] 'OD'
'OD';
'FOR' J 'FROM' L2 'BY' U2-L2 'TO' U2
'DO' 'FOR' I 'FROM' L1 'TO' U1
'DO' U[I,J]:= F[I,J]/M[I,J,0] 'OD'
'OD'
'ELSE' 'INT' L = LEVEL'OF'UG;
'REAL' H = XH0/2**L, K = YH0/2**L;
'INT' I:= 'ROUND'((OMEGA[1,1]-X0)/H),
J:= 'ROUND'((OMEGA[1,2]-Y0)/K), IN,JN,DI,DJ;
'FOR' S 'FROM' 2 'TO' 1'UPB'OMEGA
'DO' IN:= 'ROUND'((OMEGA[S,1]-X0)/H); DI:= 'SIGN'(IN-I);
JN:= 'ROUND'((OMEGA[S,2]-Y0)/K); DJ:= 'SIGN'(JN-J);
'WHILE' U[I,J]:= F[I,J]/M[I,J,0];
I /= IN 'AND' J /= JN
'DO' I+:= DI; J+:=DJ 'OD'
'OD'
'FI'
#$E#;
#E$#
# POINT RELAXATION PROCEDURES # 'PR' EJECT 'PR'
#'BOOL' SYMMETRIC:= 'FALSE', BACKWARD:= 'FALSE',
REVERSE := 'FALSE', ZEBRA := 'FALSE';
#
'PROC' PGS RELAX = ('NETMAT'M, 'NET' U,F)'VOID':
#$B RELPGS B$#
'BEGIN' # ## POINT GAUSS SEIDEL (PGS) # ##
'INT' L1:= 1'LWB'U, U1:= 1'UPB'U, START1, STEP1, STOP1,
L2:= 2'LWB'U, U2:= 2'UPB'U, START2, STEP2, STOP2;
'TO' ( SYMMETRIC ! 2 ! 1)
'DO' ( BACKWARD ! START1:= U1; STEP1:= -1; STOP1:= L1
! START1:= L1; STEP1:= 1; STOP1:= U1 );
( REVERSE ! START2:= U2; STEP2:= -1; STOP2:= L2
! START2:= L2; STEP2:= 1; STOP2:= U2 );
'FOR' I 'FROM' START1 'BY' STEP1 'TO' STOP1
'DO' 'REF'[] 'REAL' FI= F[I,@L2], UIM= U[(I>L1!I-1!I),@L2],
UI= U[I,@L2], UIP= U[(I<U1!I+1!I),@L2];
'REF'[,]'REAL' MI= M[I,@L2,@-3];
'FOR' J 'FROM' START2 'BY' STEP2 'TO' STOP2
'DO' 'INT' JM= (J>L2!J-1!J), JP= (J<U2!J+1!J);
'REF'[]'REAL' MIJ = MI[J,@-3];
UI[J]:= ( MIJ[-3]*UIM[J]+MIJ[-2]*UIM[JP]+
MIJ[-1]*UI [JM] - FI[J]+MIJ[ 1]*UI [JP]+
MIJ[ 2]*UIP[JM]+MIJ[ 3]*UIP[J] )/ -MIJ[ 0]
'OD'
'OD' ;
( SYMMETRIC! REVERSE:='NOT'REVERSE; BACKWARD:='NOT'BACKWARD)
'OD'
'END'
#$E#;
#E$#
# LINE RELAXATION PROCEDURES # 'PR' EJECT 'PR'
'PROC' LGS RELAX = ('NETMAT'M, 'NET' U,F)'VOID':
#$B RELLGS B$#
'BEGIN' # ## LINE GAUSS SEIDEL (LGS) # ##
'INT' ST = ( ZEBRA ! 2 ! 1 );
'INT' L1:= 1'LWB'U, U1:= 1'UPB'U, START, STEP, STOP;
'FOR' K 'TO' ( SYMMETRIC 'OR' ZEBRA ! 2 ! 1)
'DO' ( BACKWARD ! START := U1; STEP := -ST; STOP := L1
! START := L1; STEP := ST; STOP := U1 );
( ZEBRA
! ( SYMMETRIC /= 'ODD'(K+START) ! START+:= 'SIGN'STEP )
# ## ( SYMMETRIC ! EVEN-ODD ! ODD-EVEN ) HALF STEP # ##);
'FOR' I 'FROM' START 'BY' STEP 'TO' STOP
'DO' LINE RELAX ( U[ (I>L1!I-1!I),], U[I,],
U[ (I<U1!I+1!I),], F[I,], M[I,,@-3] )
'OD';
(SYMMETRIC ! BACKWARD:= 'NOT' BACKWARD )
'OD'
'END'
#$E#;
#E$#
'PROC' LINE RELAX = ( 'VEC'UM,U,UP,F,'REF'[,]'REAL'M )'VOID':
#$B RELLINE B$#
'BEGIN' 'VEC' B= M[, 1],# ##B[L:K], SUP# ## N = M[,-3], NE= M[,-2],
A= M[, 0],# ##A[L:K], DIA# ## S = M[, 3], SW= M[, 2],
C= M[,-1];# ##C[L:K], SUB# ##
# ##NOT EXISTING MATRIX ELEMENTS: C[L]= B[K]= 0 !!!!!!# ##
'INT' L='LWB'F, K='UPB'F; [L:K]'REAL' AA;
( 'UPB'A /= K 'OR' 'UPB'U /= K ! ERROR );
'INT' I:=L; 'REAL' G:= 0, P; AA[L]:= 1.0;
'FOR' J 'FROM' L 'TO' K
'DO' AA[J]:= A[J] - B[I]* ( P:= C[J]/AA[I] );
G := F[J] - N[J]*UM[J] -
SW[J]*UP[I] - S[J]*UP[J] - G*P;
( J<K ! G -:= NE[J]*UM[J+1] );
U[ J]:= G; I:= J
'OD';
'FOR' J 'FROM' K 'BY' -1 'TO' L
'DO' U [J]:= G := ( U[J] - B[J]*G )/AA[J] 'OD'
'END'
#$E#;
#E$#
# DAMPED JACOBI RELAXATION # 'PR' EJECT 'PR'
'PROC' JAC RELAX = ('NETMAT'M, 'NET' U,F)'VOID':
#$B RELJAC B$#
'BEGIN' # ## APPROX.OF JACOBIAN = 2 * DIAGONAL OF JACOBIAN # ##
'INT' L1= 1'LWB'U, L2= 2'LWB'U,
U1= 1'UPB'U, U2= 2'UPB'U;
'HEAP'[L1:U1,L2:U2]'REAL' UNEW;
'REF'[ ]'REAL' UIM:= U[L1,@L2], UI, UIP:= U[L1,@L2];
'FOR' I 'FROM' L1 'TO' U1
'DO' ( UI:= UIP; I = U1 ! 'SKIP' ! UIP:= U[I+1,@L2] );
'REF' []'REAL' FI = F[I,@L2];
'REF'[,]'REAL' MI:= M[I,@L2,@-3];
'INT' JM:= L2, JJ, JP:= L2;
'FOR' J 'FROM' L2 'TO' U2
'DO' ( JJ:= JP; J=U2 ! 'SKIP' ! JP+:= 1 );
'REF'[]'REAL' MIJ = MI[JJ,@-3];
UNEW[I,JJ]:= ( MIJ[-3]*UIM[JJ] + MIJ[-2]*UIM[JP] +
MIJ[-1]*UI [JM] - MIJ[ 0]*UI [JJ] + MIJ[ 1]*UI [JP] +
MIJ[ 2]*UIP[JM] + MIJ[ 3]*UIP[JJ]
- FI [JJ])/(-2.0*MIJ[0]) ;
JM:= JJ
'OD';
UIM:= UI
'OD' ;
U:= UNEW
'END'
#$E#;
#E$#
# ILU RELAXATION # 'PR' EJECT 'PR'
'PROC' ILU RELAX = ('REF''NETMAT' M,'NETMAT' A,'NET' U,F)'VOID':
#$B RELILUX B$#
'BEGIN' ( M :=: 'NETMAT'('NIL') !
M := ILUDECOMP( A ) );
# ## THE NOT EFFICIENT VERSION !!! # ##
'NET' X = RESIDUAL (A,U,F);
'INT' L1= 1'LWB'X, L2= 2'LWB'X,
U1= 1'UPB'X, U2= 2'UPB'X;
'INT' L1P:=L1, U1M:=U1, L2P:=L2, U2M:=U2;
'REF'[]'REAL' XIM:= X[L1,@L2],XI,XIP:=X[U1,@L2];
'FOR' I 'FROM' L1P 'TO' U1M
'DO' XI:= X[I,@L2];
'REF'[,]'REAL' MI= M[I,@L2,@-3];
'INT' JM:= L2, JP:=L2P;
'FOR' J 'FROM' L2P 'TO' U2M
'DO' ( J<U2 ! JP+:=1 );
'REF'[]'REAL' MIJ= MI[J,@-3];
XI[J]:= ( MIJ[-3]*XIM[J]+MIJ[-2]*XIM[JP]+
MIJ[-1]*XI [JM] - XI [J] )/ -MIJ[0];
JM:= J
'OD' ;
XIM:= XI
'OD' ;
'FOR' I 'FROM' U1M 'BY' -1 'TO' L1P
'DO' XI:= X[I,@L2];
'REF'[,]'REAL' MI= M[I,@L2,@-3];
'INT' JM:= U2M, JP:= U2;
'FOR' J 'FROM' U2M 'BY' -1 'TO' L2P
'DO' ( J>L2 ! JM-:=1 );
'REF'[]'REAL' MIJ= MI[J,@-3];
XI[J]:= -( -XI [J]+MIJ[ 1]*XI [JP]+
MIJ[ 2]*XIP[JM]+MIJ[ 3]*XIP[J] );
JP:= J
'OD' ;
XIP:= XI
'OD' ;
'FOR' I 'FROM' L1P 'TO' U1M 'DO'
'FOR' J 'FROM' L2P 'TO' U2M 'DO'
U[I,J]+:= X[I,J] 'OD' 'OD'
'END'
#$E#;
#E$#
# PREPARATION OF ILU RELAXATION # 'PR' EJECT 'PR'
'PROC' ILU DECOMP = ('NETMAT' A )'NETMAT':
#$B ILUDEC B$#
'BEGIN' 'NETMAT' M = A[@1,@1,@-3];
'INT' II = 1'UPB'M, JJ= 2'UPB'M;
'HEAP'[1:II,1:JJ,-3:3]'REAL' AA; 'ZERO' AA;
'REF' [,]'REAL' MM3 = M[,,-3], L3 = AA[,,-3],
MM2 = M[,,-2], L2 = AA[,,-2],
MM1 = M[,,-1], L1 = AA[,,-1],
M0 = M[,, 0], L0 = AA[,, 0],
M1 = M[,, 1], U1 = AA[,, 1],
M2 = M[,, 2], U2 = AA[,, 2],
M3 = M[,, 3], U3 = AA[,, 3];
'REF' [ ]'REAL' MM11 = MM1[1,], L11 = L1[1,],
M01 = M0 [1,], L01 = L0[1,],
M11 = M1 [1,@2], U11 = U1[1,@2];
'REAL' L01JM:= L01[1]:= M01[1];
'FOR' J 'FROM' 2 'TO' JJ
'DO' L01JM := L01[J]:= M01[J] -
( L11[J]:= MM11[J] ) *
( U11[J]:= M11 [J] / L01JM # ## I.E. L01[J-1]# ## )
'OD';
'FOR' I 'FROM' 2 'TO' II
'DO''INT' IM = I-1;
'REF' [ ]'REAL'
M0I = M0 [I,], L0I = L0[I,], L0IM = L0[IM,],
MM1I = MM1[I,], L1I = L1[I,], L1IM = L1[IM,@0],
MM2I = MM2[I,], L2I = L2[I,],
MM3I = MM3[I,], L3I = L3[I,],
M1I = M1[I ,],
M1IM = M1[IM,], U1I = U1[I,], U1IM = U1[IM, ],
M2IM = M2[IM,@0], U2IM = U2[IM,@0],
M3IM = M3[IM, ], U3IM = U3[IM, ];
'FOR' J 'TO' JJ-1
'DO' L2I[J] := MM2I[J] - U1IM[J] *
( L3I[J] := MM3I[J] );
U2IM[J] := ( M2IM[J] - L1IM[J] *
( U3IM[J] := M3IM[J] / L0IM[J] ) )/L0IM[J+1]
'OD';
L3I [JJ]:= MM3I[JJ];
U3IM[JJ]:= M3IM[JJ]/ L0IM[JJ];
L0I [ 1]:= M0I [ 1] - L3I[1]*U3IM[1] - L2I[1]*U2IM[1];
'FOR' J 'FROM' 2 'TO' JJ
'DO' 'INT' JM = J - 1;
L0I[J] := M0I [J] -
( L1I[J] := MM1I[J] - L3I[J]*U2IM[JM] )*
( U1I[JM] :=(M1I[JM] - L2I[JM]*U3IM[J])/L0I[JM] )
- ( J=JJ ! 0.0 ! L2I[J]*U2IM[J] )
- L3I[J]*U3IM[J]
'OD'
'OD';
AA[@(1'LWB'A),@(2'LWB'A),@-3]
'END'
#$E#;
#E$#
# ILLU RELAXATION # 'PR' EJECT 'PR'
'PROC' SOLL = ('INT' I, 'NETMAT' DEC, 'NET' R)'VOID':
#$B SOLL B$#
'BEGIN' 'REF'[]'REAL' L = DEC[I, ,-1],
D = DEC[I, , 0],
U = DEC[I, , 1],
Z = R [I, ];
'INT' LL = 'LWB' Z, UU = 'UPB' Z;
'FOR' J 'FROM' LL+1 'TO' UU 'DO' Z[J]+:= L[J]*Z[J-1] 'OD';
'FOR' J 'FROM' LL 'TO' UU 'DO' Z[J]/:= D[J] 'OD';
'FOR' J 'FROM' UU-1 'BY' -1 'TO' LL 'DO' Z[J]+:= U[J]*Z[J+1] 'OD'
'END'
#$E#;
#E$#
'PROC' ILLU RELAX = ('REF''NETMAT' DEC,'NETMAT'JAC, 'NET' U,F)'VOID':
#$B ILLUR B$#
'BEGIN' 'INT' L1= 1'LWB'U , U1= 1'UPB'U, L2= 2'LWB'U, U2= 2'UPB'U;
( 'NETMAT'(DEC) :=: 'NETMAT'('NIL') ! ILLUDEC (JAC,DEC) );
[L1:U1,L2:U2]'REAL' DU,RH;
RH:= RESIDUAL(JAC,U,F);
SOLL(L1,DEC,RH);
'FOR' I 'FROM' L1+1 'TO' U1
'DO' 'FOR' J 'FROM' L2 'TO' U2
'DO' RH[I,J]-:= JAC[I,J,-3]*RH[I-1,J ] +
( J<U2 ! JAC[I,J,-2]*RH[I-1,J+1] ! 0.0 )
'OD';
SOLL(I,DEC,RH)
'OD';
DU[U1,]:=RH[U1,];
'FOR' I 'FROM' U1-1 'BY' -1 'TO' L1
'DO' 'FOR' J 'FROM' L2 'TO' U2
'DO' DU[I,J] := JAC[I,J, 3]*DU[I+1,J ] +
( J>L2 ! JAC[I,J, 2]*DU[I+1,J-1] ! 0.0 )
'OD';
SOLL(I,DEC,DU);
'FOR' J 'FROM' L2 'TO' U2
'DO' DU[I,J] := RH[I,J] - DU[I,J] 'OD'
'OD';
'FOR' I 'FROM' L1 'TO' U1 'DO'
'FOR' J 'FROM' L2 'TO' U2 'DO'
U[I,J]+:= DU[I,J]
'OD' 'OD'
'END'
#$E#;
#E$#
'PROC' ILLUDEC = ('NETMAT' JAC, 'REF''NETMAT' DECOMP )'VOID':
#$B ILLUD B$#
'BEGIN' 'INT' L1= 1'LWB'JAC, U1= 1'UPB'JAC,
L2= 2'LWB'JAC, U2= 2'UPB'JAC;
'INT' IP;
'REAL' DD,LL,II,L DINV U;
[L2:U2,-1:+1]'REAL' D;
[L2:U2,-2:+2]'REAL' DINV;
[L2:U2,-1:+2]'REAL' L DINV;
'HEAP'[L1:U1,L2:U2,-1:+1]'REAL' DEC;
D[L2:U2,-1:+1]:= JAC[L1,L2:U2,-1:+1];
DD:= DEC[L1,L2,0]:= D[L2,0];
'FOR' J 'FROM' L2 'TO' U2-1
'DO' DEC[L1,J ,+1]:= -D[J ,+1]/DD;
DEC[L1,J+1,-1]:= LL:=-D[J+1,-1]/DD;
DEC[L1,J+1, 0]:= DD:= D[J+1, 0] + D[J,1]*LL
'OD';
'FOR' I 'FROM' L1 'TO' U1-1
'DO' IP:= I+1;
DINV[U2,0]:= II:= 1.0/DEC[I,U2,0];
'FOR' J 'FROM' U2-1 'BY' -1 'TO' L2
'DO' DINV[ J,0]:= II:= 1.0/DEC[I, J,0] +
II * DEC[I,J,1]*DEC[I,J+1,-1]
'OD';
'FOR' K 'TO' 2 'DO'
'FOR' J 'FROM' U2 'BY'-1 'TO' L2+K 'DO'
DINV[J ,-K]:= DINV[J ,1-K]*DEC[I,J-K+1,-1];
DINV[J-K, K]:= DINV[J-K+1,K-1]*DEC[I,J-K ,+1]
'OD' 'OD';
'FOR' K 'FROM' -1 'TO' 2 'DO'
'FOR' J 'FROM' L2+(K=-1!1!0) 'TO' U2-(K=2!2!1)
'DO' L DINV[J ,K]:= JAC[IP,J ,-3]*DINV[J ,K ] +
JAC[IP,J ,-2]*DINV[J+1,K-1]
'OD';
( K<1 !
L DINV[U2,K]:= JAC[IP,U2,-3]*DINV[U2 ,K ] )
'OD';
'FOR' K 'FROM' -1 'TO' 1 'DO'
'FOR' J 'FROM' L2+(K=-1!1!0) 'TO' U2-(K=1!1!0)
'DO' L DINV U := L DINV[J,K ]*JAC[I,J+K ,3];
(J+K<U2 !
L DINV U+:= L DINV[J,K+1]*JAC[I,J+K+1,2] );
D[J,K] := JAC[IP,J,K] - L DINV U
'OD' 'OD';
DD:= DEC[IP,L2,0]:= D[L2,0];
'FOR' J 'FROM' L2 'TO' U2-1
'DO' DEC[IP,J ,+1]:= -D[J ,+1]/DD;
DEC[IP,J+1,-1]:= LL:=-D[J+1,-1]/DD;
DEC[IP,J+1, 0]:= DD:= D[J+1, 0] + D[J,1]*LL
'OD' 'OD';
DECOMP:= DEC
'END'
#$E#;
#E$#
# OTHER POINTWISE RELAXATION # 'PR' EJECT 'PR'
#
'INT' TH RELAX STRATEGY := 1;
#
'PROC' TH RELAX = ('NETMAT'M, 'NET' U,F)'VOID':
#$B RELTH B$#
'BEGIN' 'INT' L1= 1'LWB'U, L2= 2'LWB'U,
U1= 1'UPB'U, U2= 2'UPB'U,
START= 'ABS' TH RELAX STRATEGY;
'FOR' COLOR 'FROM' START
'BY' 'SIGN' TH RELAX STRATEGY
'TO' START + 2
'DO''REF'[ ]'REAL' UIM:= U[L1,@L2], UI, UIP:= U[L1,@L2];
'FOR' I 'FROM' L1 'TO' U1
'DO' ( UI:= UIP; I = U1 ! 'SKIP' ! UIP:= U[I+1,@L2] );
'REF' []'REAL' FI = F[I,@L2];
'REF'[,]'REAL' MI:= M[I,@L2,@-3];
'INT' JM:= L2, JP:= U2;
'FOR' J 'FROM' L2+((COLOR+I)'MOD'3)'BY' 3 'TO' U2
'DO' ( J=L2 ! 'SKIP' ! JM:= J-1 );
( J=U2 ! 'SKIP' ! JP:= J+1 );
'REF'[]'REAL' MIJ = MI[J,@-3];
UI[J] := ( MIJ[-3]*UIM[J] + MIJ[-2]*UIM[JP] +
MIJ[-1]*UI [JM] - FI [J] + MIJ[ 1]*UI [JP] +
MIJ[ 2]*UIP[JM] + MIJ[ 3]*UIP[J] )/ -MIJ[ 0]
'OD';
UIM:= UI
'OD''OD'
'END'
#$E#;
#E$#
'PR' EJECT 'PR'
# 'PROC' FOUR COLOR RELAX = 'REF'[,]'INT':
'HEAP'[1:4,1:2]'INT' := ((0,1),(1,0),(1,1),(0,0));
'REF' [,]'INT' CH RELAX STRATEGY:= FOUR COLOR RELAX;
#
'PROC' CH RELAX = ('NETMAT'M, 'NET' U,F)'VOID':
#$B RELCH B$#
'BEGIN' 'INT' L1= 1'LWB'U, L2= 2'LWB'U,
U1= 1'UPB'U, U2= 2'UPB'U;
[L2:U2]'REAL' DUMMY; 'ZERO' DUMMY;
'FOR' S 'TO' 1'UPB' CH RELAX STRATEGY
'DO' 'INT' C1 = CH RELAX STRATEGY[S,1],
C2 = CH RELAX STRATEGY[S,2];
'REF'[ ]'REAL' UIM, UI, UIP;
'FOR' I 'FROM' L1+C1 'BY' 2 'TO' U1
'DO' UI:= U[I,@L2];
UIM:= ( I=L1 ! DUMMY ! U[I-1,@L2]);
UIP:= ( I=U1 ! DUMMY ! U[I+1,@L2]);
'REF' []'REAL' FI = F[I,@L2];
'REF'[,]'REAL' MI:= M[I,@L2,@-3];
'INT' JM:= L2, JP:= L2;
'FOR' J 'FROM' L2+C2 'BY' 2 'TO' U2
'DO' ( J=U2 ! 'SKIP' ! JP:= J+1 );
'REF'[]'REAL' MIJ = MI[J,@-3];
UI[J] := ( MIJ[-3]*UIM[J] + MIJ[-2]*UIM[JP] +
MIJ[-1]*UI [JM] - FI [J] + MIJ[ 1]*UI [JP] +
MIJ[ 2]*UIP[JM] + MIJ[ 3]*UIP[J] )/ -MIJ[ 0] ;
JM:= JP
'OD';
UIM:= UIP
'OD''OD'
'END'
#$E#;
#E$#
# GRID STYLE PROCEDURES # 'PR' EJECT 'PR'
'PROC' SOLVE DIRECTLY = ('REF''GRID' U, 'GRID' F)'VOID':
#$B SOLVED B$#
SOLVE (JACOBIAN[LEVEL'OF'F],NET'OF'U,NET'OF'F)
#$E#;
#E$#
'OP' 'NORM' = ( 'GRID'A) []'REAL':
#$B NRMG B$#
'TRUE''NORM'A
#$E#;
#E$#
'OP' 'NORM' = ('BOOL' TRUE, 'GRID'A) []'REAL':
#$B NRMGRID B$#
TRUE 'NORM'NET'OF'A
#$E#;
#E$#
'OP''ERRORNORM' = ('GRID' U )[]'REAL':
#$B NRMGERR B$#
'TRUE' 'NORM' (U- SOLUTION)
#$E#;
#E$#
'OP''RESIDUALNORM' = ('GRID' FG, UG )[]'REAL':
#$B NRMRESD B$#
RESIDUALNORM(JACOBIAN[LEVEL'OF'FG],NET'OF'FG,NET'OF'UG)
#$E#;
#E$#
'OP' * = ('NETMAT' M, 'GRID' U)'NET':
#$B NET** B$#
M*NET'OF'U
#$E#;
#E$#
'PROC' RESIDUAL GRID= ('GRID' U,F )'GRID':
#$B RESGR B$#
F - 'LH' U
#$E#;
#E$#
'OP' 'LH' = ('GRID' U)'GRID':
#$B LH B$#
(LEVEL'OF'U, JACOBIAN[LEVEL'OF'U]*(NET'OF'U))
#$E#;
#E$#
# CENTRAL PROCEDURES GRID STYLE # 'PR' EJECT 'PR'
# ************************************************************** #
# ********* GLOBAL PARAMETERS FOR GRID STYLE ******** #
# ************************************************************** #
#
'PROC' PROLONGATE := ('GRID' A)'GRID':
(LEVEL'OF'A+1,LIN INT POL(NET'OF'A)) ;
'PROC' INTERPOLATE:= ('GRID' A)'GRID':
(LEVEL'OF'A+1,SQR INT POL(NET'OF'A)) ;
'PROC' RESTRICTBAR:= ('GRID' A)'GRID':
(LEVEL'OF'A-1,LIN WEIGHT(NET'OF'A)) ;
'PROC' RESTRICT := ('GRID' A)'GRID':
(LEVEL'OF'A-1,INJECTION (NET'OF'A)) ;
'PROC' RELAX := ('REF''GRID' U, 'GRID' F)'VOID':
('INT' L = LEVEL'OF' F; ( L /= LEVEL'OF'U ! ERROR);
ILU! ILURELAX (DECOMP[L], JACOBIAN[L],
NET'OF'U , NET'OF'F )
! LGSRELAX ( JACOBIAN[L],NET'OF'U,NET'OF'F) );
#
'PROC' MG FAS = ('INT' L, 'REF'[]'GRID' U,F)'VOID':
#$B MGFAS B$#
'IF' L = 0
'THEN' SOLVE DIRECTLY (U[0],F[0])
;REP(0,0)
'ELSE' 'TO' P 'DO' RELAX(U[L],F[L]) 'OD';
REP(1,L);
'GRID' Y = RESTRICT ( U[L] );
# ## OR Y = 'COPY' U[L-1]; # ##
F[L-1]:= 'LH' Y + RESTRICTBAR(
RESIDUAL GRID ( U[L], F[L] ))/MU;
REP(2,L-1);
'TO' ( L=1 ! 1 ! S ) 'DO' MG FAS (L-1,U,F) 'OD';
U[L] +:= MU * PROLONGATE (U[L-1]-Y);
REP(3,L);
'TO' Q 'DO' RELAX(U[L],F[L]) 'OD'
;REP(4,L)
'FI'
#$E#;
#E$#
'PROC' MG CS= ('INT' L, 'REF'[]'GRID' U,F)'VOID':
#$B MGCS B$#
'IF' L = 0
'THEN' SOLVE DIRECTLY (U[0],F[0])
;REP(0,0)
'ELSE' 'TO' P 'DO' RELAX(U[L],F[L]) 'OD';
REP(1,L);
'IF' S > 0
'THEN' F[L-1]:= RESTRICTBAR( RESIDUAL GRID( U[L], F[L] ));
'ZERO' U[L-1];
REP(2,L-1);
'TO' (L=1 ! 1 ! S ) 'DO' MG CS(L-1,U,F) 'OD';
U[L] +:= PROLONGATE (U[L-1])
;REP(3,L)
'FI';
'TO' Q 'DO' RELAX(U[L],F[L]) 'OD'
;REP(4,L)
'FI'
#$E#;
#E$#
'PROC' FMGG =('PROBLEM' PROBLEM, 'REF'[]'GRID' U,F)'VOID':
#$B FMGG B$#
'BEGIN' # ## EXPECTS [0:M]'GRID' U; U[0] = STARTAPPROX. # ##
'BOOL' LIN = LINEAR'OF'PROBLEM;
'INT' L = 'UPB'U;
[1:3]'REAL' RES;
JACOBIAN:= 'HEAP'[0:L]'NETMAT';
DECOMP := 'HEAP'[0:L]'NETMAT';
'FOR' I 'FROM' 0 'TO' L
'DO' LEVEL'OF'F[I]:= I;
DECOMP [I]:= 'NETMAT'('NIL')
'OD';
DISCRETIZATION (PROBLEM, U[0],NET'OF'F[0],JACOBIAN[0]);
SOLVE DIRECTLY (U[0],F[0]);
GOON FMGG ( U[0], F[0], RES );
REP (12,0);
'FOR' K 'TO' L
'DO' U[K]:= INTERPOLATE (U[K-1]);
REP (11,K);
DISCRETIZATION (PROBLEM,U[K],NET'OF'F[K],JACOBIAN[K]);
( PR >= 0
! DIRICHLET RELAX(PROBLEM,U[K],JACOBIAN[K],NET'OF'F[K]);
'TO' PR 'DO' RELAX(U[K],F[K]) 'OD'
);
REP (13,K);
'TO' R 'WHILE'
( LIN ! MG CS (K,U,F) ! MG FAS (K,U,F));
GOON FMGG ( U[K], F[K], RES )
'DO' 'SKIP' 'OD'
;REP (12,K)
'OD'
'END'
#$E#;
#E$#
'PROC' GOON FMGG1 = ('GRID' U,F, 'REF'[]'REAL' RES) 'BOOL':
#$B FMGG1 B$#
'BEGIN' 'INT' K = LEVEL'OF'U;
'GRID' RESID = RESIDUAL GRID (U,F);
[]'REAL' RESN = 'FALSE''NORM' RESID;
PRINT((NEWLINE," LEVEL ",WHOLE(K,0)));
FL('ERRORNORM'U);FL(RESN);
( K > 0 ! FX(RES,RESN) );
RES:= RESN; RESN[3] > 1.0E-12
'END'
#$E#;
#E$#
# CENTRAL PROCEDURES NEW STYLE # 'PR' EJECT 'PR'
# ************************************************************** #
# ********* GLOBAL PARAMETERS FOR THE MG-ALGORITHM ******** #
# ************************************************************** #
#
'BOOL' ILU:= 'FALSE';
'INT' P:= 1, S:= 2, Q:= 1, R, PR:= 0; R:= S; 'REAL'MU:= 1.0;
'OP' 'I' = ('NET' N)'NET': SQR INT POL(N);
'OP' 'P' = ('NET' N)'NET': LIN INT POL(N);
'OP' 'R' = ('NET' N)'NET': INJECTION (N);
'OP' 'RBAR'= ('NET' N)'NET': LIN WEIGHT (N);
'PROC' RELAX NET :=
('REF''NETMAT' D,'NET' M, 'NET'U,F )'VOID': LGS RELAX(M,U,F);
#
'PROC' MLA FAS= ('INT'L,'REF'[]'NETMAT'D,J,'REF'[]'NET' U,F)'VOID':
#$B MLAFAS B$#
'IF' L = 0
'THEN' SOLVE (J[0],U[0],F[0]) ;REP(0,0)
'ELSE' 'TO'P 'DO' RELAX NET(D[L],J[L],U[L],F[L])'OD';REP(1,L);
'NET' Y = 'R' U[L]; # ## OR 'COPY' U[L-1] # ##
F[L-1]:= J[L-1]*U[L-1] +
'RBAR' (RESIDUAL (J[L],U[L],F[L])/MU);
REP(2,L-1);
( L=1
! SOLVE (J[0],U[0],F[0])
! 'TO' S 'DO' MLA FAS(L-1,D,J,U,F) 'OD'
);
U[L] +:= MU * 'P' (U[L-1]-Y); REP(3,L);
'TO'Q 'DO' RELAX NET(D[L],J[L],U[L],F[L])'OD';REP(4,L)
'FI'
#$E#;
#E$#
'PROC' MLA CS = ('INT'L,'REF'[]'NETMAT'D,J,'REF'[]'NET' U,F)'VOID':
#$B MLACS B$#
'IF' L = 0
'THEN' SOLVE (J[0],U[0],F[0]) ;REP(0,0)
'ELSE' 'TO'P 'DO' RELAX NET(D[L],J[L],U[L],F[L])'OD';REP(1,L);
F[L-1]:= 'RBAR' RESIDUAL (J[L],U[L],F[L]);
'ZERO' U[L-1]; REP(2,L-1);
( L=1
! SOLVE (J[0],U[0],F[0])
! 'TO' S 'DO' MLA CS (L-1,D,J,U,F) 'OD'
);
U[L] +:= 'P' U[L-1]; REP(3,L);
'TO'Q 'DO' RELAX NET(D[L],J[L],U[L],F[L])'OD';REP(4,L)
'FI'
#$E#;
#E$#
'PR' EJECT 'PR'
'PROC' FMG =('PROBLEM' PROBLEM, 'REF'[]'GRID' UG)'VOID':
#$B FMG B$#
'BEGIN' # ## EXPECTS [0:M]'GRID' UG; UG[0]= STARTAPPROX. # ##
'BOOL' LIN = LINEAR'OF'PROBLEM;
'INT' L = 'UPB'UG;
[1:3]'REAL' RES;
[0:L]'NETMAT' JAC, DEC;
'FOR' I 'FROM' 0 'TO' L
'DO' DECOMP[I]:= 'NETMAT'('NIL') 'OD';
[0:L]'NET' F; 'REF'[]'NET' U= NET'OF'UG;
DISCRETIZATION(PROBLEM, UG[0],F[0],JAC[0]);
SOLVE (JAC[0],U[0],F[0]);
GOON FMG (JAC[0], U[0], F[0], RES );
REP (12,0);
'FOR' K 'TO' L
'DO' UG[K]:= 'GRID'(K,'I' U[K-1]);
REP (11,K);
DISCRETIZATION (PROBLEM,UG[K],F[K],JAC[K]);
( PR >= 0
! DIRICHLET RELAX (PROBLEM,UG[K],JAC[K],F[K]);
'TO'PR'DO' RELAX NET(DEC[K],JAC[K],U[K],F[K])'OD'
);
REP (13,K);
'TO' R 'WHILE'
( LIN ! MLA CS (K,DEC,JAC,U,F)
! MLA FAS(K,DEC,JAC,U,F));
GOON FMG (JAC[K], U[K], F[K], RES )
'DO' 'SKIP' 'OD'
;REP (12,K)
'OD'
'END'
#$E#;
#E$#
# NAG CONTRIBUTION : THE MORE EFFICIENT VERSION # 'PR' EJECT 'PR'
'PROC' ADD INT POL = ('NET' FINE,NET )'VOID':
#$B AIP B$#
'BEGIN' 'INT' L1 = 1'LWB'NET, L2 = 2'LWB'NET,
B1 = 1'UPB'NET, B2 = 2'UPB'NET;
'INT' JJ; 'REAL' U2,U3,U4;
'REF'[]'REAL' UIP= NET[L1,@L2],
UPP= FINE[2*L1,@2*L2];
JJ:= 2*L2; UPP[JJ]+:=(U4:=UIP[L2]);
'FOR' JP 'FROM' L2+1 'TO' B2
'DO' U3:= U4; U4:= UIP[JP];
UPP[JJ+:=1]+:= (U3+U4)/2;
UPP[JJ+:=1]+:= U4
'OD';
'FOR' IP 'FROM' L1+1 'TO' B1
'DO' 'REF'[]'REAL' UI = NET [IP-1 ,@ L2],
UIP = NET [IP ,@ L2],
UMM = FINE[2*IP-1,@2*L2],
UPP = FINE[2*IP ,@2*L2];
JJ:= 2*L2; U2:= UI[L2]; U4:= UIP[L2];
UMM[JJ]:= (U2+U4)/2; UPP[JJ]:= U4;
'FOR' JP 'FROM' L2+1 'TO' B2
'DO' JJ+:= 1; U2:= UI [JP];
U3:= U4; U4:= UIP[JP];
UMM[JJ]+:= (U2+U3)/2;
UPP[JJ]+:= (U3+U4)/2;
JJ+:= 1;
UMM[JJ]+:= (U2+U4)/2;
UPP[JJ]+:= U4
'OD'
'OD'
'END'
#$E#;
#E$#
'PR' EJECT 'PR'
'PROC' WEIGHT RES = ('NETMAT'M,'NET'U,F,COARSE)'VOID':
#$B WR B$#
'BEGIN' 'INT' L1 = (1'LWB'F), L2 = (2'LWB'F),
U1 = (1'UPB'F), U2 = (2'UPB'F);
'INT' NL1 = L1'OVER'2, NL2 = L2'OVER'2,
NU1 = U1'OVER'2, NU2 = U2'OVER'2;
'BOOL' U IS ZERO = 'FALSE';
'INT' II,IM,IP, JJ,JM,JP;
[L2:U2]'REAL' SIM,SII,SIP,ZERO; 'ZERO' ZERO; 'ZERO' SIM;
'FOR' I 'FROM' NL1 'TO' NU1
'DO' II:= 2*I; IP:= II+1; IM:= II-1;
'REF'[]'REAL' CI = COARSE[I,@NL2];
'IF' U IS ZERO
'THEN' SII:= F[II,@L2]; SIP:= F[IP,@L2]
'ELSE'
'FOR' I 'FROM' II 'TO' IP
'WHILE' I <= U1
'DO' 'REF' []'REAL' UIM = (I=L1!ZERO!U[I-1,@L2]),
UI = U[I ,@L2] ,
UIP = (I=U1!ZERO!U[I+1,@L2]),
SI = (I=II!SII!SIP),
FI = F[I,@L2];
'REF'[,]'REAL' MI = M[I,@L2,@-3];
'INT' JM, JJ:= L2, JP:= L2;
'FOR' J 'FROM' L2 'TO' U2
'DO' (JM:= JJ; JJ:= JP; J=U2 ! 'SKIP' ! JP+:= 1 );
'REF'[]'REAL' MIJ = MI[JJ,@-3];
SI[JJ]:=
FI[JJ] - ( MIJ[-3]*UIM[JJ]+MIJ[-2]*UIM[JP]+
MIJ[-1]*UI [JM]+MIJ[ 0]*UI [JJ]+MIJ[ 1]*UI [JP]+
MIJ[ 2]*UIP[JM]+MIJ[ 3]*UIP[JJ] )
'OD'
'OD'
'FI';
JM:= JP:= L2+1;
CI[NL2]:= ( SIM[L2]+SIM[JP]
+SII[JP]
+SIP[L2] )*0.5 + SII[L2];
'FOR' J 'FROM' NL2+1 'TO' NU2-1
'DO' JM:= JP; JJ:= JM+1; JP:= JJ+1;
CI[J]:= ( SIM[JJ]+SIM[JP]
+SII[JM] +SII[JP]
+SIP[JM]+SIP[JJ] )*0.5 + SII[JJ]
'OD';
CI[NU2]:= ( SIM[U2]
+SII[JM]
+SIP[JM]+SIP[U2] )*0.5 + SII[U2];
SIM:= SIP
'OD'
'END'
#$E#;
#E$#
'PR' EJECT 'PR'
'PROC' REDUCE = ('NETMAT' AFI,'REF''NETMAT' ACO,
'NET' FFI,'REF''NET' FCO,UCO)'VOID':
#$B REDUCE B$#
'BEGIN' 'INT' L1 = (1'LWB'AFI)'OVER'2, U1 = (1'UPB'AFI)'OVER'2,
L2 = (2'LWB'AFI)'OVER'2, U2 = (2'UPB'AFI)'OVER'2;
ACO:='HEAP' [L1:U1,L2:U2,-3:3]'REAL';
FCO:='HEAP' [L1:U1,L2:U2 ]'REAL';
UCO:='HEAP' [L1:U1,L2:U2 ]'REAL'; 'ZERO' UCO;
'REAL' Q= 0.25;
'INT' TI,TIP,TK,TKP;
'REAL' FFB,FFD,FFE;
[1:3,1:3,-3:3]'REAL' FINE;
'REF'[]'REAL'
A = FINE[1,1,@-3], B = FINE[1,2,@-3], C = FINE[1,3,@-3],
D = FINE[2,1,@-3], E = FINE[2,2,@-3], F = FINE[2,3,@-3],
G = FINE[3,1,@-3], H = FINE[3,2,@-3], J = FINE[3,3,@-3];
'REF''REAL'
AA =A[ 0], AB =A[ 1], AD =A[ 3],
BA =B[-1], BB =B[ 0], BC =B[ 1], BD =B[ 2], BE =B[ 3],
CB =C[-1], CC =C[ 0], CE =C[ 2], CF =C[ 3],
DA =D[-3], DB =D[-2], DD =D[ 0], DE =D[ 1], DG =D[ 3],
EB =E[-3], EC =E[-2], ED =E[-1], EE =E[ 0],
EF =E[ 1], EG =E[ 2], EH =E[ 3],
FC =F[-3], FE =F[-1], FF =F[ 0], FH =F[ 2], FJ =F[ 3],
GD =G[-3], GE =G[-2], GG =G[ 0], GH =G[ 1],
HE =H[-3], HF =H[-2], HG =H[-1], HH =H[ 0], HJ =H[ 1],
JF =J[-3], JH =J[-1], JJ =J[ 0];
'ZERO' FCO[L1,]; 'ZERO' ACO[ L1, ,];
'FOR' I 'FROM' L1 'TO' U1-1
'DO' TI:= I+I; TIP:= TI+2;
FCO[I+1,L2]:= 0; 'ZERO' ACO[I+1,L2,];
'FOR' K 'FROM' L2 'TO' U2-1
'DO' TK:= K+K; TKP:= TK+2;
FFD:= FFI[TI+1,TK ];
FFB:= FFI[TI ,TK+1];
FFE:= FFI[TI+1,TK+1];
((FCO[I ,K ]+:= FFD+FFB)*:= 0.5)+:= FFI[TI,TK];
FCO[I+1,K ] := FFD+FFE;
FCO[I ,K+1]+:= FFB+FFE;
FINE[1:3,1:3,]:= AFI[TI:TIP,TK:TKP,];
'REF'[]'REAL' A = ACO[I ,K ,@-3],
C = ACO[I ,K+1,@-3],
G = ACO[I+1,K ,@-3],
J = ACO[I+1,K+1,@-3];
# ##AA# ##((A[ 0]+:=
(AB+BA+AD+DA)*2+ BB+DD+BD+DB )*:=Q)+:=AA;
# ##CC# ## C[ 0]+:= (CE+EC+CB+BC)*2+ EE+BB+BE+EB+EF+FE;
# ##GG# ## G[ 0]+:= (GE+EG+GD+DG)*2+ EE+DD+DE+ED+EH+HE;
# ##JJ# ## J[ 0] := FH+HF;
# ##AC# ##( A[ 1]+:= (AB+BC)*2+ BB+BE+DB+DE)*:=Q;
# ##CA# ##( C[-1]+:= (BA+CB)*2+ BB+EB+BD+ED)*:=Q;
# ##AG# ##( A[ 3]+:= (AD+DG)*2 +DD+BD+DE+BE)*:=Q;
# ##GA# ##( G[-3]+:= (DA+GD)*2 +DD+DB+ED+EB)*:=Q;
# ##GC# ##( G[-2] := (GE+EC)*2 + EE+HE+DE+HF+DB+EF+EB)*:=Q;
# ##CG# ##( C[ 2] := (EG+CE)*2 + EE+EH+ED+FH+BD+FE+BE)*:=Q;
# ##GJ# ## G[ 1] := EH+HF+EF;
# ##JG# ## J[-1] := HE+FH+FE;
# ##CJ# ## C[ 3] := EH+EF+FH;
# ##JC# ## J[-3] := HE+FE+HF
'OD';
FFD:= FFI[TI+1,TKP ];
((FCO[I ,U2]+:= FFD)*:=0.5)+:= FFI[TI,2*U2];
FCO[I+1,U2] := FFD;
FINE[1:3,1,]:= AFI[TI:TIP,TKP,];
'REF'[]'REAL' A = ACO[I ,U2,@-3],
G = ACO[I+1,U2,@-3];
# ##AA# ##((A[ 0]+:= (AD+DA)*2 + DD)*:= Q)+:=AA;
# ##GG# ## G[ 0]+:= (GD+DG)*2 + DD;
# ##GA# ##( G[-3]+:= (GD+DA)*2 + DD)*:=Q;
# ##AG# ##( A[ 3]+:= (AD+DG)*2 + DD)*:=Q;
G[-2] := G[ 1]:= 0.0
'OD';
'FOR' K 'FROM' L2 'TO' U2-1
'DO' TK:= K+K; TKP:= TK+2;
FFB:= FFI[2*U1,TK+1];
((FCO[U1,K ]+:= FFB)*:=0.5)+:= FFI[2*U1,TK];
FCO[U1,K+1]+:= FFB;
FINE[1,1:3,]:= AFI[TIP,TK:TKP,];
'REF'[]'REAL' A = ACO[U1,K ,@-3],
C = ACO[U1,K+1,@-3];
# ##AA# ##((A[ 0]+:= (AB+BA)*2 + BB)*:= Q)+:=AA;
# ##CC# ## C[ 0]+:= (CB+BC)*2 + BB;
# ##CA# ##( C[-1]+:= (CB+BA)*2 + BB)*:=Q;
# ##AC# ##( A[ 1]+:= (AB+BC)*2 + BB)*:=Q;
C[ 2] := C[ 3]:= 0.0
'OD';
(FCO[U1,U2] *:=0.5)+:=FFI[2*U1,2*U2];
# ##AA# ##(ACO[U1,U2,0]*:=Q)+:=AFI[2*U1,2*U2,0]
'END'
#$E#;
#E$#
# PRELIMINARY NAG MGM ROUTINE # 'PR' EJECT 'PR'
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #
# EVERYWHERE THE MATRIX DOES NOT DEFINE THE NETMAT M, #
# !!!!!!!!!!! M SHOULD CONTAIN ZEROES !!!!!!!!!!!!!!! #
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #
'PROC' MGM = ('REF'[]'NETMAT' LH, 'REF'[]'NET' UH,FH,
'INT' ITMAX,P,Q,S,T, 'REF'[]'NETMAT' ILLU,
'REF''INT' ITUSED,
'PROC'('INT','NETMAT','NET','NET')'BOOL' GO ON MGM,
'PROC'('INT','STRING')'VOID' FAIL)'VOID':
#$B MGM B$#
'BEGIN' 'INT' L='UPB'UH, R = S;
'REF'[]'NETMAT' ILLUD =
( ILLU :=: 'REF'[]'NETMAT'('NIL')
! 'LOC'[0:L]'NETMAT' ! ILLU );
'PROC' MG = ('INT'L)'VOID':
'IF' L = 0
'THEN' ILLU RELAX(ILLUD[0],LH[0],UH[0],FH[0])
'ELSE' 'TO' P 'DO' ILLU RELAX(ILLUD[L],LH[L],UH[L],FH[L])'OD';
FH[L-1]:= LIN WEIGHT( RESIDUAL (LH[L],UH[L],FH[L]) );
# ## WEIGHT RES (LH[L],UH[L],FH[L],FH[L-1]); # ##
'ZERO' UH[L-1];
'TO' (L=1!T!S) 'DO' MG (L-1) 'OD';
UH[L] +:= LIN INT POL ( UH[L-1]);
# ## ADD INT POL (UH[L],UH[L-1]); # ##
'TO' Q 'DO' ILLU RELAX(ILLUD[L],LH[L],UH[L],FH[L])'OD'
'FI';
'INT' ERR =
( 'LWB' UH /= 0 'OR' 'LWB'FH /= 0 'OR' 'LWB'LH /= 0
'OR' 'UPB'FH /= L 'OR' 'UPB'LH /= L ! 1
!: 'NETMAT' LL = LH[L];
3'LWB'LL /=-3 'OR' 3'UPB'LL /= 3 ! 2
!: S <= 0 'OR' S > 3 'OR' T <= 0 ! 3
!: ITMAX<0 'OR' P<0 'OR' Q<0 ! 4
!: 'LWB'ILLUD /= 0 'OR' 'UPB'ILLUD /=L ! 5
! 0 );
( ERR>0 ! FAIL ( ERR," MGM "));
'IF' ITUSED < 0
'THEN' # ## CREATE GALERKIN APPROXIMATIONS # ##
'FOR' I 'FROM' L 'BY' -1 'TO' 1
'DO' REDUCE(LH[I],LH[I-1],FH[I],FH[I-1],UH[I-1])
'OD';
ITUSED:= 0
'FI';
'IF' ITUSED = 0
'THEN' 'FOR'K 'FROM' 0 'TO' L
'DO' ILLUDEC (LH[K],ILLUD[K]) 'OD';
# ## APPLY FULL MULTIGRID # ##
'TO' T 'DO' MG(0) 'OD';
'FOR' K 'TO' L-1
'DO' UH[K]:= SQR INT POL (UH[K-1]);
'TO' R 'DO' MG (K) 'OD'
'OD'; UH[L]:= SQR INT POL (UH[L-1]);
GOON MGM (ITUSED,LH[L],UH[L],FH[L])
'FI';
'TO' ITMAX
'WHILE' MG (L); ITUSED +:= 1;
GOON MGM (ITUSED, LH[L], UH[L], FH[L])
'DO' 'SKIP' 'OD'
'END'
#$E#;
#E$#
# BLACK-BOX SOLUTION PROCEDURE # 'PR' EJECT 'PR'
'PROC' RAP FMGG =('PROBLEM' PROBLEM, 'REF'[]'GRID' U,F)'VOID':
#$B RAPFMGG B$#
'BEGIN'# ## EXPECTS [0:M]'GRID' U,F; # ##
# ## AND U[0]='GRID'( 0,'REF'[L1:U1,L2:U2]'REAL') # ##
'INT' L = 'UPB'U;
( 'NOT' LINEAR'OF'PROBLEM ! FAIL ( 1, "RAP FMGG") );
( 'LWB'U /= 0 'OR' 'LWB'F /= 0 'OR' 'UPB'F /= L
! FAIL ( 2, "RAP FMGG") );
[0:L]'NETMAT' JAC;
'NETMAT' N0 = NET'OF'(U[0]); 'INT' TTL = 2**L;
DISCRETIZATION (PROBLEM,
'GRID'( L,'LOC'[(1'LWB'N0)*TTL:(1'UPB'N0)*TTL,
(2'LWB'N0)*TTL:(2'UPB'N0)*TTL]'REAL'),
NET'OF'F[L], JACOBIAN[L]);
MGM(JAC,NET'OF'U,NET'OF'F,
100,1,0,1,1,'NIL','LOC''INT',MGM GOON,FAIL)
'END'
#$E#;
#E$#
# OPERATORS WITH ASYMMETRIC WEIGHTS # 'PR' EJECT 'PR'
# FOR EXPERIMENTAL PURPOSES. PWH/830314 #
# [-3:3]'REAL' WASYM,PASYM,RASYM;
PASYM[@1]:= (0.5,0.5,0.5,1.0,0.5,0.5,0.5);
RASYM := WASYM := PASYM;
#
'PROC' WEIGHT = ('NET' NET)'NET':
#$B RESTRW B$#
'BEGIN' 'INT' L1 = (1'LWB'NET), L2 = (2'LWB'NET),
U1 = (1'UPB'NET), U2 = (2'UPB'NET);
'INT' NL1 = L1'OVER'2, NL2 = L2'OVER'2,
NU1 = U1'OVER'2, NU2 = U2'OVER'2;
'REAL'M3=WASYM[ 3],M2=WASYM[ 2],M1=WASYM[ 1],
W3=WASYM[-3],W2=WASYM[-2],W1=WASYM[-1],W0=WASYM[ 0];
'HEAP'[NL1:NU1,NL2:NU2]'REAL' COARSE;
'INT' II,IM,IP, JJ,JM,JP;
[L2:U2]'REAL' ZERO; 'ZERO' ZERO;
'FOR' I 'FROM' NL1 'TO' NU1
'DO' II:= 2*I; IP:= II+1; IM:= II-1;
'REF'[]'REAL' CI = COARSE[I,@NL2],
UIM= ( IM<L1 ! ZERO ! NET[IM,@L2]),
UII= NET[II,@L2] ,
UIP= ( IP>U1 ! ZERO ! NET[IP,@L2]);
JP:= L2+1;
CI[NL2]:= W3*UIM[L2]+W2*UIM[JP]
+W0*UII[L2]+M1*UII[JP]
+M3*UIP[L2] ;
'FOR' J 'FROM' NL2+1 'TO' NU2-1
'DO' JJ:= 2*J; JP:= JJ+1; JM:= JJ-1;
CI[J]:= W3*UIM[JJ]+W2*UIM[JP]
+W1*UII[JM]+W0*UII[JJ]+M1*UII[JP]
+M2*UIP[JM]+M3*UIP[JJ]
'OD';
JM:=U2-1;
CI[NU2]:= W3*UIM[U2]
+W1*UII[JM]+W0*UII[U2]
+M2*UIP[JM]+M3*UIP[U2]
'OD'; COARSE
'END'
#$E#;
#E$#
# GALERKIN OPERATOR CONSTRUCTION # 'PR' EJECT 'PR'
'OP' 'RASP' = ('NETMAT' AFI )'NETMAT':
#$B RAPW B$#
'BEGIN' 'INT' L1 = (1'LWB'AFI)'OVER'2, U1 = (1'UPB'AFI)'OVER'2,
L2 = (2'LWB'AFI)'OVER'2, U2 = (2'UPB'AFI)'OVER'2;
'HEAP' [L1:U1,L2:U2,-3:3]'REAL' ACO; 'ZERO' ACO;
'INT' TI,TK,TIM,TKM;
'REAL'M3=PASYM[-3],M2=PASYM[-2],M1=PASYM[-1],
W3=PASYM[ 3],W2=PASYM[ 2],W1=PASYM[ 1],W0=PASYM[ 0],
Q3=RASYM[-3],Q2=RASYM[-2],Q1=RASYM[-1],
P3=RASYM[ 3],P2=RASYM[ 2],P1=RASYM[ 1],P0=RASYM[ 0];
[1:3,1:3,-3:3]'REAL' FINE;
'REF'[]'REAL'
A = FINE[1,1,@-3], B = FINE[1,2,@-3], C = FINE[1,3,@-3],
D = FINE[2,1,@-3], E = FINE[2,2,@-3], F = FINE[2,3,@-3],
G = FINE[3,1,@-3], H = FINE[3,2,@-3], J = FINE[3,3,@-3];
'REF''REAL'
AA =A[ 0], AB =A[ 1], AD =A[ 3],
BA =B[-1], BB =B[ 0], BC =B[ 1], BD =B[ 2], BE =B[ 3],
CB =C[-1], CC =C[ 0], CE =C[ 2], CF =C[ 3],
DA =D[-3], DB =D[-2], DD =D[ 0], DE =D[ 1], DG =D[ 3],
EB =E[-3], EC =E[-2], ED =E[-1], EE =E[ 0],
EF =E[ 1], EG =E[ 2], EH =E[ 3],
FC =F[-3], FE =F[-1], FF =F[ 0], FH =F[ 2], FJ =F[ 3],
GD =G[-3], GE =G[-2], GG =G[ 0], GH =G[ 1],
HE =H[-3], HF =H[-2], HG =H[-1], HH =H[ 0], HJ =H[ 1],
JF =J[-3], JH =J[-1], JJ =J[ 0];
TI:= 2*L1; TK:= 2*L2;
FINE[3,3,]:= AFI[TI,TK,];
# ##JJ# ##ACO[L1,L2,0]:= P0*JJ*W0;
'FOR' K 'FROM' L2+1 'TO' U2
'DO' TKM:= TK; TK:= 2*K;
FINE[3,1:3,]:= AFI[TI,TKM:TK,];
# ##GG# ##ACO[L1,K-1,0] +:= P0*GH*M1+P1*HG*W0+P1*HH*M1;
# ##JJ# ##ACO[L1,K ,0] +:= P0*JH*W1+Q1*HJ*W0+Q1*HH*W1
+P0*JJ*W0;
# ##GJ# ##ACO[L1,K-1,1] +:= P0*GH*W1+P1*HJ*W0+P1*HH*W1;
# ##JG# ##ACO[L1,K ,-1] +:= Q1*HG*W0+P0*JH*M1+Q1*HH*M1
'OD';
'FOR' I 'FROM' L1+1 'TO' U1
'DO' TIM:= TI; TI:= 2*I; TK:= 2*L2;
FINE[1:3,3,]:= AFI[TIM:TI,TK,];
# ##CC# ##ACO[I-1,L2,0] +:= P0*CF*M3+P3*FC*W0+P3*FF*M3;
# ##JJ# ##ACO[I ,L2,0] +:= P0*JF*W3+Q3*FJ*W0
+Q3*FF*W3+P0*JJ*W0;
# ##CJ# ##ACO[I-1,L2,3] +:= P0*CF*W3+P3*FJ*W0+P3*FF*W3;
# ##JC# ##ACO[I ,L2,-3] +:= Q3*FC*W0+P0*JF*M3+Q3*FF*M3;
'FOR' K 'FROM' L2+1 'TO' U2
'DO' TKM:= TK; TK:= 2*K;
FINE[1:3,1:3,]:= AFI[TIM:TI,TKM:TK,];
'REF'[]'REAL' A = ACO[I-1,K-1,@-3],
C = ACO[I-1,K ,@-3],
G = ACO[I ,K-1,@-3],
J = ACO[I ,K ,@-3];
# ##AA# ## A[ 0]+:= P1*BD*M3+P3*DB*M1;
# ##CC# ## C[ 0]+:= P0*CE*M2+P2*EC*W0+P0*CF*M3
+P3*FC*W0+P2*EE*M2+P3*FF*M3
+Q1*BE*M2+P2*EB*W1+P2*EF*M3+P3*FE*M2;
# ##GG# ## G[ 0]+:= P0*GE*W2+Q2*EG*W0+P0*GH*M1
+P1*HG*W0+Q2*EE*W2+P1*HH*M1
+Q3*DE*W2+Q2*ED*W3+Q2*EH*M1+P1*HE*W2;
# ##JJ# ## J[ 0]+:= P0*JJ*W0+P0*JF*W3+Q3*FJ*W0
+P0*JH*W1+Q1*HJ*W0+Q3*FF*W3
+Q1*HH*W1+Q3*FH*W1+Q1*HF*W3;
# ##AC# ## A[ 1]+:= P1*BE*M2+P3*DB*W1+P3*DE*M2;
# ##CA# ## C[-1]+:= P2*EB*M1+Q1*BD*M3+P2*ED*M3;
# ##AG# ## A[ 3]+:= P1*BD*W3+P3*DE*W2+P1*BE*W2;
# ##GA# ## G[-3]+:= Q3*DB*M1+Q2*ED*M3+Q2*EB*M1;
# ##GC# ## G[-2]+:= P0*GE*M2+Q2*EC*W0+Q2*EE*M2
+P1*HE*M2+Q3*DE*M2+P1*HF*M3
+Q3*DB*W1+Q2*EF*M3+Q2*EB*W1;
# ##CG# ## C[ 2]+:= P2*EG*W0+P0*CE*W2+P2*EE*W2
+P2*EH*M1+P2*ED*W3+P3*FH*M1
+Q1*BD*W3+P3*FE*W2+Q1*BE*W2;
# ##GJ# ## G[ 1]+:= P0*GH*W1+P1*HJ*W0+P1*HH*W1
+Q2*EH*W1+P1*HF*W3+Q2*EF*W3;
# ##JG# ## J[-1]+:= Q1*HG*W0+P0*JH*M1+Q1*HH*M1
+Q1*HE*W2+Q3*FH*M1+Q3*FE*W2;
# ##CJ# ## C[ 3]+:= P0*CF*W3+P3*FJ*W0+P3*FF*W3
+P2*EH*W1+P2*EF*W3+P3*FH*W1;
# ##JC# ## J[-3]+:= Q3*FC*W0+P0*JF*M3+Q3*FF*M3
+Q1*HE*M2+Q3*FE*M2+Q1*HF*M3
'OD' 'OD'; ACO
'END'
#$E#;
#E$#
# ORIENTATION:
ACO = COARSE K-1 K
--------------------------> Y
!
! FINE 1 2 3
! -3 -2
! I-1 1 A -- B -- C ! /
! ! / ! / ! -1 - 0 - 1
! 2 D -- E -- F / !
! ! / ! / ! 2 3
! I 3 G -- H -- J
X V
#
# COLLECTION OF PRIMARY TESTPROBLEMS # 'PR' EJECT 'PR'
# 'PROBLEM' STANDARD = 'PROBLEM' ('TRUE','NIL',
('REF''REAL' AXX,AXY,AYX,AYY, 'REAL' X,Y)'VOID':
(AXX:=A11; AXY:=A12; AYX:=A21; AYY:=A22 ),
('REF''REAL' A,B,C,F, 'REAL' X,Y,U)'VOID':
(A := B1; B := B2; C := C0;
F:= STANDARD RHS(X,Y) ),
('REF''REAL' A, 'REAL' X,Y)'VOID': (A:= 0.0 ),
('REF''REAL' B,C, 'REAL' X,Y,U)'VOID':
(B:= 1.0E+40; C:= B * DIRICHLET(X,Y) ) );
'PROC' STANDARD RHS := ('REAL' X,Y)'REAL': S0 + X*S1 + Y*S2;
'PROC' DIRICHLET := ('REAL' X,Y)'REAL': X*X + Y*Y;
'PROC' SOLUTION := ('REAL' X,Y)'REAL': X*X + Y*Y;
'REAL' A11:=1,A12:=0,A21:=0,A22:=1,B1:= 0,B2:=0,C0:=0,
S0:=-4,S1:=0,S2:=0;
#
'PROC' ANISOTROPIC = ('REAL' EPSILON,ALFA)'VOID':
#$B ANISOP B$#
('REAL' C = COS(ALFA), S = SIN(ALFA);
A11:= EPSILON*C*C + S*S; A12:= A21:= (EPSILON-1)*C*S;
A22:= EPSILON*S*S + C*C;
B1:= B2:= C0:= S1:= S2:= 0; S0:= -2*(EPSILON+1) )
#$E#;
#E$#
'PROC' CONVECTION = ('REAL' EPSILON,ALFA)'VOID':
#$B CONVEC B$#
(A11:= A22:= EPSILON; A12:= A21:= 0;
B1 := COS(ALFA); B2 := SIN(ALFA); C0:= 0;
S0:= -4*EPSILON; S1:= 2*B1; S2:= 2*B2 )
#$E#;
#E$#
'PROC' POISSON = 'VOID':
#$B POISON B$#
ANISOTROPIC(1,0)
#$E#;
#E$#
# OUTPUT PROCEDURES # 'PR' EJECT 'PR'
'OP' 'PRINT' = ('NET' NET)'VOID':
#$B PRNET B$#
('INT' L1 = 1'LWB'NET, L2 = 2'LWB'NET,
B1 = 1'UPB'NET, B2 = 2'UPB'NET;
PRINT((NEWLINE,NEWLINE,L1,B1,L2,B2));
'FOR' I 'FROM' L1 'TO' B1
'DO' PRINT(NEWLINE);
'FOR' J 'FROM' L2 'TO' B2
'DO' PRINT((FLOAT(NET[I,J],11,4,3)," "))
'OD' 'OD'; PRINT(NEWLINE))
#$E#;
#E$#
'OP' 'PRINT' = ('NETMAT' NET)'VOID':
#$B PRNETM B$#
('INT' L1 = 1'LWB'NET, L2 = 2'LWB'NET,
B1 = 1'UPB'NET, B2 = 2'UPB'NET;
PRINT((NEWLINE,NEWLINE, L1,B1,L2,B2));
'FOR' I 'FROM' L1 'TO' B1
'DO' PRINT(NEWLINE);
'FOR' J 'FROM' L2 'TO' B2
'DO' PRINT(NEWLINE);
'FOR' K 'FROM' -3 'TO' +3
'DO' PRINT((FLOAT(NET[I,J,K],12,6,2)," "))
'OD''OD''OD';PRINT(NEWLINE))
#$E#;
#E$#
'OP' 'PRINT' = ('GRID' A)'VOID':
#$B PRGRID B$#
( PRINT((NEWLINE,NEWLINE," LEVEL = ",LEVEL'OF'A));
'PRINT' NET'OF'A )
#$E#;
#E$#
'PROC' FL = ([]'REAL' R)'VOID':
#$B FL B$#
('FOR' I 'TO' 'UPB' R
'DO' PRINT((" ",FLOAT(R[I],12,6,2))) 'OD')
#$E#;
#E$#
'PROC' FX = ([]'REAL' R,S)'VOID':
#$B FX B$#
('FOR' I 'TO' 'UPB' R
'DO' PRINT((" ",FIXED(R[I]/S[I],6,2))) 'OD')
#$E#;
#E$#
'PROC' SPC = 'VOID':
#$B SPC B$#
(COLLECT GARBAGE; PRINT((
NEWLINE,"GARBAGE:",COLLECTIONS,GARBAGE,COLLECT SECONDS,
NEWLINE,"SPACE :",PROGSIZE,STACKSIZE,HEAPSIZE,AVAILABLE,
PROGSIZE+STACKSIZE+HEAPSIZE+AVAILABLE,
NEWLINE,"TIME :",CLOCK )))
#$E#;
#E$#
'OP''ERRORPRINT' = ('GRID' U ) 'VOID':
#$B PRERR B$#
'PRINT'(U-SOLUTION)
#$E#;
#E$#
# GLOBAL PARAMETERS # 'PR' EJECT 'PR'
# ************************************************************** #
# ********* GLOBAL PARAMETERS FOR THE MG-ALGORITHM ******** #
# ************************************************************** #
'PROC' PROLONGATE := ('GRID' A)'GRID': (ERROR; A);
'PROC' INTERPOLATE:= ('GRID' A)'GRID': (ERROR; A);
'PROC' RESTRICTBAR:= ('GRID' A)'GRID': (ERROR; A);
'PROC' RESTRICT := ('GRID' A)'GRID': (ERROR; A);
'PROC' RELAX := ('REF''GRID' U, 'GRID' F)'VOID': (ERROR);
# ************************************************************** #
'BOOL' ILU:= 'FALSE';
'INT' P:= 1, S:= 2, Q:= 1, R, PR:= 0; R:= S; 'REAL'MU:= 1.0;
'OP' 'I' = ('NET' N)'NET': SQR INT POL(N);
'OP' 'P' = ('NET' N)'NET': LIN INT POL(N);
'OP' 'R' = ('NET' N)'NET': INJECTION (N);
'OP' 'RBAR'= ('NET' N)'NET': LIN WEIGHT (N);
'PROC' RELAX NET :=
('REF''NETMAT' D,'NETMAT' M, 'NET'U,F )'VOID': 'SKIP';
# ************************************************************** #
# ********* GLOBAL PARAMETERS FOR THE RELAXATION ******** #
# ************************************************************** #
'BOOL' SYMMETRIC:= 'FALSE', BACKWARD:= 'FALSE',
REVERSE := 'FALSE', ZEBRA := 'FALSE';
'INT' TH RELAX STRATEGY := 1;
[,]'INT' FOUR COLOR RELAX = ((0,1),(1,0),(1,1),(0,0));
'FLEX'[1:1,1:2]'INT' CH RELAX STRATEGY:= FOUR COLOR RELAX;
# ************************************************************** #
# ********* GLOBAL PARAMETERS FOR ASYMM.RESTRICTONS ******* #
# ************************************************************** #
[-3:3]'REAL' WASYM,PASYM,RASYM;
PASYM[@1]:= (0.5,0.5,0.5,1.0,0.5,0.5,0.5);
RASYM := WASYM := PASYM;
# ************************************************************** #
# ********* COMMON GLOBAL PARAMETERS ******************** #
# ************************************************************** #
'REF'[]'NETMAT' JACOBIAN,DECOMP;
'PROC' DISCRETIZATION := ('PROBLEM' PROBLEM, 'GRID' U,
'REF''NET' RHS, 'REF''NETMAT' JAC )'VOID':
'SKIP';
'PROC' BOUNDARY CONDITIONS := ('PROBLEM' PROBLEM, 'GRID' U,
'NET' RHS, 'NETMAT' JAC )'VOID':
'SKIP';
'REAL' Q0:= 0, Q1:= 1, Q2:= 1;
'PROC' LUMP = ('BOOL' B)'VOID':
( B ! Q0:=0; Q1:=1; Q2:=1 ! Q0:=0.25; Q1:=0.50; Q2:=0.25 );
'REAL' HB FACTOR:= 1/3;
'PROC' HB:= ('REF''REAL' AXX,AXY,AYX,AYY,'REAL' BX,BY,H,K)'REAL':0.0;
'PROC' GOON FMGG := ('GRID' U,F, 'REF'[]'REAL' RES)'BOOL': 'TRUE';
'PROC' GOON FMG := ('NETMAT' JAC, 'NET' U,F, 'REF'[]'REAL' RES)
'BOOL': 'TRUE';
'PROC' MGM GOON:=('INT' ITNUM,'NETMAT'LH,'NET'UH,FH)'BOOL':'TRUE';
'PROC' REP := ('INT' NUMBER, LEVEL)'VOID': 'SKIP';
'PROC' FAIL= ('INT' N,'STRING' TEXT)'VOID':
( PRINT((NEWLINE,TEXT,N,NEWLINE)); ERROR);
'BOOL' MONIT:= 'FALSE';
'PR' EJECT 'PR'
# ************************************************************** #
# **************** GLOBAL PARAMETERS ************************* #
# *********** FOR PRIMARY TESTPROBLEMS ********************* #
# ************************************************************** #
'PROBLEM' STANDARD = 'PROBLEM' ('TRUE','NIL',
('REF''REAL' AXX,AXY,AYX,AYY, 'REAL' X,Y)'VOID':
(AXX:=A11; AXY:=A12; AYX:=A21; AYY:=A22 ),
('REF''REAL' A,B,C,F, 'REAL' X,Y,U)'VOID':
(A := B1; B := B2; C := C0;
F:= STANDARD RHS(X,Y) ),
('REF''REAL' A, 'REAL' X,Y)'VOID': (A:= 0.0 ),
('REF''REAL' B,C, 'REAL' X,Y,U)'VOID':
(B:= 1.0E+40; C:= B * DIRICHLET(X,Y) ) );
'PROC' STANDARD RHS := ('REAL' X,Y)'REAL': S0 + X*S1 + Y*S2;
'PROC' DIRICHLET := ('REAL' X,Y)'REAL': X*X + Y*Y;
'PROC' SOLUTION := ('REAL' X,Y)'REAL': X*X + Y*Y;
'REAL' A11:=1,A12:=0,A21:=0,A22:=1,B1:= 0,B2:=0,C0:=0,
S0:=-4,S1:=0,S2:=0;
'PR' PROG 'PR' 'SKIP'
'END'