FEM+DDM applied to scalar Helmholtz - necessity of lagrange multipliers?

Question

I am seeking to understand DDM's and their application to Maxwell's equations, though I am settling for the scalar Helmholtz as a baby step. Unfortunately I have hit a conceptual snag that I could use help with. My model problem is a 2D rectangle [-1,1]x[-1/2,1/2] split at x = 0 into left/right domains ($\Omega_L$ and $\Omega_R$). After expanding $u_L$ with scalar nodal interpolants, and testing with $v_L$ (galerkin, same space for $v_L$), and integration by parts over $\Omega_L$, I am left with the usual:

$\int \left(\nabla u_L \cdot \nabla v_L - k^2 u_L \cdot v_L \right) d\Omega_L - \int \left(v_L \frac{\partial u}{\partial n_L} \right) d\Gamma_L$ = 0

Where $\Gamma_L$ is the boundary of $\Omega_L$ at x = 0. For simplicity, I've omitted all the forcing terms (a point/line source) and boundary terms on surfaces other than $\Gamma_L$ (they're all ABC's or dirichlet scatterers). Because it's a wave equation, a first order Robin TC is used at the DDM interface:

$\frac{\partial u_L}{\partial n_L} + jk\cdot u_L = -\frac{\partial u_R}{\partial n_R} + jk\cdot u_R$

From here, there are two implementations I've pursued, (1) and (2). In strategy (1), which is adapted from a reference document (old colleagues thesis), a lagrange multiplier/auxiliary variable is introduced for each of the fluxes (surface currents) $\frac{\partial u_L}{\partial n_L}$ and $\frac{\partial u_R}{\partial n_R}$:

$\frac{\partial u_L}{\partial n_L} = g_L$, $\frac{\partial u_R}{\partial n_R} = g_R$

The lagrange multipliers $g_L$ are discretized in space over $\Gamma$ with 1D nodal interpolants. Furthermore, the TC itself is tested over $\Gamma$ using the same space for tests $h_L$. This yields:

$\int \left(\nabla u_L \cdot \nabla v_L - k^2 u_L \cdot v_L \right) d\Omega_L - \int \left(v_L \cdot g_L \right) d\Gamma_L$ = 0

and

$\int \left(h_L\cdot g_L\right)d\Gamma_L + jk \int \left(h_L\cdot u_L\right) d\Gamma_L + \int \left(h_L\cdot g_R\right)d\Gamma_L - jk \int \left(h_L\cdot u_R\right) d\Gamma_L = 0$

Repeating the same testing process over the other domain $\Omega_R$ gives you the desired 4 equation/4 unknown system:

$\left[ \begin{array}{cccc} \mathbf H_{vL,uL} & +k \mathbf T_{vL,gL} & \mathbf 0 & \mathbf 0\\ +k\mathbf T_{hL,uL} & +jk \mathbf T_{hL,gL} & -k \mathbf T_{hL,uR}&+jk \mathbf T_{hL,gR}\\ \mathbf 0 & \mathbf 0 & \mathbf H_{vR,uR} & +k \mathbf T_{vR,gR}\\ -k\mathbf T_{hR,uL} & +jk \mathbf T_{hR,gL} & +k \mathbf T_{hR,uR}&+jk \mathbf T_{hR,gR}\\ \end{array} \right] \left[ \begin{array}{c} u_L \\ g_L \\ u_R \\ g_R \\ \end{array} \right] = \mathbf 0 $

Where the $\mathbf H$'s are the 2D integrals of the $\nabla^2-k^2$ operator, and all the $\mathbf T$'s are basically 1D mass/gramian matrices. Note there has been a little bit of tricky rescaling applied to the TC's to make the 2x2 L/R subdomain problems complex-symmetric, to ease the memory requirements to direct solve them. This system can be solved w/ either jacobi splitting or krylov/GMRES method (with block diagonal preconditioning, because the L and R subproblems resemble FEM's with low order ABC's enforced via multiplier's, well posed for all k). This method works as expected, in the following ways:

• the eigenvalues of the $\mathbf M^{-1} \mathbf N$ splitting are all on or within the unit circle, which matches the reflection coefficient behavior of the "continuous" ddm algorithm, where you inject plane wave "error" solutions for $u_L$, $u_R$ and analyze how these errors decay to zero as a function of $k_x$ and $k_y$.

• on a conformal mesh, the fully discrete algorithm yields the exact same answer (to machine precision, or at least the residual tolerance of your iterative solver) as the non-ddm/monolithic FEM approach.

Bottom line is, this (1) approach works and mimics the results from my reference thesis.

My confusion arises because I tried an alternate approach, (2), that I think should yield the same results but in fact does not. In strategy (2), lagrange multipliers are not used, and instead we solve the TC for the $\frac{\partial u_L}{\partial n_L}$...

$\frac{\partial u_L}{\partial n_L} = -jk\cdot u_L -\frac{\partial u_R}{\partial n_R} + jk\cdot u_R$

... and then plug it into the boundary term in the volume test to yield:

$\int \left(\nabla u_L \cdot \nabla v_L - k^2 u_L \cdot v_L \right) d\Omega_L + jk \int \left(v_L \cdot u_L \right) d\Gamma_L - \int \left(v_L \cdot \frac{\partial u_R}{\partial n_R}\right) d\Gamma_L + jk\int \left(v_L \cdot u_R\right) d\Gamma_L$

When you repeat the process on $\Omega_R$, you get the following system:

$\left[ \begin{array}{cccc} \mathbf H & \mathbf H & \mathbf 0 & \mathbf 0\\ \mathbf H & \mathbf H+jk \mathbf T & \mathbf 0 & \mathbf D+jk \mathbf T\\ \mathbf 0 & \mathbf 0 & \mathbf H & \mathbf H\\ \mathbf D+jk \mathbf T & \mathbf 0 & \mathbf H & \mathbf H+jk \mathbf T \\ \end{array} \right] \left[ \begin{array}{c} u_L^{\Omega} \\ u_L^{\Gamma} \\ u_R^{\Omega} \\ u_R^{\Gamma} \\ \end{array} \right] = \mathbf 0 $

This is actually more like what I expected to see, because the 2x2 subdomain blocks are now clearly classic FEM's with ABC's on their boundaries (the $\mathbf H+jk\mathbf T$ stuff). The $\mathbf D$ "derivative" operator discretizes the $\int \left(v_L \cdot \frac{\partial u_R}{\partial n_R}\right) d\Gamma_L$, kinda like you're using it for forcing (Robin) data. This approach works in that it converges and yields an answer that appears correct, but there are some things I don't like:

• both the preconditioned krylov and block jacobi iterations converge - but actually faster than they should. In brief, the first order Robin TC should yield unit reflection coefficients for evanescent fourier modes - which motivates the quest for higher order tc's. But for this algorithm all the eigenvalues of $\mathbf M^{-1} \mathbf N$ are considerably within the unit circle - NONE are actually on it, while this is good, it's better than theory suggests, which makes me distrust it/think it's wrong.
• On a conformal mesh, the converged DMM result is quite close to the monolithic FEM result, but the error is far above machine precision. I drive GMRES to 10 digits, the backwards residual of the DDM Ax=b is 10 digits down, but the DDM and the monolithic FEM only agree to 3-4 digits. This is giving me heartburn / trustworthiness issues.

So the question (at long last) is this: have I commited some crime in strategy (2) when I forgoed the use lagrange multipliers to replace the flux of u, and instead just differentiated the primal variable directly via $\mathbf D$? If strategy (1) is the only way that really works, and (2) is junk, that's fine, I'll do (1) - it's more just understanding why (2) doesn't work, so that I won't commit the same crime again unknowingly.

Unfortunately, I am not super familiar with a lot of DDM literature/techniques, so this issue might be covered in a pretty basic text on a simpler model problem (e.g. laplace/poisson?). Thanks in advance - I know this question is very long but just wanted to be as thorough as possible. I am pretty careful (normally) about implementation correctness, so I suspect there is some misunderstanding in the (2) formulation, not it's implementation (but no promises)

Answer 1

As noone else is replying so far I'll try to make remarks and question that could check helf for your search for deeper insight.

• One question from my side concerning your discretization within the subdomains: You say your ansatz functions are "nodal interpolants". Does that imply that the ansatz functions are continuous in each subdomain or are you using piecewise polynomials which are discontinuous? Which order within each element are you using? If you're using continuous functions, that means that in your DD-solution you should have $u_L|_{\Gamma} = u_R|_{\Gamma}$ if the solution coincides with the monolithic approach, right? Did you check for that?
• Which degrees of freedom do you mean with $u_L^{\Gamma}$? Only the nodes at the interface or all degrees of freedom which are on elements at the boundary?
• Have you tried the problem in 1D? If you do it in 1D it is perhaps even simpler to investigate certain properties of your discretization, e.g. the system matrix. I suggest you try a simple 1D problem. Than you solve both problems with a direct solver (just to exclude another error source at the beginning). You could further check the eigenvalues of the system matrices for both cases and compare them.

Answer 2

I have not been able to resolve completely why Strategy (2) doesn't work.

My best intuition is that for a given PDE/mono-fem-solution, there always exists some choice of boundary dirichlet/neumann/robin data that will give you the same (restricted) solution in each subdomain. But it's not necessarily computable from the $u$ values of adjacent domains, or even the entire $u$ state over every other domain. The correct robin data could be something totally different, it's unknown. That's what the addition of new unknowns $\lambda$ (lagrange multipliers) represents - the "correct" but unknown robin data. They must be solved concurrently with $u$, and should have similar "galerkin-ness" property as $u$ (in the sense that error in $\lambda$ should be constrained orthogonal to the space that represents it).

This answer has a lot of squish, but I have pressed on. Happy to accept something better if it ever appears!

Answer 3

I am trying to do something similar to what you are writing here. In short what I understand is that, in the first approach, you do not explicitly compute the $\partial u /\partial n$ operator while in the second you do this.

My question here is about the signs of $gL$ and $gR$ in your first formulation. How do you decide on this ? If I am using a CG-FEM discretization is this the same as FETI-H approach ? Do you have some paper or some work where you explain this in greater depth ?

My current project is precisely working on TC for DDT problems in frequency domain for Helmholtz problems.