How to use Adaptive Nonlinear Variational Solver with mixed nonlinear problem?

variational formulation
1

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.

2

You can resolve this by replacing

w = TrialFunction(W)

with

w = Function(W)

and then removing the line w = Function(W) that appears later in the script, above the definition of the goal functional. (You also need to uncomment

parameters["refinement_algorithm"] = "plaza_with_parent_facets"

for the adaptivity to work.)

3

Thank you so much, that does the trick.

For anyone else having these difficulties, I also neglected to write the Jacobian. In addition to the modifications suggested by @kamensky , you also need to add

dw = TrialFunction(W)

J = derivative(F, w, dw)

and define the problem with

problem = NonlinearVariationalProblem(F, w, bcs, J)

Works like a charm.