# 1 D Diffusion equation FDM with different layers

I'm trying to solve this particular equation

$\frac{\partial u}{\partial t} = \frac{\partial}{\partial x} \big[D_{i}(x)\frac{\partial u}{\partial x} \big] + S(x,t)$

where the $i$ index denotes different layers. Let's say that $i \in \{1,2,3\}$.

The interface conditions are:

$-D_1(x)\frac{\partial u}{\partial x}\bigg|_{x=x_1^-} = -D_2(x)\frac{\partial u}{\partial x}\bigg|_{x=x_1^+}$

$-D_2(x)\frac{\partial u}{\partial x}\bigg|_{x=x_2^-} = -D_3(x)\frac{\partial u}{\partial x}\bigg|_{x=x_2^+}$

and at the boundaries:

$\frac{\partial u}{\partial x}\bigg|_{x=0}=0$

$\frac{\partial u}{\partial x}\bigg|_{x=L}=0$

with initial conditions: $u(x,t_0) = u_0$

where $[0;L]$ denotes the spatial domain over which I want to solve my problem. Every single layer has its own length, let's define with $L_1,L_2$ and $L_3$ and $L = \sum_i L_i$.

I'm using the MOL approach in order to obtain a continous time-discrete space formulation.

Moreover let's assume that the source term $S(x,t)$ is not zero within the layers 1 and 3, but it's null within the second layer.

My current workflow has been:

1 - Divide every signel layer into $N_i$ points with the FDM method

2 - Treat the interfaces points as "common" points; therefore the total number of discretized points will be $N_1+N_2+N_3-2$

3 - At the interfaces points ($N_1$ and $N_1+N_2-1$) enforce the interface BCs. In $N_1$ I am treating that point as part of the layer 1, while in $N_1+N_2-1$ I'm treating that point as part of the layer 2 (is that correct?)

4 - At the end points I'm enforcing the null-flux conditions.

I am able to solve the equation, but I'm experiencing an increase of the overall mass in the system. As far as I know this is a common problem when dealing with conservation of mass problems and FDM; theoretically I could adopt the FVM method in order to avoid this behaviour, but I was courious to understand if I was doing right in this way or wrong. Understand if it was up to the numerical approach or up to some kind of mine bug.

Moreover: the term $S(x,t)$ I'm assuming to be non-zero in $[1;N_1]$ and $[N_1+N_2;N_1+N_2+N_3-2]$

Thanks

• You can always evaluate your solution at the boundary / interfaces and check your conditions if they hold. And of course check the dgl itself with the solution. – Bort Mar 19 '15 at 13:27
• Hi, thanks for your reply! How can I check a posteriori the boundary/interface conditions? I'm using the built-in solvers like ode15i/s to solve these equations and therefore I was implementing everything in the fnction specs. Is there a way to separate the domains and check the conditions as you mentioned or do I have to define a "custom" sover for this? Moreover: is that important to consider the interface points as parto only of one domain? – FancyPants Mar 29 '15 at 13:04
• Use the method of manufactured solutions to do code verification, you'll find the bugs if there are any. – chris Mar 31 '15 at 6:54
• As mentioned in other works I'm just enforcing the boundary interface conditions using different discretization scheme. Let's say that I'm using some backward scheme once I am in the j-th point (looking from the left domain) and I use a forward approximation in the same point but now as if I was in the right domain. Summing the two equations leads to obtain a single equation for that point. I'm not sure that it is correct . – FancyPants Mar 31 '15 at 9:39