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
Please login to add an answer/comment or follow this question.

Similar posts:
Search »