### snes solver diverges for different initial guess.

68
views
0
6 weeks ago by
Howdy all,

So I'm a bit worried. Fenics defaults all functions to a value of 0, if I run the snes solver with this guess I get a good solution. My question though is when I assign the value of w to create a starting guess of .5, we get divergence of the linear solve. I'm just curious if I'm doing this correctly. What worries me is that if I set the value of w to be 0, i get it and even for negative values i.e -.1 I can also solve it (which is confusing because ln(-.1) is not defined in the real numbers. Thoughts?

from dolfin import *
import numpy as np
import math
import random

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

class Obstaclethree(SubDomain):
def inside(self, x, on_boundary):
return (((between(x[0], (-1.0, -(5.0/8.0))) or between(x[0], ((5.0/8.0),1.0 )) ) and between(x[1], (-1.0, 1.0))) or (between(x[0], (-1.0, 1.0)) and (between(x[1], (-1.0, -(5.0/8.0)) or between(x[1], ((5.0/8.0), 1.0))))))

def helmholtz_state(w,x_size,y_size,alpha_val):
#set_log_active(False)
left = Left()
top = Top()
right = Right()
bottom = Bottom()
obstacle = Obstacle()
obstacletwo = Obstacletwo()
obstaclethree = Obstaclethree()

# 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('CG', triangle, degree=1, dim=5)
V = FunctionSpace(mesh, element)
W= FunctionSpace(mesh, "CG", 1)
K = FunctionSpace(mesh,"DG",0)

# 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)
obstaclethree.mark(domains, 3)

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

w_index_list=[1]

mu_val = alpha_val
for w_val in w_index_list:

u=Function(V)
guess = Expression("0.5",degree=1)
assign(u.sub(2),interpolate(guess,V.sub(2).collapse()))
myfile = File("guess.pvd")
myfile << u.sub(2)

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

beta_one = Constant(str(w_val))

mu = Constant(str(mu_val))
L =0.5*((u_r -cos_f )**2 + (u_i - sin_f)**2)*dx(2) +.5*beta_one*(w-.5)**2*dx(1) -  mu*ln(w)*dx(1) - mu*ln(1.0-w)*dx(1)+  inner(grad(u_r),grad(lambo))*dx + ko*u_i*lambo*ds- 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 - kosquared*u_i*lambt*dx - kosquared*q*w*u_i*lambt*dx(1) - uo_i*q*w*lambt*dx(1)

first = derivative(L,u,TestFunction(V))
second = derivative(first,u,TrialFunction(V))
problem = NonlinearVariationalProblem(first, u,J=second)

lower = Function(V)
upper = Function(V)

w_lower = Constant(0.01)
w_upper = Constant(0.99)

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": {"line_search":"basic","linear_solver": "lu","maximum_iterations": 2000,"report": False,"error_on_nonconvergence": False,"absolute_tolerance":1E-6}}
solver.parameters.update(snes_solver_parameters)
print(solver.solve())

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

size=128
helmholtz_state(None,size,size,.1)​
Community: FEniCS Project

0
24 days ago by
Hello,

Finding a suitable initial guess is part of solving the system. Just for illustration purposes, the PETSc manual offers the reason for the SNES solver convergence, and it states clearly that at some cases: "All one can do at this point is change the initial guess and try again."  http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/SNES/SNESConvergedReason.html

I would suggest you to run the solver again but setting the set_log_level(50) (critical) and evaluate properly the differences on the residual of the solver based on the initial condition, for example..

Hope it helps somehow,

Regards, Leonardo