I’m struggling to find resources (other than some things using Lagrange multipliers) to enforce a normalisation condition to the solution of a PDE. For PDE’s where the solution is a probability density function
\int_{-\infty}^{\infty}P(x,t)\,dx = 1
must hold. How would one implement this with use in FeniCS? So far I have only been able to ensure this holds for my problem by choosing a slow diffusion coefficient in the PDE so that over a large time the integral is preserved to 1. This becomes problematic when one wants to investigate the diffusion coefficient. Any help/comments appreciated.
Using Galerkin’s method with zero-flux BCs should be globally-conservative. Using the notation of the reference linked in your post, the strong form of the PDE is
\partial_tP + \nabla\cdot\mathbf{S} = 0\text{ ,}
where I’ve generalized slightly to multiple dimensions and a vector-valued flux. If one takes the weak form
\int_\Omega(\partial_tP) Q\,d\Omega - \int_\Omega\mathbf{S}\cdot\nabla Q\,d\Omega = 0~~~\forall Q\text{ ,}
the absence of a boundary term from the integration by parts weakly enforces the BC \mathbf{S}\cdot\mathbf{n} = 0 on \partial\Omega, where \mathbf{n} is the outward-facing surface normal to \partial\Omega. As long as the test space that Q is in contains a constant function over \Omega (i.e., Q\equiv C\in\mathbb{R}), which is true in most standard finite element spaces (when no Dirichlet BCs are present), then the weak form holds for Q such that \nabla Q = \mathbf{0}, so
\int_\Omega\partial_tP\,d\Omega = 0\text{ ,}
i.e., the integral of P remains constant at its initial value. This would remain true (exactly) with standard time integrators, like Euler’s method, midpoint rule, etc.
From your other posts, it looks like you’ve expanded \mathbf{S} into many terms and the system is no longer clearly in this simple conservation form. I would recommend trying to offload as much algebra as possible to FEniCS’s UFL, to avoid mistakes and keep the abstract structure of the problem clearly visible.
2 Likes
Thanks David. Appreciate your help I will work on this.
Hi David I thought some more about this. l think I understand that discretising the time derivative on the left hand side causes the loss of conservation? I didn’t think that the weak form you have written could be handled without approximating the derivative like I have done previously?
Discretizing the time derivative should not hurt conservation. For example, with Euler’s method, the choice of Q \equiv 1 gives
\int_\Omega \frac{1}{\Delta t}\left(P^{n+1} - P^n\right)\,d\Omega = 0
which simplifies to
\int_\Omega P^{n+1}\,d\Omega = \int_\Omega P^n\,d\Omega\text{ .}
1 Like
It looks like this is the approach I have used in my problem in recent posts? Isn’t this the standard way of finding the bilinear form and approximating the derivative. If that’s true my problem is the same but with the exception that my \mathbf{S} is more complicated. The code
from __future__ import print_function
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
L = 20 #Final L
num_steps = 1000 #Number of length steps
dl = L / num_steps #Length step size
#Define positive constants
A_1 = .009
A_2 = .008
A_3 = .0096
A_4 = .002
A_5 = .001
#Intial mesh values
x0 = 0.0001
y0 = 0.0001
z0 = 1
#End mesh values
x1 = 2000
y1 = 2000
z1 = 2000
#Number of cells in each direction of the mesh
nx = 10
ny = 10
nz = 10
#Define mesh
mesh = BoxMesh(Point(x0, y0, z0),Point( x1, y1, z1), nx, ny, nz)
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)
#u_D = Expression("(7./pow(cosh(7.*(x[0])),2))*(7./pow(cosh(7.*(x[1])),2))*(7./pow(cosh(7.*(x[2]-1.)),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_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
#Starting to build length steps
u = Function(V)
l = 0
#Create lists
L_vals = []
prob_volume = []
for n in range(num_steps):
#Update current lengthj
l += dl
#Solve
solve(a_form == L_form, u)
u_file = File('u.pvd')
u_file << u, l
#Store lengths
L_vals.append(l)
#Store probability volume
prob_volume.append(assemble(u*dx))
#Update previous solution
u_n.assign(u)
#Print probability volume at each length step (should be preserved =1)
print("Area under PDF is:")
print(assemble(u*dx))
#Plot probability volume vs length
plt.xlabel('$L$')
axes = plt.gca()
axes.set_ylim([0,1.1])
plt.ylabel('$\\int_{\\Omega}P(L,x,y,z)d\\mathbf{x}$')
plt.scatter(L_vals,prob_volume)
plt.show()
however when I integrate for larger L (time essentially) the integral drops off
Please tell me if I have completely missed your point!?
It looks like the formulation in the code is not in the conservative form of a flux weighted against the gradient of the test function. For instance, the test function v appears in some terms within the dl*(...) part of a_form without a spatial derivative. For clarity (and to avoid algebra errors or typos), I would try to first have line of the form
S = ...
defining the flux, then define
a_form = (u*v - dl*dot(S,grad(v)))*dx
(What is presumably \nabla\cdot\mathbf{S} in the PDE is expanded into many terms, so it’s not easy for me to pick out what S should be, although it should follow from the physics of the problem; again, I’ve not studied the Fokker–Planck equation myself.)
HI David thanks for your response. I tried this with a simpler PDE. Suppose
\begin{align}\partial_{t}P(x,t) &= \partial_{x}\bigg((x^2-1)\partial_{x}P(x,t)\bigg) \\
&= \nabla\cdot\mathbf{S},
\end{align}
where \mathbf{S} = (x^2-1)\partial_x P. The weak form is then
\int_{\Omega}(\partial_{t}P^{n+1})Q\,d\Omega = \int_{\Omega}\bigg(P^{n}Q-\Delta L \mathbf{S}\partial_{x}Q\, \bigg)d\Omega
so (where u=P and Q=v)
S = ((x[0]**2-1)*u.dx(0))
a_form = (u*v - dl*S*v.dx(0))*dx
L_form = u_n*v*dx
If that is true, I essentially need to repeat this form for a more complicated \mathbf{S}?
Yes, that sounds right, with the nitpick that I’d multiply the scalar (x^2-1)\partial_xP by \mathbf{e}_1 for a vector of dimension one, since FEniCS’s UFL distinguishes between scalars and vectors of dimension one. So, if you defined S = as_vector(((x[0]**2-1)*u.dx(0),)) (where, in general, the argument to as_vector would be a tuple of expressions for its components), you could use the direct notation of dot(S,grad(v)), which would be more concise with more space dimensions. (Also, I’ll mention that, for \vert x\vert > 1, your given flux would become unstable, like heat flowing to higher temperatures.)
1 Like
Fascinating. I’m more confused now. If I don’t use the Galerkin bilinear form of the PDE I wrote above (with Dirac mass initial condition) I obtain
a_form = (u*v + dt*(x[0]*x[0]-1)*u.dx(0)*v.dx(0))*dx
L_form = u_n*v*dx
Integrating this over x\in [1,500) the quantity \int_{1}^{500}P(t,x)dx is conserved at each timestep to be equal to 1 (less than 0.0001% error). Why is this? Because in some problems I consider (the PDE I wrote above) where x\in[1,\infty)?
First, there is a sign difference in the diffusion term of your second formulation, so now it is diffusive for \vert x\vert > 1 and anti-diffusive for \vert x\vert < 1. Second, the formulation can be unstable but conservative; running the heat equation backwards in time on an insulated domain still conserves the total energy, but is highly sensitive to perturbations in the initial condition.
2 Likes
Thank you for your reply. I apologise for my error, that makes sense. Kudos.
