Support Floquet periodic boundary conditions

CHANGELOG.md

@@ -48,6 +48,9 @@ The format of this changelog is based on
   - Added an option to use the complex-valued system matrix for the coarse level solve (sparse
     direct solve) instead of the real-valued approximation. This can be specified with
     `config["Solver"]["Linear"]["ComplexCoarseSolve"]`.
+  - Added support for Floquet periodic boundary conditions with phase-delay constraints.
+    The Floquet wave vector can be specified along with periodic boundaries in the
+    `config["Boundaries"]["Periodic"]` configuration.
 
 ## [0.13.0] - 2024-05-20
 

docs/src/config/boundaries.md

@@ -60,6 +60,10 @@
     [
         ...
     ],
+    "FloquetWaveVector":
+    [
+        ...
+    ],
     "Postprocessing":
     {
         "SurfaceFlux":
@@ -125,7 +129,10 @@ electrostatic simulations. Entries of the capacitance matrix are extracted corre
 each terminal boundary.
 
 `"Periodic"` :  Array of objects for configuring periodic boundary conditions for surfaces
-with meshes that are identical after a specified translation.
+with meshes that are identical after translation and/or rotation.
+
+`"FloquetWaveVector"` :  Array for specifying Floquet wave vector for
+meshes generated with built-in periodicity.
 
 `"Postprocessing"` :  Top-level object for configuring boundary postprocessing.
 
@@ -527,7 +534,9 @@ boundary.
     {
         "DonorAttributes": [<int array>],
         "ReceiverAttributes": [<int array>],
-        "Translation": [<float array>]
+        "Translation": [<float array>],
+        "AffineTransformation": [<float array>],
+        "FloquetWaveVector": [<float array>]
     },
     ...
 ]
@@ -541,8 +550,27 @@ attributes for this periodic boundary.
 `"ReceiverAttributes" [None]` :  Integer array of the receiver attributes of the mesh boundary
 attributes for this periodic boundary.
 
-`"Translation" [None]` :  Defines the distance between the donor and receiver attributes in
-mesh units.
+`"Translation" [None]` :  Optional floating point array defining the distance from the donor
+attribute to the receiver attribute in mesh units. If neither `"Translation"` or
+`"AffineTransformation"` are specified, the transformation between donor and receiver boundaries
+is automatically detected.
+
+`"AffineTransformation" [None]` :  Optional floating point array of size 16 defining the
+three-dimensional (4 x 4) affine transformation matrix (in row major format) from the donor attribute
+to the receiver attribute in mesh units. If neither `"Translation"` or `"AffineTransformation"` are
+specified, the transformation between donor and receiver boundaries is automatically detected.
+
+`"FloquetWaveVector" [None]` :  Optional floating point array defining the phase delay between
+this pair of donor and receiver periodic boundaries in the X/Y/Z directions in radians per mesh
+unit. If multiple periodic boundary pairs are used, the Floquet wave vector will be summed over
+the periodic boundary pairs.
+
+## `boundaries["FloquetWaveVector"]`
+
+Optional floating point array defining the phase delay between the periodic boundaries in the X/Y/Z
+directions in radians per mesh unit, for meshes generated with built-in periodicity. This should not
+be used for non-periodic meshes, or for meshes generated without built-in periodicity. In the latter
+case, the Floquet wave vector should be specified via `"boundaries["Periodic"]["FloquetWaveVector"]"`.
 
 ## `boundaries["Postprocessing"]["SurfaceFlux"]`
 

docs/src/config/solver.md

@@ -432,10 +432,14 @@ iterative solver choices, and the default choice depends on the iterative solver
   - `"Default"`
 
 `"DivFreeTol" [1.0e-12]` :  Relative tolerance for divergence-free cleaning used in the
-eigenmode simulation type.
+eigenmode simulation type. Ignored if non-zero Floquet wave vector is specified in
+[`config["Boundaries"]["Periodic"]["FloquetWaveVector"]`](boundaries.md##boundaries%5B%%22Periodic%22%5D%22FloquetWaveVector%22%5D)
+or [`config["Boundaries"]["FloquetWaveVector"]`](boundaries.md##boundaries%5B%%22FloquetWaveVector%22%5D).
 
 `"DivFreeMaxIts" [1000]` :  Maximum number of iterations for divergence-free cleaning use in
-the eigenmode simulation type.
+the eigenmode simulation type. Ignored if non-zero Floquet wave vector is specified in
+[`config["Boundaries"]["Periodic"]["FloquetWaveVector"]`](boundaries.md##boundaries%5B%%22Periodic%22%5D%22FloquetWaveVector%22%5D)
+or [`config["Boundaries"]["FloquetWaveVector"]`](boundaries.md##boundaries%5B%%22FloquetWaveVector%22%5D).
 
 `"EstimatorTol" [1.0e-6]` :  Relative tolerance for flux projection used in the
 error estimate calculation.

docs/src/examples/cylinder.md

@@ -254,15 +254,18 @@ such modes to higher frequencies. The relevant modes are tabulated as
 
 For this problem, we use curved tetrahedral elements from the mesh file
 [`mesh/cavity_tet.msh`](https://github.com/awslabs/palace/blob/main/examples/cylinder/mesh/cavity_tet.msh),
-and the configuration file
-[`waveguide.json`](https://github.com/awslabs/palace/blob/main/examples/cylinder/waveguide.json).
+and the configuration files
+[`waveguide.json`](https://github.com/awslabs/palace/blob/main/examples/cylinder/waveguide.json) and
+[`floquet.json`](https://github.com/awslabs/palace/blob/main/examples/cylinder/floquet.json).
 
-The main difference between this configuration file and those used in the cavity example is
+The main difference between these configuration files and those used in the cavity example is
 in the `"Boundaries"` object: `waveguide.json` specifies a perfect electric conductor
 (`"PEC"`) boundary condition for the exterior surface and a periodic boundary condition
 (`"Periodic"`) on the cross-sections of the cylinder (in the $z-$ direction). The periodic
 attribute pairs are defined by `"DonorAttributes"` and `"ReceiverAttributes"`, and the
-distance between them is given by the `"Translation"` vector in mesh units.
+distance between them is given by the `"Translation"` vector in mesh units. In `floquet.json`,
+an additional `"FloquetWaveVector"` specifies the phase delay between the donor and receiver
+boundaries in the X/Y/Z directions.
 
 The file `postpro/waveguide/eig.csv` contains information about the computed eigenfrequencies and
 associated quality factors:

docs/src/guide/boundaries.md

@@ -74,7 +74,7 @@ Periodic boundary conditions on an existing mesh can be specified using the
 ["Periodic"](../config/boundaries.md#boundaries%5B%22Periodic%22%5D) boundary keyword. This
 boundary condition enforces that the solution on the specified boundaries be exactly equal,
 and requires that the surface meshes on the donor and receiver boundaries be identical up to
-translation. Periodicity in *Palace* is also supported through meshes generated
+translation or rotation. Periodicity in *Palace* is also supported through meshes generated
 incorporating periodicity as part of the meshing process.
 
 ## Lumped and wave port excitation

docs/src/index.md

@@ -62,5 +62,4 @@ the frequency or time domain, using the
   - Improved adaptive mesh refinement (AMR) support for transient simulations and numeric
     wave ports on nonconformal meshes
   - Perfectly matched layer (PML) boundaries
-  - Periodic boundaries with phase delay constraints
   - Automatic mesh generation and optimization

examples/cylinder/floquet.json

@@ -0,0 +1,69 @@
+{
+  "Problem":
+  {
+    "Type": "Eigenmode",
+    "Verbose": 2,
+    "Output": "postpro/floquet"
+  },
+  "Model":
+  {
+    "Mesh": "mesh/cylinder_tet.msh",
+    "L0": 1.0e-2,  // cm
+  },
+  "Domains":
+  {
+    "Materials":
+    [
+      {
+        "Attributes": [1],
+        "Permeability": 1.0,
+        "Permittivity": 2.08,
+        "LossTan": 0.0004
+      }
+    ],
+    "Postprocessing":
+    {
+      "Energy":
+      [
+        {
+          "Index": 1,
+          "Attributes": [1]
+        }
+      ]
+    }
+  },
+  "Boundaries":
+  {
+    "Periodic":
+    [
+      {
+        "DonorAttributes": [2],
+        "ReceiverAttributes": [3],
+        "FloquetWaveVector": [0.0, 0.0, 0.2]
+      }
+    ],
+    "PEC":
+    {
+      "Attributes": [4]
+    }
+  },
+  "Solver":
+  {
+    "Order": 4,
+    "Device": "CPU",
+    "Eigenmode":
+    {
+      "N": 15,
+      "Tol": 1.0e-8,
+      "Target": 2.0,  // TE f111 ~ 2.9 GHz
+      "Save": 15
+    },
+    "Linear":
+    {
+      "Type": "Default",
+      "KSPType": "GMRES",
+      "Tol": 1.0e-8,
+      "MaxIts": 100
+    }
+  }
+}

examples/cylinder/waveguide.json

@@ -39,7 +39,7 @@
       {
         "DonorAttributes": [2],
         "ReceiverAttributes": [3],
-        "Translation": [0.0, 0.0, 5.48] // in L0 units
+        "Translation": [0.0, 0.0, -5.48] // in L0 units
       }
     ],
     "PEC":

palace/drivers/drivensolver.cpp

@@ -8,6 +8,7 @@
 #include "fem/errorindicator.hpp"
 #include "fem/mesh.hpp"
 #include "linalg/errorestimator.hpp"
+#include "linalg/floquetcorrection.hpp"
 #include "linalg/ksp.hpp"
 #include "linalg/operator.hpp"
 #include "linalg/vector.hpp"
@@ -117,17 +118,17 @@ ErrorIndicator DrivenSolver::SweepUniform(SpaceOperator &space_op, PostOperator
   auto C = space_op.GetDampingMatrix<ComplexOperator>(Operator::DIAG_ZERO);
   auto M = space_op.GetMassMatrix<ComplexOperator>(Operator::DIAG_ZERO);
   auto A2 = space_op.GetExtraSystemMatrix<ComplexOperator>(omega0, Operator::DIAG_ZERO);
+  auto FP = space_op.GetFloquetMatrix<ComplexOperator>(Operator::DIAG_ZERO);
   const auto &Curl = space_op.GetCurlMatrix();
 
   // Set up the linear solver and set operators for the first frequency step. The
   // preconditioner for the complex linear system is constructed from a real approximation
   // to the complex system matrix.
   auto A = space_op.GetSystemMatrix(std::complex<double>(1.0, 0.0), 1i * omega0,
                                     std::complex<double>(-omega0 * omega0, 0.0), K.get(),
-                                    C.get(), M.get(), A2.get());
+                                    C.get(), M.get(), A2.get(), FP.get());
   auto P = space_op.GetPreconditionerMatrix<ComplexOperator>(1.0, omega0, -omega0 * omega0,
                                                              omega0);
-
   ComplexKspSolver ksp(iodata, space_op.GetNDSpaces(), &space_op.GetH1Spaces());
   ksp.SetOperators(*A, *P);
 
@@ -147,6 +148,15 @@ ErrorIndicator DrivenSolver::SweepUniform(SpaceOperator &space_op, PostOperator
       iodata.solver.linear.estimator_mg);
   ErrorIndicator indicator;
 
+  // If using Floquet BCs, a correction term (kp x E) needs to be added to the B field.
+  std::unique_ptr<FloquetCorrSolver<ComplexVector>> floquet_corr;
+  if (FP)
+  {
+    floquet_corr = std::make_unique<FloquetCorrSolver<ComplexVector>>(
+        space_op.GetMaterialOp(), space_op.GetPeriodicOp(), space_op.GetNDSpace(),
+        space_op.GetRTSpace(), iodata.solver.linear.tol, iodata.solver.linear.max_it, 0);
+  }
+
   // Main frequency sweep loop.
   int step = step0;
   double omega = omega0;
@@ -164,7 +174,7 @@ ErrorIndicator DrivenSolver::SweepUniform(SpaceOperator &space_op, PostOperator
       A2 = space_op.GetExtraSystemMatrix<ComplexOperator>(omega, Operator::DIAG_ZERO);
       A = space_op.GetSystemMatrix(std::complex<double>(1.0, 0.0), 1i * omega,
                                    std::complex<double>(-omega * omega, 0.0), K.get(),
-                                   C.get(), M.get(), A2.get());
+                                   C.get(), M.get(), A2.get(), FP.get());
       P = space_op.GetPreconditionerMatrix<ComplexOperator>(1.0, omega, -omega * omega,
                                                             omega);
       ksp.SetOperators(*A, *P);
@@ -179,6 +189,12 @@ ErrorIndicator DrivenSolver::SweepUniform(SpaceOperator &space_op, PostOperator
     Curl.Mult(E.Real(), B.Real());
     Curl.Mult(E.Imag(), B.Imag());
     B *= -1.0 / (1i * omega);
+    if (FP)
+    {
+      // Calculate B field correction for Floquet BCs.
+      // B = -1/(iω) ∇ x E - 1/ω kp x E
+      floquet_corr->AddMult(E, B, -1.0 / omega);
+    }
     post_op.SetEGridFunction(E);
     post_op.SetBGridFunction(B);
     post_op.UpdatePorts(space_op.GetLumpedPortOp(), space_op.GetWavePortOp(), omega);

palace/drivers/eigensolver.cpp

@@ -9,6 +9,7 @@
 #include "linalg/arpack.hpp"
 #include "linalg/divfree.hpp"
 #include "linalg/errorestimator.hpp"
+#include "linalg/floquetcorrection.hpp"
 #include "linalg/ksp.hpp"
 #include "linalg/operator.hpp"
 #include "linalg/slepc.hpp"
@@ -36,6 +37,10 @@ EigenSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
   auto K = space_op.GetStiffnessMatrix<ComplexOperator>(Operator::DIAG_ONE);
   auto C = space_op.GetDampingMatrix<ComplexOperator>(Operator::DIAG_ZERO);
   auto M = space_op.GetMassMatrix<ComplexOperator>(Operator::DIAG_ZERO);
+  auto FP = space_op.GetFloquetMatrix<ComplexOperator>(Operator::DIAG_ZERO);
+  auto A2 = space_op.GetExtraSystemMatrix<ComplexOperator>(1.0, Operator::DIAG_ZERO);
+  A2 = nullptr;
+
   const auto &Curl = space_op.GetCurlMatrix();
   SaveMetadata(space_op.GetNDSpaces());
 
@@ -124,11 +129,25 @@ EigenSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
                                           : EigenvalueSolver::ScaleType::NONE;
   if (C)
   {
-    eigen->SetOperators(*K, *C, *M, scale);
+    if (FP)
+    {
+      eigen->SetOperators(*K, *C, *M, *FP, scale);
+    }
+    else
+    {
+      eigen->SetOperators(*K, *C, *M, scale);
+    }
   }
   else
   {
-    eigen->SetOperators(*K, *M, scale);
+    if (FP)
+    {
+      eigen->SetOperators(*K, *M, *FP, scale);
+    }
+    else
+    {
+      eigen->SetOperators(*K, *M, scale);
+    }
   }
   eigen->SetNumModes(iodata.solver.eigenmode.n, iodata.solver.eigenmode.max_size);
   eigen->SetTol(iodata.solver.eigenmode.tol);
@@ -154,7 +173,7 @@ EigenSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
   // Construct a divergence-free projector so the eigenvalue solve is performed in the space
   // orthogonal to the zero eigenvalues of the stiffness matrix.
   std::unique_ptr<DivFreeSolver<ComplexVector>> divfree;
-  if (iodata.solver.linear.divfree_max_it > 0)
+  if (iodata.solver.linear.divfree_max_it > 0 && !FP)
   {
     Mpi::Print(" Configuring divergence-free projection\n");
     constexpr int divfree_verbose = 0;
@@ -165,6 +184,15 @@ EigenSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
     eigen->SetDivFreeProjector(*divfree);
   }
 
+  // If using Floquet BCs, a correction term (kp x E) needs to be added to the B field.
+  std::unique_ptr<FloquetCorrSolver<ComplexVector>> floquet_corr;
+  if (FP)
+  {
+    floquet_corr = std::make_unique<FloquetCorrSolver<ComplexVector>>(
+        space_op.GetMaterialOp(), space_op.GetPeriodicOp(), space_op.GetNDSpace(),
+        space_op.GetRTSpace(), iodata.solver.linear.tol, iodata.solver.linear.max_it, 0);
+  }
+
   // Set up the initial space for the eigenvalue solve. Satisfies boundary conditions and is
   // projected appropriately.
   if (iodata.solver.eigenmode.init_v0)
@@ -242,7 +270,7 @@ EigenSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
   // to the complex system matrix.
   auto A = space_op.GetSystemMatrix(std::complex<double>(1.0, 0.0), 1i * target,
                                     std::complex<double>(-target * target, 0.0), K.get(),
-                                    C.get(), M.get());
+                                    C.get(), M.get(), A2.get(), FP.get());
   auto P = space_op.GetPreconditionerMatrix<ComplexOperator>(1.0, target, -target * target,
                                                              target);
   auto ksp = std::make_unique<ComplexKspSolver>(iodata, space_op.GetNDSpaces(),
@@ -306,6 +334,12 @@ EigenSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
     Curl.Mult(E.Real(), B.Real());
     Curl.Mult(E.Imag(), B.Imag());
     B *= -1.0 / (1i * omega);
+    if (FP)
+    {
+      // Calculate B field correction for Floquet BCs.
+      // B = -1/(iω) ∇ x E - 1/ω kp x E.
+      floquet_corr->AddMult(E, B, -1.0 / omega);
+    }
     post_op.SetEGridFunction(E);
     post_op.SetBGridFunction(B);
     post_op.UpdatePorts(space_op.GetLumpedPortOp(), omega.real());