Elliptic-Parabolic System of Equations

Question

I am trying to solve a coupled system of equations made up of a parabolic and elliptic equation. For some reason, when I use a mixed space and solve the equations in a single solve call, the parabolic equation remains constant. Is there a way to solve these equations as a single system? Below is a minimal example.

Import necessary packages

from fenics import *
from dolfin import *
import numpy as np

Define parameters and constants

tol = 1E-14
density = 1.0
specificHeat = 1.0
ambAlpha = 1.0
ambTemp = 1.0
coolAlpha = 1.0
coolTemp = 0.5
thermalConductivity = 1.0
perfusionRate = 1.0
bloodTemp = 1.0
radDiffusion = 1.0
absCoefficient = 5.0
appliedPower = 18.0
appliedArea = 0.2

Set time-stepping parameters

T = 10.0
num_steps = 100
dt = T / num_steps

Create mesh and define function space

nx = ny = 10
mesh = UnitSquareMesh(nx,ny)

boundaryMarkers = FacetFunction('size_t', mesh)

Set the ambient boundary

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

bx1 = Boundary1()
bx1.mark(boundaryMarkers, 1)

Set the radiative boundary

class Boundary2(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0, tol) and (x[1] >= 0.25 and x[1] <= 0.75)

bx2 = Boundary2()
bx2.mark(boundaryMarkers, 2)

Set the cooled boundary

class Boundary3(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0, tol) and (x[1] < 0.25 or x[1] > 0.75)

bx3 = Boundary3()
bx3.mark(boundaryMarkers, 3)

P1 = FiniteElement("Lagrange", triangle, 1)
element = MixedElement([P1,P1])
V = FunctionSpace(mesh, element)

Define boundary conditions

bcs = [] # No Dirichlet boundary conditions

Redefine the measure ds for boundary integrals

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

Define variational problem

u = Function(V)
u1, u2 = split(u)

u_n = Function(V)
u1_n, u2_n = split(u_n)

v1, v2 = TestFunctions(V)

Define initial conditions

u_n = interpolate(Constant((0.5,0.0)), V)

Create progress bar

progress = Progress('Time-stepping')
set_log_level(PROGRESS)

Save solution to file in VTK format

vtkfile1 = File('paraview/temperature.pvd')
vtkfile2 = File('paraview/radiation.pvd')

Time-stepping

t = 0

k = Constant(dt)

for n in range(num_steps):
# Update time
t += dt

# Solve equations
F = density * specificHeat * ((u1 - u1_n) / k) * v1 * dx \
    - ambAlpha * (ambTemp - u1) * v1 * ds(1) \
    - coolAlpha * (coolTemp - u1) * v1 * ds(2) \
    - coolAlpha * (coolTemp - u1) * v1 * ds(3) \
    + dot(thermalConductivity * grad(u1), grad(v1)) * dx \
    - absCoefficient * u2 * v1 * dx \
    - perfusionRate * (bloodTemp - u1) * v1 * dx

F += dot(grad(u2), grad(v2)) * dx             \
    - appliedPower/appliedArea * v2 * ds(3)  \
    + 1/2 * u2 * v2 * ds(1)                   \
    + absCoefficient * u2 * v2 * dx

solve(F == 0, u)

# Write solution to VTK format
_u1, _u2 = u.split()
vtkfile1 << (_u1,t)
vtkfile2 << (_u2,t)

# Update previous solution
u_n.assign(u)

# Update progress boundary
progress.update(t / T)

Comments

Hi, did you include the entire code? Because I guess that u_n is the previous step, and it is not being updated.

---

It appears that a portion of the code got cut off when copying. I've added the rest, the update for u_n was cut off, but it is in the version that is giving me trouble.

PS: Is there a trick to getting the code to appear properly? It was complete in the original post but wasn't displayed.

---

No idea, and also no idea why it is not working. Maybe you could try splitting u_n in each iteration because the splitted functions could be copies and not references, but no idea honestly.