Same time dependent temperature problem in two different ways: one works, one doesn't. Why?

variational formulation
1

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

T_solo_t00001
T_solo_t00001801×355 14.9 KB

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)

T_coup_t00002
T_coup_t00002816×355 39.6 KB

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.

2

I just figured it out! I didn’t define the InitialConditions class, therefore my u0, p0, T0, were not of the right form. I am still new at FEniCS so unfortunately I don’t know why it works one way and not the other… But I do know that adding

class InitialConditions(UserExpression):
    def __init__(self, **kwargs):
        super().__init__(**kwargs)

    def eval(self, values, x):
        values[0] = 0.0
        values[1] = 0.0
        values[2] = 0.0
        values[3] = 0.0
    def value_shape(self):
        return (4,)

to the above, and then

w0 = Function(W)
w = Function(W)
u0, p0, T0 = split(w0)
u, p, T    = split(w)
#Interpolate the initial conditions
w_init = InitialConditions(degree=1)
w.interpolate(w_init)
w0.interpolate(w_init)

for w and w0 (instead of what I wrote above), everything works! I don’t know why… I will mark this subject as closed to not waste anyone’s time, but if you read this and know why that fixed it, I would be really grateful if you could reply.

1 Like