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
