This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

wrong result in subdomain and a plot is "not a number"

0 votes

Dear fenics-support,

i have a question :(

The following code is an example i found:

HEAT EQAUTION

from fenics import *
from dolfin import *
from dolfin_adjoint import *
from math import *
import sys
import scipy
import numpy

n = 200
mesh = RectangleMesh(Point(-1, -1), Point(1, 1), n, n)

V = FunctionSpace(mesh, "CG", 1)
W = FunctionSpace(mesh, "DG", 1)
u = Function(V, name = "State")
m = Function(W, name = "Control")
v = TestFunction(V)

Run the forward model once to create the simulation record

F = (inner(grad(u), grad(v)) - mv)dx
bc = DirichletBC(V, 0.0, "on_boundary")
solve(F == 0, u, bc)

The functional of interest is the normed difference between desired

and simulated temperature profile

x = triangle.x
u_desired = Expression("exp(-1/(1-x[0]x[0])-1/(1-x[1]x[1]))")
J = Functional((0.5inner(u-u_desired, u-u_desired))dx*dt[FINISH_TIME]) #Maybe *dx(1)

Run the optimisation

J_Tilde = ReducedFunctional(J, Control(m, value=m))

Make sure you have scipy >= 0.11 installed

m_opt = minimize(J_Tilde, method = "L-BFGS-B", tol=2e-12, bounds = (0.0, 0.5), options = {"disp": True})

F_opt = replace(F, {m: m_opt})
solve(F_opt == 0, u, bc)

plot(m_opt, mesh)
plot(mesh)
u_d_plot = interpolate(u_desired, V)
plot(u_d_plot, mesh)
plot(u, mesh)

interactive()

It runs!
I found the example called "Poisson equation with multiple subdomains" from fenics webpage.
I want to create an "great" domain and solve the problem from example "HEAT EQUATION" in an obstacle of the great domain. And in the rest of the domain should be zero!

So I have written the following code:

#

Start of Code

#
#

import all what is possible and of interest maybe

from fenics import *
from dolfin import *
from dolfin_adjoint import *
from math import *

another imports

import sys
import scipy
import numpy

#
#

subdomain

class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], -3.0)

class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 3.0)

class Bottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], -3.0)

class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 3.0)

class Obstacle(SubDomain):
def inside(self, x, on_boundary):
return (between(x[0], (-1.0, 1.0)) and between(x[1], (-1.0, 1.0)))

#
#

Initialize sub-domain instances

left = Left()
top = Top()
right = Right()
bottom = Bottom()
obstacle = Obstacle()

#
#

define mesh

n = 600
mesh = RectangleMesh(Point(-3.0, -3.0), Point(3.0, 3.0), n, n)

mesh = UnitSquareMesh(600, 600)

#
#

Initialize mesh function for interior domains

domains = CellFunction("size_t", mesh)
domains.set_all(0)
obstacle.mark(domains, 1)

#
#

Initialize mesh function for boundary domains

boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)

#
#

Define input data

a0 = Constant(1.0)
a1 = Constant(1.0)
f = Constant(0)

#
#

define functions and function spaces

V = FunctionSpace(mesh, "CG", 1)
W = FunctionSpace(mesh, "DG", 1)
u = Function(V, name = "State")
m = Function(W, name = "Control")
v = TestFunction(V)

#
#

dirichlet boundary conditions

bc = DirichletBC(V, 0.0, "on_boundary")

bc = [DirichletBC(V, 0.0, boundaries, 1),
DirichletBC(V, 0.0, boundaries, 2),
DirichletBC(V, 0.0, boundaries, 3),
DirichletBC(V, 0.0, boundaries, 4)]

#
#

Define new measures associated with the interior domains and

exterior boundaries

dx = Measure("dx")[domains]

#
#

Run the forward model once to create the simulation record

F = (inner(a0grad(u), grad(v)) - mv)dx(1) + (inner(a1grad(u), grad(v)) - fv)dx(0)
solve(F == 0, u, bc)

#
#

The functional of interest is the normed difference between desired

and simulated temperature profile

x = triangle.x
u_desired = Expression("exp(-1/(1-x[0]x[0])-1/(1-x[1]x[1]))")
J = Functional((0.5inner(u-u_desired, u-u_desired))dx*dt[FINISH_TIME])

#
#

Run the optimisation

J_Tilde = ReducedFunctional(J, Control(m, value = m))

Make sure you have scipy >= 0.11 installed

m_opt = minimize(J_Tilde, method = "L-BFGS-B", tol = 2e-12, bounds = (0.0, 0.5),
options = {"disp": True})

replace

F_opt = replace(F, {m: m_opt})
solve(F_opt == 0, u, bc)

#
#

plots

plot(m_opt, title = "control")
plot(mesh)
u_d_plot = interpolate(u_desired, V)
plot(u_d_plot, title = "u_desired")
plot(u, title = "u_real")

option for plots

interactive()

#
#

end of code

#

Anywhere must be a mistake in the code and my thinking. Because it must be the same result. But it isnt.
And also my u_desired is "nan" (not a number). What is the mistake?
And is it possible to draw a line or is it possible to make a distinction of the two domains in the plots?
Please help me to fix my problem!

Best regards

Christian
Student

asked Dec 9, 2015 by Christian FEniCS Novice (150 points)
...