[2696] | 1 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 2 | ! Driver to test the Second Order Adjoint (SOA) model |
---|
| 3 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 4 | |
---|
| 5 | PROGRAM KPP_ROOT_Driver |
---|
| 6 | |
---|
| 7 | USE KPP_ROOT_Model |
---|
| 8 | USE KPP_ROOT_Initialize, ONLY: Initialize |
---|
| 9 | |
---|
| 10 | KPP_REAL :: T, DVAL(NSPEC), DV(NVAR) |
---|
| 11 | KPP_REAL :: ADJ1(NVAR), ADJ2(NVAR) |
---|
| 12 | KPP_REAL :: FinalRef(NVAR), FinalRef_tlm(NVAR) |
---|
| 13 | KPP_REAL :: FinalPert(NVAR), FinalPert_tlm(NVAR) |
---|
| 14 | KPP_REAL :: Hs(NVAR,NVAR), AdjTlm(NVAR,NVAR), SOA(NVAR) |
---|
| 15 | INTEGER :: i, j, ind_1 = 1, ind_2 = 2, ind_COST |
---|
| 16 | LOGICAL :: LHessian, LAdjTlm |
---|
| 17 | |
---|
| 18 | !~~~> Number of second order adjoints |
---|
| 19 | ! i.e., number of vectors U_i s.t. Sigma_i = Hess*U_i |
---|
| 20 | ! 1 <= i <= NSOA |
---|
| 21 | |
---|
| 22 | INTEGER, PARAMETER :: NSOA = 1 |
---|
| 23 | |
---|
| 24 | KPP_REAL :: Y_tlm(NVAR,NSOA) |
---|
| 25 | KPP_REAL :: Y_adj(NVAR,1) |
---|
| 26 | KPP_REAL :: Y_soa(NVAR,NSOA) |
---|
| 27 | |
---|
| 28 | STEPMIN = 0.0d0 |
---|
| 29 | STEPMAX = 0.0d0 |
---|
| 30 | |
---|
| 31 | DO i=1,NVAR |
---|
| 32 | RTOL(i) = 1.0d-4 |
---|
| 33 | ATOL(i) = 1.0d-3 |
---|
| 34 | END DO |
---|
| 35 | |
---|
| 36 | |
---|
| 37 | !~~~> The cost function is 0.5*VAR(ind_COST, tF)**2 |
---|
| 38 | ind_COST = ind_O3 |
---|
| 39 | |
---|
| 40 | |
---|
| 41 | !~~~> Compute the full Hessian NVARxNVAR (or not) |
---|
| 42 | LHessian = .TRUE. |
---|
| 43 | |
---|
| 44 | |
---|
| 45 | !~~~> Compute the full Adjoint*Tlm NVARxNVAR (or not) |
---|
| 46 | LAdjTlm = .TRUE. |
---|
| 47 | |
---|
| 48 | !~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 49 | PRINT*,'FIRST RUN, Reference' |
---|
| 50 | CALL Initialize() |
---|
| 51 | DV = 0.0d0 |
---|
| 52 | DV(ind_NO2) = 0.02*VAR(ind_NO2) |
---|
| 53 | DV(ind_NO) = 0.04*VAR(ind_NO) |
---|
| 54 | ! DV(ind_PAN) = -0.2*VAR(ind_PAN) |
---|
| 55 | ! DV(ind_O3) = 0.1*VAR(ind_O3) |
---|
| 56 | ! DO i=1,NVAR |
---|
| 57 | ! DV(i) = rand()*VAR(i) |
---|
| 58 | ! END DO |
---|
| 59 | Y_tlm(1:NVAR,1) = DV |
---|
| 60 | Y_adj(1:NVAR,1) = 0.0d0 |
---|
| 61 | Y_soa(1:NVAR,1:NSOA) = 0.0d0 |
---|
| 62 | CALL INTEGRATE_SOA( NSOA, VAR, Y_tlm, Y_adj, Y_soa, TSTART, TEND, & |
---|
| 63 | ATOL, RTOL) |
---|
| 64 | FinalRef = VAR |
---|
| 65 | FinalRef_tlm(1:NVAR) = Y_tlm(1:NVAR,1) |
---|
| 66 | !~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 67 | |
---|
| 68 | !~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 69 | PRINT*,'FIRST RUN, Perturbed' |
---|
| 70 | CALL Initialize() |
---|
| 71 | VAR = VAR + DV |
---|
| 72 | Y_tlm(1:NVAR,1) = DV |
---|
| 73 | Y_adj(1:NVAR,1) = 0.0d0 |
---|
| 74 | Y_soa(1:NVAR,1:NSOA) = 0.0d0 |
---|
| 75 | CALL INTEGRATE_SOA( NSOA, VAR, Y_tlm, Y_adj, Y_soa, TSTART, TEND, & |
---|
| 76 | ATOL, RTOL) |
---|
| 77 | FinalPert = VAR |
---|
| 78 | FinalPert_tlm(1:NVAR) = Y_tlm(1:NVAR,1) |
---|
| 79 | !~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 80 | |
---|
| 81 | !~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 82 | PRINT*,'PERTURBED RUN' |
---|
| 83 | CALL Initialize() |
---|
| 84 | VAR = VAR + DV |
---|
| 85 | Y_tlm(1:NVAR,1) = DV |
---|
| 86 | Y_adj(1:NVAR,1) = 0.0d0 |
---|
| 87 | Y_adj(ind_COST,1) = FinalPert(ind_COST) |
---|
| 88 | Y_soa(1:NVAR,1) = 0.0d0 |
---|
| 89 | Y_soa(ind_COST,1) = FinalPert_tlm(ind_COST) |
---|
| 90 | CALL INTEGRATE_SOA( NSOA, VAR, Y_tlm, Y_adj, Y_soa, TSTART, TEND, & |
---|
| 91 | ATOL, RTOL) |
---|
| 92 | F2 = VAR(ind_COST)**2/2 |
---|
| 93 | ADJ2 = Y_adj(:,1) |
---|
| 94 | PRINT*,'PERTURBED FUN = ',F2 |
---|
| 95 | !~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 96 | |
---|
| 97 | !~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 98 | PRINT*,'REFERENCE RUN' |
---|
| 99 | CALL Initialize() |
---|
| 100 | Y_tlm(1:NVAR,1) = DV |
---|
| 101 | Y_adj(1:NVAR,1) = 0.0d0 |
---|
| 102 | Y_adj(ind_COST,1) = FinalRef(ind_COST) |
---|
| 103 | Y_soa(1:NVAR,1) = 0.0d0 |
---|
| 104 | Y_soa(ind_COST,1) = FinalRef_tlm(ind_COST) |
---|
| 105 | CALL INTEGRATE_SOA( NSOA, VAR, Y_tlm, Y_adj, Y_soa, TSTART, TEND, & |
---|
| 106 | ATOL, RTOL) |
---|
| 107 | F1 = VAR(ind_COST)**2/2 |
---|
| 108 | ADJ1 = Y_adj(:,1) |
---|
| 109 | SOA = Y_soa(:,1) |
---|
| 110 | PRINT*,'REFERENCE FUN = ',F1 |
---|
| 111 | !~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 112 | |
---|
| 113 | !~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 114 | IF (LHessian) THEN |
---|
| 115 | DO i = 1,NVAR |
---|
| 116 | PRINT*,'.... FULL HESSIAN column ',i,' of ',NVAR |
---|
| 117 | CALL Initialize() |
---|
| 118 | Y_tlm(1:NVAR,1) = 0.0d0 |
---|
| 119 | Y_tlm(i,1) = 1.0d0 |
---|
| 120 | Y_adj(1:NVAR,1) = 0.0d0 |
---|
| 121 | Y_adj(ind_COST,1) = FinalRef(ind_COST) |
---|
| 122 | Y_soa(1:NVAR,1) = 0.0d0 |
---|
| 123 | Y_soa(ind_COST,1) = FinalRef_tlm(ind_COST) |
---|
| 124 | CALL INTEGRATE_SOA( NSOA, VAR, Y_tlm, Y_adj, Y_soa, TSTART, TEND, & |
---|
| 125 | ATOL, RTOL) |
---|
| 126 | CALL Initialize() |
---|
| 127 | Y_soa(1:NVAR,1) = 0.0d0 |
---|
| 128 | Y_soa(ind_COST,1) = Y_tlm(ind_COST,1) |
---|
| 129 | Y_tlm(1:NVAR,1) = 0.0d0 |
---|
| 130 | Y_tlm(i,1) = 1.0d0 |
---|
| 131 | Y_adj(1:NVAR,1) = 0.0d0 |
---|
| 132 | Y_adj(ind_COST,1) = FinalRef(ind_COST) |
---|
| 133 | CALL INTEGRATE_soa(Nsoa, VAR, Y_tlm, Y_adj, Y_soa, TSTART, TEND, & |
---|
| 134 | ATOL, RTOL) |
---|
| 135 | Hs(1:NVAR,i) = Y_soa(1:NVAR,1) |
---|
| 136 | END DO |
---|
| 137 | ! Write Hessian in file |
---|
| 138 | OPEN(150,FILE='KPP_ROOT_Hess.m') |
---|
| 139 | DO i=1,NVAR |
---|
| 140 | WRITE(150,224) (Hs(i,j),j=1,NVAR) |
---|
| 141 | END DO |
---|
| 142 | CLOSE(150) |
---|
| 143 | END IF |
---|
| 144 | !~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 145 | |
---|
| 146 | !~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 147 | IF (LAdjTlm) THEN |
---|
| 148 | DO i = 1,NVAR |
---|
| 149 | PRINT*,'.... FULL ADJ*TLM column ',i,' of ',NVAR |
---|
| 150 | CALL Initialize() |
---|
| 151 | Y_tlm(1:NVAR,1) = 0.0d0 |
---|
| 152 | Y_tlm(i,1) = 1.0d0 |
---|
| 153 | Y_adj(1:NVAR,1) = 0.0d0 |
---|
| 154 | Y_adj(ind_COST,1) = FinalRef(ind_COST) |
---|
| 155 | Y_soa(1:NVAR,1) = 0.0d0 |
---|
| 156 | Y_soa(ind_COST,1) = FinalRef_tlm(ind_COST) |
---|
| 157 | CALL INTEGRATE_SOA( NSOA, VAR, Y_tlm, Y_adj, Y_soa, TSTART, TEND, & |
---|
| 158 | ATOL, RTOL) |
---|
| 159 | CALL Initialize() |
---|
| 160 | Y_adj(1:NVAR,1) = Y_tlm(1:NVAR) |
---|
| 161 | Y_soa(1:NVAR,1) = 0.0d0 |
---|
| 162 | Y_tlm(1:NVAR,1) = 0.0d0 |
---|
| 163 | CALL INTEGRATE_soa(Nsoa, VAR, Y_tlm, Y_adj, Y_soa, TSTART, TEND, & |
---|
| 164 | ATOL, RTOL) |
---|
| 165 | AdjTlm(1:NVAR,i) = Y_adj(1:NVAR,1) |
---|
| 166 | END DO |
---|
| 167 | ! Write Hessian in file |
---|
| 168 | OPEN(150,FILE='KPP_ROOT_AdjTlm.m') |
---|
| 169 | DO i=1,NVAR |
---|
| 170 | WRITE(150,224) (AdjTlm(i,j),j=1,NVAR) |
---|
| 171 | END DO |
---|
| 172 | CLOSE(150) |
---|
| 173 | END IF |
---|
| 174 | !~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 175 | |
---|
| 176 | WRITE(6,*) 'FD versus SOA:' |
---|
| 177 | DO i=1,NVAR |
---|
| 178 | WRITE(6,222) ADJ2(i)-ADJ1(i), SOA(i), & |
---|
| 179 | ABS(ADJ2(i)-ADJ1(i)-SOA(i))/ & |
---|
| 180 | MAX(ABS(ADJ2(i)-ADJ1(i)),ABS(SOA(i)),1.d-16), & |
---|
| 181 | TRIM(SPC_NAMES(i)) |
---|
| 182 | END DO |
---|
| 183 | 222 FORMAT('FD=',E12.5,' soa=',E12.5,' Diff=',E12.5,' (',A,')') |
---|
| 184 | |
---|
| 185 | WRITE(6,FMT="(A14,1X,E14.7)") 'F2-F1 = ',F2-F1 |
---|
| 186 | WRITE(6,FMT="(A14,1X,E14.7)") 'First = ',DOT_PRODUCT(ADJ1,DV) |
---|
| 187 | WRITE(6,FMT="(A14,1X,E14.7)") 'F2-F1-First = ',F2-F1-DOT_PRODUCT(ADJ1,DV) |
---|
| 188 | WRITE(6,FMT="(A14,1X,E14.7)") 'Second = ',DOT_PRODUCT(SOA,DV)/2 |
---|
| 189 | |
---|
| 190 | IF (LHessian) THEN |
---|
| 191 | WRITE(6,*) 'The Hessian:' |
---|
| 192 | DO i=1,NVAR |
---|
| 193 | WRITE(6,223) (Hs(i,j),j=1,NVAR) |
---|
| 194 | END DO |
---|
| 195 | 223 FORMAT(100(E12.5,1X)) |
---|
| 196 | END IF |
---|
| 197 | |
---|
| 198 | |
---|
| 199 | OPEN(150,FILE='KPP_ROOT_Cost.m') |
---|
| 200 | WRITE(150,FMT="(A14,1X,E14.7)") 'F2-F1 = ',F2-F1 |
---|
| 201 | WRITE(150,FMT="(A14,1X,E14.7)") 'First = ',DOT_PRODUCT(ADJ1,DV) |
---|
| 202 | WRITE(150,FMT="(A14,1X,E14.7)") 'F2-F1-First = ',F2-F1-DOT_PRODUCT(ADJ1,DV) |
---|
| 203 | WRITE(150,FMT="(A14,1X,E14.7)") 'Second = ',DOT_PRODUCT(SOA,DV)/2 |
---|
| 204 | CLOSE(150) |
---|
| 205 | |
---|
| 206 | OPEN(150,FILE='KPP_ROOT_FD.m') |
---|
| 207 | DO i=1,NVAR |
---|
| 208 | WRITE(150,224) ADJ2(i)-ADJ1(i), SOA(i), & |
---|
| 209 | ABS(ADJ2(i)-ADJ1(i)-SOA(i))/ & |
---|
| 210 | MAX(ABS(ADJ2(i)-ADJ1(i)),ABS(SOA(i)),1.d-16) |
---|
| 211 | END DO |
---|
| 212 | CLOSE(150) |
---|
| 213 | |
---|
| 214 | OPEN(150,FILE='KPP_ROOT_Sol.m') |
---|
| 215 | DO i=1,NVAR |
---|
| 216 | WRITE(150,224) FinalRef(i) |
---|
| 217 | END DO |
---|
| 218 | CLOSE(150) |
---|
| 219 | |
---|
| 220 | OPEN(150,FILE='KPP_ROOT_Adj.m') |
---|
| 221 | DO i=1,NVAR |
---|
| 222 | WRITE(150,224) ADJ1(i) |
---|
| 223 | END DO |
---|
| 224 | CLOSE(150) |
---|
| 225 | |
---|
| 226 | 224 FORMAT(10000(E24.14,1X)) |
---|
| 227 | |
---|
| 228 | |
---|
| 229 | END PROGRAM KPP_ROOT_Driver |
---|
| 230 | |
---|
| 231 | |
---|