Output global monitoring quantities to a csv file

docs/source/user/post-processing/index.rst

@@ -7,4 +7,5 @@ simulation output.
 .. toctree::
    :maxdepth: 2
 
+   monitoring.rst
    paraview.rst

docs/source/user/post-processing/monitoring.rst

@@ -0,0 +1,93 @@
+Monitoring runs
+===============
+
+During a simulation x3d2 writes global scalar quantities to a
+``monitoring.csv`` file at every output step (controlled by ``n_output``
+in the input file). This file is useful for tracking the evolution of
+the flow and for comparison with reference data.
+
+Output format
+-------------
+
+The file is a comma-separated CSV with a header line:
+
+.. code-block:: text
+
+   # time, enstrophy, ke, eps, div_u_max, div_u_mean
+
+Each subsequent row contains one record per output step.
+
+Quantities
+----------
+
+.. list-table::
+   :header-rows: 1
+   :widths: 20 40 40
+
+   * - Column
+     - Description
+     - Formula
+   * - ``time``
+     - Simulation time
+     - :math:`t`
+   * - ``enstrophy``
+     - Spatially-averaged enstrophy
+     - :math:`\mathcal{E} = \frac{1}{2N} \sum |\nabla \times \mathbf{u}|^2`
+   * - ``ke``
+     - Spatially-averaged kinetic energy
+     - :math:`E_k = \frac{1}{2N} \sum (u^2 + v^2 + w^2)`
+   * - ``eps``
+     - Dissipation rate (2nd derivative form)
+     - :math:`\varepsilon = -\frac{\nu}{N} \sum \mathbf{u} \cdot \nabla^2 \mathbf{u}`
+   * - ``div_u_max``
+     - Maximum of :math:`|\nabla \cdot \mathbf{u}|`
+     - Divergence-free check
+   * - ``div_u_mean``
+     - Mean of :math:`|\nabla \cdot \mathbf{u}|`
+     - Divergence-free check
+
+.. note::
+
+   The dissipation rate uses the second-derivative form, which is
+   equivalent to the strain-rate form for periodic flows.
+
+.. note::
+
+   The divergence columns (``div_u_max`` and ``div_u_mean``) measure how
+   well the incompressibility constraint :math:`\nabla \cdot \mathbf{u} = 0`
+   is satisfied.  Since x3d2 enforces the incompressibility constraint via a pressure
+   projection, these values should remain
+   close to machine precision (typically :math:`\sim 10^{-14}` in double
+   precision).  A growing divergence indicates a problem with the pressure
+   solver or time-stepping stability.
+
+Example: plotting with Python
+-----------------------------
+
+.. code-block:: python
+
+   import pandas as pd
+   import matplotlib.pyplot as plt
+
+   df = pd.read_csv("monitoring.csv", comment="#")
+   df.columns = df.columns.str.strip()
+
+   fig, axes = plt.subplots(2, 2, figsize=(10, 6))
+
+   axes[0, 0].plot(df["time"], df["enstrophy"])
+   axes[0, 0].set_ylabel("Enstrophy")
+
+   axes[0, 1].plot(df["time"], df["ke"])
+   axes[0, 1].set_ylabel("Kinetic energy")
+
+   axes[1, 0].plot(df["time"], df["eps"])
+   axes[1, 0].set_ylabel("Dissipation rate")
+
+   axes[1, 1].plot(df["time"], df["div_u_max"])
+   axes[1, 1].set_ylabel("div(u) max")
+
+   for ax in axes.flat:
+       ax.set_xlabel("Time")
+
+   plt.tight_layout()
+   plt.savefig("monitoring.png")

src/CMakeLists.txt

@@ -7,7 +7,8 @@ set(SRC
   ordering.f90
   poisson_fft.f90
   solver.f90
-  postprocess.f90
+  postprocess/postprocess.f90
+  postprocess/monitoring.f90
   tdsops.f90
   time_integrator.f90
   ordering.f90

src/case/base_case.f90

@@ -5,7 +5,8 @@ module m_base_case
 
   use m_allocator, only: allocator_t
   use m_base_backend, only: base_backend_t
-  use m_common, only: dp, DIR_X, DIR_Z, DIR_C, VERT
+  use m_common, only: dp, DIR_X, DIR_C, VERT
+  use m_monitoring, only: monitoring_t
   use m_field, only: field_t, flist_t
   use m_mesh, only: mesh_t
   use m_solver, only: solver_t, init
@@ -19,6 +20,7 @@ module m_base_case
   type, abstract :: base_case_t
     class(solver_t), allocatable :: solver
     type(io_manager_t) :: io_mgr
+    type(monitoring_t) :: monitoring
   contains
     procedure(boundary_conditions), deferred :: boundary_conditions
     procedure(initial_conditions), deferred :: initial_conditions
@@ -29,8 +31,6 @@ module m_base_case
     procedure :: case_finalise
     procedure :: set_init
     procedure :: run
-    procedure :: print_enstrophy
-    procedure :: print_div_max_mean
   end type base_case_t
 
   abstract interface
@@ -109,13 +109,16 @@ subroutine case_init(self, backend, mesh, host_allocator)
       call self%initial_conditions()
     end if
 
+    call self%monitoring%init(self%solver)
+
   end subroutine case_init
 
   subroutine case_finalise(self)
     class(base_case_t) :: self
 
     if (self%solver%mesh%par%is_root()) print *, 'run end'
 
+    call self%monitoring%finalise()
     call self%io_mgr%finalise()
   end subroutine case_finalise
 
@@ -158,57 +161,6 @@ end function field_func
 
   end subroutine set_init
 
-  subroutine print_enstrophy(self, u, v, w)
-    !! Reports the enstrophy
-    implicit none
-
-    class(base_case_t), intent(in) :: self
-    class(field_t), intent(in) :: u, v, w
-
-    class(field_t), pointer :: du, dv, dw
-    real(dp) :: enstrophy
-
-    du => self%solver%backend%allocator%get_block(DIR_X, VERT)
-    dv => self%solver%backend%allocator%get_block(DIR_X, VERT)
-    dw => self%solver%backend%allocator%get_block(DIR_X, VERT)
-
-    call self%solver%curl(du, dv, dw, u, v, w)
-
-    enstrophy = 0.5_dp*(self%solver%backend%scalar_product(du, du) &
-                        + self%solver%backend%scalar_product(dv, dv) &
-                        + self%solver%backend%scalar_product(dw, dw)) &
-                /self%solver%ngrid
-
-    if (self%solver%mesh%par%is_root()) print *, 'enstrophy:', enstrophy
-
-    call self%solver%backend%allocator%release_block(du)
-    call self%solver%backend%allocator%release_block(dv)
-    call self%solver%backend%allocator%release_block(dw)
-
-  end subroutine print_enstrophy
-
-  subroutine print_div_max_mean(self, u, v, w)
-    !! Reports the div(u) at cell centres
-    implicit none
-
-    class(base_case_t), intent(in) :: self
-    class(field_t), intent(in) :: u, v, w
-
-    class(field_t), pointer :: div_u
-    real(dp) :: div_u_max, div_u_mean
-
-    div_u => self%solver%backend%allocator%get_block(DIR_Z)
-
-    call self%solver%divergence_v2p(div_u, u, v, w)
-
-    call self%solver%backend%field_max_mean(div_u_max, div_u_mean, div_u)
-    if (self%solver%mesh%par%is_root()) &
-      print *, 'div u max mean:', div_u_max, div_u_mean
-
-    call self%solver%backend%allocator%release_block(div_u)
-
-  end subroutine print_div_max_mean
-
   subroutine run(self)
     !! Runs the solver forwards in time from t=t_0 to t=T, performing
     !! postprocessing/IO and reporting diagnostics.

src/case/channel.f90

@@ -158,8 +158,8 @@ subroutine postprocess_channel(self, iter, t)
       print *, 'time =', t, 'iteration =', iter
     end if
 
-    call self%print_enstrophy(self%solver%u, self%solver%v, self%solver%w)
-    call self%print_div_max_mean(self%solver%u, self%solver%v, self%solver%w)
+    call self%monitoring%write_step( &
+      self%solver, t, self%solver%u, self%solver%v, self%solver%w)
 
   end subroutine postprocess_channel
 

src/case/cylinder.f90

@@ -123,8 +123,8 @@ subroutine postprocess_cylinder(self, iter, t)
       print *, 'time =', t, 'iteration =', iter
     end if
 
-    call self%print_enstrophy(self%solver%u, self%solver%v, self%solver%w)
-    call self%print_div_max_mean(self%solver%u, self%solver%v, self%solver%w)
+    call self%monitoring%write_step( &
+      self%solver, t, self%solver%u, self%solver%v, self%solver%w)
 
   end subroutine postprocess_cylinder
 

src/case/generic.f90

@@ -90,8 +90,8 @@ subroutine postprocess_generic(self, iter, t)
       print *, 'time =', t, 'iteration =', iter
     end if
 
-    call self%print_enstrophy(self%solver%u, self%solver%v, self%solver%w)
-    call self%print_div_max_mean(self%solver%u, self%solver%v, self%solver%w)
+    call self%monitoring%write_step( &
+      self%solver, t, self%solver%u, self%solver%v, self%solver%w)
 
   end subroutine postprocess_generic
 

src/case/tgv.f90

@@ -109,8 +109,8 @@ subroutine postprocess_tgv(self, iter, t)
       print *, 'time =', t, 'iteration =', iter
     end if
 
-    call self%print_enstrophy(self%solver%u, self%solver%v, self%solver%w)
-    call self%print_div_max_mean(self%solver%u, self%solver%v, self%solver%w)
+    call self%monitoring%write_step( &
+      self%solver, t, self%solver%u, self%solver%v, self%solver%w)
 
   end subroutine postprocess_tgv
 

src/postprocess/monitoring.f90

@@ -0,0 +1,134 @@
+module m_monitoring
+  !! Computes and logs global scalar monitoring quantities to a file.
+  !!
+  !! Tracks the following quantities:
+  !!
+  !! - Kinetic energy: \( E_k = \frac{1}{2N} \sum (u^2 + v^2 + w^2) \)
+  !! - Enstrophy: \( \mathcal{E} = \frac{1}{2N} \sum |\nabla \times \mathbf{u}|^2 \)
+  !! - Dissipation rate: \( \varepsilon = -\frac{\nu}{N} \sum \mathbf{u} \cdot \nabla^2 \mathbf{u} \)
+  !! - Divergence: \( \max |\nabla \cdot \mathbf{u}| \) and mean (divergence-free check)
+
+  use m_common, only: dp, DIR_X, DIR_Z, VERT
+  use m_field, only: field_t
+  use m_solver, only: solver_t
+
+  implicit none
+
+  type :: monitoring_t
+    integer, private :: file_unit = -1
+    logical, private :: is_root = .false.
+  contains
+    procedure :: init
+    procedure :: write_step
+    procedure :: finalise
+  end type monitoring_t
+
+contains
+
+  subroutine init(self, solver)
+    class(monitoring_t), intent(inout) :: self
+    class(solver_t), intent(in) :: solver
+
+    self%is_root = solver%mesh%par%is_root()
+
+    if (self%is_root) then
+      open (newunit=self%file_unit, file='monitoring.csv', &
+            status='replace', action='write')
+      write (self%file_unit, '(A)') &
+        '# time, enstrophy, ke, eps, div_u_max, div_u_mean'
+    end if
+
+  end subroutine init
+
+  subroutine write_step(self, solver, t, u, v, w)
+    !! Computes and reports monitoring quantities, writing to both stdout
+    !! and the monitoring file.
+    class(monitoring_t), intent(inout) :: self
+    class(solver_t), intent(inout) :: solver
+    real(dp), intent(in) :: t
+    class(field_t), intent(in) :: u, v, w
+
+    class(field_t), pointer :: du, dv, dw, div_u, lapl
+    real(dp) :: enstrophy, ke, eps, div_u_max, div_u_mean
+
+    !! Kinetic energy: \( E_k = \frac{1}{2N} \sum (u^2 + v^2 + w^2) \)
+    ke = 0.5_dp*(solver%backend%scalar_product(u, u) &
+                 + solver%backend%scalar_product(v, v) &
+                 + solver%backend%scalar_product(w, w)) &
+         /solver%ngrid
+
+    !! Enstrophy: \( \mathcal{E} = \frac{1}{2N} \sum |\nabla \times \mathbf{u}|^2 \)
+    du => solver%backend%allocator%get_block(DIR_X, VERT)
+    dv => solver%backend%allocator%get_block(DIR_X, VERT)
+    dw => solver%backend%allocator%get_block(DIR_X, VERT)
+
+    call solver%curl(du, dv, dw, u, v, w)
+
+    enstrophy = 0.5_dp*(solver%backend%scalar_product(du, du) &
+                        + solver%backend%scalar_product(dv, dv) &
+                        + solver%backend%scalar_product(dw, dw)) &
+                /solver%ngrid
+
+    call solver%backend%allocator%release_block(du)
+    call solver%backend%allocator%release_block(dv)
+    call solver%backend%allocator%release_block(dw)
+
+    !! Dissipation rate (2nd derivative form):
+    !! \( \varepsilon = -\frac{\nu}{N} \sum \mathbf{u} \cdot \nabla^2 \mathbf{u} \)
+    !! Equivalent to the strain-rate form for periodic flows.
+    lapl => solver%backend%allocator%get_block(DIR_X, VERT)
+
+    call solver%vector_calculus%laplacian( &
+      lapl, u, solver%xdirps%der2nd, solver%ydirps%der2nd, &
+      solver%zdirps%der2nd &
+      )
+    eps = solver%backend%scalar_product(u, lapl)
+
+    call solver%vector_calculus%laplacian( &
+      lapl, v, solver%xdirps%der2nd, solver%ydirps%der2nd, &
+      solver%zdirps%der2nd &
+      )
+    eps = eps + solver%backend%scalar_product(v, lapl)
+
+    call solver%vector_calculus%laplacian( &
+      lapl, w, solver%xdirps%der2nd, solver%ydirps%der2nd, &
+      solver%zdirps%der2nd &
+      )
+    eps = eps + solver%backend%scalar_product(w, lapl)
+
+    call solver%backend%allocator%release_block(lapl)
+
+    eps = -solver%nu*eps/solver%ngrid
+
+    !! Divergence: \( \nabla \cdot \mathbf{u} \) max and mean (should be ~0)
+    div_u => solver%backend%allocator%get_block(DIR_Z)
+
+    call solver%divergence_v2p(div_u, u, v, w)
+
+    call solver%backend%field_max_mean(div_u_max, div_u_mean, div_u)
+
+    call solver%backend%allocator%release_block(div_u)
+
+    ! Print to stdout and write to file (root only)
+    if (self%is_root) then
+      print *, 'enstrophy:', enstrophy
+      print *, 'div u max mean:', div_u_max, div_u_mean
+
+      write (self%file_unit, '(ES20.12,5(",",ES20.12))') &
+        t, enstrophy, ke, eps, div_u_max, div_u_mean
+      flush (self%file_unit)
+    end if
+
+  end subroutine write_step
+
+  subroutine finalise(self)
+    class(monitoring_t), intent(inout) :: self
+
+    if (self%is_root .and. self%file_unit /= -1) then
+      close (self%file_unit)
+      self%file_unit = -1
+    end if
+
+  end subroutine finalise
+
+end module m_monitoring