Hi everyone,
I am trying to implement a coupled flow and temperature time dependent problem. Since I am new to these problems, I have started with simple steps. I have one code that solves the time dependent problem for the temperature, with a given velocity. This code works perfectly!
However, when I try to run the coupled example, suddenly I get the wrong temperature solution (for the same velocity). I think it’s some issue to do with the “split” function, as the second more complicated example is a mixed space problem. However, no matter how much I have tried, I can’t seem to fix it.
The simple only temperature code with purely Dirichlet conditions is
from dolfin import *
import warnings
import numpy as np
# Choose a mesh.
mesh = RectangleMesh(Point(0, 0), Point(1, 1), 10, 10, "crossed")
# Define mesh normal
n = FacetNormal(mesh)
# Prevent warning
warnings.simplefilter(action='ignore', category=FutureWarning)
# define the full boundary
def full_boundary(x, on_boundary):
return on_boundary
# Class for assembling the matrix A and the vector b, as well as apply bcs
class GoverningEquations(NonlinearProblem):
def __init__(self, a, L, bcs):
NonlinearProblem.__init__(self)
self.L = L
self.a = a
self.bcs = bcs
def F(self, b, x):
assemble(self.L, tensor=b)
for bc in self.bcs:
bc.apply(b, x)
def J(self, A, x):
assemble(self.a, tensor=A)
for bc in self.bcs:
bc.apply(A)
# Define source term for govn equation
class QExactDef(UserExpression):
def __init__(self, t, **kwargs):
super().__init__(kwargs)
self._t = t
# self.D = D
def eval_cell(self, value, x, ufc_cell):
value[0] = -2*self._t(x) + x[0]*x[0] + 2*self._t(x)*self._t(x)*x[0]*(np.sin(2*np.pi*x[1]))
def value_shape(self):
return ()
# Define velocity field
class uExactDef(UserExpression):
def __init__(self, t, **kwargs):
super().__init__(kwargs)
self._t = t
def eval_cell(self, value, x, ufc_cell):
value[0] = -self._t(x) * np.sin(2 * np.pi * (1 - x[1]))
value[1] = Constant(0.0)
def value_shape(self):
return (2,)
# Exact solution for the temperature
class TExactDef(UserExpression):
def __init__(self, t, **kwargs):
super().__init__(kwargs)
self._t = t
def eval_cell(self, value, x, ufc_cell):
value[0] = self._t(x)*(x[0]*x[0])
def value_shape(self):
return()
# Define function space
S = FiniteElement("CG", triangle, 1)
W = FunctionSpace(mesh, S)
# Define Function for the previously converged step and TestFunction(s)
dw = TrialFunction(W) # used for the bilinear form
s1 = TestFunction(W) # used for the linear form
theta = 0.5 # parameter for the theta scheme in time
dt = 0.0001 # timestep
t = Constant(0) # t that is evaluated in the functions
# Form compiler options
parameters["form_compiler"]["optimize"] = True
parameters["form_compiler"]["cpp_optimize"] = True
x = SpatialCoordinate(mesh)
# Define the variational form
facet_f = MeshFunction("size_t", mesh, mesh.topology().dim() - 1) # FACET function
ds = Measure("ds", domain=mesh, subdomain_data=facet_f)
# Evaluate Q, u and T
Qexact = QExactDef(t)
uexact = uExactDef(t)
Texact = TExactDef(t)
T0 = Function(FunctionSpace(mesh, S))
T0.interpolate(Texact)
T = Function(W)
# Define timestepping for T
T_mid = (1.0 - theta)*T0 + theta*T
L0 = (1/dt) * (T*s1 - T0*s1)*dx() + (dot(uexact, grad(T_mid)) * s1 + inner(grad(s1), grad(T_mid))) * dx() \
- s1 * dot(grad(T_mid), n) * ds()
L1 = - Qexact*s1 * dx()
L = L0 + L1
J2 = derivative(L, T, dw)
# Output file
file_T = File("Results/temperature_solo.pvd", "compressed")
# Step in time
t2 = 0.0 # this is the parameter that goes in the while loop
T_final = 10*dt
file_T << (T0, t2)
bcs = [DirichletBC(W, Texact, full_boundary)]
# Evaluating the problem
problem = GoverningEquations(J2, L, bcs)
# solver parameters
solver = NewtonSolver()
solver.parameters["linear_solver"] = "lu"
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6
while(t2 < T_final):
t2 += dt
T0.vector()[:] = T.vector()
# advance Q and u in time
Qexact._t.assign(t2)
uexact._t.assign(t2)
Texact._t.assign(t2)
solver.solve(problem, T.vector())
file_T << (T, t2)
This gives the following solution for T at t= 0.0001
which is the correct answer! However, when I do the same temperature problem, coupled with Stokes flow, it doesn’t work. In the second problem, the coupling is only one way, so the velocity field is unaffected by the temperature, but the temperature needs the velocity solution. The code is
from dolfin import *
from mshr import *
import matplotlib.pyplot as plt
import numpy as np
import warnings
# Choose a mesh. You can also import a previously adapted mesh
# mesh = Mesh('Meshes/CoupledRefinedMeshGamma5Cylindric.xml')
mesh = RectangleMesh(Point(0, 0), Point(1, 1), 10, 10, "crossed")
n = FacetNormal(mesh)
warnings.simplefilter(action='ignore', category=FutureWarning)
tol = 1E-14
def full_boundary(x, on_boundary):
return on_boundary
def top_boundary(x, on_boundary):
return on_boundary and abs(x[1]) < tol
# Define Taylor--Hood function space W
V = VectorElement("CG", triangle, 2)
Q = FiniteElement("CG", triangle, 1)
S = FiniteElement("CG", triangle, 1)
W = FunctionSpace(mesh, MixedElement([V, Q, S]))
# Define Function for the previously converged step and TestFunction(s)
# w0 = Function(W)
dw = TrialFunction(W) # used for the bilinear form
(v, q, s1) = TestFunctions(W) # used for the linear form
theta = 0.5 # parameter for the theta scheme in time
dt = 0.0001
t = Constant(0)
# Classes for the time dependent exact solutions and boundary conditions
class uExactDef(UserExpression):
def __init__(self, t, **kwargs):
super().__init__(kwargs)
self._t = t
def eval_cell(self, value, x, ufc_cell):
value[0] = -self._t(x)*np.sin(2*np.pi*(1-x[1]))
value[1] = Constant(0.0)
def value_shape(self):
return(2,)
class fExactDef(UserExpression):
def __init__(self, t, **kwargs):
super().__init__(kwargs)
self._t = t
self.mu = 1
def eval_cell(self, value, x, ufc_cell):
value[0] = 4*np.pi*np.pi*self._t(x)*self.mu*np.sin(2*np.pi*x[1])
value[1] = 0.0
def value_shape(self):
return(2,)
class pExactDef(UserExpression):
def __init__(self, t, **kwargs):
super().__init__(kwargs)
self._t = t
def eval_cell(self, value, x, ufc_cell):
value[0] = 0.0
def value_shape(self):
return()
class TExactDef(UserExpression):
def __init__(self, t, **kwargs):
super().__init__(kwargs)
self._t = t
def eval_cell(self, value, x, ufc_cell):
value[0] = self._t(x)*(x[0]*x[0])
def value_shape(self):
return()
class QExactDef(UserExpression):
def __init__(self, t, **kwargs):
super().__init__(kwargs)
self._t = t
# self.D = D
def eval_cell(self, value, x, ufc_cell):
value[0] = -2*self._t(x) + x[0]*x[0] + 2*self._t(x)*self._t(x)*x[0]*(np.sin(2*np.pi*x[1]))
def value_shape(self):
return ()
class GoverningEquations(NonlinearProblem):
def __init__(self, a, L, bcs):
NonlinearProblem.__init__(self)
self.L = L
self.a = a
self.bcs = bcs
def F(self, b, x):
assemble(self.L, tensor=b)
for bc in self.bcs:
bc.apply(b, x)
def J(self, A, x):
assemble(self.a, tensor=A)
for bc in self.bcs:
bc.apply(A)
# Form compiler options
parameters["form_compiler"]["optimize"] = True
parameters["form_compiler"]["cpp_optimize"] = True
# Evaluate the boundary condition expressions
uexact = uExactDef(t)
pexact = pExactDef(t)
Texact = TExactDef(t)
Qexact = QExactDef(t)
fexact = fExactDef(t)
x = SpatialCoordinate(mesh)
mu = Constant(1.0)
def epsilon(v):
return sym(as_tensor([[v[0].dx(0), v[0].dx(1)],
[v[1].dx(0), v[1].dx(1)]]))
# stress tensor
def sigma(v, p):
return 2*mu*epsilon(v)-p*Identity(2)
# Define the variational form
facet_f = MeshFunction("size_t", mesh, mesh.topology().dim() - 1) # FACET function
CompiledSubDomain('near(x[1], 0.0)').mark(facet_f, 0)
CompiledSubDomain('near(x[0], 1.0)').mark(facet_f, 1) # wall
CompiledSubDomain('near(x[1], 1.0)').mark(facet_f, 2) # outflow
CompiledSubDomain('near(x[0], 0.0)').mark(facet_f, 3)
# Create the measure
ds = Measure("ds", domain=mesh, subdomain_data=facet_f)
print("Evaluating the steady state solution")
u0 = Function(FunctionSpace(mesh, V))
p0 = Function(FunctionSpace(mesh, Q))
T0 = Function(FunctionSpace(mesh, S))
w0 = Function(W)
# Interpolate the initial conditions
#
u0.interpolate(uexact)
p0.interpolate(pexact)
T0.interpolate(Texact)
print("Beginning the unsteady part of the problem")
w = Function(W)
### I THINK THE PROBLEM IS HERE ###
u, p, T = split(w)
# u^{k+theta}
u_mid = (1.0 - theta)*u0 + theta*u
T_mid = (1.0 - theta)*T0 + theta*T
p_mid = (1.0 - theta)*p0 + theta*p
# staying with the notation from Cahn Hilliard (for future questions)
L0 = (1/dt) * (T*s1 - T0*s1)*dx() + (dot(u_mid, grad(T_mid)) * s1 + inner(grad(s1), grad(T_mid))) * dx()\
- Qexact*s1 * dx() - s1 * dot(grad(T_mid), n) * ds()
L1 = inner(sigma(u, p), epsilon(v)) * dx() - div(u) * q * dx() \
- dot(fexact, v) * dx() - dot(dot(sigma(u, p), v), n) * ds(1) - dot(dot(sigma(u, p), v), n) * ds(3)
L = L0 + L1
J2 = derivative(L, w, dw)
# Output file
file_u = File("Results/velocity_coupled.pvd", "compressed")
file_p = File("Results/pressure_coupled.pvd", "compressed")
file_T = File("Results/temperature_coupled.pvd", "compressed")
# Step in time
t2 = 0.0
T_final =5*dt
file_u << (u0, t2)
file_p << (p0, t2)
file_T << (T0, t2)
bcu = DirichletBC(W.sub(0), uexact, full_boundary)
bcP_inflow = DirichletBC(W.sub(1), pexact, top_boundary)
bcT = DirichletBC(W.sub(2), Texact, full_boundary)
bcs = [bcu, bcP_inflow, bcT]
problem = GoverningEquations(J2, L, bcs)
solver = NewtonSolver()
solver.parameters["linear_solver"] = "lu"
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6
while(t2 < T_final):
t2 += dt
print("t2 is ", t2)
### OR MAYBE THE PROBLEM IS HERE ###
w0.vector()[:] = w.vector()
uexact._t.assign(t2)
pexact._t.assign(t2)
Texact._t.assign(t2)
Qexact._t.assign(t2)
fexact._t.assign(t2)
solver.solve(problem, w.vector())
file_u << (w.split()[0], t2)
file_p << (w.split()[1], t2)
file_T << (w.split()[2], t2)
I apologise for the length. The class definitions are the same in both examples, apart from a top boundary (for specifying pressure) and the velocity source term (fexact). The only equation with a time derivative is the equation involving temperature. The velocity solution I get is correct, however the temperature looks correct only for t = 0.0001 (the first timestep), but then it looks like (for t = 0.0002)
It looks like the time-stepping only occurs for the Dirichlet boundary conditions. I think that either T_mid or T0 is not being updated at each time-step, but I don’t know why (since I am doing the same thing as in the first example that did work) .
I would really appreciate any advice.