Hi David, thanks for your reply. I have rewritten the PDE as
\frac{\partial P}{\partial L}(L,x,y) = \nabla\cdot \mathbf{a}^\star\nabla + \mathbf{c}^\star\cdot\nabla P(L,x,y) + QP(L,x,y)
where
\begin{align}
\boldsymbol{a}^\star = \begin{bmatrix}
\frac{A_3^2y^2}{2} & \frac{A_3^2xy}{2} \\
\frac{A_3^2xy}{2} & \frac{A_3^2x^2}{2}
\end{bmatrix},
\quad
\boldsymbol{c}^\star = \begin{bmatrix}
-\bigg(\frac{A_3^2y^2}{2x} + A_4x\bigg) - \frac{A_3^2y}{2} \\
-\bigg(\frac{A_3^2x^2}{2y} - A_4y\bigg) - \frac{A_3^2x}{2}
\end{bmatrix},\quad
Q = \frac{A_3^2}{2}(y^2/x^2 + x^2/y^2 - 2) \end{align}
Approximating the length derivative we have
\begin{align}
\frac{P^{n+1} - P^n}{\Delta L} \approx (\nabla\cdot \mathbf{a}^\star\nabla + \mathbf{c}^\star\cdot\nabla)P^{n+1}.
\end{align}
Multiply by test function v(\mathbf{x}) and integrate over our domain \Omega to obtain
\begin{align}
\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}
\int_{\Omega}\bigg( v(\mathbf{x})P^{n+1} + \Delta L\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 we 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}
so
\begin{align}
\int_{\Omega}\bigg( v(\mathbf{x})P^{n+1} + \Delta L\mathbf{S}\bigg)d\mathbf{x} = \int_{\Omega}v(\mathbf{x})P^{n}d\mathbf{x},
\end{align}
so
\begin{align}
a(P,v) = L_{n+1}(v),
\end{align}
where
\begin{align}
a(P,v) &= \int_{\Omega}\bigg( v(\mathbf{x})P + \Delta L\mathbf{S}\bigg)d\mathbf{x}, \\
L_{n+1}(v) &= \int_{\Omega}v(\mathbf{x})P^{n}d\mathbf{x}.&
\end{align}
My code for solving this is
from __future__ import print_function
from fenics import *
import numpy as np
#constants
A_3 = .2
A_4 = .1
#Final L
L = 1
#Number of length steps
num_steps = 10
#Length step size
dl = L / num_steps
#Intial mesh values
x0 = .1
y0 = 1
# #End mesh values
x1 = 10
y1 = 10
#Cell numbers in mesh
nx = 100
ny = 100
#Define mesh
mesh = RectangleMesh(Point(x0,y0), Point(x1,y1), nx, ny, "right")
V = FunctionSpace(mesh, 'P', 1)
x = SpatialCoordinate(mesh)
mpx = np.array(mesh.coordinates())
dx = Measure('dx',mesh)
#Dirac delta intitial condition
u_D = Expression("(1./(pow(x[0],2) + pow(eps, 2)))*(1./(pow(x[1]-1.,2) + pow(eps, 2)))", eps=1e-6, degree=2)
#Define initial value
u_n = interpolate(u_D, V)
C = assemble(u_n*dx)
u_n.assign(u_n/C)
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
#Bilinear form full
a = as_matrix([[.5*A_3**2*x[1]**2, .5*A_3**2*x[0]*x[1]],
[.5*A_3**2*x[0]*x[1], .5*A_3**2*x[0]**2]])
cee = as_vector([-(.5*A_3**2*(x[1]**2/x[0]) + A_4*x[0]) -.5*A_3**2*x[1],
-(.5*A_3**2*(x[0]**2/x[1])) - .5*A_3**2*x[0] ])
#a_grad_u = dot(a,grad(u))
a_grad_u = a*grad(u)
c_dot_grad_u = dot(cee,grad(u))
S = dot(grad(v),a_grad_u) - v*c_dot_grad_u -.5*A_3**2*(x[0]**2/x[1]**2 + x[1]**2/x[0]**2 - 2)*u*v
a_form = (u*v + dl*S)*dx
L_form = u_n*v*dx
#Starting to build length steps
u = Function(V)
l = 0
#Solve
for n in range(num_steps):
#Update current length
l += dl
print('L =',l)
#Solve
solve(a_form == L_form, u)
u_file = File('u_sol_.pvd')
u_file << u, l
#Probability volume
print("Area under PDF is:")
print(assemble(u*dx))
#Compute transmission coefficient at each length step
print("Transmission coefficient is:")
print(assemble((u/x[1]**2)*dx ) )
#Update previous solution
u_n.assign(u)
#Print messages
print("We are at step",n+1,"out of",num_steps , "step(s)." )
I am interested in computing the quantity in my code assemble((u/x[1]**2)*dx )
but this blows up. Any thoughts? Does my formulation look correct?