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 25 dof system to machine precision using arnoldi solver
solve > solver returned with residual ...
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 (if defined)
- 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 (if defined)
- nutils.solver.solve_linear(target, residual, *, constrain=None, lhs0=None, arguments={}, **kwargs)¶
solve linear problem
- Parameters:
target (
str
) – Name of the target: anutils.function.Argument
inresidual
.residual (
nutils.evaluable.AsEvaluableArray
) – Residual integral, depends ontarget
constrain (
numpy.ndarray
with dtypefloat
) – Defines the fixed entries of the coefficient vectorarguments (
collections.abc.Mapping
) – Defines the values fornutils.function.Argument
objects in residual. Thetarget
should not be present inarguments
. Optional.
- Returns:
Array of
target
values for whichresidual == 0
- Return type:
- nutils.solver.newton(target, residual, *, jacobian=None, lhs0=None, relax0=1.0, constrain=None, linesearch='__legacy__', 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
andD
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: anutils.function.Argument
inresidual
.residual (
nutils.evaluable.AsEvaluableArray
) –relax0 (
float
) – Initial relaxation value.constrain (
numpy.ndarray
with dtypebool
orfloat
) – Masks the free vector entries asFalse
(boolean) or NaN (float). In the remaining positions the values oflhs0
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 fornutils.function.Argument
objects in residual. Iftarget
is present inarguments
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
andF
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: anutils.function.Argument
inresidual
.residual (
nutils.evaluable.AsEvaluableArray
) –constrain (
numpy.ndarray
with dtypebool
orfloat
) – Masks the free vector entries asFalse
(boolean) or NaN (float). In the remaining positions the values oflhs0
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 fornutils.function.Argument
objects in residual. Iftarget
is present inarguments
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:
target (
str
) – Name of the target: anutils.function.Argument
inresidual
.residual (
nutils.evaluable.AsEvaluableArray
) –inertia (
nutils.evaluable.AsEvaluableArray
) –timestep (
float
) – Initial time step, will scale up as residual decreasesconstrain (
numpy.ndarray
with dtypebool
orfloat
) – Masks the free vector entries asFalse
(boolean) or NaN (float). In the remaining positions the values oflhs0
are returned unchanged (boolean) or overruled by the values in constrain (float).arguments (
collections.abc.Mapping
) – Defines the values fornutils.function.Argument
objects in residual. Iftarget
is present inarguments
then it is used as the initial guess for the iterative procedure.
- Yields:
numpy.ndarray
with dtypefloat
– Tuple of coefficient vector and residual norm
- nutils.solver.thetamethod(target, residual, inertia, timestep, theta, *, lhs0=None, target0=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:
target (
str
) – Name of the target: anutils.function.Argument
inresidual
.residual (
nutils.evaluable.AsEvaluableArray
) –inertia (
nutils.evaluable.AsEvaluableArray
) –timestep (
float
) – The time step.theta (
float
) – Theta value (theta=1 for implicit Euler, theta=0.5 for Crank-Nicolson)residual0 (
nutils.evaluable.AsEvaluableArray
) – Optional additional residual component evaluated in previous timestepconstrain (
numpy.ndarray
with dtypebool
orfloat
) – Masks the free vector entries asFalse
(boolean) or NaN (float). In the remaining positions the values oflhs0
are returned unchanged (boolean) or overruled by the values in constrain (float).newtontol (
float
) – Residual tolerance of individual timestepsarguments (
collections.abc.Mapping
) – Defines the values fornutils.function.Argument
objects in residual. Iftarget
is present inarguments
then it is used as the initial condition.timetarget (
str
) – Name of thenutils.function.Argument
that represents time. Optional.time0 (
float
) – The intial time. Default:0.0
.
- Yields:
numpy.ndarray
– Coefficient vector for all timesteps after the initial condition.
- nutils.solver.impliciteuler(target, residual, inertia, timestep, *, theta=1, lhs0=None, target0=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:
target (
str
) – Name of the target: anutils.function.Argument
inresidual
.residual (
nutils.evaluable.AsEvaluableArray
) –inertia (
nutils.evaluable.AsEvaluableArray
) –timestep (
float
) – The time step.theta (
float
) – Theta value (theta=1 for implicit Euler, theta=0.5 for Crank-Nicolson)residual0 (
nutils.evaluable.AsEvaluableArray
) – Optional additional residual component evaluated in previous timestepconstrain (
numpy.ndarray
with dtypebool
orfloat
) – Masks the free vector entries asFalse
(boolean) or NaN (float). In the remaining positions the values oflhs0
are returned unchanged (boolean) or overruled by the values in constrain (float).newtontol (
float
) – Residual tolerance of individual timestepsarguments (
collections.abc.Mapping
) – Defines the values fornutils.function.Argument
objects in residual. Iftarget
is present inarguments
then it is used as the initial condition.timetarget (
str
) – Name of thenutils.function.Argument
that represents time. Optional.time0 (
float
) – The intial time. Default:0.0
.
- Yields:
numpy.ndarray
– Coefficient vector for all timesteps after the initial condition.
- nutils.solver.cranknicolson(target, residual, inertia, timestep, *, theta=0.5, lhs0=None, target0=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:
target (
str
) – Name of the target: anutils.function.Argument
inresidual
.residual (
nutils.evaluable.AsEvaluableArray
) –inertia (
nutils.evaluable.AsEvaluableArray
) –timestep (
float
) – The time step.theta (
float
) – Theta value (theta=1 for implicit Euler, theta=0.5 for Crank-Nicolson)residual0 (
nutils.evaluable.AsEvaluableArray
) – Optional additional residual component evaluated in previous timestepconstrain (
numpy.ndarray
with dtypebool
orfloat
) – Masks the free vector entries asFalse
(boolean) or NaN (float). In the remaining positions the values oflhs0
are returned unchanged (boolean) or overruled by the values in constrain (float).newtontol (
float
) – Residual tolerance of individual timestepsarguments (
collections.abc.Mapping
) – Defines the values fornutils.function.Argument
objects in residual. Iftarget
is present inarguments
then it is used as the initial condition.timetarget (
str
) – Name of thenutils.function.Argument
that represents time. Optional.time0 (
float
) – The intial time. Default:0.0
.
- 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=None, failrelax=1e-06, **kwargs)¶
find the minimizer of a given functional
- Parameters:
target (
str
) – Name of the target: anutils.function.Argument
inresidual
.functional (scalar
nutils.evaluable.AsEvaluableArray
) – The functional the should be minimized by varying targettol (
float
) – Target residual norm.arguments (
collections.abc.Mapping
) – Defines the values fornutils.function.Argument
objects in residual. Iftarget
is present inarguments
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 dtypefloat
) – Defines the fixed entries of the coefficient vectorrelax0 (
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