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)
