source: CIVL/examples/fortran/nek5000/core/cvode_driver.f

main
Last change on this file was ea777aa, checked in by Alex Wilton <awilton@…>, 3 years ago

Moved examples, include, build_default.properties, common.xml, and README out from dev.civl.com into the root of the repo.

git-svn-id: ​svn://vsl.cis.udel.edu/civl/trunk@5704 fb995dde-84ed-4084-dfe6-e5aef3e2452c
  • Property mode set to 100755
File size: 17.7 KB
Line 
1#ifndef CVODE
2 subroutine cv_setsize
3
4 include 'SIZE'
5 if(nio.eq.0) write(6,*) 'ABORT: not compiled with CVODE support!'
6 call exitt
7
8 return
9 end
10c----------------------------------------------------------------------
11 subroutine cv_init
12
13 include 'SIZE'
14 if(nio.eq.0) write(6,*) 'ABORT: not compiled with CVODE support!'
15 call exitt
16
17 return
18 end
19c----------------------------------------------------------------------
20 subroutine cdscal_cvode
21
22 include 'SIZE'
23 if(nio.eq.0) write(6,*) 'ABORT: not compiled with CVODE support!'
24 call exitt
25
26 return
27 end
28c----------------------------------------------------------------------
29#else
30 subroutine cv_setsize
31
32 include 'SIZE'
33 include 'TOTAL'
34 include 'CVODE'
35
36 integer*8 i8glsum
37
38 nxyz = lx1*ly1*lz1
39
40 ! set local ODE size
41 cv_nlocal = 0
42 do i = 2,nfield
43 if (ifcvfld(i)) then
44 ntot = nxyz*nelfld(i)
45 cv_nlocal = cv_nlocal + ntot
46 endif
47 enddo
48
49 ! check array size is large enough
50 if (cv_nlocal .gt. cv_lysize) then
51 if(nio.eq.0) write(6,*)
52 & 'ABORT cv_setsize(): workspace too small, check SIZE! '
53 call exitt
54 endif
55
56 ! determine global ODE size
57 cv_nglobal = i8glsum(cv_nlocal,1)
58
59 return
60 end
61c----------------------------------------------------------------------
62 subroutine cv_init
63
64 include 'SIZE'
65 include 'TOTAL'
66 include 'CVODE'
67
68 common /CVWRK1/ y0(cv_lysize)
69
70 ! cvode will not allocate these arrays
71 common /cv_rout/ rout(6),rpar(1)
72
73 integer*8 iout,ipar
74 integer cvcomm
75 common /cv_iout/ iout(21),ipar(1),cvcomm
76
77 real cv_atol_(lx1,ly1,lz1,lelt,ldimt)
78 common /CV_YDOT/ cv_atol_ ! used as scratch
79
80 integer cv_meth
81
82 real*8 etime1
83 character*15 txt_meth,txt_itask
84
85 real atol_t(ldimt)
86
87
88 nxyz = lx1*ly1*lz1
89 ifcvodeinit = .false.
90
91 if(nio.eq.0) write(*,*) 'Initializing CVODE ...'
92 etime1=dnekclock_sync()
93
94 cv_timel = -1
95 call cv_rstat
96
97 ! set solver parameters
98 cv_itask = param(160)
99 cv_meth = param(161) ! AM or BDF
100 cv_rtol = param(163)
101 cv_dtmax = param(164)
102 cv_sigs = param(165)
103 cv_delt = param(166)
104 cv_ipretype = param(167) ! 0: no, 1:left, 2: right
105 cv_maxl = 20 ! max dimension of Krylov subspace
106 cv_iatol = 2 ! 1: scalar 2: vector
107
108 ! setup absolute tolerances
109 if (cv_iatol.eq.1) then
110 cv_atol(1) = param(162)
111 else if (cv_iatol.eq.2) then
112 do i = 2,nfield
113 if (ifcvfld(i)) then
114 ntot = nxyz*nelfld(i)
115 call cfill(cv_atol_(1,1,1,1,i-1),atol(i),ntot)
116 endif
117 enddo
118 call cvpack(cv_atol,cv_atol_,.false.)
119 endif
120
121 ! initialize vector module
122 call create_comm(cvcomm)
123#ifdef MPI
124 call fnvinitp(cvcomm, 1, cv_nlocal, cv_nglobal, ier)
125#else
126 call fnvinits(1, cv_nlocal, ier)
127#endif
128 if (ier.ne.0) then
129 write(*,'(a,i3)') 'ABORT cv_init(): fnvinitp ier=', ier
130 call exitt
131 endif
132
133 ! initialize cvode
134 call cvpack(y0,t,.false.)
135 call fcvmalloc(time, y0, cv_meth, itmeth, cv_iatol,
136 & cv_rtol, cv_atol, iout, rout, ipar, rpar, ier)
137 if (ier.ne.0) then
138 write(*,'(a,i3)') 'ABORT cv_init(): fcvmalloc ier=', ier
139 call exitt
140 endif
141
142 ! initialize linear solver
143 call fcvspGMR(cv_ipretype,igstype,cv_maxl,cv_delt,ier)
144 if (ier .ne. 0) then
145 write(*,'(a,i3)') 'ABORT cv_init(): fcvspgmr ier=', ier
146 call exitt
147 endif
148
149 ! preconiditoner
150 if(cv_ipretype.gt.0) call fcvspilssetprec(1, ier)
151
152 ! enable user-supplied Jacobian routine
153 call fcvspilssetjac(1, ier)
154
155 ! print solver settings
156 etime1 = dnekclock_sync() - etime1
157 if(nio.eq.0) then
158 if(cv_meth.eq.1) txt_meth='AM'
159 if(cv_meth.eq.2) txt_meth='BDF'
160 write(6,'(A,A)') ' integration method : ',
161 & txt_meth
162 if(cv_itask.eq.1) txt_itask='NORMAL'
163 if(cv_itask.eq.3) txt_itask='NORMAL_TSTOP'
164 write(6,'(A,A)') ' integration mode : ',
165 & txt_itask
166 write(6,'(A,i10)') ' Nglobal : ',
167 & cv_nglobal
168 write(6,'(A,1pe8.1)') ' relative tolerance : ',
169 & cv_rtol
170
171 do i = 2,nfield
172 if (ifcvfld(i)) then
173 if (cv_iatol.eq.1) then
174 dd = cv_atol(1)
175 else
176 ntot = nxyz*nelfld(i)
177 dd = vlmax(cv_atol_(1,1,1,1,i-1),ntot)
178 endif
179
180 write(6,1000) i, dd
181 1000 format(3x,'absolute tolerance field ',i3,' : ',1pe8.1)
182 endif
183 enddo
184
185 write(6,'(A,i3)') ' krylov dimension : ',
186 & cv_maxl
187 write(6,'(A,f5.3)') ' ratio linear/non-linear tol : ',
188 & cv_delt
189 write(6,'(A,i3)') ' preconditioner : ',
190 & int(param(167))
191 write(6,'(A,1pe8.1)') ' dt_max : ',
192 & cv_dtmax
193 write(6,'(A,f5.3)') ' increment factor DQJ : ',
194 & cv_sigs
195 write(6,'(A,g15.3,A,/)') ' done :: initializing CVODE',
196 & etime1, ' sec'
197 endif
198
199 ifcvodeinit = .true.
200
201 return
202 end
203c----------------------------------------------------------------------
204 subroutine cdscal_cvode
205c
206c Top level driver for CVODE
207c webpage: https://computation.llnl.gov/casc/sundials
208c
209c Integrate the IVP d/dt[y] = f(y(t),t); y(t=t0) := f0
210c using BDF(stiff) or AM(non-stiff).
211c f denotes the RHS function and is evaluated in fcvfun().
212c
213 include 'SIZE'
214 include 'TOTAL'
215 include 'CVODE'
216
217 common /CVWRK1/ y(cv_lysize)
218
219 common /CVWRK2/ vx_ (lx1,ly1,lz1,lelv)
220 & ,vy_ (lx1,ly1,lz1,lelv)
221 & ,vz_ (lx1,ly1,lz1,lelv)
222 & ,xm1_(lx1,ly1,lz1,lelv)
223 & ,ym1_(lx1,ly1,lz1,lelv)
224 & ,zm1_(lx1,ly1,lz1,lelv)
225 & ,wx_ (lx1,ly1,lz1,lelv)
226 & ,wy_ (lx1,ly1,lz1,lelv)
227 & ,wz_ (lx1,ly1,lz1,lelv)
228
229 integer*8 iout,ipar
230 integer cvcomm
231 common /cv_iout/ iout(21),ipar(1),cvcomm
232
233 nxyz = lx1*ly1*lz1
234 ntot = nxyz * nelv
235
236 if (.not.ifcvodeinit) then
237 write(6,*) 'ABORT: CVODE was not initialized!'
238 call exitt
239 endif
240
241 ! save old output values
242 if (cv_itask.eq.1) then
243 do i=1,21
244 iout_save(i) = iout(i)
245 enddo
246 endif
247
248 ! save coord and mesh vel
249 if(ifmvbd) then
250 call copy(xm1_,xm1,ntot)
251 call copy(ym1_,ym1,ntot)
252 if (if3d) call copy(zm1_,zm1,ntot)
253 call copy(wx_,wx,ntot)
254 call copy(wy_,wy,ntot)
255 if (if3d) call copy(wz_,wz,ntot)
256 endif
257
258 ! save velocities
259 call copy(vx_,vx,ntot)
260 call copy(vy_,vy,ntot)
261 if (if3d) call copy(vz_,vz,ntot)
262
263 ! MAIN solver call
264 time_ = time
265 call fcvsetrin('MAX_STEP',cv_dtmax,ier)
266c call fcvsetiin('MAX_ORD' ,3 ,ier)
267
268 if (cv_itask.eq.3) then
269 if(istep.gt.1) then
270 call cvpack(y,t,.false.)
271 call fcvreinit(timef,y,cv_iatol,cv_rtol,cv_atol,ier)
272 if (ier .ne. 0) then
273 write(*,'(a,i3)') 'ABORT: fcvsetrin ier=', ier
274 call exitt
275 endif
276 cv_timel = 0
277 call cv_rstat
278 endif
279 call fcvsetrin('STOP_TIME',time,ier)
280 endif
281
282 call fcvode(time,tout,y,cv_itask,ier)
283 time = time_
284
285 if (ier.lt.0) then
286 if (nid.eq.0) then
287 write(*,'(a)') ' Restart integrator and try again ...'
288 endif
289 call cvpack(y,t,.false.)
290 call fcvreinit(timef,y,cv_iatol,cv_rtol,cv_atol,ier)
291 cv_timel = 0
292 call cv_rstat
293 call fcvode(time,tout,y,cv_itask,ier)
294 time = time_
295 endif
296
297 if (ier.lt.0) then
298 if (nid.eq.0) then
299 write(*,'(a,i3)') 'ABORT: fcvode returned ier =', ier
300 endif
301 call exitt
302 endif
303
304 cv_istep = cv_istep + 1
305 call cvunpack(t,y)
306
307 ! restore velocities
308 call copy(vx,vx_,ntot)
309 call copy(vy,vy_,ntot)
310 if (if3d) call copy(vz,vz_,ntot)
311 if (param(99).gt.0) call set_convect_new(vxd,vyd,vzd,vx,vy,vz)
312
313 ! restore coord and mesh vel
314 if(ifmvbd) then
315 call copy(xm1,xm1_,ntot)
316 call copy(ym1,ym1_,ntot)
317 if (if3d) call copy(zm1,zm1_,ntot)
318 call copy(wx,wx_,ntot)
319 call copy(wy,wy_,ntot)
320 if (if3d) call copy(wz,wz_,ntot)
321
322 call cv_eval_geom
323 endif
324
325 ! print solver statistics
326 if (nid.eq.0) call cv_pstat
327
328 return
329 end
330c----------------------------------------------------------------------
331 subroutine cv_upd_v
332c
333 include 'SIZE'
334 include 'TSTEP'
335 include 'DEALIAS'
336 include 'SOLN'
337 include 'INPUT'
338 include 'CVODE'
339
340 common /CVWRK2/ vx_ (lx1,ly1,lz1,lelv)
341 & ,vy_ (lx1,ly1,lz1,lelv)
342 & ,vz_ (lx1,ly1,lz1,lelv)
343 & ,xm1_(lx1,ly1,lz1,lelv)
344 & ,ym1_(lx1,ly1,lz1,lelv)
345 & ,zm1_(lx1,ly1,lz1,lelv)
346 & ,wx_ (lx1,ly1,lz1,lelv)
347 & ,wy_ (lx1,ly1,lz1,lelv)
348 & ,wz_ (lx1,ly1,lz1,lelv)
349
350 ntot = lx1*ly1*lz1*nelv
351
352 call sumab(vx,vx_,vxlag,ntot,cv_ab,nbd)
353 call sumab(vy,vy_,vylag,ntot,cv_ab,nbd)
354 if(if3d) call sumab(vz,vz_,vzlag,ntot,cv_ab,nbd)
355
356 return
357 end
358c----------------------------------------------------------------------
359 subroutine cv_upd_w
360c
361 include 'SIZE'
362 include 'TSTEP'
363 include 'MVGEOM'
364 include 'INPUT'
365 include 'CVODE'
366
367 common /CVWRK2/ vx_ (lx1,ly1,lz1,lelv)
368 & ,vy_ (lx1,ly1,lz1,lelv)
369 & ,vz_ (lx1,ly1,lz1,lelv)
370 & ,xm1_(lx1,ly1,lz1,lelv)
371 & ,ym1_(lx1,ly1,lz1,lelv)
372 & ,zm1_(lx1,ly1,lz1,lelv)
373 & ,wx_ (lx1,ly1,lz1,lelv)
374 & ,wy_ (lx1,ly1,lz1,lelv)
375 & ,wz_ (lx1,ly1,lz1,lelv)
376
377 ntot = lx1*ly1*lz1*nelv
378
379 call sumab(wx,wx_,wxlag,ntot,cv_ab,nbd)
380 call sumab(wy,wy_,wylag,ntot,cv_ab,nbd)
381 if(if3d) call sumab(wz,wz_,wzlag,ntot,cv_ab,nbd)
382
383 return
384 end
385c----------------------------------------------------------------------
386 subroutine cv_upd_coor
387c
388 include 'SIZE'
389 include 'TSTEP'
390 include 'GEOM'
391 include 'MVGEOM'
392 include 'INPUT'
393 include 'CVODE'
394
395 common /CVWRK2/ vx_ (lx1,ly1,lz1,lelv)
396 & ,vy_ (lx1,ly1,lz1,lelv)
397 & ,vz_ (lx1,ly1,lz1,lelv)
398 & ,xm1_(lx1,ly1,lz1,lelv)
399 & ,ym1_(lx1,ly1,lz1,lelv)
400 & ,zm1_(lx1,ly1,lz1,lelv)
401 & ,wx_ (lx1,ly1,lz1,lelv)
402 & ,wy_ (lx1,ly1,lz1,lelv)
403 & ,wz_ (lx1,ly1,lz1,lelv)
404
405 COMMON /SCRSF/ dtmp(lx1*ly1*lz1*lelv)
406
407 ntot = lx1*ly1*lz1*nelv
408
409 call sumab(dtmp,wx_,wxlag,ntot,cv_abmsh,nabmsh)
410 call add3 (xm1,xm1_,dtmp,ntot)
411
412 call sumab(dtmp,wy_,wylag,ntot,cv_abmsh,nabmsh)
413 call add3 (ym1,ym1_,dtmp,ntot)
414
415 if(if3d) then
416 call sumab(dtmp,wz_,wzlag,ntot,cv_abmsh,nabmsh)
417 call add3 (zm1,zm1_,dtmp,ntot)
418 endif
419
420 return
421 end
422c----------------------------------------------------------------------
423 subroutine fcvfun (time_, y, ydot, ipar, rpar, ier)
424c
425c Compute RHS function f (called within cvode)
426c CAUTION: never touch y!
427c
428 include 'SIZE'
429 include 'TOTAL'
430 include 'CTIMER'
431 include 'CVODE'
432
433 real time_,y(*),ydot(*),rpar(*)
434 integer*8 ipar(*)
435
436 real w1(lx1,ly1,lz1,lelt),
437 $ w2(lx1,ly1,lz1,lelt),
438 $ w3(lx1,ly1,lz1,lelt)
439
440 real ydott(lx1,ly1,lz1,lelt,ldimt)
441 common /CV_YDOT/ ydott
442
443 ifcvfun = .true.
444 etime1 = dnekclock()
445 time = time_
446 nxyz = lx1*ly1*lz1
447 ntotv = nxyz*nelv
448
449 if (time.ne.cv_timel) then
450 call cv_settime
451
452 if(nio.eq.0) write(6,10) istep,time,time-cv_timel
453 10 format(4x,i7,2x,'t=',1pE14.7,' stepsize=',1pE13.4)
454
455 call cv_upd_v
456 call copy(w1,vx,ntotv)
457 call copy(w2,vy,ntotv)
458 if (if3d) call copy(w3,vz,ntotv)
459
460 if (ifmvbd) then
461 call cv_upd_coor
462 call cv_eval_geom
463 call cv_upd_w
464 call sub2(vx,wx,ntotv)
465 call sub2(vy,wy,ntotv)
466 if (if3d) call sub2(vz,wz,ntotv)
467 endif
468
469 if (param(99).gt.0) call set_convect_new(vxd,vyd,vzd,vx,vy,vz)
470
471 call copy(vx,w1,ntotv)
472 call copy(vy,w2,ntotv)
473 if (if3d) call copy(vz,w3,ntotv)
474
475 cv_timel = time
476 endif
477
478 call cvunpack(t,y)
479
480 if(nid.eq.0 .and. loglevel.gt.2) write(6,*) 'fcvfun',
481 $ ifdqj
482
483 ifield = 1
484 call vprops ! we may use fluid properties somewhere
485 do ifield=2,nfield
486 if (ifcvfld(ifield)) call vprops
487 enddo
488
489 do ifield=2,nfield
490 if (ifcvfld(ifield)) then
491 ntot = nxyz*nelfld(ifield)
492 call makeq
493
494 if (iftmsh(ifield)) then
495 call dssum(bq(1,1,1,1,ifield-1),lx1,ly1,lz1)
496 call col2(bq(1,1,1,1,ifield-1),bintm1,ntot)
497
498 call col3(w1,vtrans(1,1,1,1,ifield),bm1,ntot)
499 call dssum(w1,lx1,ly1,lz1)
500 call col2(w1,bintm1,ntot)
501 else
502 call copy(w1,vtrans(1,1,1,1,ifield),ntot)
503 endif
504
505 call invcol3(ydott(1,1,1,1,ifield-1),bq(1,1,1,1,ifield-1),
506 & w1,ntot)
507 endif
508 enddo
509
510 if (ifgsh_fld_same) then ! all fields are on the v-mesh
511 istride = lx1*ly1*lz1*lelt
512 call nvec_dssum(ydott,istride,nfield-1,gsh_fld(1))
513 else
514 do ifield = 2,nfield
515 if (ifcvfld(ifield) .and. gsh_fld(ifield).ge.0) then
516 if(.not.iftmsh(ifield))
517 & call dssum(ydott(1,1,1,1,ifield-1),lx1,ly1,lz1)
518 endif
519 enddo
520 endif
521
522 do ifield = 2,nfield
523 if (ifcvfld(ifield)) then
524 ntot = nxyz*nelfld(ifield)
525 if (.not.iftmsh(ifield) .and. gsh_fld(ifield).ge.0) then
526 call col2(ydott(1,1,1,1,ifield-1),binvm1,ntot)
527 endif
528 endif
529 enddo
530
531 call cvpack(ydot,ydott,.true.)
532
533 tcvf = tcvf + dnekclock()-etime1
534 ncvf = ncvf + 1
535
536 ier = 0
537 ifcvfun = .false.
538
539 return
540 end
541
542#endif
543c----------------------------------------------------------------------
544 subroutine cv_eval_geom
545
546 call glmapm1
547 call geodat1
548 call geom2
549 call volume
550 call setinvm
551 call setdef
552
553 return
554 end
555
556c----------------------------------------------------------------------
557 subroutine cv_settime
558
559 include 'SIZE'
560 include 'TSTEP'
561 include 'CVODE'
562
563 cv_dtNek = time - timef ! stepsize between nek and cvode
564
565 cv_dtlag(1) = cv_dtNek
566 cv_dtlag(2) = dtlag(2)
567 cv_dtlag(3) = dtlag(3)
568
569 call rzero(cv_abmsh,3)
570 call setabbd(cv_abmsh,cv_dtlag,nabmsh,1)
571 do i = 1,3
572 cv_abmsh(i) = cv_dtNek*cv_abmsh(i)
573 enddo
574
575 call rzero(cv_ab,3)
576 call setabbd(cv_ab,cv_dtlag,nbd,nbd)
577
578 call rzero(cv_bd,4)
579 call setbd(cv_bd,cv_dtlag,nbd)
580
581 return
582 end
583c----------------------------------------------------------------------
584 subroutine cv_pstat
585
586 include 'SIZE'
587 include 'CVODE'
588
589 integer*8 iout,ipar
590 integer cvcomm
591 common /cv_iout/ iout(21),ipar(1),cvcomm
592
593 real nli_nni
594
595 nstp = iout(3) - iout_save(3)
596 nfe = iout(4) - iout_save(4) + iout(20)-iout_save(20)
597 nni = iout(7) - iout_save(7)
598 nli = iout(20)- iout_save(20)
599 nli_nni = float(nli)/max(1,nni)
600
601 nfe_avg = (1.-1./cv_istep)*nfe_avg + nfe/cv_istep
602 nli_nni_avg = (1.-1./cv_istep)*nli_nni_avg + nli_nni/cv_istep
603
604 write(*,'(13x,a8,i8,3x,a8,i8,3x,a10,f5.1)')
605 & 'nsteps: ', nstp ,
606 & 'nfe: ', nfe ,
607 & '<nfe>: ', nfe_avg
608 write(*,'(13x,a8,i8,3x,a8,i8,3x,a10,f5.1)')
609 & 'nni: ', nni ,
610 & 'nli: ', nli ,
611 & '<nli/nni>:', nli_nni_avg
612 write(*,'(13x,a8,i8,3x,a8,i8,3x,a10,i5)')
613 & 'ncfn: ' ,iout(6)-iout_save(6),
614 & 'ncfl: ' ,iout(21)-iout_save(21),
615 & 'netf: ',iout(5)-iout_save(5)
616
617 return
618 end
619c----------------------------------------------------------------------
620 subroutine cv_rstat
621
622 include 'SIZE'
623 include 'CVODE'
624
625 do i=1,21
626 iout_save(i) = 0.0
627 enddo
628
629 nfe_avg = 0
630 nli_nni_avg = 0
631 cv_istep = 0
632
633 return
634 end
635c----------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.