Exception: Unable to extract all indices when trying to solve a coupled PDE (Stokes + adv-diff eq)

On Tue, 30 Apr 2013 23:21:16 -0000
Noemi Petra <email address hidden> wrote:
> New question #227910 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/227910
>
> Hi,
>
> I am having difficulties solving a coupled Stokes with an adv-diff
> equation. In particular the problem comes in when trying to define
> week forms that involve both variables i.e., (u, grad(T))*s, where u
> is the velocity solution of the Stokes and T is the temperature
> solution of the adv-diff equation. I posted below a minimalistic code
> that shows the problem. When I run it I get "Exception: Unable to
> extract all indices." I am wondering if you could give a feedback or
> if you have an idea of what the problem could be? Thanks in advance!
>
> Noemi
>
> ------------------------------------------------------
> from dolfin import *
>
> #Define Mesh and Elements
> mesh = UnitSquare(16,16)
> V = VectorFunctionSpace(mesh,"Lagrange", 2)
> Q = FunctionSpace(mesh, "Lagrange", 1)
> S = FunctionSpace(mesh, "Lagrange", 1)
> ME = MixedFunctionSpace([V, Q, S])
>
> # Define Dirichlet boundary conditions
> def all_boundary(x, on_boundary):
> return on_boundary
>
> v0 = Expression(("0","0"))
> bc1 = DirichletBC(ME.sub(0), v0, all_boundary)
> bc2 = DirichletBC(ME.sub(1), 0, all_boundary)
> bc3 = DirichletBC(ME.sub(2), 0, all_boundary)
> bcs = [bc1, bc2, bc3]
>
> # body forces
> f = Expression (("1","1"))
> g = Expression(("1"))
>
> (u, p, T) = TrialFunctions(ME)
> (v, q, s) = TestFunctions(ME)
>
> # Define Variational problems
> a = inner(grad(u),grad(v))*dx - p*div(v)*dx + div(u)*q*dx
> a = a + inner(grad(T), grad(s))*dx + inner(u,grad(T))*s*dx

This is apparent non-linearity: inner(u,grad(T))*s*dx

> L = inner(f,v)*dx + g*s*dx
>
> w = Function(ME)
> solve(a == L, w, bcs, solver_parameters={'newton_solver':
> {'relative_tolerance': 1e-6}})
>

On Tue, 30 Apr 2013 23:21:16 -0000
Noemi Petra <email address hidden> wrote:
> New question #227910 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/227910
>
> Hi,
>
> I am having difficulties solving a coupled Stokes with an adv-diff
> equation. In particular the problem comes in when trying to define
> week forms that involve both variables i.e., (u, grad(T))*s, where u
> is the velocity solution of the Stokes and T is the temperature
> solution of the adv-diff equation. I posted below a minimalistic code
> that shows the problem. When I run it I get "Exception: Unable to
> extract all indices." I am wondering if you could give a feedback or
> if you have an idea of what the problem could be? Thanks in advance!

solve(a == L, ...) works with with a bilinear and L linear. But your
form a is not bilinear, it's quadratic in trial functions which is not
allowed.

To solve nonlinear problems: solve(F==0, ...) with F rank 1
(linear in test functions) form. Nonlinearity of F is covered by
arbitrary dependence on some function (not trial function). Check
corrections below.

>
> Noemi
>
> ------------------------------------------------------
> from dolfin import *
>
> #Define Mesh and Elements
> mesh = UnitSquare(16,16)
> V = VectorFunctionSpace(mesh,"Lagrange", 2)
> Q = FunctionSpace(mesh, "Lagrange", 1)
> S = FunctionSpace(mesh, "Lagrange", 1)
> ME = MixedFunctionSpace([V, Q, S])
>
> # Define Dirichlet boundary conditions
> def all_boundary(x, on_boundary):
> return on_boundary
>
> v0 = Expression(("0","0"))
> bc1 = DirichletBC(ME.sub(0), v0, all_boundary)
> bc2 = DirichletBC(ME.sub(1), 0, all_boundary)
> bc3 = DirichletBC(ME.sub(2), 0, all_boundary)
> bcs = [bc1, bc2, bc3]
>
> # body forces
> f = Expression (("1","1"))
> g = Expression(("1"))
>
> (u, p, T) = TrialFunctions(ME)

w = Function(ME)
(u, p, T) = split(w)

> (v, q, s) = TestFunctions(ME)
>
> # Define Variational problems
> a = inner(grad(u),grad(v))*dx - p*div(v)*dx + div(u)*q*dx
> a = a + inner(grad(T), grad(s))*dx + inner(u,grad(T))*s*dx
> L = inner(f,v)*dx + g*s*dx

F = a - L
J = derivative(F, w)

>
> w = Function(ME)
> solve(a == L, w, bcs, solver_parameters={'newton_solver':
> {'relative_tolerance': 1e-6}})

solve(F == 0, w, bcs, J, ...)

Ah, thanks very much. I thought the Newton solver should work with a == L as well. Yes, your suggestion did the trick. I appreciate your fast assistance. For those who are interested in such coupling, here is the complete code for the Stokes - adv-diff coupled PDE solve:

Thanks again!
Noemi

------------------
from dolfin import *

 #Define Mesh and Elements
mesh = UnitSquare(32,32)
V = VectorFunctionSpace(mesh,"Lagrange", 2)
Q = FunctionSpace(mesh, "Lagrange", 1)
S = FunctionSpace(mesh, "Lagrange", 1)
ME = MixedFunctionSpace([V, Q, S])

# Define Dirichlet boundary conditions
def all_boundary(x, on_boundary):
 return on_boundary

v0 = Expression(("0","0"))
bc1 = DirichletBC(ME.sub(0), v0, all_boundary)
bc2 = DirichletBC(ME.sub(1), 0, all_boundary)
bc3 = DirichletBC(ME.sub(2), 0, all_boundary)
bcs = [bc1, bc2, bc3]

# body forces
f = Expression (("1","1"))
g = Expression(("1"))

w = Function(ME)
(u, p, T) = split(w)
(v, q, s) = TestFunctions(ME)

# Define Variational problems
a = inner(grad(u),grad(v))*dx - p*div(v)*dx + div(u)*q*dx
a = a + inner(grad(T), grad(s))*dx + inner(u,grad(T))*s*dx
L = inner(f,v)*dx + g*s*dx
F = a - L

# Compute solution
solve(F == 0, w, bcs, solver_parameters={"newton_solver":
                                        {"relative_tolerance": 1e-6}})

# alternative way to solve the nl problem with provided Jacobian
#J = derivative(F, w)
#problem = NonlinearVariationalProblem(F,w,bcs,J)
#solver = NonlinearVariationalSolver(problem)
#solver.solve()

# Get sub-functions
u, p, T = w.split()

# Save solution in VTK format
ufile_pvd = File("velocity.pvd")
ufile_pvd << u
pfile_pvd = File("pressure.pvd")
pfile_pvd << p
Tfile_pvd = File("temperature.pvd")
Tfile_pvd << T

# Plot solution
plot(u)
plot(p)
plot(T)
interactive()