From 5deb07627f88cfa85e5912735d5fb1bd55364871 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Mon, 6 Nov 2017 13:21:06 +0100 Subject: [PATCH] Performance improvement for transfinite interpolation manifold. Pass initial point along to the push_forward_gradient call rather than re-computing it. Use quasi-Newton method (Broyden's method) rather than full Newton with finite differences to reduce number of calls to the compute_transfinite_interpolation. --- include/deal.II/grid/manifold_lib.h | 22 +++++--- source/grid/manifold_lib.cc | 86 +++++++++++++++++++++++------ 2 files changed, 84 insertions(+), 24 deletions(-) diff --git a/include/deal.II/grid/manifold_lib.h b/include/deal.II/grid/manifold_lib.h index fd8d8ccf42..0cf05a34b5 100644 --- a/include/deal.II/grid/manifold_lib.h +++ b/include/deal.II/grid/manifold_lib.h @@ -780,12 +780,12 @@ private: * Push forward operation. * * @note This internal function is currently not compatible with the - * ChartManifold::pull_back() function because the given class represents an - * atlas of charts, not a single chart. Thus, the pull_back() operation is - * only valid with the additional information of the chart, given by a cell - * on the coarse grid. An alternative implementation could shift the index - * depending on the coarse cell for a 1-to-1 relation between the chart space - * and the image space. + * ChartManifold::push_forward() function because the given class represents + * an atlas of charts, not a single chart. Thus, the push_forward() + * operation is only valid with the additional information of the chart, + * given by a cell on the coarse grid. An alternative implementation could + * shift the index depending on the coarse cell for a 1-to-1 relation + * between the chart space and the image space. */ Point push_forward(const typename Triangulation::cell_iterator &cell, @@ -793,10 +793,18 @@ private: /** * Gradient of the push_forward method. + * + * @note This internal function is not compatible with the + * ChartManifold::pull_back() function because the given class represents an + * atlas of charts, not a single chart. Furthermore, this private function + * also requires the user to provide the result of the push_forward() call + * on the chart point for the use case of this function, namely inside a + * Newton iteration where the gradient is computed by finite differences. */ DerivativeForm<1,dim,spacedim> push_forward_gradient(const typename Triangulation::cell_iterator &cell, - const Point &chart_point) const; + const Point &chart_point, + const Point &pushed_forward_chart_point) const; /** * The underlying triangulation. diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc index 44a0a1c503..ebac5f2764 100644 --- a/source/grid/manifold_lib.cc +++ b/source/grid/manifold_lib.cc @@ -1107,6 +1107,8 @@ namespace } // next subtract the contributions of the lines + const auto weights_view_line = make_array_view(weights.begin(), weights.begin()+2); + const auto points_view_line = make_array_view(points.begin(), points.begin()+2); for (unsigned int line=0; line::lines_per_cell; ++line) { const double line_point = (line < 8 ? chart_point[1-(line%4)/2] : chart_point[2]); @@ -1142,8 +1144,8 @@ namespace weights[0] = 1. - line_point; weights[1] = line_point; new_point -= my_weight * - cell.line(line)->get_manifold().get_new_point(points_view, - weights_view); + cell.line(line)->get_manifold().get_new_point(points_view_line, + weights_view_line); } } @@ -1181,11 +1183,10 @@ template DerivativeForm<1,dim,spacedim> TransfiniteInterpolationManifold ::push_forward_gradient(const typename Triangulation::cell_iterator &cell, - const Point &chart_point) const + const Point &chart_point, + const Point &pushed_forward_chart_point) const { // compute the derivative with the help of finite differences - Point point = compute_transfinite_interpolation(*cell, chart_point, - coarse_cell_is_flat[cell->index()]); DerivativeForm<1,dim,spacedim> grad; for (unsigned int d=0; d modified[d] += step; Tensor<1,spacedim> difference = compute_transfinite_interpolation(*cell, modified, - coarse_cell_is_flat[cell->index()]) - point; + coarse_cell_is_flat[cell->index()]) - + pushed_forward_chart_point; for (unsigned int e=0; e Point chart_point = GeometryInfo::project_to_unit_cell(cell->real_to_unit_cell_affine_approximation(point)); - // run Newton iteration. As opposed to the various mapping implementations, - // this class does not throw exception upon failure as those are relatively - // expensive and failure occurs quite regularly in the implementation of the - // compute_chart_points method. + // run quasi-Newton iteration with a combination of finite differences for + // the exact Jacobian and "Broyden's good method". As opposed to the various + // mapping implementations, this class does not throw exception upon failure + // as those are relatively expensive and failure occurs quite regularly in + // the implementation of the compute_chart_points method. Tensor<1,spacedim> residual = point - compute_transfinite_interpolation(*cell, chart_point, coarse_cell_is_flat[cell->index()]); const double tolerance = 1e-21 * cell->diameter() * cell->diameter(); double residual_norm_square = residual.norm_square(); + DerivativeForm<1,dim,spacedim> inv_grad; for (unsigned int i=0; i<100; ++i) { if (residual_norm_square < tolerance) - return chart_point; + { + // do a final update of the point with the last available Jacobian + // information. The residual is close to zero due to the check + // above, but me might improve some of the last digits by a final + // Newton-like step with step length 1 + Tensor<1,dim> update; + for (unsigned int d=0; d grad = push_forward_gradient(cell, chart_point); - if (grad.determinant() <= 0.0) - return outside; - DerivativeForm<1,dim,spacedim> inv_grad = grad.covariant_form(); + // every 8 iterations, including the first time around, we create an + // approximation of the Jacobian with finite differences. Broyden's + // method usually does not need more than 5-8 iterations, but sometimes + // we might have had a bad initial guess and then we can accelerate + // convergence considerably with getting the actual Jacobian rather than + // using secant-like methods (one gradient calculation in 3D costs as + // much as 3 more iterations). this usually happens close to convergence + // and one more step with the finite-differenced Jacobian leads to + // convergence + if (i%8 == 0) + { + // if the determinant is zero, the mapping is not invertible and we are + // outside the valid chart region + DerivativeForm<1,dim,spacedim> grad + = push_forward_gradient(cell, chart_point, + Point(point-residual)); + if (grad.determinant() <= 0.0) + return outside; + inv_grad = grad.covariant_form(); + } Tensor<1,dim> update; for (unsigned int d=0; d alpha > 1e-7) alpha *= 0.5; + const Tensor<1,spacedim> old_residual = residual; while (alpha > 1e-7) { Point guess = chart_point + alpha*update; @@ -1270,6 +1299,29 @@ TransfiniteInterpolationManifold } if (alpha < 1e-7) return outside; + + // update the inverse Jacobian with "Broyden's good method" and + // Sherman-Morrison formula for the update of the inverse, see + // https://en.wikipedia.org/wiki/Broyden%27s_method + const Tensor<1,dim> delta_x = alpha*update; + + // switch sign in residual as compared to the wikipedia article because + // we use a negative definition of the residual with respect to the + // Jacobian + const Tensor<1,spacedim> delta_f = old_residual - residual; + + Tensor<1,dim> Jinv_deltaf; + for (unsigned int d=0; d factor = (delta_x - Jinv_deltaf)/(delta_x * Jinv_deltaf); + Tensor<1,spacedim> jac_update; + for (unsigned int d=0; d