 # Simple time dependent problem, how to *stop* updating temperature dependent expressions?

Hi everyone.
It’s a bit of a weird question. I have a simple time dependent problem that works well (which I also described here). I am trying to write code that is flexible, so that other people can use it for more complicated user expressions.

For illustration purposes, consider the source term

``````class QExactDef(UserExpression):
def __init__(self, t, **kwargs):
super().__init__(kwargs)
self._t = t

def eval_cell(self, value, x, ufc_cell):
value = -2*self._t(x) + x*x + 2*self._t(x)*self._t(x)*x*(np.sin(2*np.pi*x))

def value_shape(self):
return ()
``````

which I initialise with

``````t = Constant(0.0)
Qexact = QExactDef(t)
``````

Then I build my variational forms and boundary conditions. The timestepping occurs in a while loop for a parameter t2 (I was worried I would get confused between t and t2 so I kept them as separate), and inside it I assign the new t2 for Qexact, namely

``````while(t2 < T_final):
t2 += dt
T0.vector()[:] = T.vector()
Qexact._t.assign(t2)
solver.solve(problem, T.vector())
file_T << (T, t2)
``````

No problems there. However, what if I want to define Q the same way as above, but I don’t want Q to advance in time? Intuitively, I assumed that it would just be the same as above, but commenting out the line

``````# Qexact._t.assign(t2)
``````

which (to me) means that Qexact is evaluated at 0.0 once, and then remains the same value at every timestep. However, somehow Qexact is still progressing in time, as the solution I get is exactly the same whether or not I add that line. I also tried

``````Qexact._t.assign(0.0)
``````

with the same results, Qexact still evolves in time. I suppose my question is, where does the function actually get evaluated at t2?

I wanted to write flexible code where all of my parameters could in principle depend on time, but then the user can choose which ones to fix as a constant and which ones to evolve. I am concerned that I don’t understand the timestepping enough to allow for that.

If you are interested, the full code for the MWE 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 InitialConditions(UserExpression):
def __init__(self, **kwargs):
super().__init__(**kwargs)

def eval(self, value, x):
value = 0.0

def value_shape(self):
return()

# 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 = -2*self._t(x) + x*x + 2*self._t(x)*self._t(x)*x*(np.sin(2*np.pi*x))

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 = -self._t(x) * np.sin(2 * np.pi * (1 - x))
value = 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 = self._t(x)*(x*x)

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)
T0 = Function(W)

T_init = InitialConditions(degree=1)
T.interpolate(T_init)
T0.interpolate(T_init)

# 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_3.pvd", "compressed")
# Step in time
t2 = 0.0 # this is the parameter that goes in the while loop
T_final = 10*dt

bcs = [DirichletBC(W, Texact, full_boundary)]
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)
``````

which is pretty much the Cahn-Hilliard example.

Thank you!

In your MWE, have you commented all of the following lines, or only the line `Qexact._t.assign(t2)`?

In your MWE, `Qexact._t`, `uexact._t`, and `Texact._t` all refer to the same object; that is:

The following function calls

have exactly the same effect as calling

``````t.assign(t2)
``````

Therefore, calling any one of `Qexact._t.assign(t2)`, `uexact._t.assign(t2)`, or `Texact._t.assign(t2)` will change the value that is “seen” by all three of `Qexact`, `uexact`, and `Texact`.

You could change this behavior by doing:

``````t = Constant(0) # t that is evaluated in the functions
tQ = Constant(0) # t that is evaluated in Qexact
Qexact = QExactDef(tQ)
uexact = uExactDef(t)
Texact = TExactDef(t)
``````

Hi!
In the text I only mentioned Qexact but I also had the other two functions that had the same issue. I should have removed them for the sake of the MWE, my bad.
I see I have completely misunderstood where the evaluation actually happens! Now that I know it’s when I call the class the first time, it makes perfect sense.

Thanks again!!