Dear all,
I am currently trying to solve a mass flow problem in 2D with an inequality constraint on the solution. The PDE it self reads:
div(f)=0;
cfunc(norm(f))f/norm(f)=-grad(u);
some boundary conditions
I was able to solve this for now (newton with relaxation parameter), and would like to further impose an inequality constraint, say:
rho = cfunc(norm(f))norm(f) between [0, 5].
According to some other threads I should be using snes as the nonlinear solver. However the attempt is not successful as the solver exits with divergence after the first iteration.
*** Error: Unable to solve nonlinear system with PETScSNESSolver.
*** Reason: Solver did not converge.
*** Where: This error was encountered inside PETScSNESSolver.cpp.
*** Process: 0
***
*** DOLFIN version: 2018.1.0
I sort of sense the reason to this problem, as the residual is increasing at the very first iteration. This is observed in Newton solver as well, because the newton step is too large to ensure descent, and is tackled by using relaxation_parameter. However, there is no such parameter for SNES solver.
I would like to ask 3 questions:
- Is there a parameter for SNES solver that resembles the relaxation_parameter in Newton solver? Or is there anyway I can manually require this?
- For Newton solver, how do I impose adaptive relaxation_parameter? i.e. I would like the stepsize to still be large for fast convergence, and decrease as the iteration # increases.
- If SNES eventual do not work out, can I impose this constraint in Newton_solver using lagrangian multiplier? The constraint is over the domain. I am not familiar with doing this in Fenics and it would be great if you can suggest a demo on this.
Below is some critical lines of the code:
sele = FiniteElement("CG", mesh.ufl_cell(), 2)
vele = VectorElement("CG", mesh.ufl_cell(), 2)
s_space = FunctionSpace(mesh,sele)
v_space = FunctionSpace(mesh,vele)
mixed_space = FunctionSpace(mesh, MixedElement([sele, vele,sele]))
U = Function(mixed_space)
u,f,rho = split(U)
V = TestFunction(mixed_space)
v,tau,p = split(V)
n = FacetNormal(mesh)
def norma(x):
return pow((pow(x[0],2)+pow(x[1],2)),0.5)
def cfunc(f):
return (alpha*pow(norma(f),gammaa)+beta*pow(norma(f),gammab))
F_gov = inner(f,grad(v))*dx+cfunc(f)*inner(f/norma(f),tau)*dx-div(tau)*u*dx\
+inner(tau,n)*u*ds\
rho_constr = (rho-norma(f)/cfunc(f))*p*dx
F = F_gov + rho_constr
J = derivative(F,U)
problem = NonlinearVariationalProblem(F, U, bc, J)
solver=NonlinearVariationalSolver(problem)
prm=solver.parameters
prm['nonlinear_solver']='snes'
prm['snes_solver']['absolute_tolerance']=1e-6
prm['snes_solver']['maximum_iterations']=2000
#imposing constraint on rho, and set the bounds for u and f(vector) to be infinity
lower = Function(mixed_space)
upper = Function(mixed_space)
ninfty = Function(s_space); ninfty.vector()[:] = -np.infty
pinfty = Function(s_space); pinfty.vector()[:] = np.infty
fa = FunctionAssigner(mixed_space, [s_space, v_space, s_space])
fa.assign(lower, [ninfty, interpolate(Constant((-np.infty,-np.infty)), v_space), interpolate(y_lower, s_space)])
fa.assign(upper, [pinfty, interpolate(Constant((-np.infty,-np.infty)), v_space), interpolate(y_upper, s_space)])
problem.set_bounds(lower, upper)
info(solver.parameters,True)
import numpy as np
U.vector().set_local(-np.random.rand(U.vector().size())) # nonzero initial guess
solver.solve()
Any help would be sincerely appreciated!
Thanks in advance.