Hello everyone,
I have a problem with a problem of reaction-diffusion-advection. i have implemented this as i shown below.
from fenics import *
import numpy as np
# Create mesh and define function space
num_steps = 5000 # number of time steps
dt = 1e-3 # time step size
D = 0.5 # diffusion coefficient
tol = 1e-6
tolb = 1e-4
impresion = 10
# Create mesh and define function space
nx = ny = 100
mesh = RectangleMesh.create([Point(0, -0.4),Point(1.2, 0.4)],[nx, ny],CellType.Type.quadrilateral)
V = FunctionSpace(mesh,'DG', 2)
#Define parameters
G_0 = 0.25
x0 = 0.2
c = 0.3
lamda = 0.04
H = 0.4
k = 0.04
xt = 1.2
x0 = 0
U_0 = 3.2
teta = 0.8
#Define the Dirichlet condition
phi_D = Expression('c*G_0/2.0*exp(-pow((2.0*x[1])/(pi*H),2.0))*sin((19.0*c/pi)*t)', degree=2,c=c,G_0=G_0,H=H,t=0)
def boundary_D(x):
return near(x[0], 0.0, tolb)
bcD = DirichletBC(V, phi_D, boundary_D)
#Define velocity
um = 'U_0*(x[1]*exp(-pow((5*x[1])/(H),2))*cos((3*pi*x[0])/(xt))+tanh(5/2*(2*x[1]+H)/(teta*H))-tanh(5/2*(2*x[1]-H)/(teta*H))-1)'
vm = 'U_0*(-(3*pi*pow(H,2))/(50*xt)*exp(-pow((5*x[1])/(H),2)))*sin((3*pi*x[0])/(xt))'
vel= Expression((um,vm),degree=2,U_0=U_0,H=H,xt=xt,teta=teta)
#Define variational problem
phi = TrialFunction(V)
w = TestFunction(V)
phi_n = interpolate(Constant(0), V)
F = ((phi - phi_n) / dt)*w*dx + dot(vel, grad(phi))*w*dx\
-D*dot(grad(phi), grad(w))*dx
a, L = lhs(F),rhs(F)
# Compute solution
phi = Function(V)
# Create VTK files for visualization output
vtkfile_phi = File('punto2cfd/phi.pvd')
# Time-stepping and increment
t = 0
z = 0
inc = 0
for n in range(num_steps):
# Update current time
t += dt
inc += 1
G.t = t
phi_D.t=t
# Solve variational problem for time step
solve(a==L, phi, bcD)
num=inc/impresion-1
if (int(num)==z) or (int(num-z)==0):
# Save solution to file (VTK)
vtkfile_phi << (phi, t)
z +=1
# diff_L2 = errornorm(phi, phi_n, 'L2')/dt
# if diff_L2<tol:
# break
# Update previous solution
phi_n.assign(phi)
# Update progress bar
# progress.update(t / T)
print('t_f=',t)
# Hold plot
The issue is that the solution is zero at every time, it looks like the dirichlet condition have not been put. I tried to implement this with Lagrangian elements and the solution is unstable. I can not find the error in the code. I would be very grateful if you help me with this. Thank you in advance.