// Newton.cc -- solve H(h) = 0 via Newton's method
// $Id$
//
// <<<prototypes for functions local to this file>>>
// Newton_solve - driver to solve H(h) = 0 via Newton's method
//

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

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

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

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

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

#include "gfns.hh"
#include "AHFinderDirect.hh"

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

//
// ***** prototypes for functions local to this file *****
//

namespace {
	  }

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

//
// This function solves H(h) = 0 for h via Newton's method.
//
bool Newton_solve(patch_system& ps,
		  const struct cactus_grid_info& cgi,
		  const struct geometry_info& gi,
		  const struct Jacobian_info& Jacobian_info,
		  const struct solver_info& solver_info,
		  Jacobian& Jac)
{
	for (int iteration = 1 ;
	     iteration <= solver_info.max_Newton_iterations ;
	     ++iteration)
	{
	CCTK_VInfo(CCTK_THORNSTRING,
		   "Newton iteration %d", iteration);

	jtutil::norm<fp> H_norms;
	horizon_function(ps, cgi, gi, true, true, &H_norms);
	CCTK_VInfo(CCTK_THORNSTRING,
		   "   H rms-norm=%.2e, infinity-norm=%.2e",
		   H_norms.rms_norm(), H_norms.infinity_norm());

	if (H_norms.infinity_norm() <= solver_info.H_norm_for_convergence)
	   then {
		CCTK_VInfo(CCTK_THORNSTRING,
			   "==> finished (||H|| < %g)",
			   double(solver_info.H_norm_for_convergence));
		return true;				// *** NORMAL RETURN ***
		}

	// take a Newton step
	horizon_Jacobian(ps, cgi, gi, Jacobian_info, Jac);
	CCTK_VInfo(CCTK_THORNSTRING,
		   "solving linear system for %d unknowns", Jac.NN());
	Jac.solve_linear_system(gfns::gfn__H,		// rhs gridfn
				gfns::gfn__Delta_h);	// soln gridfn
	CCTK_VInfo(CCTK_THORNSTRING,
		   "done solving linear system");
	jtutil::norm<fp> Delta_h_norms;
		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)
		{
		const fp Delta_h = p.gridfn(gfns::gfn__Delta_h, irho,isigma);
		Delta_h_norms.data(Delta_h);

		p.ghosted_gridfn(gfns::gfn__h, irho,isigma) -= Delta_h;
		}
		}
		}
	CCTK_VInfo(CCTK_THORNSTRING,
		   "moved h by Delta_h (rhs-norm=%.2e, infinity-norm=%.2e)",
		   Delta_h_norms.rms_norm(), Delta_h_norms.infinity_norm());

	if (Delta_h_norms.infinity_norm() <= solver_info
					     .Delta_h_norm_for_convergence)
	   then {
		CCTK_VInfo(CCTK_THORNSTRING,
			   "==> finished (||Delta_h|| < %g)",
			   double(solver_info.Delta_h_norm_for_convergence));
		return true;				// *** NORMAL RETURN ***
		}
	}

CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
	   "Newton_solve: no convergence in %d iterations!\n",
	   solver_info.max_Newton_iterations);
return false;						// *** ERROR RETURN ***
}