# Function with small values is turned into all 0 when call vector() function

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
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)

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

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)
# plt.colorbar(lol)
# plt.show()
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()

tao = PETSc.TAO().create(comm=None)
print(dir(tao))
#print(dir(tao))
tao.setInitial(w.vec())
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())``````

You do not copy your upper and lower bounds when you create them. You have to do a deepcopy:

``````w=w.vector()
lb=w.vec().copy()
ub=w.vec().copy()

vec_size = lb.size
for i in range(0,vec_size):
lb.setValue(i,0.0)
ub.setValue(i,10.0)

``````

You could also use the built in functionality in dolfin adjoint for optimization.
Only change you have to make to the code is to not send the jacobian into solve.
Then you could do the following:

``````    Jhat = ReducedFunctional(J, Control(w))
s_opt = minimize(Jhat, bounds=[0,1], options={"tol":1e-6, "disp":True, "maxiter":50})
File("dJ.pvd").write(Jhat.derivative())
File("optimal.pvd").write(s_opt)
print(min(s_opt.vector().get_local()), max(s_opt.vector().get_local()))

``````