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

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