Hi everybody,
I am trying to solve a variational inequality equation which is nothing else but an elasticity problem with state constraint on the state variable, the displacement field.
In light of the “NonlinearVariationalProblem” and “NonlinearVariationalSolver”, we can solve this problem easily with snes nonlinear_solver of PETSc. What I wanna do is to set_options_prefix for the snes object associating to the “NonlinearVariationalSolver“ so that I can get the object and view the snes.
My understanding is that we are not allowed get the snes object of the “NonlinearVariationalSolver”, but we can view the snes object by setting the parameter i.e. “report” :True. We can also update parameters of the “NonlinearVariationalSolver” by using “NonlinearVariationalSolver().parameters.update()” and print these information to the screen by using “NonlinearVariationalSolver().parameters.str(True)”.
Is that all we can do to see the details of the snes_solver that is used by the NonlinearVariationalSolver?
In stead of using the “NonlinearVariationalSolver”, I am trying to use the PETScSNESSolver directly. I see that in the class PETScSNESSolver, there is a method “set_options_prefix”. However, I do not know why when I set a prefix for my PETScSNESSolver, I got this error
“AttributeError: ‘dolfin.cpp.nls.PETScSNESSolver’ object has no attribute 'set_options_prefix’”
Furthermore, I want to solve the problem using the snes-vi method, but I could not set_bounds for the problem as well as solving it using the method “vinewtonrsls”.
However, I do not know how to set bounds for the solver PETScSNESSolver. I do see that PETScSNESSolver allow us to call with bounds in the input arguments, but when I tried to run the problem with bounds and vinewtonrsls method, I got this error
“TypeError: solve(): incompatible function arguments. The following argument types are supported:
- (self: dolfin.cpp.nls.PETScSNESSolver, arg0: dolfin::NonlinearProblem, arg1: dolfin.cpp.la.GenericVector) -> Tuple[int, bool]”
Could anyone explain to me if I understand these solvers correctly or not and please give me some advice? Thank you so much.
I am using dolfin for python version 2018.1.0.
I am sorry. I do not know how to attach a python file. A minimum working example without state inequality constraint on the displacement field is below.
**"""
The program is to solve a variational inequality problem for elasticity.
The linear elastic problem:
-div((h + delta)^p * sigma(u)) = f
u = uD on the boundary.
sigma (u) . n = g on the Neumann boundary
delta > 0 i.e. delta:= 1e-8
“”"
from future import print_function
import matplotlib.pyplot as plt
from fenics import *
from dolfin import *
import numpy as np
import sys, os, sympy, math
from mshr import *
import petsc4py as PETSc
set_log_level(1000)
try:
from petsc4py import PETSc
except:
print ("\n*** Warning: you need to have petsc4py installed for this module to run\n")
print ("\nExiting.\n")
exit()
if not has_petsc4py():
print ("\n*** Warning: Dolfin is not compiled with petsc4py support\n")
print ("\nExiting.\n")
exit()
def main():
# Optimization options for the form compiler
parameters[“form_compiler”][“cpp_optimize”] = True
parameters[“form_compiler”][“representation”] = “uflacs”
parameters[“linear_algebra_backend”] = “PETSc”
global delta, ninfty, pinfty;
delta = 1e-8; # a small number to prevent singularity
pinfty = np.inf; ninfty = -np.inf;
h0Val = 1.; p = 4.; q = 2.;
ubVal_x = 1000.; lbVal_x = -1000; ubVal_y = 1000.; lbVal_y = -1000;
Lx = 1.; Ly = .1; Nx = 150; Ny = 30;
fVal = 0. ; f = Constant((0., -fVal)) # body force
# Traction vectors on the Neumann B.C.
gL = Constant((0., 0.)); gT = Constant((0.0, -1.))
gR = Constant((0., 0.)); gB = Constant((0.0, 0.0))
bcs = []
# Elasticity parameters. Note: scale u by changing E
E, nu = Constant(1e3), Constant(0.3)
mu, lambda_ = E/2/(1 + nu), E*nu/(1 + nu)/(1 - 2*nu)
model = "plane_stress"
if model == "plane_stress":
lambda_ = 2*mu*lambda_/(lambda_+2*mu)
savedir = "QA_snes_solver/Lx%s_Ly%s_Nx%s_Ny%s"%(Lx,Ly,Nx,Ny)
# Create mesh and define functional spaces
#------------------------------------
mesh = RectangleMesh(Point(0.,0.), Point(Lx, Ly), Nx, Ny)
Vu = VectorFunctionSpace(mesh, 'P', 2)
Vh = FunctionSpace(mesh, 'P', 1)
# Define the boundaries
class Left(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0., DOLFIN_EPS)
class Top(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], Ly, DOLFIN_EPS)
class Right(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], Lx, DOLFIN_EPS)
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 0., DOLFIN_EPS)
sub_left = Left(); sub_right = Right();
sub_top = Top(); sub_bottom = Bottom();
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim() - 1) # "size_t", not "double"
boundary_markers.set_all(4)
sub_left.mark(boundary_markers, 0) # marked as 0
sub_top.mark(boundary_markers, 1) # marked as 1
sub_right.mark(boundary_markers, 2) # marked as 2
sub_bottom.mark(boundary_markers, 3) # marked as 3
# Redefine boundary integration measure
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
# Dirichlet B.C.
bcl = DirichletBC(Vu, Constant((0.0, 0.0)), sub_left); bcs = [bcl]
# Define variational problem for PDE
# -------------------------------------------
u_1= Function(Vu)
u_2 = Function(Vu)
h = Function(Vh)
d = u_1.geometric_dimension() # space dimension
# Define strain and stress
def epsilon(u):
return sym(grad(u))
def sigma(u):
return lambda_*tr(epsilon(u))*Identity(d) + 2.0*mu*epsilon(u)
# Initial guess for h
h0 = Expression('h0Val', degree=1, h0Val = h0Val)
h.assign(interpolate(h0, Vh))
ubVal_Ex = Expression(("xmax","ymax"), degree=1, xmax=ubVal_x, ymax=ubVal_y)
lbVal_Ex = Expression(("xmin","ymin - x[1]"), degree=1, xmin=lbVal_x, ymin=lbVal_y)
lb_u = interpolate(lbVal_Ex, Vu); ub_u = interpolate(ubVal_Ex, Vu)
outputFile = File(savedir + "/ub_u.pvd"); outputFile << ub_u
outputFile = File(savedir + "/lb_u.pvd"); outputFile << lb_u
# Solve the PDE problem using Snes solver
u_1 = nls_solver(u_1, h, f, p, bcs, ub_u, lb_u, gL, gT, gR, gB, lambda_, mu, ds)
u_2 = snes_solver(u_2, h, f, p, bcs, ub_u, lb_u, gL, gT, gR, gB, lambda_, mu, ds)
print("\n The difference between the two method is \n", errornorm(u_1, u_2))
# Save solution to file in VTK format
File(savedir + "/thickness.pvd") << h
return
def nls_solver(u, h, f, p, bcs, ub, lb, gL, gT, gR, gB, lambda_, mu, ds):
Vu = u.function_space(); mesh = Vu.mesh();
# Define variational problem
du = TrialFunction(Vu); v = TestFunction(Vu)
d = u.geometric_dimension() # space dimension
# Define strain and stress
def epsilon(u):
return sym(grad(u))
def sigma(u):
return lambda_*tr(epsilon(u))*Identity(d) + 2*mu*epsilon(u)
# The weak form
PDE = (h**p + delta)* inner(sigma(u), epsilon(v))*dx - dot(f, v)*dx - dot(gL, v)*ds(0) - dot(gT, v)*ds(1)- dot(gR, v)*ds(2) - dot(gB, v)*ds(3)
# Differentiate to get Jacobian.
PDE_u = derivative(PDE, u, du)
snes_solver_parameters_bounds = {"nonlinear_solver": "snes", "snes_solver": {"linear_solver": "umfpack", "maximum_iterations": 30, "preconditioner": "none", "report": True, "line_search": "basic","sign": "default","method":"vinewtonrsls", "absolute_tolerance":1e-8, "relative_tolerance":1e-8, "solution_tolerance":1e-8, "krylov_solver": {"report" : True, "monitor_convergence" : True,"absolute_tolerance":1e-10,"relative_tolerance":1e-10}}}
# Set up the non-linear problem
problem_pde = NonlinearVariationalProblem(PDE, u, bcs, J=PDE_u)
problem_pde.set_bounds(lb.vector(), ub.vector())
nls_solver = NonlinearVariationalSolver(problem_pde)
nls_solver.parameters.update(snes_solver_parameters_bounds)
print("parameters of the nls_solver is \n", nls_solver.parameters.str(True))
(iter, converged) = nls_solver.solve()
if not converged:
warning("Convergence is not guaranteed when modifying some parameters or using PETSc.")
return u
def snes_solver(u, h, f, p, bcs, ub, lb, gL, gT, gR, gB, lambda_, mu, ds):
Vu = u.function_space(); mesh = Vu.mesh();
# Define variational problem
du = TrialFunction(Vu); v = TestFunction(Vu)
d = u.geometric_dimension() # space dimension
# Define strain and stress
def epsilon(u):
return sym(grad(u))
def sigma(u):
return lambda_*tr(epsilon(u))*Identity(d) + 2*mu*epsilon(u)
# The weak form
PDE = (h**p + delta)* inner(sigma(u), epsilon(v))*dx - dot(f, v)*dx - dot(gL, v)*ds(0) - dot(gT, v)*ds(1)- dot(gR, v)*ds(2) - dot(gB, v)*ds(3)
class GeneralProblem(NonlinearProblem):
def __init__(self, F, u, bcs):
NonlinearProblem.__init__(self)
self.fform = F
self.u = u
self.bcs = bcs
self.jacobian = derivative(F, u)
def F(self, b, x):
assemble(self.fform, tensor=b)
[bc.apply(b, x) for bc in self.bcs]
def J(self, A, x):
assemble(self.jacobian, tensor=A)
[bc.apply(A) for bc in self.bcs]
problem_pde = GeneralProblem(PDE, u, bcs)
my_snes_solver = PETScSNESSolver()
my_snes_solver.set_options_prefix(“mypde”)
solver_parameters = {"linear_solver": "umfpack", "maximum_iterations": 30, "preconditioner": "none", "report": True, "line_search": "basic","sign": "default","method":"newtonls", "absolute_tolerance":1e-8, "relative_tolerance":1e-8, "solution_tolerance":1e-8, "krylov_solver": {"report" : True, "monitor_convergence" : True,"absolute_tolerance":1e-10,"relative_tolerance":1e-10, "error_on_nonconvergence" : True, "maximum_iterations": 1000}, "lu_solver" : {"report": True, "symmetric": False, "verbose": False}}
# method: vinewtonrsls, newtonls,
my_snes_solver.parameters.update(solver_parameters)
snes = my_snes_solver.snes()
print("\n view the snes: Start\n")
snes.view()
print("\n view the snes: End\n")
print("\n, In the snes_solver, the parameters of snes_solver are \n", my_snes_solver.parameters.str(True))
(iter, converged) = my_snes_solver.solve(problem_pde, u.vector(), lb.vector(), ub.vector())
(iter, converged) = my_snes_solver.solve(problem_pde, u.vector())
if not converged:
warning("Convergence is not guaranteed when modifying some parameters or using PETSc.")
return u
if name == ‘main’:
main()
**