Hello all. I am experimenting with a solution of the Fokker Planck equation of the following first-order form (arising from a noisy carrier-tracking problem in radio science):
\begin{eqnarray}
\frac{\partial p}{\partial t} &=& p + \mathbf{a}\cdot \mathbf{J} - \nabla\cdot (\overline{\mathbf{B}}\cdot\mathbf{J}) \\
\mathbf{J} &=& -\nabla p ,\end{eqnarray}
where \mathbf{J} is a three component probability current, p is the joint probability function in x, y and z and
\mathbf{a} = \left\lbrace\begin{array}{c} (y - x + \sin z)\alpha \\
(x-y)\beta \\
Kx - \omega_0\end{array}\right\rbrace,
with
\overline{\mathbf{B}} = \left\lbrack\begin{array}{ccc}
\alpha^2 N_0 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0\end{array}\right \rbrack.
The solution domain \Omega is the rectangular parallelpiped volume (-0.5, -0.5, -\pi) to (0.5, 0.5, \pi). Neumann boundary conditions are enforced on x = -0.5, 0.5 and y=-0.5, 0.5 surfaces and an absorbing boundary p = 0 is assumed on z = -\pi, \pi surfaces. Time discretization is using the implicit Euler method:
p^n = p^{n+1} -\Delta t (p^{n+1} - \mathbf{a}\cdot \mathbf{J}^{n+1} + \nabla\cdot (\overline{\mathbf{B}}\cdot\mathbf{J}^{n+1})).
For solving this, I am using the least-squares finite element approach (similar to what I did here). The residual is constructed as follows:
R = \vert\!\vert p^{n+1}-p^n-\Delta t (p^{n+1} - \mathbf{a}\cdot\mathbf{J}^{n+1}+\nabla\cdot(\overline{\mathbf{B}}\cdot\mathbf{J}^{n+1}))\vert\!\vert^2_\Omega +
\qquad\qquad\vert\!\vert\mathbf{J}^{n+1}+\nabla p^{n+1}\vert\!\vert^2_\Omega + \vert\!\vert\mathbf{n\cdot J}^{n+1}\vert\!\vert^2_{\partial\Omega_N} + \vert\!\vert p^{n+1}\vert\!\vert^2_{\partial\Omega_A}.
Minimization of the residual is achieved using the usual Ritz process on the unknown variables (p^{n+1}, \mathbf{J}^{n+1}). The initial conditions are set in p^n.
I am doing something incorrectly, most likely, in setting up the initial conditions or the way I do the multiplication of the coefficients. The code gets as far as the assembly process and dies as the result of an “arity” error. I’ve had a look at other posts that describe similar errors, but have had no luck in making my code work. The error output is:
Generating mesh with CGAL 3D mesh generator
<Mesh of topological dimension 3 (tetrahedra) with 7421 vertices and 37627 cells, ordered>
Solution vector size = 221284
Calling FFC just-in-time (JIT) compiler, this may take some time.
Traceback (most recent call last):
File "FP1.py", line 103, in <module>
A, b = assemble_system(a, L, bcs)
File "/usr/lib/python3/dist-packages/dolfin/fem/assembling.py", line 360, in assemble_system
A_dolfin_form = _create_dolfin_form(A_form, form_compiler_parameters)
File "/usr/lib/python3/dist-packages/dolfin/fem/assembling.py", line 58, in _create_dolfin_form
function_spaces=function_spaces)
File "/usr/lib/python3/dist-packages/dolfin/fem/form.py", line 44, in __init__
mpi_comm=mesh.mpi_comm())
File "/usr/lib/python3/dist-packages/dolfin/jit/jit.py", line 47, in mpi_jit
return local_jit(*args, **kwargs)
File "/usr/lib/python3/dist-packages/dolfin/jit/jit.py", line 97, in ffc_jit
return ffc.jit(ufl_form, parameters=p)
File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 217, in jit
module = jit_build(ufl_object, module_name, parameters)
File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 133, in jit_build
generate=jit_generate)
File "/usr/lib/python3/dist-packages/dijitso/jit.py", line 165, in jit
header, source, dependencies = generate(jitable, name, signature, params["generator"])
File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 66, in jit_generate
prefix=module_name, parameters=parameters, jit=True)
File "/usr/lib/python3/dist-packages/ffc/compiler.py", line 143, in compile_form
prefix, parameters, jit)
File "/usr/lib/python3/dist-packages/ffc/compiler.py", line 185, in compile_ufl_objects
analysis = analyze_ufl_objects(ufl_objects, kind, parameters)
File "/usr/lib/python3/dist-packages/ffc/analysis.py", line 90, in analyze_ufl_objects
for form in forms)
File "/usr/lib/python3/dist-packages/ffc/analysis.py", line 90, in <genexpr>
for form in forms)
File "/usr/lib/python3/dist-packages/ffc/analysis.py", line 174, in _analyze_form
do_apply_restrictions=True)
File "/usr/lib/python3/dist-packages/ufl/algorithms/compute_form_data.py", line 418, in compute_form_data
check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode) # Currently testing how fast this is
File "/usr/lib/python3/dist-packages/ufl/algorithms/check_arities.py", line 177, in check_form_arity
check_integrand_arity(itg.integrand(), arguments, complex_mode)
File "/usr/lib/python3/dist-packages/ufl/algorithms/check_arities.py", line 159, in check_integrand_arity
arg_tuples = map_expr_dag(rules, expr, compress=False)
File "/usr/lib/python3/dist-packages/ufl/corealg/map_dag.py", line 37, in map_expr_dag
result, = map_expr_dags(function, [expression], compress=compress)
File "/usr/lib/python3/dist-packages/ufl/corealg/map_dag.py", line 86, in map_expr_dags
r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
File "/usr/lib/python3/dist-packages/ufl/algorithms/check_arities.py", line 48, in sum
raise ArityMismatch("Adding expressions with non-matching form arguments {0} vs {1}.".format(_afmt(a), _afmt(b)))
ufl.algorithms.check_arities.ArityMismatch: Adding expressions with non-matching form arguments () vs ('v_0',).
My model code is:
from dolfin import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt
import os, sys, traceback
pi = 3.1415926536
tol = 1.0e-12
# Problem parameters
alfa = 1.0
beta = 0.1
N0 = 0.1
K = 1.0
omega = 1.0
sig = 0.3
# time stepping params
dt = 0.01
NumSteps = 100
class Wall(SubDomain):
def inside(self, x, on_boundary):
return on_boundary
class AbsorbingBC(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (near(x[2], -pi, tol) or near(x[2], pi, tol))
box = Box(dolfin.Point(-0.5,-0.5, -pi), dolfin.Point(0.5, 0.5, pi))
mesh = generate_mesh(box, 70)
info(mesh)
# Mark boundaries
sub_domains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
sub_domains.set_all(2)
wall = Wall()
wall.mark(sub_domains, 0)
absorb = AbsorbingBC()
absorb.mark(sub_domains, 1)
File("BoxSubDomains.pvd").write(sub_domains)
# PDE Coefficients , C1 is vector, C2 is 3x3 matrix
C1 = Expression(('(x[1] - x[0] + sin(x[3])) * alfa', '(x[0] - x[1]) * beta', 'K * x[0] - omega'), degree = 3, alfa = alfa, beta = beta, K = K, omega = omega)
C2 = Constant(((N0, 0.0, 0.0), (0.0, 0.0, 0.0), (0.0, 0.0, 0.0)))
# Set up function spaces
cell = tetrahedron
ele = FiniteElement('Lagrange', cell, 2)
vec_ele = VectorElement('Lagrange', cell, 2)
V4 = FunctionSpace(mesh, MixedElement([ele, vec_ele]))
V = FunctionSpace(mesh, ele)
u1 = Function(V4)
un = Function(V4)
pv, Jv = TestFunctions(V4)
#surface integral definitions from boundaries
ds = Measure('ds', domain = mesh, subdomain_data = sub_domains)
# set initial conditions
g = Expression(('exp(-(x[0]*x[0]+x[1]*x[1]+x[2]*x[2])/(4.0 * pi * pi)) / (2.0 * sig * sig)', '0.0', '0.0', '0.0'), degree = 3, sig = sig)
pp0 = interpolate(g, V4)
un.assign(pp0)
#Boundary condition dictionary
boundary_conditions = {0: {'Wall' : 0.0},
1: {'AbsorbingBC': 0.0}}
n = FacetNormal(mesh)
pr, Jr = split(u1) #Split for forms
pn, Jn = split(un)
#Build boundary conditions
bcs = []
integrals_load =[]
#Neumann BC on impermeable walls
for i in boundary_conditions:
if 'Wall' in boundary_conditions[i]:
bb1 = inner(dot(n, Jv), dot(n, Jr)) * ds(i)
integrals_load.append(bb1)
# Build absorbing BC
for i in boundary_conditions:
if 'AbsorbingBC' in boundary_conditions[i]:
bb = inner(pv, pr) * ds(i)
integrals_load.append(bb)
a = (inner(pv - dt * (pv + dot(C1, Jv) - div(C2 * Jv)), pr - dt * (pr + dot(C1, Jr) - div(C2 * Jv))) + inner(Jv + grad(pv), Jr + grad(pr))) * dx + sum(integrals_load)
L = inner(pv - dt * (pv + dot(C1, Jv) - div(C2 * Jv)), pn) * dx
vdim = u1.vector().size()
print ('Solution vector size = ',vdim)
u1.vector().zero()
A, b = assemble_system(a, L, bcs) #XXXXX Seems to die here XXXX
solver = PETScKrylovSolver("cg","sor")
solver.parameters['absolute_tolerance'] = 1e-5
solver.parameters['relative_tolerance'] = 1e-12
solver.parameters['maximum_iterations'] = 10000
solver.parameters['monitor_convergence'] = True
solver.parameters['error_on_nonconvergence'] = False
solver.parameters['nonzero_initial_guess'] = True
solver.parameters['report'] = True
solver.set_operators(A, A)
fpp = File('Prob.pvd')
fpJ = File('Curr.pvd')
p, J = u1.split()
# Do time stepping
t = 0
for nn in range(NumSteps):
t += dt
#Solve
solver.solve(u1.vector(),b)
#Update solution
un.assign(u1)
fpp << (p, t)
fpJ << (J, t)
sys.exit(0)Any hints as to where I am going wrong would be very much appreciated! I am using Fenics 2019 on an Ubuntu 18 PC where Fenics was installed from PPA. This version of Fenics has worked well with other problems, so the installation seems OK.