Turn on operator partial assembly by default

CHANGELOG.md

@@ -13,6 +13,12 @@ The format of this changelog is based on
 
 ## In progress
 
+  - Changed default value of `config["Solver"]["PartialAssemblyOrder"]` in order to activate
+    operator partial assembly by default for all operators in all simulation types.
+  - Changed the normalization of computed eigenmodes for consistency across different domain
+    decompositions. Eigenvectors are now normalized with respect to the mass matrix for unit
+    domain electric field energy.
+
 ## [0.12.0] - 2023-12-21
 
   - Added support for operator partial assembly for high-order finite element spaces based

docs/src/config/solver.md

@@ -43,9 +43,11 @@ with
 
 `"Order" [1]` :  Finite element order (degree). Arbitrary high-order spaces are supported.
 
-`"PartialAssemblyOrder" [100]` :  Order at which to switch from full assembly of finite
+`"PartialAssemblyOrder" [1]` :  Order at which to switch from full assembly of finite
 element operators to [partial assembly](https://mfem.org/howto/assembly_levels/). Setting
-this parameter equal to 1 will fully activate operator partial assembly on all levels.
+this parameter equal to 1 will fully activate operator partial assembly on all levels, while
+setting it to some large number (greater than the finite element order) will result in
+fully assembled operators as sparse matrices.
 
 `"Device" ["CPU"]` :  The runtime device configuration passed to
 [MFEM](https://mfem.org/howto/assembly_levels/) in order to activate different options

docs/src/examples/cavity.md

@@ -70,9 +70,8 @@ and
 [`cavity_impedance.json`](https://github.com/awslabs/palace/blob/main/examples/cavity/cavity_impedance.json).
 
 In both, the [`config["Problem"]["Type"]`](../config/problem.md#config%5B%22Problem%22%5D)
-field is set to `"Eigenmode"`, and we use the mesh shown above with a single level of
-uniform mesh refinement (`"UniformLevels": 1`). The material properties for Teflon are
-entered under
+field is set to `"Eigenmode"`, and we use the mesh shown above. The material properties for
+Teflon are entered under
 [`config["Domains"]["Materials"]`](../config/domains.md#domains%5B%22Materials%22%5D). The
 [`config["Domains"]["Postprocessing"]["Energy]"`](../config/domains.md#domains%5B%22Postprocessing%22%5D%5B%22Energy%22%5D)
 object is used to extract the quality factor due to bulk dielectric loss; in this problem
@@ -90,10 +89,10 @@ In both cases, we configure the eigenvalue solver to solve for the ``15`` lowest
 modes above ``2.0\text{ GHz}`` (the dominant mode frequencies for both the
 ``\text{TE}`` and ``\text{TM}`` cases fall around ``2.9\text{ GHz}`` frequency for this
 problem). A sparse direct solver is used for the solutions of the linear system resulting
-from the spatial discretization of the governing equations, using in this case a second-
+from the spatial discretization of the governing equations, using in this case a fourth-
 order finite element space.
 
-The frequencies for the lowest order ``\text{TE}`` and ``\text{TM}`` modes computed using
+The frequencies for the lowest-order ``\text{TE}`` and ``\text{TM}`` modes computed using
 the above formula for this problem are listed in the table below.
 
 | ``(n,m,l)`` | ``f_{\text{TE}}``       | ``f_{\text{TM}}``       |

docs/src/examples/cpw.md

@@ -80,9 +80,8 @@ sweep algorithm is given a tolerance of ``1\times10^{-3}`` for choosing the samp
 points; the simulation with uniform ports uses ``9`` frequency samples and that with wave
 ports uses ``10``. Despite the much finer frequency resolution, the adaptive frequency
 sweep simulations take roughly the same amount of time as the uniform ones where the
-resulting resolution is worse by a factor of ``20``. Lastly, for all simulations, a single
-level of uniform mesh refinement is applied to the initial mesh and a first-order finite
-element approximation for the solution is used.
+resulting resolution is worse by a factor of ``20``. Lastly, for all simulations, a second-
+order finite element approximation for the solution is used.
 
 The results from the four different simulations are presented in the plots below. Note that
 here, ``\text{dB}`` means ``20\log_{10}(|S_{ij}|)``:
@@ -105,8 +104,8 @@ ports, namely the lumped port excitation exhibits much higher reflection that th
 ports. This can be attributed to the coarse meshes used for these examples. Indeed,
 refining the mesh or increasing the order of the solution approximation resolves this issue
 and leads to better agreement between the lumped port and wave port results. See below for
-the results with again a single level of mesh refinement but ``p = 2`` for the order of the
-solution space.
+the results with again ``p = 2`` for the order of the solution space but with a single level
+of mesh refinement as well.
 
 ```@raw html
 <br/><p align="center">

examples/cavity/cavity_impedance.json

@@ -9,10 +9,7 @@
   {
     "Mesh": "mesh/cavity.msh",
     "L0": 1.0e-2,  // cm
-    "Refinement":
-    {
-      "UniformLevels": 1
-    }
+    "Refinement": {}
   },
   "Domains":
   {
@@ -48,7 +45,7 @@
   },
   "Solver":
   {
-    "Order": 2,
+    "Order": 4,
     "Eigenmode":
     {
       "N": 15,

examples/cavity/cavity_pec.json

@@ -9,10 +9,7 @@
   {
     "Mesh": "mesh/cavity.msh",
     "L0": 1.0e-2,  // cm
-    "Refinement":
-    {
-      "UniformLevels": 1
-    }
+    "Refinement": {}
   },
   "Domains":
   {
@@ -45,7 +42,7 @@
   },
   "Solver":
   {
-    "Order": 2,
+    "Order": 4,
     "Eigenmode":
     {
       "N": 15,

examples/cpw/cpw_lumped_adaptive.json

@@ -9,10 +9,7 @@
   {
     "Mesh": "mesh/cpw_lumped.msh",
     "L0": 1.0e-6,  // μm
-    "Refinement":
-    {
-      "UniformLevels": 1
-    }
+    "Refinement": {}
   },
   "Domains":
   {
@@ -150,7 +147,7 @@
   },
   "Solver":
   {
-    "Order": 1,
+    "Order": 2,
     "Driven":
     {
       "MinFreq": 2.0,  // GHz

examples/cpw/cpw_lumped_uniform.json

@@ -9,10 +9,7 @@
   {
     "Mesh": "mesh/cpw_lumped.msh",
     "L0": 1.0e-6,  // μm
-    "Refinement":
-    {
-      "UniformLevels": 1
-    }
+    "Refinement": {}
   },
   "Domains":
   {
@@ -150,7 +147,7 @@
   },
   "Solver":
   {
-    "Order": 1,
+    "Order": 2,
     "Driven":
     {
       "MinFreq": 2.0,  // GHz

examples/cpw/cpw_wave_adaptive.json

@@ -9,10 +9,7 @@
   {
     "Mesh": "mesh/cpw_wave.msh",
     "L0": 1.0e-6,  // μm
-    "Refinement":
-    {
-      "UniformLevels": 1
-    }
+    "Refinement": {}
   },
   "Domains":
   {
@@ -114,7 +111,7 @@
   },
   "Solver":
   {
-    "Order": 1,
+    "Order": 2,
     "Driven":
     {
       "MinFreq": 2.0,  // GHz

examples/cpw/cpw_wave_uniform.json

@@ -9,10 +9,7 @@
   {
     "Mesh": "mesh/cpw_wave.msh",
     "L0": 1.0e-6,  // μm
-    "Refinement":
-    {
-      "UniformLevels": 1
-    }
+    "Refinement": {}
   },
   "Domains":
   {
@@ -114,7 +111,7 @@
   },
   "Solver":
   {
-    "Order": 1,
+    "Order": 2,
     "Driven":
     {
       "MinFreq": 2.0,  // GHz

palace/drivers/eigensolver.cpp

@@ -148,13 +148,13 @@ EigenSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
   std::unique_ptr<Operator> KM;
   if (iodata.solver.eigenmode.mass_orthog)
   {
-    // Mpi::Print(" Basis uses M-inner product\n");
-    // KM = spaceop.GetInnerProductMatrix(0.0, 1.0, nullptr, M.get());
-    // eigen->SetBMat(*KM);
-
-    Mpi::Print(" Basis uses (K + M)-inner product\n");
-    KM = spaceop.GetInnerProductMatrix(1.0, 1.0, K.get(), M.get());
+    Mpi::Print(" Basis uses M-inner product\n");
+    KM = spaceop.GetInnerProductMatrix(0.0, 1.0, nullptr, M.get());
     eigen->SetBMat(*KM);
+
+    // Mpi::Print(" Basis uses (K + M)-inner product\n");
+    // KM = spaceop.GetInnerProductMatrix(1.0, 1.0, K.get(), M.get());
+    // eigen->SetBMat(*KM);
   }
 
   // Construct a divergence-free projector so the eigenvalue solve is performed in the space
@@ -272,6 +272,14 @@ EigenSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
       spaceop.GetMaterialOp(), spaceop.GetNDSpace(), iodata.solver.linear.estimator_tol,
       iodata.solver.linear.estimator_max_it, 0);
   ErrorIndicator indicator;
+  if (!KM)
+  {
+    // Normalize the finalized eigenvectors with respect to mass matrix (unit electric field
+    // energy) even if they are not computed to be orthogonal with respect to it.
+    KM = spaceop.GetInnerProductMatrix(0.0, 1.0, nullptr, M.get());
+    eigen->SetBMat(*KM);
+    eigen->RescaleEigenvectors(num_conv);
+  }
   for (int i = 0; i < iodata.solver.eigenmode.n; i++)
   {
     eigen->GetEigenvector(i, E);