Hi, I have a PDE
\begin{align} \frac{\partial P}{\partial t}(L,x,y) = \nabla\cdot \mathbf{a}^\star\nabla + \mathbf{c}^\star\cdot\nabla P(L,x,y) \end{align}
with Neumann boundary conditions and Dirac delta initial condition P(t=0,x,y) = \delta(x)\delta(y-1).
Approximating the time derivative we have
\begin{align}
\frac{P^{n+1} - P^n}{\Delta t} \approx (\nabla\cdot \mathbf{a}^\star\nabla + \mathbf{c}^\star\cdot\nabla + Q)P^{n+1}.
\end{align}
where Q = A_3^2x/(2y^2) - A_3^2 - A_4 and the A terms are order 1, also
\begin{align}
\mathbf{a}^\star = \begin{bmatrix}
2A_3^2xy^2 & A_3^2xy \\
A_3^2xy & \frac{A_3^2}{2}x
\end{bmatrix},\quad
\mathbf{c}^\star = \begin{bmatrix}
-A_3^2x - 2x A_4 \\
A_3^2y - \frac{A_3^2}{2}\bigg(2y + \frac{x}{y}\bigg) + A_4y
\end{bmatrix}.
\end{align}
Multiply by test function v(\mathbf{x}) and integrate over our domain \Omega to obtain
\begin{align}\nonumber
&\int_{\Omega} \bigg(v(\mathbf{x})P^{n+1} - \Delta L\bigg(v(\mathbf{x})\nabla\cdot\mathbf{a}^\star\nabla P^{n+1} + v(\mathbf{x})\mathbf{c}^\star\cdot\nabla P^{n+1} + v(\mathbf{x})QP^{n+1}\bigg)\bigg)d\mathbf{x} \\
&= \int_{\Omega}v(\mathbf{x})P^{n}d\mathbf{x}.
\end{align}
To formulate the weak form of the problem I integrate by parts the second term to obtain
\begin{align}
\int_{\Omega} v(\mathbf{x})\nabla\cdot\mathbf{a}^\star\nabla P d\mathbf{x} = \int_{\Omega} v(\mathbf{x})\mathbf{n}\cdot \mathbf{a}^\star\nabla P ds - \int_{\Omega} \nabla v(\mathbf{x})\cdot \mathbf{a}^\star\nabla P d\mathbf{x},
\end{align}
with Neumann boundary conditions the surface integral vanishes to obtain
\begin{align}
\int_{\Omega} v(\mathbf{x})\nabla\cdot\mathbf{a}^\star\nabla P d\mathbf{x} = - \int_{\Omega} \nabla v(\mathbf{x})\cdot \mathbf{a}^\star\nabla P d\mathbf{x}.
\end{align}
Hence weak form is written as
\begin{align}\nonumber
&\int_{\Omega}\bigg( v(\mathbf{x})P^{n+1} + \Delta t\bigg(\nabla v(\mathbf{x})\cdot \mathbf{a}^\star\nabla P^{n+1} - v(\mathbf{x})\mathbf{c}^\star\cdot\nabla P^{n+1} - v(\mathbf{x})QP^{n+1}\bigg)\bigg)d\mathbf{x} \\
&= \int_{\Omega}v(\mathbf{x})P^{n}d\mathbf{x},
\end{align}
where I define
\begin{align}
\mathbf{S} = \nabla v(\mathbf{x})\cdot \mathbf{a}^\star\nabla P^{n+1} - v(\mathbf{x})\mathbf{c}^\star\cdot\nabla P^{n+1} - v(\mathbf{x})QP^{n+1},
\end{align}
My problem is that when I integrate over time, the solution P decays, I.e. the initial condition is not preserved. Is there a way to change this formulation to conserve the initial Dirac delta condition?