Howdy,
So without getting into the details of PDE optimization,
I’ve solved my state equation fine - when I use my state equation and plug it into adjoint (you can think the adjoint is the same PDE in the adjoint variable, just with a different source term than the state had) I get 0 always for some reason. I basically used the collections count to show that after calling the adjoint solve, all elements are 0. When I compare this to the values in the solution to the state equation, it looks to have solved fine. Can I get some help trouble shooting this? Here is a minimal code.
from dolfin import *
import math
import collections
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 zBottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], -1.0)
class zTop(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))) and between(x[2], (-(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)) and between(x[2], (.7, 1.0)))
def direction_three(psi,theta):
return (math.sin(psi)*math.cos(theta),math.sin(psi)*math.sin(theta),math.cos(psi))
class ObstacleFour(SubDomain):
def inside(self, x, on_boundary):
return (between(x[0], (.7, 1.0)) and between(x[1], (.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 + (x[2]-.85)**2 < pow(.1,2)
set_log_level(LogLevel.ERROR)
#set_log_active(False)
left = Left()
top = Top()
right = Right()
bottom = Bottom()
topz = zTop()
bottomz=zBottom()
obstacle = Obstacle()
obstacletwo = Obstacletwo()
obstaclefour= ObstacleFour()
circle = Circle()
weights=[1.0]
p1 = Point(-1.0,-1.0,-1.0)
p2 = Point(1.0,1.0,1.0)
size=128
mesh = BoxMesh(p1,p2,size,size,int(size/16))
print("building function space")
element = VectorElement('CG', mesh.ufl_cell(), 1, dim=2)
V = FunctionSpace(mesh, element)
objfunc=FunctionSpace(mesh, "CG", 1)
dg = FunctionSpace(mesh,'DG',0)
domains = MeshFunction("size_t", mesh,mesh.topology().dim())
domains.set_all(0)
obstacle.mark(domains, 1)
obstacletwo.mark(domains, 3)
obstaclefour.mark(domains,4)
circle.mark(domains,5)
domain_number=5
dx = Measure("dx")(subdomain_data=domains)
kosquared= Constant(36*pi*pi)
ko= Constant(6*pi)
q= Constant(.75)
psi = 0.0
theta = math.pi/4.0
direction = direction_three(psi,theta)
w=Constant(1.0)
k0=6.0*math.pi
exp = str(k0)+"*(" + str(direction[0]) + "*x[0]+" + str(direction[1]) + "*x[1]+" + str(direction[2]) + "*x[2])"
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=Function(V)
(u_r,u_i) = split(u)
print(type(u))
print((type(u_r),type(u_i)))
tests = TestFunction(V)
v_r = tests[0]
v_i =tests[1]
state=inner(grad(u_r),grad(v_r))*dx + inner(grad(u_i),grad(v_i))*dx +ko*u_i*v_r*ds -ko*u_r*v_i*ds -v_r*(kosquared*(1+q*w)*u_r +kosquared*q*w*cos_f)*dx -v_i*(kosquared*(1+q*w)*u_i +kosquared*q*w*sin_f)*dx
Jac = derivative(state, u, TrialFunction(V))
solve(state == 0, u, J=Jac,solver_parameters={'newton_solver': {'linear_solver': 'mumps'}})
#u_list.append((u.sub(0,True),u.sub(1,True)))
all_states = split(u)
u_r = all_states[0]
u_i = all_states[1]
print((type(u_r),type(u_i)))
print(collections.Counter(u.sub(0,True).vector()[:]))
print(collections.Counter(u.sub(1,True).vector()[:]))
theta_counter=0
multiplers = Function(V)
tests = TestFunction(V)
v_r = tests[0]
v_i =tests[1]
print("adjoint")
print(weights[theta_counter])
print(domain_number)
(lambda_1,lambda_2) = split(multiplers)
adjoint_eq = weights[theta_counter]*(u_r - cos_f)*v_r*dx(domain_number) + weights[theta_counter]*(u_i - sin_f)*v_i*dx(domain_number) +inner(grad(v_r),grad(lambda_1))*dx +inner(grad(v_i),grad(lambda_2))*dx - lambda_1*kosquared*(1.0+q*w)*v_r*dx - lambda_2*kosquared*(1+q*w)*v_i*dx - ko*v_r*lambda_2*ds + ko*v_i*lambda_1*ds
Jac = derivative(adjoint_eq, multiplers, TrialFunction(V))
solve(adjoint_eq == 0, multiplers, J=Jac,solver_parameters={'newton_solver': {'linear_solver': 'mumps'}})
(lambda_1,lambda_2) = split(multiplers)
print(collections.Counter(multiplers.sub(0,True).vector()[:]))
print(collections.Counter(multiplers.sub(1,True).vector()[:]))