Merge pull request #233 from awslabs/sjg/error-estimator-dev

Upgrade error estimation for time-dependent problems

CHANGELOG.md

@@ -37,6 +37,9 @@ The format of this changelog is based on
     disable the associated boundary condition and only use the surface for postprocessing.
   - Changed the smooth flux space for the electrostatic error estimator to fix performance
     on problems with material interfaces.
+  - Fixed error estimation bug affecting time-dependent simulation types (driven, transient,
+    eigenmode) where the recovery of the electric flux density also needs to be taken into
+    account in addition to the magnetic field.
   - Fixed a bug related to mesh cleaning for unspecified domains and mesh partitioning.
   - Change computation of domain energy postprocessing for electrostatic and magnetostatic
     simulation types in order to improve performance.

docs/src/reference.md

@@ -10,7 +10,7 @@
 The solver computes a finite element approximation to the three-dimensional, time-harmonic
 Maxwell's equations in second-order form. The nondimensionalized, source-free, boundary
 value problem for ``\bm{E}(\bm{x})\in\mathbb{C}^3``, ``\bm{x}\in\Omega``,
-``\partial\Omega=\Gamma``, where
+``\partial\Omega = \Gamma``, where
 ``\bm{\mathscr{E}}(\bm{x},t) = \text{Re}\{\bm{E}(\bm{x})e^{i\omega t}\}`` denotes the
 electric field, is written as
 
@@ -53,21 +53,21 @@ complex-valued quantity:
 ```
 
 where ``\varepsilon_r'`` is the real relative permittivity and ``\tan{\delta}`` is the loss
-tangent. Alternatively, conductor loss is modeled by Ohm's law ``\bm{J}=\sigma\bm{E}`` with
-electrical conductivity ``\sigma>0.0``. For a superconducting domain, the constitive
+tangent. Alternatively, conductor loss is modeled by Ohm's law ``\bm{J} = \sigma\bm{E}``
+with electrical conductivity ``\sigma > 0.0``. For a superconducting domain, the constitive
 current-field relationship given by Ohm's law is replaced by that given by the London
 equations:
 
 ```math
-\frac{\partial \bm{J}}{\partial t}=\frac{1}{\mu_r\lambda_L^2}\bm{E}
+\frac{\partial \bm{J}}{\partial t} = \frac{1}{\mu_r\lambda_L^2}\bm{E}
 ```
 
 where ``\lambda_L = \sqrt{m/\mu n_s e^2}/L_0`` is the nondimensionalized London penetration
 depth. In this case, the term ``+i\omega\sigma \bm{E}`` arising for a normal conductor in
 the time-harmonic Maxwell's equations becomes ``+(\mu_r \lambda_L^2)^{-1}\bm{E}``.
 
-The domain boundary ``\Gamma=\Gamma_{PEC}\cup\Gamma_{PMC}\cup\Gamma_{Z}``, is separated into
-perfect electric conductor (PEC), perfect magnetic conductor (PMC), and impedance
+The domain boundary ``\Gamma = \Gamma_{PEC}\cup\Gamma_{PMC}\cup\Gamma_{Z}``, is separated
+into perfect electric conductor (PEC), perfect magnetic conductor (PMC), and impedance
 boundaries, respectively. The PEC boundary condition is a homogeneous Dirichlet condition,
 while the PMC boundary condition is the natural boundary condition for the problem and is
 satisfied at all exterior boundaries by the finite element formulation. Impedance
@@ -201,12 +201,12 @@ In the time domain, the time histories of the port voltages can be Fourier-trans
 get their frequency domain representation for scattering parameter calculation.
 
 Numeric wave ports assume a field with known normal-direction dependence
-``\bm{E}(\bm{x})=\bm{e}(\bm{x}_t)e^{ik_n x_n}`` where ``k_n`` is the propagation constant.
+``\bm{E}(\bm{x}) = \bm{e}(\bm{x}_t)e^{ik_n x_n}`` where ``k_n`` is the propagation constant.
 For each operating frequency ``\omega``, a two-dimensional eigenvalue problem is solved on
 the port yielding the mode shapes ``\bm{e}_m`` and associated propagation constants
 ``k_{n,m}``. These are used in the full 3D model where the Robin port boundary condition has
-coefficient ``\gamma=i\text{Re}\{k_{n,m}\}/\mu_r`` and the computed mode is used to compute
-the incident field in the source term ``\bm{U}^{inc}`` at excited ports. Scattering
+coefficient ``\gamma = i\text{Re}\{k_{n,m}\}/\mu_r`` and the computed mode is used to
+compute the incident field in the source term ``\bm{U}^{inc}`` at excited ports. Scattering
 parameter postprocessing takes the same form as the lumped port counterpart using the
 computed modal solutions. Since the propagation constants are known for each wave port,
 scattering parameter de-embedding can be performed by specifying an offset distance ``d``
@@ -229,7 +229,7 @@ condition, is a special case of the general impedance boundary condition describ
 ```
 
 This is also known as the Sommerfeld radiation condition, and one can recognize the
-dependence on the impedance of free space ``Z_0^{-1}=\sqrt{\mu_r^{-1}\varepsilon_r}``. The
+dependence on the impedance of free space ``Z_0^{-1} = \sqrt{\mu_r^{-1}\varepsilon_r}``. The
 second-order absorbing boundary condition is
 
 ```math
@@ -238,7 +238,7 @@ second-order absorbing boundary condition is
     - \beta\nabla\times[(\nabla\times\bm{E})_n\bm{n}] = 0
 ```
 
-where assuming an infinite radius of curvature ``\beta=\mu_r^{-1}c_0/(2i\omega)``, and the
+where assuming an infinite radius of curvature ``\beta = \mu_r^{-1}c_0/(2i\omega)``, and the
 contribution depending on ``(\nabla\cdot\bm{E}_t)`` has been neglected.
 
 Additionally, while metals with finite conductivity can be modeled using an impedance
@@ -249,7 +249,7 @@ account the frequency dependence of the skin depth is
 Z_s = \frac{1+i}{\delta\sigma}
 ```
 
-where ``\delta=\sqrt{2/\mu_r\sigma\omega}`` is the skin depth and ``\sigma`` is the
+where ``\delta = \sqrt{2/\mu_r\sigma\omega}`` is the skin depth and ``\sigma`` is the
 conductivity of the metal. Another model, which takes into account finite thickness effects,
 is given by
 
@@ -258,7 +258,7 @@ Z_s = \frac{1}{\delta\sigma}\left(\frac{\sinh{\nu}+\sin{\nu}}{\cosh{\nu}+\cos{\n
     + i\frac{\sinh{\nu}-\sin{\nu}}{\cosh{\nu}+\cos{\nu}}\right)
 ```
 
-where ``\nu=h/\delta`` and ``h`` is the layer thickness. This model correctly produces the
+where ``\nu = h/\delta`` and ``h`` is the layer thickness. This model correctly produces the
 DC limit when ``h\ll\delta``.
 
 ## Energy-participation ratios
@@ -294,8 +294,8 @@ Finally, the total electric energy in mode ``m`` is
     + \sum_j \frac{1}{2} \, C_jV_{mj}^2
 ```
 
-where ``\bm{D}_m=\varepsilon_r\bm{E}_m`` is the electric flux density for mode ``m`` and the
-second term on the right-hand side accounts for any lumped capacitive boundaries with
+where ``\bm{D}_m = \varepsilon_r\bm{E}_m`` is the electric flux density for mode ``m`` and
+the second term on the right-hand side accounts for any lumped capacitive boundaries with
 nonzero circuit capacitance ``C_j``.
 
 The EPR can also be used to estimate mode quality factors due to input-output(I-O) line
@@ -373,7 +373,7 @@ quality factor for interface ``j`` is given by
 ```
 
 where ``\bm{E}_n`` denotes the normal field to the interface and
-``\bm{E}_t=\bm{E}-\bm{E}_n`` denotes the tangential field.
+``\bm{E}_t = \bm{E}-\bm{E}_n`` denotes the tangential field.
 
 ## Lumped parameter extraction
 
@@ -424,7 +424,7 @@ specified surfaces of the model. The magnetic field energy associated with any s
 \mathcal{E}(\bm{A}_i) = \frac{1}{2}\int_\Omega\mu_r^{-1}\bm{B}_i\cdot\bm{B}_i\,dV
 ```
 
-where ``\bm{B}_i=\nabla\times\bm{A}_i`` is the magnetic flux density. Then, the entries of
+where ``\bm{B}_i = \nabla\times\bm{A}_i`` is the magnetic flux density. Then, the entries of
 the inductance matrix, ``\bm{M}``, are given by
 
 ```math
@@ -435,6 +435,26 @@ the inductance matrix, ``\bm{M}``, are given by
 where ``I_i`` is the excitation current for port ``i``, computed by integrating the source
 surface current ``\bm{J}_s^{inc}`` over the surface of the port.
 
+## Error estimation and adaptive mesh refinement (AMR)
+
+Error estimation is used to provide element-wise error estimates for AMR, as well as a
+global error indicator used to terminate AMR iterations or provide an estimate for solution
+accuracy. A Zienkiewicz–Zhu (ZZ) error estimator based on [[5]](#References) is
+implemented, which measures the error in the recovered magnetic field and electric flux
+density. On element ``K``, we have
+
+```math
+\eta^2_K = \eta_{m,2}^2+\eta_{e,K}^2 =
+    \|\mu_r^{1/2}\bm{R}_{ND}(\mu^{-1}\bm{B})
+    - (\mu_r^{-1/2}\bm{B})\|_{L^2(\Omega_K)}^2
+    + \|\varepsilon_r^{-1/2}\bm{R}_{RT}(\varepsilon_r\bm{E})
+    - (\varepsilon_r^{1/2}\bm{E})\|_{L^2(\Omega_K)}^2
+```
+
+where ``\bm{R}_{ND}`` and ``\bm{R}_{RT}`` are the smooth-space recovery operators which
+orthogonally project their argument onto ``H(\text{curl})`` and ``H(\text{div})``,
+discretized by Nédélec and Raviart-Thomas elements, respectively.
+
 ## References
 
 [1] J.-M. Jin, _The Finite Element Method in Electromagnetics_, Wiley-IEEE Press, Hoboken,
@@ -444,4 +464,6 @@ Oxford, 2003.\
 [3] L. Vardapetyan and L. Demkowicz, Full-wave analysis of dielectric waveguides at a given
 frequency, _Mathematics of Computation_ 72 (2003) 105-129.\
 [4] J. Wenner, R. Barends, R. C. Bialczak, et al., Surface loss of superconducting coplanar
-waveguide resonators, _Applied Physics Letters_ 99, 113513 (2011).
+waveguide resonators, _Applied Physics Letters_ 99, 113513 (2011).\
+[5] S. Nicaise, On Zienkiewicz-Zhu error estimators for Maxwell’s equations, _Comptes Rendus
+Mathematique_ 340 (2005) 697-702.

palace/drivers/drivensolver.cpp

@@ -141,9 +141,10 @@ ErrorIndicator DrivenSolver::SweepUniform(SpaceOperator &spaceop, PostOperator &
   B = 0.0;
 
   // Initialize structures for storing and reducing the results of error estimation.
-  CurlFluxErrorEstimator<ComplexVector> estimator(
-      spaceop.GetMaterialOp(), spaceop.GetNDSpaces(), iodata.solver.linear.estimator_tol,
-      iodata.solver.linear.estimator_max_it, 0, iodata.solver.linear.estimator_mg);
+  TimeDependentFluxErrorEstimator<ComplexVector> estimator(
+      spaceop.GetMaterialOp(), spaceop.GetNDSpaces(), spaceop.GetRTSpaces(),
+      iodata.solver.linear.estimator_tol, iodata.solver.linear.estimator_max_it, 0,
+      iodata.solver.linear.estimator_mg);
   ErrorIndicator indicator;
 
   // Main frequency sweep loop.
@@ -181,8 +182,8 @@ ErrorIndicator DrivenSolver::SweepUniform(SpaceOperator &spaceop, PostOperator &
     postop.SetEGridFunction(E);
     postop.SetBGridFunction(B);
     postop.UpdatePorts(spaceop.GetLumpedPortOp(), spaceop.GetWavePortOp(), omega);
-    double E_elec = postop.GetEFieldEnergy();
-    double E_mag = postop.GetHFieldEnergy();
+    const double E_elec = postop.GetEFieldEnergy();
+    const double E_mag = postop.GetHFieldEnergy();
     Mpi::Print(" Sol. ||E|| = {:.6e} (||RHS|| = {:.6e})\n",
                linalg::Norml2(spaceop.GetComm(), E),
                linalg::Norml2(spaceop.GetComm(), RHS));
@@ -194,7 +195,7 @@ ErrorIndicator DrivenSolver::SweepUniform(SpaceOperator &spaceop, PostOperator &
 
     // Calculate and record the error indicators.
     Mpi::Print(" Updating solution error estimates\n");
-    estimator.AddErrorIndicator(E, indicator);
+    estimator.AddErrorIndicator(E, B, E_elec + E_mag, indicator);
 
     // Postprocess S-parameters and optionally write solution to disk.
     Postprocess(postop, spaceop.GetLumpedPortOp(), spaceop.GetWavePortOp(),
@@ -238,9 +239,10 @@ ErrorIndicator DrivenSolver::SweepAdaptive(SpaceOperator &spaceop, PostOperator
   B = 0.0;
 
   // Initialize structures for storing and reducing the results of error estimation.
-  CurlFluxErrorEstimator<ComplexVector> estimator(
-      spaceop.GetMaterialOp(), spaceop.GetNDSpaces(), iodata.solver.linear.estimator_tol,
-      iodata.solver.linear.estimator_max_it, 0, iodata.solver.linear.estimator_mg);
+  TimeDependentFluxErrorEstimator<ComplexVector> estimator(
+      spaceop.GetMaterialOp(), spaceop.GetNDSpaces(), spaceop.GetRTSpaces(),
+      iodata.solver.linear.estimator_tol, iodata.solver.linear.estimator_max_it, 0,
+      iodata.solver.linear.estimator_mg);
   ErrorIndicator indicator;
 
   // Configure the PROM operator which performs the parameter space sampling and basis
@@ -262,7 +264,18 @@ ErrorIndicator DrivenSolver::SweepAdaptive(SpaceOperator &spaceop, PostOperator
   {
     // Add the HDM solution to the PROM reduced basis.
     promop.UpdatePROM(omega, E);
-    estimator.AddErrorIndicator(E, indicator);
+
+    // Compute B = -1/(iω) ∇ x E on the true dofs, and set the internal GridFunctions in
+    // PostOperator for energy postprocessing and error estimation.
+    BlockTimer bt0(Timer::POSTPRO);
+    Curl.Mult(E.Real(), B.Real());
+    Curl.Mult(E.Imag(), B.Imag());
+    B *= -1.0 / (1i * omega);
+    postop.SetEGridFunction(E, false);
+    postop.SetBGridFunction(B, false);
+    const double E_elec = postop.GetEFieldEnergy();
+    const double E_mag = postop.GetHFieldEnergy();
+    estimator.AddErrorIndicator(E, B, E_elec + E_mag, indicator);
   };
   promop.SolveHDM(omega0, E);
   UpdatePROM(omega0);
@@ -349,8 +362,8 @@ ErrorIndicator DrivenSolver::SweepAdaptive(SpaceOperator &spaceop, PostOperator
     postop.SetEGridFunction(E);
     postop.SetBGridFunction(B);
     postop.UpdatePorts(spaceop.GetLumpedPortOp(), spaceop.GetWavePortOp(), omega);
-    double E_elec = postop.GetEFieldEnergy();
-    double E_mag = postop.GetHFieldEnergy();
+    const double E_elec = postop.GetEFieldEnergy();
+    const double E_mag = postop.GetHFieldEnergy();
     Mpi::Print(" Sol. ||E|| = {:.6e}\n", linalg::Norml2(spaceop.GetComm(), E));
     {
       const double J = iodata.DimensionalizeValue(IoData::ValueType::ENERGY, 1.0);
@@ -397,8 +410,8 @@ void DrivenSolver::Postprocess(const PostOperator &postop,
   // The internal GridFunctions for PostOperator have already been set from the E and B
   // solutions in the main frequency sweep loop.
   const double freq = iodata.DimensionalizeValue(IoData::ValueType::FREQUENCY, omega);
-  double E_cap = postop.GetLumpedCapacitorEnergy(lumped_port_op);
-  double E_ind = postop.GetLumpedInductorEnergy(lumped_port_op);
+  const double E_cap = postop.GetLumpedCapacitorEnergy(lumped_port_op);
+  const double E_ind = postop.GetLumpedInductorEnergy(lumped_port_op);
   PostprocessCurrents(postop, surf_j_op, step, omega);
   PostprocessPorts(postop, lumped_port_op, step, omega);
   if (surf_j_op.Size() == 0)

palace/drivers/eigensolver.cpp

@@ -275,9 +275,10 @@ EigenSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
 
   // Calculate and record the error indicators, and postprocess the results.
   Mpi::Print("\nComputing solution error estimates and performing postprocessing\n");
-  CurlFluxErrorEstimator<ComplexVector> estimator(
-      spaceop.GetMaterialOp(), spaceop.GetNDSpaces(), iodata.solver.linear.estimator_tol,
-      iodata.solver.linear.estimator_max_it, 0, iodata.solver.linear.estimator_mg);
+  TimeDependentFluxErrorEstimator<ComplexVector> estimator(
+      spaceop.GetMaterialOp(), spaceop.GetNDSpaces(), spaceop.GetRTSpaces(),
+      iodata.solver.linear.estimator_tol, iodata.solver.linear.estimator_max_it, 0,
+      iodata.solver.linear.estimator_mg);
   ErrorIndicator indicator;
   if (!KM)
   {
@@ -314,13 +315,13 @@ EigenSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
     postop.SetEGridFunction(E);
     postop.SetBGridFunction(B);
     postop.UpdatePorts(spaceop.GetLumpedPortOp(), omega.real());
-    double E_elec = postop.GetEFieldEnergy();
-    double E_mag = postop.GetHFieldEnergy();
+    const double E_elec = postop.GetEFieldEnergy();
+    const double E_mag = postop.GetHFieldEnergy();
 
     // Calculate and record the error indicators.
     if (i < iodata.solver.eigenmode.n)
     {
-      estimator.AddErrorIndicator(E, indicator);
+      estimator.AddErrorIndicator(E, B, E_elec + E_mag, indicator);
     }
 
     // Postprocess the mode.
@@ -339,8 +340,8 @@ void EigenSolver::Postprocess(const PostOperator &postop,
 {
   // The internal GridFunctions for PostOperator have already been set from the E and B
   // solutions in the main loop over converged eigenvalues.
-  double E_cap = postop.GetLumpedCapacitorEnergy(lumped_port_op);
-  double E_ind = postop.GetLumpedInductorEnergy(lumped_port_op);
+  const double E_cap = postop.GetLumpedCapacitorEnergy(lumped_port_op);
+  const double E_ind = postop.GetLumpedInductorEnergy(lumped_port_op);
   PostprocessEigen(i, omega, error_bkwd, error_abs, num_conv);
   PostprocessPorts(postop, lumped_port_op, i);
   PostprocessEPR(postop, lumped_port_op, i, omega, E_elec + E_cap);

palace/drivers/electrostaticsolver.cpp

@@ -47,7 +47,7 @@ ElectrostaticSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
 
   // Initialize structures for storing and reducing the results of error estimation.
   GradFluxErrorEstimator estimator(
-      laplaceop.GetMaterialOp(), laplaceop.GetH1Space(), laplaceop.GetRTSpaces(),
+      laplaceop.GetMaterialOp(), laplaceop.GetNDSpace(), laplaceop.GetRTSpaces(),
       iodata.solver.linear.estimator_tol, iodata.solver.linear.estimator_max_it, 0,
       iodata.solver.linear.estimator_mg);
   ErrorIndicator indicator;
@@ -75,7 +75,7 @@ ElectrostaticSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
     Grad.AddMult(V[step], E, -1.0);
     postop.SetVGridFunction(V[step]);
     postop.SetEGridFunction(E);
-    double E_elec = postop.GetEFieldEnergy();
+    const double E_elec = postop.GetEFieldEnergy();
     Mpi::Print(" Sol. ||V|| = {:.6e} (||RHS|| = {:.6e})\n",
                linalg::Norml2(laplaceop.GetComm(), V[step]),
                linalg::Norml2(laplaceop.GetComm(), RHS));
@@ -86,7 +86,7 @@ ElectrostaticSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
 
     // Calculate and record the error indicators.
     Mpi::Print(" Updating solution error estimates\n");
-    estimator.AddErrorIndicator(V[step], indicator);
+    estimator.AddErrorIndicator(E, E_elec, indicator);
 
     // Postprocess field solutions and optionally write solution to disk.
     Postprocess(postop, step, idx, E_elec, (step == nstep - 1) ? &indicator : nullptr);

palace/drivers/magnetostaticsolver.cpp

@@ -47,8 +47,8 @@ MagnetostaticSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
   std::vector<double> I_inc(nstep);
 
   // Initialize structures for storing and reducing the results of error estimation.
-  CurlFluxErrorEstimator<Vector> estimator(
-      curlcurlop.GetMaterialOp(), curlcurlop.GetNDSpaces(),
+  CurlFluxErrorEstimator estimator(
+      curlcurlop.GetMaterialOp(), curlcurlop.GetRTSpace(), curlcurlop.GetNDSpaces(),
       iodata.solver.linear.estimator_tol, iodata.solver.linear.estimator_max_it, 0,
       iodata.solver.linear.estimator_mg);
   ErrorIndicator indicator;
@@ -77,7 +77,7 @@ MagnetostaticSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
     Curl.Mult(A[step], B);
     postop.SetAGridFunction(A[step]);
     postop.SetBGridFunction(B);
-    double E_mag = postop.GetHFieldEnergy();
+    const double E_mag = postop.GetHFieldEnergy();
     Mpi::Print(" Sol. ||A|| = {:.6e} (||RHS|| = {:.6e})\n",
                linalg::Norml2(curlcurlop.GetComm(), A[step]),
                linalg::Norml2(curlcurlop.GetComm(), RHS));
@@ -89,7 +89,7 @@ MagnetostaticSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
 
     // Calculate and record the error indicators.
     Mpi::Print(" Updating solution error estimates\n");
-    estimator.AddErrorIndicator(A[step], indicator);
+    estimator.AddErrorIndicator(B, E_mag, indicator);
 
     // Postprocess field solutions and optionally write solution to disk.
     Postprocess(postop, step, idx, I_inc[step], E_mag,