Help with writing variational form in python

1

Hi, I have a problem with trying to get my variational problem to work correctly. I have

\partial_{L} P(L, \mathbf{x})=(\nabla \cdot \mathbf{a} \nabla+\mathbf{c} \cdot \nabla) P(L, \mathbf{x})

with

\mathbf{a}=\left[\begin{array}{ccc}\frac{A_{5}^{2} z^{2}}{2} & 0 & \frac{A_{5}^{2} x z}{2} \\ 0 & \frac{1}{2}\left(A_{1}^{2}+A_{2}^{2}+A_{3}^{2} \frac{z^{2}}{x^{2}}\right) & 0 \\ \frac{A_{3}^{2} x z}{2} & 0 & \frac{A_{2}^{2} x^{2}}{2}\end{array}\right], \quad \mathbf{c} = \begin{bmatrix}-A_4x \\ 0 \\ A_4z\end{bmatrix}.

My variational problem is then

\frac{P^{n+1}-P^{n}}{\Delta L} \approx(\nabla \cdot \mathbf{a} \nabla+\mathbf{c} \cdot \nabla) P^{n+1},

so

\int_{\Omega}\left(v(\mathrm{x}) P^{n+1}+\Delta L \mathrm{S}\right) d \mathrm{x}=\int_{\Omega} v(\mathrm{x}) P^{n} d \mathrm{x},

with

\mathrm{S}=\nabla v(\mathrm{x}) \cdot \mathrm{a} \nabla P^{n+1}-v(\mathrm{x}) \mathrm{c} \cdot \nabla P^{n+1}.

Furthermore

\begin{aligned} a(P, v) &=\int_{\Omega}(v(\mathbf{x}) P+\Delta L \mathbf{S}) d \mathbf{x} \\ L_{n+1}(v) &=\int_{\Omega} v(\mathbf{x}) P^{n} d \mathbf{x} \end{aligned}.

I try to formulate the problem in python as

a = np.array([[0.5*A_3**2*x[2]**2,0,0.5*A_3**2*x[0]*x[2]],[0,0.5*(A_1**2+A_2**2+A_3**2*(x[2]/x[0])**2),0],[0.5*A_3**2*x[0]*x[2],0,0.5*A_3**2*x[0]**2]])
cee = np.array([-0.5*A_3**2*x[0]+(A_3**2/2-A_4)*x[0], -A_5, -0.5*A_3**2*x[2]+(A_3**2/2+A_4)*x[2]])
a_grad_u = np.multiply(a,grad(u))
c_dot_grad_u = cee.dot(grad(u))
S = np.dot(grad(v),a_grad_u) - np.multiply(v,c_dot_grad_u)
a_form = (u*v + dl*(S))*dx
L_form = u_n*v*dx

Next I have

L = .2   #Final L
num_steps = 10  #Number of length steps
dl = L / num_steps #Length step size

then I solve

for n in range(num_steps):

#Update current lengthj
    l += dl
    

    #Solve
    solve(a_form == L_form, u)

#Update previous solution
    u_n.assign(u)

but I get the error

line 227, in solve
return dolfin.la.solver.solve(*args)
TypeError: solve() missing 1 required positional argument: ‘b’

Have I input my \mathbf{S} incorrectly? I feel like I’m making a really dumb mistake.

2

Hi,
you should post the whole code with proper formatting so that people can help.
At first glance, you should not use Numpy dot or multiply functions but UFL dot or inner.

1 Like
3

Apologies. Here is my code

from __future__ import print_function
from fenics import *
import numpy as np

L = .2   #Final L
num_steps = 10  #Number of length steps
dl = L / num_steps #Length step size

#Define positive constants
A_1 = 6
A_2 = 7
A_3 = 5.4
A_4 = 6
A_5 = 7

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

#End mesh values
x1 = 3
y1 = 30
z1 = 3

#Number of cells in each direction of the mesh
nxyz = 20

#Define mesh
mesh = BoxMesh(Point(x0, y0, z0),Point( x1, y1, z1), nxyz, nxyz, nxyz)

V = FunctionSpace(mesh, 'P', 1)
x = SpatialCoordinate(mesh)
mpx=np.array(mesh.coordinates())
dx = Measure('dx',mesh)

#Dirac Delta (mass) initial condition with epsilon definition
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)

#Bilinear form
a = np.array([[0.5*A_3**2*x[2]**2,0,0.5*A_3**2*x[0]*x[2]],[0,0.5*(A_1**2+A_2**2+A_3**2*(x[2]/x[0])**2),0],[0.5*A_3**2*x[0]*x[2],0,0.5*A_3**2*x[0]**2]])
cee = np.array([-0.5*A_3**2*x[0]+(A_3**2/2-A_4)*x[0], -A_5, -0.5*A_3**2*x[2]+(A_3**2/2+A_4)*x[2]])
a_grad_u = np.multiply(a,grad(u))
c_dot_grad_u = cee.dot(grad(u))

S = dot(grad(v),a_grad_u) - np.multiply(v,c_dot_grad_u)

#Form
a_form = (u*v + dl*(S))*dx
L_form = u_n*v*dx

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

for n in range(num_steps):
    print("We have",num_steps - n, "steps remaining." )
    

    #Update current lengthj
    l += dl
    

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

  

    #Update previous solution
    u_n.assign(u)
4

I may say, originally I had my \mathbf{S} written as

a = np.array([[0.5*A_3**2*x[2]**2,0,0.5*A_3**2*x[0]*x[2]],[0,0.5*(A_1**2+A_2**2+A_3**2*(x[2]/x[0])**2),0],[0.5*A_3**2*x[0]*x[2],0,0.5*A_3**2*x[0]**2]])
cee = np.array([-0.5*A_3**2*x[2]+(A_3**2/2-A_4)*x[0], -A_5, -0.5*A_3**2*x[0]+(A_3**2/2+A_4)*x[2]])
a_dot_grad_u = a.dot(grad(u))
c_dot_grad_u = cee.dot(grad(u))
S = np.dot(grad(v),a_dot_grad_u) - np.dot(v,c_dot_grad_u)

so

a_form = (u*v + dl*(S))*dx
L_form = u_n*v*dx

and the code ran. I’m not sure it’s correct however since I use dot instead of multiply?

5

I think my formulation is fine, no? Since \nabla P is a scalar, i can just write dot(a,grad(u)) where u=P?

6

In your code, just change to

a_grad_u = dot(as_matrix(a),grad(u))