Unable to define linear variational problem in a splitting method

Errors
1

Hello,

I’m trying to solve a system of two equations on a hollow cylinder. The first equation is the time-varying 3D heat diffusion equation with a constant source term for cylinder temperature (Tc). The other is a steady 1D equation for temperature of fluid (Tw) flowing through the “pipe” along the axis of the hollow cylinder; this equation only lives on the surface of the hole through the cylinder.

My source term is highly nonlinear in cylinder temperature Tc and so in the past I’ve used FEniCS with a splitting method for problems that don’t have “pipes” in them. For the situation with pipes, I’ve managed to solve a simpler version of this problem (one with a constant source term and no splitting used during time stepping) using mixed function spaces. In trying to use the actual nonlinear source term, it seemed easier to forego mixed function spaces and solve with separate function spaces.

I’ve followed the Navier-Stokes example in the tutorial, but I’m getting an error "
Unable to define linear variational problem a(u, v) = L(v) for all v" similar to this but I can’t see that the source of my error is the same.

The basic model is

\eqalign{ {{\partial T_c} \over {\partial t}} &= \nabla^2T_c + \dot e \\ {{d T_w} \over {d z}} &= \big( T_c |_{r=r_p} - T_w \big) }

with Robin conditions all around the cylinder and a constant inlet fluid temperature. My minimal example is (and I hope this isn’t too long):

from fenics import *
from mshr import *


# ...........................
# Parameters

# domain geometry
rc = 1              # radius of the cylinder
rp = 0.2            # radius of the pipe
L = 1                # length of the cylinder
# note that here the center of the domain is (0,0)
cx = 0
cy = 0

# mshr parameters
segments = 64
resolution = 30        

fe_degree_c = 2      # degree of the finite element for cylinder
fe_degree_w = 1      # degree of the finite element for water

# boundary and initial condtions re: cylinder
Tc_i  = 0           # initial temperature of cylinder
T_amb = 0.0          # ambient temperature
h_amb  = 0.001       # ambient convection coefficient at outer radius of cylinder

# boundary and initial conditions re: pipe
Tw_e = 0             # entering water temperature
h_pipe  = 100        # convection coefficient in the pipe

# source term
dot_e = 1

# for time-stepping
t_f = 10             # final time
dt = 0.25        # time step size

XDMF_file_Tc = XDMFFile('Tc.xdmf')
XDMF_file_Tw = XDMFFile('Tw.xdmf')



# ...........................
# Problem setup

# mesh
domain = Cylinder(Point(cx, cy, L), Point(cx, cy, 0), rc, rc, segments) - Cylinder(Point(cx, cy, L), Point(cx, cy, 0), rp, rp)
mesh = generate_mesh(domain, resolution)

# boundaries
boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
boundary_markers.set_all(0)

class OnBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary

class OnPipe(SubDomain):
    def inside(self, x, on_boundary):
        # use a small increase over the radius to provide a little
        # "robustness" in finding the inner pipe surface
        return pow(x[0] - cx, 2) + pow(x[1] - cy, 2) <= pow(rp, 2)*(1+0.01*rp)

# set outer boundary marker to "2"
OnBoundary = OnBoundary()
OnBoundary.mark(boundary_markers, 2)

# set inner pipe marker to "3"
OnPipe = OnPipe()
OnPipe.mark(boundary_markers, 3)

# for the boundary
ds = Measure('ds', subdomain_data=boundary_markers)
# for outer boundary
ds = ds(2)
# for the pipe surface
dp = ds(3)



# function, trial, and test spaces for cylinder temperature
V = FunctionSpace(mesh, 'P', int(fe_degree_c))
Tc = TrialFunction(V)
v = TestFunction(V)

# trial functions for the intermediate steps in the splitting scheme for cylinder temperature
Tc_a = TrialFunction(V)
Tc_b = TrialFunction(V)

# function, trial, and test spaces for water temperature
Q = FunctionSpace(mesh, 'P', int(fe_degree_w))
Tw = TrialFunction(Q)
q = TestFunction(Q)



# initial condition
# Tc_n represents the solution at the previous (and initial) timestep
Tn = interpolate(Constant(Tc_i), V)   
Tw = interpolate(Constant(Tw_e), Q)


# Robin boundary condition for the outer surfaces of the cylinder
# in the 'real' code there are multiple Robin and Neumann boundary conditions
# so I keep to that basic structure here; this is applied in the middle
# splitting step so Tc_b is the one to use here
cylinder_bc = dt*(h_amb)*(Tc_b-T_amb)*v*ds(2)


# Dirichlet boundary condition for water inlet temperature
class AtInlet(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[2], height_m)

inlet = AtInlet()
pipe_bc = DirichletBC(Q, Constant(Tw_e), inlet)



# functions and variational formulations needed in the splitting method
# use Strang RDR (reaction-diffusion-reaction) splitting
#   at each timestep:
#       - reaction: advance temperature "adiabatically": no heat transfer, but use egen for 1/2 timestep
#       - diffusion: using this new temperature field, simulate a purely diffusive problem (i.e. only heat transfer) for a full timestep
#       - reaction: with the newly diffused temperature field, advance temperature adiabatically again for a second 1/2 timestep

# first step is a simple Euler advance of (adiabatic) temperature over half a time step
def TadvAdb(delt, T):
    "advance (adiabatic) temperature using Cervera egen by Euler method"
    V = T.function_space()
    Tmesh = V.mesh()
    degree = V.ufl_element().degree()
    W = FunctionSpace(Tmesh, 'P', degree)
    Tnext = project(T + (delt/2)*dot_e, W, solver_type='iterative')
    return Tnext


# second step is a pure diffusion step; the Robin boundary condition is inherently included here
#  Tn here will be the temperature from the first step: Tn = TadvAdb()
F2 = Tc_b*v*dx + dt*dot(grad(Tc_b), grad(v))*dx - Tn*v*dx + cylinder_bc + dt*h_pipe*(Tn-Tw)*v*dp
a2 = lhs(F2)
L2 = rhs(F2)

# third step is another call to function TadvAdb(...)


# for the water temperature: a fourth step
F4 = (Tw.dx(2))*q*dp - (Tc - Tw)*q*dp
a4 = lhs(F4)
L4 = rhs(F4)


A2 = assemble(a2)
A4 = assemble(a4)




# ...........................
# Time stepping loop
t = 0

Tc = Function(V)
Tc_a = Function(V)
Tc_b = Function(V)
Tw = Function(Q)

while t < t_f/dt:
    
    # udate current time and iteration number
    t += dt

    # first splitting step: adiabatic advance over 1/2 timestep
    Tc_a = TadvAdb(dt, Tn)
    Tn.assign(Tc_a)

    # second splitting step: pure diffusion over full timestep
    b2 = assemble(L2)
    solve(A2, Tc_b.vector(), b2, 'bicgstab', 'hypre_amg')
    Tn.assign(Tc_b)

    # third splitting step: adiabatic advance over 1/2 timestep
    Tc =  TadvAdb(dt, Tn)
    Tn.assign(Tc)

    # solve for water
    solve(a4 == L4, Tw, pipe_bc)


    XDMF_file_Tc.write(Tc, t)
    XDMF_file_Tw.write(Tw, t)

print(' ')
print('Done')

The error occurs at

solve(a4 == L4, Tw, pipe_bc)

*** Error: Unable to define linear variational problem a(u, v) = L(v) for all v.
*** Reason: Expecting the solution variable u to be a member of the trial space.
*** Where: This error was encountered inside LinearVariationalProblem.cpp.
*** Process: 0

*** DOLFIN version: 2019.1.0

2

It would be best if you could reduce this to a MWE since few of us have time to debug entire research projects.

At first glance, do you think this is a valid weak formulation?

F4 = (Tw.dx(2))*q*dp - (Tc - Tw)*q*dp
a4 = lhs(F4)
L4 = rhs(F4)

And what do you expect will happen to the first term in F4 when you’ve prescribed fe_degree_w = 1?

Furthermore, you’ve assigned:

Tw = TrialFunction(Q)

followed almost immediately by

Tw = interpolate(Constant(Tw_e), Q)
3

My apologies for the long example. I’ll pare it down next time.

And what do you expect will happen to the first term in F4 when you’ve prescribed fe_degree_w = 1 ?

Tw = interpolate(Constant(Tw_e), Q)

There were my problems, I’ve got it to work. Thanks.

1 Like