From 5265ca4f0c9f13c437b0c7d56b6385701ae47a28 Mon Sep 17 00:00:00 2001 From: Irufan Ahmed Date: Sat, 28 Mar 2026 19:23:34 +0000 Subject: [PATCH 1/2] output monitors to csv file --- src/CMakeLists.txt | 3 +- src/case/base_case.f90 | 60 ++---------- src/case/channel.f90 | 4 +- src/case/cylinder.f90 | 4 +- src/case/generic.f90 | 4 +- src/case/tgv.f90 | 4 +- src/postprocess/monitoring.f90 | 134 ++++++++++++++++++++++++++ src/{ => postprocess}/postprocess.f90 | 19 ++-- 8 files changed, 160 insertions(+), 72 deletions(-) create mode 100644 src/postprocess/monitoring.f90 rename src/{ => postprocess}/postprocess.f90 (93%) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 028e9d9b6..3d6f4b616 100644 --- a/src/CMakeLists.txt +++ b/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 diff --git a/src/case/base_case.f90 b/src/case/base_case.f90 index 8464db500..f12e9e6a3 100644 --- a/src/case/base_case.f90 +++ b/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,6 +109,8 @@ 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) @@ -116,6 +118,7 @@ subroutine case_finalise(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. diff --git a/src/case/channel.f90 b/src/case/channel.f90 index 846dcf93e..8d4b876b1 100644 --- a/src/case/channel.f90 +++ b/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 diff --git a/src/case/cylinder.f90 b/src/case/cylinder.f90 index a15130e26..9a6afc52d 100644 --- a/src/case/cylinder.f90 +++ b/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 diff --git a/src/case/generic.f90 b/src/case/generic.f90 index 193767999..66927aead 100644 --- a/src/case/generic.f90 +++ b/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 diff --git a/src/case/tgv.f90 b/src/case/tgv.f90 index 971ee29fb..dfb4b9e83 100644 --- a/src/case/tgv.f90 +++ b/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 diff --git a/src/postprocess/monitoring.f90 b/src/postprocess/monitoring.f90 new file mode 100644 index 000000000..a47ab3774 --- /dev/null +++ b/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 diff --git a/src/postprocess.f90 b/src/postprocess/postprocess.f90 similarity index 93% rename from src/postprocess.f90 rename to src/postprocess/postprocess.f90 index fe0a54eb2..719f3cf42 100644 --- a/src/postprocess.f90 +++ b/src/postprocess/postprocess.f90 @@ -39,14 +39,13 @@ subroutine compute_derived_fields(solver, output_vorticity, & v_y, dvdy_y, v_z, dvdz_z, & w_y, dwdy_y, w_z, dwdz_z - ! Lazily allocate output fields + ! Allocate output fields on first use if (output_vorticity .and. .not. associated(solver%vort)) & solver%vort => solver%backend%allocator%get_block(DIR_X, VERT) if (output_qcriterion .and. .not. associated(solver%qcrit)) & solver%qcrit => solver%backend%allocator%get_block(DIR_X, VERT) - ! --- Compute all 9 velocity gradient components --- - + ! Compute all 9 velocity gradient components ! du/dx (x-derivative in x-pencil, no reorder needed) dudx => solver%backend%allocator%get_block(DIR_X) call solver%backend%tds_solve(dudx, solver%u, & @@ -133,8 +132,10 @@ subroutine compute_derived_fields(solver, output_vorticity, & ! dvdx, dvdy_x, dvdz_x ! dwdx, dwdy_x, dwdz_x - ! Vorticity magnitude: - ! |w| = sqrt((dw/dy-dv/dz)^2+(du/dz-dw/dx)^2+(dv/dx-du/dy)^2) + !! Vorticity magnitude: + !! \( |\boldsymbol{\omega}| = \sqrt{(\partial w/\partial y - \partial v/\partial z)^2 + !! + (\partial u/\partial z - \partial w/\partial x)^2 + !! + (\partial v/\partial x - \partial u/\partial y)^2} \) if (output_vorticity) then solver%vort%data = sqrt( & (dwdy_x%data - dvdz_x%data)**2 + & @@ -143,9 +144,9 @@ subroutine compute_derived_fields(solver, output_vorticity, & ) end if - ! Q-criterion: - ! Q = -0.5*(dudx^2+dvdy^2+dwdz^2) - ! - dudy*dvdx - dudz*dwdx - dvdz*dwdy + !! Q-criterion: + !! \( Q = -\frac{1}{2}(u_{x}^2 + v_{y}^2 + w_{z}^2) + !! - u_{y}v_{x} - u_{z}w_{x} - v_{z}w_{y} \) if (output_qcriterion) then solver%qcrit%data = & -0.5_dp*(dudx%data**2 & @@ -181,7 +182,7 @@ subroutine compute_pressure_vert(solver) error stop 'compute_pressure_vert: pressure not yet computed' end if - ! Lazy allocate pressure_vert on first call + ! Allocate pressure_vert on first use if (.not. associated(solver%pressure_vert)) then solver%pressure_vert => & solver%backend%allocator%get_block(DIR_X, VERT) From 0b6dd17665be07f7eb9753be0d55e1d77beb1c1c Mon Sep 17 00:00:00 2001 From: Irufan Ahmed Date: Sat, 28 Mar 2026 23:27:08 +0000 Subject: [PATCH 2/2] docs: update docs to include run monitoring --- docs/source/user/post-processing/index.rst | 1 + .../user/post-processing/monitoring.rst | 93 +++++++++++++++++++ 2 files changed, 94 insertions(+) create mode 100644 docs/source/user/post-processing/monitoring.rst diff --git a/docs/source/user/post-processing/index.rst b/docs/source/user/post-processing/index.rst index 0f650de71..9dad19a51 100644 --- a/docs/source/user/post-processing/index.rst +++ b/docs/source/user/post-processing/index.rst @@ -7,4 +7,5 @@ simulation output. .. toctree:: :maxdepth: 2 + monitoring.rst paraview.rst \ No newline at end of file diff --git a/docs/source/user/post-processing/monitoring.rst b/docs/source/user/post-processing/monitoring.rst new file mode 100644 index 000000000..006fd7a1b --- /dev/null +++ b/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")