source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/drv/general_hsv.f90 @ 4282

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

Merge of branch palm4u into trunk

File size: 9.3 KB
Line 
1!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2!   Driver to test Hessian Singular Vectors
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 :: Hs(NVAR,NVAR), SOA(NVAR)
12      INTEGER :: i, j, ind_1 = 1, ind_2 = 2, ind_COST
13      LOGICAL :: LHessian
14      INTEGER, PARAMETER :: NSOA = 1
15      KPP_REAL :: Y_tlm(NVAR,NSOA)
16      KPP_REAL :: Y_adj(NVAR,1)
17      KPP_REAL :: Y_soa(NVAR,NSOA)
18 
19      integer          lwork, n
20      parameter        (n = NVAR)
21     
22      ! Number of desired eigenvectors
23      INTEGER, PARAMETER :: NEIG = 4
24      double complex   zwork(2000), alpha(NEIG), beta(NEIG), &
25                       eivec(n,NEIG), tmp(n), residu(n)
26
27      integer          kmax, jmax, jmin, maxstep, method, m, l, maxnmv,  &
28                       order, testspace
29      double precision tol, lock, dznrm2
30      logical          wanted
31      double complex   target
32      real             elapse
33      real             etime, tarray(2)
34
35      target = (0.0,0.0)
36      tol = 1.d-9
37      kmax = NEIG           ! Number of wanted solutions
38      jmin = 10; jmax = 20  ! min/max size of search space
39     
40      maxstep = 30  ! Max number of Jacobi-Davidson iterations
41      lock = 1.d-9
42!...     order =  0: nearest to target
43!...     order = -1: smallest real part
44!...     order =  1: largest real part
45!...     order = -2: smallest complex part
46!...     order =  2: largest complex part
47      order = 1
48
49      method = 2  ! 1=gmres(m), 2=cgstab(l)
50      m = 30      ! for gmres(m):
51      l = 2       ! for cgstab(l):
52
53      maxnmv = 20 !...     maximum number of matvecs in cgstab or gmres
54     
55      testspace = 3
56!...     Testspace 1: w = "Standard Petrov" * v (Section 3.1.1)
57!...     Testspace 2: w = "Standard 'variable' Petrov" * v (Section 3.1.2)
58!...     Testspace 3: w = "Harmonic Petrov" * v (Section 3.5.1)
59
60      if ( method .eq. 1 ) then
61         lwork =  4 +  m  + 5*jmax + 3*kmax
62      else
63         lwork = 10 + 6*l + 5*jmax + 3*kmax
64      end if
65      wanted = .true.
66
67      call jdqz(alpha, beta, eivec, wanted, n, target, tol, &
68          kmax, jmax, jmin, &
69          method, m, l, maxnmv, maxstep, &
70          lock, order, testspace, zwork, lwork )
71
72      elapse = etime( tarray )
73
74      print*,'Number of converged eigenvalues = ',kmax
75     
76!...     Compute the norms of the residuals:
77      do j = 1, kmax
78         call amul  ( n, eivec(1,j), residu )
79         call zscal ( n, beta(j), residu, 1)
80         call bmul  ( n, eivec(1,j), tmp )
81         call zaxpy( n, -alpha(j), tmp, 1, residu, 1 )
82         print '("lambda(",i2,"): (",1p,e11.4,",",e11.4, &
83             " )")', j,alpha(j)/beta(j)
84         print '(a30,d13.6)', '||beta Ax - alpha Bx||:', &
85               dznrm2( n, residu, 1 )
86      end do
87      write(*,10) tarray(1), elapse
88
89   10 format(1x,'END JDQZ AFTER ',f6.2,' SEC. CPU-TIME AND ', f6.2, &
90            ' SEC. ELAPSED TIME' )
91
92      open(120,file='KPP_ROOT_vec.m')
93      do i=1,NVAR
94        write(120,120) (eivec(i,j), j=1,NEIG)
95      end do
96      close(120)       
97
98      open(120,file='KPP_ROOT_val.m')
99      do i=1,NEIG
100        write(120,120) ALPHA(i), BETA(i), ALPHA(i)/BETA(i)
101      end do
102      close(120)       
103     
104120   format(10000(E14.7,1X))
105
106END PROGRAM KPP_ROOT_Driver
107
108
109!==========================================================================
110!     Dummy preconditioner subroutine
111!==========================================================================
112      subroutine PRECON( neq, q )
113!...............................................
114!...     Subroutine to compute q = K^-1 q
115!...............................................
116      integer        neq, i
117      double complex q(neq)
118      end subroutine PRECON
119
120
121!==========================================================================
122!     matrix vector multiplication subroutine
123!
124!
125!     Compute the matrix vector multiplication y<---A*x
126!     where A = adjoint(Model)*tlm(Model)
127!
128!     Note: we multiply the real and imaginaray parts separately
129!           x = (xr,xi)
130!           y = (yr,yi)
131!           We run the code twice such that:
132!           yr <-- A*xr and yi <-- A*xi
133!
134!==========================================================================
135      subroutine AMUL (n, u, v)
136  USE KPP_ROOT_Model
137  USE KPP_ROOT_Initialize, ONLY: Initialize
138      integer           n, j
139      Double Complex u(n), v(n)
140      Double precision xr(n), yr(n),  xi(n), yi(n), three, two 
141      parameter         (three = 3.0D+0, two = 2.0D+0)
142      INTEGER, PARAMETER :: NSOA = 1
143      KPP_REAL :: Y_tlm(NVAR,NSOA)
144      KPP_REAL :: Y_adj(NVAR,1)
145      KPP_REAL :: Y_soa(NVAR,NSOA)
146      integer, save :: indexA = 0
147     
148      indexA = indexA + 1
149      print*,'AMUL #',indexA
150
151! Real and imaginary parts of complex inputs
152      xr(1:n) = DBLE( u(1:n) )
153      xi(1:n) = DIMAG( u(1:n) )
154
155
156      DO i=1,NVAR
157        RTOL(i) = 1.0d-4
158        ATOL(i) = 1.0d-3
159      END DO
160!~~~>  The cost function is 0.5*VAR(ind_COST, tF)**2
161      ind_COST = ind_O3
162
163! First do a FWD and a TLM integration
164      CALL Initialize()
165      Y_tlm(1:NVAR,1) = xr(1:NVAR)
166      Y_adj(1:NVAR,1) = 0.0d0
167      Y_soa(1:NVAR,1) = 0.0d0
168      CALL INTEGRATE_SOA( NSOA, VAR, Y_tlm, Y_adj, Y_soa, TSTART, TEND, &
169                            ATOL, RTOL)
170
171! Set Lambda(t_final) = Y_tlm(t_final) and do backward integration
172      CALL Initialize()
173      Y_adj(1:NVAR,1) = Y_tlm(1:NVAR)
174      Y_tlm(1:NVAR,1) = 0.0d0
175      Y_soa(1:NVAR,1) = 0.0d0
176      CALL INTEGRATE_SOA( NSOA, VAR, Y_tlm, Y_adj, Y_soa, TSTART, TEND, &
177                            ATOL, RTOL)
178      yr(1:NVAR) = Y_adj(1:NVAR)
179
180! First do a FWD and a TLM integration
181      CALL Initialize()
182      Y_tlm(1:NVAR,1) = xi(1:NVAR)
183      Y_adj(1:NVAR,1) = 0.0d0
184      Y_soa(1:NVAR,1) = 0.0d0
185      CALL INTEGRATE_SOA( NSOA, VAR, Y_tlm, Y_adj, Y_soa, TSTART, TEND, &
186                            ATOL, RTOL)
187
188! Set Lambda(t_final) = Y_tlm(t_final) and do backward integration
189      CALL Initialize()
190      Y_adj(1:NVAR,1) = Y_tlm(1:NVAR)
191      Y_tlm(1:NVAR,1) = 0.0d0
192      Y_soa(1:NVAR,1) = 0.0d0
193      CALL INTEGRATE_SOA( NSOA, VAR, Y_tlm, Y_adj, Y_soa, TSTART, TEND, &
194                            ATOL, RTOL)
195      yi(1:NVAR) = Y_adj(1:NVAR)
196
197! Convert double outputs to complex
198      DO i=1,n
199        v(i) = DCMPLX( yr(i), yi(i) )
200      END DO   
201
202
203      print*,'   (END AMUL #',indexA,')'
204     
205      end  subroutine AMUL
206
207!==========================================================================
208!     matrix vector multiplication subroutine
209!
210!     Compute the matrix vector multiplication y<---B*x
211!     where B = Hessian, and y = SOA
212!
213!     Note: we multiply the real and imaginaray parts separately
214!           x = (xr,xi)
215!           y = (yr,yi)
216!           We run the code twice such that:
217!           yr <-- B*xr and yi <-- B*xi
218!==========================================================================
219
220      subroutine BMUL (n, u, v)
221  USE KPP_ROOT_Model
222  USE KPP_ROOT_Initialize, ONLY: Initialize
223      integer           n, j
224      Double Complex u(n), v(n)
225      Double precision xr(n), yr(n), xi(n), yi(n)
226      INTEGER, PARAMETER :: NSOA = 1
227      KPP_REAL :: Y_tlm(NVAR,NSOA)
228      KPP_REAL :: Y_adj(NVAR,1)
229      KPP_REAL :: Y_soa(NVAR,NSOA)
230      KPP_REAL :: Final(NVAR), FinalTlm(NVAR)
231      integer, save :: indexB = 0
232     
233      indexB = indexB + 1
234      print*,'BMUL #',indexB
235
236! Real and imaginary parts of complex inputs
237      xr(1:n) = DBLE( u(1:n) )
238      xi(1:n) = DIMAG( u(1:n) )
239
240      DO i=1,NVAR
241        RTOL(i) = 1.0d-4
242        ATOL(i) = 1.0d-3
243      END DO
244!~~~>  The cost function is 0.5*VAR(ind_COST, tF)**2
245      ind_COST = ind_O3
246
247! Forward & TLM integration
248      CALL Initialize()
249      Y_tlm(1:NVAR,1) = xr
250      Y_adj(1:NVAR,1) = 0.0d0
251      Y_soa(1:NVAR,1) = 0.0d0
252      CALL INTEGRATE_SOA( NSOA, VAR, Y_tlm, Y_adj, Y_soa, TSTART, TEND, &
253                            ATOL, RTOL)
254      Final(1:NVAR)    = VAR(1:NVAR)                       
255      FinalTlm(1:NVAR) = Y_tlm(1:NVAR)                     
256                           
257! Adjoint and SOA integration, real part
258      CALL Initialize()
259      Y_adj(1:NVAR,1) = 0.0d0; Y_adj(ind_COST,1) = Final(ind_COST)
260      Y_soa(1:NVAR,1) = 0.0d0; Y_soa(ind_COST,1) = FinalTlm(ind_COST)
261      Y_tlm(1:NVAR,1) = xr
262      CALL INTEGRATE_SOA( NSOA, VAR, Y_tlm, Y_adj, Y_soa, TSTART, TEND, &
263                            ATOL, RTOL)
264                           
265      yr(1:NVAR) = Y_soa(1:NVAR,1)                         
266
267
268! Forward & TLM integration
269      CALL Initialize()
270      Y_tlm(1:NVAR,1) = xi
271      Y_adj(1:NVAR,1) = 0.0d0
272      Y_soa(1:NVAR,1) = 0.0d0
273      CALL INTEGRATE_SOA( NSOA, VAR, Y_tlm, Y_adj, Y_soa, TSTART, TEND, &
274                            ATOL, RTOL)
275      Final(1:NVAR)    = VAR(1:NVAR)                       
276      FinalTlm(1:NVAR) = Y_tlm(1:NVAR)                     
277
278! Adjoint and SOA integration, imaginary part
279      CALL Initialize()
280      Y_adj(1:NVAR,1) = 0.0d0; Y_adj(ind_COST,1) = Final(ind_COST)
281      Y_soa(1:NVAR,1) = 0.0d0; Y_soa(ind_COST,1) = FinalTlm(ind_COST)
282      Y_tlm(1:NVAR,1) = xi
283      CALL INTEGRATE_SOA( NSOA, VAR, Y_tlm, Y_adj, Y_soa, TSTART, TEND, &
284                            ATOL, RTOL)
285                           
286      yi(1:NVAR) = Y_soa(1:NVAR,1)                         
287
288! Convert double outputs to complex
289      DO i=1,n
290        v(i) = DCMPLX( yr(i), yi(i) )
291      END DO   
292
293
294      print*,'   (END BMUL #',indexB,')'
295     
296      end subroutine BMUL 
297! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
298
299
Note: See TracBrowser for help on using the repository browser.