Add FieldPoyntingFlux reduced diagnostic (#5475)

This adds a reduced diagnostic that calculates the Poynting flux on the
surfaces of the domain, providing the power flow into and out of the
domain. This also includes the time integrated data.

When using the implicit evolve scheme, to get the energy accounting
correct, the flux needs to be calculated at the mid step. For this
reason, the `ComputeDiagsMidStep` was added which is called directly at
the appropriate times.

Because of the time integration, there are two main differences of this
reduced diagnostic compared to the others. The first is that it is
calculated every time step in order to get the full resolution in time.
The intervals parameter still controls how often the diagnostic data is
written out. The second is that a facility is added to write out the
values of the time integration to a file when a checkpoint is made, so
on a restart the integration can continue with the previous values. The
facility was written in a general way so that other reduced diagnostics
can also do this.

The CI test using the implicit solver is dependent on PR #5498 and PR
#5489.

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

Docs/source/usage/parameters.rst

@@ -3182,6 +3182,12 @@ This shifts analysis from post-processing to runtime calculation of reduction op
         Note that the fields are averaged on the cell centers before their maximum values are
         computed.
 
+    * ``FieldPoyntingFlux``
+        Integrates the normal Poynting flux over each domain boundary surface and also integrates the flux over time.
+        This provides the power and total energy loss into or out of the simulation domain.
+        The output columns are the flux for each dimension on the lower boundaries, then the higher boundaries,
+        then the integrated energy loss for each dimension on the the lower and higher boundaries.
+
     * ``FieldProbe``
         This type computes the value of each component of the electric and magnetic fields
         and of the Poynting vector (a measure of electromagnetic flux) at points in the domain.

Examples/Tests/pec/CMakeLists.txt

@@ -40,3 +40,23 @@ add_warpx_test(
     "analysis_default_regression.py --path diags/diag1000010"  # checksum
     OFF  # dependency
 )
+
+add_warpx_test(
+    test_2d_pec_field_insulator_implicit  # name
+    2  # dims
+    2  # nprocs
+    inputs_test_2d_pec_field_insulator_implicit  # inputs
+    "analysis_pec_insulator_implicit.py diags/diag1000020"  # analysis
+    "analysis_default_regression.py --path diags/diag1000020"  # checksum
+    OFF  # dependency
+)
+
+add_warpx_test(
+    test_2d_pec_field_insulator_implicit_restart  # name
+    2  # dims
+    2  # nprocs
+    inputs_test_2d_pec_field_insulator_implicit_restart  # inputs
+    "analysis_pec_insulator_implicit.py diags/diag1000020"  # analysis
+    "analysis_default_regression.py --path diags/diag1000020"  # checksum
+    test_2d_pec_field_insulator_implicit  # dependency
+)

Examples/Tests/pec/analysis_pec_insulator_implicit.py

@@ -0,0 +1,57 @@
+#!/usr/bin/env python3
+
+#
+#
+# This file is part of WarpX.
+#
+# License: BSD-3-Clause-LBNL
+#
+# This is a script that analyses the simulation results from
+# the scripts `inputs_test_2d_pec_field_insulator_implicit` and
+# `inputs_test_2d_pec_field_insulator_implicit_restart`.
+# The scripts model an insulator boundary condition on part of the
+# upper x boundary that pushes B field into the domain. The implicit
+# solver is used, converging to machine tolerance. The energy accounting
+# should be exact to machine precision, so that the energy is the system
+# should be the same as the amount of energy pushed in from the boundary.
+# This is checked using the FieldEnergy and FieldPoyntingFlux reduced
+# diagnostics.
+import sys
+
+import matplotlib
+
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import numpy as np
+
+# this will be the name of the plot file
+fn = sys.argv[1]
+
+EE = np.loadtxt(f"{fn}/../reducedfiles/fieldenergy.txt", skiprows=1)
+SS = np.loadtxt(f"{fn}/../reducedfiles/poyntingflux.txt", skiprows=1)
+SSsum = SS[:, 2:6].sum(1)
+EEloss = SS[:, 7:].sum(1)
+
+dt = EE[1, 1]
+
+fig, ax = plt.subplots()
+ax.plot(EE[:, 0], EE[:, 2], label="field energy")
+ax.plot(SS[:, 0], -EEloss, label="-flux*dt")
+ax.legend()
+ax.set_xlabel("time (s)")
+ax.set_ylabel("energy (J)")
+fig.savefig("energy_history.png")
+
+fig, ax = plt.subplots()
+ax.plot(EE[:, 0], (EE[:, 2] + EEloss) / EE[:, 2].max())
+ax.set_xlabel("time (s)")
+ax.set_ylabel("energy difference/max energy (1)")
+fig.savefig("energy_difference.png")
+
+tolerance_rel = 1.0e-13
+
+energy_difference_fraction = np.abs((EE[:, 2] + EEloss) / EE[:, 2].max()).max()
+print(f"energy accounting error = {energy_difference_fraction}")
+print(f"tolerance_rel = {tolerance_rel}")
+
+assert energy_difference_fraction < tolerance_rel

Examples/Tests/pec/inputs_test_2d_pec_field_insulator_implicit

@@ -0,0 +1,73 @@
+# Maximum number of time steps
+max_step = 20
+
+# number of grid points
+amr.n_cell = 32 32
+amr.blocking_factor = 16
+
+# Maximum level in hierarchy (for now must be 0, i.e., one level in total)
+amr.max_level = 0
+
+# Geometry
+geometry.dims = 2
+geometry.prob_lo = 0.     2.e-2 # physical domain
+geometry.prob_hi = 1.e-2  3.e-2
+
+# Boundary condition
+boundary.field_lo = neumann      periodic
+boundary.field_hi = pec_insulator periodic
+
+insulator.area_x_hi(y,z) = (2.25e-2 <= z and z <= 2.75e-2)
+insulator.By_x_hi(y,z,t) = min(t/1.0e-12,1)*1.e1*3.3e-4
+
+warpx.serialize_initial_conditions = 1
+
+# Implicit setup
+# Note that this is the CFL step size for the explicit simulation, over 2.
+# This value allows quick convergence of the Picard solver.
+warpx.const_dt = 7.37079480234276e-13/2.
+
+algo.maxwell_solver = Yee
+algo.evolve_scheme = "theta_implicit_em"
+#algo.evolve_scheme = "semi_implicit_em"
+
+implicit_evolve.theta = 0.5
+#implicit_evolve.max_particle_iterations = 21
+#implicit_evolve.particle_tolerance = 1.0e-12
+
+implicit_evolve.nonlinear_solver = "picard"
+picard.verbose = true
+picard.max_iterations = 25
+picard.relative_tolerance = 0.0
+picard.absolute_tolerance = 0.0
+picard.require_convergence = false
+
+#implicit_evolve.nonlinear_solver = "newton"
+#newton.verbose = true
+#newton.max_iterations = 20
+#newton.relative_tolerance = 1.0e-20
+#newton.absolute_tolerance = 0.0
+#newton.require_convergence = false
+
+#gmres.verbose_int = 2
+#gmres.max_iterations = 1000
+#gmres.relative_tolerance = 1.0e-20
+#gmres.absolute_tolerance = 0.0
+
+# Verbosity
+warpx.verbose = 1
+
+# Diagnostics
+diagnostics.diags_names = diag1 chk
+diag1.intervals = 20
+diag1.diag_type = Full
+
+chk.intervals = 10
+chk.diag_type = Full
+chk.format = checkpoint
+
+warpx.reduced_diags_names = fieldenergy poyntingflux
+poyntingflux.type = FieldPoyntingFlux
+poyntingflux.intervals = 1
+fieldenergy.type = FieldEnergy
+fieldenergy.intervals = 1

Examples/Tests/pec/inputs_test_2d_pec_field_insulator_implicit_restart

@@ -0,0 +1,5 @@
+# base input parameters
+FILE = inputs_test_2d_pec_field_insulator_implicit
+
+# test input parameters
+amr.restart = "../test_2d_pec_field_insulator_implicit/diags/chk000010"

Examples/Tests/reduced_diags/inputs_test_3d_reduced_diags

@@ -68,7 +68,7 @@ photons.uz_th = 0.2
 #################################
 ###### REDUCED DIAGS ############
 #################################
-warpx.reduced_diags_names = EP NP EF PP PF MF MR FP FP_integrate FP_line FP_plane FR_Max FR_Min FR_Integral Edotj
+warpx.reduced_diags_names = EP NP EF PP PF MF PX MR FP FP_integrate FP_line FP_plane FR_Max FR_Min FR_Integral Edotj
 EP.type = ParticleEnergy
 EP.intervals = 200
 EF.type = FieldEnergy
@@ -79,6 +79,8 @@ PF.type = FieldMomentum
 PF.intervals = 200
 MF.type = FieldMaximum
 MF.intervals = 200
+PX.type = FieldPoyntingFlux
+PX.intervals = 200
 FP.type = FieldProbe
 FP.intervals = 200
 #The probe is placed at a cell center to match the value in the plotfile

Python/pywarpx/picmi.py

@@ -4074,6 +4074,7 @@ def __init__(
             "FieldEnergy",
             "FieldMomentum",
             "FieldMaximum",
+            "FieldPoyntingFlux",
             "RhoMaximum",
             "ParticleNumber",
             "LoadBalanceCosts",

Regression/Checksum/benchmarks_json/test_2d_pec_field_insulator_implicit.json

@@ -0,0 +1,14 @@
+{
+  "lev=0": {
+    "Bx": 0.0,
+    "By": 0.35907571934346943,
+    "Bz": 0.0,
+    "Ex": 36840284.366667606,
+    "Ey": 0.0,
+    "Ez": 107777138.0847348,
+    "jx": 0.0,
+    "jy": 0.0,
+    "jz": 0.0
+  }
+}
+

Regression/Checksum/benchmarks_json/test_2d_pec_field_insulator_implicit_restart.json

@@ -0,0 +1,14 @@
+{
+  "lev=0": {
+    "Bx": 0.0,
+    "By": 0.35907571934346943,
+    "Bz": 0.0,
+    "Ex": 36840284.366667606,
+    "Ey": 0.0,
+    "Ez": 107777138.0847348,
+    "jx": 0.0,
+    "jy": 0.0,
+    "jz": 0.0
+  }
+}
+

Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.H

@@ -35,6 +35,8 @@ class FlushFormatCheckpoint final : public FlushFormatPlotfile
                               const amrex::Vector<ParticleDiag>& particle_diags) const;
 
     void WriteDMaps (const std::string& dir, int nlev) const;
+
+    void WriteReducedDiagsData (std::string const & dir) const;
 };
 
 #endif // WARPX_FLUSHFORMATCHECKPOINT_H_

Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp

@@ -5,6 +5,7 @@
 #   include "BoundaryConditions/PML_RZ.H"
 #endif
 #include "Diagnostics/ParticleDiag/ParticleDiag.H"
+#include "Diagnostics/ReducedDiags/MultiReducedDiags.H"
 #include "Fields.H"
 #include "Particles/WarpXParticleContainer.H"
 #include "Utils/TextMsg.H"
@@ -174,6 +175,8 @@ FlushFormatCheckpoint::WriteToFile (
 
     WriteDMaps(checkpointname, nlev);
 
+    WriteReducedDiagsData(checkpointname);
+
     VisMF::SetHeaderVersion(current_version);
 
 }
@@ -263,3 +266,12 @@ FlushFormatCheckpoint::WriteDMaps (const std::string& dir, int nlev) const
         }
     }
 }
+
+void
+FlushFormatCheckpoint::WriteReducedDiagsData (std::string const & dir) const
+{
+    if (ParallelDescriptor::IOProcessor()) {
+        auto & warpx = WarpX::GetInstance();
+        warpx.reduced_diags->WriteCheckpointData(dir);
+    }
+}

Source/Diagnostics/ReducedDiags/CMakeLists.txt

@@ -9,6 +9,7 @@ foreach(D IN LISTS WarpX_DIMS)
         FieldEnergy.cpp
         FieldMaximum.cpp
         FieldMomentum.cpp
+        FieldPoyntingFlux.cpp
         FieldProbe.cpp
         FieldProbeParticleContainer.cpp
         FieldReduction.cpp

Source/Diagnostics/ReducedDiags/FieldPoyntingFlux.H

@@ -0,0 +1,63 @@
+/* Copyright 2019-2020
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#ifndef WARPX_DIAGNOSTICS_REDUCEDDIAGS_FIELDPOYTINGFLUX_H_
+#define WARPX_DIAGNOSTICS_REDUCEDDIAGS_FIELDPOYTINGFLUX_H_
+
+#include "ReducedDiags.H"
+
+#include <string>
+
+/**
+ * \brief This class mainly contains a function that computes the field Poynting flux,
+ *        S = E cross B, integrated over each face of the domain.
+ */
+class FieldPoyntingFlux : public ReducedDiags
+{
+public:
+
+    /**
+     * \brief Constructor
+     *
+     * \param[in] rd_name reduced diags names
+     */
+    FieldPoyntingFlux (const std::string& rd_name);
+
+    /**
+     * \brief Call the routine to compute the Poynting flux if needed
+     *
+     * \param[in] step current time step
+     */
+    void ComputeDiags (int step) final;
+
+    /**
+     * \brief Call the routine to compute the Poynting flux at the mid step time level
+     *
+     * \param[in] step current time step
+     */
+    void ComputeDiagsMidStep (int step) final;
+
+    /**
+     * \brief This function computes the electromagnetic Poynting flux,
+     * obtained by integrating the electromagnetic Poynting flux density g = eps0 * (E x B)
+     * on the surface of the domain.
+     *
+     * \param[in] step current time step
+     */
+    void ComputePoyntingFlux ();
+
+    void WriteCheckpointData (std::string const & dir) final;
+
+    void ReadCheckpointData (std::string const & dir) final;
+
+private:
+
+    bool use_mid_step_value = false;
+
+};
+
+#endif