Algebraic solver component in Fenics

Question

Hi, I'm trying to identify what component in Fenics is in charge of solving system of algebraic equations, i.e: when there are nonlinear terms in the PDE, the system of equations will have polynomial terms on the coefficients. Since Fenics seems to be capable of solving iteratively some nonlinear PDE, I presume this capability is lying somewhere. I want to know if this module can be used stand-alone

To give a bit more of context, I'm trying to solve steady-state heat equation on 2D surfaces with radiative terms (prop. to T^4). The mass matrix itself is relatively straightforward to compute as the powers of the node basis will be zero for nodes that are not neighbours, etc.
so, really the tricky part is the iterative solver.

I considered writing the problem as a variational expression, but the solution I'm finally after involves adding radiative terms from other surfaces, which in Fenics would be multiple SquareUnitMesh objects or their rectangular equivalents being coupled via view factor matrix terms. But I'm hoping to use the solution to the uncoupled nonlinear heat problem as an initial-guess to the fully coupled problem (i.e with cross-terms between multiple surfaces)

Edit I couldn't find a way to separate and reuse the implementation of the solver, so I decided to begin implementing my own: https://bitbucket.org/CharlesJQuarra/tensor-ais/wiki/Home

Answer 1

In FEniCS the iterative procedure of solving a non-linear problem is
by default Newton's method. Here, the problem is linearized, the
Jacobian is computed and the corresponding linear systems solved.

You may take more control over this process either by implementing your
own solver or constructing your own variant of an approximate Jacobian.
The following example shows a Picard iteration:

"""
Program to solve x**4 - 4 with Picard
The solution x = pluss/minus sqrt(2)

We assume a problem on variational form:
F(u)vdx = (a(u)u - b)vdx = 0
Then, if up is the previous solution,
we find u as:
a(up)uvdx = b
and update:
u = alphau + (1-alpha)up
where alpha is the regularization parameter
"""

from dolfin import *
import math

mesh = UnitSquareMesh(1,1)
V = FunctionSpace(mesh, "DG", 0)

u = TrialFunction(V)
U = Function(V)
U.vector()[:] = 2.7
v = TestFunction(V)

a = U**3uvdx
L = 4v*dx

F = action(a, U) - L
picard_problem = NonlinearVariationalProblem(F, U, J=a)
solver = NonlinearVariationalSolver(picard_problem)

solver.parameters["newton_solver"]["relaxation_parameter"] = 0.4
solver.parameters["newton_solver"]["maximum_iterations"] = 40
solver.parameters["newton_solver"]["absolute_tolerance"] = 1.0e-6

solver.solve()

print "The true solution sqrt(2)=%e and the numerical solution is %e " % (math.sqrt(2), max(U.vector().array()[:]))