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

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)

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.
• Sorry for the lack of reply, but here it goes: – rchilton1980 Jan 30 '14 at 16:57
• * functions are continuous in the domain, so they can absorb a derivative/integrate by parts. The solutions do match in gamma in strategy (1), they are only close in strategy (2) – rchilton1980 Jan 30 '14 at 16:59
• * means DoF's on elements on the boundary - not just the dof's right on the interface. To see why, strategy (2) requires computing a normal derivative of u, so you have a reach "a little inside" the "other domain" in order to evaluate normal behavior – rchilton1980 Jan 30 '14 at 17:01
• * unfortunately in this case the 1D problem does not capture the essential behavior, because there is no such thing as oblique or evanescent waves on the boundary (which is a big part of understanding the behavior of the robin TC / higher order TC's) This still might be worth doing, I will comment further if I ever get around to that experiment. – rchilton1980 Jan 30 '14 at 17:04
• I do understand why you want to consider the problem in at least 2 dimensions. Nevertheless for debugging I'd always try to simplify the problem as mich as possible. – hrothgar grumpen Jan 31 '14 at 22:48

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!

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.

• At the risk of subverting the structure of the site (you are Q'ing in your A, and I am A'ing in a comment), I'll try to help. The signs of the lagrange multipliers are "induced" by your normal convention (it's d/dn, right?). I always choose the normal to be out of the domain. I think always-inward would be fine too, as long as you are consistent. Regarding references, I think this warrants another question. I work mainly on the Maxwell system, I only framed my question within Helmholtz so that more posters would understand it. So most of my knowledge is drawn from CEM literature .. – rchilton1980 Nov 8 '17 at 13:53
• .. that said, I do have some familiarity with DDM for Helmholtz and will happily supply an answer if you issue a new question for reference requests (I will keep my eyes peeled, and answer, even if not immediately). A new question with more visibility will also get you higher quality answers from other experts on the site. – rchilton1980 Nov 8 '17 at 13:55