source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/drv/general_soa_check.f90 @ 3698

Last change on this file since 3698 was 2696, checked in by kanani, 7 years ago

Merge of branch palm4u into trunk

File size: 7.1 KB
Line 
1!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2!   Driver to test the Second Order Adjoint (SOA) model
3!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
4
5PROGRAM 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 
183222   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 
195223   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
226224   FORMAT(10000(E24.14,1X))     
227
228
229END PROGRAM KPP_ROOT_Driver
230
231
Note: See TracBrowser for help on using the repository browser.