dolfin plot with colorbar - white patches in solution (PDE constrained optimization problem)


94
views
0
10 weeks ago by
Howdy folks,

So when I plot one of my solutions : I get these white patches  - which seems to be some sort of weird out of the colorbar range. I'm doing a constrained optimization problem so on this solution w, I have 0 <= w <=1



I was curious if someone knew how to adjust the color bar so I don't have these patches? Is there just some issue with my code that I'm not seeing?

from dolfin import *
import matplotlib.pyplot as plt
import numpy as np
import math
from random import randint
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)
	u=Function(W)

	# 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
	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)
	u=Function(V)
	u_r,u_i,w,lambo,lambt = split(u)


	w_list = list(range(100,101))
	lagrangian = []

	
	alpha_one = Constant(".0001")
	alpha_two = Constant(".0001")
	beta_one = Constant("1.0")
	L =0.5*((u_r -cos_f )**2 + (u_i - sin_f)**2)*dx(2) + beta_one*(w**2)*dx +  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)
	obj = 0.5*((u_r -cos_f )**2 + (u_i - sin_f)**2)*dx(2)
	first = derivative(L,u,TrialFunction(V))
	second = derivative(first,u)
	bc=[]
	problem = NonlinearVariationalProblem(first, u, bc, J=second)
	
	lower = Function(V)
	upper = Function(V)

	w_lower = Constant(0.0)
	w_upper = Constant(1.0) 

	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)
	snes_solver_parameters = {"nonlinear_solver": "snes","snes_solver": {"linear_solver": "lu","maximum_iterations": 100,"report": False,"error_on_nonconvergence": False}}
	solver.parameters.update(snes_solver_parameters)
	solver.solve()


	r=dolfin.plot(w)
	plt.colorbar(r)
	plt.show()




for i in range(5,6):
	print(i)
	helmholtz_state(None,2**i,2**i)
Community: FEniCS Project
Use paraview or something sensible. Check for NaN in your solution vector.
written 10 weeks ago by Nate  
Thank you Nate for the suggestion. That helped a lot. Everything is fine - seems to be just a plotting issue in plt (these white spots are not in the paraview and everything is between 0 and 1). Thank you.
written 10 weeks ago by Ryan Vogt  

1 Answer


0
10 weeks ago by
Just to answer the question - seems like there is some sort of issue with the plt plot. When I just project w, and plot in paraview everything looks great.


As with my code I just write this at the bottom after the solving i.e

problem.set_bounds(lower, upper)
	solver  = NonlinearVariationalSolver(problem)
	snes_solver_parameters = {"nonlinear_solver": "snes","snes_solver": {"linear_solver": "lu","maximum_iterations": 100,"report": False,"error_on_nonconvergence": False}}
	solver.parameters.update(snes_solver_parameters)
	solver.solve()
	u_r,u_i,w,lambo,lambt = split(u)
	w = project(w, W)
	file = File("wopt.pvd")
	file << w
Please login to add an answer/comment or follow this question.

Similar posts:
Search »