Solvers

Using topologies, bases and integrals, we now have the tools in place to start performing some actual functional-analytical operations. We start with what is perhaps the simplest of its kind, the least squares projection, demonstrating the different implementations now available to us and working our way up from there.

Taking the geometry component \( x_0 \) as an example, to project it onto the basis \( {φ_n} \) means finding the coefficients \( \hat{u}_n \) such that

\[ \left(\int_Ω φ_n φ_m \ dV\right) \hat u_m = \int_Ω φ_n x_0 \ dV \]

for all \( φ_n \), or \( A_{nm} \hat{u}_m = f_n \). This is implemented as follows:

A = topo.integral('basis_m basis_n dV' @ ns, degree=2).eval()
f = topo.integral('basis_n x_0 dV' @ ns, degree=2).eval()
A.solve(f)
# solve > solving 5 dof system to machine precision using arnoldi solver
# solve > solver returned with residual 3e-17±1e-15
# array([0.  , 0.25, 0.5 , 0.75, 1.  ])±1e-15

Alternatively, we can write this in the slightly more general form

\[ R_n := \int_Ω φ_n (u - x_0) \ dV = 0. \]

res = topo.integral('basis_n (u - x_0) dV' @ ns, degree=2)

Taking the derivative of \( R_n \) to \( \hat{u}m \) gives the above matrix \( A{nm} \), and substituting for \( \hat{u} \) the zero vector yields \( -f_n \). Nutils can compute those derivatives for you, using the method Array.derivative() to compute the derivative with respect to an nutils.function.Argument, returning a new nutils.function.Array.

A = res.derivative('lhs').eval()
f = -res.eval(lhs=numpy.zeros(5))
A.solve(f)
# solve > solving 5 dof system to machine precision using arnoldi solver
# solve > solver returned with residual 3e-17±1e-15
# array([0.  , 0.25, 0.5 , 0.75, 1.  ])±1e-15

The above three lines are so common that they are combined in the function nutils.solver.solve_linear:

solver.solve_linear('lhs', res)
# solve > solving 5 dof system to machine precision using arnoldi solver
# solve > solver returned with residual 3e-17±1e-15
# array([0.  , 0.25, 0.5 , 0.75, 1.  ])±1e-15

We can take this formulation one step further. Minimizing

\[ S := \int_Ω (u - x_0)^2 \ dV \]

for \( \hat{u} \) is equivalent to the above two variants. The derivative of \( S \) to \( \hat{u}_n \) gives \( 2 R_n \):

sqr = topo.integral('(u - x_0)^2 dV' @ ns, degree=2)
solver.solve_linear('lhs', sqr.derivative('lhs'))
# solve > solving 5 dof system to machine precision using arnoldi solver
# solve > solver returned with residual 6e-17±1e-15
# array([0.  , 0.25, 0.5 , 0.75, 1.  ])±1e-15

The optimization problem can also be solved by the nutils.solver.optimize function, which has the added benefit that \( S \) may be nonlinear in \( \hat{u} \) --- a property not used here.

solver.optimize('lhs', sqr)
# optimize > solve > solving 5 dof system to machine precision using arnoldi solver
# optimize > solve > solver returned with residual 0e+00±1e-15
# optimize > optimum value 0.00e+00±1e-15
# array([0.  , 0.25, 0.5 , 0.75, 1.  ])±1e-15

Nutils also supports solving a partial optimization problem. In the Laplace problem stated above, the Dirichlet boundary condition at \( Γ_\text{left} \) minimizes the following functional:

sqr = topo.boundary['left'].integral('(u - 0)^2 dS' @ ns, degree=2)

By passing the droptol argument, nutils.solver.optimize returns an array with nan ('not a number') for every entry for which the optimization problem is invariant, or to be precise, where the variation is below droptol:

cons = solver.optimize('lhs', sqr, droptol=1e-15)
# optimize > constrained 1/5 dofs
# optimize > optimum value 0.00e+00
cons
# array([ 0., nan, nan, nan, nan])±1e-15

Consider again the Laplace problem stated above. The residual is implemented as

res = topo.integral('∇_i(basis_n) ∇_i(u) dV' @ ns, degree=0)
res -= topo.boundary['right'].integral('basis_n dS' @ ns, degree=0)

Since this problem is linear in argument lhs, we can use the nutils.solver.solve_linear method to solve this problem. The constraints cons are passed via the keyword argument constrain:

lhs = solver.solve_linear('lhs', res, constrain=cons)
# solve > solving 4 dof system to machine precision using arnoldi solver
# solve > solver returned with residual 9e-16±1e-15
lhs
# array([0.  , 0.25, 0.5 , 0.75, 1.  ])±1e-15

For nonlinear residuals you can use nutils.solver.newton.

The Nutils Book

Solvers