Hi,
I was solving a simple parabolic heat equation with step initial condition in the middle of the domain and no flux BC at the boundaries. For simplicity, take 1D interval as the domain. When I solve it, for smaller number of elements, say 5, no flux BC does not hold and also I get overshoot at the nodes next to discontinuous edge (at least at early time steps). I also did manual FE discretization and I saw the same behavior. However, when I increase the number of the elements (like > 11), all of the mentioned issues are resolved. I was wondering why is it like that? Also, is there any workaround for it? I appreciate if anyone can help me.
Here is my code:
Governing equation: time-dependent heat equation
ut = D grad.grad(u) + f
grad_n(u).n = 0 at x = 0,1
u(x,0) = step at x=0.5
from dolfin import *
from numpy import *
import matplotlib.pyplot as plt
set_log_active(False)
mesh = UnitIntervalMesh(5)
x_mesh = mesh.coordinates()[:,0]
V = FunctionSpace(mesh,"CG",1)
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0.0)
g = Constant(0.0)
D = Constant(0.01)
dt = 0.1
t = 0.
T = 1
u_ic = Expression('(x[0]<0.5)')
u_1 = interpolate(u_ic,V)
theta = 1.0 # theta in [0.5,1]
a = vudx + thetadtDinner(nabla_grad(v), nabla_grad(u))dx
L = vu_1dx - (1.-theta)dtDinner(nabla_grad(v), nabla_grad(u_1))dx + dtDvgds + dtvf*dx
A = assemble(a)
b = None
solver = LUSolver(A)
solver.parameters["reuse_factorization"]=True
w = Function(V)
count = 0
fn = 0
while (t<T):
if mod(count,1)==0:
fn += 1
plt.figure(fn)
plt.plot(x_mesh,u_1.vector().array(),'-o',markersize=8,lw=2)
plt.xlabel('Domain')
plt.ylabel('Solution')
plt.ylim(-0.1,1.1)
plt.grid()
b = assemble(L, tensor=b)
solver.solve(w.vector(),b)
u_1.assign(w)
t += dt
count += 1
plt.show()
print A.array()
