Import latest changes (Nek5000#605)

* direct solve for p0th * fix token parser to get restart file in par * add ifvvisp logical
fdzu · Mar 22, 2019 · c181c74 · c181c74
1 parent c2f51ef
commit c181c74
Show file tree

Hide file tree

Showing 16 changed files with 66 additions and 91 deletions.
diff --git a/bin/visnek b/bin/visnek
@@ -1,17 +1,14 @@
-#!/bin/sh
-#
-#  Shell script that generated VisIt metafile blah.nek5000 either from SESSION.NAME or from the first command line argument
-#
+#!/bin/bash
 
-nfld=0					# # of field files
+nfld=0
 
 if [ $# -eq 0 ]; then
     if [ -f SESSION.NAME ]; then
         base=$(head -1 SESSION.NAME)	# case name
     else
        nfld=-1				# exit
        echo
-       echo '  No SESSION.NAME file -- Use "visnek case_name"'
+       echo '  No SESSION.NAME file found -- use "visnek case_name"'
        echo
     fi
 else
@@ -21,13 +18,6 @@ fi
 
 if [ $nfld -ne -1 ]; then
 
-    echo
-    echo "  Generating $base.nek5000 file ..."
-    echo
-
-    #echo 'endian little'			>  $base.nek5000
-    #echo 'endian big'			>  $base.nek5000
-
     if [ -e "$base.fld01" ]; then
         nfld=$(ls -1 "$base".fld[0-9][0-9]* | wc -l)
         echo "filetemplate: ${base}.fld%02d" > "$base.nek5000"
@@ -55,18 +45,13 @@ if [ $nfld -ne -1 ]; then
     fi
 fi
 
+echo " Found $nfld field file(s)"
+echo " Generating metadata file $base.nek5000 ..."
 
-echo 'firsttimestep: 1' >> "$base.nek5000"
+echo 'firsttimestep: 1'    >> "$base.nek5000"
 echo "numtimesteps: $nfld" >> "$base.nek5000"
+#echo 'meshcoords: 1'      >> "$base.nek5000"
 
-#echo 'meshcoords: 1'			>> "$base.nek5000"
-
-
-echo
-echo " Assuming that coordinates are in the first file ${base}*.f*01 -- otherwise edit $base.nek5000 file ..."
-echo
-echo ' Also one may need to correct for an endian by inserting the first line "endian big" or "endian little" in case of, e.g., files transfered from other systems'
-echo
 echo
-echo " Total $nfld file(s) are in the generated file $base.nek5000 ... Done!"
+echo " Now run visit -o $base.nek5000 or paraview --data=$base.nek5000"
 echo
diff --git a/core/3rd_party/finiparser.c b/core/3rd_party/finiparser.c
@@ -224,7 +224,6 @@ void finiparser_findTokens(char *key, char *delim, int *icounter,int key_len,int
            token[++i] = strtok(NULL,d);
     }
     *icounter = i;
-    free(newstr);
     free(d);  
 
     return;

diff --git a/core/CVODE b/core/CVODE
@@ -2,7 +2,7 @@ c
 c     Variables for the method of characteristics
 c
       integer cv_lysize
-      parameter (cv_lysize=lcvx1*lcvy1*lcvz1*lcvelt*ldimt+1)
+      parameter (cv_lysize=lcvx1*lcvy1*lcvz1*lcvelt*ldimt)
 
       integer         cv_nfld, cv_iatol
       integer         cv_maxl, cv_itask, cv_ipretype
@@ -16,11 +16,10 @@ c
       common /lcvode/ ifcvodeinit, ifdqj, ifcvfun
 
       real cv_atol(cv_lysize)
-
-      real cv_dtlag(lorder),cv_abmsh(lorder),cv_ab(lorder)
+      real cv_dtlag(3),cv_abmsh(3),cv_ab(3),cv_bd(4)
       real cv_rtol, cv_sigs, cv_delt
       real cv_time,cv_timel,cv_dtnek,cv_dtmax
-      common /rcvode/ cv_atol, cv_dtlag, cv_abmsh, cv_ab,
+      common /rcvode/ cv_atol, cv_dtlag, cv_abmsh, cv_ab, cv_bd,
      &                cv_rtol, cv_sigs, cv_delt, 
      &                cv_time, cv_timel, cv_dtnek, cv_dtmax
 

diff --git a/core/INPUT b/core/INPUT
@@ -49,7 +49,7 @@ c
      $               ,ifcyclic,ifmoab,ifcoup, ifvcoup, ifusermv,ifreguo
      $               ,ifxyo_,ifaziv,ifneknek,ifneknekm,ifneknekc
      $               ,ifcvfld(ldimt1),ifdp0dt
-     $               ,ifmpiio,ifrich
+     $               ,ifmpiio,ifrich,ifvvisp
 
       common /input3/ if3d,ifflow,ifheat,iftran,ifaxis,ifstrs,ifsplit
      $               ,ifmgrid 
@@ -65,7 +65,7 @@ c
      $               ,ifcyclic,ifmoab,ifcoup, ifvcoup, ifusermv,ifreguo
      $               ,ifxyo_,ifaziv,ifneknek,ifneknekm,ifneknekc
      $               ,ifcvfld,ifdp0dt
-     $               ,ifmpiio,ifrich
+     $               ,ifmpiio,ifrich,ifvvisp
 
       logical         ifnav
       equivalence    (ifnav, ifadvc(1))

diff --git a/core/PARALLEL.patch b/core/PARALLEL.patch
@@ -16,6 +16,3 @@ index 212b297..1b6d207 100644
       $       ,nelgv,nelgt
 -      common /hcglb/ nvtot,nelg,lglel,gllel,gllnid,nelgv,nelgt
 +      common /hcglb/ nvtot,nelg,lglel,nelgv,nelgt
-
-       logical         ifgprnt
-       common /diagl/  ifgprnt
diff --git a/core/SOLN b/core/SOLN
@@ -88,8 +88,8 @@ c     Solution data for magnetic field
       real qtl(lx2,ly2,lz2,lelt), usrdiv(lx2,ly2,lz2,lelt)
       common /diverg/ qtl, usrdiv
 
-      real p0th, dp0thdt, gamma0
-      common /p0therm/ p0th, dp0thdt, gamma0
+      real p0th, dp0thdt, gamma0, p0thn, p0thlag(2)
+      common /p0therm/ p0th, dp0thdt, gamma0, p0thn, p0thlag
 
       real  v1mask (lx1,ly1,lz1,lelv)
      $     ,v2mask (lx1,ly1,lz1,lelv)

diff --git a/core/cvode_aux.h b/core/cvode_aux.h
@@ -1,13 +1,12 @@
-      subroutine cvunpack(w1,w1p,y)
+      subroutine cvunpack(w1,y)
 c
-c     copy the internal cvode vector y to nek array w1 and w1p 
+c     copy the internal cvode vector y to nek array w1 
 c
       include 'SIZE'
       include 'TOTAL'
       include 'CVODE'
 
       real w1(lx1,ly1,lz1,lelt,*)
-      real w1p
       real y(*)
 
       nxyz = lx1*ly1*lz1
@@ -22,17 +21,12 @@ c
            endif
         enddo
 
-      if (ifdp0dt) then
-         w1p = y(j)
-         j = j + 1 
-      endif
-
       return
       end
 c----------------------------------------------------------------------
-      subroutine cvpack(y,w1,w1p,ifrhs)
+      subroutine cvpack(y,w1,ifrhs)
 c
-c     copy the nek array w1 and w1p to the internal cvode vector y
+c     copy the nek array w1 to the internal cvode vector y
 c     note: assumes temperature is stored in ifield=2 (only for ifdp0dt)
 c
       include 'SIZE'
@@ -41,7 +35,6 @@ c
 
       real y(*)
       real w1(lx1,ly1,lz1,lelt,*)
-      real w1p
       logical ifrhs
 
       common /scrsf/ dtmp(lx1,ly1,lz1,lelt)
@@ -69,11 +62,5 @@ c
          endif
       enddo
 
-      if (ifdp0dt) then
-         y(j) = w1p
-         if (ifrhs) y(j) = w1p
-      endif
-
       return
       end
-
diff --git a/core/cvode_driver.f b/core/cvode_driver.f
@@ -45,7 +45,6 @@ subroutine cv_setsize
            cv_nlocal = cv_nlocal + ntot
          endif
       enddo
-      if (ifdp0dt) cv_nlocal = cv_nlocal + 1
 
       ! check array size is large enough
       if (cv_nlocal .gt. cv_lysize) then
@@ -116,7 +115,7 @@ subroutine cv_init
                call cfill(cv_atol_(1,1,1,1,i-1),atol(i),ntot)
             endif
          enddo
-         call cvpack(cv_atol,cv_atol_,atol(2),.false.)
+         call cvpack(cv_atol,cv_atol_,.false.)
       endif
 
       ! initialize vector module
@@ -132,7 +131,7 @@ subroutine cv_init
       endif
 
       ! initialize cvode
-      call cvpack(y0,t,p0th,.false.)
+      call cvpack(y0,t,.false.)
       call fcvmalloc(time, y0, cv_meth, itmeth, cv_iatol,
      &               cv_rtol, cv_atol, iout, rout, ipar, rpar, ier)
       if (ier.ne.0) then
@@ -268,7 +267,7 @@ subroutine cdscal_cvode
 
       if (cv_itask.eq.3) then
         if(istep.gt.1) then 
-          call cvpack(y,t,p0th,.false.)
+          call cvpack(y,t,.false.)
           call fcvreinit(timef,y,cv_iatol,cv_rtol,cv_atol,ier)
           if (ier .ne. 0) then
             write(*,'(a,i3)') 'ABORT: fcvsetrin ier=', ier
@@ -287,7 +286,7 @@ subroutine cdscal_cvode
          if (nid.eq.0) then
             write(*,'(a)') ' Restart integrator and try again  ...'
          endif
-         call cvpack(y,t,p0th,.false.)
+         call cvpack(y,t,.false.)
          call fcvreinit(timef,y,cv_iatol,cv_rtol,cv_atol,ier)
          cv_timel = 0
          call cv_rstat
@@ -303,7 +302,7 @@ subroutine cdscal_cvode
       endif
 
       cv_istep = cv_istep + 1 
-      call cvunpack(t,p0th,y)
+      call cvunpack(t,y)
 
       ! restore velocities               
       call copy(vx,vx_,ntot)             
@@ -476,7 +475,7 @@ subroutine fcvfun (time_, y, ydot, ipar, rpar, ier)
          cv_timel = time          
       endif
 
-      call cvunpack(t,p0th,y)          
+      call cvunpack(t,y)          
 
       if(nid.eq.0 .and. loglevel.gt.2) write(6,*) 'fcvfun',
      $                                 ifdqj
@@ -529,7 +528,7 @@ subroutine fcvfun (time_, y, ydot, ipar, rpar, ier)
          endif
       enddo
 
-      call cvpack(ydot,ydott,dp0thdt,.true.)
+      call cvpack(ydot,ydott,.true.)
 
       tcvf = tcvf + dnekclock()-etime1 
       ncvf = ncvf + 1 
@@ -568,13 +567,16 @@ subroutine cv_settime
       cv_dtlag(3) = dtlag(3)
 
       call rzero(cv_abmsh,3)
-      call setabbd (cv_abmsh,cv_dtlag,nabmsh,1)
+      call setabbd(cv_abmsh,cv_dtlag,nabmsh,1)
       do i = 1,3
          cv_abmsh(i) = cv_dtNek*cv_abmsh(i) 
       enddo
 
       call rzero(cv_ab,3)
-      call setabbd (cv_ab,cv_dtlag,nbd,nbd)
+      call setabbd(cv_ab,cv_dtlag,nbd,nbd)
+
+      call rzero(cv_bd,4)
+      call setbd(cv_bd,cv_dtlag,nbd)
 
       return
       end

diff --git a/core/drive1.f b/core/drive1.f
@@ -145,7 +145,8 @@ subroutine nek_init(comm)
 
       call dofcnt
 
-      jp = 0            ! Set perturbation field count to 0 for baseline flow
+      jp = 0  ! Set perturbation field count to 0 for baseline flow
+      p0thn = p0th
 
       call in_situ_init()
 

diff --git a/core/drive2.f b/core/drive2.f
@@ -45,18 +45,16 @@ subroutine initdat
 
 c     Set default logicals
 
-      IFDOIT    = .FALSE.
-      IFCVODE   = .false.
-      IFEXPLVIS = .false.
+      ifdoit    = .false.
+      ifcvode   = .false.
+      ifexplvis = .false.
+      ifvvisp   = .true.
 
       ifsplit = .false.
       if (lx1.eq.lx2) ifsplit=.true.
 
       if_full_pres = .false.
 
-c     Turn off (on) diagnostics for communication
-      IFGPRNT= .FALSE.
-
       CALL RZERO (PARAM,200)
 C
 C     The initialization of CBC is done in READAT
@@ -84,6 +82,8 @@ subroutine initdat
       CALL RZERO(USRDIV,NTOT)
       CALL RZERO(QTL,NTOT)
 
+      NSTEPS = 0
+
       RETURN
       END
 C

diff --git a/core/makeq_aux.f b/core/makeq_aux.f
@@ -10,6 +10,12 @@ subroutine makeq_aux
 
       if (if_conv_std) call makeuq ! this will zero bq first
 
+      if (ifdp0dt .and.  ifield.eq.2 .and. .not.ifcvfld(ifield)) then
+            dd = dp0thdt * (gamma0 - 1.)/gamma0
+            ntot = nxyz*nelv
+            call add2s2 (bq(1,1,1,1,ifield-1),bm1,dd,ntot) 
+      endif
+
       if(filterType.eq.2) call make_hpf
 
       return

diff --git a/core/navier1.f b/core/navier1.f
@@ -1743,7 +1743,6 @@ subroutine setbd (bd,dtbd,nbd)
       INTEGER IR(ldim),IC(ldim)
       REAL BD(1),DTBD(1)
 C
-      CALL RZERO (BD,10)
       IF (NBD.EQ.1) THEN
          BD(1) = 1.
          BDF   = 1.

diff --git a/core/plan4.f b/core/plan4.f
@@ -150,7 +150,7 @@ subroutine crespsp (respr)
       call opadd2cm (wa1,wa2,wa3,ta1,ta2,ta3,scale)
 
 c compute stress tensor for ifstrs formulation - variable viscosity Pn-Pn
-      if (ifstrs) then
+      if (ifstrs .and. ifvvisp) then
          call opgrad   (ta1,ta2,ta3,vdiff)
          call invcol2  (ta1,vdiff,ntot1)
          call invcol2  (ta2,vdiff,ntot1)
@@ -596,20 +596,34 @@ subroutine userqtl_scig
          call copy(w2,w1,ntot)
          call col2(w1,bm1,ntot)
 
-         p0alph1 = p0th / glsum(w1,ntot)
+         p0alph1 = 1. / glsum(w1,ntot)
 
          call copy   (w1,qtl,ntot)
          call col2   (w1,bm1,ntot)
 
          termQ = glsum(w1,ntot)
          if (ifcvfun) then
             termV = glcflux(vx,vy,vz)
+            prhs  = p0alph1*(termQ - termV)
+            pcoef =(cv_bd(1) - cv_dtNek*prhs)
+            call add3s2(Saqpq,p0thn,p0thlag(1),cv_bd(2),cv_bd(3),1)
+            if(nbd.eq.3) call add2s2(Saqpq,p0thlag(2),cv_bd(4),1)
+            p0th  = Saqpq / pcoef
          else
             termV = glcflux(vx_e,vy_e,vz_e)
+            prhs  = p0alph1*(termQ - termV)
+            pcoef =(bd(1) - dt*prhs)
+            call add3s2(Saqpq,p0thn,p0thlag(1),bd(2),bd(3),1)
+            if(nbd.eq.3) call add2s2(Saqpq,p0thlag(2),bd(4),1)
+            p0th  = Saqpq / pcoef
+            p0thlag(2) = p0thlag(1)
+            p0thlag(1) = p0thn
+            p0thn      = p0th
          endif
-         dp0thdt = p0alph1*(termQ - termV)
 
-         dd =-dp0thdt/p0th
+         dp0thdt= prhs*p0th
+
+         dd =-prhs
          call cmult(w2,dd,ntot)
          call add2 (qtl,w2,ntot)
       endif