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
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 - cx, 2) + pow(x - 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, 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