Solution Blows Up for Linear PDE

1

Hi, everyone. I’m relatively new to FEniCS, but I’ve gone through the textbook. I’m trying to solve an equation:
\partial(h)/\partial(t) = \nabla(\nabla(h)) - Da*h
x = 0: \nabla(h) = 0
x =1: \nabla(h) = Bi*(h - eta0)
t = 0: h = 0
where Da, eta0, and Bi are constants and h = h(x,t) from 0 to 1. The time goes from 0 to 1, and I know this problem has been solved before; I’m trying to reprove it. However, when I compute the solution and view in Paraview, it oscillates and continues to grow larger when I know it should stay between 0 and 1. Here’s my code:

T = 1            
num_steps = 20   
dt = T / num_steps 
Bi = 2		  
Da = 2		  
eta0 = -1	  

mesh = IntervalMesh(100,0,1)
Q = FunctionSpace(mesh,'P',1)

# Define test and regular functions
h = TrialFunction(Q)
ht = TestFunction(Q)

# Define boundary conditions

tol = 1E-14 
class BoundaryX0(SubDomain): 
	def inside(self, x, on_boundary): 
		return on_boundary and near(x[0], 0, tol)

class BoundaryX1(SubDomain): 
	def inside(self, x, on_boundary): 
		return on_boundary and near(x[0], 1, tol)

boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1, 0)
bx0 = BoundaryX0()
bx0.mark(boundary_markers, 0) 
bx1 = BoundaryX1()
bx1.mark(boundary_markers, 1) 

ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)

r = Bi
s = eta0
g = 0

boundary_conditions = {0: {'Robin': (r, s)},
			 1: {'Neumann':	0}}

integrals_N = []
for i in boundary_conditions:
	if 'Neumann' in boundary_conditions[i]: 
		if boundary_conditions[i]['Neumann'] != 0: 
			g = boundary_conditions[i]['Neumann'] 
			integrals_N.append(g*ht*ds(i)) 

integrals_R = []
for i in boundary_conditions:
	if 'Robin' in boundary_conditions[i]: 
		r, s = boundary_conditions[i]['Robin']
		integrals_R.append(r*(h-s)*ht*ds(i))


# Define initial values
h_0 = Constant(0) 
h_n = project(h_0, Q) 

# Define expressions used in variational forms
k = Constant(dt)


F1 =    (h-h_n)*(ht/k)*dx - dot(grad(h), grad(ht))*dx \
	+ Da*h*ht*dx \
	+sum(integrals_R) + sum(integrals_N)

a, L = lhs(F1), rhs(F1)

# Create VTK files for visualization output
vtkfile_h = File('early_bazant/h.pvd')

# Create progress bar
# progress = Progress('Time-stepping')
progress = Progress('Time-stepping', num_steps)
# set_log_level(PROGRESS)

# Time-stepping for potential
h = Function(Q)
t = 0
for n in range(num_steps):

    # Update current time
    t += dt

    # Solve variational problem for time step
    solve(a == L, h)

    # Save solution to file (VTK)
    vtkfile_h << (h, t)

    # Update previous solution
    h_n.assign(h)

I’m running on Windows using Ubuntu 20.04 LTS (GNU/Linux 4.4.0-19041-Microsoft x86_64) and Python3.

2

Also, I realize defining Robin boundaries gets confusing so I’ve tried it using Bi = -2 as well, and the solution still blows up.

3

Please format your latex using $ for encapsulation, and ``` to encapsulate code, such that it is easy to read for other, than can then in turn provide suggestions.
Note that the diffusive term of the heat equation should change sign when doing integration by parts.

This should therefore be positive and not negative

1 Like
4

Sorry for the rough format; should be fixed.
Your note on the correct sign of that term fixed it; thank you so much!