# *****************************************************************# # # # # # 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= ( IMU1 ! 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= ( IMU1 ! 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= ( IMU1 ! 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' IL1 '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' IL1 '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' IL1 '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[(IL2!J-1!J), JP= (JL1!I-1!I),], U[I,], U[ (IL2 ! 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 ] + ( JL2 ! 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 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= ( IMU1 ! 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'