Greetings to everyone.
As a part of my study I need to solve convection-diffusion-reaction equation:
\begin{cases}
\frac{\partial \varphi}{\partial t} - div{(\mu grad{\varphi})} + div{({U}\varphi)} + \beta\varphi = f, \; (x, y, t) \in \Omega \times (0, T), \\
U^{(-)}_n \varphi + \mu \ \frac{\partial \varphi}{\partial n} = G, \; (x, y, t) \in \Gamma \times (0, T), \\
\varphi = \varphi_0, \; (x, y) \in \Omega, \; t = 0, \\
\end{cases}
where \mu, b = const > 0, U= (u_x, u_y), U^{(\pm)}_n = [|({U}, {n})| \pm ({U}, {n})]/2.
I wrote the following function to solve it:
import fenics as fe
from fenics import dx, ds
import ufl
import numpy as np
import pylab as plt
import mshr
def StraightProblemSolver(mesh, V, t, T, num_steps, mu, beta, U, f, G, phi_0):
dt = (T - t) / num_steps
phi_n = fe.Function(V)
phi_n.assign(phi_0)
n = fe.FacetNormal(mesh)
U_n_p = (ufl.algebra.Abs(fe.dot(U, n)) + fe.dot(U, n))/2.0
phi, phic = fe.TrialFunction(V), fe.TestFunction(V)
F = phi*phic*dx + dt*mu*fe.dot(fe.grad(phi), fe.grad(phic))*dx - dt*fe.dot(fe.grad(phic), U)*phi*dx + \
dt*beta*phi*phic*dx + dt*U_n_p*phi*phic*ds - (dt*f*phic*dx + phi_n*phic*dx + dt*G*phic*ds)
a = fe.lhs(F)
L = fe.rhs(F)
tcurrent = t
phi = fe.Function(V)
for time_step in range(num_steps):
tcurrent += dt
f.t = tcurrent
G.t = tcurrent
U.t = tcurrent
fe.solve(a == L, phi)
phi_n.assign(phi)
return phi
To test it I also generated some test cases, for example:
def test1_StrainghtProblemSolver():
all_N = [10, 14, 20, 44]
all_num_steps = [50, 100, 200, 1000]
for (N, num_steps) in zip(all_N, all_num_steps):
mesh = fe.UnitSquareMesh(N, N)
V = fe.FunctionSpace(mesh, "P", 1)
# Problem definition here:
mu = 0.000005
beta = 0.01
t = 0.0
T = 1.0
phi_e = fe.Expression("t*(7*x[0]*x[0]-5*x[1]*x[1])", degree=2, domain=mesh, t=np.nan)
f = fe.Expression("(7*x[0]*x[0]-5*x[1]*x[1]) - 0.1*exp(-t)*t*(14.0*x[0]*x[0]+10.0*x[1]*x[1])"+\
"- 4.0*mu*t + beta*t*(7*x[0]*x[0]-5*x[1]*x[1])", \
degree=2, domain = mesh, mu=mu, beta=beta, t=np.nan)
U = fe.Expression(("-0.1*x[0]*exp(-t)", "0.1*x[1]*exp(-t)"), degree=2, domain=mesh, t=np.nan)
class BoundaryExp(fe.UserExpression):
def __init__(self, t, mu, **kwargs):
self.t = t
self.mu = mu
super().__init__(**kwargs)
def eval(self, value, x):
if fe.near(x[0], 0.0, 1e-6):
value[0] = 0.0
elif fe.near(x[0], 1.0, 1e-6):
value[0] = 14.0*self.mu*self.t + 0.1*fe.exp(-self.t)*(7.0 - 5.0*x[1]*x[1])*self.t
elif fe.near(x[1], 0.0, 1e-6):
value[0] = 0.0
elif fe.near(x[1], 1.0, 1e-6):
value[0] = -10.0*mu*self.t
else:
value[0] = 0.0
def value_shape(self):
return ()
G = BoundaryExp(element=V.ufl_element(), degree=2, domain=mesh, mu=mu, t=np.nan)
# End problem definition
phi_e.t = t
phi_0 = fe.interpolate(phi_e, V)
phi = StraightProblemSolver(mesh, V, t, T, num_steps, mu, beta, U, f, G, phi_0)
phi_e.t = T
L2error = fe.errornorm(phi_e, phi, "L2")
Cerror = np.max(np.abs(phi.vector()[:] - fe.interpolate(phi_e, V).vector()[:]))
Cerrorat = np.argmax(np.abs(phi.vector()[:] - fe.interpolate(phi_e, V).vector()[:]))
Cerrorat = V.tabulate_dof_coordinates()[Cerrorat]
print("N:", N, " h^2:", mesh.hmax() ** 2.0, " dt:", (T-t)/num_steps)
print("L2:", L2error, " C: ", Cerror, " at", Cerrorat)
graph = fe.plot(ufl.algebra.Abs(phi - phi_e))
plt.colorbar(graph)
plt.show()
def test4_StrainghtProblemSolver():
all_N = [18, 24, 35, 70]
all_num_steps = [50, 100, 200, 800]
for (N, num_steps) in zip(all_N, all_num_steps):
domain = mshr.Ellipse(fe.Point(1,0.5),1,.5,N)
mesh = mshr.generate_mesh(domain, int(N * .75))
V = fe.FunctionSpace(mesh, "P", 1)
# Problem definition here:
mu = 0.00001
beta = 0.01
t = 0.0
T = 1.0
phi_e = fe.Expression("t*(7*cos(pi*x[0]) + 5*cos(pi*x[1]))", degree=2, domain=mesh, t=np.nan)
f = fe.Expression("(7*cos(pi*x[0]) + 5*cos(pi*x[1]))" + \
" - pi*t*(x[1]*(1-x[1])*(1-2*x[0])*5*sin(pi*x[1]) - x[0]*(1-x[0])*(1-2*x[1])*7*sin(pi*x[0]))" + \
" - mu*pi*pi*t*(7*cos(pi*x[0]) + 5*cos(pi*x[1]))" + \
" + beta*t*(7*cos(pi*x[0]) + 5*cos(pi*x[1]))", degree=2, domain=mesh, mu=mu, beta=beta, t=np.nan)
U = fe.Expression(("-x[0]*(1-x[0])*(1-2*x[1])", "x[1]*(1-x[1])*(1-2*x[0])"), degree=2, domain=mesh, t=np.nan)
class BoundaryExp(fe.UserExpression):
def __init__(self, t, mu, **kwargs):
self.t = t
self.mu = mu
super().__init__(**kwargs)
def eval(self, value, x):
nx = 2.0*(x[0] - 1.0)
ny = 8.0*(x[1] - 0.5)
l = fe.sqrt(nx*nx + ny*ny)
nx = nx / l
ny = ny / l
U_n = -nx*x[0]*(1-x[0])*(1-2*x[1]) + ny*x[1]*(1-x[1])*(1-2*x[0])
U_n_m = (abs(U_n) - U_n)/2.0
Der = -np.pi*self.t*(7.0*np.sin(np.pi*x[0])*nx + 5.0*np.sin(np.pi*x[1])*ny)
value[0] = U_n_m*self.t*(7.0*np.cos(np.pi*x[0]) + 5.0*np.cos(np.pi*x[1])) + self.mu*Der
def value_shape(self):
return ()
G = BoundaryExp(element=V.ufl_element(), degree=2, domain=mesh, mu=mu, t=np.nan)
# End problem definition
phi_e.t = t
phi_0 = fe.interpolate(phi_e, V)
phi = StraightProblemSolver(mesh, V, t, T, num_steps, mu, beta, U, f, G, phi_0)
phi_e.t = T
L2error = fe.errornorm(phi_e, phi, "L2")
Cerror = np.max(np.abs(phi.vector()[:] - fe.interpolate(phi_e, V).vector()[:]))
Cerrorat = np.argmax(np.abs(phi.vector()[:] - fe.interpolate(phi_e, V).vector()[:]))
Cerrorat = V.tabulate_dof_coordinates()[Cerrorat]
print("N:", N, " h^2:", mesh.hmax() ** 2.0, " dt:", (T-t)/num_steps)
print("L2:", L2error, " C: ", Cerror, " at", Cerrorat)
graph = fe.plot(ufl.algebra.Abs(phi - phi_e))
plt.colorbar(graph)
plt.show()
The problem is that I’m getting large error at the points on the boundary. I’m new to fenics and also don’t have enough experience with finite element method, so I can’t figure it out on my own if I made a mistake or if this is a expected result.
I would appreciate any help. Thanks in advance.
