Namespace for black-box FD Newton solver. More...

Typedefs
typedef void(*	ResidualFctPt) (const Vector< double > &parameters, const Vector< double > &unknowns, Vector< double > &residuals)

typedef void(*	JacobianFctPt) (const Vector< double > &, const Vector< double > &, DenseDoubleMatrix &jacobian)

Functions
void	black_box_fd_newton_solve (ResidualFctPt residual_fct, const Vector< double > &params, Vector< double > &unknowns, JacobianFctPt jacobian_fct)

void	line_search (const Vector< double > &x_old, const double half_residual_squared_old, const Vector< double > &gradient, ResidualFctPt residual_fct, const Vector< double > &params, Vector< double > &newton_dir, Vector< double > &x, double &half_residual_squared, const double &stpmax)
	Line search helper for globally convergent Newton method. More...

Variables
unsigned	Max_iter = 20
	Max. # of Newton iterations. More...

unsigned	N_iter_taken = 0
	Number of Newton iterations taken in most recent invocation. More...

bool	Doc_Progress = false

double	FD_step = 1.0e-8
	FD step. More...

double	Tol = 1.0e-8
	Tolerance. More...

bool	Use_step_length_control = false
	Use steplength control do make globally convergent (default false) More...

Detailed Description

Namespace for black-box FD Newton solver.

/////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////

Typedef Documentation

◆ JacobianFctPt

typedef void(* oomph::BlackBoxFDNewtonSolver::JacobianFctPt) (const Vector< double > &, const Vector< double > &, DenseDoubleMatrix &jacobian)

◆ ResidualFctPt

typedef void(* oomph::BlackBoxFDNewtonSolver::ResidualFctPt) (const Vector< double > &parameters, const Vector< double > &unknowns, Vector< double > &residuals)

Function Documentation

◆ black_box_fd_newton_solve()

void oomph::BlackBoxFDNewtonSolver::black_box_fd_newton_solve	(	ResidualFctPt	residual_fct,
		const Vector< double > &	params,
		Vector< double > &	unknowns,
		JacobianFctPt	jacobian_fct
	)

Black-box FD Newton solver: Calling sequence for residual function is

residual_fct(parameters,unknowns,residuals)

where all arguments are double Vectors. unknowns.size() = residuals.size()

Reset number of Newton iterations taken in most recent invocation

     {
       // Jacobian, current and advanced residual Vectors
       unsigned ndof = unknowns.size();
       DenseDoubleMatrix jacobian(ndof);
       Vector<double> residuals(ndof);
       Vector<double> residuals_pls(ndof);
       Vector<double> dx(ndof);
       Vector<double> gradient(ndof);
       Vector<double> newton_direction(ndof);
  
       double half_residual_squared = 0.0;
       double max_step = 0.0;
  
       N_iter_taken = 0;
  
       // Newton iterations
       for (unsigned iloop = 0; iloop < Max_iter; iloop++)
       {
         // Evaluate current residuals
         residual_fct(params, unknowns, residuals);
  
         // Get half of squared residual and find maximum step length
         // for step length control
         if (Use_step_length_control)
         {
           half_residual_squared = 0.0;
           double sum = 0.0;
           for (unsigned i = 0; i < ndof; i++)
           {
             sum += unknowns[i] * unknowns[i];
             half_residual_squared += residuals[i] * residuals[i];
           }
           half_residual_squared *= 0.5;
           max_step = 100.0 * std::max(sqrt(sum), double(ndof));
         }
  
  
         // Check max. residuals
         double max_res = std::fabs(*std::max_element(
           residuals.begin(), residuals.end(), AbsCmp<double>()));
  
  
         // Doc progress?
         if (Doc_Progress)
         {
           oomph_info << "\nNewton iteration iter=" << iloop
                      << "\ni residual[i] unknown[i] " << std::endl;
           for (unsigned i = 0; i < ndof; i++)
           {
             oomph_info << i << " " << residuals[i] << " " << unknowns[i]
                        << std::endl;
           }
         }
  
  
         // Converged?
         if (max_res < Tol)
         {
           return;
         }
  
         // Next iteration...
         N_iter_taken++;
  
         // ...and how would Sir like his Jacobian?
         if (jacobian_fct == 0)
         {
           // FD loop for Jacobian
           for (unsigned i = 0; i < ndof; i++)
           {
             double backup = unknowns[i];
             unknowns[i] += FD_step;
  
             // Evaluate advanced residuals
             residual_fct(params, unknowns, residuals_pls);
  
             // Do FD
             for (unsigned j = 0; j < ndof; j++)
             {
               jacobian(j, i) = (residuals_pls[j] - residuals[j]) / FD_step;
             }
  
             // Reset fd step
             unknowns[i] = backup;
           }
         }
         // Analytical Jacobian
         else
         {
           jacobian_fct(params, unknowns, jacobian);
         }
  
  
         if (Doc_Progress)
         {
           oomph_info << "\n\nJacobian: " << std::endl;
           jacobian.sparse_indexed_output(*(oomph_info.stream_pt()));
           oomph_info << std::endl;
         }
  
         // Get gradient
         if (Use_step_length_control)
         {
           for (unsigned i = 0; i < ndof; i++)
           {
             double sum = 0.0;
             for (unsigned j = 0; j < ndof; j++)
             {
               sum += jacobian(j, i) * residuals[j];
             }
             gradient[i] = sum;
           }
         }
  
         // Solve
         jacobian.solve(residuals, newton_direction);
  
         // Update
         if (Use_step_length_control)
         {
           for (unsigned i = 0; i < ndof; i++)
           {
             newton_direction[i] *= -1.0;
           }
           // Update with steplength control
           Vector<double> unknowns_old(unknowns);
           double half_residual_squared_old = half_residual_squared;
           line_search(unknowns_old,
                       half_residual_squared_old,
                       gradient,
                       residual_fct,
                       params,
                       newton_direction,
                       unknowns,
                       half_residual_squared,
                       max_step);
         }
         else
         {
           // Direct Newton update:
           for (unsigned i = 0; i < ndof; i++)
           {
             unknowns[i] -= newton_direction[i];
           }
         }
       }
  
  
       // Failed to converge
       std::ostringstream error_stream;
       error_stream << "Newton solver did not converge in " << Max_iter
                    << " steps " << std::endl;
  
       throw OomphLibError(
         error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
     }

References Doc_Progress, boost::multiprecision::fabs(), FD_step, i, j, line_search(), max, Max_iter, N_iter_taken, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, Particles2023AnalysisHung::params, oomph::DoubleMatrixBase::solve(), oomph::Matrix< T, MATRIX_TYPE >::sparse_indexed_output(), sqrt(), oomph::OomphInfo::stream_pt(), Tol, and Use_step_length_control.

Referenced by AxisymOscillatingDisk::AxisymOscillatingDisk(), oomph::SarahBL::exact_soln(), SarahBL::exact_soln(), oomph::SarahBL::full_exact_soln(), SarahBL::full_exact_soln(), and main().

◆ line_search()

void oomph::BlackBoxFDNewtonSolver::line_search	(	const Vector< double > &	x_old,
		const double	half_residual_squared_old,
		const Vector< double > &	gradient,
		ResidualFctPt	residual_fct,
		const Vector< double > &	params,
		Vector< double > &	newton_dir,
		Vector< double > &	x,
		double &	half_residual_squared,
		const double &	stpmax
	)

Line search helper for globally convergent Newton method.

     {
       const double min_fct_decrease = 1.0e-4;
       double convergence_tol_on_x = 1.0e-16;
       double f_aux = 0.0;
       double lambda_aux = 0.0;
       double proposed_lambda = 0.0;
       unsigned n = x_old.size();
       double sum = 0.0;
       for (unsigned i = 0; i < n; i++)
       {
         sum += newton_dir[i] * newton_dir[i];
       }
       sum = sqrt(sum);
       if (sum > stpmax)
       {
         for (unsigned i = 0; i < n; i++)
         {
           newton_dir[i] *= stpmax / sum;
         }
       }
       double slope = 0.0;
       for (unsigned i = 0; i < n; i++)
       {
         slope += gradient[i] * newton_dir[i];
       }
       if (slope >= 0.0)
       {
         std::ostringstream error_stream;
         error_stream << "Roundoff problem in lnsrch: slope=" << slope << "\n";
         throw OomphLibError(
           error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
       }
       double test = 0.0;
       for (unsigned i = 0; i < n; i++)
       {
         double temp =
           std::fabs(newton_dir[i]) / std::max(std::fabs(x_old[i]), 1.0);
         if (temp > test) test = temp;
       }
       double lambda_min = convergence_tol_on_x / test;
       double lambda = 1.0;
       while (true)
       {
         for (unsigned i = 0; i < n; i++)
         {
           x[i] = x_old[i] + lambda * newton_dir[i];
         }
  
         // Evaluate current residuals
         Vector<double> residuals(n);
         residual_fct(params, x, residuals);
         half_residual_squared = 0.0;
         for (unsigned i = 0; i < n; i++)
         {
           half_residual_squared += residuals[i] * residuals[i];
         }
         half_residual_squared *= 0.5;
  
         if (lambda < lambda_min)
         {
           for (unsigned i = 0; i < n; i++) x[i] = x_old[i];
  
           // Create an Oomph Lib warning
           OomphLibWarning("Warning: Line search converged on x only!",
                           "BlackBoxFDNewtonSolver::line_search()",
                           OOMPH_EXCEPTION_LOCATION);
           return;
         }
         else if (half_residual_squared <=
                  half_residual_squared_old + min_fct_decrease * lambda * slope)
         {
           return;
         }
         else
         {
           if (lambda == 1.0)
           {
             proposed_lambda =
               -slope / (2.0 * (half_residual_squared -
                                half_residual_squared_old - slope));
           }
           else
           {
             double r1 = half_residual_squared - half_residual_squared_old -
                         lambda * slope;
             double r2 = f_aux - half_residual_squared_old - lambda_aux * slope;
             double a_poly =
               (r1 / (lambda * lambda) - r2 / (lambda_aux * lambda_aux)) /
               (lambda - lambda_aux);
             double b_poly = (-lambda_aux * r1 / (lambda * lambda) +
                              lambda * r2 / (lambda_aux * lambda_aux)) /
                             (lambda - lambda_aux);
             if (a_poly == 0.0)
             {
               proposed_lambda = -slope / (2.0 * b_poly);
             }
             else
             {
               double discriminant = b_poly * b_poly - 3.0 * a_poly * slope;
               if (discriminant < 0.0)
               {
                 proposed_lambda = 0.5 * lambda;
               }
               else if (b_poly <= 0.0)
               {
                 proposed_lambda =
                   (-b_poly + sqrt(discriminant)) / (3.0 * a_poly);
               }
               else
               {
                 proposed_lambda = -slope / (b_poly + sqrt(discriminant));
               }
             }
             if (proposed_lambda > 0.5 * lambda)
             {
               proposed_lambda = 0.5 * lambda;
             }
           }
         }
         lambda_aux = lambda;
         f_aux = half_residual_squared;
         lambda = std::max(proposed_lambda, 0.1 * lambda);
       }
     }