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[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 ()
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.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[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)
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!