Periodic boundary conditions

[kpobrien]

Thanks for developing and maintaining Palace. It's awesome to have an open-source EM solver for superconducting circuits. I have a few questions about periodic boundary conditions in Palace:

1. Does Palace support periodic boundary conditions? I searched the code and docs and did not see it, but may have missed it.
2. Are periodic boundary conditions a feature you plan to add?
3. If not, would you potentially accept a PR adding this feature?

A bit of motivation for the request: For periodic structures (eg. metamaterials, traveling wave parametrric amplifiers (TWPAs), antenna arrays, etc) it is important to have periodic boundary conditions, so that one can simulate a single unit cell and extract properties of the larger system. In addition, for many of these applications, it is important to be able to apply a phase between the periodic boundaries in order to simulate wave propagation along that axis. This is sometimes called Floquet or Bloch periodicity.

[sebastiangrimberg]

Hi @kpobrien, thanks a lot for your interest! The short answer is that currently Palace does not have a periodic BC support but we would absolutely welcome it. We've gotten away with just manually building models which contain many repeated unit cells and using the fact that Palace scales quite well to just simulate the larger problem (without needing periodic BC). This doesn't work all the time once you have very large arrays, clearly, and I'd be happy to work with you or whomever else to try to design and implement it.

The longer answer: MFEM (the finite element library on which Palace is based), does have native support for periodic meshes and thus traditional periodic BCs (no phase delay). See for example this comment: https://github.com/mfem/mfem/issues/649#issuecomment-432773499. If you are sure that your mesh is periodic, we could add the necessary vertex matching as a preprocessing step and everything else should just work.

For the phase delay, this requires a more general approach since instead of constraining the dofs on the slave boundary which map 1:1 with counterparts on the master, there is an added coefficient (dof_slave = alpha * dof_master). I think for this we likely will have to not use MFEM's periodic mesh approach and instead will have to do the elimination on the assembled linear systems by augmenting the current way of handling Dirichlet BC.

Maybe this was too in the weeds but hopefully it is helpful in starting down this path. Thanks again for your interest and willingness to contribute this feature!

[kpobrien]

Thank you @sebastiangrimberg! Your going into the weeds was very helpful.

If you are sure that your mesh is periodic, we could add the necessary vertex matching as a preprocessing step and everything else should just work.

I think it is a good idea to enforce a periodic mesh for both the periodic and periodic with phase delay cases, especially given how MFEM works. In my experience as a FEM user, even if the software supports periodic boundary conditions with an arbitrary mesh, a periodic mesh results in a lower error.

I did some digging and found this issue https://github.com/mfem/mfem/issues/1553 where the MFEM developers recommend solving a periodic system with phase delay by incorporating the phase delay into the PDE so that the modified PDE is periodic and one can use the native periodic boundary conditions. In other words, applying Bloch's theorem. What do you think of this approach?

[sebastiangrimberg]

Sorry for the slow response here! This (https://github.com/mfem/mfem/issues/1553) makes total sense and makes the implementation nicer too I think. If I follow correctly, for the actual BC Palace sees we only will have to deal with the exactly periodic case which is handled already by MFEM and will require probably one or two lines of code in Palace to set up the periodic mesh boundaries from the configuration file data.

Then, for enforcing a particular phase delay, it seems like it will just take someone going through the frequency domain formulation (https://awslabs.github.io/palace/stable/reference/#Mathematical-background) with the correct expansion for the E-field phasor with the delay term and see where the coefficients or extra terms show up as functions of the wave number, propagation direction, and phase delay. I can help with adding these to the frequency domain operators in Palace once we get there. It might take a bit to settle on a nice way to describe these boundary pairs in the configuration file (probably OK to restrict to periodic boundary pairs exactly aligned with the X/Y/Z axes for now, so that the vector wavenumber in the e^{ik \cdot x} simplifies when plugged into the PDE).

[kpobrien]

Thanks @sebastiangrimberg!

If I follow correctly, for the actual BC Palace sees we only will have to deal with the exactly periodic case which is handled already by MFEM and will require probably one or two lines of code in Palace to set up the periodic mesh boundaries from the configuration file data.

That's my understanding as well. Great to hear this won't require major changes to Palace.

Then, for enforcing a particular phase delay, it seems like it will just take someone going through the frequency domain formulation (https://awslabs.github.io/palace/stable/reference/#Mathematical-background) with the correct expansion for the E-field phasor with the delay term and see where the coefficients or extra terms show up as functions of the wave number, propagation direction, and phase delay.

I will happily volunteer to do this and post the results here (unless someone beats me to it). I'll try to find time to look at this in the next few days.

I can help with adding these to the frequency domain operators in Palace once we get there.

Perfect! That sounds like the hard part to me, so I greatly appreciate your help.

It might take a bit to settle on a nice way to describe these boundary pairs in the configuration file (probably OK to restrict to periodic boundary pairs exactly aligned with the X/Y/Z axes for now, so that the vector wavenumber in the e^{ik \cdot x} simplifies when plugged into the PDE).

Yes, I agree we might have to iterate a bit to find the right path, even at the level of whether the user should specify the k vector or phases. That restriction of pairs aligned with the axes would be reasonable, as long as we can apply X,Y,Z phases simultaneously. For example, for mapping out the band structure of a photonic crystal you may want to simultaneously create phase differences along different axes. Let's see if it turns out to be simpler to add these restrictions or solve the more general case!

[EMinsight]

I am wondering how mfem to enforce the two side of Bloch boundary have exact same mesh, this usually done in the mesh engine.

[kpobrien]

Here is a first pass (following the conventions and formatting in the docs except switching from \bm to \mathbf due to limitations of Github's markdown):

We assume that the electric field is a periodic field $\mathbf{E}_p$ multiplied by a phase factor $\mathbf{E}(\mathbf{x})=\mathbf{E}_p(\mathbf{x})e^{-i \mathbf{k}_p \cdot \mathbf{x}}$ where $\mathbf{k}_p$ is the user specified Bloch wavevector. The curl operator acting on $\mathbf{E}(\mathbf{x})$ gives:

\begin{aligned}
\nabla\times\mathbf{E} &= (\nabla\times\mathbf{E_p})e^{-i \mathbf{k}_p \cdot \mathbf{x}} -i\mathbf{k}_p\times \mathbf{E}_p e^{-i \mathbf{k}_p \cdot \mathbf{x}} \\  &= ((\nabla-i\mathbf{k}_p)\times\mathbf{E_p})e^{-i \mathbf{k}_p \cdot \mathbf{x}}
\end{aligned}

We can either modify the curl operator $\nabla \to \nabla -i\mathbf{k}_p$ and multiply the electric field by $e^{-i \mathbf{k}_p \cdot \mathbf{x}}$ (outside the curl) or substitute it in and include the extra terms:

\begin{aligned}
\nabla\times\mu_r^{-1}\nabla\times\mathbf{E}_p
    -i\mathbf{k}_p\times\mu_r^{-1}\nabla\times\mathbf{E}_p
    -i\nabla\times(\mu_r^{-1}\mathbf{k}_p\times\mathbf{E}_p)
    -\mathbf{k}_p\times\mu_r^{-1}\mathbf{k}_p\times\mathbf{E}_p
    + i\omega\sigma\mathbf{E}_p
    - \omega^2\varepsilon_r\mathbf{E}_p &= 0 \,,\; \mathbf{x}\in\Omega \\
\mathbf{n}\times\mathbf{E}_p &= 0 \,,\; \mathbf{x}\in\Gamma_{PEC} \\
\mathbf{n}\times(\mu_r^{-1}\nabla\times\mathbf{E}_p)
    -i\mathbf{n}\times(\mu_r^{-1}\mathbf{k}_p\times\mathbf{E}_p) &= 0 \,,\; \mathbf{x}\in\Gamma_{PMC} \\
\mathbf{n}\times(\mu_r^{-1}\nabla\times\mathbf{E}_p)
    -i\mathbf{n}\times(\mu_r^{-1}\mathbf{k}_p\times\mathbf{E}_p)
    + \gamma\mathbf{n}\times(\mathbf{n}\times\mathbf{E}_p) &= \mathbf{U}^{inc}e^{i \mathbf{k}_p \cdot \mathbf{x}} \,,\; \mathbf{x}\in\Gamma_{Z}
\end{aligned}

[sebastiangrimberg]

Nice!! Any chance you want to derive the weak form too? Or is this something I should work on 😁 (It's not immediately clear the what the best way is to handle the discretization yet, need to either find it in literature or derive it ourselves).

[sebastiangrimberg]

I am wondering how mfem to enforce the two side of Bloch boundary have exact same mesh, this usually done in the mesh engine.

@EMinsight, yes, MFEM requires the mesh is actually periodic before performing a matching between pairs of faces: https://github.com/mfem/mfem/blob/master/mesh/mesh.cpp#L4973C44-L4973C44. If the user instead encodes the periodicity directly into the mesh file, this is handled "automagically" by the mesh parser in MFEM.

[kpobrien]

@sebastiangrimberg, do you have a writeup or could you point me to a reference or the location in the code where you have the weak form of the unmodified Maxwell's equations as used in Palace? It's getting to a busy point in the semester for me, but I'll help out as I can.

[sebastiangrimberg]

Here's one: https://arxiv.org/pdf/1910.10069.pdf (and associated textbook).

Anyway, I've gone ahead and started looking at the extra terms. Changing the convention to $\mathbf{E}=\mathbf{E}_pe^{+i\mathbf{k}_p\cdot\mathbf{x}}$ to be consistent with the phasor definition we use ($\mathcal{E}(t)=\mathbf{E}e^{+i\omega t}$), the weak form for the curl-curl term becomes (for some vector test function $\mathbf{F}$, and with $\hat{\nabla} = \nabla + i\mathbf{k}_p$):

\begin{aligned}
(\hat{\nabla}\times(\mu^{-1}\hat{\nabla}\times\mathbf{E}_p), \mathbf{F})_{\Omega} &= (\mu^{-1}\nabla\times\mathbf{E}_p,\nabla\times\mathbf{F})_{\Omega}  \\
&\quad+ i(\mu^{-1}\mathbf{k}_p\times\mathbf{E}_p,\nabla\times\mathbf{F})_{\Omega}  \\
&\quad- i(\mu^{-1}\nabla\times\mathbf{E}_p,\mathbf{k}_p\times\mathbf{F})_{\Omega}  \\
&\quad+ (\mu^{-1}\mathbf{k}_p\times\mathbf{E}, \mathbf{k}_p\times\mathbf{F})_{\Omega}  \\
&\quad- (\mu^{-1}\hat{\nabla}\times\mathbf{E}_p,\mathbf{n}\times\mathbf{F})_{\partial\Omega}
\end{aligned}

The first term on the RHS is the standard stiffness term, and fourth is just a contribution to the standard mass term with coefficient $[\mathbf{k}_p\times]^T\mu^{-1}[\mathbf{k}_p\times]$ where $[\mathbf{k}_p\times]$ is the matrix representation of the cross-product with $\mathbf{k}_p$.

The last term comes from the integration by parts and is the natural boundary condition.

The second and third terms are new, and should be fine to implement. One thing that worries me is that

i(\mu^{-1}\mathbf{k}_p\times\mathbf{E}_p,\nabla\times\mathbf{F})_{\Omega} -i(\mu^{-1}\nabla\times\mathbf{E}_p,\mathbf{k}_p\times\mathbf{F})_{\Omega}

is skew symmetric, which is the first symmetry-breaking term we have encountered. This might have implications for linear solver performance, and it should probably be double checked that I haven't messed up a minus sign somewhere in my algebra.

Anyway, I think this is a reassuring starting point if we can double check that this is correct. The $\mathbf{k}_p$ vector should be specified in the configuration file and will end up in a new PeriodicBoundaryOperator class, which SpaceOperator should handle for the addition of the new terms in the governing equations (the new entries will be described by MixedVectorWeakCurl and MixedVectorCurl integrators, I believe).

Then we'll have to figure out what will need to be done about the asymmetry of the operator and if the current linear solver will need any adjustments.

Lastly, for my own understanding, what kind of driving excitations do you need for these types of simulations? I imagine what is currently in Palace might not be sufficient. Relatedly, what are the quantities of interest you care about postprocessing from the computed fields?

[kpobrien]

Thanks for the references and for starting on the weak form derivations! Regarding the sign convention, the argument for my choice of $\mathbf{E}(\mathbf{x}) = \mathbf{E}_p e^{-i \mathbf{k}_p \cdot \mathbf{x}}$ is that using the phasor convention for Palace, a plane wave propagating in the direction $\mathbf{k}_p/|\mathbf{k}_p|$ would be $\mathbf{E}_0 e^{-i \mathbf{k}_p \cdot \mathbf{x}+i\omega t}$. Imagining solving this in a box with periodic boundary conditions on all sides with user specified wavevector $\mathbf{k}_p$, the periodic field $\mathbf{E}_p$ should be constant. If we drop the $i\omega t$ term, then $\mathbf{E}_0 e^{-i \mathbf{k}_p \cdot \mathbf{x}} = \mathbf{E}_p e^{-i \mathbf{k}_p \cdot \mathbf{x}}$ which implies $\mathbf{E}_0 = \mathbf{E}_p$ or in other words, the field we would solve in MFEM would be constant as desired in this case.

Anyway, I think this is a reassuring starting point if we can double check that this is correct. The \mathbf{k}_p vector should be specified in the configuration file and will end up in a new PeriodicBoundaryOperator class, which SpaceOperator should handle for the addition of the new terms in the governing equations (the new entries will be described by MixedVectorWeakCurl and MixedVectorCurl integrators, I believe).

I'll try to double check that and this all sounds good.

The use case I had in mind was eigenmode simulations where there is no excitation. The output of that simulation would be the fields / eigenmode for the user defined $\mathbf{k}_p$ and the (potentially complex) eigenfrequency. In postprocessing, assuming the mode is TEM, I would compute line integrals of E and H to calculate voltage and current and/or flux with the end goal of calculating the characteristic impedance and wavevector as a function of frequency. For different applications like antennas, driven simulations from a lumped port or current source or with a background field (eg. a plane wave) which is scattered by the antenna are typical excitation conditions.

[pocdn]

This would be a great addition I would be using now if it existed.

An example of using eigenmode solver with swept periodic boundary phase in competing HFSS (Master+Slave GUI definition) to compute dispersion diagram is added in links below. These links give more detail on how this periodic boundary change is used and would be a welcome addition to palace.

Tutorial on Dispersion Diagram I: Parallel Plate

Youtube: (1/3) Dispersion diagram plot-HFSS (Part-1: Simulation setup for periodic boundary condition) by RF Design Basics