Dear FEniCS community,
I have a diffusion equation for unit vector \mathbf{n}(\mathbf{x}, t)
\frac{\partial}{\partial t} \mathbf{n}(\mathbf{x}, t) = \nabla^2 \mathbf{n} (\mathbf{x}, t), \quad |\mathbf{n}| = 1,
with periodic boundary condition and homogeneous initial condition
\mathbf{n}(\mathbf{x},0) = (0, 1).
The solution should not evolve due to the homogeneous initial condition. Thus the Laplacian of \mathbf{n} should be (0, 0).
However, when implemented using FEniCS, \mathbf{n} is not exactly (0,1) due to numerical round-off error. The laplacian of \mathbf{n}
ufl.nabla_div(ufl.nabla_grad(n))
seems to exaggerate the round-off error of \mathbf{n}, and diverges a lot from (0,0).
What is the right way to calculate Laplacian in this situation?
Please enlighten me on this. Thanks.
Below are the whole code can be run and sample output.
PS: I use fenics-2019.1.0-py39hf3d152e_18 installed via Conda if it helps.
code:
import dolfin as df
import ufl
import numpy as np
num_steps = 500 # number of time steps
dt = 1e-3 # time step size
# set mesh size
nn = 80
nx = nn
ny = nn
class ICn(df.UserExpression):
def eval(self, values, x):
values[0] = 0
values[1] = 1
def value_shape(self):
return (2,)
class PeriodicBoundary(df.SubDomain):
def inside(self, x, on_boundary):
return bool((df.near(x[0], 0) or df.near(x[1], 0)) and
(not ((df.near(x[0], 0) and df.near(x[1], 1)) or
(df.near(x[0], 1) and df.near(x[1], 0)))) and
on_boundary)
def map(self, x, y):
if df.near(x[0], 1) and df.near(x[1], 1):
y[0] = x[0] - 1.
y[1] = x[1] - 1.
elif df.near(x[0], 1):
y[0] = x[0] - 1.
y[1] = x[1]
else: # near(x[1], 1)
y[0] = x[0]
y[1] = x[1] - 1.
def laplacian(n):
"""Calculate laplacian of director."""
return ufl.nabla_div(ufl.nabla_grad(n))
def normalize(w, dim=2, epsilon=0.01):
"""Check and normalize director to norm 1."""
wv = w.vector()[:]
wa = np.reshape(wv, (-1, dim))
norms = np.linalg.norm(wa, axis=1, keepdims=True)
if(np.max(np.abs(norms-1.)) > epsilon):
wa = wa/norms
for i in range(dim):
wv[i::dim] = wa[:, i]
w.vector().set_local(wv)
w.vector().apply('')
return w
# Create mesh and define function spaces
mesh = df.UnitSquareMesh(nx, ny)
# apply periodic boundary condition for director field
pbc = PeriodicBoundary()
W = df.VectorFunctionSpace(mesh, 'P', 3, constrained_domain=pbc)
DG = df.VectorFunctionSpace(mesh, 'DG', 1)
# Define trial and test functions
n = df.TrialFunction(W)
m = df.TestFunction(W)
# Define functions for solutions at previous and current time steps
n_n = df.Function(W)
n_ = df.Function(W)
ntmp = df.Function(W)
n_init = ICn()
n_n.interpolate(n_init)
k = df.Constant(dt)
F0 = ufl.dot((n-n_n)/k, m)*ufl.dx \
+ ufl.inner(ufl.nabla_grad(n_n), ufl.nabla_grad(m))*ufl.dx \
a0 = ufl.lhs(F0)
L0 = ufl.rhs(F0)
# Assemble matrices
A0 = df.assemble(a0)
# Time-stepping
t = 0
for i in range(num_steps):
# Update current time
t += dt
b0 = df.assemble(L0)
df.solve(A0, ntmp.vector(), b0)
# normalize unit director
ntmp = normalize(ntmp) # check and normalize to norm 1
n_.vector()[:] = ntmp.vector()
laplace = df.project(laplacian(n_), DG)
print('------ n ----------------------')
print(n_.vector()[:])
print('------ laplacian of n -----------')
print(laplace.compute_vertex_values())
# Update previous solution
n_n.assign(n_)
Sample output:
------ n ----------------------
[0. 1. 0. ... 1. 0. 1.]
------ laplacian of n -----------
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 ... 1.97993207e-07 1.98205422e-07 -1.97993207e-07]
------ n ----------------------
[0. 1. 0. ... 1. 0. 1.]
------ laplacian of n -----------
[ 0. 0. 0. ... -0.00022009 -0.00022037 0.00031604]
------ n ----------------------
[0. 1.000001 0. ... 1.000001 0. 1.000001]
------ laplacian of n -----------
[ 0. 0. 0. ... 0.34403651 0.34418522 -0.56959499]
------ n ----------------------
[0. 0.99798558 0. ... 0.99798608 0. 0.99798836]
------ laplacian of n -----------
[0. 0. 0. ... -624.23276296 -624.17340055 1037.42146611]
------ n ----------------------
[0. 1. 0. ... 1. 0. 1.]
------ laplacian of n -----------
[ 0. 0. 0. ... 230399.99999999 230400. -806400.00000001]