Howdy,
So In my code I have this gradient that I calculate with dolfin adjoint. It has small entries, here I plot it
When I call the vector() function on it it makes them all zeros (which obvious when I do optimization makes it only do 1 step and call victory cause the gradient is the zero vector)
Can someone help?
from petsc4py import PETSc
from dolfin import *
import numpy as np
from fenics_adjoint import *
import math
import pylab as plt
class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], -1.0)
class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 1.0)
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], -1.0)
class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 1.0)
class Obstacle(SubDomain):
def inside(self, x, on_boundary):
return (between(x[1], (-(5.0/8.0), (5.0/8.0))) and between(x[0], (-(5.0/8.0), (5.0/8.0))))
class Obstacletwo(SubDomain):
def inside(self, x, on_boundary):
return (between(x[0], (-.5, .5)) and between(x[1], (.7, 1.0)))
class Obstaclethree(SubDomain):
def inside(self, x, on_boundary):
return (((between(x[0], (-1.0, -(5.0/8.0))) or between(x[0], ((5.0/8.0),1.0 )) ) and between(x[1], (-1.0, 1.0))) or (between(x[0], (-1.0, 1.0)) and (between(x[1], (-1.0, -(5.0/8.0)) or between(x[1], ((5.0/8.0), 1.0))))))
class ObstacleFour(SubDomain):
def inside(self, x, on_boundary):
return (between(x[0], (.7, 1.0)) and between(x[1], (.7, 1.0)))
class Circle(SubDomain):
def inside(self, x, on_boundary):
return (x[0]-.85)**2 + (x[1]-.85)**2 < pow(.1,2)
def objectivegradient(tao,X,G):
w=X
size=128
left = Left()
top = Top()
right = Right()
bottom = Bottom()
obstacle = Obstacle()
obstacletwo = Obstacletwo()
obstaclethree = Obstaclethree()
obstaclefour= ObstacleFour()
circle = Circle()
p1 = Point(-1.0,-1.0)
p2 = Point(1.0,1.0)
mesh = RectangleMesh(p1,p2,size,size)
element = VectorElement('CG', triangle, 1, dim=2)
V = FunctionSpace(mesh, element)
singlespace = FunctionSpace(mesh,'CG',1)
domains = MeshFunction("size_t", mesh,mesh.topology().dim())
domains.set_all(0)
obstacle.mark(domains, 1)
obstacletwo.mark(domains, 2)
obstaclethree.mark(domains, 3)
obstaclefour.mark(domains,4)
circle.mark(domains,5)
boundaries = MeshFunction("size_t", mesh,mesh.topology().dim()-1)
boundaries.set_all(0)
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)
dx = Measure("dx")(subdomain_data=domains)
ds = Measure("ds")(subdomain_data=boundaries)
u = Function(V)
all_states = split(u)
Lagrangian=Expression("0.0",degree=2)*dx(domain=mesh)
Remainder = Expression("0.0",degree=2)*dx(domain=mesh)
mycon=Expression("0.5",degree=2)*dx(domain=mesh)
F=Expression("0.0",degree=2)*dx(domain=mesh)
kosquared= Constant(36*pi*pi)
ko= Constant(6*pi)
q= Constant(.75)
objective=None
function_place = 0
direction = [math.cos(math.pi/4.0),math.sin(math.pi/4.0)]
print("Building Piece: " + str(i))
exp = "6*pi*(" + str(direction[0]) + "*x[0] + " + str(direction[1]) + "*x[1])"
uo_r = Expression("36*pi*pi*cos(" +str(exp) + ")",degree=2)
uo_i = Expression("36*pi*pi*sin(" +str(exp) +")",degree=2)
cos_f = Expression("cos("+str(exp) + ")",degree=2)
sin_f = Expression("sin("+str(exp) + ")",degree=2)
u_r = all_states[function_place]
u_i = all_states[function_place +1]
tests = TestFunction(V)
v_r = tests[0]
v_i =tests[1]
#wp = all_states[-1]
wp = Function(singlespace)
#w=PETScVector(w)
#print(type(wp))
wp.vector().set_local(PETScVector(w))
w=wp
domain_number=5
F=inner(grad(u_r),grad(v_r))*dx + inner(grad(u_i),grad(v_i))*dx + ko*u_i*v_r*ds(1) -ko*u_r*v_i*ds(1)+ ko*u_i*v_r*ds(2) -ko*u_r*v_i*ds(2)+ ko*u_i*v_r*ds(3) -ko*u_r*v_i*ds(3)+ ko*u_i*v_r*ds(4) -ko*u_r*v_i*ds(4) - kosquared*u_r*v_r*dx -kosquared*q*w*u_r*v_r*dx(1) -kosquared*u_i*v_i*dx -kosquared*q*w*u_i*v_i*dx(1) - uo_r*q*w*v_r*dx(1) - uo_i*q*w*v_i*dx(1)
Jac = derivative(F, u, TrialFunction(V))
solve(F == 0, u, J=Jac)
all_states = split(u)
u_r = all_states[0]
u_i = all_states[1]
# objective = .5*((u_r -cos_f)**2 + (u_i -sin_f)**2 )*dx(domain_number) + (1.0/2.0)*(w-.5)**2*dx(domain_number)
objective = .5*((u_r -cos_f)**2 + (u_i -sin_f)**2 )*dx(domain_number) + (1.0/2.0)*(w-.5)**2*dx(domain_number)
J =assemble(objective)
print(J)
gradient = compute_gradient(J,Control(w))
# lol=plot(gradient)
# plt.colorbar(lol)
# plt.show()
print(type(gradient))
gradient = gradient.vector()
print(gradient[:])
gradient = as_backend_type(gradient)
gradient = gradient.vec()
#print(gradient.view())
G=gradient
return J
def monitor(tao):
print("-------")
print("running")
print("-------")
size=128
p1 = Point(-1.0,-1.0)
p2 = Point(1.0,1.0)
mesh = RectangleMesh(p1,p2,size,size)
element = VectorElement('CG', triangle, 1, dim=5)
singlespace = FunctionSpace(mesh,'CG',1)
w=Function(singlespace)
w=w.vector()
print(type(w))
lb=w.vec()
ub=w.vec()
vec_size = lb.size
for i in range(0,vec_size):
lb.setValue(i,0.0)
ub.setValue(i,1.0)
# x = PETScVector(ub)
# v = Function(singlespace)
# v.vector().set_local(x)
# print("success")
# lol=plot(v)
# plt.colorbar(lol)
# plt.show()
#objectivegradient(w)
tao = PETSc.TAO().create(comm=None)
print(dir(tao))
#print(dir(tao))
tao.setInitial(w.vec())
tao.setObjectiveGradient(objectivegradient)
tao.setVariableBounds(lb,ub)
tao.setType('bncg')
tao.setFromOptions()
tao.setMonitor(monitor)
tao.solve()
opt_w=tao.solution
w = Function(singlespace)
opt_w = PETScVector(opt_w)
w.vector().set_local(opt_w)
print(tao.getConvergedReason())
# lol=plot(w)
# plt.colorbar(lol)
# plt.show()
#print(opt_w.view())
# print(lb.size)
# print(opt_w.size)
# print(lul.size())
# print(w.vector().size())
# print(dir(w.vector()))
# lol=plot(w)
# plt.colorbar(lol)
# plt.show()
#print(opt_w.view())
# print(dir(tao))
# print(p.view())