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

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 = VectorFunctionS

Q = FunctionSpace(mesh, "Lagrange", 1)

S = FunctionSpace(mesh, "Lagrange", 1)

ME = MixedFunctionSp

# Define Dirichlet boundary conditions

def all_boundary(x, on_boundary):

return on_boundary

v0 = Expression(

bc1 = DirichletBC(

bc2 = DirichletBC(

bc3 = DirichletBC(

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(

a = a + inner(grad(T), grad(s))*dx + inner(u,

L = inner(f,v)*dx + g*s*dx

w = Function(ME)

solve(a == L, w, bcs, solver_

## Question information

- Language:
- English Edit question

- Status:
- Solved

- For:
- DOLFIN Edit question

- Assignee:
- No assignee Edit question

- Solved by:
- Noemi Petra

- Solved:
- 2013-05-01

- Last query:
- 2013-05-01

- Last reply:
- 2013-05-01

Jan Blechta (blechta) said : | #1 |

On Tue, 30 Apr 2013 23:21:16 -0000

Noemi Petra <email address hidden> wrote:

> New question #227910 on DOLFIN:

> https:/

>

> 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 = VectorFunctionS

> Q = FunctionSpace(mesh, "Lagrange", 1)

> S = FunctionSpace(mesh, "Lagrange", 1)

> ME = MixedFunctionSp

>

> # Define Dirichlet boundary conditions

> def all_boundary(x, on_boundary):

> return on_boundary

>

> v0 = Expression(

> bc1 = DirichletBC(

> bc2 = DirichletBC(

> bc3 = DirichletBC(

> 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(

> a = a + inner(grad(T), grad(s))*dx + inner(u,

This is apparent non-linearity: inner(u,

> L = inner(f,v)*dx + g*s*dx

>

> w = Function(ME)

> solve(a == L, w, bcs, solver_

> {'relative_

>

Jan Blechta (blechta) said : | #2 |

On Tue, 30 Apr 2013 23:21:16 -0000

Noemi Petra <email address hidden> wrote:

> New question #227910 on DOLFIN:

> https:/

>

> 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 = VectorFunctionS

> Q = FunctionSpace(mesh, "Lagrange", 1)

> S = FunctionSpace(mesh, "Lagrange", 1)

> ME = MixedFunctionSp

>

> # Define Dirichlet boundary conditions

> def all_boundary(x, on_boundary):

> return on_boundary

>

> v0 = Expression(

> bc1 = DirichletBC(

> bc2 = DirichletBC(

> bc3 = DirichletBC(

> 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(

> a = a + inner(grad(T), grad(s))*dx + inner(u,

> L = inner(f,v)*dx + g*s*dx

F = a - L

J = derivative(F, w)

>

> w = Function(ME)

> solve(a == L, w, bcs, solver_

> {'relative_

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

Noemi Petra (noemi-a) said : | #3 |

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 = VectorFunctionS

Q = FunctionSpace(mesh, "Lagrange", 1)

S = FunctionSpace(mesh, "Lagrange", 1)

ME = MixedFunctionSp

# Define Dirichlet boundary conditions

def all_boundary(x, on_boundary):

return on_boundary

v0 = Expression(

bc1 = DirichletBC(

bc2 = DirichletBC(

bc3 = DirichletBC(

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(

a = a + inner(grad(T), grad(s))*dx + inner(u,

L = inner(f,v)*dx + g*s*dx

F = a - L

# Compute solution

solve(F == 0, w, bcs, solver_

# alternative way to solve the nl problem with provided Jacobian

#J = derivative(F, w)

#problem = NonlinearVariat

#solver = NonlinearVariat

#solver.solve()

# Get sub-functions

u, p, T = w.split()

# Save solution in VTK format

ufile_pvd = File("velocity.

ufile_pvd << u

pfile_pvd = File("pressure.

pfile_pvd << p

Tfile_pvd = File("temperatu

Tfile_pvd << T

# Plot solution

plot(u)

plot(p)

plot(T)

interactive()