Writing equation in conservative flux form

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?