Coupling Boundary Condition of one PDE with source term of another PDE

Question

We have a system of equations, wherein the BC of one PDE is coupled with the source term of another PDE.

We have a regular 2D unit grid in x and y.

There are two PDEs to be solved

• The first PDE (elliptic diffusion problem) is defined only at $y = 1$, acting along the x-axis (i.e. it acts in the x-direction and only along the top of the cartesian grid). This x-axis is discretized with a fixed grid-spacing, generating a finite number of nodes. Let this set of nodes/co-ordinates be represented by $‘X’$.
• The second PDE is a time-varying diffusion problem. This is defined only along the y-axis, but for all x-nodes (i.e. for all $‘X’$), where the 1st PDE is being solved.

PDE1:

$$\nabla.(S \nabla a) = f(\omega,a)$$ BC1: $\frac{\partial a}{\partial x} = 1$, at $x=1$ (Neumann)

BC2: $a = 0$, at $x=0$ (Dirichelet)

PDE2:

$$\frac{\partial B}{\partial t} = \nabla.( \left( \begin{array}{cc} 0 & 0 \\ 0 & D \end{array} \right) \nabla B)$$

BC1: $B = 0$, at $y=0$ for all $'X'$, i.e. along the bottom face

BC2: $\frac{\partial B}{\partial y} = g(\omega, a)$, at $y=1$ (Neumann)

$f$ and $g$ are linear functions in $\omega(x,y,t)$ and $a(x,y,t)$. $B$ is defined in the 2-D grid as $B(x,y,t)$.

Importantly, $$\omega=B_{y=1} \text{ for all } 'X'$$

i.e., the BC2 of the 2nd PDE couples with the Implicit Source Term of the 1st PDE along the top face of the Cartesian mesh.

How would one approach this problem in order to obtain a numerical solution via a FD/FV implementation?

Answer 1

As written, it's not clear to me that it's actually fully coupled. That is, the solution to $B$ depends (via its BC2) on the solution to $a$, but assuming $\omega$ is simply some specified function of position and time, it looks like $a$ doesn't depend on the solution to $B$. If that's the case, then simply solve for $a$ on a 1D mesh first, then you have the (constant in time) BC for $B$ and you can solve it as though there were no coupling. This can be done using a number of techniques for applying BC's within the FD/FV method, but it didn't appear that this was what you were asking. If $\omega$ is a function of $B$, then the system would be fully coupled. I'll assume that they're fully coupled.

If you're confident that the diffusion in PDE2 is not going to ever have a component in the $x$ direction, and you think FD/FV is good for your problem (I so no reason why they wouldn't be), I would consider casting this as two separate, coupled 1D grids to avoid needlessly calculating zero-fluxes in the $x$ direction in the diffusing region. The first grid would be 1D in the $x$ direction, over which PDE1 can be defined. Then, you could repeat a 1D $y$ direction grid at each discretization point in $x$, and define PDE2 over each one of those. Then it's more like an equivalent 1+1D system. To solve that in a coupled system, as you suggested, you can use FD/FV to discretize in space, then step forward using the method of lines. Because $a$ doesn't appear with a time-derivative anywhere, the discretized problem will be a system of (index-1) differential algebraic equations (DAE's).

Once you've done the FD/FV, you would end up with an equation like this $$0 = f_1(\omega, \{a_i\})$$ at each interior mesh point in the $x$ mesh, and $$\frac{\partial B}{\partial t} = f_2(\{B_i\})$$ at each interior mesh point on the $y$ meshes where $f_1$ and $f_2$ are functions which depend on your discretization choice.

There are a number of ways to solve index-1 DAE's including (1) simple implicit Euler time stepping, (2) the idas module of the SUNDIALS suite, which has been wrapped for some other languages (Python, Julia...), or perhaps (3) Matlab's ode15s, which also enables solution of index-1 DAE's by passing in a singular mass matrix.

Then, as you said, the coupling would come through the boundary conditions, which you can implement a number of ways. One common approach is to use "ghost points," or fictitious points added on each side of the mesh. So for example, for the right side of the $x$ mesh, say it has index $N_x$, you could add an additional equation for the fictitious point $N_{x}+1$, $$a_{N_{x}+1} = a_{N_x}$$ and simply write the grid point equation for $a_{N_x}$ as you would for any interior point (typically as a function of $a_{N_{x}-1}, a_{N_x}, a_{N_{x}+1}$). The above is a case that's trivial enough that you could choose to substitute this algebraic equation into the discretization for point $N_x$ on the $x$ mesh, but when things get uglier, you can simply leave them as extra equations. For example, looking at a top point for the $B$ mesh(es) at $x$-index $x_i$, you could have the $y$-index point $N_{y}+1$ defined with the equation $$\left(B_{x_i,N_{y}+1} - B_{x_i,N_y}\right)/\Delta y = g(\omega,a_{x_i})$$ where $\Delta y$ is the grid spacing in the $y$ direction. Then, again, if you can solve for $B_{x_i,N_y}$, you could substitute it into the $x_i, N_{y}$ equation, but if $\omega$ were some non-linear function of $B$, you could solve this algebraic equation along with the others in your DAE system. You can certainly change the discretization scheme to get higher order accuracy, but hopefully this demonstrates the idea.