Newton solver

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@651 f88db872-0e4f-0410-b76b-b9085cfa78c5

1 files changed, 123 insertions, 0 deletions

diff --git a/src/gr/Newton.cc b/src/gr/Newton.cc
new file mode 100644
index 0000000..011f34c
--- /dev/null
+++ b/src/gr/Newton.cc
@@ -0,0 +1,123 @@
+// 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 "gfn.hh"
+#include "AHFinderDirect.hh"
+
+//******************************************************************************
+
+//
+// ***** prototypes for functions local to this file *****
+//
+
+namespace {
+	  }
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
+// This function solves H(h) = 0 for h via Newton's method.
+// At the moment the NP Jacobian is used.
+//
+void Newton_solve(patch_system& ps,
+		  const struct cactus_grid_info& cgi,
+		  const struct geometry_interpolator_info& gii,
+		  const char Jacobian_type[],
+		  fp perturbation_amplitude,
+		  int max_Newton_iterations,
+		  fp H_norm_for_convergence)
+{
+Jacobian& Jac = create_Jacobian(ps, Jacobian_type);
+
+	for (int iteration = 1 ;
+	     iteration <= max_Newton_iterations ;
+	     ++iteration)
+	{
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "Newton iteration %d", iteration);
+
+	jtutil::norm<fp> H_norms;
+	horizon_function(ps, cgi, gii, 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() <= H_norm_for_convergence)
+	   then {
+		CCTK_VInfo(CCTK_THORNSTRING, "==> finished!");
+		return;					// *** NORMAL RETURN ***
+		}
+
+	// take a Newton step
+	horizon_Jacobian_NP(ps, cgi, gii, Jac, perturbation_amplitude);
+	jtutil::array_scale(Jac.NN(),
+			    -1.0,
+			    ps.gridfn_data(nominal_gfns::gfn__H));
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "solving linear system");
+	Jac.solve_linear_system(nominal_gfns::gfn__H,		// rhs gridfn
+				nominal_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(nominal_gfns::gfn__Delta_h,
+					    irho,isigma);
+		Delta_h_norms.data(Delta_h);
+
+		p.ghosted_gridfn(ghosted_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());
+	}
+
+CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
+	   "Newton_solve: failed to converge in %d iterations!",
+	   max_Newton_iterations);
+}