Merge pull request #42 from chowland/time-averaging

Add on-the-fly time averaging for 3D fields and energy spectra

.fortls

@@ -0,0 +1,5 @@
+{
+  "source_dirs": ["src/**"],
+  "excl_paths": ["src/old/**", "tools/**"],
+  "include_dirs": ["obj/**"]
+}

Makefile

@@ -10,7 +10,7 @@ FLAVOUR=GNU
 # MARENOSTRUM (Intel): fabric intel mkl impi hdf5 fftw szip
 # SUPERMUC (Intel): fftw hdf5
 # DISCOVERER:
-#	GNU: hdf5/1/1.14/1.14.0-gcc-openmpi fftw/3/latest-gcc-openmpi lapack
+#	GNU: hdf5/1/1.14/latest-gcc-openmpi fftw/3/latest-gcc-openmpi lapack
 #	Intel: hdf5/1/1.14/latest-intel-openmpi fftw/3/latest-gcc-openmpi mkl
 
 #=======================================================================
@@ -36,13 +36,14 @@ ifeq ($(MACHINE),PC)
 # Intel Debug Flags
 # FC += -O0 -g -traceback -check bounds
 	ifeq ($(FLAVOUR),GNU)
-		LDFLAGS = -L$(HOME)/fftw-install/lib -lfftw3 -llapack -ldl
+		LDFLAGS = -lfftw3 -llapack -ldl
 	else
 		LDFLAGS = -lfftw3 -qmkl=sequential
 	endif
 endif
 ifeq ($(MACHINE),DISCOVERER)
 	ifeq ($(FLAVOUR),GNU)
+		FC = h5pfc -cpp -fdefault-real-8 -fdefault-double-8
 		LDFLAGS += -lfftw3 -llapack -ldl
 	else
 		LDFLAGS += -lfftw3 -qmkl=sequential
@@ -94,7 +95,7 @@ OBJS = obj/main.o obj/CalcMaxCFL.o \
 	obj/MakeMovieYCut.o obj/MakeMovieZCut.o obj/MpiAuxRoutines.o \
 	obj/QuitRoutine.o obj/ReadInputFile.o obj/ResetLogs.o \
 	obj/SetTempBCs.o obj/SolveImpEqnUpdate_Temp.o obj/SolveImpEqnUpdate_X.o \
-	obj/SolveImpEqnUpdate_YZ.o obj/SpecRoutines.o \
+	obj/SolveImpEqnUpdate_YZ.o \
 	obj/TimeMarcher.o obj/WriteFlowField.o obj/WriteGridInfo.o \
 	obj/CalcWriteQ.o obj/GlobalQuantities.o obj/ReadFlowInterp.o
 
@@ -119,14 +120,15 @@ OBJS += obj/mean_zplane.o
 # Module object files
 MOBJS = obj/param.o obj/decomp_2d.o obj/AuxiliaryRoutines.o obj/decomp_2d_fft.o \
 	obj/pressure.o obj/HermiteInterpolations.o obj/grid.o obj/h5_tools.o obj/means.o \
-	obj/ibm_param.o obj/IBMTools.o obj/moisture.o obj/salinity.o obj/phasefield.o
+	obj/ibm_param.o obj/IBMTools.o obj/moisture.o obj/salinity.o obj/phasefield.o \
+	obj/time_averaging.o obj/spectra.o
 
 #=======================================================================
 #  Files that create modules:
 #=======================================================================
 MFILES = param.F90 decomp_2d.F90 AuxiliaryRoutines.F90 decomp_2d_fft.F90 \
 	pressure.F90 HermiteInterpolations.F90 grid.F90 ibm_param.F90 IBMTools.F90 \
-	moisture.F90 salinity.F90 phasefield.F90
+	moisture.F90 salinity.F90 phasefield.F90 time_averaging.F90 spectra.F90
 
 #============================================================================ 
 #  make PROGRAM   
@@ -169,6 +171,10 @@ $(OBJDIR)/phasefield.o: src/phasefield.F90 obj/salinity.o
 	$(FC) -c -o $@ $< $(LDFLAGS)
 $(OBJDIR)/moisture.o: src/moisture.F90
 	$(FC) -c -o $@ $< $(LDFLAGS)
+$(OBJDIR)/time_averaging.o: src/time_averaging.F90
+	$(FC) -c -o $@ $< $(LDFLAGS)
+$(OBJDIR)/spectra.o: src/spectra.F90 obj/time_averaging.o obj/pressure.o
+	$(FC) -c -o $@ $< $(LDFLAGS)
 $(OBJDIR)/%.o: src/%.F90 $(MOBJS)
 	$(FC) -c -o $@ $< $(LDFLAGS)
 $(OBJDIR)/%.o: src/flow_solver/%.F90 $(MOBJS)

docs/spectra.md

@@ -0,0 +1,23 @@
+# Power spectra and time-averaging
+
+## Spectrum calculation (`afid_spectra`)
+
+The module [`afid_spectra`](https://github.com/chowland/AFiD-MuRPhFi/blob/main/src/spectra.F90) is used to compute and store power spectra from the simulation fields.
+Since the computation of two-dimensional Fourier transforms involves the expensive operation of transposing the global 3-D arrays, spectra are not calculated by default.
+This option is activated by including a file `spectra.in` in the simulation directory that specifies the time from which spectra should be calculated.
+
+Fourier transforms are performed in the `y` and `z` directions in `CalcFourierCoef` to produce the discretized equivalent of the three-dimensional array
+
+$$
+\hat{f}(x,k_y, k_z) = \iint f(x,y,z) e^{-i (k_y y + k_z z)} \,\mathrm{d}y \,\mathrm{d}z
+$$
+
+Currently, the power spectra $\Phi_{ff}=|\hat{f}|^2$ of the three velocity components `vx`, `vy`, `vz` and temperature `temp` are computed and stored along with the co-spectrum $\Phi_{u\theta}=\hat{S}^* \hat{u}$ of the wall-normal concentration flux.
+Each spectrum is computed at each time step and integrated in time to provide a final time-averaged three-dimensional spectrum that is written to `outputdir` by `WriteSpectra`.
+
+## Time-averaged fields (`afid_averaging`)
+
+A similar procedure is followed by the module `afid_averaging` to provide three-dimensional fields of time-averaged variables.
+The same input file `spectra.in` is needed to activate this operation, and the averaging occurs from the time specified in the input file.
+The subroutine `UpdateTemporalAverages` integrates the four basic variables in time at each time step to construct the temporal averages.
+The time-averaged fields are eventually written to `outputdir/fields/xx_mean.h5`.

examples/RayleighBenardConvection/spectra.in

@@ -1,3 +1,2 @@
-0.00132	0.01321	0.04028	0.012731165	0.5
-! This file contains 5 values of wall-normal
-! position from which to generate 1D spectra
+! Sim time at which to start averaging
+100.0

mkdocs.yml

@@ -59,4 +59,5 @@ nav:
         # - 'buoyancy.md'
         - 'interpolation.md'
         - 'stretched_grids.md'
-        - 'ibm.md'
+        - 'ibm.md'
+        - 'spectra.md'

src/flow_solver/CheckDivergence.F90

@@ -13,7 +13,7 @@ subroutine CheckDivergence(qmax,qmaxr)
     use local_arrays, only: vy,vx,vz
     use afid_salinity, only: vyr,vxr,vzr
     use mpih
-    use decomp_2d, only: xstart,xend,xstartr,xendr,nrank
+    use decomp_2d, only: xstart,xend,xstartr,xendr!,nrank
     implicit none
     real,intent(out) :: qmax,qmaxr
     integer :: jc,kc,kp,jp,ic,ip

src/flow_solver/CreateGrid.F90

@@ -14,10 +14,7 @@ subroutine CreateGrid
     use GridModule
     implicit none
 
-    real :: x1,x2,x3
-
-    integer :: i, j, kc, km, kp
-    logical :: fexist
+    integer :: kc
 
     do kc=1,nxm
         kmv(kc)=kc-1

src/flow_solver/CreateInitialConditions.F90

@@ -16,11 +16,10 @@ subroutine CreateInitialConditions
     use afid_salinity, only: RayS
     use afid_phasefield, only: pf_eps, read_phase_field_params
     implicit none
-    integer :: j,k,i,kmid, io
+    integer :: j,k,i,kmid
     real :: xxx,yyy,zzz,eps,varptb,amp
     real :: h0,t0,Lambda,r, x0, A, B, alpha
     real, dimension(11) :: yh, zh
-    logical :: exists
 
     call random_seed()
 
@@ -37,15 +36,29 @@ subroutine CreateInitialConditions
     end if
 
     if (gAxis == 1) then
-        if ((RayT < 0) .and. (RayS <0)) then
-            !CJH: Stratified shear layer initial condition
-            do i=xstart(3),xend(3)
-                do j=xstart(2),xend(2)
-                    do k=1,nxm
-                        vy(k,j,i) = tanh(xm(k) - alx3/2.0)
+        if (RayT < 0) then
+            if (RayS <0) then
+                !CJH: Stratified shear layer initial condition
+                do i=xstart(3),xend(3)
+                    do j=xstart(2),xend(2)
+                        do k=1,nxm
+                            vy(k,j,i) = tanh(xm(k) - alx3/2.0)
+                            ! vz(k,j,i) = 1.0/cosh(xm(k) - alx3/2.0)
+                        end do
                     end do
                 end do
-            end do
+            else
+                !CJH: Salt-fingering initial condition
+                do i=xstart(3),xend(3)
+                    do j=xstart(2),xend(2)
+                        do k=1,nxm
+                            vx(k,j,i) = 0.0
+                            vy(k,j,i) = 0.0
+                            vz(k,j,i) = 0.0
+                        end do
+                    end do
+                end do
+            end if
         else
             !CJH: RBC initial condition as used in AFiD 1.0
             eps = 0.01
@@ -178,8 +191,8 @@ subroutine CreateInitialConditions
                 do k=1,nxm
                     call random_number(varptb)
                     if (abs(xm(k)-0.5) + eps > 0.5) then
-                    amp = 0.5 - abs(xm(k)-0.5) ! CJH Prevent values of |T| exceeding 0.5
-                    temp(k,j,i) = temp(k,j,i) + amp*(2.d0*varptb - 1.d0)
+                        amp = 0.5 - abs(xm(k)-0.5) ! CJH Prevent values of |T| exceeding 0.5
+                        temp(k,j,i) = temp(k,j,i) + amp*(2.d0*varptb - 1.d0)
                     else
                     temp(k,j,i) = temp(k,j,i) + eps*(2.d0*varptb - 1.d0)
                     end if

src/flow_solver/ExplicitTermsTemp.F90

@@ -14,7 +14,7 @@ subroutine ExplicitTermsTemp
     use decomp_2d, only: xstart,xend
     implicit none
     integer :: jc,kc,ic
-    integer :: km,kp,jm,jp,im,ip
+    integer :: jm,jp,im,ip
     real    :: htx,hty,htz,udy,udz
     real    :: udzq,udyq
     real    :: dyyt,dzzt

src/flow_solver/ImplicitAndUpdateTemp.F90

@@ -17,9 +17,7 @@ subroutine ImplicitAndUpdateTemp
     use ibm_param
     implicit none
     integer :: jc,kc,ic
-    integer :: km,kp,ke
     real    :: alpec,dxxt
-    real    :: app,acc,amm
 
     alpec=al/pect
 

src/flow_solver/ImplicitAndUpdateVX.F90

@@ -18,7 +18,7 @@ subroutine ImplicitAndUpdateVX
     use ibm_param
     implicit none
     integer :: jc,kc
-    integer :: km,kp,ic,ke
+    integer :: km,kp,ic
     real    :: alre,udx3
     real    :: amm,acc,app,dxp,dxxvx
 
@@ -67,7 +67,7 @@ subroutine ImplicitAndUpdateVX
     end do
 !$OMP END PARALLEL DO
 
-    iF(ANY(IsNaN(rhs))) write(*,*) 'NaN in rhs'
+    ! iF(ANY(IsNaN(rhs))) write(*,*) 'NaN in rhs'
 
 !  Solve equation and update velocity
 

src/flow_solver/ImplicitAndUpdateVY.F90

@@ -18,7 +18,7 @@ subroutine ImplicitAndUpdateVY
     use ibm_param
     implicit none
     integer :: kc,jmm,jc,ic
-    integer :: kpp,kmm,ke
+    integer :: kpp,kmm
     real    :: alre,udy
     real    :: amm,acc,app
     real    :: dyp,dxxvy

src/flow_solver/ImplicitAndUpdateVZ.F90

@@ -18,7 +18,7 @@ subroutine ImplicitAndUpdateVZ
     use ibm_param
     implicit none
     integer :: kc,jc,ic,imm
-    integer :: kmm,kpp,ke
+    integer :: kmm,kpp
     real    :: alre,amm,acc,app,udz
     real    :: dxxvz,dzp
 

src/flow_solver/QuitRoutine.F90

@@ -17,6 +17,8 @@ subroutine QuitRoutine(tin,normalexit,errorcode)
       use afid_moisture, only: DeallocateMoistVariables
       use afid_phasefield, only: DeallocatePFVariables
       use afid_salinity, only: DeallocateSalVariables
+      use afid_averaging
+      use afid_spectra
       implicit none
       logical, intent(in) :: normalexit
       integer :: errorcode
@@ -33,6 +35,10 @@ subroutine QuitRoutine(tin,normalexit,errorcode)
         if(nrank.eq.0) write(6,'(a,f10.2,a)') '  Total Iteration Time = ',(tin(3) -tin(2))/3600.0,' h.'
         ! if (statcal) call WriteStats
         call WriteFlowField(.true.)
+        if (specwrite) then
+          call WriteTemporalAverages
+          call WriteSpectra
+        end if
       else
         call MPI_Abort(MPI_COMM_WORLD, 1, ierr)
       endif
@@ -44,6 +50,10 @@ subroutine QuitRoutine(tin,normalexit,errorcode)
       if (phasefield) call DeallocatePFVariables
       if (IBM) call DeallocateIBMVariables
       if (moist) call DeallocateMoistVariables
+      if (specwrite) then
+        call DeallocateAveragingVariables
+        call DeallocateSpectra
+      end if
       call HdfClose
       call decomp_2d_fft_finalize
       call decomp_2d_finalize

src/flow_solver/ReadInputFile.F90

@@ -13,10 +13,9 @@ subroutine ReadInputFile
     use afid_phasefield, only: pf_A, pf_D, pf_eps, pf_Lambda, pf_S, pf_Tm
     implicit none
     character(len=4) :: dummy
-    integer flagstat,flagbal,flagmelt
+    integer :: flagmelt
     integer :: flagMR, flagsal, flagPF
-    integer :: FFscaleS!, pf_IBM
-    logical fexist
+    integer :: FFscaleS
 
     open(unit=15,file='bou.in',status='old')
         read(15,301) dummy

src/grid.F90

@@ -227,7 +227,7 @@ subroutine scallop_grid(c_grd, m_grd, Nm, grd_len, Retau, dw)
         real, intent(in) :: grd_len, Retau, dw
 
         integer :: k, ks
-        real :: alpha, kb, dyw, sig, dxlo, dxup, dxsmooth
+        real :: alpha, kb, sig, dxlo, dxup, dxsmooth
 
         ! Index of roughness height
         kb = 0.2*Retau/dw

src/h5tools/MakeMovieXCut.F90

@@ -10,13 +10,13 @@
 subroutine Mkmov_xcut
     use mpih
     use hdf5
-    use decomp_2d, only: xstart,xend,xstartr,xendr,DECOMP_2D_COMM_CART_X
+    use decomp_2d, only: xstart,xend,xstartr,xendr!,DECOMP_2D_COMM_CART_X
     use local_arrays, only: vz,vy,vx,temp
-    use afid_salinity, only: sal
-    use afid_phasefield, only: phi
+    use afid_salinity, only: sal, RayS
+    ! use afid_phasefield, only: phi
     use afid_moisture, only: humid
     use h5_tools
-    use param, only: nxm, nxmr, IBM
+    use param, only: nxm, nxmr, IBM, RayT
     implicit none
     character(70) :: filename
     character(4) :: varname
@@ -28,6 +28,7 @@ subroutine Mkmov_xcut
 
     ! Select plane - plane next to lower wall
     ic = 1
+    if (RayS < 0 .and. RayT < 0) ic = nxm/2
     if (IBM) ic = nxm/2
 
     ! Record filename as string
@@ -64,6 +65,7 @@ subroutine Mkmov_xcut
         !! Repeat on refined grid to save salinity
         ! Select wall plane index for refined grid
         ic = 1
+        if (RayS < 0 .and. RayT < 0) ic = nxmr/2
         if (IBM) ic = nxmr/2
 
         call h5_open_or_create(file_id, filename, comm, fexist)

src/h5tools/MakeMovieYCut.F90

@@ -11,7 +11,7 @@ subroutine Mkmov_ycut
     use param, only: nym, nymr
     use mpih
     use hdf5
-    use decomp_2d, only: xstart,xend,xstartr,xendr,DECOMP_2D_COMM_CART_X
+    use decomp_2d, only: xstart,xend,xstartr,xendr!,DECOMP_2D_COMM_CART_X
     use local_arrays, only: vz,vy,vx,temp
     use afid_salinity, only: sal
     use afid_phasefield, only: phi

src/h5tools/MakeMovieZCut.F90

@@ -11,11 +11,10 @@ subroutine Mkmov_zcut
     use param, only: nzm, nzmr
     use mpih
     use hdf5
-    use decomp_2d, only: xstart,xend,xstartr,xendr,DECOMP_2D_COMM_CART_X
-    use local_arrays, only: vz,vy,vx,temp, pr
-    ! use mgrd_arrays, only: sal, phi, phic, tempr
+    use decomp_2d, only: xstart,xend,xstartr,xendr!,DECOMP_2D_COMM_CART_X
+    use local_arrays, only: vz,vy,vx,temp!, pr
     use afid_salinity, only: sal
-    use afid_phasefield, only: phi, phic, tempr
+    use afid_phasefield, only: phi!, phic, tempr
     use afid_moisture, only: humid
     use h5_tools
     implicit none