dolfin team mailing list archive

[Marie Rognes <meg@xxxxxxxxx>]

Achim Schroll wrote:
> Question #105989 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/105989
>
>     Status: Answered => Open
>
> Achim Schroll is still having a problem:
> I still appreciate any help to set up function spaces to automatically differentiate and solve a nonlinear system of the following type:
> dot(h*grad(s), grad(phi))*dx + f*phi*dx = 0
> dot(s*grad(h), grad(psi))*dx + g*psi*dx = 0
> Unknowns are u=(s,h), while (f,g) are given.
> Best, Achim
>
>   


Example attached. (The example is just made up, so the Newton solver 
does not converge and so the answer is just rubbish, but it should 
illustrate how to go about.)

You are probably particularily interested in line 25.



--
Marie

from dolfin import *

# Define mesh
n = 8
mesh = UnitSquare(n, n)

# Define sources
f = Expression("sin(pi*x[0])")
g = Expression("2.0")

# Define function spaces at will
V = FunctionSpace(mesh, "CG", 1)
Q = FunctionSpace(mesh, "CG", 1)
M = V * Q # Older versions M = V + Q

# Define boundary condition(s)
boundary = lambda x, on_boundary: on_boundary
bcs = [DirichletBC(M.sub(0), Constant(1.0), boundary),
       DirichletBC(M.sub(1), Constant(0.0), boundary)]

# Define function for solution
u = Function(M)

# Take components for easy form-writing. (Do NOT use split.)
(s, h) = as_vector((u[0], u[1]))

# Define test functions
(phi, psi) = TestFunctions(M)

# Define (nonlinear) form
F = (dot(h*grad(s), grad(phi)) + f*phi + dot(s*grad(h), grad(psi)) + g*psi)*dx

# Differentiate F with respect to u in the direction of du:
du = TrialFunction(M)
dF = derivative(F, u, du)

# Define non-linear variational problem (dF(v, du; u) = - F(v; u))
pde = VariationalProblem(dF, F, bcs, nonlinear=True)
(s, h) = pde.solve().split()

# Plot solution
plot(s)
plot(h)
interactive()