### PDE Optimization : Zero Jacobian and Zero Hessian of Lagrangian Function

69
views
-1
11 weeks ago by
Hey all,

I was very curious why would I be getting a zero jacobian and zero hessian? I'm not sure if I'm displaying them right, but when I assemble and print, they appear to be 0.  When I try to solve the problem via nonlinearvariationalproblem, it explodes, which would be consistent with having a singular jacobian/ hessian.

J = [ 0.  0.  0. ...,  0.  0.  0.]

H =[[ 0.  0.  0. ...,  0.  0.  0.]

[ 0.  0.  0. ...,  0.  0.  0.]

[ 0.  0.  0. ...,  0.  0.  0.]

...,

[ 0.  0.  0. ...,  0.  0.  0.]

[ 0.  0.  0. ...,  0.  0.  0.]

[ 0.  0.  0. ...,  0.  0.  0.]]

Some insight would be much appreciated.
from dolfin import *
import matplotlib.pyplot as plt
import numpy as np
import math
from random import randint
from square import square_class
from numpy.random import rand
from petsc4py import PETSc
# from scipy import *
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)))

def helmholtz_state(w,x_size,y_size):
left = Left()
top = Top()
right = Right()
bottom = Bottom()
obstacle = Obstacle()
obstacletwo = Obstacletwo()

# Create mesh and define function space
p1 = Point(-1.0,-1.0)
p2 = Point(1.0,1.0)
mesh = RectangleMesh(p1,p2,x_size,y_size)
element = VectorElement('P', triangle, 1, dim=5)
V = FunctionSpace(mesh, element)
W= FunctionSpace(mesh, "CG", 1)

# Initialize mesh function for interior domains
domains = MeshFunction("size_t", mesh,mesh.topology().dim())
domains.set_all(0)
obstacle.mark(domains, 1)
obstacletwo.mark(domains, 2)

# Initialize mesh function for boundary domains
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)

# Define variational problem
u = Function(V)
u_r,u_i,w,lambo,lambt = split(u)
uo_r = Expression("36*pi*pi*cos(6*pi*x[1])",degree=2)
uo_i = Expression("36*pi*pi*sin(6*pi*x[1])",degree=2)
ko= Constant(6*pi)
kosquared= Constant(36*pi*pi)
q= Constant(.75)

cos_f = Expression("cos(6*pi*x[1])",degree=2)
sin_f = Expression("sin(6*pi*x[1])",degree=2)

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

L =0.5*((u_r -cos_f )**2 + (u_i - sin_f)**2)*dx(2) + inner(grad(u_r),grad(lambo))*dx + ko*u_i*lambo*ds(1)+ ko*u_i*lambo*ds(2) + ko*u_i*lambo*ds(3) + ko*u_i*lambo*ds(4) - kosquared*u_r*lambo*dx - kosquared*q*w*u_r*lambo*dx(1) - uo_r*q*w*lambo*dx(1) + inner(grad(u_i),grad(lambt))*dx - ko*u_r*lambt*ds(1)- ko*u_r*lambt*ds(2) - ko*u_r*lambt*ds(3) - ko*u_r*lambt*ds(4) - kosquared*u_i*lambt*dx - kosquared*q*w*u_i*lambt*dx(1) - uo_i*q*w*lambt*dx(1)

J = derivative(L, u)

# second derivative
H = derivative(J, u)

JMat= assemble(J)
HMat= assemble(H)
print(JMat.get_local())
print(HMat.array())

# snes_solver_parameters = {"nonlinear_solver": "snes","snes_solver": {"linear_solver": "lu","maximum_iterations": 100,"report": False,"error_on_nonconvergence": False}}
# bcs = []
# problem = NonlinearVariationalProblem(J,u,bcs,H)
# lower = Function(V)
# upper = Function(V)

# w_lower = Constant(0.0) # lower bound for state
# w_upper = Constant(1.0) # upper bound

# ninfty = Function(W);
# ninfty.vector()[:] = -np.infty
# pinfty= Function(W)
# pinfty.vector()[:] =  np.infty

# fa = FunctionAssigner(V, [W, W,W,W,W])
# fa.assign(lower, [ninfty,ninfty, interpolate(w_lower, W), ninfty,ninfty])
# fa.assign(upper, [pinfty,pinfty, interpolate(w_upper, W), pinfty,pinfty])

# #problem.set_bounds(lower, upper)

# solver  = NonlinearVariationalSolver(problem)
# solver.parameters.update(snes_solver_parameters)
# info(solver.parameters, True)

# (iter, converged) = solver.solve()

# u_r, u_i, w,lambo,lambt = u.split()

# file = File("opur" + str(x_size)+".pvd")
# file << u_r
# file = File("opw" + str(x_size)+".pvd")
# file << w

# file = File("uix.pvd")
# file << u_i

# dolfin,plot(u_r)
# plt.show()
# dolfin,plot(u_i)
# plt.show()

# dolfin.plot(w)
# plt.show()

# Plot solution
# plot(u, interactive=True)

for i in range(5,6):
print(i)
helmholtz_state(None,2**i,2**i)

#multiple of 60 for mesh.​

Community: FEniCS Project