Hello,
I’m sorry if this is not the right place to post this since I dont know if it’s really an inherent problem of fenics. Seeking help elsewhere is quite tricky since not everyone is used to fenics.
I am trying to solve successively two time dependent ADEs (which differ only by their boundary conditions) with the discontinuous galerkin method.
I don’t encounter any problem when I solve only one of the two problems.
However, when I solve both, I get the error, in attachment,
I have checked several times for functions and variables with the same name because it is essentially what I found on other threads that i could explain this error and found nothing.
Do you have any idea to solve this issue?
Thank you in advance for your help
Here is the code that I did not truncate in order to have the maximum information on a possible error
from dolfin import *
from fenics import *
#mesh
my_mesh = UnitSquareMesh (50, 50, 'right/left' )
# functional space
my_degree=1
V_dg = FunctionSpace ( my_mesh, 'DG', my_degree )
V_u = VectorFunctionSpace ( my_mesh, 'CG', my_degree + 1 )
# Create velocity function.
u = interpolate ( Constant ( ( "1.0", "0.0" ) ), V_u )
v = interpolate ( Constant ( ( "0.0", "1.0" ) ), V_u )
# Mesh-related functions.
n = FacetNormal ( my_mesh )
h = CellDiameter ( my_mesh )
# Parametres
kappa = Constant (0.1) # diffusion
f = Constant ( 0.0 ) # Source term.
alpha = Constant ( 5.0 ) # Penalty term.
# BCs
class DirichletBoundaryX( SubDomain ):
def inside ( self, x, on_boundary ):
return abs ( x[0] - 1.0 ) < DOLFIN_EPS and on_boundary
bc = DirichletBC ( V_dg, 1.0, DirichletBoundaryX(), "geometric" )
class DirichletBoundaryY( SubDomain ):
def inside ( self, x, on_boundary ):
return abs ( x[1] - 1.0 ) < DOLFIN_EPS and on_boundary
bcy = DirichletBC ( V_dg, 1.0, DirichletBoundaryY(), "geometric" )
# function from chapter 31 fenicsbook
# Advective and convective terms
def advection(phi, psi, uT, n, theta=0.5):
# Define |u * n|
unT = abs(dot(uT("+"), n("+")))
# Contributions from cells
a_cell = - dot(uT*phi, grad(psi))*dx
# Contributions from interior facets
a_int = (dot(uT("+"), jump(psi, n))*avg(phi) + 0.5*unT*dot(jump(phi, n), jump(psi, n)))*dS
return theta*(a_cell + a_int)
def diffusion(phi, psi, k_c, alpha, n, h, theta=0.5):
# Contribution from the cells
a_cell = k_c*dot(grad(phi), grad(psi))*dx
# Contribution from the interior facets
tmp = (alpha("+")/h("+")*dot(jump(psi, n), jump(phi, n))\
- dot(avg(grad(psi)), jump(phi, n))\
- dot(jump(psi, n), avg(grad(phi))))*dS
a_int = k_c("+")*tmp
return theta*(a_cell + a_int)
# Time dependent problem
T = 1.0 # final time
num_steps = 5 # number of time steps
dt = T/num_steps # time step size
t = 0
k = Constant(dt)
# function (trial and test for the two problem)
phi = TrialFunction ( V_dg )
psi = TestFunction ( V_dg )
##############################
# 1st ADE with "bc" as BCs
###############################
# null initial condition
phix=interpolate(Constant(0.0),V_dg)
# Bilinear form.
a1 = ((phi - phix)/k)*psi*dx\
+ advection(phi, psi, u, n, theta=0.5)\
+ diffusion(phi, psi, kappa, alpha, n, h, theta=0.5)\
+ advection(phix, psi, u, n, theta=0.5)\
+ diffusion(phix, psi, kappa, alpha, n, h, theta=0.5)
# Linear form.
L1 = psi*f*dx
form = a1 - L1
a = lhs(form)
L = rhs(form)
phi_h = Function ( V_dg )
t = 0
for n in range(num_steps):
t += dt
solve(a == L, phi_h, bc)
phix.assign(phi_h)
##############################
# 2nd ADE with "bcy" as BCs
###############################
# null initial condition
phiy=interpolate(Constant(0.0),V_dg)
# Bilinear form.
a1y = ((phi - phiy)/k)*psi*dx\
+ advection(phi, psi, v, n, theta=0.5)\
+ diffusion(phi, psi, kappa, alpha, n, h, theta=0.5)\
+ advection(phiy, psi, v, n, theta=0.5)\
+ diffusion(phiy, psi, kappa, alpha, n, h, theta=0.5)
# Linear form.
L1y = psi*f*dx
formy = a1y - L1y
ay = lhs(formy)
Ly = rhs(formy)
phi_hy = Function ( V_dg )
t = 0
for n in range(num_steps):
t += dt
solve(ay == Ly, phi_hy, bcy)
phiy.assign(phi_hy)