Support Floquet periodic boundary conditions

CHANGELOG.md

@@ -45,6 +45,9 @@ The format of this changelog is based on
   - Added adaptive time-stepping capability for transient simulations. The new ODE integrators
     rely on the SUNDIALS library and can be specified by setting the
     `config["Solver"]["Transient"]["Type"]` option to `"CVODE"` or `"ARKODE"`.
+  - 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

@@ -125,7 +125,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"` :  Object for configuring Floquet periodic boundary conditions for
+meshes generated with built-in periodicity.
 
 `"Postprocessing"` :  Top-level object for configuring boundary postprocessing.
 
@@ -527,7 +530,9 @@ boundary.
     {
         "DonorAttributes": [<int array>],
         "ReceiverAttributes": [<int array>],
-        "Translation": [<float array>]
+        "Translation": [<float array>],
+        "AffineTransformation": [<float array>],
+        "FloquetWaveVector": [<float array>]
     },
     ...
 ]
@@ -541,8 +546,26 @@ 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 between the donor
+and receiver attributes 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) between the donor
+and receiver attributes 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/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

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());

palace/fem/libceed/coefficient.cpp

@@ -17,7 +17,7 @@ namespace
 
 inline auto CoeffDim(int dim)
 {
-  return dim * (dim + 1) / 2;
+  return dim * dim;
 }
 
 inline void MakeDiagonalCoefficient(int dim, CeedIntScalar *mat_coeff, CeedScalar a,
@@ -30,7 +30,7 @@ inline void MakeDiagonalCoefficient(int dim, CeedIntScalar *mat_coeff, CeedScala
   }
   for (int di = 0; di < dim; ++di)
   {
-    const int idx = (di * dim) - (((di - 1) * di) / 2);
+    const int idx = di * (dim + 1);
     mat_coeff[coeff_dim * k + idx].second = a;
   }
 }
@@ -86,8 +86,7 @@ PopulateCoefficientContext(int dim, const MaterialPropertyCoefficient *Q, double
     AttrMat(ctx.data())[i].first = (k < 0) ? zero_mat : k;
   }
 
-  // Copy material properties: Matrix-valued material properties are always assumed to be
-  // symmetric and we store only the lower triangular part.
+  // Copy material properties
   ctx[1 + attr_mat.Size()].first = mat_coeff.SizeK() + 1;
   for (int k = 0; k < mat_coeff.SizeK(); k++)
   {
@@ -100,10 +99,10 @@ PopulateCoefficientContext(int dim, const MaterialPropertyCoefficient *Q, double
     {
       for (int dj = 0; dj < dim; ++dj)
       {
-        for (int di = dj; di < dim; ++di)
+        for (int di = 0; di < dim; ++di)
         {
           // Column-major ordering.
-          const int idx = (dj * dim) - (((dj - 1) * dj) / 2) + di - dj;
+          const int idx = (dj * dim) + di;
           MatCoeff(ctx.data())[coeff_dim * k + idx].second = a * mat_coeff(di, dj, k);
         }
       }

palace/fem/libceed/integrator.cpp

@@ -204,7 +204,7 @@ std::vector<CeedInt> QuadratureDataSetup(unsigned int ops, Ceed ceed,
   PalaceCeedCall(ceed, CeedBasisGetNumQuadraturePoints(basis, &num_qpts));
   for (auto size : active_input_sizes)
   {
-    q_data_size += size * (size + 1) / 2;
+    q_data_size += size * size;
   }
 
   PalaceCeedCall(

palace/fem/libceed/operator.hpp

@@ -42,6 +42,9 @@ class Operator : public palace::Operator
   ~Operator() override;
 
   CeedOperator operator[](std::size_t i) const { return op[i]; }
+
+  CeedOperator GetTranspose(std::size_t i) const { return op_t[i]; }
+
   auto Size() const { return op.size(); }
 
   void AddSubOperator(CeedOperator sub_op, CeedOperator sub_op_t = nullptr);