Add FieldPoyntingFlux reduced diagnostic

Docs/source/usage/parameters.rst

@@ -3150,6 +3150,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(
     diags/diag1000010  # output
     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  # analysis
+    diags/diag1000020  # output
+    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  # analysis
+    diags/diag1000020  # output
+    OFF  # dependency
+)

Examples/Tests/pec/analysis_pec_insulator_implicit.py

@@ -0,0 +1,67 @@
+#!/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 os
+import sys
+
+import matplotlib
+
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import numpy as np
+
+sys.path.insert(1, "../../../../warpx/Regression/Checksum/")
+from checksumAPI import evaluate_checksum
+
+# 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
+
+# compare checksums
+evaluate_checksum(
+    test_name=os.path.split(os.getcwd())[1],
+    output_file=sys.argv[1],
+)

Examples/Tests/pec/inputs_test_2d_pec_field_insulator_implicit

@@ -0,0 +1,71 @@
+# 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
+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

@@ -3996,6 +3996,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

@@ -34,6 +34,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"
@@ -171,6 +172,8 @@ FlushFormatCheckpoint::WriteToFile (
 
     WriteDMaps(checkpointname, nlev);
 
+    WriteReducedDiagsData(checkpointname);
+
     VisMF::SetHeaderVersion(current_version);
 
 }
@@ -242,3 +245,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,62 @@
+/* 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.
+ */
+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