DOLFIN

handling a term involving composition

Asked by Doug Arnold

Is it possible in Dolfin to solve a problem (given here in
1D for simplicity) like

 -u''(x) + u(x^2) = f(x)

on an interval, with boundary conditions. (For example,
if f(x) = x^4 - 2 and the BC are u(0)=0, u(1)=1, then
u(x)=x^2 is a solution.) The oddity is, of course, that the
lower order term is u(x^2), not u(x). In this simple
example, of course, there are analytic things I can do,
but the real problem is much more complicated.

Actually I want to solve a nonlinear problem, where the
lower order term is u(...u(x)...) with the argument
to u involving u(x). When I iterate this will lead
to terms like u(...uold(x)...), where uold will be
a known finite element function, and u a new finite element
function to be determined.

Is there any way to define the bilinear form "u(x^2)*v(x)*dx"
or "u(...uold(x)...)*v(x)*dx" in Dolfin?

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Marie Rognes
Solved:
Last query:
Last reply:
Revision history for this message
Best Marie Rognes (meg-simula) said : 		#1

Try the below:

# Last changed: 2011-11-23

"""
Problem:

- u''(x) + u(u(x)) = f(x)
u(0) = 0, u(1) = 1
f(x) = x**4 - 2
"""

from dolfin import *

# Turn on -O2 optimization
parameters["form_compiler"]["cpp_optimize"] = True

# Allow extrapolation from outside domain into domain (might be
# fragile for this purpose, but seems to work in this example)
parameters["allow_extrapolation"] = True

# Defining w(x) = u(u(x)) dependence here:
class Composition(Expression):
    def eval(self, values, x):
        y = self.u(x[0])
        values[0] = self.u(y)

# Define source
f = Expression("pow(x[0], 4) - 2", degree=4)

# Define mesh, function space and basis functions
mesh = UnitInterval(32)
V = FunctionSpace(mesh, "CG", 1)
du = TrialFunction(V)
v = TestFunction(V)

# Define boundary conditions
bcs = [DirichletBC(V, 0.0, "near(x[0], 0.0)"),
       DirichletBC(V, 1.0, "near(x[0], 1.0)")]

# Define initial guess for u_
u_ = Function(V)
u = Function(V)

# Define relation w(x) = u(u(x)) cf above
w = Composition(degree=2)

# Iterate: solve for u in terms of w, and update w
max_iterations = 100
tolerance = 1.e-12
for i in range(max_iterations):

    # Set previous u in w
    w.u = u_

    # Solve for u in terms of w
    F = inner(grad(du), grad(v))*dx + w*v*dx - f*v*dx
    a, L = system(F)
    solve(a == L, u, bcs)

    # Check || u - u_|| as convergence criterion
    error = errornorm(u, u_)
    print "error = ", error
    if error < tolerance:
        break

    # Update u_
    u_.assign(u)

# Define exact solution for this test case
exact = Expression("pow(x[0], 2)", degree=2)

# Compute actual error (seems to converge as n -> \infty)
actual_error = errornorm(u, exact)
print "actual_error = ", actual_error

# Inspect solution in eye-norm
plot(u, title="Approximated")
plot(exact, title="Exact", mesh=mesh)
interactive()

Revision history for this message
Doug Arnold (dnarnold) said : 		#2

Super! Thanks a lot. -- Doug

On 11/22/2011 08:35 PM, Marie Rognes wrote:
> """
> Problem:
>
> - u''(x) + u(u(x)) = f(x)
> u(0) = 0, u(1) = 1
> f(x) = x**4 - 2
> """
>
> from dolfin import *
>
> # Turn on -O2 optimization
> parameters["form_compiler"]["cpp_optimize"] = True
>
> # Allow extrapolation from outside domain into domain (might be
> # fragile for this purpose, but seems to work in this example)
> parameters["allow_extrapolation"] = True
>
> # Defining w(x) = u(u(x)) dependence here:
> class Composition(Expression):
> def eval(self, values, x):
> y = self.u(x[0])
> values[0] = self.u(y)
>
> # Define source
> f = Expression("pow(x[0], 4) - 2", degree=4)
>
> # Define mesh, function space and basis functions
> mesh = UnitInterval(32)
> V = FunctionSpace(mesh, "CG", 1)
> du = TrialFunction(V)
> v = TestFunction(V)
>
> # Define boundary conditions
> bcs = [DirichletBC(V, 0.0, "near(x[0], 0.0)"),
> DirichletBC(V, 1.0, "near(x[0], 1.0)")]
>
> # Define initial guess for u_
> u_ = Function(V)
> u = Function(V)
>
> # Define relation w(x) = u(u(x)) cf above
> w = Composition(degree=2)
>
> # Iterate: solve for u in terms of w, and update w
> max_iterations = 100
> tolerance = 1.e-12
> for i in range(max_iterations):
>
> # Set previous u in w
> w.u = u_
>
> # Solve for u in terms of w
> F = inner(grad(du), grad(v))*dx + w*v*dx - f*v*dx
> a, L = system(F)
> solve(a == L, u, bcs)
>
> # Check || u - u_|| as convergence criterion
> error = errornorm(u, u_)
> print "error = ", error
> if error< tolerance:
> break
>
> # Update u_
> u_.assign(u)
>
> # Define exact solution for this test case
> exact = Expression("pow(x[0], 2)", degree=2)
>
> # Compute actual error (seems to converge as n -> \infty)
> actual_error = errornorm(u, exact)
> print "actual_error = ", actual_error
>
> # Inspect solution in eye-norm
> plot(u, title="Approximated")
> plot(exact, title="Exact", mesh=mesh)
> interactive(

Revision history for this message
Doug Arnold (dnarnold) said : 		#3

Thanks Marie Rognes, that solved my question.