Merge pull request #27 from matt-frey/domain-averaged

Calculate domain-averaged kinetic energy, potential energy and enstrophy
matt-frey · Dec 4, 2022 · 7ce2372 · 7ce2372
2 parents 6217fc6 + 4e56705
commit 7ce2372
Show file tree

Hide file tree

Showing 8 changed files with 22 additions and 44 deletions.
diff --git 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>", ke[0])
     #print("initial <EN>", 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
@@ -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>", ke[0])
     #print("initial <EN>", en[0])
 

diff --git 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>", ke[0])
     #print("initial <EN>", en[0])
 

diff --git 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
@@ -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
@@ -78,19 +78,19 @@ 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)
 
             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
@@ -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
@@ -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