// Newton.cc -- solve Theta(h) = 0 via Newton's method
// $Header$
//
// Newton_solve - driver to solve Theta(h) = 0 via Newton's method
//

#include <stdio.h>
#include <assert.h>
#include <math.h>

#include "util_Table.h"
#include "cctk.h"
#include "cctk_Arguments.h"

#include "config.h"
#include "stdc.h"
#include "../jtutil/util.hh"
#include "../jtutil/array.hh"
#include "../jtutil/cpm_map.hh"
#include "../jtutil/linear_map.hh"
using jtutil::error_exit;

#include "../patch/coords.hh"
#include "../patch/grid.hh"
#include "../patch/fd_grid.hh"
#include "../patch/patch.hh"
#include "../patch/patch_edge.hh"
#include "../patch/patch_interp.hh"
#include "../patch/ghost_zone.hh"
#include "../patch/patch_system.hh"

#include "../elliptic/Jacobian.hh"

#include "../gr/gfns.hh"
#include "../gr/gr.hh"

#include "driver.hh"

//******************************************************************************

//
// This function solves Theta(h) = 0 for h via Newton's method.
//
// Arguments:
// initial_find_flag = true ==> This is the first attempt to find this
//				horizon, or this is a subsequent attempt
//				and the immediately previous attempt failed.
//		       false ==> This isn't the first attempt to find this
//				 horizon, and we found it successfully on
//				 the immediately previous attempt.
// hn = the horizon number (used only in naming output files)
//
// Results:
// This function returns a success flag: this is true if (and only if)
// it found an h satisfying Theta(h) = 0 to within the error tolerances,
// otherwise it's false.
//
bool Newton_solve(patch_system& ps,
		  Jacobian& Jac,
		  const struct cactus_grid_info& cgi,
		  const struct geometry_info& gi,
		  const struct Jacobian_info& Jacobian_info,
		  const struct solver_info& solver_info,
		  bool initial_find_flag,
		  const struct IO_info& IO_info,
		  int hn, const struct verbose_info& verbose_info)
{
const int max_Newton_iterations
    = initial_find_flag ? solver_info.max_Newton_iterations__initial
			: solver_info.max_Newton_iterations__subsequent;

	for (int iteration = 1 ;
	     iteration <= max_Newton_iterations ;
	     ++iteration)
	{
	if (verbose_info.print_algorithm_details)
	   then CCTK_VInfo(CCTK_THORNSTRING,
			   "Newton iteration %d/%d",
			   iteration, max_Newton_iterations);
	if (solver_info.debugging_output_at_each_Newton_iteration)
	   then output_gridfn(ps, gfns::gfn__h,
			      IO_info, IO_info.h_base_file_name,
			      hn, verbose_info.print_algorithm_details,
			      iteration);

	//
	// evaluate the expansion Theta(h)
	//
	jtutil::norm<fp> Theta_norms;
	if (! expansion(ps,
			cgi, gi,
			true,		// yes, we want the Jacobian coeffs
			verbose_info.print_algorithm_details,
			&Theta_norms))
	   then return false;				// *** ERROR RETURN ***
	if (verbose_info.print_algorithm_highlights)
	   then CCTK_VInfo(CCTK_THORNSTRING,
			  "   iteration %d: Theta ||rms||=%.1e, ||infty||=%.1e",
			   iteration,
			   Theta_norms.rms_norm(), Theta_norms.infinity_norm());
	if (solver_info.debugging_output_at_each_Newton_iteration)
	   then output_gridfn(ps, gfns::gfn__Theta,
			      IO_info, IO_info.Theta_base_file_name,
			      hn, verbose_info.print_algorithm_details,
			      iteration);

	//
	// convergence test on ||Theta||
	//
	if (Theta_norms.infinity_norm() <= solver_info
					   .Theta_norm_for_convergence)
	   then {
		if (verbose_info.print_algorithm_details)
		   then CCTK_VInfo(CCTK_THORNSTRING,
				   "==> finished (||Theta|| < %g)",
				   double(solver_info
					  .Theta_norm_for_convergence));
		return true;				// *** NORMAL RETURN ***
		}

	//
	// compute the Newton step
	//
	if (! expansion_Jacobian(ps, Jac,
				 cgi, gi, Jacobian_info,
				 verbose_info.print_algorithm_details))
	   then return false;				// *** ERROR RETURN ***
	if (verbose_info.print_algorithm_details)
	   then CCTK_VInfo(CCTK_THORNSTRING,
			   "solving linear system for Delta_h (%d unknowns)",
			   Jac.NN());
	const fp rcond = Jac.solve_linear_system(gfns::gfn__Theta,
						 gfns::gfn__Delta_h);
	if	(rcond == 0.0)
	   then {
		CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
		   "Newton_solve: Jacobian matrix is numerically singular!");
		return false;				// *** ERROR RETURN ***
		}
	if (verbose_info.print_algorithm_details)
	   then {
		if (rcond > 0.0)
		   then CCTK_VInfo(CCTK_THORNSTRING,
			   "done solving linear system (condition number %.1e)",
				   double(rcond));
		   else CCTK_VInfo(CCTK_THORNSTRING,
				   "done solving linear system");
		}

	if (solver_info.debugging_output_at_each_Newton_iteration)
	   then output_gridfn(ps, gfns::gfn__Delta_h,
			      IO_info, IO_info.Delta_h_base_file_name,
			      hn, verbose_info.print_algorithm_details,
			      iteration);

	//
	// if the Newton step is too large, scale it down
	//
	jtutil::norm<fp> h_norms;
	ps.ghosted_gridfn_norms(gfns::gfn__h, h_norms);
	const fp max_allowable_Delta_h
		= solver_info.max_Delta_h_over_h * h_norms.mean();
	jtutil::norm<fp> Delta_h_norms;
	ps.gridfn_norms(gfns::gfn__Delta_h, Delta_h_norms);
	const fp max_Delta_h = Delta_h_norms.infinity_norm();
	const fp scale = (max_Delta_h > max_allowable_Delta_h)
			 ? max_allowable_Delta_h / max_Delta_h
			 : 1.0;

	//
	// take the Newton step (scaled if need be)
	//
		for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
		{
		patch& p = ps.ith_patch(pn);

			for (int irho = p.min_irho() ;
			     irho <= p.max_irho() ;
			     ++irho)
			{
			for (int isigma = p.min_isigma() ;
			     isigma <= p.max_isigma() ;
			     ++isigma)
			{
			p.ghosted_gridfn(gfns::gfn__h, irho,isigma)
			   -= scale * p.gridfn(gfns::gfn__Delta_h, irho,isigma);
			}
			}
		}
	if (verbose_info.print_algorithm_details)
	   then {
		if (scale == 1.0)
		   then CCTK_VInfo(CCTK_THORNSTRING,
			   "h += Delta_h (rms-norm=%.1e, infinity-norm=%.1e)",
				   Delta_h_norms.rms_norm(),
				   Delta_h_norms.infinity_norm());
		   else CCTK_VInfo(CCTK_THORNSTRING,
			   "h += %g * Delta_h (infinity-norm clamped to %.2g)",
				   scale,
				   scale * Delta_h_norms.infinity_norm());
		}

	//
	// convergence test on ||Delta_h||
	//
	if ( scale * Delta_h_norms.infinity_norm()
	     <= solver_info.Delta_h_norm_for_convergence )
	   then {
		if (verbose_info.print_algorithm_details)
		   then CCTK_VInfo(CCTK_THORNSTRING,
			   "==> finished (||Delta_h|| < %g)",
				   double(solver_info
					  .Delta_h_norm_for_convergence));
		if (solver_info.final_Theta_update_if_Delta_h_converged)
		   then {
			if (verbose_info.print_algorithm_details)
			   then CCTK_VInfo(CCTK_THORNSTRING,
					 "    doing final Theta(h) evaluation");
			if (! expansion(ps,
					cgi, gi,
					false,		// no Jacobian coeffs
					verbose_info.print_algorithm_details))
			   then return false;		// *** ERROR RETURN ***
			}
		return true;				// *** NORMAL RETURN ***
		}
	}

if (solver_info.final_Theta_update_if_Delta_h_converged)
   then {
	if (verbose_info.print_algorithm_details)
	   then CCTK_VInfo(CCTK_THORNSTRING,
			   "    doing final Theta(h) evaluation");
	if (! expansion(ps,
			cgi, gi,
			false,		// no Jacobian coeffs
			verbose_info.print_algorithm_details))
	   then return false;				// *** ERROR RETURN ***
	}
return false;						// *** ERROR RETURN ***
}