diff --git a/.fortls b/.fortls new file mode 100644 index 00000000..7bfcbc46 --- /dev/null +++ b/.fortls @@ -0,0 +1,5 @@ +{ + "source_dirs": ["src/**"], + "excl_paths": ["src/old/**", "tools/**"], + "include_dirs": ["obj/**"] +} \ No newline at end of file diff --git a/Makefile b/Makefile index 08767ec8..5a01f543 100644 --- a/Makefile +++ b/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) diff --git a/docs/spectra.md b/docs/spectra.md new file mode 100644 index 00000000..3a477237 --- /dev/null +++ b/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`. \ No newline at end of file diff --git a/examples/RayleighBenardConvection/spectra.in b/examples/RayleighBenardConvection/spectra.in index 8d350088..e82b57eb 100644 --- a/examples/RayleighBenardConvection/spectra.in +++ b/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 \ No newline at end of file +! Sim time at which to start averaging +100.0 \ No newline at end of file diff --git a/mkdocs.yml b/mkdocs.yml index cd127a53..3656d3e3 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -59,4 +59,5 @@ nav: # - 'buoyancy.md' - 'interpolation.md' - 'stretched_grids.md' - - 'ibm.md' \ No newline at end of file + - 'ibm.md' + - 'spectra.md' \ No newline at end of file diff --git a/src/flow_solver/CheckDivergence.F90 b/src/flow_solver/CheckDivergence.F90 index 91458532..3ace7230 100644 --- a/src/flow_solver/CheckDivergence.F90 +++ b/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 diff --git a/src/flow_solver/CreateGrid.F90 b/src/flow_solver/CreateGrid.F90 index 477dc367..e5fdc0c2 100644 --- a/src/flow_solver/CreateGrid.F90 +++ b/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 diff --git a/src/flow_solver/CreateInitialConditions.F90 b/src/flow_solver/CreateInitialConditions.F90 index e0582e8e..14e2c294 100644 --- a/src/flow_solver/CreateInitialConditions.F90 +++ b/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 diff --git a/src/flow_solver/ExplicitTermsTemp.F90 b/src/flow_solver/ExplicitTermsTemp.F90 index e0de8c4f..4919c197 100644 --- a/src/flow_solver/ExplicitTermsTemp.F90 +++ b/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 diff --git a/src/flow_solver/ImplicitAndUpdateTemp.F90 b/src/flow_solver/ImplicitAndUpdateTemp.F90 index 80400bb0..6bd613cb 100644 --- a/src/flow_solver/ImplicitAndUpdateTemp.F90 +++ b/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 diff --git a/src/flow_solver/ImplicitAndUpdateVX.F90 b/src/flow_solver/ImplicitAndUpdateVX.F90 index 741a147c..f85ed39f 100644 --- a/src/flow_solver/ImplicitAndUpdateVX.F90 +++ b/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 diff --git a/src/flow_solver/ImplicitAndUpdateVY.F90 b/src/flow_solver/ImplicitAndUpdateVY.F90 index bffd435b..049d6ee2 100644 --- a/src/flow_solver/ImplicitAndUpdateVY.F90 +++ b/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 diff --git a/src/flow_solver/ImplicitAndUpdateVZ.F90 b/src/flow_solver/ImplicitAndUpdateVZ.F90 index b8862a98..c2c05e7f 100644 --- a/src/flow_solver/ImplicitAndUpdateVZ.F90 +++ b/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 diff --git a/src/flow_solver/QuitRoutine.F90 b/src/flow_solver/QuitRoutine.F90 index 1dccc9ab..b5608e21 100644 --- a/src/flow_solver/QuitRoutine.F90 +++ b/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 diff --git a/src/flow_solver/ReadInputFile.F90 b/src/flow_solver/ReadInputFile.F90 index 470a8716..4e0e1765 100644 --- a/src/flow_solver/ReadInputFile.F90 +++ b/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 diff --git a/src/grid.F90 b/src/grid.F90 index 67c09b96..6c273abb 100644 --- a/src/grid.F90 +++ b/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 diff --git a/src/h5tools/MakeMovieXCut.F90 b/src/h5tools/MakeMovieXCut.F90 index de690e8a..ad5a5aa8 100644 --- a/src/h5tools/MakeMovieXCut.F90 +++ b/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) diff --git a/src/h5tools/MakeMovieYCut.F90 b/src/h5tools/MakeMovieYCut.F90 index 8596f5af..a3f6e9bd 100644 --- a/src/h5tools/MakeMovieYCut.F90 +++ b/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 diff --git a/src/h5tools/MakeMovieZCut.F90 b/src/h5tools/MakeMovieZCut.F90 index 36e7a33f..34e97ffd 100644 --- a/src/h5tools/MakeMovieZCut.F90 +++ b/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 diff --git a/src/h5tools/h5_tools.F90 b/src/h5tools/h5_tools.F90 index 33466dd5..fbfd1744 100644 --- a/src/h5tools/h5_tools.F90 +++ b/src/h5tools/h5_tools.F90 @@ -1,8 +1,8 @@ module h5_tools use HDF5 - use mpih, only: MPI_INFO_NULL + use mpih, only: MPI_INFO_NULL, MPI_COMM_WORLD use param, only: salinity, phasefield, nx, nxm, nym, nzm, time, tframe, nxmr, nymr, nzmr, IBM, moist - use decomp_2d, only: xstart, xstartr + use decomp_2d, only: xstart, xstartr, xend implicit none integer(HID_T) :: group_id, plist_id @@ -71,7 +71,7 @@ subroutine write_H5_plane(file_id, varname, var, axis) character, intent(in) :: axis integer, parameter :: ndims = 2 - integer(HID_T) :: filespace, slabspace, memspace, dset + integer(HID_T) :: filespace, memspace, dset integer(HSIZE_T), dimension(2) :: dims, data_count, data_offset character(len=5) :: frame character(len=10) :: dsetname @@ -150,4 +150,155 @@ subroutine InitSliceCommunicators end subroutine InitSliceCommunicators +!> Write out a 3D array that does not have a halo +subroutine write_3D_array(filename, arr) + character(len=30), intent(in) :: filename + real, dimension(1:nx,xstart(2):xend(2),xstart(3):xend(3)), intent(in) :: arr + + integer(hid_t) :: file_id, filespace, slabspace, memspace, dset + integer(hsize_t), dimension(3) :: dims, data_count, data_offset + integer :: comm, ndims + + comm = MPI_COMM_WORLD + info = MPI_INFO_NULL + + ndims = 3 + + ! Set the size of the global array + dims(1)=nx + dims(2)=nym + dims(3)=nzm + + ! Create a dataspace for the (global) array in the file + call h5screate_simple_f(ndims, dims, filespace, hdf_error) + + ! Set the size of the data local to this process + data_count(1) = nx + data_count(2) = xend(2)-xstart(2)+1 + data_count(3) = xend(3)-xstart(3)+1 + + ! Set the offset position of the local data in the global array + data_offset(1) = 0 + data_offset(2) = xstart(2)-1 + data_offset(3) = xstart(3)-1 + + ! Create a property list and set it up with the MPI comm + call h5pcreate_f(H5P_FILE_ACCESS_F, plist_id, hdf_error) + call h5pset_fapl_mpio_f(plist_id, comm, info, hdf_error) + + ! Create the HDF5 file + call h5fcreate_f(filename, H5F_ACC_TRUNC_F, file_id, & + hdf_error, access_prp=plist_id) + ! Close the property list used to set the file access property + call h5pclose_f(plist_id, hdf_error) + + ! Create a dataset with name 'var' + call h5dcreate_f(file_id, 'var', H5T_NATIVE_DOUBLE, & + filespace, dset, hdf_error) + + ! Create a dataspace for the array in the (local) memory + call h5screate_simple_f(ndims, data_count, memspace, hdf_error) + + ! Get dataspace of file and select slab (subsection) from offsets + call h5dget_space_f(dset, slabspace, hdf_error) + call h5sselect_hyperslab_f(slabspace, H5S_SELECT_SET_F, & + data_offset, data_count, hdf_error) + + ! Create the property defining the collective write + call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdf_error) + call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, & + hdf_error) + + ! Write the dataset + call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, & + arr, dims, & + hdf_error, file_space_id = slabspace, mem_space_id = memspace, & + xfer_prp = plist_id) + + ! Close the property list, dataset, dataspace, and file + call h5pclose_f(plist_id, hdf_error) + call h5dclose_f(dset, hdf_error) + call h5sclose_f(filespace, hdf_error) + call h5sclose_f(memspace, hdf_error) + call h5sclose_f(slabspace, hdf_error) + call h5fclose_f(file_id, hdf_error) + +end subroutine write_3D_array + +!> Write out a 3D spectrum decomposed using spectral space +subroutine write_3D_spectrum(filename, arr) + use decomp_2d_fft + character(len=30), intent(in) :: filename + real, dimension(1:nxm,sp%xst(2):sp%xen(2),sp%xst(3):sp%xen(3)), intent(in) :: arr + + integer(hid_t) :: file_id, filespace, slabspace, memspace, dset + integer(hsize_t), dimension(3) :: dims, data_count, data_offset + integer :: comm, ndims + + comm = MPI_COMM_WORLD + info = MPI_INFO_NULL + + ndims = 3 + + ! Set the size of the global array + dims(1)=nxm + dims(2)=nym + dims(3)=nzm + + ! Create a dataspace for the (global) array in the file + call h5screate_simple_f(ndims, dims, filespace, hdf_error) + + ! Set the size of the data local to this process + data_count(1) = nxm + data_count(2) = sp%xen(2)-sp%xst(2)+1 + data_count(3) = sp%xen(3)-sp%xst(3)+1 + + ! Set the offset position of the local data in the global array + data_offset(1) = 0 + data_offset(2) = sp%xst(2)-1 + data_offset(3) = sp%xst(3)-1 + + ! Create a property list and set it up with the MPI comm + call h5pcreate_f(H5P_FILE_ACCESS_F, plist_id, hdf_error) + call h5pset_fapl_mpio_f(plist_id, comm, info, hdf_error) + + ! Create the HDF5 file + call h5fcreate_f(filename, H5F_ACC_TRUNC_F, file_id, & + hdf_error, access_prp=plist_id) + ! Close the property list used to set the file access property + call h5pclose_f(plist_id, hdf_error) + + ! Create a dataset with name 'var' + call h5dcreate_f(file_id, 'var', H5T_NATIVE_DOUBLE, & + filespace, dset, hdf_error) + + ! Create a dataspace for the array in the (local) memory + call h5screate_simple_f(ndims, data_count, memspace, hdf_error) + + ! Get dataspace of file and select slab (subsection) from offsets + call h5dget_space_f(dset, slabspace, hdf_error) + call h5sselect_hyperslab_f(slabspace, H5S_SELECT_SET_F, & + data_offset, data_count, hdf_error) + + ! Create the property defining the collective write + call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdf_error) + call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, & + hdf_error) + + ! Write the dataset + call h5dwrite_f(dset, H5T_NATIVE_DOUBLE, & + arr, dims, & + hdf_error, file_space_id = slabspace, mem_space_id = memspace, & + xfer_prp = plist_id) + + ! Close the property list, dataset, dataspace, and file + call h5pclose_f(plist_id, hdf_error) + call h5dclose_f(dset, hdf_error) + call h5sclose_f(filespace, hdf_error) + call h5sclose_f(memspace, hdf_error) + call h5sclose_f(slabspace, hdf_error) + call h5fclose_f(file_id, hdf_error) + +end subroutine write_3D_spectrum + end module h5_tools \ No newline at end of file diff --git a/src/h5tools/mean_zplane.F90 b/src/h5tools/mean_zplane.F90 index 3d329317..72a7dbf6 100644 --- a/src/h5tools/mean_zplane.F90 +++ b/src/h5tools/mean_zplane.F90 @@ -3,7 +3,7 @@ subroutine mean_zplane use param, only: nx 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 @@ -132,7 +132,7 @@ subroutine mean_yplane use param, only: nx 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 diff --git a/src/ibm/IBMTools.F90 b/src/ibm/IBMTools.F90 index 53ad72e8..1135f3df 100644 --- a/src/ibm/IBMTools.F90 +++ b/src/ibm/IBMTools.F90 @@ -12,7 +12,7 @@ subroutine calc_interface_height(ph, h) !< input variable real, intent(out) :: h(xstartr(2)-lvlhalo:,xstartr(3)-lvlhalo:) real :: dxx(nxmr) - integer :: i,j,k,kp + integer :: i,j,k h(:,:) = 0.0 diff --git a/src/ibm/topogr_ibm.F90 b/src/ibm/topogr_ibm.F90 index 825bf1ea..26f40500 100644 --- a/src/ibm/topogr_ibm.F90 +++ b/src/ibm/topogr_ibm.F90 @@ -421,7 +421,7 @@ subroutine topogr end do end do end do - iF(ANY(IsNaN(distb))) write(*,*) 'NaN in distb' + ! iF(ANY(IsNaN(distb))) write(*,*) 'NaN in distb' if(l.eq.1) then if(n.gt.mpun) & diff --git a/src/main.F90 b/src/main.F90 index 063be8cf..fee66d09 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -11,6 +11,8 @@ program AFiD use afid_moisture use afid_salinity use afid_phasefield + use afid_averaging + use afid_spectra use h5_tools, only: InitSliceCommunicators ! use stat_arrays, only: nstatsamples,vx_global,vy_global,vz_global @@ -24,7 +26,7 @@ program AFiD integer :: prow=0,pcol=0 integer :: lfactor,lfactor2 character(100) :: arg - logical :: nanexist, write_mean_planes=.true. + logical :: write_mean_planes=.true.!, nanexist ! real,allocatable,dimension(:,:) :: dummy,dscan,dbot ! integer :: comm,ierror,row_id,row_coords(2),ic,jc,kc @@ -109,7 +111,7 @@ program AFiD call WriteGridInfo inquire(file=trim("spectra.in"),exist=specwrite) - if (specwrite) call InitPowerSpec + ! if (specwrite) call InitPowerSpec !m=================================== !m=================================== @@ -184,6 +186,11 @@ program AFiD ! if (phasefield) call UpdateIBMLocation end if + if (specwrite) then + call InitAveragingVariables + call InitSpectra + end if + !EP Update all relevant halos call update_halo(vx,lvlhalo) call update_halo(vy,lvlhalo) @@ -218,7 +225,7 @@ program AFiD end if call CalcMeanProfiles - if (specwrite) call WritePowerSpec + ! if (specwrite) call WritePowerSpec if(ismaster) write(6,*) 'Write plane slices' call Mkmov_xcut call Mkmov_ycut @@ -298,6 +305,11 @@ program AFiD call TimeMarcher + if (specwrite .and. time > tav_start) then + call UpdateTemporalAverages + call UpdateSpectra + end if + time=time+dt if(mod(time,tout).lt.dt) then @@ -306,11 +318,11 @@ program AFiD write(6,'(a,ES11.4,a,i9,a,ES11.4)') ' T = ',time,' NTIME = ',ntime,' DT = ',dt endif call CalcMeanProfiles - if (specwrite) then - if (ismaster) write(*,*) "Writing power spectra" - call WritePowerSpec - if (ismaster) write(*,*) "Done writing power spectra" - end if + ! if (specwrite) then + ! if (ismaster) write(*,*) "Writing power spectra" + ! call WritePowerSpec + ! if (ismaster) write(*,*) "Done writing power spectra" + ! end if if(ismaster) then open(96,file='outputdir/cfl.out',status='unknown',position='append',access='sequential') write(96,769) ntime,time,dt,instCFL*dt!,vx_global,vy_global,vz_global diff --git a/src/moisture.F90 b/src/moisture.F90 index c5af3885..c40ed41e 100644 --- a/src/moisture.F90 +++ b/src/moisture.F90 @@ -150,7 +150,7 @@ end subroutine SetHumidityBCs subroutine CreateInitialHumidity use mpih integer :: ic, jc, kc - real :: rnum, r, r2, amp + real :: rnum, r2, amp real :: bz(nxm), qz(nxm) logical :: exists character(len=30) :: dsetname, filename diff --git a/src/flow_solver/SpecRoutines.F90 b/src/old/SpecRoutines.F90 similarity index 100% rename from src/flow_solver/SpecRoutines.F90 rename to src/old/SpecRoutines.F90 diff --git a/src/phasefield.F90 b/src/phasefield.F90 index 32f44595..b1390a82 100644 --- a/src/phasefield.F90 +++ b/src/phasefield.F90 @@ -77,9 +77,8 @@ end subroutine DeallocatePFVariables !> Set the initial state of the phase field subroutine CreateInitialPhase - use afid_salinity, only: sal, PraS - real :: A, B, alpha, t0, x0, h0 - integer :: i, j, k + ! use afid_salinity, only: sal, PraS + real :: h0 if (salinity) then !! Ice above salty water (1D_DDMelting example) @@ -101,6 +100,9 @@ subroutine CreateInitialPhase elseif (pf_IC==3) then call set_flat_interface(0.5, .true.) call add_temperature_mode(amp=0.1, ymode=2, zmode=2) + + else if (pf_IC==4) then + call set_ice_sphere(r0=0.1) !! 1D supercooling example elseif (pf_IC==5) then @@ -298,6 +300,40 @@ subroutine set_ice_disc(r0) end do end subroutine set_ice_disc +!> Add an ice sphere at the centre of the domain of radius r0 +!! with a corresponding temperature profile +subroutine set_ice_sphere(r0) + use local_arrays, only: temp + real, intent(in) :: r0 !! Initial disc radius + + real :: r + integer :: i, j, k + + ! Temperature field (0 in disc, 1 out of disc, tanh interface width approx 1e-2) + do i=xstart(3),xend(3) + do j=xstart(2),xend(2) + do k=1,nxm + r = sqrt((xm(k) - 0.5*alx3)**2 + (ym(j) - 0.5*ylen)**2 + (zm(i) - 0.5*zlen)**2) + ! temp(k,j,i) = 0.5*(1.0 + tanh(100.0*(r - r0))) + if (r > r0) then + temp(k,j,i) = 1.0 + else + temp(k,j,i) = 0.0 + end if + end do + end do + end do + ! Phase-field (0 out of disc, 1 in disc) + do i=xstartr(3),xendr(3) + do j=xstartr(2),xendr(2) + do k=1,nxmr + r = sqrt((xmr(k) - 0.5*alx3)**2 + (ymr(j) - 0.5*ylen)**2 + (zmr(i) - 0.5*zlen)**2) + phi(k,j,i) = 0.5*(1.0 - tanh(0.5*(r - r0)/pf_eps)) + end do + end do + end do +end subroutine set_ice_sphere + !> Compute the explicit terms for the phase-field evolution !! and store the result in `hphi` subroutine ExplicitPhase diff --git a/src/salinity.F90 b/src/salinity.F90 index ab82ad41..1a98ad00 100644 --- a/src/salinity.F90 +++ b/src/salinity.F90 @@ -143,7 +143,7 @@ subroutine CreateInitialSalinity !! Rayleigh-Taylor setup for pore-scale simulation if (IBM) then call SetSaltTwoLayer(h0=0.5*alx3, eps=1e-7, stable=.false.) - call AddSalinityNoise(amp=0.1, localised=.true., h0=0.5*alx3) + call AddSalinityNoise(amp=0.1, localised=.true., h0=0.5*alx3, extent=0.01) !! Bounded double-diffusive convection (begin with small amplitude noise + BLs) else if ((active_S==1) .and. (active_T==1) .and. (gAxis==1)) then @@ -160,7 +160,7 @@ subroutine CreateInitialSalinity !! Stratified shear layer setup else if ((RayS < 0) .and. (RayT < 0)) then call SetSaltTwoLayer(h0=0.5*alx3, eps=1.0, stable=.true.) - call AddSalinityNoise(amp=1e-2, localised=.true., h0=0.5*alx3) + call AddSalinityNoise(amp=1e-2, localised=.true., h0=0.5*alx3, extent=1.0) !! Default: linear profile + small noise (e.g. RBC, VC) else @@ -226,10 +226,11 @@ end subroutine SetSaltTwoLayer !> Add random noise to the salinity field, either locally at an interface !! or uniformly. In both cases, noise is limited such that the absolute value !! of salinity does not exceed 0.5 -subroutine AddSalinityNoise(amp, localised, h0) +subroutine AddSalinityNoise(amp, localised, h0, extent) real, intent(in) :: amp !! Amplitude of random noise logical, intent(in) :: localised !! Flag determining whether to add noise around interface or everywhere real, intent(in), optional :: h0 !! Height of interface if using localised noise + real, intent(in), optional :: extent !! Width of localised noise region integer :: i, j, k real :: a2, varptb @@ -242,7 +243,7 @@ subroutine AddSalinityNoise(amp, localised, h0) call random_number(varptb) !! Add noise locally if (localised) then - sal(k,j,i) = sal(k,j,i) + amp/cosh((xmr(k) - h0)/0.01)**2*varptb + sal(k,j,i) = sal(k,j,i) + amp/cosh((xmr(k) - h0)/extent)**2*varptb ! Restrict initial salinity field to [-0.5,0.5] sal(k,j,i) = min(0.5, sal(k,j,i)) sal(k,j,i) = max(-0.5, sal(k,j,i)) diff --git a/src/spectra.F90 b/src/spectra.F90 new file mode 100644 index 00000000..1922ea6e --- /dev/null +++ b/src/spectra.F90 @@ -0,0 +1,207 @@ +!> Module for calculating energy/power spectra of various quantities +!! Spectra are stored as time averages and written out at the end of a simulation +module afid_spectra + use fftw_params + use afid_pressure + use decomp_2d_fft + + real, allocatable, dimension(:,:,:) :: vx_spec !! Power spectrum of vx + real, allocatable, dimension(:,:,:) :: vy_spec !! Power spectrum of vy + real, allocatable, dimension(:,:,:) :: vz_spec !! Power spectrum of vz + real, allocatable, dimension(:,:,:) :: te_spec !! Power spectrum of temp + + real, allocatable, dimension(:,:,:) :: wSr_spec !! Cospectrum of wall-normal S-flux (real part) + real, allocatable, dimension(:,:,:) :: wSi_spec !! Cospectrum of wall-normal S-flux (imaginary part) + +contains + +!> Allocate memory for 3D arrays of energy spectra +subroutine InitSpectra + use AuxiliaryRoutines, only: AllocateReal3DArray + + allocate(fouvar1(sp%xst(1):sp%xen(1), & + sp%xst(2):sp%xen(2),sp%xst(3):sp%xen(3))) + allocate(fouvar2(sp%xst(1):sp%xen(1), & + sp%xst(2):sp%xen(2), sp%xst(3):sp%xen(3))) + + call AllocateReal3DArray(vx_spec,1,nxm,sp%xst(2),sp%xen(2),sp%xst(3),sp%xen(3)) + call AllocateReal3DArray(vy_spec,1,nxm,sp%xst(2),sp%xen(2),sp%xst(3),sp%xen(3)) + call AllocateReal3DArray(vz_spec,1,nxm,sp%xst(2),sp%xen(2),sp%xst(3),sp%xen(3)) + call AllocateReal3DArray(te_spec,1,nxm,sp%xst(2),sp%xen(2),sp%xst(3),sp%xen(3)) + + if (salinity) then + call AllocateReal3DArray(wSr_spec,1,nxm,sp%xst(2),sp%xen(2),sp%xst(3),sp%xen(3)) + call AllocateReal3DArray(wSi_spec,1,nxm,sp%xst(2),sp%xen(2),sp%xst(3),sp%xen(3)) + + wSr_spec = 0.0 + wSi_spec = 0.0 + end if + + vx_spec = 0.0 + vy_spec = 0.0 + vz_spec = 0.0 + te_spec = 0.0 + +end subroutine InitSpectra + +!> Free memory from arrays used to store energy spectra +!! (at end of simulation) +subroutine DeallocateSpectra + use AuxiliaryRoutines, only: DestroyReal3DArray + + call DestroyReal3DArray(vx_spec) + call DestroyReal3DArray(vy_spec) + call DestroyReal3DArray(vz_spec) + call DestroyReal3DArray(te_spec) + call DestroyReal3DArray(wSr_spec) + call DestroyReal3DArray(wSi_spec) + + deallocate(fouvar1, fouvar2) + +end subroutine DeallocateSpectra + +!> Perform a Fourier transform of the variable `var` in y and z, +!! returning the transformed complex array `fouvar` +subroutine CalcFourierCoef(var, fouvar) + real, dimension(1:nxm, xstart(2):xend(2), xstart(3):xend(3)), intent(in) :: var + complex, dimension(sp%xst(1):sp%xen(1), & + sp%xst(2):sp%xen(2), sp%xst(3):sp%xen(3)), intent(out) :: fouvar + + real :: inyzm + + inyzm = 1.0/(nzm*nym) + + ! Transpose to y-pencil to perform FFT in y locally + call transpose_x_to_y(var, ry1, ph) + ! Plan the FFTs if this is not yet done (should be done by init pressure) + if (.not. planned) call PlanFourierTransform + ! Perform the FFT in y + call dfftw_execute_dft_r2c(fwd_guruplan_y, ry1, cy1) + + ! Transpose to z-pencil to perform FFT in z locally + call transpose_y_to_z(cy1, cz1, sp) + ! Perform the FFT in z + call dfftw_execute_dft(fwd_guruplan_z, cz1, cz1) + ! Normalize the result (FFTW does not do this automatically) + cz1 = cz1*inyzm + + !CJH: Is this really necessary? Could we work with z-pencils? + call transpose_z_to_x(cz1, fouvar, sp) + +end subroutine CalcFourierCoef + +!> Calculate the power spectrum of the variable `var` as a function of +!! height and horizontal wavenumbers +subroutine AddRealSpectrum(var, spectrum) + use param, only: dt + real, dimension(1:nxm,xstart(2):xend(2),xstart(3):xend(3)), intent(in) :: var + real, dimension(1:nxm,sp%xst(2):sp%xen(2),sp%xst(3):sp%xen(3)), intent(inout) :: spectrum + integer :: i, j, k + + ! Compute FFT on the variable + call CalcFourierCoef(var, fouvar1) + + do i=sp%xst(3),sp%xen(3) + do j=sp%xst(2),sp%xen(2) + do k=1,nxm + spectrum(k,j,i) = spectrum(k,j,i) + dt*abs(fouvar1(k,j,i))**2 + end do + end do + end do + +end subroutine AddRealSpectrum + +!> Calculate the cospectrum of the variables `var1` and `var2` as a function of +!! height and horizontal wavenumbers and add the real part to `rspec` and the +!! imaginary part to ispec (scaled by dt) +subroutine AddCospectrum(var1, var2, rspec, ispec) + use param, only: dt + real, dimension(1:nxm,xstart(2):xend(2),xstart(3):xend(3)), intent(in) :: var1, var2 + real, dimension(1:nxm,sp%xst(2):sp%xen(2),sp%xst(3):sp%xen(3)), intent(inout) :: rspec, ispec + complex :: cspec + integer :: i, j, k + + ! Compute FFT on the variable + call CalcFourierCoef(var1, fouvar1) + call CalcFourierCoef(var2, fouvar2) + + do i=sp%xst(3),sp%xen(3) + do j=sp%xst(2),sp%xen(2) + do k=1,nxm + cspec = fouvar1(k,j,i)*conjg(fouvar2(k,j,i)) + rspec(k,j,i) = rspec(k,j,i) + dt*real(cspec) + ispec(k,j,i) = ispec(k,j,i) + dt*aimag(cspec) + end do + end do + end do + +end subroutine AddCospectrum + +!> Update the spectra +subroutine UpdateSpectra + use local_arrays, only: vx, vy, vz, temp, dph, dq + use afid_salinity, only: salc + + call AddRealSpectrum(vx(1:nxm,xstart(2):xend(2),xstart(3):xend(3)), vx_spec) + call AddRealSpectrum(vy(1:nxm,xstart(2):xend(2),xstart(3):xend(3)), vy_spec) + call AddRealSpectrum(vz(1:nxm,xstart(2):xend(2),xstart(3):xend(3)), vz_spec) + call AddRealSpectrum(temp(1:nxm,xstart(2):xend(2),xstart(3):xend(3)), te_spec) + + ! Interpolate vx and salc to the cell centre, temporarily storing in dph + do i=xstart(3),xend(3) + do j=xstart(2),xend(2) + do k=1,nxm + dph(k,j,i) = 0.5*(vx(k,j,i) + vx(k+1,j,i)) + dq(k,j,i) = 0.5*(salc(k,j,i) + salc(k,j+1,i)) + end do + end do + end do + + if (salinity) then + call AddCospectrum(dq(1:nxm,xstart(2):xend(2),xstart(3):xend(3)), & + dph(1:nxm,xstart(2):xend(2),xstart(3):xend(3)), & + wSr_spec, wSi_spec) + end if + +end subroutine UpdateSpectra + +!> Normalize the time-averaged spectra and write them out to files +subroutine WriteSpectra + use h5_tools, only: write_3D_spectrum + use afid_averaging, only: tav_start + real :: tav_length !! Duration of time averaging interval + real :: i_tav + character(len=30) :: filename + + tav_length = time - tav_start + i_tav = 1.0/tav_length + + vx_spec = vx_spec*i_tav + vy_spec = vy_spec*i_tav + vz_spec = vz_spec*i_tav + te_spec = te_spec*i_tav + if (salinity) then + wSr_spec = wSr_spec*i_tav + wSi_spec = wSi_spec*i_tav + end if + + ! Save the arrays + filename = 'outputdir/vx_spectrum.h5' + call write_3D_spectrum(filename, vx_spec) + filename = 'outputdir/vy_spectrum.h5' + call write_3D_spectrum(filename, vy_spec) + filename = 'outputdir/vz_spectrum.h5' + call write_3D_spectrum(filename, vz_spec) + filename = 'outputdir/te_spectrum.h5' + call write_3D_spectrum(filename, te_spec) + + if (salinity) then + filename = 'outputdir/wSr_spectrum.h5' + call write_3D_spectrum(filename, wSr_spec) + filename = 'outputdir/wSi_spectrum.h5' + call write_3D_spectrum(filename, wSi_spec) + end if + +end subroutine WriteSpectra + +end module afid_spectra \ No newline at end of file diff --git a/src/time_averaging.F90 b/src/time_averaging.F90 new file mode 100644 index 00000000..42ba6e6c --- /dev/null +++ b/src/time_averaging.F90 @@ -0,0 +1,111 @@ +!> Module providing the ability to store time-averaged data that +!! is updated at each time step of the simulation +module afid_averaging + use param, only: nx, nxm, time, dt + use local_arrays, only: vx, vy, vz, temp + use decomp_2d, only: xstart, xend + + real, allocatable, dimension(:,:,:) :: vx_tav !! Time average of x-component of velocity + real, allocatable, dimension(:,:,:) :: vy_tav !! Time average of y-component of velocity + real, allocatable, dimension(:,:,:) :: vz_tav !! Time average of z-component of velocity + real, allocatable, dimension(:,:,:) :: te_tav !! Time average of temperature field + + real :: tav_start !! Time at which time-averaging begins + +contains + +!> Allocate the memory for the arrays storing temporal averages +subroutine InitAveragingVariables + use AuxiliaryRoutines, only: AllocateReal3DArray + logical :: exists + integer :: io + + ! Note we don't need ghost cells for these array + call AllocateReal3DArray(vx_tav,1,nx,xstart(2),xend(2),xstart(3),xend(3)) + call AllocateReal3DArray(vy_tav,1,nx,xstart(2),xend(2),xstart(3),xend(3)) + call AllocateReal3DArray(vz_tav,1,nx,xstart(2),xend(2),xstart(3),xend(3)) + call AllocateReal3DArray(te_tav,1,nx,xstart(2),xend(2),xstart(3),xend(3)) + + vx_tav = 0.0 + vy_tav = 0.0 + vz_tav = 0.0 + te_tav = 0.0 + + ! Read in time at which to start averaging + inquire(file="spectra.in", exist=exists) + if (exists) then + open(file="spectra.in", newunit=io, status="old", action="read") + read(io,*) + read(io,*) tav_start + end if + + ! If starting a simulation already past the specified averaging + ! time, set tav_start as current time + if (time > tav_start) then + tav_start = time + end if + +end subroutine InitAveragingVariables + +!> Free memory from arrays used to store temporal averages +!! (at end of simulation) +subroutine DeallocateAveragingVariables + use AuxiliaryRoutines, only: DestroyReal3DArray + + call DestroyReal3DArray(vx_tav) + call DestroyReal3DArray(vy_tav) + call DestroyReal3DArray(vz_tav) + call DestroyReal3DArray(te_tav) + +end subroutine DeallocateAveragingVariables + +!> Add dt*f(t) to each array used for time-averaging +subroutine UpdateTemporalAverages + integer :: i, j, k + + do i=xstart(3),xend(3) + do j=xstart(2),xend(2) + do k=1,nxm + vx_tav(k,j,i) = vx_tav(k,j,i) + dt*vx(k,j,i) + vy_tav(k,j,i) = vy_tav(k,j,i) + dt*vy(k,j,i) + vz_tav(k,j,i) = vz_tav(k,j,i) + dt*vz(k,j,i) + te_tav(k,j,i) = te_tav(k,j,i) + dt*temp(k,j,i) + end do + end do + end do + +end subroutine UpdateTemporalAverages + +!> Normalize the arrays to complete the time-averaging +!! and write them out to fields/var_mean.h5 +subroutine WriteTemporalAverages + use h5_tools, only: write_3D_array + real :: tav_length !! Duration of time averaging interval + real :: i_tav !! Inverse of averaging duration + character(len=30) :: filename + + tav_length = time - tav_start + i_tav = 1.0/tav_length + + ! Only write out the array if we've started averaging + if (tav_length > 0) then + ! Normalize the arrays + vx_tav = vx_tav*i_tav + vy_tav = vy_tav*i_tav + vz_tav = vz_tav*i_tav + te_tav = te_tav*i_tav + + ! Save the arrays + filename = 'outputdir/fields/vx_mean.h5' + call write_3D_array(filename, vx_tav) + filename = 'outputdir/fields/vy_mean.h5' + call write_3D_array(filename, vy_tav) + filename = 'outputdir/fields/vz_mean.h5' + call write_3D_array(filename, vz_tav) + filename = 'outputdir/fields/te_mean.h5' + call write_3D_array(filename, te_tav) + end if + +end subroutine WriteTemporalAverages + +end module afid_averaging \ No newline at end of file diff --git a/submit_scripts/submit_discoverer.sh b/submit_scripts/submit_discoverer.sh index 79f4e47c..fe55dc94 100644 --- a/submit_scripts/submit_discoverer.sh +++ b/submit_scripts/submit_discoverer.sh @@ -10,7 +10,7 @@ #SBATCH --exclusive module purge -module load hdf5/1/1.14/1.14.0-gcc-openmpi +module load hdf5/1/1.14/latest-gcc-openmpi module load lapack/latest-gcc module load fftw/3/latest-gcc-openmpi diff --git a/tools/afidtools/moisture.py b/tools/afidtools/moisture.py new file mode 100644 index 00000000..dd4b4aff --- /dev/null +++ b/tools/afidtools/moisture.py @@ -0,0 +1,17 @@ +# Tools specific to post-processing moist simulations + +class MoistParams: + """ + Class containing input parameters for moist convection, + reading the input from the file `humid.in` + """ + def __init__(self, folder): + fname = folder+"/humid.in" + with open(fname, "r") as f: + flist = f.readlines() + self.alpha = float(flist[4]) + self.beta = float(flist[7]) + self.gamma = float(flist[10]) + self.tau = float(flist[13]) + self.qfixN, self.qfixS = [int(x) for x in flist[20].split()] + self.Sm = float(flist[23]) \ No newline at end of file