I’m trying to solve
\begin{align}
\frac{\partial P}{\partial L}(L, x, y, z)&=\frac{A_{3}^{2} z^{2}}{2} \frac{\partial^{2} P}{\partial x^{2}}(L, x, y, z)
\\ &+A_{3}^{2} x z \frac{\partial^{2} P}{\partial x \partial z}(L, x, y, z)+\left(A_{1}^{2}+A_{2}^{2}+A_{3}^{2} \frac{z^{2}}{x^{2}}\right) \frac{\partial^{2} P}{\partial y^{2}}(L, x, y, z)
\\ &+A_{3}^{2} x^{2} \frac{\partial^{2} P}{\partial z^{2}}(L, x, y, z) \\ &+\left(\frac{A_{3}^{2}}{2}-A_{4}\right) x \frac{\partial P}{\partial x}(L, x, y, z)
\\ &+\left(\frac{A_{3}^{2}}{2}+A_{4}\right) z \frac{\partial P}{\partial z}(L, x, y, z) \\ &-A_{5} \frac{\partial P}{\partial y}(L, x, y, z).
\end{align}
The initial condition is P(L=0,x,y,z) = \delta(x)\delta(y)\delta(z-1) where \delta denotes a Dirac delta function (Dirac mass) at L=0. We also have pure Neumann boundary conditions. I think I have correctly obtained the variational form
\begin{align}\nonumber
a(P,v) &= \int_{\Omega}\bigg[Pv + \Delta L\bigg(\frac{A_3^2z^2}{2}\frac{\partial v}{\partial x}\frac{\partial P}{\partial x} + \bigg(A_1^2 + A_2^2 + A_3^2\frac{z^2}{x^2}\bigg)\frac{\partial v}{\partial y}\frac{\partial P}{\partial y} + A_3^2x^2\frac{\partial v}{\partial x}\frac{\partial P}{\partial z} \\ \nonumber
&+ A_3^2\bigg(zv + xz\frac{\partial v}{\partial x}\bigg)\frac{\partial P}{\partial z} - \bigg(\frac{A_3^2}{2} - A_4\bigg)x\frac{\partial P}{\partial x}v \\
&-\bigg(\frac{A_3^2}{2} + A_4\bigg)z\frac{\partial P}{\partial z}v + A_5\frac{\partial P}{\partial y}v\bigg)\bigg]\,d\mathbf x, \\
L_{n+1}(v) &= \int_{\Omega}P^{n}v\,d\mathbf x,&
\end{align}
where \mathbf x = [x,y,z]^{T}. I have tried implementing this via
from __future__ import print_function
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate
L = 1 # final L
num_steps = 10 # number of length steps
dl = L / num_steps # length step size
theta = 1 #Lloc constant term
#Define material constants
A_1 = 1
A_2 = 1
A_3 = 1
A_4 = 1
# Create mesh and define function space
nx = 1000
mesh = IntervalMesh(nx, 1., 10)
#mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, 'P', 1)
x = SpatialCoordinate(mesh)
# Dirac-delta initial condition
u_D = Expression("(1./(pow(x[0],2) + pow(eps, 2)))*(1./(pow(x[1],2) + pow(eps, 2)))*(1./(pow(x[2]-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)
a_form = (u*v + dl*((A_3**2/2)*x[2]**2*v.dx(0)*u.dx(0)
+ (A_1**2 + A_2**2 + A_3**2*(x[2]**2/x[0]**2))*v.dx(1)*u.dx(1)
+ A_3**2*x[0]**2*v.dx(0)*u.dx(2)
+ A_3**2*(x[2]*v +x[0]*x[2]*v.dx(0))*u.dx(2)
- (A_3**2/2 + A_4)*x[2]*u.dx(2)*v
+ A_5*u.dx(1)*v
))*dx
L_form = u_n*v*dx
# Time-stepping
u = Function(V)
l = 0
#List for values of L
L_vals = []
for n in range(num_steps):
# Update current time
l += dl
# Compute solution
solve(a_form == L_form, u)
vtkfile = File('ab_sol.pvd')
vtkfile << u
#Plot solution
plot(u)
# Update previous solution
u_n.assign(u)
#Print area under pdf at each length step
print(assemble(u*dx))
#Plotting PDF solution, comment out for viewing mean moment integral solutions
plt.xlabel('$?$')
plt.ylabel('$P$')
plt.show()
but I get the error code
Assertion failed: (index >= 0 && index < size()), function operator[], file /opt/anaconda3/envs/fenicsproject/include/eigen3/Eigen/src/Core/DenseCoeffsBase.h, line 162. Abort trap: 6
Would appreciate some tips. I’m not used to solving PDE’s with these higher dimensional terms. Thanks in advance.