Viscous dissipation becomes zero in post-processing mode with load_fld

[Chinthaka Jacob]

Dear Nek Communtiy,

  I am encountering a persistent issue when computing viscous dissipation (based on comp_gije --> comp_sije --> mag_tensor_e) in post-processing mode using checkpoint files generated in the Pn–Pn-2 formulation. This is for time-averaging viscous dissipation across many checkpoint files. The goal is to avoid rerunning the simulation and to perform consistent post-processing.

  In post-processing mode, I loop over multiple .f0000* files using load_fld() inside userchk. The velocity and scalar fields are correctly loaded (non-zero). However, all gradient-based quantities, such as SNRM, collapse to zero.

  I observed that the failure appears at the level of gradient evaluation, not data loading. Given that comp_gije internally relies on geometric factors such as rxm1, sxm1, txm1, and Jacobians, it seems likely that these are not properly initialized in post-processing mode when using load_fld(). I presume that load_fld() correctly populates the field variables (vx, vy, vz, t), but does not reinitialize geometry/metric operators, leading to:

\nabla{u}=0 ==> Sij​=0 ==> S:S=0.

Would adding a call to geom_reset(1) after each load_fld() call resolve this issue? 

1. Is it expected that load_fld() does not reinitialize geometric factors required by comp_gije?
2. Is call geom_reset(1) the correct and sufficient fix in this context?
3. Are there additional routines (e.g. related to derivative operators or mapping) that must be explicitly called in post-processing mode? 

Any clarification of minimal working examples would be highly appreciated. I suspect others may encounter the same issue when working with gradient-based diagnostics in post-processing mode.

Chinthaka.

Here is the snippet of my userchk:

c-----------------------------------------------------------------------
      subroutine userchk()

      include 'SIZE'
      include 'TOTAL'
      include 'RESTART'
     
      parameter(lt=lx1*ly1*lz1*lelv)     
      character*80 filename(9999)
      integer nfiles, n
      real mean_temp, mean_KExy, mean_KEz

c    work arrays for viscous dissipation
c    sij : strain tensor storage per element
c    snrm: scalar field holding S:S
c    one : vector of ones for inner product with bm1

      real sij(lx1*ly1*lz1,ldim,ldim,lelv)
      real snrm(lx1*ly1*lz1,lelv)
      real eta_sum(lx1*ly1*lz1,lelv)
      real one(lt)

c    scalar outputs
      real eps_int, eta_mean_int
      real scale, smin, smax

      ifxyo = .true.  
      n = nx1*ny1*nz1*nelv

c    initialize
      call cfill(one, 1.0, n)
     
      if (nid.eq.0) then
         open(unit=199,file='file.list',form='formatted',status='old')
         read(199,*) nfiles
         read(199,'(A80)') (filename(i),i=1,nfiles)
         do k = 1,nfiles
            print*, filename(k)
      enddo
         close(199)
      endif

      call bcast(nfiles,isize)
      call bcast(filename,nfiles*80)

      call rzero(eta_sum,n)
      do i = 1,nfiles
         call load_fld(filename(i))

c        Compute avg quantities
         mean_temp = glsc2(t,bm1,n)/volvm1
         mean_KExy = 0.5D0*(glsc3(vx,bm1,vx,n) +
     $                      glsc3(vy,bm1,vy,n))/volvm1
         mean_KEz = 0.5D0*glsc3(vz,bm1,vz,n)/volvm1
       
c        Save them
         if (nid .eq. 0) then !root process
           write(51,'(8E15.7)') time,mean_temp,mean_KExy,mean_KEz
         endif

        vol = glsc3(one,bm1,one,n)
        call geom_reset(1) ! reset geometric factors

c------------------- compute S:S pointwise---------------------

         do e=1,nelv
            call comp_gije(sij(1,1,1,e),
     $           vx(1,1,1,e), vy(1,1,1,e), vz(1,1,1,e), e)
            call comp_sije(sij(1,1,1,e))
            call mag_tensor_e(snrm(1,e), sij(1,1,1,e))
         enddo

         call cmult(snrm,2.0,n)

         smin = glmin(snrm,n)
         smax = glmax(snrm,n)

         eps_int = glsc3(snrm,bm1,snrm,n)

         if (nid .eq. 0) then
            write(6,'(A,I6)') 'Processed file ',i
            write(6,'(A,1PE15.7,A,1PE15.7)')
     $      '   snrm min/max = ',smin,'  ',smax
            write(6,'(A,1PE15.7)') '   eps_int     = ',eps_int
         endif

c        accumulate for time-average
         call add2(eta_sum,snrm,n)

      enddo

c------------------- compute time-average ------------------------

         scale = 1.0/(nfiles)
         call cmult(eta_sum,scale,n)

c        diagnostic averaged field
         eta_mean_int = glsc3(eta_sum,bm1,one,n)/volvm1
     
      return
      end
c-----------------------------------------------------------------------  

[YuHsiang Lan]

Hi Chinthaka,

Does the code behave when you run with nstep=1?

I typically just call exitt manually in userchk and never run under post-processing mode...

Also, does your mesh coordinates come from the checkpoint files?

If so, you will need to manually call geom_reset afterward.

The code doesn't re-compute them even the load_fld might have updated the mesh coordinates.

However, I'd still expect the mesh and those geometric factors has been computed based on the rea/re2.

Twice, actually. The first one is generated in gengeom and the second one is called from fix_geom after usrdat2.

The logfile should have something like these:

    generate geometry data

    regenerate geometry data

In fact, I just tested with eddy example that nstep=0 still have rx, ry, sx and sy at the first userchk.

      subroutine userchk
      include 'SIZE'    
      include 'TOTAL'

      n = lx1*ly1*lz1*nelv
      a1 = glamax(rxm1, n)
      a2 = glamax(rym1, n)
      a3 = glamax(sxm1, n)
      a4 = glamax(sym1, n)

      if (nio.eq.0) write(*,*)'dbg',istep,nsteps,a1,a2,a3,a4

      return

      end

log:

 nsteps=0 -> skip time loop
 running solver in post processing mode

 call userchk
 dbg           0           0  0.19634954084939693        5.8175686490358203E-014   4.4853010194856324E-014  0.19634954084939693

BTW, it's usually fine but accumulating fields like this can lose accuracy if you have, says, 1000 files. 

         call add2(eta_sum,snrm,n)

Unless the numbers are iid to zero mean, it's basically doing the round-off "big_number += small number" and you will lose significant digits from small number.

A safer approach is to follow avg_all

         beta  = dtime/atime
         alpha = 1.-beta
         ! compute averages E(X)
         call avg1    (uavg,vx,alpha,beta,ntot ,'um  ',ifverbose)

Hope this helps,

Yu-Hsiang

--

[Chinthaka Jacob]

Many thanks, Yu-Hsiang, for the detailed explanation! This resolved the issue cleanly, and the note on accumulation accuracy was particularly helpful going forward.

Chinthaka.