Error with Dirchlet condition

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.

Your formulation isn’t sufficient for elements with a discontinuous basis. See the standard paper and the dolfin DG demo.

1 Like

I very much appreciate that.

1 Like