Using BoxMesh for 2d problem

mesh
1

Hi, I have a PDE \partial_L P(L,\boldsymbol{x}) = (\mathbf{a}\nabla + \mathbf{b})\cdot\nabla P(L,\mathbf{x}),

where

\mathbf{a} = \begin{bmatrix} \frac{A_3^2z^2}{2} & 0 & \frac{A_3^2xz}{2} \\ 0 &\frac12\bigg(A_1^2+A_2^2+A_3^2\frac{z^2}{x^2}\bigg) & 0 \\ \frac{A_3^2xz}{2} & 0 & \frac{A_3^2x^2}{2} \end{bmatrix}, \quad \mathbf{b} = \begin{bmatrix} (\frac{A_3^2}{2}-A_4)x \\ -A_5 \\ (\frac{A_3^2}{2}+A_4)z \end{bmatrix},

with x \in [0,x_\infty), \, y \in [0,y_\infty) and z \in [1,z_\infty) and impose natural Neumann boundary conditions.

The variational problem is then

a(P,v) = L_{n+1}(v),
where
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}, \\ \mathbf{S} = \nabla v(\mathbf{x})\cdot \mathbf{a}\nabla P^{n+1} - v(\mathbf{x})\mathbf{c}\cdot\nabla P^{n+1}.

If I want to remove the y variable from the PDE, can I simply use a 3D boxmesh but remove the parts in the \mathbf{a} and \mathbf{c} which contain y derivates? Code below.

from __future__ import print_function
from fenics import *
import numpy as np
import matplotlib.pyplot as plt

A_1 = 0.6667979098259997

A_2 = -0.15822994796144957

A_3 = -0.3953106541204352

A_4 = 0.10550739857257128

A_5 = 0.1562705132611264

#Final L
L = 3
#Number of length steps
num_steps = 200
#Length step size
dl = L / num_steps

#Intial mesh values
x0 = .01
y0 = .01
z0 = 1

#End mesh values
x1 = 1
y1 = 1
z1 = 4


n_x = 1
n_y = 10
n_z = 500

mesh = BoxMesh(Point(x0, y0, z0),Point(x1, y1, z1), n_x, n_y, n_z)
V = FunctionSpace(mesh, 'P', 1)

x = SpatialCoordinate(mesh)
mpx = np.array(mesh.coordinates())
dx = Measure('dx',mesh)


#Multiple Dirac Delta (mass) initial condition with epsilon definition
#Dirac delta without h variable
u_D = Expression("(7./pow(cosh(7.*(x[0])),2))*(7./pow(cosh(7.*(x[2]-1.)),2))", eps=1e-3, 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
a = as_matrix([[0.5*A_3**2*x[2]**2,0,0.5*A_3**2*x[0]*x[2]],
[0,0,0],
[0.5*A_3**2*x[0]*x[2],0,0.5*A_3**2*x[0]**2]])

cee = as_vector([-A_4*x[0],
0,
A_4*x[2]])

a_grad_u = dot(a,grad(u))
c_dot_grad_u = dot(cee,grad(u))
S = dot(grad(v),a_grad_u) - v*c_dot_grad_u
a_form = (u*v + dl*(S))*dx
L_form = u_n*v*dx

#Starting to build length steps
u = Function(V)
l = 0

#Create empty lists
mean_t = []
std_t = []
L_vals = []
prob_volume = []
u_vals = [] 

for n in range(num_steps):
    #Update current length
    l += dl
    

    #Solve
    solve(a_form == L_form, u)
    
    u_file = File('u.pvd')
    u_file << u, l
    
    #Store lengths
    L_vals.append(l)

    print("Area under PDF is:")
    print(assemble(u*dx))
    #Compute transmission coefficient at each length step
    print("Transmission coefficient is:")
    print(assemble((u/x[2]**2)*dx ) )
    #Store power transmission coefficient in list at each l step
    mean_t.append( (assemble( (u/x[2]**2 )*dx) ) )
    
    #Store std of power transmission at each l
    std_t.append(assemble( (u/x[2]**4 )*dx)  - assemble( (u/x[2]**2 )*dx)**2 )
    #Store probability volume
    prob_volume.append(assemble(u*dx))



    #Update previous solution
    u_n.assign(u)

#Plot statistical moment

fig = plt.figure()
ax1 = fig.add_subplot(111)
#Just E[tau]
plt.suptitle('Mean Power Transmission Versus Length')
plt.xlabel('$L$')
plt.ylabel('$\\mathbf{E}[\\tau(L)]$')
ax1.scatter(L_vals,mean_t,s=2,c='red')
ax1.set_ylim([0,1.1])
plt.savefig('mean_t_L.pdf')
plt.close()