solver

The solver module defines solvers for problems of the kind res = 0 or ∂inertia/∂t + res = 0, where res is a nutils.evaluable.AsEvaluableArray. To demonstrate this consider the following setup:

>>> from nutils import mesh, function, solver
>>> ns = function.Namespace()
>>> domain, ns.x = mesh.rectilinear([4,4])
>>> ns.basis = domain.basis('spline', degree=2)
>>> cons = domain.boundary['left,top'].project(0, onto=ns.basis, geometry=ns.x, ischeme='gauss4')
project > constrained 11/36 dofs, error 0.00e+00/area
>>> ns.u = 'basis_n ?lhs_n'

Function u represents an element from the discrete space but cannot not evaluated yet as we did not yet establish values for ?lhs. It can, however, be used to construct a residual functional res. Aiming to solve the Poisson problem u_,kk = f we define the residual functional res = v,k u,k + v f and solve for res == 0 using solve_linear:

>>> res = domain.integral('(basis_n,i u_,i + basis_n) d:x' @ ns, degree=2)
>>> lhs = solver.solve_linear('lhs', residual=res, constrain=cons)
solve > solving for argument lhs (36) using direct method
solve > residual norm: 1.5e-15

The coefficients lhs represent the solution to the Poisson problem.

In addition to solve_linear the solver module defines newton and pseudotime for solving nonlinear problems, as well as impliciteuler for time dependent problems.

class nutils.solver.NormBased(minscale=0.01, acceptscale=0.6666666666666666, maxscale=2.0)

Bases: object

Line search abstraction for Newton-like iterations, computing relaxation values that correspond to greatest reduction of the residual norm.

Parameters:

  • minscale (float) – Minimum relaxation scaling per update. Must be strictly greater than zero.

  • acceptscale (float) – Relaxation scaling that is considered close enough to optimality to to accept the current Newton update. Must lie between minscale and one.

  • maxscale (float) – Maximum relaxation scaling per update. Must be greater than one, and therefore always coincides with acceptance, determining how fast relaxation values rebound to one if not bounded by optimality.

__weakref__

list of weak references to the object

class nutils.solver.MedianBased(minscale=0.01, acceptscale=0.6666666666666666, maxscale=2.0, quantile=0.5)

Bases: object

Line search abstraction for Newton-like iterations, computing relaxation values such that half (or any other configurable quantile) of the residual vector has its optimal reduction beyond it. Unline the NormBased approach this is invariant to constant scaling of the residual items.

Parameters:

  • minscale (float) – Minimum relaxation scaling per update. Must be strictly greater than zero.

  • acceptscale (float) – Relaxation scaling that is considered close enough to optimality to to accept the current Newton update. Must lie between minscale and one.

  • maxscale (float) – Maximum relaxation scaling per update. Must be greater than one, and therefore always coincides with acceptance, determining how fast relaxation values rebound to one if not bounded by optimality.

  • quantile (float) – Fraction of the residual vector that is aimed to have its optimal reduction at a smaller relaxation value. The default value of one half corresponds to the median. A value close to zero means tighter control, resulting in strong relaxation.

__weakref__

list of weak references to the object

class nutils.solver.System(residual, /, trial, test=None)

Bases: object

System of one or more variables.

The System class represents a linear or nonlinear problem of one or more variables, and offers several methods to solve it. The main solution method is solve, which serves as a general entry point that offloads the problem to specialized solvers. The other prominent method is step, which adds functionality to deal with time dependent problems.

__init__(self, residual, /, trial, test=None)

Construct a system.

Parameters:

  • residual – Any object that supports conversion to a scalar evaluable Array via the AsEvaluableArray protocol. This is the function that the trial arguments will seek to make constant in all the test arguments, or zero in all the derivatives. If provided as a tuple of vectors then these are already the derivatives to be made zero, and no test arguments may be specified.

  • trial – Names of the trial arguments, provided as a string (‘arg1,arg2’) or tuple of strings ((‘arg1’, ‘arg2’)).

  • test – Names of the test arguments (optional) provided as a string or tuple of strings. Must be omitted if residual is a tuple of vectors; otherwise, may be omitted if test and trial are equal.

solve(self, *, arguments={}, constrain={}, tol=0.0, miniter=0, maxiter=None, method=None)

Solve the system.

Determines the trial arguments for which the derivatives of the system’s scalar functional are zero in all the test arguments.

Parameters:

  • arguments – Arguments required to evaluate the system. Arguments that coincide with the system’s trial arguments serve as initial guess.

  • constrain – Constrained values for any of the trial arguments, supplied as float vectors with NaN entries or boolean vectors, the latter holding values at those of the initial guess.

  • tol – Maximum residual norm, stopping criterion for the iterative solver. Required for nonlinear problems.

  • miniter – Minimum number of iterations for the iterative solver, overruling to tol based stopping criterion.

  • maxiter – Maximum number of iterations, after which a SolverError is raised.

  • method – Iterative solution method, in the form of a callable that returns an iterator of (arguments, resnorm, value) triplets.

step(self, *, arguments, suffix, timearg=None, timesteparg=None, timestep=None, maxretry=2, **solveargs)

Advance a time step.

This method is best described by an example. Let timearg equal ‘t’ and suffix ‘0’, and let our system’s trial arguments be ‘u’ and ‘v’. This method then creates argument copies ‘t0’, ‘u0’, ‘v0’, advances ‘t’ by timestep, and solves for the new ‘u’ and ‘v’.

Parameters:

  • arguments – Arguments required to evaluate the system. Arguments that coincide with the system’s trial arguments serve as initial value.

  • suffix – String suffix to add to argument names to denote their value at the beginning of the time step.

  • timearg – Name of the scalar argument that tracks the time (optional).

  • timesteparg – Name of the scalar argument that tracks the timestep (optional).

  • timestep – Size of the time increment (required if either timearg or timesteparg are specified).

  • maxretry – If either timearg or timesteparg are specified and affecting the system, then this positive integer determines how many levels of timestep bisections will be considered to recover from solver or matrix errors.

  • solveargs – Remaining keyword arguments are passed on to the solver method.

solve_constraints(self, *, droptol, arguments={}, constrain={}, linargs={})

Solve for Dirichlet constraints.

This method is similar to solve, but with two key differences.

The method is limited to linear systems, but adds the ability to solve a limited class of singular systems. It does so by isolating the subset of arguments that contribute (up to droptol) to the residual, and solving the corresponding submatrix. The remaining argument values are returned as NaN (not a number).

The second key difference with solve is that the returned dictionary is augmented with the remaining _constrain_ items, rather than those from _arguments_, reflecting the method’s main utility of forming Dirichlet constraints. This allows for the aggregation of constraints by calling the method multiple times in series.

Parameters:

  • arguments – Arguments required to evaluate the system. Any arguments that coincide with the system’s trial arguments serve as initial guess.

  • constrain – Constrained values for any of the trial arguments, supplied as float vectors with NaN entries or boolean vectors, the latter holding values at those of the initial guess.

  • linargs – Keyword arguments to be passed on to the linear solver.

  • droptol – Minimum absolute value if matrix entries for rows and columns to participate in the optimization problem.

__weakref__

list of weak references to the object

class nutils.solver.Direct(**linargs)

Bases: object

Direct solution of a linear system.

__weakref__

list of weak references to the object

class nutils.solver.Newton(**linargs)

Bases: object

Newton solver.

This solver generates the iterates x_n+1 = x_n - J(x_n)^-1 r(x_n) consistent with a standard Newton process.

__weakref__

list of weak references to the object

class nutils.solver.ReuseNewton(require=0.5, **linargs)

Bases: object

Newton with approximate Jacobians.

This is a modification of the vanilla Newton process, in which a new Jacobian matrix is NOT created for every update, but an existing matrix is reused instead for as long as |res(x - J^-1 res(x))| < require |res(x)|.

When the condition no longer holds, a new Jacobian matrix is assembled and a standard Newton update is performed without regard for the require parameter.

__weakref__

list of weak references to the object

class nutils.solver.LinesearchNewton(strategy=NormBased(minscale=0.01, acceptscale=0.6666666666666666, maxscale=2.0), failrelax=1e-06, relax0=1.0, **linargs)

Bases: object

Newton solver with automatic relaxation.

This solver generates the iterates x_n+1 = x_n - α J(x_n)^-1 r(x_n), where 0 < α <= 1 is a relaxation value. The value for α is initially relax0, and subsequently determined by a configurable strategy. If the value falls below failrelax then the solver fails with a SolverError.

The default strategy is NormBased, which employs a polynomial interpolation of the residual norm between the current value of x and the prospective next value to determine at what value of α the residual norm is likely minimal.

__weakref__

list of weak references to the object

class nutils.solver.Minimize(rampup=0.5, rampdown=-1.0, failrelax=-10.0, **linargs)

Bases: object

Direct minimization of scalar functional.

This solver constructs a trajectory of steepest descent to find the minimizer of a locally convex scalar function.

__weakref__

list of weak references to the object

class nutils.solver.Arnoldi(maxiter=2, **linargs)

Bases: object

Arnoldi iterations for linear systems.

Upon first use, the Arnoldi solver behaves identically to Direct, in that it returns the solution to the linear system A x = b.

The difference arises when the method is reused. If the matrix and residual depend on a second set of arguments, as in A(y) x = b(y), then the Arnoldi solver will aim to benefit from an existing factorization of P = A(y_previous)^-1 by projecting x on the growing subspace [P b(y), P^2 b(y), P^3 b(y), ..).

If the desired tolerance is not reached at a configurable subspace dimension of maxiter, then the system is solved with A(y) and the cached matrix updated for subsequent reuse.

__weakref__

list of weak references to the object

class nutils.solver.Pseudotime(inertia, timestep, **linargs)

Bases: object

Inertia assisted Newton.

This solver is a variant of the Newton solver for problems that do not converge along the standard Newton path. For problems that represent the steady state solution of a dynamic system, the pseudo-time solver uses a provided inertia matrix to progress towards the steady state along the (approximately) physical path using ever larger time steps.

Concretely, rather than forming updates of the form J(x) dx = r(x), the pseudo-time solver solves (J(x) + dJ / dt) dx = r(x), where dJ is the provided inertia matrix and dt is the configured initial timestep scaled by the ratio of |r(x0)| / |r(x)|.

__weakref__

list of weak references to the object

nutils.solver.solve_linear(target, residual, *, constrain=None, lhs0=None, arguments={}, **kwargs)

solve linear problem

Parameters:
Returns:

Array of target values for which residual == 0

Return type:

numpy.ndarray

nutils.solver.newton(target, residual, *, jacobian=None, lhs0=None, relax0=1.0, constrain=None, linesearch=NormBased(minscale=0.01, acceptscale=0.6666666666666666, maxscale=2.0), failrelax=1e-06, arguments={}, **kwargs)

iteratively solve nonlinear problem by gradient descent

Generates targets such that residual approaches 0 using Newton procedure with line search based on the residual norm. Suitable to be used inside solve.

An optimal relaxation value is computed based on the following cubic assumption:

|res(lhs + r * dlhs)|^2 = A + B * r + C * r^2 + D * r^3

where A, B, C and D are determined based on the current residual and tangent, the new residual, and the new tangent. If this value is found to be close to 1 then the newton update is accepted.

Parameters:

  • target (str) – Name of the target: a nutils.function.Argument in residual.

  • residual (nutils.evaluable.AsEvaluableArray)

  • relax0 (float) – Initial relaxation value.

  • constrain (numpy.ndarray with dtype bool or float) – Masks the free vector entries as False (boolean) or NaN (float). In the remaining positions the values of lhs0 are returned unchanged (boolean) or overruled by the values in constrain (float).

  • linesearch (Callable[[float, float, float, float], Tuple[float, bool]]) – Callable that defines relaxation logic. The callable takes four arguments: the current residual and directional derivative, and the candidate residual and directional derivative, with derivatives normalized to unit length; and returns the optimal scaling and a boolean flag that marks whether the candidate should be accepted.

  • failrelax (float) – Fail with exception if relaxation reaches this lower limit.

  • arguments (collections.abc.Mapping) – Defines the values for nutils.function.Argument objects in residual. If target is present in arguments then it is used as the initial guess for the iterative procedure.

Yields:

numpy.ndarray – Coefficient vector that approximates residual==0 with increasing accuracy

nutils.solver.minimize(target, energy, *, lhs0=None, constrain=None, rampup=0.5, rampdown=-1.0, failrelax=-10.0, arguments={}, **kwargs)

iteratively minimize nonlinear functional by gradient descent

Generates targets such that residual approaches 0 using Newton procedure with line search based on the energy. Suitable to be used inside solve.

An optimal relaxation value is computed based on the following assumption:

energy(lhs + r * dlhs) = A + B * r + C * r^2 + D * r^3 + E * r^4 + F * r^5

where A, B, C, D, E and F are determined based on the current and new energy, residual and tangent. If this value is found to be close to 1 then the newton update is accepted.

Parameters:

  • target (str) – Name of the target: a nutils.function.Argument in residual.

  • residual (nutils.evaluable.AsEvaluableArray)

  • constrain (numpy.ndarray with dtype bool or float) – Masks the free vector entries as False (boolean) or NaN (float). In the remaining positions the values of lhs0 are returned unchanged (boolean) or overruled by the values in constrain (float).

  • rampup (float) – Value to increase the relaxation power by in case energy is decreasing.

  • rampdown (float) – Value to decrease the relaxation power by in case energy is increasing.

  • failrelax (float) – Fail with exception if relaxation reaches this lower limit.

  • arguments (collections.abc.Mapping) – Defines the values for nutils.function.Argument objects in residual. If target is present in arguments then it is used as the initial guess for the iterative procedure.

Yields:

numpy.ndarray – Coefficient vector that approximates residual==0 with increasing accuracy

nutils.solver.pseudotime(target, residual, inertia, timestep, *, lhs0=None, constrain=None, arguments={}, **kwargs)

iteratively solve nonlinear problem by pseudo time stepping

Generates targets such that residual approaches 0 using hybrid of Newton and time stepping. Requires an inertia term and initial timestep. Suitable to be used inside solve.

Parameters:
Yields:

numpy.ndarray with dtype float – Tuple of coefficient vector and residual norm

nutils.solver.thetamethod(target, residual, inertia, timestep, theta, *, lhs0=None, constrain=None, newtontol=1e-10, arguments={}, newtonargs={}, timetarget='_thetamethod_time', time0=0.0, historysuffix='0')

solve time dependent problem using the theta method

Parameters:
Yields:

numpy.ndarray – Coefficient vector for all timesteps after the initial condition.

nutils.solver.impliciteuler(target, residual, inertia, timestep, *, theta=1, lhs0=None, constrain=None, newtontol=1e-10, arguments={}, newtonargs={}, timetarget='_thetamethod_time', time0=0.0, historysuffix='0')

solve time dependent problem using the theta method

Parameters:
Yields:

numpy.ndarray – Coefficient vector for all timesteps after the initial condition.

nutils.solver.cranknicolson(target, residual, inertia, timestep, *, theta=0.5, lhs0=None, constrain=None, newtontol=1e-10, arguments={}, newtonargs={}, timetarget='_thetamethod_time', time0=0.0, historysuffix='0')

solve time dependent problem using the theta method

Parameters:
Yields:

numpy.ndarray – Coefficient vector for all timesteps after the initial condition.

nutils.solver.optimize(target, functional, *, tol=0.0, arguments={}, droptol=None, constrain=None, lhs0=None, relax0=1.0, linesearch=NormBased(minscale=0.01, acceptscale=0.6666666666666666, maxscale=2.0), failrelax=1e-06, **kwargs)

find the minimizer of a given functional

Parameters:

  • target (str) – Name of the target: a nutils.function.Argument in residual.

  • functional (scalar nutils.evaluable.AsEvaluableArray) – The functional the should be minimized by varying target

  • tol (float) – Target residual norm.

  • arguments (collections.abc.Mapping) – Defines the values for nutils.function.Argument objects in residual. If target is present in arguments then it is used as the initial guess.

  • droptol (float) – Threshold for leaving entries in the return value at NaN if they do not contribute to the value of the functional.

  • constrain (numpy.ndarray with dtype float) – Defines the fixed entries of the coefficient vector

  • relax0 (float) – Initial relaxation value.

  • linesearch (Callable[[float, float, float, float], Tuple[float, bool]]) – Callable that defines relaxation logic. The callable takes four arguments: the current residual and directional derivative, and the candidate residual and directional derivative, with derivatives normalized to unit length; and returns the optimal scaling and a boolean flag that marks whether the candidate should be accepted.

  • failrelax (float) – Fail with exception if relaxation reaches this lower limit.

Yields:

numpy.ndarray – Coefficient vector corresponding to the functional optimum