Hello,

I want to use the AdaptiveNonlinearVariationalSolver to adapt my mesh automatically, but my problem is nonlinear and in a mixed space. I can solve my problem generally with solve(), but now that I want to put it in the correct format for AdaptiveNonlinearVariationalSolver, it doesn’t work. Currently the error is: ValueError: too many values to unpack (expected 2). I gathered that it is because I define it as a Nonlinear problem (instead of MixedNonlinear), and that’s not correct. However, I can’t figure out the right way to fix it.

This is my MWE:

```
from dolfin import *
from mshr import *
import matplotlib.pyplot as plt
import numpy as np
mesh = RectangleMesh(Point(0, 0), Point(0.2, 1), 8, 8)
n = FacetNormal(mesh)
# 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 and TestFunction(s)
w = TrialFunction(W)
(u, p, T) = split(w)
(v, q, s1) = split(TestFunction(W))
# Define the parameters and bcs
mu = Constant(1.0)
Qc2 = Constant(1.0) # this will be a function of x[0]
Ta = Constant(1.0) # this will be a function of x[1]
Pe = Constant(27.0)
Bi = Constant(11.6)
u_in = Constant(-2.0)
u_c = Constant(-1.0)
inflow = 'near(x[1], 1.0) && x[0]<=0.1'
wall = 'near(x[0], 0.2)'
centre = 'near(x[0], 0.0)'
outflow = 'near(x[1], 0.0)'
bcu_inflow = DirichletBC(W.sub(0), (0.0, u_in), inflow)
bcu_wall = DirichletBC(W.sub(0), (0.0, u_c), wall)
bcu_outflow = DirichletBC(W.sub(0), (0.0, u_c), outflow)
bcu_symmetry = DirichletBC(W.sub(0).sub(0), Constant(0.0), centre)
bcT_inflow = DirichletBC(W.sub(2), 0.0, inflow)
bcs = [bcu_inflow, bcu_wall, bcu_outflow, bcu_symmetry, bcT_inflow]
x = SpatialCoordinate(mesh)
def epsilon(v):
return sym(as_tensor([[v[0].dx(0), v[0].dx(1), 0],
[v[1].dx(0), v[1].dx(1), 0],
[0, 0, 0]]))
# stress tensor
def sigma(v, p):
return 2*mu*epsilon(v)-Id(p)
def Id(p):
return as_tensor([[p, 0, 0],
[0, p, 0],
[0, 0, p]])
def cond(v):
return sym(as_tensor([[v[0].dx(0), v[0].dx(1)],
[v[1].dx(0), v[1].dx(1)]]))
def sigmabc(v, p):
return 2*mu*cond(v) - p*Identity(2)
facet_f = MeshFunction("size_t", mesh, mesh.topology().dim() - 1) # FACET function
CompiledSubDomain('near(x[1], 1.0) && x[0]<=0.1').mark(facet_f, 0)
CompiledSubDomain('near(x[1], 1.0) && x[0]>=0.1').mark(facet_f, 1)
CompiledSubDomain('near(x[0], 0.2)').mark(facet_f, 2) # wall
CompiledSubDomain('near(x[1], 0.0)').mark(facet_f, 3) # outflow
CompiledSubDomain('near(x[0], 0.0)').mark(facet_f, 4)
CompiledSubDomain('near(x[1], 1.0) && x[0]>0.1 && x[0]<0.2').mark(facet_f, 5)
# Create the measure
ds = Measure("ds", domain=mesh, subdomain_data=facet_f)
# Define the variational form
f = Constant((0, -1))
a = (inner(sigma(u, p), epsilon(v)) - div(u) * q + (dot(u, grad(T)) * s1 + (
1 / Pe) * inner(grad(s1), grad(T)))) * dx() - (1 / Pe) * (-Bi * s1 * T * ds(2))
L = (- dot(f, v) - (1/Pe) * Qc2*s1) * dx() - (1/Pe) * (s1 * Bi * Ta * ds(2))
F = a + L # For nonlinear problem, write F = 0
w = Function(W)
# Define goal functional (quantity of interest)
M = inner(w[0], w[0])*dx()
# Define error tolerance
tol = 1.e-6
# Solve equation a = L with respect to u and the given boundary
# conditions, such that the estimated error (measured in M) is less
# than tol
problem = NonlinearVariationalProblem(F, w, bcs) # The problem is here
#problem = MixedNonlinearVariationalProblem(F, w, bcs) # This doesn't work
solver = AdaptiveNonlinearVariationalSolver(problem, M)
solver.parameters["error_control"]["dual_variational_solver"]["linear_solver"] = "umfpack"
#parameters["refinement_algorithm"] = "plaza_with_parent_facets"
solver.solve(tol)
solver.summary()
# Plot solutions
(u0, p0, T0) = w.root_node().split()
(u1, p1, T1) = w.leaf_node().split()
R = w.leaf_node().function_space()
plot(R.mesh())
```

I think the error is because of the term (dot(u, grad(T)) * s1 , but I don’t know how else to write the variational formulation. Any help would be greatly appreciated.

NOTE: I am aware that my problem for u and p is decoupled from my problem for T, so I could just solve for u and p, and then use that solution to solve for T separately. However, I am working my way up to the fully coupled problem, which is why I have written in that way.