diff --git a/plotting/plot_grid_resolution.py b/plotting/plot_grid_resolution.py index 36991fe..a64c6ed 100644 --- a/plotting/plot_grid_resolution.py +++ b/plotting/plot_grid_resolution.py @@ -64,12 +64,6 @@ t, ke, en = np.loadtxt(os.path.join(fpath, 'beltrami_' + str(grid) + '_ecomp.asc'), skiprows=1, unpack=True) - voli = 1.0 / np.pi ** 3 - - # calculate mean KE and mean EN - ke *= voli - en *= voli - #print("initial ", ke[0]) #print("initial ", en[0]) @@ -115,11 +109,6 @@ os.path.join(fpath, 'beltrami_' + str(grid) + '_vorticity.asc'), skiprows=1, unpack=True) - voli = 1.0 / np.pi ** 3 - - # calculate mean KE and mean EN - en *= voli - maxen[i] = en.max() vmax[i] = vormax.max() diff --git a/plotting/plot_ke_en_varying_grid_resolution.py b/plotting/plot_ke_en_varying_grid_resolution.py index 92e2194..2e05bf5 100644 --- a/plotting/plot_ke_en_varying_grid_resolution.py +++ b/plotting/plot_ke_en_varying_grid_resolution.py @@ -54,12 +54,6 @@ t, ke, en = np.loadtxt(os.path.join(fpath, 'beltrami_' + str(grid) + '_ecomp.asc'), skiprows=1, unpack=True) - voli = 1.0 / np.pi ** 3 - - # calculate mean KE and mean EN - ke *= voli - en *= voli - #print("initial ", ke[0]) #print("initial ", en[0]) diff --git a/plotting/plot_ke_en_varying_prediss.py b/plotting/plot_ke_en_varying_prediss.py index 45c264e..2b0fd36 100644 --- a/plotting/plot_ke_en_varying_prediss.py +++ b/plotting/plot_ke_en_varying_prediss.py @@ -63,12 +63,6 @@ t, ke, en = np.loadtxt(os.path.join(fpath, 'beltrami_' + str(grid) + '_' + pred + '_ecomp.asc'), skiprows=1, unpack=True) - voli = 1.0 / np.pi ** 3 - - # calculate mean KE and mean EN - ke *= voli - en *= voli - #print("initial ", ke[0]) #print("initial ", en[0]) diff --git a/plotting/plot_ke_with_volume_render.py b/plotting/plot_ke_with_volume_render.py index f2ea24c..8b8b35c 100644 --- a/plotting/plot_ke_with_volume_render.py +++ b/plotting/plot_ke_with_volume_render.py @@ -93,9 +93,6 @@ def add_inset(ax, bounds, t, ke, step): t, ke, en = np.loadtxt(os.path.join(filepath, 'beltrami_' + str(grid) + '_ecomp.asc'), skiprows=1, unpack=True) -voli = 1.0 / np.pi ** 3 -ke = ke * voli - # calculate mean KE and mean EN ax.plot(t, ke, color='blue', linewidth=1.0) ax.set_xlabel(r'time, $t$') diff --git a/plotting/plot_peaks.py b/plotting/plot_peaks.py index fb30cf5..412a44d 100644 --- a/plotting/plot_peaks.py +++ b/plotting/plot_peaks.py @@ -55,11 +55,6 @@ _, vormax, _, _, _, _, _ = np.loadtxt(os.path.join(fpath, 'beltrami_' + str(grid) + '_vorticity.asc'), skiprows=1, unpack=True) - voli = 1.0 / np.pi ** 3 - - # calculate mean KE and mean EN - en *= voli - maxen[i] = en.max() vmax[i] = vormax.max() diff --git a/src/fields/field_diagnostics_netcdf.f90 b/src/fields/field_diagnostics_netcdf.f90 index 7a9bc6c..c4951cd 100644 --- a/src/fields/field_diagnostics_netcdf.f90 +++ b/src/fields/field_diagnostics_netcdf.f90 @@ -78,9 +78,9 @@ subroutine create_netcdf_field_stats_file(basename, overwrite) call define_netcdf_dataset( & ncid=ncid, & name='ke', & - long_name='kinetic energy', & + long_name='domain-averaged kinetic energy', & std_name='', & - unit='m^5/s^2', & + unit='m^2/s^2', & dtype=NF90_DOUBLE, & dimids=(/t_dim_id/), & varid=ke_id) @@ -88,9 +88,9 @@ subroutine create_netcdf_field_stats_file(basename, overwrite) call define_netcdf_dataset( & ncid=ncid, & name='en', & - long_name='enstrophy', & + long_name='domain-averaged enstrophy', & std_name='', & - unit='m^3/s^2', & + unit='1/s^2', & dtype=NF90_DOUBLE, & dimids=(/t_dim_id/), & varid=en_id) @@ -99,9 +99,9 @@ subroutine create_netcdf_field_stats_file(basename, overwrite) call define_netcdf_dataset( & ncid=ncid, & name='pe', & - long_name='potential energy', & + long_name='domain-averaged potential energy', & std_name='', & - unit='m^5/s^2', & + unit='m^2/s^2', & dtype=NF90_DOUBLE, & dimids=(/t_dim_id/), & varid=pe_id) diff --git a/src/fields/fields.f90 b/src/fields/fields.f90 index 6a4edd1..d7878d4 100644 --- a/src/fields/fields.f90 +++ b/src/fields/fields.f90 @@ -3,7 +3,7 @@ ! and functions. ! ============================================================================= module fields - use parameters, only : nx, ny, nz, vcell, ncelli, dx, lower, extent, ngrid + use parameters, only : nx, ny, nz, vcell, ncelli, dx, lower, extent, ngrid, vdomaini use constants, only : zero, f12, f14, one, two use merge_sort use inversion_utils, only : fftxys2p, diffx, diffy, central_diffz & @@ -115,9 +115,13 @@ subroutine calculate_peref zmean = zmean + gam * two * vcell peref = peref - b(k) * vcell * zmean enddo + + ! divide by domain volume to get domain-averaged peref + peref = peref * vdomaini #endif end subroutine calculate_peref + ! domain-averaged potential energy function get_potential_energy() result(pe) double precision :: pe #ifdef ENABLE_BUOYANCY @@ -135,14 +139,13 @@ function get_potential_energy() result(pe) enddo enddo - pe = pe * vcell - - pe = pe - peref + pe = pe * ncelli - peref #else pe = zero #endif end function get_potential_energy + ! domain-averaged kinetic energy function get_kinetic_energy() result(ke) double precision :: ke @@ -156,10 +159,11 @@ function get_kinetic_energy() result(ke) + vel(nz, :, :, 2) ** 2 & + vel(nz, :, :, 3) ** 2) - ke = ke * vcell + ke = ke * ncelli end function get_kinetic_energy + ! domain-averaged enstrophy function get_enstrophy() result(en) double precision :: en @@ -173,11 +177,12 @@ function get_enstrophy() result(en) + vor(nz, :, :, 2) ** 2 & + vor(nz, :, :, 3) ** 2) - en = en * vcell + en = en * ncelli end function get_enstrophy #ifdef ENABLE_BUOYANCY + ! domain-averaged grad(b) integral function get_gradb_integral() result(enb) double precision :: enb double precision :: ds(0:nz, 0:nx-1, 0:ny-1) @@ -207,7 +212,7 @@ function get_gradb_integral() result(enb) + f12 * sum(mag(0, :, :)) & + f12 * sum(mag(nz, :, :)) - enb = enb * vcell + enb = enb * ncelli end function get_gradb_integral #endif diff --git a/src/utils/parameters.f90 b/src/utils/parameters.f90 index 1170f80..b5defbb 100644 --- a/src/utils/parameters.f90 +++ b/src/utils/parameters.f90 @@ -37,6 +37,8 @@ module parameters ! domain size double precision :: extent(3) + double precision :: vdomaini + ! domain centre double precision :: center(3) @@ -87,6 +89,8 @@ subroutine update_parameters endif endif + vdomaini = one / product(extent) + vcell = product(dx) vcelli = one / vcell