Newton Solver Benchmark Set - Nonlinear Channel Flow#
The files in this directory
can be used to recreate the nonlinear channel flow figures from
Fraters et al. [2019].
Start by compiling the simple_nonlinear plugin in this directory as described (at the example of another benchmark)
in Running benchmarks that require code .
After compiling the plugin you can run (or adjust) the script run.sh, which will run all the models necessary
to recreate Fig. 1 and Fig. 2 of Fraters et al. [2019] and the other figures on this page
(by default the script will assume your processor has 4 compute cores, but you can adjust the number in the script).
After running the model series you can use the gnuplot files plot_v.gnuplot and plot_t.gnuplot to recreate the following figures.
Fig. 288 Convergence history for several methods for a rheology with n = 3 where in- and outflow are described by prescribing the traction. Top row: Computations in which we switch abruptly from Picard iterations to Newton iterations. Bottom row: Same as top, but using the residual scaling method in which we gradually introduce the Newton method. Left-hand column: Unmodified Newton iterations. Right-hand column: Results where we applied the SPD stabilizations to the Newton matrix. Horizontal axes: number of the non-linear (outer) iteration. Vertical axes: non-linear residual.#
Fig. 289 Convergence history for several methods for a rheology with n = 3 where in- and outflow are described by prescribing the velocity. Top row: Computations in which we switch abruptly from Picard iterations to Newton iterations. Bottom row: Same as top, but using the residual scaling method in which we gradually introduce the Newton method. Left-hand column: Unmodified Newton iterations. Right-hand column: Results where we applied the SPD stabilizations to the Newton matrix. Horizontal axes: number of the non-linear (outer) iteration. Vertical axes: non-linear residual.#
The nonlinear convergence behavior is shown in Fig. 288 and Fig. 289.
Note, that the linear residual is stalling in some cases. This is related to the choice of boundary velocity (in Fig. 289)
and boundary traction (in Fig. 288). You can adjust the base input files input_t.prm and input_v.prm to
test what happens when varying the boundary tractions or velocities.
Also note that compared to the original paper we have added a new model series to the figure that uses the default Picard solver single Advection, iterated Stokes (purple lines). You can see that this solver converges at about the same rate as the defect correction Picard solver, however its residual eventually stagnates at a fixed value. This is the expected behavior of a solver that computes the actual solution vector instead of an update to the solution vector. The accuracy of the linear solver in these cases limits the accuracy that the nonlinear solver can reach and changing the linear solver tolerance controls the exact residual at which the nonlinear solver stagnates. At the same time this illustrates the advantage of defect-correction solvers (using either Picard or Newton iterations). Even with a coarse linear solver tolerance these nonlinear solvers continue to converge to the actual solution.
Since the publication of Fraters et al. [2019] we have continued to work on improving and optimizing the Newton solver implementation and parameters. The result of this work are default parameters that can solve many models faster than using the parameters originally used. In particular two choices influence the solver convergence and are discussed here:
We can make use of a coarser maximum linear solver tolerance than the one used for the publications (here
1e-20). See the parameter Parameter: Maximum linear Stokes solver tolerance. For the example with velocity boundary conditions this change cuts the cost for the linear Stokes solver in half, while not affecting the number of nonlinear iterations significantly.Choosing less (but not 0) line search iterations actually benefits the stability of the Newton solver. This is because too many line search iterations will prevent the Newton solver from modifying the solution to a value that has a higher residual. In other words if the solver gets stuck in a local minimum there is no way for it to leave this minimum behind. Therefore by default we now limit the number of line-search iterations. See the parameter Parameter: Max Newton line search iterations for more details.
Applying these modifications here is a version of the figures above that was run with the current (as of version 3.0) default parameters of ASPECT:
Fig. 290 As above for the model with pressure boundary conditions, but using the updated default solver parameters.#
Fig. 291 As above for the model with velocity boundary conditions, but using the updated default solver parameters.#
You can easily recreate these figures by replacing mLT_1e-8 with mLT_1e-2, and by replacing LS_100 with LS_5 in the plotting scripts.