Problem implementing 4d PDE

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.

You are working on a 1D space with \texttt{IntervalMesh} thus \texttt{x[2]} and \texttt{x[1]} ain’t defined. The same goes with \texttt{UnitSquare} try to create if you can a 3D mesh knowing its limits.

The following documentation might help you if your domain is a polygon

https://fenicsproject.org/docs/dolfin/1.4.0/python/demo/documented/mesh-generation/python/documentation.html

1 Like

I tried using a boxmesh since I know that x,y \in[0,10] and z\in [1,10] via

#Intial mesh values
x0 = 0
y0 = 0
z0 = 1

#End mesh values
x1 = 10
y1 = 10
z1 = 10 

#Number of cells in each direction of the mesh
nx = 100
ny = 100
nz = 100

#Using 3d boxmesh
mesh = BoxMesh(x0, y0, z0, x1, y1, z1, nx, ny, nz)
plot(mesh)
#mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, 'P', 1)

x = SpatialCoordinate(mesh)

Is this along the right lines? Still error codes.

You have to declare point and not just values.
Try the following code

mesh = BoxMesh(Point(x0, y0, z0),Point( x1, y1, z1), nx, ny, nz)

Thanks. It runs now but I encounter a memory error

Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Solving linear variational problem.

UMFPACK V5.7.8 (Nov 9, 2018): ERROR: out of memory

Traceback (most recent call last):
File “abrahams_rot.py”, line 85, in
solve(a_form == L_form, u)
File “/opt/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/dolfin/fem/solving.py”, line 220, in solve
_solve_varproblem(*args, **kwargs)
File “/opt/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/dolfin/fem/solving.py”, line 247, in _solve_varproblem
solver.solve()
RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at


*** fenics-support@googlegroups.com


*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.


*** -------------------------------------------------------------------------
*** Error: Unable to successfully call PETSc function ‘KSPSolve’.
*** Reason: PETSc error code is: 76 (Error in external library).
*** Where: This error was encountered inside /usr/local/miniconda/conda-bld/fenics-pkgs_1582286161246/work/dolfin/dolfin/la/PETScKrylovSolver.cpp.
*** Process: 0


*** DOLFIN version: 2019.1.0
*** Git changeset:
*** -------------------------------------------------------------------------

Okay, if I reduce the number of cells in each direction it actually compiles. Problem is now how to plot this solution. Say, how to plot P(L=.5,x=1,y=1,z) on the vertical axis and z on the horizontal axis? When I try plot(u) in the loop I get the error

NotImplementedError: It is not currently possible to manually set the aspect on 3D axes

Try using paraview for the advanced visualization. To export your solution (which I assume has the name u) you can run the following commands, in the for loop

u_file = File('u.pvd')
u_file << u, l

which generates the file u.pvd that can be opened with Paraview. There, the timesteps correspond to your first coordinate, i.e., L

Hi, I tried this. It simply gives me a cuboid with constant values lol. Not sure what to make of it.

Maybe you can try to project on the subdomain created by the line x=1,y=1

Try using different filters, e.g., the slice filter.
From the outside you can only see the boundary values…

Yes, I realised that quickly. It’s the same story. It’s constant throughout…

You mean, plot P(L_{fixed},x_{fixed},z) on the vertical axis and z on the horizontal axis?

If you want. But first the most complicated part is to project … with Fenics it’s may be possible by creating the subdomain wanted (the line x=1,y=1 obviously L is fixed when you have a solution).

But what I have done for some of my previous works is to select the values of P I wanted to plot with numpy arrays.

import numpy as np
#x are the cordinates on the mesh
n=len(P.vector())
P_arr=np.array(P.vector())
x_arr=np.reshape(np.array(x.vector()),n//3,3)

eps=1e-3
P_line=P_arr[np.where((np.abs(x_arr[:,0]-1)<eps) & (np.abs(x_arr[:,1]-1)<eps))]
z_line=x_arr[np.where((np.abs(x_arr[:,0]-1)<eps) & (np.abs(x_arr[:,1]-1)<eps)),2]
import matplotlib.pyplot as plt
plt.scatter(z_line,P_line)

Try with different values of eps

Hi, I tried this. I understood this your “x” as being my “mpx” term

mesh = BoxMesh(Point(x0, y0, z0),Point( x1, y1, z1), nx, ny, nz)
V = FunctionSpace(mesh, 'P', 1)
x = SpatialCoordinate(mesh)
mpx=mesh.coordinates()

so that

m=len(u.vector())
P_arr=np.array(u.vector())
x_arr=np.reshape(np.array(mpx.vector()),m//3,3)

eps=1e-3
P_line=P_arr[np.where((np.abs(x_arr[:,0]-1)<eps) & (np.abs(x_arr[:,1]-1)<eps))]
z_line=x_arr[np.where((np.abs(x_arr[:,0]-1)<eps) & (np.abs(x_arr[:,1]-1)<eps)),2]

plt.scatter(z_line,P_line)

but I get the error

x_arr=np.reshape(np.array(mpx.vector()),m//3,3)

AttributeError: ‘numpy.ndarray’ object has no attribute ‘vector’

My x was a vector function similar to u.
Here it seems that mpx is already an array and I’m not sure how this array is organized can you print the begining of it.
If it has the same organization as my x then just remove the .vector()
Otherwise:

V2=VectorFunctionSpace(mesh,'P',1)
X=Expression(('x[0]','x[1]','x[2]'),element = V2.ufl_element())
x=interpolate(X,V2)
x_arr=np.reshape(np.array(x.vector()),m//3,3)

mpx is an array

(1331, 3)
[[ 0.  0.  1.]
 [ 2.  0.  1.]
 [ 4.  0.  1.]
 ...
 [16. 20. 20.]
 [18. 20. 20.]
 [20. 20. 20.]]

with shape (1331, 3)

I try

m=len(u.vector())
print(m//3)
P_arr=np.array(u)
x_arr=np.reshape(np.array(mpx),m//3,3)

but I get the error

ValueError: Non-string object detected for the array ordering. Please pass in ‘C’, ‘F’, ‘A’, or ‘K’ instead…

It seems that

x_arr=mpx 

I hope that it is the same order for mpx and P

Hi, this works now. However I get an empty array for z_line and P_line has too many indices for the array. Perhaps a problem with my solution to this PDE…

It’s weird indeed… did you use the method I proposed ?

Yes. I will update my solution further once I have formatted it correctly. I am still having some trouble with trusting this solution…