KTH framework for Nek5000 toolboxes; testing version  0.0.1
cvode_driver.f
Go to the documentation of this file.
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
10 c----------------------------------------------------------------------
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
19 c----------------------------------------------------------------------
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
28 c----------------------------------------------------------------------
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
61 c----------------------------------------------------------------------
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
203 c----------------------------------------------------------------------
204  subroutine cdscal_cvode
205 c
206 c Top level driver for CVODE
207 c webpage: https://computation.llnl.gov/casc/sundials
208 c
209 c Integrate the IVP d/dt[y] = f(y(t),t); y(t=t0) := f0
210 c using BDF(stiff) or AM(non-stiff).
211 c f denotes the RHS function and is evaluated in fcvfun().
212 c
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)
266 c 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
330 c----------------------------------------------------------------------
331  subroutine cv_upd_v
332 c
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
358 c----------------------------------------------------------------------
359  subroutine cv_upd_w
360 c
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
385 c----------------------------------------------------------------------
386  subroutine cv_upd_coor
387 c
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
422 c----------------------------------------------------------------------
423  subroutine fcvfun (time_, y, ydot, ipar, rpar, ier)
424 c
425 c Compute RHS function f (called within cvode)
426 c CAUTION: never touch y!
427 c
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
543 c----------------------------------------------------------------------
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 
556 c----------------------------------------------------------------------
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
583 c----------------------------------------------------------------------
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
619 c----------------------------------------------------------------------
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
635 c----------------------------------------------------------------------
void exitt()
Definition: comm_mpi.f:604
subroutine setinvm
Definition: coef.f:1334
subroutine glmapm1
Definition: coef.f:575
subroutine volume
Definition: coef.f:1003
subroutine geom2
Definition: coef.f:823
subroutine geodat1
Definition: coef.f:670
subroutine create_comm(inewcomm)
Definition: comm_mpi.f:441
real *8 function dnekclock()
Definition: comm_mpi.f:393
real *8 function dnekclock_sync()
Definition: comm_mpi.f:401
subroutine set_convect_new(cr, cs, ct, ux, uy, uz)
Definition: convect.f:846
subroutine cvunpack(w1, y) c c copy the internal cvode vector y to nek array w1 c include 'SIZE' include 'TOTAL' include 'CVODE' real w1(lx1
subroutine cv_init
Definition: cvode_driver.f:12
subroutine cv_eval_geom
Definition: cvode_driver.f:545
subroutine fcvfun(time_, y, ydot, ipar, rpar, ier)
Definition: cvode_driver.f:424
subroutine cv_settime
Definition: cvode_driver.f:558
subroutine cv_upd_v
Definition: cvode_driver.f:332
subroutine cv_pstat
Definition: cvode_driver.f:585
subroutine cv_upd_w
Definition: cvode_driver.f:360
subroutine cv_rstat
Definition: cvode_driver.f:621
subroutine cv_upd_coor
Definition: cvode_driver.f:387
subroutine cv_setsize
Definition: cvode_driver.f:3
subroutine cdscal_cvode
Definition: cvode_driver.f:21
subroutine nvec_dssum(u, stride, n, gs_handle)
Definition: dssum.f:261
subroutine dssum(u, nx, ny, nz)
Definition: dssum.f:34
subroutine setdef
Definition: genxyz.f:352
subroutine makeq
Definition: makeq.f:3
subroutine col3(a, b, c, n)
Definition: math.f:598
subroutine sub2(a, b, n)
Definition: math.f:164
subroutine col2(a, b, n)
Definition: math.f:564
real function vlmax(vec, n)
Definition: math.f:396
subroutine copy(a, b, n)
Definition: math.f:260
subroutine add3(a, b, c, n)
Definition: math.f:644
subroutine cfill(a, b, n)
Definition: math.f:244
subroutine invcol3(a, b, c, n)
Definition: math.f:98
subroutine rzero(a, n)
Definition: math.f:208
subroutine setbd(bd, dtbd, nbd)
Definition: navier1.f:1736
subroutine setabbd(ab, dtlag, nab, nbd)
Definition: navier1.f:1663
subroutine sumab(v, vv, vvlag, ntot, ab_, nab_)
Definition: plan4.f:526
subroutine vprops
Definition: vprops.f:2