Computation of hydraulic balance-ODE fails with non-linear solver

1

Hello,

I am trying to figure out why a non-linear ordinary differential equation with analytical solution does not converge numerically (with FEM) under some conditions. This simple example could help me to understand some of the weaknesses of the newton solver for solving non-linear problems.

The equation results from a pure hydraulic balance:

\frac{du}{dt} + \sqrt{u} = f , and the variational formulation is as follows:

F=\int\limits_{\Omega}(\frac{du}{dt})v\,dx\,+\,\int\limits_{\Omega}(\sqrt{u})v\,dx\,-\,\int\limits_{\Omega}(f)v\,dx=0

There are three lines of code which may be determinant to achieve an accurate result without divergence:
Mesh definition: mesh = IntervalMesh(3,0,5)
u function initialization: u.vector()[:] = 1.0
Source therm value: f1 = Constant(0.5)

from fenics import *
import matplotlib.pyplot as plt

class BoundaryLeft(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and (near(x[0], 0.))

#Define mesh
mesh=IntervalMesh(3,0,5)

#Define function space
V=FunctionSpace(mesh,'CG',1)

#Define boundary conditions
boundaryLeft=BoundaryLeft()
u_0 = Constant(0.0)
bc = DirichletBC(V,u_0,boundaryLeft)

#Define variational problem
u = Function(V)
u.vector()[:] = 1.0

#u = interpolate(Expression("x[0] + 1.0", element=V.ufl_element()), V)
v = TestFunction(V)

f1 = Constant(0.5)
F = v*u.dx(0)*dx+v*sqrt(u)*dx-v*f1*dx

solve(F == 0, u, bc,solver_parameters={"newton_solver": {
        "relative_tolerance": 1e-7, "maximum_iterations": 50,
        "convergence_criterion": "residual",
        "linear_solver":"gmres", "preconditioner": "jacobi",
        "krylov_solver":{"nonzero_initial_guess":True}}})
plot(u)
plot(mesh)
plt.show()
error_L2 = errornorm(u_0,u,'L2')
print(error_L2)

E.g. : Changing the source therm value below 0.5 leads to solution-failure.
Is it a problem of the finite element method? Solver implementation? Should I linearize the variational problem in these cases?

Thank you in advance,
Santiago

2

You sometimes need to manually guard against nans in the variational form. Using a PETSc SNES solver with backtracking ("bt") line search can also improve robustness. For example, the following modifications run without error:

#f1 = Constant(0.5)
f1 = Constant(0.01)

#F = v*u.dx(0)*dx+v*sqrt(u)*dx-v*f1*dx
def safeSqrt(x):
    return conditional(gt(x,DOLFIN_EPS),sqrt(x),Constant(0.0))
F = v*u.dx(0)*dx+v*safeSqrt(u)*dx-v*f1*dx

#solve(F == 0, u, bc,solver_parameters={"newton_solver": {
#        "relative_tolerance": 1e-7, "maximum_iterations": 50,
#        "convergence_criterion": "residual",
#        "linear_solver":"gmres", "preconditioner": "jacobi",
#        "krylov_solver":{"nonzero_initial_guess":True}}})
solve(F==0,u,bc,
      solver_parameters
      ={"nonlinear_solver":"snes",
        "snes_solver":{"line_search":"bt"}})
1 Like
3

Hello kamensky

I thank you very much for the example; I think it can be pretty handy to avoid problems of divergence.

I still haven’t been able to find extensive information about the solvers in FEniCS and some alternatives to improve them according to the problem nature.

Is it possible to find this kind of documentation about solvers in FEniCS?, or is it the user who must discover it as he acquires experience using the software?

All the best,
Santiago

4

If you create a solver object, you can use the info() function to list all the different parameter options:

from dolfin import *

# Create some trivial nonlinear solver instance, for demonstration purposes:
mesh = UnitIntervalMesh(1)
V = FunctionSpace(mesh,"CG",1)
problem = NonlinearVariationalProblem(Function(V)*TestFunction(V)*dx,
                                      Function(V))
solver = NonlinearVariationalSolver(problem)

# Use the info() function to print out all the options for its parameters:
info(solver.parameters,True)
1 Like
5

If you want more control over the solver, you can also use solver = PETScSNESSolver() and (solve with solver.solve(problem, u.vector()), problem being an inherited class of NonlinearProblem, see the Cahn-Hilliard demo). PETSc SNES is documented here, see section 5.2 for options for the nonlinear solvers.

You can set any PETSc options with, for instance, PETScOption.set('snes_linesearch_monitor').
edit: to get the option accepted, you need then to do:

solver = PETScSNESSolver()
PETScOptions.set('snes_linesearch_monitor')    # for example some monitoring
PETScOptions.set('snes_view')                            # print the solver config
snes = solver.snes().create()
snes.setFromOptions()

solver.solve(problem, u.vector())

With solver.snes() and solver.snes().ksp you can access the underlying PETSc nonlinear and linear (ksp) solvers (petsc4py objects) and manipulate them, extract convergence history and such.

1 Like
6

Thank you very much dajuno. I haven’t tried any of those solution strategies; I will check it out!

7

Do you solve hydraulic problems by FEniCS?Is it convenient for you to share reference example?