C C AUTHOR: ALFIO BORZI' C INSTITUT FÜR MATHEMATIK C KARL-FRANZENS-UNIVRSITÄT GRAZ C HEINRICHSTR. 36, A-8010 GRAZ C AUSTRIA C E-MAIL: ALFIO.BORZI@UNI-GRAZ.AT C C DATE: 2002 C PROJECT: FULL MULTIGRID FOR OPTIMALITY SYSTEMS C LANGUAGE: FORTRAN 77 C CONSTRAINTS: NONE - PUBLIC DOMAIN SOFTWARE C C FOR THE PRESENTATION OF THIS CODE SEE: C (1) ALFIO BORZI' AND KARL KUNISCH, C The Numerical Solution of the Steady State Solid C Fuel Ignition Model and Its Optimal Control, C SIAM J. Sci. Comp., 22(1) (2000), pp. 263-284. C FOR THE ANALYSIS OF THE CONVERGENCE PROPERTIES OF THIS CODE C AND FOR ACCURACY ESTIMATES (IN THE LINEAR CASE) SEE: C (2) ALFIO BORZI', KARL KUNISCH, AND DO YOUNG KWAK, C Accuracy and Convergence Properties of the Finite C Difference Multigrid Solution of an Optimal Control C Optimality System, SIAM J. Control Optim., (2002). C------------------------------------------------------------------------- C THIS CODE IS DEVELOPED ON TECHNIQUES AND SOFTWARE GIVEN IN: C ACHI BRANDT, C MULTI-LEVEL ADAPTIVE SOLUTIONS TO BOUNDARY-VALUE PROBLEMS, C MATHEMATICS OF COMPUTATION 31 (1977), pp. 333-390. C-------------------------------------------------------------------------- C PROGRAM CONTROLLA C C FULL MULTI-GRID ALGORITHM FOR SOLVING C THE OPTIMALITY SYSTEM FOR THE OPTIMAL CONTROL PROBLEM C C DELTA(Q) + DLT*EXP(Q) = F AND C C min |Q - Z|^2/2 + BL*|EXP(Q)-EXP(Z)|^2/2 + RNU*|F|^2/2 C C WITH DIRICHLET BOUNDARY CONDITIONS ON A RECTANGLE. C RELAXATION IS BY COLLECTIVE GAUSS-SEIDEL-NEWTON. C C THE ALGORITHM CAN WORK WITH THE FIXED OPTION AS WELL AS C WITH THE ACCOMODATIVE OPTION,WITH THE FULL APPROXIMATION SCHEME C C EXPLANATIONS OF PARAMETERS: C ---------------------------- C C NX0- NUMBER OF INTERVALS IN X-DIRECTION C NY0- NUMBER OF INTERVALS IN Y-DIRECTION C H0- LENGTH OF EACH INTERVAL C M- NUMBER OF LEVELS C LIN- THE LEVEL FOR WHICH THE FMG MULTI-GRID PROCESS IS TO BEGIN C C THE PARAMETERS NR1, NR2, ETA AND DELTA ARE USED FOR SWITCHING C BETWEEN WORK LEVEL C C IVW- NUMBER OF 2H-CYCLES WITHIN H-CYCLE. C IVW=1 ENTAILS V-CYCLES C IVW=2 ENTAILS W-CYCLES C C NR1- NUMBER OF RELAXATION SWEEPS IN EACH CYCLE BEFORE TRAN- C SFER TO THE COARSER GRID C C NR2- NUMBER OF SWEEPS IN EACH CYCLE AFTER COMING BACK FROM C THE COARSER GRID C C NR2L- NR2 FOR THE CURRENTLY FINEST LEVEL L C C FOR FULLY ACCOMODATIVE ALGORITHM USE LARGE NR1 AND NR2. C C ETA- ON EACH RELAXATION THE RESIDUALS ARE CALCULATED. C IF THE RATIO BETWEEN THE RESIDUAL IN A CERTAIN C RELAXATION AND THE RESIDUAL IN PREVIOUS RELAXATION C IS GREATER THAN ETA,THRE IS A TRANSFER TO A COARSER C GRID C C DELTA- THE CONVERGENCE CRITERION FOR LEVEL K IS DETERMINED C BY DELTA TIMES THE RESIDUAL NORM OF THE LAST C RELAXATION ON LEVEL K+1 C C ETA=1.E40, DELTA=0 CAN BE USED FOR STRICTLY FIXED ALGORITHM C C FOR THE MULTI-GRID ALGORITHM,THE STOPPING C CRITERIA FOR CURRENTLY FINEST LEVEL L ARE DETERMINED BY THE C PARAMETERS:TOL,RATIO,PREC,PRECM,WMAX,WMAXM,NCYCL,NCYCM C (WHICHEVER IS MET FIRST) C C TOL- TOLERANCE C RATIO- THE TOLERANCE FOR LEVEL L IS TOL*((RATIO)**L) C PREC- THE WORK ON LEVEL L IS STOPPED WHEN THE RESIDUAL C NORM OF RELAXATION ON LEVEL L IS SMALLER THAN PREC C TIMES TAU FOR LEVEL L-1 C PRECM- USED INSTEAD OF PREC WHEN L=M C WMAX- MAXIMUM ALLOWED WORK UNITS FOR WORKING ON LEVEL C L OF THE FULL MULTI-GRID ALGORITHM C WMAXM- USED INSTEAD OF WMAX WHEN WORKING ON LEVEL M C NCYCL- NUMBER OF CYCLES ON CURRENTLY FINEST GRID L C (EXCEPT FOR L=LIN AND L=M) C NCYCIN- NCYC FOR L=LIN C NCYCM- NCYC FOR L=M. C C THE STATE VAR. IS STORED IN Q(.) C THE ADJOINT VAR. IS STORED IN U(.) C C UQ(X,Y), UU(X,Y) - C FUNCTIONS FOR INITIALIZATION OF STATE AND ADJOINT VAR. C FQ(X,Y), FU(X,Y) - C RIGHT HAND SIDES OF THE OPT. SYSTEM C C ZF(X,Y) OBJECTIVE FUNCTION C IMPLICIT REAL *8 (A-H,O-Z) C EXTERNAL UQ,UU,UP,FQ,FU,FP,ZF C C GRID PARAMETERS C NX0=4 NY0=4 H0=0.25 M=7 C C MLAT PARAMETERS C LIN=3 IVW=1 NR1=2 NR2=2 NR2L=2 ETA=1.D40 DELTA=0.0 TOL=1.0D-10 RATIO=.25 PREC=1.0D-10 PRECM=1.0D-10 WMAX=50. WMAXM=50. NCYCL=5 NCYCIN=5 NCYCM=5 C C PROBLEM'S PARAMETER C DLT=6.8 BLT=6.8 RNUT=1.0E-4 C WRITE(*,5) WRITE(*,6) 5 FORMAT(25X,'O U T P U T O F C O N T R O L L A') 6 FORMAT(1H+,25X,25(1H_),//) WRITE(*,2)NX0,NY0,H0,M,LIN,IVW WRITE(*,3)NR1,NR2,NR2L,ETA,DELTA,TOL,RATIO,PREC,PRECM WRITE(*,4)WMAX,WMAXM,NCYCL,NCYCIN,NCYCM 2 FORMAT(1X,'NX0=',I2,1X,',NY0=',I2,1X,',H0=',F5.2,1X,',M=',I2, *1X,',LIN=',I2,1X,',IVW=',I1,',') 3 FORMAT(/,1X,'NR1=',I2,1X,',NR2=',I2,1X,',NR2L=',I2,1X,',ETA=', *F5.2,1X,',DELTA=',F5.2,1X,',TOL=',F5.2,1X,',RATIO=',F4.3,1X, *',PREC=',E11.3,1X,',PRECM=',E11.3,',') 4 FORMAT(/,1X,'WMAX=',F5.1,1X,',WMAXM=',F5.1,1X,',NCYCL=',I2,1X, *',NCYCIN=',I2,1X,',NCYCM=',I2,//) C C MLAT C CALL MULTIG(NX0,NY0,H0,M,LIN,IVW,NR1,NR2,NR2L,ETA,DELTA,RATIO 1,TOL,PREC,PRECM,WMAX,WMAXM,NCYCIN,NCYCL,NCYCM 2,UQ,UU,UP,FQ,FU,FP,ZF,DLT,BLT,RNUT) C CALL PRINTE(M,2*M) STOP END C REAL *8 FUNCTION FQ(X,Y) IMPLICIT REAL *8 (A-H,O-Z) FQ=0. RETURN END C REAL *8 FUNCTION UQ(X,Y) IMPLICIT REAL *8 (A-H,O-Z) UQ=0.0 RETURN END C REAL *8 FUNCTION FU(X,Y) IMPLICIT REAL *8 (A-H,O-Z) FU=0. RETURN END C REAL *8 FUNCTION UU(X,Y) IMPLICIT REAL *8 (A-H,O-Z) UU=0.0 RETURN END C REAL *8 FUNCTION FP(X,Y) IMPLICIT REAL *8 (A-H,O-Z) FP=0. RETURN END C REAL *8 FUNCTION UP(X,Y) IMPLICIT REAL *8 (A-H,O-Z) UP=0. RETURN END C REAL *8 FUNCTION ZF(X,Y) IMPLICIT REAL *8 (A-H,O-Z) A=1.0 PI=4.0D0 * DATAN(1.0D0) PI2=PI*PI C ZF=A*DSIN(PI*X)*DSIN(PI*Y)/PI2 ZF=A/PI2 RETURN END C SUBROUTINE MULTIG(NX0,NY0,H0,M,LIN,IVW,NR1,NR2,NR2L,ETA,DELTA *,RATIO,TOL,PREC,PRECM,WMAX,WMAXM,NCYCIN,NCYCL,NCYCM *,UQ,UU,UP,FQ,FU,FP,ZF,DLT,BLT,RNUT) IMPLICIT REAL *8 (A-H,O-Z) EXTERNAL UQ,UU,UP,FQ,FU,FP,ZF DIMENSION EPS(10),ICYC(10) COMMON/IIQQ/IQ COMMON/PPARA/DL,BL,RNU DL=DLT BL=BLT RNU=RNUT IQ=1 C DO 1 K=1,M K2=2**(K-1) CALL GRDFN(K,NX0*K2+1,NY0*K2+1,H0/K2) 1 CALL GRDFN(K+M,NX0*K2+1,NY0*K2+1,H0/K2) M1=M-1 DO 50 K=LIN,M 50 CALL PUTFQ(K+M,FQ,2) DO 51 K=LIN,M 51 CALL PUTFU(K+M,FU,2) C CALL PUTFQ(LIN,UQ,0) CALL PUTFU(LIN,UU,0) CALL PUTBQ(UQ,LIN) CALL PUTBU(UU,LIN) C DO 55 K=1,M 55 CALL PUTAR(K,ZF) C WU=0. C DO 10 L=LIN,M C C BEGIN NEW FINEST LEVEL L. C WRITE(*,6) L 6 FORMAT(/,20(1H.),I3,2X,20(1H.)/) TOLL=TOL*((RATIO)**L) EPS(L)=TOLL WU=.25*WU ICYC(L)=NCYCL ICYC(M)=NCYCM ICYC(LIN)=NCYCIN IRE=NR1 IBA=1 WMAXL=WMAX IF(L.EQ.M)WMAXL=WMAXM PRECL=PREC IF(L.EQ.M)PRECL=PRECM K=L C C BEGIN A NEW WORK LEVEL K C ERRQO=100.0 ERRUO=100.0 5 ERR=1.E20 3 ERRP=ERR IF(IRE.GT.0) GOTO 9 IF(IBA.EQ.2) GO TO 14 IRE=NR2 IBA=2 GOTO 4 14 ICYC(K)=ICYC(K)-1 IF(K.EQ.L) THEN WRITE(*,*) '---- CONV. FACT.: ', ERRQ/ERRQO, ERRU/ERRUO ERRQO=ERRQ ERRUO=ERRU ENDIF IF(ICYC(K).LE.0) GOTO 2 IRE=NR1 IBA=1 9 CALL RELAX(K,K+M,ERRQ,ERRU) ERR=ERRQ+ERRU IF((ERR/ERRP).GT.1.0) IRE=1 WU=WU+4.**(K-L) IRE=IRE-1 WRITE(*,40)K,ERRQ,ERRU,WU 40 FORMAT(' LEV',I2,' RES.N. Q =', E10.3,1X, * ' RES.N. U =', E10.3, * ' WORK=', F7.3) IF(ERR.LE.EPS(K)) GO TO 2 IF(WU.GE.WMAXL) GOTO 20 IF((ERR/ERRP).LT.ETA) GOTO 3 IF(ERR.LE.EPS(K)) GO TO 2 C C TRANSFER TO COARSER LEVEL C 4 IF(K.EQ.1) GO TO 3 CALL RESCAL(K,K+M,K+M-1) EPS(K-1)=DELTA*ERR K=K-1 CALL PUTU(K+1,K) CALL CRSRES(K,K+M) IF(K.EQ.(L-1)) EPS(L)=MAX(PRECL,TOLL) ICYC(K)=IVW ICYC(1)=1 IRE=NR1 IBA=1 GOTO 5 C C TRANSFER TO FINER GRID C 2 IF(K.EQ.L) GOTO 20 CALL SUBTRT(K+1,K) CALL INTADD(K,K+1) IRE=NR2 IBA=2 K=K+1 IF(K.EQ.L) IRE=NR2L GOTO 5 20 IF(L.EQ.M) GOTO 12 CALL INTRP3(L,L+1) CALL PUTBQ(UQ,L+1) CALL PUTBU(UU,L+1) 12 CALL DIFMX(L,ZF) 10 CONTINUE RETURN END C SUBROUTINE PUTU(KF,KC) IMPLICIT REAL *8 (A-H,O-Z) COMMON IUF(300),IUC(300),IKR0(300) COMMON /SOLUTIO/ Q(200000) COMMON /LAGRMUL/ U(200000) COMMON /TARGET / Z(200000) CALL KEY(KF,IUF,IIF,JJF,HF) CALL KEY(KC,IUC,IIC,JJC,HC) DO 1 IC=1,IIC IF=2*IC-1 IFO=IUF(IF) ICO=IUC(IC) JF=-1 DO 1 JC=1,JJC JF=JF+2 U(ICO+JC)=U(IFO+JF) Q(ICO+JC)=Q(IFO+JF) 1 CONTINUE RETURN END C SUBROUTINE SUBTRT(KF,KC) IMPLICIT REAL *8 (A-H,O-Z) COMMON IUF(300),IUC(300),IKR0(300) COMMON /SOLUTIO/ Q(200000) COMMON /LAGRMUL/ U(200000) COMMON /TARGET / Z(200000) CALL KEY(KF,IUF,IIF,JJF,HF) CALL KEY(KC,IUC,IIC,JJC,HC) DO 1 IC=1,IIC IF=2*IC-1 IFO=IUF(IF) ICO=IUC(IC) JF=-1 DO 1 JC=1,JJC JF=JF+2 U(ICO+JC)=U(ICO+JC)-U(IFO+JF) Q(ICO+JC)=Q(ICO+JC)-Q(IFO+JF) 1 CONTINUE RETURN END C SUBROUTINE CRSRES(K,KRHS) IMPLICIT REAL *8 (A-H,O-Z) COMMON IST(300),IRHS(300),IKR0(300) COMMON /SOLUTIO/ Q(200000) COMMON /LAGRMUL/ U(200000) COMMON /TARGET / Z(200000) COMMON/PPARA/DL,BL,RNU CALL KEY(K,IST,II,JJ,H) CALL KEY(KRHS,IRHS,II,JJ,H) H2=H*H I1=II-1 J1=JJ-1 DO 1 I=2,I1 IR=IRHS(I) IO=IST(I) IM=IST(I-1) IP=IST(I+1) DO 1 J=2,J1 B=-U(IR+J)-U(IO+J+1)-U(IO+J-1)-U(IM+J)-U(IP+J) U(IR+J)=-B-4.*U(IO+J) + DL*H2 * EXP(Q(IO+J)) * U(IO+J) X +H2*(Q(IO+J)-Z(IO+J)) + H2*BL*EXP(Q(IO+J)) X *(EXP(Q(IO+J))-EXP(Z(IO+J))) A=-Q(IR+J)-Q(IO+J+1)-Q(IO+J-1)-Q(IM+J)-Q(IP+J) Q(IR+J)=-A-4.*Q(IO+J) + DL*H2 * EXP(Q(IO+J)) X -H2 * U(IO+J)/RNU 1 CONTINUE RETURN END C SUBROUTINE GRDFN(N,IMAX,JMAX,HH) IMPLICIT REAL *8 (A-H,O-Z) COMMON/GRD/NST(21),IMX(21),JMX(21),H(21) COMMON/IIQQ/IQ NST(N)=IQ IMX(N)=IMAX JMX(N)=JMAX H(N)=HH IQ=IQ+IMAX*JMAX RETURN END C SUBROUTINE KEY(K,IST,IMAX,JMAX,HH) IMPLICIT REAL *8 (A-H,O-Z) COMMON/GRD/NST(21),IMX(21),JMX(21),H(21) DIMENSION IST(1) IMAX=IMX(K) JMAX=JMX(K) IS=NST(K)-JMAX-1 DO 1 I=1,IMAX IS=IS + JMAX 1 IST(I)=IS HH=H(K) RETURN END C SUBROUTINE PUTFQ(K,F,NH) IMPLICIT REAL *8 (A-H,O-Z) EXTERNAL F COMMON IST(300),IKR0(300),IKO(300) COMMON /SOLUTIO/ Q(200000) COMMON /LAGRMUL/ U(200000) COMMON /TARGET / Z(200000) CALL KEY (K,IST,II,JJ,H) H2=H**NH DO 1 I=1,II DO 1 J=1,JJ X=(I-1)*H Y=(J-1)*H 1 Q(IST(I)+J)=F(X,Y)*H2 RETURN END C SUBROUTINE PUTFU(K,F,NH) IMPLICIT REAL *8 (A-H,O-Z) EXTERNAL F COMMON IST(300),IKR0(300),IKO(300) COMMON /SOLUTIO/ Q(200000) COMMON /LAGRMUL/ U(200000) COMMON /TARGET / Z(200000) CALL KEY (K,IST,II,JJ,H) H2=H**NH DO 1 I=1,II DO 1 J=1,JJ X=(I-1)*H Y=(J-1)*H 1 U(IST(I)+J)=F(X,Y)*H2 RETURN END C SUBROUTINE PUTBQ(F,K) IMPLICIT REAL *8 (A-H,O-Z) EXTERNAL F COMMON IST(300),IKR0(300),IKO(300) COMMON /SOLUTIO/ Q(200000) COMMON /LAGRMUL/ U(200000) COMMON /TARGET / Z(200000) CALL KEY(K,IST,II,JJ,H) II1=II-1 DO 1 J=1,JJ X=0. Y=(J-1)*H Q(IST(1)+J)=F(X,Y) X=(II-1)*H Q(IST(II)+J)=F(X,Y) 1 CONTINUE DO 2 I=2,II1 Y=0. X=(I-1)*H Q(IST(I)+1)=F(X,Y) Y=(JJ-1)*H Q(IST(I)+JJ)=F(X,Y) 2 CONTINUE RETURN END C SUBROUTINE PUTBU(F,K) IMPLICIT REAL *8 (A-H,O-Z) EXTERNAL F COMMON IST(300),IKR0(300),IKO(300) COMMON /SOLUTIO/ Q(200000) COMMON /LAGRMUL/ U(200000) COMMON /TARGET / Z(200000) CALL KEY(K,IST,II,JJ,H) II1=II-1 DO 1 J=1,JJ X=0. Y=(J-1)*H U(IST(1)+J)=F(X,Y) X=(II-1)*H U(IST(II)+J)=F(X,Y) 1 CONTINUE DO 2 I=2,II1 Y=0. X=(I-1)*H U(IST(I)+1)=F(X,Y) Y=(JJ-1)*H U(IST(I)+JJ)=F(X,Y) 2 CONTINUE RETURN END C SUBROUTINE PUTAR(K,F) IMPLICIT REAL *8 (A-H,O-Z) EXTERNAL F COMMON IST(300),IKR0(300),IKO(300) COMMON /SOLUTIO/ Q(200000) COMMON /LAGRMUL/ U(200000) COMMON /TARGET / Z(200000) CALL KEY (K,IST,II,JJ,H) DO 1 I=1,II DO 1 J=1,JJ X=(I-1)*H Y=(J-1)*H 1 Z(IST(I)+J)=F(X,Y) RETURN END C SUBROUTINE RELAX(K,KRHS,ERRQ,ERRU) C COLLECTIVE GS NEWTON IMPLICIT REAL *8 (A-H,O-Z) COMMON IST(300),IRHS(300),IKR0(300) COMMON /SOLUTIO/ Q(200000) COMMON /LAGRMUL/ U(200000) COMMON /TARGET / Z(200000) COMMON/PPARA/DL,BL,RNU CALL KEY(K,IST,II,JJ,H) CALL KEY(KRHS,IRHS,II,JJ,H) H2=H*H I1=II-1 J1=JJ-1 ERRQ=0. ERRU=0. C DO 1 I=2,I1 IR=IRHS(I) IO=IST(I) IM=IST(I-1) IP=IST(I+1) DO 1 J=2,J1 RDEN=16.*RNU-8.*DL*EXP(Q(IO+J))*H2*RNU X +H2*H2*(1.+DL*DL*EXP(2.*Q(IO+J))*RNU X +DL*EXP(Q(IO+J))*U(IO+J)+BL*EXP(Q(IO+J)) X *(2.*EXP(Q(IO+J))-EXP(Z(IO+J))) ) QN=(Q(IO+J+1)+Q(IO+J-1)+Q(IM+J)+Q(IP+J)) UN=(U(IO+J+1)+U(IO+J-1)+U(IM+J)+U(IP+J)) TNQ1=(-4.+DL*EXP(Q(IO+J))*H2)*RNU*(-Q(IR+J) X +DL*EXP(Q(IO+J))*H2+QN-4.*Q(IO+J)-H2*U(IO+J)/RNU) TNQ2=H2*(-U(IR+J)+BL*EXP(2.*Q(IO+J))*H2-BL*EXP(Q(IO+J) X +Z(IO+J))*H2+UN+H2*Q(IO+J)-4.*U(IO+J) X +DL*EXP(Q(IO+J))*H2*U(IO+J)-H2*Z(IO+J)) TNU1=(-4.+DL*EXP(Q(IO+J))*H2)*RNU*(-U(IR+J) X +BL*EXP(2.*Q(IO+J))*H2-BL*EXP(Q(IO+J)+Z(IO+J))*H2 X +UN-4.*U(IO+J)+H2*Q(IO+J) X +DL*EXP(Q(IO+J))*U(IO+J)*H2-H2*Z(IO+J)) TNU2=-RNU*H2*(-Q(IR+J)+DL*EXP(Q(IO+J))*H2 X +QN-4.*Q(IO+J)-H2*U(IO+J)/RNU) X *(1.+DL*EXP(Q(IO+J))*U(IO+J)+BL*EXP(Q(IO+J)) X *(2.*EXP(Q(IO+J))-EXP(Z(IO+J))) ) TNQ=(TNQ1+TNQ2)/RDEN TNU=(TNU1+TNU2)/RDEN QX=Q(IO+J)-TNQ UX=U(IO+J)-TNU ER=Q(IO+J)-QX ERRQ=ERRQ+ER*ER Q(IO+J)=QX ER=U(IO+J)-UX ERRU=ERRU+ER*ER U(IO+J)=UX 1 CONTINUE ERRQ=SQRT(ERRQ)/(4.*H) ERRU=SQRT(ERRU)/(4.*H) RETURN END C SUBROUTINE INTADD(KC,KF) IMPLICIT REAL *8 (A-H,O-Z) COMMON ISTC(300),ISTF(300),IKR0(300) COMMON /SOLUTIO/ Q(200000) COMMON /LAGRMUL/ U(200000) COMMON /TARGET / Z(200000) CALL KEY(KC,ISTC,IIC,JJC,HC) CALL KEY(KF,ISTF,IIF,JJF,HF) DO 1 IC=2,IIC IF=2*IC-1 JF=1 IFO=ISTF(IF) IFM=ISTF(IF-1) ICO=ISTC(IC) ICM=ISTC(IC-1) DO 1 JC=2,JJC JF=JF+2 C A=.5*(Q(ICO+JC)+Q(ICO+JC-1)) AM=.5*(Q(ICM+JC)+Q(ICM+JC-1)) Q(IFO+JF) = Q(IFO+JF)+Q(ICO+JC) Q(IFM+JF) = Q(IFM+JF)+.5*(Q(ICO+JC)+Q(ICM+JC)) Q(IFO+JF-1)=Q(IFO+JF-1)+A Q(IFM+JF-1)=Q(IFM+JF-1)+.5*(A+AM) C B=.5*(U(ICO+JC)+U(ICO+JC-1)) BM=.5*(U(ICM+JC)+U(ICM+JC-1)) U(IFO+JF) = U(IFO+JF)+U(ICO+JC) U(IFM+JF) = U(IFM+JF)+.5*(U(ICO+JC)+U(ICM+JC)) U(IFO+JF-1)=U(IFO+JF-1)+B U(IFM+JF-1)=U(IFM+JF-1)+.5*(B+BM) 1 CONTINUE RETURN END C SUBROUTINE RESCAL(KF,KRF,KRC) IMPLICIT REAL *8 (A-H,O-Z) COMMON IUF(300),IRF(300),IRC(300) COMMON /SOLUTIO/ Q(200000) COMMON /LAGRMUL/ U(200000) COMMON /TARGET / Z(200000) COMMON/PPARA/DL,BL,RNU DIMENSION RESQ(257,257), RESU(257,257) CALL KEY(KF,IUF,IIF,JJF,HF) CALL KEY(KRF,IRF,IIF,JJF,HF) CALL KEY(KRC,IRC,IIC,JJC,HC) H2F=HF*HF I1=IIF-1 J1=JJF-1 DO 1 IF=2,I1 IFR=IRF(IF) IFO=IUF(IF) IFM=IUF(IF-1) IFP=IUF(IF+1) DO 1 JF=2,J1 C S=Q(IFO+JF+1)+Q(IFO+JF-1)+Q(IFM+JF)+Q(IFP+JF) RESQ(IF,JF)=4.*(Q(IFR+JF)-S+4.*Q(IFO+JF)) X -4.*H2F*(DL*EXP(Q(IFO+JF)) - U(IFO+JF)/RNU) C R=U(IFO+JF+1)+U(IFO+JF-1)+U(IFM+JF)+U(IFP+JF) RESU(IF,JF)=4.*(U(IFR+JF)-R+4.*U(IFO+JF)) X -4.*H2F* ( DL*EXP(Q(IFO+JF))*U(IFO+JF) X + (Q(IFO+JF)-Z(IFO+JF)) X + BL*EXP(Q(IFO+JF)) X * (EXP(Q(IFO+JF))-EXP(Z(IFO+JF))) ) 1 CONTINUE C TRANSFER IIC1=IIC-1 JJC1=JJC-1 DO 2 IC=2,IIC1 ICR=IRC(IC) IF=2*IC-1 JF=1 DO 2 JC=2,JJC1 JF=JF+2 C Half-weighting SQ=RESQ(IF,JF+1)+RESQ(IF,JF-1)+RESQ(IF-1,JF)+RESQ(IF+1,JF) SU=RESU(IF,JF+1)+RESU(IF,JF-1)+RESU(IF-1,JF)+RESU(IF+1,JF) Q(ICR+JC) = 0.5 * RESQ(IF,JF) + 0.125 * SQ U(ICR+JC) = 0.5 * RESU(IF,JF) + 0.125 * SU 2 CONTINUE RETURN END C SUBROUTINE INTRP3(KC,KF) IMPLICIT REAL *8 (A-H,O-Z) COMMON IUF(300),IUC(300),IKR0(300) COMMON /SOLUTIO/ Q(200000) COMMON /LAGRMUL/ U(200000) COMMON /TARGET / Z(200000) CALL KEY(KC,IUC,IIC,JJC,HC) CALL KEY(KF,IUF,IIF,JJF,HF) DO 20 IC=1,IIC IC0 = IUC(IC) IF0 = IUF(2*IC-1) Q(IF0+1) = Q(IC0+1) Q(IF0+2) = (5*Q(IC0+1)+15*Q(IC0+2)-5*Q(IC0+3)+Q(IC0+4))/16 JJC2 = JJC-2 DO 10 JC=2,JJC2 JF = 2*JC-1 Q(IF0+JF) = Q(IC0+JC) Q(IF0+JF+1) = (-Q(IC0+JC-1)+9.*Q(IC0+JC) + +9.*Q(IC0+JC+1)-Q(IC0+JC+2))/16. 10 CONTINUE Q(IF0+JJF-2) = Q(IC0+JJC-1) Q(IF0+JJF-1) = (Q(IC0+JJC-3)-5.*Q(IC0+JJC-2) + +15.*Q(IC0+JJC-1)+5.*Q(IC0+JJC))/16. Q(IF0+JJF) = Q(IC0+JJC) 20 CONTINUE IM1 = IUF(1) I0 = IUF(2) IP1 = IUF(3) IP3 = IUF(5) IP5 = IUF(7) DO 30 J=1,JJF Q(I0+J) = (5.*Q(IM1+J)+15.*Q(IP1+J) + -5.*Q(IP3+J)+Q(IP5+J))/16. 30 CONTINUE IIF3 = IIF-3 DO 40 I=4,IIF3,2 IM3 = IUF(I-3) IM1 = IUF(I-1) I0 = IUF(I) IP1 = IUF(I+1) IP3 = IUF(I+3) DO 40 J=1,JJF Q(I0+J) = (-Q(IM3+J)+9.*Q(IM1+J) + +9.*Q(IP1+J)-Q(IP3+J))/16. 40 CONTINUE IM5 = IUF(IIF-6) IM3 = IUF(IIF-4) IM1 = IUF(IIF-2) I0 = IUF(IIF-1) IP1 = IUF(IIF) DO 50 J=1,JJF Q(I0+J) = (Q(IM5+J)-5.*Q(IM3+J) + +15.*Q(IM1+J)+5.*Q(IP1+J))/16. 50 CONTINUE C DO 21 IC=1,IIC IC0 = IUC(IC) IF0 = IUF(2*IC-1) U(IF0+1) = U(IC0+1) U(IF0+2) = (5*U(IC0+1)+15*U(IC0+2)-5*U(IC0+3)+U(IC0+4))/16 JJC2 = JJC-2 DO 11 JC=2,JJC2 JF = 2*JC-1 U(IF0+JF) = U(IC0+JC) U(IF0+JF+1) = (-U(IC0+JC-1)+9.*U(IC0+JC) + +9.*U(IC0+JC+1)-U(IC0+JC+2))/16. 11 CONTINUE U(IF0+JJF-2) = U(IC0+JJC-1) U(IF0+JJF-1) = (U(IC0+JJC-3)-5.*U(IC0+JJC-2) + +15.*U(IC0+JJC-1)+5.*U(IC0+JJC))/16. U(IF0+JJF) = U(IC0+JJC) 21 CONTINUE IM1 = IUF(1) I0 = IUF(2) IP1 = IUF(3) IP3 = IUF(5) IP5 = IUF(7) DO 31 J=1,JJF U(I0+J) = (5.*U(IM1+J)+15.*U(IP1+J) + -5.*U(IP3+J)+U(IP5+J))/16. 31 CONTINUE IIF3 = IIF-3 DO 41 I=4,IIF3,2 IM3 = IUF(I-3) IM1 = IUF(I-1) I0 = IUF(I) IP1 = IUF(I+1) IP3 = IUF(I+3) DO 41 J=1,JJF U(I0+J) = (-U(IM3+J)+9.*U(IM1+J) + +9.*U(IP1+J)-U(IP3+J))/16. 41 CONTINUE IM5 = IUF(IIF-6) IM3 = IUF(IIF-4) IM1 = IUF(IIF-2) I0 = IUF(IIF-1) IP1 = IUF(IIF) DO 51 J=1,JJF U(I0+J) = (U(IM5+J)-5.*U(IM3+J) + +15.*U(IM1+J)+5.*U(IP1+J))/16. 51 CONTINUE RETURN END C SUBROUTINE DIFMX(K,F) IMPLICIT REAL *8 (A-H,O-Z) EXTERNAL F COMMON IST(300),IRHS(300),IKO(300) COMMON /SOLUTIO/ Q(200000) COMMON /LAGRMUL/ U(200000) COMMON /TARGET / Z(200000) COMMON/PPARA/DL,BL,RNU CALL KEY(K,IST,II,JJ,H) I1=II-1 J1=JJ-1 DIFTR=0. RFNORM=0.0 DO 1 I=2,I1 IO=IST(I) DO 1 J=2,J1 RFNORM=RFNORM + U(IST(I)+J)*U(IST(I)+J) 1 DIFTR=DIFTR+(Z(IST(I)+J)-Q(IST(I)+J))**2 RFNORM=H*SQRT(RFNORM)/RNU DIFTR =H*SQRT(DIFTR) WRITE(*,70) K, DIFTR,RFNORM 70 FORMAT(5X,' LEVEL ',I2,2X,'Tracking err. |Q-Z|_L2 =',E13.5 * ,3X,'Control |F|_L2 =',E13.5) RETURN END C SUBROUTINE PRINTE(K,KRHS) IMPLICIT REAL *8 (A-H,O-Z) COMMON IST(300),IRHS(300),IKR0(300) COMMON /SOLUTIO/ Q(200000) COMMON /LAGRMUL/ U(200000) COMMON /TARGET / Z(200000) COMMON/PPARA/DL,BL,RNU CALL KEY(K,IST,II,JJ,H) CALL KEY(KRHS,IRHS,II,JJ,H) C C STATIC RESIDUAL C H2=H*H I1=II-1 J1=JJ-1 ERRQ=0. ERRU=0. DO 1 I=2,I1 IR=IRHS(I) IO=IST(I) IM=IST(I-1) IP=IST(I+1) DO 1 J=2,J1 C S=Q(IO+J+1)+Q(IO+J-1)+Q(IM+J)+Q(IP+J) ER=Q(IR+J)-S+4.*Q(IO+J) X -H2*DL*EXP(Q(IO+J)) + H2*U(IO+J)/RNU ERRQ=ERRQ+ER*ER C R=U(IO+J+1)+U(IO+J-1)+U(IM+J)+U(IP+J) ER=U(IR+J)-R+4.*U(IO+J) X -H2*DL*EXP(Q(IO+J))*U(IO+J) X -H2*(Q(IO+J)-Z(IO+J)) X -H2*BL*EXP(Q(IO+J)) X * (EXP(Q(IO+J))-EXP(Z(IO+J))) ERRU=ERRU+ER*ER 1 CONTINUE ERRQ=SQRT(ERRQ)/H ERRU=SQRT(ERRU)/H WRITE(*,*)' STATIC RESIDUAL STATE & ADJOINT:', ERRQ, ERRU C C PRINT THE SOLUTION (FINEST GRID) C JSTEP=1 ISTEP=1 IF(JJ.GT.9) JSTEP=INT(JJ/27.)+1 IF(II.GT.9) ISTEP=INT(II/27.)+1 OPEN(UNIT=10,FILE='solq.dat',STATUS='UNKNOWN') C WRITE(10,*)'STATE' DO 90 I=1,II,ISTEP 90 WRITE(10,111) (Q(IST(I)+J),J=1,JJ,JSTEP) CLOSE(10) OPEN(UNIT=11,FILE='solu.dat',STATUS='UNKNOWN') C WRITE(11,*)'ADJOIT' DO 91 I=1,II,ISTEP 91 WRITE(11,111) (U(IST(I)+J),J=1,JJ,JSTEP) CLOSE(11) 111 FORMAT(1X,128(1X,E8.2)) RETURN END C