# Numerical fluxes for parabolic interface problems with aligned mesh Xiaoyu Wei <xywei@illinois.edu> $\Omega = \Omega_1 \cup \Gamma \cup \Omega_2$, consider the heat eqn $$ u_t = \nabla\cdot(\alpha(x) \nabla u), $$ subject to homogeneous Neumann BCs on $\partial \Omega$. Here, $\alpha$ is piecewise-smooth, having a jump across $\Gamma$. WLOG, denote $\alpha(x) = \alpha_1(x)\chi_1(x)+ \alpha_2(x)\chi_2(x)$, where $\alpha_{1,2}$ are smooth. **There are two jump conditions** a time stepping scheme needs to preserve. Firt of all, existing [parabolic regularity results](https://doi.org/10.1007/s40879-017-0168-y) show that $u\in W^{1,p}$, meaning $$ u_1(x) = u_2(x), \quad \forall x \in \Gamma \tag{jump-1} $$ Consider a region $D\subset \Omega$, such that $D\cap\Omega_1=D_1$, $D\cap\Omega_2=D_2$, $D\cap\Gamma = \Pi$. Using test function $v$, integration-by-parts over $D_1$ and $D_2$ yields $$0 = \int_{D_1}[u_t v + \alpha_1(x) \nabla u\cdot\nabla v] dx - \int_{\partial D_1 \setminus \Pi} \alpha_1 \nabla u_1 \cdot n_1 v ds - \int_\Pi \alpha_1 \nabla u_1 \cdot n_1 v ds + \\ \int_{D_2}[u_t v + \alpha_2(x) \nabla u\cdot\nabla v] dx - \int_{\partial D_2 \setminus \Pi} \alpha_2 \nabla u_2 \cdot n_2 v ds - \int_\Pi \alpha_2 \nabla u_2 \cdot n_2 v ds \\ = \int_{D}[u_t v + \alpha(x) \nabla u\cdot\nabla v] dx - \int_{\partial D} \alpha \nabla u \cdot n v ds - \int_\Pi [\alpha_1 \nabla u_1 \cdot n_1 + \alpha_2 \nabla u_2 \cdot n_2] v ds \\ = \int_\Pi [\alpha_1 \nabla u_1 \cdot n_1 + \alpha_2 \nabla u_2 \cdot n_2] v ds \\ $$ where the last equality uses the fact that the PDE holds over the whole region $D$. Noticing $n_1 = - n_2$, we have $$ 0= \int_\Pi (\alpha_1 \nabla u_1 - \alpha_2\nabla u_2) \cdot n_1 v ds $$ This identity should hold for any region $D$ that intersects $\Gamma$, so $$ \alpha_1(x) \partial_n u_1(x) = \alpha_2(x) \partial_n u_2(x)\quad \forall x\in \Gamma. \tag{jump-2} $$ Implications on numerical flux choices? Fluxes are linear in updates added to $u$ at each time step. So to preserve jump conditions at the next time step, when deciding what traces of $u$, $\nabla u$ and $\alpha$ to use, need - $\hat{u^+} = \hat{u^-}$ - $\hat{\alpha^+}\hat{u_n^+} = \hat{\alpha^-}\hat{u_n^-}$ **Simple constructions satisfying the conditions:** - use averaged $\alpha$ as well as $u$ and $\nabla u$. - use averaged $u$ and $q=\alpha \nabla u$. (Use $q=\alpha \nabla u$ in the LDG formulation to avoid evaluating $u_n$ and $\alpha$ at separate places.) If we were to use internal (single sided) values of $\alpha$, then we should calculate $u_n$ on the same side as $$ \hat{u_n} = [[\alpha u_n]] / \hat{\alpha}, $$ where $[[.]]$ means some kind of average function. Remark: [Cockburn, Shu 1998](https://doi.org/10.1137/S0036142997316712) used $q=\sqrt{\alpha}\nabla u$ because it allows the proof for $L^2$ stability. Can that result be generalized to interface problems if numerical fluxes are chosen as above to respect jump conditions? (high order with presence of interface). For example, using simple arithmetic average, we have $$ \hat{u_n^+} = \frac{1}{2} u_n^+ + \frac{\alpha^-}{2\alpha^+} u_n^- , $$ which is different from a simple average. Meaning, without such treatment, the jump condidtions may be violated. In fact, the error from using a wrong formula (say, arithmetic mean) to compute $\hat{u_n}$ here is O(1) if $\alpha^+\neq\alpha^-$, because: $$ \alpha^+ (u_n^+ + u_n^-)/2 - \alpha^- (u_n^+ + u_n^-)/2 = (\alpha^+ - \alpha^-) (u_n^+ + u_n^-)/2, $$ which leads to an O(1) term in the weak form (let $D=\Omega$ in the derivation of jump-2, this imbalance in flux is integrated over the whole interface $\Gamma$).