aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--README13
-rw-r--r--doc/documentation.tex144
-rw-r--r--interface.ccl6
-rw-r--r--param.ccl73
-rw-r--r--schedule.ccl7
-rw-r--r--src/CoordTransf.c159
-rw-r--r--src/Equations.c186
-rw-r--r--src/FuncAndJacobian.c744
-rw-r--r--src/Newton.c435
-rw-r--r--src/TP_utilities.c951
-rw-r--r--src/TP_utilities.h91
-rw-r--r--src/TwoPunctures.c128
-rw-r--r--src/TwoPunctures.h95
-rw-r--r--src/make.code.defn8
-rw-r--r--test/twopunctures.par53
15 files changed, 3093 insertions, 0 deletions
diff --git a/README b/README
new file mode 100644
index 0000000..5180e00
--- /dev/null
+++ b/README
@@ -0,0 +1,13 @@
+CVS info : $Header$
+
+Cactus Code Thorn TwoPunctures
+Thorn Author(s) : Marcus Ansorg <marcus.ansorg@aei.mpg.de>
+ : Erik Schnetter <schnetter@aei.mpg.de>
+Thorn Maintainer(s) : Marcus Ansorg <marcus.ansorg@aei.mpg.de>
+ : Erik Schnetter <schnetter@aei.mpg.de>
+--------------------------------------------------------------------------
+
+Purpose of the thorn:
+
+Create initial for two puncture black holes using a single domain
+spectral method.
diff --git a/doc/documentation.tex b/doc/documentation.tex
new file mode 100644
index 0000000..60ff50f
--- /dev/null
+++ b/doc/documentation.tex
@@ -0,0 +1,144 @@
+% *======================================================================*
+% Cactus Thorn template for ThornGuide documentation
+% Author: Ian Kelley
+% Date: Sun Jun 02, 2002
+% $Header$
+%
+% Thorn documentation in the latex file doc/documentation.tex
+% will be included in ThornGuides built with the Cactus make system.
+% The scripts employed by the make system automatically include
+% pages about variables, parameters and scheduling parsed from the
+% relevant thorn CCL files.
+%
+% This template contains guidelines which help to assure that your
+% documentation will be correctly added to ThornGuides. More
+% information is available in the Cactus UsersGuide.
+%
+% Guidelines:
+% - Do not change anything before the line
+% % START CACTUS THORNGUIDE",
+% except for filling in the title, author, date, etc. fields.
+% - Each of these fields should only be on ONE line.
+% - Author names should be separated with a \\ or a comma.
+% - You can define your own macros, but they must appear after
+% the START CACTUS THORNGUIDE line, and must not redefine standard
+% latex commands.
+% - To avoid name clashes with other thorns, 'labels', 'citations',
+% 'references', and 'image' names should conform to the following
+% convention:
+% ARRANGEMENT_THORN_LABEL
+% For example, an image wave.eps in the arrangement CactusWave and
+% thorn WaveToyC should be renamed to CactusWave_WaveToyC_wave.eps
+% - Graphics should only be included using the graphicx package.
+% More specifically, with the "\includegraphics" command. Do
+% not specify any graphic file extensions in your .tex file. This
+% will allow us to create a PDF version of the ThornGuide
+% via pdflatex.
+% - References should be included with the latex "\bibitem" command.
+% - Use \begin{abstract}...\end{abstract} instead of \abstract{...}
+% - Do not use \appendix, instead include any appendices you need as
+% standard sections.
+% - For the benefit of our Perl scripts, and for future extensions,
+% please use simple latex.
+%
+% *======================================================================*
+%
+% Example of including a graphic image:
+% \begin{figure}[ht]
+% \begin{center}
+% \includegraphics[width=6cm]{MyArrangement_MyThorn_MyFigure}
+% \end{center}
+% \caption{Illustration of this and that}
+% \label{MyArrangement_MyThorn_MyLabel}
+% \end{figure}
+%
+% Example of using a label:
+% \label{MyArrangement_MyThorn_MyLabel}
+%
+% Example of a citation:
+% \cite{MyArrangement_MyThorn_Author99}
+%
+% Example of including a reference
+% \bibitem{MyArrangement_MyThorn_Author99}
+% {J. Author, {\em The Title of the Book, Journal, or periodical}, 1 (1999),
+% 1--16. {\tt http://www.nowhere.com/}}
+%
+% *======================================================================*
+
+% If you are using CVS use this line to give version information
+% $Header$
+
+\documentclass{article}
+
+% Use the Cactus ThornGuide style file
+% (Automatically used from Cactus distribution, if you have a
+% thorn without the Cactus Flesh download this from the Cactus
+% homepage at www.cactuscode.org)
+\usepackage{../../../../doc/latex/cactus}
+
+\begin{document}
+
+% The author of the documentation
+\author{Marcus Ansorg \textless marcus.ansorg@aei.mpg.de\textgreater \\ Erik Schnetter \textless schnetter@aei.mpg.de\textgreater}
+
+% The title of the document (not necessarily the name of the Thorn)
+\title{}
+
+% the date your document was last changed, if your document is in CVS,
+% please use:
+% \date{$ $Date$ $}
+\date{May 18 2004}
+
+\maketitle
+
+% Do not delete next line
+% START CACTUS THORNGUIDE
+
+% Add all definitions used in this documentation here
+% \def\mydef etc
+
+% Add an abstract for this thorn's documentation
+\begin{abstract}
+
+\end{abstract}
+
+% The following sections are suggestive only.
+% Remove them or add your own.
+
+\section{Introduction}
+
+\section{Physical System}
+
+\section{Numerical Implementation}
+
+\section{Using This Thorn}
+
+\subsection{Obtaining This Thorn}
+
+\subsection{Basic Usage}
+
+\subsection{Special Behaviour}
+
+\subsection{Interaction With Other Thorns}
+
+\subsection{Examples}
+
+\subsection{Support and Feedback}
+
+\section{History}
+
+\subsection{Thorn Source Code}
+
+\subsection{Thorn Documentation}
+
+\subsection{Acknowledgements}
+
+
+\begin{thebibliography}{9}
+
+\end{thebibliography}
+
+% Do not delete next line
+% END CACTUS THORNGUIDE
+
+\end{document}
diff --git a/interface.ccl b/interface.ccl
new file mode 100644
index 0000000..427ac85
--- /dev/null
+++ b/interface.ccl
@@ -0,0 +1,6 @@
+# Interface definition for thorn TwoPunctures
+# $Header$
+
+IMPLEMENTS: TwoPunctures
+
+INHERITS: ADMBase StaticConformal grid
diff --git a/param.ccl b/param.ccl
new file mode 100644
index 0000000..0d006df
--- /dev/null
+++ b/param.ccl
@@ -0,0 +1,73 @@
+# Parameter definitions for thorn TwoPunctures
+# $Header$
+
+SHARES: ADMBase
+
+USES KEYWORD metric_type
+
+EXTENDS KEYWORD initial_data
+{
+ "twopunctures" :: "two puncture black holes"
+}
+
+
+
+SHARES: StaticConformal
+
+USES KEYWORD conformal_storage
+
+
+
+PRIVATE:
+
+INT npoints_A "Number of coefficients in the compactified radial direction"
+{
+ 4:* :: ""
+} 30
+
+INT npoints_B "Number of coefficients in the angular direction"
+{
+ 4:* :: ""
+} 30
+
+INT npoints_phi "Number of coefficients in the phi direction"
+{
+ 4:*:2 :: ""
+} 16
+
+
+
+REAL par_b "x coordinate of the m+ puncture"
+{
+ (0.0:*) :: ""
+} 1.0
+
+REAL par_m_plus "mass of the m+ puncture"
+{
+ (0.0:* :: ""
+} 1.0
+
+REAL par_m_minus "mass of the m- puncture"
+{
+ (0.0:* :: ""
+} 1.0
+
+REAL par_P_plus[3] "momentum of the m+ puncture"
+{
+ (*:*) :: ""
+} 0.0
+
+REAL par_P_minus[3] "momentum of the m- puncture"
+{
+ (*:*) :: ""
+} 0.0
+
+REAL par_S_plus[3] "spin of the m+ puncture"
+{
+ (*:*) :: ""
+} 0.0
+
+REAL par_S_minus[3] "spin of the m- puncture"
+{
+ (*:*) :: ""
+} 0.0
diff --git a/schedule.ccl b/schedule.ccl
new file mode 100644
index 0000000..f625a55
--- /dev/null
+++ b/schedule.ccl
@@ -0,0 +1,7 @@
+# Schedule definitions for thorn TwoPunctures
+# $Header$
+
+SCHEDULE TwoPunctures IN ADMBase_InitialData
+{
+ LANG: C
+} "Create puncture black hole initial data"
diff --git a/src/CoordTransf.c b/src/CoordTransf.c
new file mode 100644
index 0000000..5b2f31a
--- /dev/null
+++ b/src/CoordTransf.c
@@ -0,0 +1,159 @@
+// TwoPunctures: File "CoordTransf.c"
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <ctype.h>
+#include <time.h>
+#include "cctk_Parameters.h"
+#include "TP_utilities.h"
+#include "TwoPunctures.h"
+
+//---------------------------------------------------------------
+void
+AB_To_XR (int nvar, double A, double B, double *X, double *R,
+ derivs U)
+// On Entrance: U.d0[]=U[]; U.d1[] =U[]_A; U.d2[] =U[]_B; U.d3[] =U[]_3;
+// U.d11[]=U[]_AA; U.d12[]=U[]_AB; U.d13[]=U[]_A3;
+// U.d22[]=U[]_BB; U.d23[]=U[]_B3; U.d33[]=U[]_33;
+// At Exit: U.d0[]=U[]; U.d1[] =U[]_X; U.d2[] =U[]_R; U.d3[] =U[]_3;
+// U.d11[]=U[]_XX; U.d12[]=U[]_XR; U.d13[]=U[]_X3;
+// U.d22[]=U[]_RR; U.d23[]=U[]_R3; U.d33[]=U[]_33;
+{
+ DECLARE_CCTK_PARAMETERS;
+ double At = 0.5 * (A + 1), A_X, A_XX, B_R, B_RR;
+ int ivar;
+
+ *X = 2 * atanh (At);
+ *R = Pih + 2 * atan (B);
+
+ A_X = 1 - At * At;
+ A_XX = -At * A_X;
+ B_R = 0.5 * (1 + B * B);
+ B_RR = B * B_R;
+
+ for (ivar = 0; ivar < nvar; ivar++)
+ {
+ U.d11[ivar] = A_X * A_X * U.d11[ivar] + A_XX * U.d1[ivar];
+ U.d12[ivar] = A_X * B_R * U.d12[ivar];
+ U.d13[ivar] = A_X * U.d13[ivar];
+ U.d22[ivar] = B_R * B_R * U.d22[ivar] + B_RR * U.d2[ivar];
+ U.d23[ivar] = B_R * U.d23[ivar];
+ U.d1[ivar] = A_X * U.d1[ivar];
+ U.d2[ivar] = B_R * U.d2[ivar];
+ }
+}
+
+//---------------------------------------------------------------
+void
+C_To_c (int nvar, double X, double R, double *x, double *r,
+ derivs U)
+// On Entrance: U.d0[]=U[]; U.d1[] =U[]_X; U.d2[] =U[]_R; U.d3[] =U[]_3;
+// U.d11[]=U[]_XX; U.d12[]=U[]_XR; U.d13[]=U[]_X3;
+// U.d22[]=U[]_RR; U.d23[]=U[]_R3; U.d33[]=U[]_33;
+// At Exit: U.d0[]=U[]; U.d1[] =U[]_x; U.d2[] =U[]_r; U.d3[] =U[]_3;
+// U.d11[]=U[]_xx; U.d12[]=U[]_xr; U.d13[]=U[]_x3;
+// U.d22[]=U[]_rr; U.d23[]=U[]_r3; U.d33[]=U[]_33;
+{
+ DECLARE_CCTK_PARAMETERS;
+ double C_c2, U_cb, U_CB;
+ dcomplex C, C_c, C_cc, c, c_C, c_CC, U_c, U_cc, U_C, U_CC, One =
+ Complex (1., 0.);
+ int ivar;
+
+ C.r = X;
+ C.i = R;
+
+ c = RCmul (par_b, Ccosh (C)); // c=b*cosh(C)
+ c_C = RCmul (par_b, Csinh (C));
+ c_CC = c;
+
+ C_c = Cdiv (One, c_C);
+ C_cc = RCmul (-1., Cmul (Cmul (C_c, C_c), Cmul (C_c, c_CC)));
+ C_c2 = C_c.r * C_c.r + C_c.i * C_c.i;
+
+ for (ivar = 0; ivar < nvar; ivar++)
+ {
+ // U_C = 0.5*(U_X3-i*U_R3)
+ // U_c = U_C*C_c = 0.5*(U_x3-i*U_r3)
+ U_C.r = 0.5 * U.d13[ivar];
+ U_C.i = -0.5 * U.d23[ivar];
+ U_c = Cmul (U_C, C_c);
+ U.d13[ivar] = 2. * U_c.r;
+ U.d23[ivar] = -2. * U_c.i;
+
+ // U_C = 0.5*(U_X-i*U_R)
+ // U_c = U_C*C_c = 0.5*(U_x-i*U_r)
+ U_C.r = 0.5 * U.d1[ivar];
+ U_C.i = -0.5 * U.d2[ivar];
+ U_c = Cmul (U_C, C_c);
+ U.d1[ivar] = 2. * U_c.r;
+ U.d2[ivar] = -2. * U_c.i;
+
+ // U_CC = 0.25*(U_XX-U_RR-2*i*U_XR)
+ // U_CB = d^2(U)/(dC*d\bar{C}) = 0.25*(U_XX+U_RR)
+ U_CC.r = 0.25 * (U.d11[ivar] - U.d22[ivar]);
+ U_CC.i = -0.5 * U.d12[ivar];
+ U_CB = 0.25 * (U.d11[ivar] + U.d22[ivar]);
+
+ // U_cc = C_cc*U_C+(C_c)^2*U_CC
+ U_cb = U_CB * C_c2;
+ U_cc = Cadd (Cmul (C_cc, U_C), Cmul (Cmul (C_c, C_c), U_CC));
+
+ // U_xx = 2*(U_cb+Re[U_cc])
+ // U_rr = 2*(U_cb-Re[U_cc])
+ // U_rx = -2*Im[U_cc]
+ U.d11[ivar] = 2 * (U_cb + U_cc.r);
+ U.d22[ivar] = 2 * (U_cb - U_cc.r);
+ U.d12[ivar] = -2 * U_cc.i;
+ }
+
+ *x = c.r;
+ *r = c.i;
+}
+
+//---------------------------------------------------------------
+void
+rx3_To_xyz (int nvar, double x, double r, double phi,
+ double *y, double *z, derivs U)
+// On Entrance: U.d0[]=U[]; U.d1[] =U[]_x; U.d2[] =U[]_r; U.d3[] =U[]_3;
+// U.d11[]=U[]_xx; U.d12[]=U[]_xr; U.d13[]=U[]_x3;
+// U.d22[]=U[]_rr; U.d23[]=U[]_r3; U.d33[]=U[]_33;
+// At Exit: U.d0[]=U[]; U.d1[] =U[]_x; U.d2[] =U[]_y; U.dz[] =U[]_z;
+// U.d11[]=U[]_xx; U.d12[]=U[]_xy; U.d1z[]=U[]_xz;
+// U.d22[]=U[]_yy; U.d2z[]=U[]_yz; U.dzz[]=U[]_zz;
+{
+ int jvar;
+ double
+ sin_phi = sin (phi),
+ cos_phi = cos (phi),
+ sin2_phi = sin_phi * sin_phi,
+ cos2_phi = cos_phi * cos_phi,
+ sin_2phi = 2 * sin_phi * cos_phi,
+ cos_2phi = cos2_phi - sin2_phi, r_inv = 1 / r, r_inv2 = r_inv * r_inv;
+
+ *y = r * cos_phi;
+ *z = r * sin_phi;
+
+ for (jvar = 0; jvar < nvar; jvar++)
+ {
+ double U_x = U.d1[jvar], U_r = U.d2[jvar], U_3 = U.d3[jvar],
+ U_xx = U.d11[jvar], U_xr = U.d12[jvar], U_x3 = U.d13[jvar],
+ U_rr = U.d22[jvar], U_r3 = U.d23[jvar], U_33 = U.d33[jvar];
+ U.d1[jvar] = U_x; // U_x
+ U.d2[jvar] = U_r * cos_phi - U_3 * r_inv * sin_phi; // U_y
+ U.d3[jvar] = U_r * sin_phi + U_3 * r_inv * cos_phi; // U_z
+ U.d11[jvar] = U_xx; // U_xx
+ U.d12[jvar] = U_xr * cos_phi - U_x3 * r_inv * sin_phi; // U_xy
+ U.d13[jvar] = U_xr * sin_phi + U_x3 * r_inv * cos_phi; // U_xz
+ U.d22[jvar] = U_rr * cos2_phi + r_inv2 * sin2_phi * (U_33 + r * U_r) // U_yy
+ + sin_2phi * r_inv2 * (U_3 - r * U_r3);
+ U.d23[jvar] = 0.5 * sin_2phi * (U_rr - r_inv * U_r - r_inv2 * U_33) // U_yz
+ - cos_2phi * r_inv2 * (U_3 - r * U_r3);
+ U.d33[jvar] = U_rr * sin2_phi + r_inv2 * cos2_phi * (U_33 + r * U_r) // U_zz
+ - sin_2phi * r_inv2 * (U_3 - r * U_r3);
+ }
+}
+
+//---------------------------------------------------------------
diff --git a/src/Equations.c b/src/Equations.c
new file mode 100644
index 0000000..b461d4e
--- /dev/null
+++ b/src/Equations.c
@@ -0,0 +1,186 @@
+// TwoPunctures: File "Equations.c"
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <ctype.h>
+#include "cctk_Parameters.h"
+#include "TP_utilities.h"
+#include "TwoPunctures.h"
+
+// U.d0[ivar] = U[ivar]; (ivar = 0..nvar-1)
+// U.d1[ivar] = U[ivar]_x;
+// U.d2[ivar] = U[ivar]_y;
+// U.d3[ivar] = U[ivar]_z;
+// U.d11[ivar] = U[ivar]_xx;
+// U.d12[ivar] = U[ivar]_xy;
+// U.d13[ivar] = U[ivar]_xz;
+// U.d22[ivar] = U[ivar]_yy;
+// U.d23[ivar] = U[ivar]_yz;
+// U.d33[ivar] = U[ivar]_zz;
+
+double
+BY_KKofxyz (double x, double y, double z)
+{
+ DECLARE_CCTK_PARAMETERS;
+ int i, j;
+ double r_plus, r2_plus, r3_plus, r_minus, r2_minus, r3_minus, np_Pp, nm_Pm,
+ Aij, AijAij, n_plus[3], n_minus[3], np_Sp[3], nm_Sm[3];
+
+ r2_plus = (x - par_b) * (x - par_b) + y * y + z * z;
+ r2_minus = (x + par_b) * (x + par_b) + y * y + z * z;
+ r_plus = sqrt (r2_plus);
+ r_minus = sqrt (r2_minus);
+ r3_plus = r_plus * r2_plus;
+ r3_minus = r_minus * r2_minus;
+
+ n_plus[0] = (x - par_b) / r_plus;
+ n_minus[0] = (x + par_b) / r_minus;
+ n_plus[1] = y / r_plus;
+ n_minus[1] = y / r_minus;
+ n_plus[2] = z / r_plus;
+ n_minus[2] = z / r_minus;
+
+ // dot product: np_Pp = (n_+).(P_+); nm_Pm = (n_-).(P_-)
+ np_Pp = 0;
+ nm_Pm = 0;
+ for (i = 0; i < 3; i++)
+ {
+ np_Pp += n_plus[i] * par_P_plus[i];
+ nm_Pm += n_minus[i] * par_P_minus[i];
+ }
+ // cross product: np_Sp[i] = [(n_+) x (S_+)]_i; nm_Sm[i] = [(n_-) x (S_-)]_i
+ np_Sp[0] = n_plus[1] * par_S_plus[2] - n_plus[2] * par_S_plus[1];
+ np_Sp[1] = n_plus[2] * par_S_plus[0] - n_plus[0] * par_S_plus[2];
+ np_Sp[2] = n_plus[0] * par_S_plus[1] - n_plus[1] * par_S_plus[0];
+ nm_Sm[0] = n_minus[1] * par_S_minus[2] - n_minus[2] * par_S_minus[1];
+ nm_Sm[1] = n_minus[2] * par_S_minus[0] - n_minus[0] * par_S_minus[2];
+ nm_Sm[2] = n_minus[0] * par_S_minus[1] - n_minus[1] * par_S_minus[0];
+ AijAij = 0;
+ for (i = 0; i < 3; i++)
+ {
+ for (j = 0; j < 3; j++)
+ { // Bowen-York-Curvature :
+ Aij =
+ +1.5 * (par_P_plus[i] * n_plus[j] + par_P_plus[j] * n_plus[i]
+ + np_Pp * n_plus[i] * n_plus[j]) / r2_plus
+ + 1.5 * (par_P_minus[i] * n_minus[j] + par_P_minus[j] * n_minus[i]
+ + nm_Pm * n_minus[i] * n_minus[j]) / r2_minus
+ - 3.0 * (np_Sp[i] * n_plus[j] + np_Sp[j] * n_plus[j]) / r3_plus
+ - 3.0 * (nm_Sm[i] * n_minus[j] + nm_Sm[i] * n_minus[j]) / r3_minus;
+ if (i == j)
+ Aij -= +1.5 * (np_Pp / r2_plus + nm_Pm / r2_minus);
+ AijAij += Aij * Aij;
+ }
+ }
+
+ return AijAij;
+}
+
+void
+BY_Aijofxyz (double x, double y, double z, double Aij[3][3])
+{
+ // Aij is a one-dimensional array with 9 elements
+ DECLARE_CCTK_PARAMETERS;
+ int i, j;
+ double r_plus, r2_plus, r3_plus, r_minus, r2_minus, r3_minus, np_Pp, nm_Pm,
+ n_plus[3], n_minus[3], np_Sp[3], nm_Sm[3];
+
+ r2_plus = (x - par_b) * (x - par_b) + y * y + z * z;
+ r2_minus = (x + par_b) * (x + par_b) + y * y + z * z;
+ r_plus = sqrt (r2_plus);
+ r_minus = sqrt (r2_minus);
+ r3_plus = r_plus * r2_plus;
+ r3_minus = r_minus * r2_minus;
+
+ n_plus[0] = (x - par_b) / r_plus;
+ n_minus[0] = (x + par_b) / r_minus;
+ n_plus[1] = y / r_plus;
+ n_minus[1] = y / r_minus;
+ n_plus[2] = z / r_plus;
+ n_minus[2] = z / r_minus;
+
+ // dot product: np_Pp = (n_+).(P_+); nm_Pm = (n_-).(P_-)
+ np_Pp = 0;
+ nm_Pm = 0;
+ for (i = 0; i < 3; i++)
+ {
+ np_Pp += n_plus[i] * par_P_plus[i];
+ nm_Pm += n_minus[i] * par_P_minus[i];
+ }
+ // cross product: np_Sp[i] = [(n_+) x (S_+)]_i; nm_Sm[i] = [(n_-) x (S_-)]_i
+ np_Sp[0] = n_plus[1] * par_S_plus[2] - n_plus[2] * par_S_plus[1];
+ np_Sp[1] = n_plus[2] * par_S_plus[0] - n_plus[0] * par_S_plus[2];
+ np_Sp[2] = n_plus[0] * par_S_plus[1] - n_plus[1] * par_S_plus[0];
+ nm_Sm[0] = n_minus[1] * par_S_minus[2] - n_minus[2] * par_S_minus[1];
+ nm_Sm[1] = n_minus[2] * par_S_minus[0] - n_minus[0] * par_S_minus[2];
+ nm_Sm[2] = n_minus[0] * par_S_minus[1] - n_minus[1] * par_S_minus[0];
+ for (i = 0; i < 3; i++)
+ {
+ for (j = 0; j < 3; j++)
+ { // Bowen-York-Curvature :
+ Aij[i][j] =
+ +1.5 * (par_P_plus[i] * n_plus[j] + par_P_plus[j] * n_plus[i]
+ + np_Pp * n_plus[i] * n_plus[j]) / r2_plus
+ + 1.5 * (par_P_minus[i] * n_minus[j] + par_P_minus[j] * n_minus[i]
+ + nm_Pm * n_minus[i] * n_minus[j]) / r2_minus
+ - 3.0 * (np_Sp[i] * n_plus[j] + np_Sp[j] * n_plus[j]) / r3_plus
+ - 3.0 * (nm_Sm[i] * n_minus[j] + nm_Sm[i] * n_minus[j]) / r3_minus;
+ if (i == j)
+ Aij[i][j] -= +1.5 * (np_Pp / r2_plus + nm_Pm / r2_minus);
+ }
+ }
+}
+
+//---------------------------------------------------------------
+//******* Nonlinear Equations **********
+//---------------------------------------------------------------
+void
+NonLinEquations (double A, double B, double X, double R,
+ double x, double r, double phi,
+ double y, double z, derivs U, double *values)
+{
+ DECLARE_CCTK_PARAMETERS;
+ double r_plus, r_minus, psi, psi2, psi4, psi7;
+ double mu;
+
+ r_plus = sqrt ((x - par_b) * (x - par_b) + y * y + z * z);
+ r_minus = sqrt ((x + par_b) * (x + par_b) + y * y + z * z);
+
+ psi =
+ 1. + 0.5 * par_m_plus / r_plus + 0.5 * par_m_minus / r_minus + U.d0[0];
+ psi2 = psi * psi;
+ psi4 = psi2 * psi2;
+ psi7 = psi * psi2 * psi4;
+
+ values[0] =
+ U.d11[0] + U.d22[0] + U.d33[0] + 0.125 * BY_KKofxyz (x, y, z) / psi7;
+
+}
+
+//---------------------------------------------------------------
+//******* Linear Equations **********
+//---------------------------------------------------------------
+void
+LinEquations (double A, double B, double X, double R,
+ double x, double r, double phi,
+ double y, double z, derivs dU, derivs U, double *values)
+{
+ DECLARE_CCTK_PARAMETERS;
+ double r_plus, r_minus, psi, psi2, psi4, psi8;
+
+ r_plus = sqrt ((x - par_b) * (x - par_b) + y * y + z * z);
+ r_minus = sqrt ((x + par_b) * (x + par_b) + y * y + z * z);
+
+ psi =
+ 1. + 0.5 * par_m_plus / r_plus + 0.5 * par_m_minus / r_minus + U.d0[0];
+ psi2 = psi * psi;
+ psi4 = psi2 * psi2;
+ psi8 = psi4 * psi4;
+
+ values[0] = dU.d11[0] + dU.d22[0] + dU.d33[0]
+ - 0.875 * BY_KKofxyz (x, y, z) / psi8 * dU.d0[0];
+}
+
+//---------------------------------------------------------------
diff --git a/src/FuncAndJacobian.c b/src/FuncAndJacobian.c
new file mode 100644
index 0000000..568ac13
--- /dev/null
+++ b/src/FuncAndJacobian.c
@@ -0,0 +1,744 @@
+// TwoPunctures: File "FuncAndJacobian.c"
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <ctype.h>
+#include <time.h>
+#include "cctk_Parameters.h"
+#include "TP_utilities.h"
+#include "TwoPunctures.h"
+
+#define FAC sin(al)*sin(be)*sin(al)*sin(be)*sin(al)*sin(be)
+//#define FAC sin(al)*sin(be)*sin(al)*sin(be)
+//#define FAC 1
+
+// -------------------------------------------------------------------------------
+int
+Index (int ivar, int i, int j, int k, int nvar, int n1, int n2, int n3)
+{
+ int i1 = i, j1 = j, k1 = k;
+
+ if (i1 < 0)
+ i1 = -(i1 + 1);
+ if (i1 >= n1)
+ i1 = 2 * n1 - (i1 + 1);
+
+ if (j1 < 0)
+ j1 = -(j1 + 1);
+ if (j1 >= n2)
+ j1 = 2 * n2 - (j1 + 1);
+
+ if (k1 < 0)
+ k1 = k1 + n3;
+ if (k1 >= n3)
+ k1 = k1 - n3;
+
+ return ivar + nvar * (i1 + n1 * (j1 + n2 * k1));
+}
+
+// -------------------------------------------------------------------------------
+void
+allocate_derivs (derivs * v, int n)
+{
+ int m = n - 1;
+ (*v).d0 = dvector (0, m);
+ (*v).d1 = dvector (0, m);
+ (*v).d2 = dvector (0, m);
+ (*v).d3 = dvector (0, m);
+ (*v).d11 = dvector (0, m);
+ (*v).d12 = dvector (0, m);
+ (*v).d13 = dvector (0, m);
+ (*v).d22 = dvector (0, m);
+ (*v).d23 = dvector (0, m);
+ (*v).d33 = dvector (0, m);
+}
+
+// -------------------------------------------------------------------------------
+void
+free_derivs (derivs * v, int n)
+{
+ int m = n - 1;
+ free_dvector ((*v).d0, 0, m);
+ free_dvector ((*v).d1, 0, m);
+ free_dvector ((*v).d2, 0, m);
+ free_dvector ((*v).d3, 0, m);
+ free_dvector ((*v).d11, 0, m);
+ free_dvector ((*v).d12, 0, m);
+ free_dvector ((*v).d13, 0, m);
+ free_dvector ((*v).d22, 0, m);
+ free_dvector ((*v).d23, 0, m);
+ free_dvector ((*v).d33, 0, m);
+}
+
+// -------------------------------------------------------------------------------
+void
+Derivatives_AB3 (int nvar, int n1, int n2, int n3, derivs v)
+{
+ int i, j, k, ivar, N, *indx;
+ double *p, *dp, *d2p, *q, *dq, *r, *dr;
+
+ N = maximum3 (n1, n2, n3);
+ p = dvector (0, N);
+ dp = dvector (0, N);
+ d2p = dvector (0, N);
+ q = dvector (0, N);
+ dq = dvector (0, N);
+ r = dvector (0, N);
+ dr = dvector (0, N);
+ indx = ivector (0, N);
+
+ for (ivar = 0; ivar < nvar; ivar++)
+ {
+ for (k = 0; k < n3; k++)
+ { // Calculation of Derivatives w.r.t. A-Dir.
+ for (j = 0; j < n2; j++)
+ { // (Chebyshev_Zeros)
+ for (i = 0; i < n1; i++)
+ {
+ indx[i] = Index (ivar, i, j, k, nvar, n1, n2, n3);
+ p[i] = v.d0[indx[i]];
+ }
+ chebft_Zeros (p, n1, 0);
+ chder (p, dp, n1);
+ chder (dp, d2p, n1);
+ chebft_Zeros (dp, n1, 1);
+ chebft_Zeros (d2p, n1, 1);
+ for (i = 0; i < n1; i++)
+ {
+ v.d1[indx[i]] = dp[i];
+ v.d11[indx[i]] = d2p[i];
+ }
+ }
+ }
+ for (k = 0; k < n3; k++)
+ { // Calculation of Derivatives w.r.t. A-Dir.
+ for (i = 0; i < n1; i++)
+ { // (Chebyshev_Zeros)
+ for (j = 0; j < n2; j++)
+ {
+ indx[j] = Index (ivar, i, j, k, nvar, n1, n2, n3);
+ p[j] = v.d0[indx[j]];
+ q[j] = v.d1[indx[j]];
+ }
+ chebft_Zeros (p, n2, 0);
+ chebft_Zeros (q, n2, 0);
+ chder (p, dp, n2);
+ chder (dp, d2p, n2);
+ chder (q, dq, n2);
+ chebft_Zeros (dp, n2, 1);
+ chebft_Zeros (d2p, n2, 1);
+ chebft_Zeros (dq, n2, 1);
+ for (j = 0; j < n2; j++)
+ {
+ v.d2[indx[j]] = dp[j];
+ v.d22[indx[j]] = d2p[j];
+ v.d12[indx[j]] = dq[j];
+ }
+ }
+ }
+ for (i = 0; i < n1; i++)
+ { // Calculation of Derivatives w.r.t. phi-Dir. (Fourier)
+ for (j = 0; j < n2; j++)
+ {
+ for (k = 0; k < n3; k++)
+ {
+ indx[k] = Index (ivar, i, j, k, nvar, n1, n2, n3);
+ p[k] = v.d0[indx[k]];
+ q[k] = v.d1[indx[k]];
+ r[k] = v.d2[indx[k]];
+ }
+ fourft (p, n3, 0);
+ fourder (p, dp, n3);
+ fourder2 (p, d2p, n3);
+ fourft (dp, n3, 1);
+ fourft (d2p, n3, 1);
+ fourft (q, n3, 0);
+ fourder (q, dq, n3);
+ fourft (dq, n3, 1);
+ fourft (r, n3, 0);
+ fourder (r, dr, n3);
+ fourft (dr, n3, 1);
+ for (k = 0; k < n3; k++)
+ {
+ v.d3[indx[k]] = dp[k];
+ v.d33[indx[k]] = d2p[k];
+ v.d13[indx[k]] = dq[k];
+ v.d23[indx[k]] = dr[k];
+ }
+ }
+ }
+ }
+ free_dvector (p, 0, N);
+ free_dvector (dp, 0, N);
+ free_dvector (d2p, 0, N);
+ free_dvector (q, 0, N);
+ free_dvector (dq, 0, N);
+ free_dvector (r, 0, N);
+ free_dvector (dr, 0, N);
+ free_ivector (indx, 0, N);
+}
+
+// -------------------------------------------------------------------------------
+void
+F_of_v (int nvar, int n1, int n2, int n3, derivs v, double *F,
+ derivs u)
+{ // Calculates the left hand sides of the non-linear equations F_m(v_n)=0
+ // and the function u (u.d0[]) as well as its derivatives
+ // (u.d1[], u.d2[], u.d3[], u.d11[], u.d12[], u.d13[], u.d22[], u.d23[], u.d33[])
+ // at interior points and at the boundaries "+/-"
+ int i, j, k, ivar, indx;
+ double al, be, A, B, X, R, x, r, phi, y, z, Am1, *values;
+ derivs U;
+
+ values = dvector (0, nvar - 1);
+ allocate_derivs (&U, nvar);
+
+ Derivatives_AB3 (nvar, n1, n2, n3, v);
+
+ for (i = 0; i < n1; i++)
+ {
+ for (j = 0; j < n2; j++)
+ {
+ for (k = 0; k < n3; k++)
+ {
+
+ al = Pih * (2 * i + 1) / n1;
+ A = -cos (al);
+ be = Pih * (2 * j + 1) / n2;
+ B = -cos (be);
+ phi = 2. * Pi * k / n3;
+
+ Am1 = A - 1;
+ for (ivar = 0; ivar < nvar; ivar++)
+ {
+ indx = Index (ivar, i, j, k, nvar, n1, n2, n3);
+ U.d0[ivar] = Am1 * v.d0[indx]; // U
+ U.d1[ivar] = v.d0[indx] + Am1 * v.d1[indx]; // U_A
+ U.d2[ivar] = Am1 * v.d2[indx]; // U_B
+ U.d3[ivar] = Am1 * v.d3[indx]; // U_3
+ U.d11[ivar] = 2 * v.d1[indx] + Am1 * v.d11[indx]; // U_AA
+ U.d12[ivar] = v.d2[indx] + Am1 * v.d12[indx]; // U_AB
+ U.d13[ivar] = v.d3[indx] + Am1 * v.d13[indx]; // U_AB
+ U.d22[ivar] = Am1 * v.d22[indx]; // U_BB
+ U.d23[ivar] = Am1 * v.d23[indx]; // U_B3
+ U.d33[ivar] = Am1 * v.d33[indx]; // U_33
+ }
+ // Calculation of (X,R) and
+ // (U_X, U_R, U_3, U_XX, U_XR, U_X3, U_RR, U_R3, U_33)
+ AB_To_XR (nvar, A, B, &X, &R, U);
+ // Calculation of (x,r) and
+ // (U, U_x, U_r, U_3, U_xx, U_xr, U_x3, U_rr, U_r3, U_33)
+ C_To_c (nvar, X, R, &x, &r, U);
+ // Calculation of (y,z) and
+ // (U, U_x, U_y, U_z, U_xx, U_xy, U_xz, U_yy, U_yz, U_zz)
+ rx3_To_xyz (nvar, x, r, phi, &y, &z, U);
+ NonLinEquations (A, B, X, R, x, r, phi, y, z, U, values);
+ for (ivar = 0; ivar < nvar; ivar++)
+ {
+ indx = Index (ivar, i, j, k, nvar, n1, n2, n3);
+ F[indx] = values[ivar] * FAC;
+ u.d0[indx] = U.d0[ivar]; // U
+ u.d1[indx] = U.d1[ivar]; // U_x
+ u.d2[indx] = U.d2[ivar]; // U_y
+ u.d3[indx] = U.d3[ivar]; // U_z
+ u.d11[indx] = U.d11[ivar]; // U_xx
+ u.d12[indx] = U.d12[ivar]; // U_xy
+ u.d13[indx] = U.d13[ivar]; // U_xz
+ u.d22[indx] = U.d22[ivar]; // U_yy
+ u.d23[indx] = U.d23[ivar]; // U_yz
+ u.d33[indx] = U.d33[ivar]; // U_zz
+ }
+ }
+ }
+ }
+ free_dvector (values, 0, nvar - 1);
+ free_derivs (&U, nvar);
+}
+
+// -------------------------------------------------------------------------------
+void
+J_times_dv (int nvar, int n1, int n2, int n3, derivs dv,
+ double *Jdv, derivs u)
+{ // Calculates the left hand sides of the non-linear equations F_m(v_n)=0
+ // and the function u (u.d0[]) as well as its derivatives
+ // (u.d1[], u.d2[], u.d3[], u.d11[], u.d12[], u.d13[], u.d22[], u.d23[], u.d33[])
+ // at interior points and at the boundaries "+/-"
+ DECLARE_CCTK_PARAMETERS;
+ int i, j, k, ivar, indx;
+ double al, be, A, B, X, R, x, r, phi, y, z, Am1, *values;
+ derivs dU, U;
+
+
+ values = dvector (0, nvar - 1);
+ allocate_derivs (&dU, nvar);
+ allocate_derivs (&U, nvar);
+
+ Derivatives_AB3 (nvar, n1, n2, n3, dv);
+
+ for (i = 0; i < n1; i++)
+ {
+ for (j = 0; j < n2; j++)
+ {
+ for (k = 0; k < n3; k++)
+ {
+
+ al = Pih * (2 * i + 1) / n1;
+ A = -cos (al);
+ be = Pih * (2 * j + 1) / n2;
+ B = -cos (be);
+ phi = 2. * Pi * k / n3;
+
+ Am1 = A - 1;
+ for (ivar = 0; ivar < nvar; ivar++)
+ {
+ indx = Index (ivar, i, j, k, nvar, n1, n2, n3);
+ dU.d0[ivar] = Am1 * dv.d0[indx]; // dU
+ dU.d1[ivar] = dv.d0[indx] + Am1 * dv.d1[indx]; // dU_A
+ dU.d2[ivar] = Am1 * dv.d2[indx]; // dU_B
+ dU.d3[ivar] = Am1 * dv.d3[indx]; // dU_3
+ dU.d11[ivar] = 2 * dv.d1[indx] + Am1 * dv.d11[indx]; // dU_AA
+ dU.d12[ivar] = dv.d2[indx] + Am1 * dv.d12[indx]; // dU_AB
+ dU.d13[ivar] = dv.d3[indx] + Am1 * dv.d13[indx]; // dU_AB
+ dU.d22[ivar] = Am1 * dv.d22[indx]; // dU_BB
+ dU.d23[ivar] = Am1 * dv.d23[indx]; // dU_B3
+ dU.d33[ivar] = Am1 * dv.d33[indx]; // dU_33
+ U.d0[ivar] = u.d0[indx]; // U
+ U.d1[ivar] = u.d1[indx]; // U_x
+ U.d2[ivar] = u.d2[indx]; // U_y
+ U.d3[ivar] = u.d3[indx]; // U_z
+ U.d11[ivar] = u.d11[indx]; // U_xx
+ U.d12[ivar] = u.d12[indx]; // U_xy
+ U.d13[ivar] = u.d13[indx]; // U_xz
+ U.d22[ivar] = u.d22[indx]; // U_yy
+ U.d23[ivar] = u.d23[indx]; // U_yz
+ U.d33[ivar] = u.d33[indx]; // U_zz
+ }
+ // Calculation of (X,R) and
+ // (dU_X, dU_R, dU_3, dU_XX, dU_XR, dU_X3, dU_RR, dU_R3, dU_33)
+ AB_To_XR (nvar, A, B, &X, &R, dU);
+ // Calculation of (x,r) and
+ // (dU, dU_x, dU_r, dU_3, dU_xx, dU_xr, dU_x3, dU_rr, dU_r3, dU_33)
+ C_To_c (nvar, X, R, &x, &r, dU);
+ // Calculation of (y,z) and
+ // (dU, dU_x, dU_y, dU_z, dU_xx, dU_xy, dU_xz, dU_yy, dU_yz, dU_zz)
+ rx3_To_xyz (nvar, x, r, phi, &y, &z, dU);
+ LinEquations (A, B, X, R, x, r, phi, y, z, dU, U, values);
+ for (ivar = 0; ivar < nvar; ivar++)
+ {
+ indx = Index (ivar, i, j, k, nvar, n1, n2, n3);
+ Jdv[indx] = values[ivar] * FAC;
+ }
+ }
+ }
+ }
+ free_dvector (values, 0, nvar - 1);
+ free_derivs (&dU, nvar);
+ free_derivs (&U, nvar);
+}
+
+// -------------------------------------------------------------------------------
+void
+JFD_times_dv (int i, int j, int k, int nvar, int n1, int n2,
+ int n3, derivs dv, derivs u, double *values)
+{ // Calculates rows of the vector 'J(FD)*dv'.
+ // First row to be calculated: row = Index(0, i, j, k; nvar, n1, n2, n3)
+ // Last row to be calculated: row = Index(nvar-1, i, j, k; nvar, n1, n2, n3)
+ // These rows are stored in the vector JFDdv[0] ... JFDdv[nvar-1].
+ DECLARE_CCTK_PARAMETERS;
+ int ivar, indx;
+ double al, be, A, B, X, R, x, r, phi, y, z, Am1;
+ double sin_al, sin_al_i1, sin_al_i2, sin_al_i3, cos_al;
+ double sin_be, sin_be_i1, sin_be_i2, sin_be_i3, cos_be;
+ double dV0, dV1, dV2, dV3, dV11, dV12, dV13, dV22, dV23, dV33,
+ ha, ga, ga2, hb, gb, gb2, hp, gp, gp2, gagb, gagp, gbgp;
+ derivs dU, U;
+
+ allocate_derivs (&dU, nvar);
+ allocate_derivs (&U, nvar);
+
+ if (k < 0)
+ k = k + n3;
+ if (k >= n3)
+ k = k - n3;
+
+ ha = Pi / n1; // ha: Stepsize with respect to (al)
+ al = ha * (i + 0.5);
+ A = -cos (al);
+ ga = 1 / ha;
+ ga2 = ga * ga;
+
+ hb = Pi / n2; // hb: Stepsize with respect to (be)
+ be = hb * (j + 0.5);
+ B = -cos (be);
+ gb = 1 / hb;
+ gb2 = gb * gb;
+ gagb = ga * gb;
+
+ hp = 2 * Pi / n3; // hp: Stepsize with respect to (phi)
+ phi = hp * j;
+ gp = 1 / hp;
+ gp2 = gp * gp;
+ gagp = ga * gp;
+ gbgp = gb * gp;
+
+
+ sin_al = sin (al);
+ sin_be = sin (be);
+ sin_al_i1 = 1 / sin_al;
+ sin_be_i1 = 1 / sin_be;
+ sin_al_i2 = sin_al_i1 * sin_al_i1;
+ sin_be_i2 = sin_be_i1 * sin_be_i1;
+ sin_al_i3 = sin_al_i1 * sin_al_i2;
+ sin_be_i3 = sin_be_i1 * sin_be_i2;
+ cos_al = -A;
+ cos_be = -B;
+
+ Am1 = A - 1;
+ for (ivar = 0; ivar < nvar; ivar++)
+ {
+ int iccc = Index (ivar, i, j, k, nvar, n1, n2, n3),
+ ipcc = Index (ivar, i + 1, j, k, nvar, n1, n2, n3),
+ imcc = Index (ivar, i - 1, j, k, nvar, n1, n2, n3),
+ icpc = Index (ivar, i, j + 1, k, nvar, n1, n2, n3),
+ icmc = Index (ivar, i, j - 1, k, nvar, n1, n2, n3),
+ iccp = Index (ivar, i, j, k + 1, nvar, n1, n2, n3),
+ iccm = Index (ivar, i, j, k - 1, nvar, n1, n2, n3),
+ icpp = Index (ivar, i, j + 1, k + 1, nvar, n1, n2, n3),
+ icmp = Index (ivar, i, j - 1, k + 1, nvar, n1, n2, n3),
+ icpm = Index (ivar, i, j + 1, k - 1, nvar, n1, n2, n3),
+ icmm = Index (ivar, i, j - 1, k - 1, nvar, n1, n2, n3),
+ ipcp = Index (ivar, i + 1, j, k + 1, nvar, n1, n2, n3),
+ imcp = Index (ivar, i - 1, j, k + 1, nvar, n1, n2, n3),
+ ipcm = Index (ivar, i + 1, j, k - 1, nvar, n1, n2, n3),
+ imcm = Index (ivar, i - 1, j, k - 1, nvar, n1, n2, n3),
+ ippc = Index (ivar, i + 1, j + 1, k, nvar, n1, n2, n3),
+ impc = Index (ivar, i - 1, j + 1, k, nvar, n1, n2, n3),
+ ipmc = Index (ivar, i + 1, j - 1, k, nvar, n1, n2, n3),
+ immc = Index (ivar, i - 1, j - 1, k, nvar, n1, n2, n3);
+ // Derivatives of (dv) w.r.t. (al,be,phi):
+ dV0 = dv.d0[iccc];
+ dV1 = 0.5 * ga * (dv.d0[ipcc] - dv.d0[imcc]);
+ dV2 = 0.5 * gb * (dv.d0[icpc] - dv.d0[icmc]);
+ dV3 = 0.5 * gp * (dv.d0[iccp] - dv.d0[iccm]);
+ dV11 = ga2 * (dv.d0[ipcc] + dv.d0[imcc] - 2 * dv.d0[iccc]);
+ dV22 = gb2 * (dv.d0[icpc] + dv.d0[icmc] - 2 * dv.d0[iccc]);
+ dV33 = gp2 * (dv.d0[iccp] + dv.d0[iccm] - 2 * dv.d0[iccc]);
+ dV12 =
+ 0.25 * gagb * (dv.d0[ippc] - dv.d0[ipmc] + dv.d0[immc] - dv.d0[impc]);
+ dV13 =
+ 0.25 * gagp * (dv.d0[ipcp] - dv.d0[imcp] + dv.d0[imcm] - dv.d0[ipcm]);
+ dV23 =
+ 0.25 * gbgp * (dv.d0[icpp] - dv.d0[icpm] + dv.d0[icmm] - dv.d0[icmp]);
+ // Derivatives of (dv) w.r.t. (A,B,phi):
+ dV11 = sin_al_i3 * (sin_al * dV11 - cos_al * dV1);
+ dV12 = sin_al_i1 * sin_be_i1 * dV12;
+ dV13 = sin_al_i1 * dV13;
+ dV22 = sin_be_i3 * (sin_be * dV22 - cos_be * dV2);
+ dV23 = sin_be_i1 * dV23;
+ dV1 = sin_al_i1 * dV1;
+ dV2 = sin_be_i1 * dV2;
+ // Derivatives of (dU) w.r.t. (A,B,phi):
+ dU.d0[ivar] = Am1 * dV0;
+ dU.d1[ivar] = dV0 + Am1 * dV1;
+ dU.d2[ivar] = Am1 * dV2;
+ dU.d3[ivar] = Am1 * dV3;
+ dU.d11[ivar] = 2 * dV1 + Am1 * dV11;
+ dU.d12[ivar] = dV2 + Am1 * dV12;
+ dU.d13[ivar] = dV3 + Am1 * dV13;
+ dU.d22[ivar] = Am1 * dV22;
+ dU.d23[ivar] = Am1 * dV23;
+ dU.d33[ivar] = Am1 * dV33;
+
+ indx = Index (ivar, i, j, k, nvar, n1, n2, n3);
+ U.d0[ivar] = u.d0[indx]; // U
+ U.d1[ivar] = u.d1[indx]; // U_x
+ U.d2[ivar] = u.d2[indx]; // U_y
+ U.d3[ivar] = u.d3[indx]; // U_z
+ U.d11[ivar] = u.d11[indx]; // U_xx
+ U.d12[ivar] = u.d12[indx]; // U_xy
+ U.d13[ivar] = u.d13[indx]; // U_xz
+ U.d22[ivar] = u.d22[indx]; // U_yy
+ U.d23[ivar] = u.d23[indx]; // U_yz
+ U.d33[ivar] = u.d33[indx]; // U_zz
+ }
+ // Calculation of (X,R) and
+ // (dU_X, dU_R, dU_3, dU_XX, dU_XR, dU_X3, dU_RR, dU_R3, dU_33)
+ AB_To_XR (nvar, A, B, &X, &R, dU);
+ // Calculation of (x,r) and
+ // (dU, dU_x, dU_r, dU_3, dU_xx, dU_xr, dU_x3, dU_rr, dU_r3, dU_33)
+ C_To_c (nvar, X, R, &x, &r, dU);
+ // Calculation of (y,z) and
+ // (dU, dU_x, dU_y, dU_z, dU_xx, dU_xy, dU_xz, dU_yy, dU_yz, dU_zz)
+ rx3_To_xyz (nvar, x, r, phi, &y, &z, dU);
+ LinEquations (A, B, X, R, x, r, phi, y, z, dU, U, values);
+ for (ivar = 0; ivar < nvar; ivar++)
+ values[ivar] *= FAC;
+
+ free_derivs (&dU, nvar);
+ free_derivs (&U, nvar);
+}
+
+// -------------------------------------------------------------------------------
+void
+SetMatrix_JFD (int nvar, int n1, int n2, int n3, derivs u,
+ int *ncols, int **cols, double **Matrix)
+{
+ DECLARE_CCTK_PARAMETERS;
+ int column, row, mcol;
+ int i, i1, i_0, i_1, j, j1, j_0, j_1, k, k1, k_0, k_1, N1, N2, N3,
+ ivar, ivar1, ntotal = nvar * n1 * n2 * n3;
+ double *values;
+ derivs dv;
+
+ values = dvector (0, nvar - 1);
+ allocate_derivs (&dv, ntotal);
+
+ N1 = n1 - 1;
+ N2 = n2 - 1;
+ N3 = n3 - 1;
+
+ for (i = 0; i < n1; i++)
+ {
+ for (j = 0; j < n2; j++)
+ {
+ for (k = 0; k < n3; k++)
+ {
+ for (ivar = 0; ivar < nvar; ivar++)
+ {
+ row = Index (ivar, i, j, k, nvar, n1, n2, n3);
+ ncols[row] = 0;
+ dv.d0[row] = 0;
+ }
+ }
+ }
+ }
+ for (i = 0; i < n1; i++)
+ {
+ for (j = 0; j < n2; j++)
+ {
+ for (k = 0; k < n3; k++)
+ {
+ for (ivar = 0; ivar < nvar; ivar++)
+ {
+ column = Index (ivar, i, j, k, nvar, n1, n2, n3);
+ dv.d0[column] = 1;
+
+ i_0 = maximum2 (0, i - 1);
+ i_1 = minimum2 (N1, i + 1);
+ j_0 = maximum2 (0, j - 1);
+ j_1 = minimum2 (N2, j + 1);
+ k_0 = k - 1;
+ k_1 = k + 1;
+/* i_0 = 0;
+ i_1 = N1;
+ j_0 = 0;
+ j_1 = N2;
+ k_0 = 0;
+ k_1 = N3;*/
+
+ for (i1 = i_0; i1 <= i_1; i1++)
+ {
+ for (j1 = j_0; j1 <= j_1; j1++)
+ {
+ for (k1 = k_0; k1 <= k_1; k1++)
+ {
+ JFD_times_dv (i1, j1, k1, nvar, n1, n2, n3,
+ dv, u, values);
+ for (ivar1 = 0; ivar1 < nvar; ivar1++)
+ {
+ if (values[ivar1] != 0)
+ {
+ row = Index (ivar1, i1, j1, k1, nvar, n1, n2, n3);
+ mcol = ncols[row];
+ cols[row][mcol] = column;
+ Matrix[row][mcol] = values[ivar1];
+ ncols[row] += 1;
+ }
+ }
+ }
+ }
+ }
+
+ dv.d0[column] = 0;
+ }
+ }
+ }
+ }
+ free_derivs (&dv, ntotal);
+ free_dvector (values, 0, nvar - 1);
+}
+
+// -------------------------------------------------------------------------------
+// Calculates the value of v at an arbitrary position (A,B,phi)
+double
+PunctEvalAtArbitPosition (double *v, double A, double B, double phi,
+ int n1, int n2, int n3)
+{
+ int i, j, k, N;
+ double *p, *values1, **values2, result;
+
+ N = maximum3 (n1, n2, n3);
+ p = dvector (0, N);
+ values1 = dvector (0, N);
+ values2 = dmatrix (0, N, 0, N);
+
+ for (k = 0; k < n3; k++)
+ {
+ for (j = 0; j < n2; j++)
+ {
+ for (i = 0; i < n1; i++)
+ p[i] = v[i + n1 * (j + n2 * k)];
+ chebft_Zeros (p, n1, 0);
+ values2[j][k] = chebev (-1, 1, p, n1, A);
+ }
+ }
+
+ for (k = 0; k < n3; k++)
+ {
+ for (j = 0; j < n2; j++)
+ p[j] = values2[j][k];
+ chebft_Zeros (p, n2, 0);
+ values1[k] = chebev (-1, 1, p, n2, B);
+ }
+
+ fourft (values1, n3, 0);
+ result = fourev (values1, n3, phi);
+
+ free_dvector (p, 0, N);
+ free_dvector (values1, 0, N);
+ free_dmatrix (values2, 0, N, 0, N);
+
+ return result;
+}
+
+// -------------------------------------------------------------------------------
+void
+calculate_derivs (int i, int j, int k, int ivar, int nvar, int n1, int n2,
+ int n3, derivs v, derivs vv)
+{
+ double al = Pih * (2 * i + 1) / n1, be = Pih * (2 * j + 1) / n2,
+ sin_al = sin (al), sin2_al = sin_al * sin_al, cos_al = cos (al),
+ sin_be = sin (be), sin2_be = sin_be * sin_be, cos_be = cos (be);
+
+ vv.d0[0] = v.d0[Index (ivar, i, j, k, nvar, n1, n2, n3)];
+ vv.d1[0] = v.d1[Index (ivar, i, j, k, nvar, n1, n2, n3)] * sin_al;
+ vv.d2[0] = v.d2[Index (ivar, i, j, k, nvar, n1, n2, n3)] * sin_be;
+ vv.d3[0] = v.d3[Index (ivar, i, j, k, nvar, n1, n2, n3)];
+ vv.d11[0] = v.d11[Index (ivar, i, j, k, nvar, n1, n2, n3)] * sin2_al
+ + v.d1[Index (ivar, i, j, k, nvar, n1, n2, n3)] * cos_al;
+ vv.d12[0] =
+ v.d12[Index (ivar, i, j, k, nvar, n1, n2, n3)] * sin_al * sin_be;
+ vv.d13[0] = v.d13[Index (ivar, i, j, k, nvar, n1, n2, n3)] * sin_al;
+ vv.d22[0] = v.d22[Index (ivar, i, j, k, nvar, n1, n2, n3)] * sin2_be
+ + v.d2[Index (ivar, i, j, k, nvar, n1, n2, n3)] * cos_be;
+ vv.d23[0] = v.d23[Index (ivar, i, j, k, nvar, n1, n2, n3)] * sin_be;
+ vv.d33[0] = v.d33[Index (ivar, i, j, k, nvar, n1, n2, n3)];
+}
+
+// -------------------------------------------------------------------------------
+double
+interpol (double a, double b, double c, derivs v)
+{
+ return v.d0[0]
+ + a * v.d1[0] + b * v.d2[0] + c * v.d3[0]
+ + 0.5 * a * a * v.d11[0] + a * b * v.d12[0] + a * c * v.d13[0]
+ + 0.5 * b * b * v.d22[0] + b * c * v.d23[0] + 0.5 * c * c * v.d33[0];
+}
+
+// -------------------------------------------------------------------------------
+// Calculates the value of v at an arbitrary position (A,B,phi)
+double
+PunctIntPolAtArbitPosition (int ivar, int nvar, int n1,
+ int n2, int n3, derivs v, double x, double y,
+ double z)
+{
+ DECLARE_CCTK_PARAMETERS;
+ double xs, ys, zs, rs2, phi, X, R, A, B, al, be, aux1, aux2, a, b, c,
+ result, Us, Ui;
+ int i, j, k;
+ derivs vv;
+ allocate_derivs (&vv, 1);
+
+ xs = x / par_b;
+ ys = y / par_b;
+ zs = z / par_b;
+ rs2 = ys * ys + zs * zs;
+ phi = atan2 (z, y);
+ if (phi < 0)
+ phi += 2 * Pi;
+
+ aux1 = 0.5 * (xs * xs + rs2 - 1);
+ aux2 = sqrt (aux1 * aux1 + rs2);
+ X = asinh (sqrt (aux1 + aux2));
+ R = asin (sqrt (-aux1 + aux2));
+ if (x < 0)
+ R = Pi - R;
+
+ A = 2 * tanh (0.5 * X) - 1;
+ B = tan (0.5 * R - Piq);
+ al = Pi - acos (A);
+ be = Pi - acos (B);
+
+ i = rint (al * n1 / Pi - 0.5);
+ j = rint (be * n2 / Pi - 0.5);
+ k = rint (0.5 * phi * n3 / Pi);
+
+ a = al - Pi * (i + 0.5) / n1;
+ b = be - Pi * (j + 0.5) / n2;
+ c = phi - 2 * Pi * k / n3;
+
+ calculate_derivs (i, j, k, ivar, nvar, n1, n2, n3, v, vv);
+ result = interpol (a, b, c, vv);
+ free_derivs (&vv, 1);
+
+// Us = (A-1)*PunctEvalAtArbitPosition(v.d0, A, B, phi, n1, n2, n3);
+ Ui = (A - 1) * result;
+// printf("%e %e %e Us=%e Ui=%e 1-Ui/Us=%e\n",x,y,z,Us,Ui,1-Ui/Us);
+
+ return Ui;
+
+/* calculate_derivs(i+1, j, k, ivar, nvar, n1, n2, n3, v, vv, 2);
+ calculate_derivs(i, j+1, k, ivar, nvar, n1, n2, n3, v, vv, 3);
+ calculate_derivs(i+1, j+1, k, ivar, nvar, n1, n2, n3, v, vv, 4);
+ calculate_derivs(i, j, k+1, ivar, nvar, n1, n2, n3, v, vv, 5);
+ calculate_derivs(i+1, j, k+1, ivar, nvar, n1, n2, n3, v, vv, 6);
+ calculate_derivs(i, j+1, k+1, ivar, nvar, n1, n2, n3, v, vv, 7);
+ calculate_derivs(i+1, j+1, k+1, ivar, nvar, n1, n2, n3, v, vv, 8);
+ v_intpol[2]=interpol(ha*(a-1),hb*b, hp*c, vv,2);
+ v_intpol[3]=interpol(ha*a, hb*(b-1),hp*c, vv,3);
+ v_intpol[4]=interpol(ha*(a-1),hb*(b-1),hp*c, vv,4);
+ v_intpol[5]=interpol(ha*a, hb*b, hp*(c-1),vv,5);
+ v_intpol[6]=interpol(ha*(a-1),hb*b, hp*(c-1),vv,6);
+ v_intpol[7]=interpol(ha*a, hb*(b-1),hp*(c-1),vv,7);
+ v_intpol[8]=interpol(ha*(a-1),hb*(b-1),hp*(c-1),vv,8);
+ for(j=1;j<=8;j++)
+ printf("Value U[V] at origin:%16.15f\tj=%d\n",-2*v_intpol[j],j);*/
+/* result = 0;
+ for(mi=i-npol; mi<=i+npol; mi++){
+ al_m = Pih*(2*mi+1)/n1;
+ for(mj=j-npol; mj<=j+npol; mj++){
+ be_m = Pih*(2*mj+1)/n2;
+ for(mk=k-npol; mk<=k+npol; mk++){
+ phi_m = 2.*Pi*mk/n3;
+ g_m = 1;
+ for(ni=i-npol; ni<=i+npol; ni++){
+ al_n = Pih*(2*ni+1)/n1;
+ if(ni != mi) g_m *= (al - al_n)/(al_m - al_n);
+ }
+ for(nj=j-npol; nj<=j+npol; nj++){
+ be_n = Pih*(2*nj+1)/n2;
+ if(nj != mj) g_m *= (be - be_n)/(be_m - be_n);
+ }
+ for(nk=k-npol; nk<=k+npol; nk++){
+ phi_n = 2.*Pi*nk/n3;
+ if(nk != mk) g_m *= (phi-phi_n)/(phi_m - phi_n);
+ }
+ result += g_m*v.d0[Index(ivar,mi,mj,mk,nvar,n1,n2,n3)];
+ }
+ }
+ } */
+}
+
+// -------------------------------------------------------------------------------
diff --git a/src/Newton.c b/src/Newton.c
new file mode 100644
index 0000000..42ceaa3
--- /dev/null
+++ b/src/Newton.c
@@ -0,0 +1,435 @@
+// TwoPunctures: File "Newton.c"
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <ctype.h>
+#include "cctk_Parameters.h"
+#include "TP_utilities.h"
+#include "TwoPunctures.h"
+
+// -----------------------------------------------------------------------------------
+void
+resid (double *res, int ntotal, double *dv, double *rhs,
+ int *ncols, int **cols, double **JFD)
+{
+ int i, m;
+ double JFDdv_i;
+
+ for (i = 0; i < ntotal; i++)
+ {
+ JFDdv_i = 0;
+ for (m = 0; m < ncols[i]; m++)
+ JFDdv_i += JFD[i][m] * dv[cols[i][m]];
+ res[i] = rhs[i] - JFDdv_i;
+ }
+}
+
+// -----------------------------------------------------------------------------------
+void
+LineRelax_al (double *dv, int j, int k, int nvar, int n1, int n2, int n3,
+ double *rhs, int *ncols, int **cols, double **JFD)
+{
+ int i, m, Ic, Ip, Im, col, ivar;
+ double *a, *b, *c, *r, *u;
+ a = dvector (0, n1 - 1);
+ b = dvector (0, n1 - 1);
+ c = dvector (0, n1 - 1);
+ r = dvector (0, n1 - 1);
+ u = dvector (0, n1 - 1);
+
+ for (ivar = 0; ivar < nvar; ivar++)
+ {
+ for (i = 0; i < n1; i++)
+ {
+ Ip = Index (ivar, i + 1, j, k, nvar, n1, n2, n3);
+ Ic = Index (ivar, i, j, k, nvar, n1, n2, n3);
+ Im = Index (ivar, i - 1, j, k, nvar, n1, n2, n3);
+ a[i] = 0;
+ b[i] = 0;
+ c[i] = 0;
+ r[i] = rhs[Ic];
+ for (m = 0; m < ncols[Ic]; m++)
+ {
+ col = cols[Ic][m];
+ if (col != Ip && col != Ic && col != Im)
+ r[i] -= JFD[Ic][m] * dv[col];
+ else
+ {
+ if (col == Im)
+ a[i] = JFD[Ic][m];
+ if (col == Ic)
+ b[i] = JFD[Ic][m];
+ if (col == Ip)
+ c[i] = JFD[Ic][m];
+ }
+ }
+ }
+ tridag (a, b, c, r, u, n1);
+ for (i = 0; i < n1; i++)
+ {
+ Ic = Index (ivar, i, j, k, nvar, n1, n2, n3);
+ dv[Ic] = u[i];
+ }
+ }
+ free_dvector (a, 0, n1 - 1);
+ free_dvector (b, 0, n1 - 1);
+ free_dvector (c, 0, n1 - 1);
+ free_dvector (r, 0, n1 - 1);
+ free_dvector (u, 0, n1 - 1);
+}
+
+// -----------------------------------------------------------------------------------
+void
+LineRelax_be (double *dv, int i, int k, int nvar, int n1, int n2, int n3,
+ double *rhs, int *ncols, int **cols, double **JFD)
+{
+ int j, m, Ic, Ip, Im, col, ivar;
+ double *a, *b, *c, *r, *u;
+ a = dvector (0, n2 - 1);
+ b = dvector (0, n2 - 1);
+ c = dvector (0, n2 - 1);
+ r = dvector (0, n2 - 1);
+ u = dvector (0, n2 - 1);
+
+ for (ivar = 0; ivar < nvar; ivar++)
+ {
+ for (j = 0; j < n2; j++)
+ {
+ Ip = Index (ivar, i, j + 1, k, nvar, n1, n2, n3);
+ Ic = Index (ivar, i, j, k, nvar, n1, n2, n3);
+ Im = Index (ivar, i, j - 1, k, nvar, n1, n2, n3);
+ a[j] = 0;
+ b[j] = 0;
+ c[j] = 0;
+ r[j] = rhs[Ic];
+ for (m = 0; m < ncols[Ic]; m++)
+ {
+ col = cols[Ic][m];
+ if (col != Ip && col != Ic && col != Im)
+ r[j] -= JFD[Ic][m] * dv[col];
+ else
+ {
+ if (col == Im)
+ a[j] = JFD[Ic][m];
+ if (col == Ic)
+ b[j] = JFD[Ic][m];
+ if (col == Ip)
+ c[j] = JFD[Ic][m];
+ }
+ }
+ }
+ tridag (a, b, c, r, u, n2);
+ for (j = 0; j < n2; j++)
+ {
+ Ic = Index (ivar, i, j, k, nvar, n1, n2, n3);
+ dv[Ic] = u[j];
+ }
+ }
+ free_dvector (a, 0, n2 - 1);
+ free_dvector (b, 0, n2 - 1);
+ free_dvector (c, 0, n2 - 1);
+ free_dvector (r, 0, n2 - 1);
+ free_dvector (u, 0, n2 - 1);
+}
+
+// -----------------------------------------------------------------------------------
+static void
+relax (double *dv, int nvar, int n1, int n2, int n3,
+ double *rhs, int *ncols, int **cols, double **JFD)
+{
+ int i, j, k, n;
+
+ for (k = 0; k < n3; k = k + 2)
+ {
+ for (n = 0; n < N_PlaneRelax; n++)
+ {
+ for (i = 2; i < n1; i = i + 2)
+ LineRelax_be (dv, i, k, nvar, n1, n2, n3, rhs, ncols, cols, JFD);
+ for (i = 1; i < n1; i = i + 2)
+ LineRelax_be (dv, i, k, nvar, n1, n2, n3, rhs, ncols, cols, JFD);
+ for (j = 1; j < n2; j = j + 2)
+ LineRelax_al (dv, j, k, nvar, n1, n2, n3, rhs, ncols, cols, JFD);
+ for (j = 0; j < n2; j = j + 2)
+ LineRelax_al (dv, j, k, nvar, n1, n2, n3, rhs, ncols, cols, JFD);
+ }
+ }
+ for (k = 1; k < n3; k = k + 2)
+ {
+ for (n = 0; n < N_PlaneRelax; n++)
+ {
+ for (i = 0; i < n1; i = i + 2)
+ LineRelax_be (dv, i, k, nvar, n1, n2, n3, rhs, ncols, cols, JFD);
+ for (i = 1; i < n1; i = i + 2)
+ LineRelax_be (dv, i, k, nvar, n1, n2, n3, rhs, ncols, cols, JFD);
+ for (j = 1; j < n2; j = j + 2)
+ LineRelax_al (dv, j, k, nvar, n1, n2, n3, rhs, ncols, cols, JFD);
+ for (j = 0; j < n2; j = j + 2)
+ LineRelax_al (dv, j, k, nvar, n1, n2, n3, rhs, ncols, cols, JFD);
+ }
+ }
+}
+
+// -----------------------------------------------------------------------------------
+void
+TestRelax (int nvar, int n1, int n2, int n3, derivs v,
+ double *dv)
+{
+ DECLARE_CCTK_PARAMETERS;
+ int ntotal = n1 * n2 * n3 * nvar, **cols, *ncols, maxcol =
+ StencilSize * nvar, j;
+ double *F, *res, **JFD;
+ derivs u;
+
+ F = dvector (0, ntotal - 1);
+ res = dvector (0, ntotal - 1);
+ allocate_derivs (&u, ntotal);
+
+ JFD = dmatrix (0, ntotal - 1, 0, maxcol - 1);
+ cols = imatrix (0, ntotal - 1, 0, maxcol - 1);
+ ncols = ivector (0, ntotal - 1);
+
+ F_of_v (nvar, n1, n2, n3, v, F, u);
+
+ SetMatrix_JFD (nvar, n1, n2, n3, u, ncols, cols, JFD);
+
+ for (j = 0; j < ntotal; j++)
+ dv[j] = 0;
+ resid (res, ntotal, dv, F, ncols, cols, JFD);
+ printf ("Davor: |F|=%20.15e\n", norm1 (res, ntotal));
+ for (j = 0; j < NRELAX; j++)
+ {
+ relax (dv, nvar, n1, n2, n3, F, ncols, cols, JFD); // solves JFD*sh = s
+ if (j % Step_Relax == 0)
+ {
+ resid (res, ntotal, dv, F, ncols, cols, JFD);
+ printf ("j=%d\t |F|=%20.15e\n", j, norm1 (res, ntotal));
+ }
+ }
+
+ resid (res, ntotal, dv, F, ncols, cols, JFD);
+ printf ("Danach: |F|=%20.15e\n", norm1 (res, ntotal));
+
+ free_dvector (F, 0, ntotal - 1);
+ free_dvector (res, 0, ntotal - 1);
+ free_derivs (&u, ntotal);
+
+ free_dmatrix (JFD, 0, ntotal - 1, 0, maxcol - 1);
+ free_imatrix (cols, 0, ntotal - 1, 0, maxcol - 1);
+ free_ivector (ncols, 0, ntotal - 1);
+}
+
+// -----------------------------------------------------------------------------------
+int
+bicgstab (int nvar, int n1, int n2, int n3, derivs v,
+ derivs dv, int output, int itmax, double tol, double *normres)
+{
+ DECLARE_CCTK_PARAMETERS;
+ int ntotal = n1 * n2 * n3 * nvar, ii, j;
+ double alpha = 0, beta = 0;
+ double rho = 0, rho1 = 1, rhotol = 1e-50;
+ double omega = 0, omegatol = 1e-50;
+ double *p, *rt, *s, *t, *r, *vv;
+ double **JFD;
+ int **cols, *ncols, maxcol = StencilSize * nvar;
+ double *F;
+ derivs u, ph, sh;
+
+ F = dvector (0, ntotal - 1);
+ allocate_derivs (&u, ntotal);
+
+ JFD = dmatrix (0, ntotal - 1, 0, maxcol - 1);
+ cols = imatrix (0, ntotal - 1, 0, maxcol - 1);
+ ncols = ivector (0, ntotal - 1);
+
+ F_of_v (nvar, n1, n2, n3, v, F, u);
+ SetMatrix_JFD (nvar, n1, n2, n3, u, ncols, cols, JFD);
+
+ /* temporary storage */
+ r = dvector (0, ntotal - 1);
+ p = dvector (0, ntotal - 1);
+ allocate_derivs (&ph, ntotal);
+// ph = dvector(0, ntotal-1);
+ rt = dvector (0, ntotal - 1);
+ s = dvector (0, ntotal - 1);
+ allocate_derivs (&sh, ntotal);
+// sh = dvector(0, ntotal-1);
+ t = dvector (0, ntotal - 1);
+ vv = dvector (0, ntotal - 1);
+
+ /* check */
+ if (output == 1)
+ printf ("bicgstab: itmax %d, tol %e\n", itmax, tol);
+
+ /* compute initial residual rt = r = F - J*dv */
+ J_times_dv (nvar, n1, n2, n3, dv, r, u);
+ for (j = 0; j < ntotal; j++)
+ rt[j] = r[j] = F[j] - r[j];
+
+ *normres = norm2 (r, ntotal);
+ if (output == 1)
+ printf ("bicgstab: %5d %10.3e\n", 0, *normres);
+
+ if (*normres <= tol)
+ return 0;
+
+ /* cgs iteration */
+ for (ii = 0; ii < itmax; ii++)
+ {
+ rho = scalarproduct (rt, r, ntotal);
+ if (fabs (rho) < rhotol)
+ break;
+
+ /* compute direction vector p */
+ if (ii == 0)
+ for (j = 0; j < ntotal; j++)
+ p[j] = r[j];
+ else
+ {
+ beta = (rho / rho1) * (alpha / omega);
+ for (j = 0; j < ntotal; j++)
+ p[j] = r[j] + beta * (p[j] - omega * vv[j]);
+ }
+
+ /* compute direction adjusting vector ph and scalar alpha */
+ for (j = 0; j < ntotal; j++)
+ ph.d0[j] = 0;
+ for (j = 0; j < NRELAX; j++) // solves JFD*ph = p by relaxation
+ relax (ph.d0, nvar, n1, n2, n3, p, ncols, cols, JFD);
+
+ J_times_dv (nvar, n1, n2, n3, ph, vv, u); // vv=J*ph
+ alpha = rho / scalarproduct (rt, vv, ntotal);
+ for (j = 0; j < ntotal; j++)
+ s[j] = r[j] - alpha * vv[j];
+
+ /* early check of tolerance */
+ *normres = norm2 (s, ntotal);
+ if (*normres <= tol)
+ {
+ for (j = 0; j < ntotal; j++)
+ dv.d0[j] += alpha * ph.d0[j];
+ if (output == 1)
+ printf ("bicgstab: %5d %10.3e %10.3e %10.3e %10.3e\n",
+ ii + 1, *normres, alpha, beta, omega);
+ break;
+ }
+
+ /* compute stabilizer vector sh and scalar omega */
+ for (j = 0; j < ntotal; j++)
+ sh.d0[j] = 0;
+ for (j = 0; j < NRELAX; j++) // solves JFD*sh = s by relaxation
+ relax (sh.d0, nvar, n1, n2, n3, s, ncols, cols, JFD);
+
+ J_times_dv (nvar, n1, n2, n3, sh, t, u); // t=J*sh
+ omega = scalarproduct (t, s, ntotal) / scalarproduct (t, t, ntotal);
+
+ /* compute new solution approximation */
+ for (j = 0; j < ntotal; j++)
+ {
+ dv.d0[j] += alpha * ph.d0[j] + omega * sh.d0[j];
+ r[j] = s[j] - omega * t[j];
+ }
+ /* are we done? */
+ *normres = norm2 (r, ntotal);
+ if (output == 1)
+ printf ("bicgstab: %5d %10.3e %10.3e %10.3e %10.3e\n",
+ ii + 1, *normres, alpha, beta, omega);
+ if (*normres <= tol)
+ break;
+ rho1 = rho;
+ if (fabs (omega) < omegatol)
+ break;
+
+ }
+
+ /* free temporary storage */
+ free_dvector (r, 0, ntotal - 1);
+ free_dvector (p, 0, ntotal - 1);
+// free_dvector(ph, 0, ntotal-1);
+ free_derivs (&ph, ntotal);
+ free_dvector (rt, 0, ntotal - 1);
+ free_dvector (s, 0, ntotal - 1);
+// free_dvector(sh, 0, ntotal-1);
+ free_derivs (&sh, ntotal);
+ free_dvector (t, 0, ntotal - 1);
+ free_dvector (vv, 0, ntotal - 1);
+
+ free_dvector (F, 0, ntotal - 1);
+ free_derivs (&u, ntotal);
+
+ free_dmatrix (JFD, 0, ntotal - 1, 0, maxcol - 1);
+ free_imatrix (cols, 0, ntotal - 1, 0, maxcol - 1);
+ free_ivector (ncols, 0, ntotal - 1);
+
+ /* iteration failed */
+ if (ii > itmax)
+ return -1;
+
+ /* breakdown */
+ if (fabs (rho) < rhotol)
+ return -10;
+ if (fabs (omega) < omegatol)
+ return -11;
+
+ /* success! */
+ return ii + 1;
+}
+
+// -------------------------------------------------------------------
+void
+Newton (int nvar, int n1, int n2, int n3,
+ derivs v, double tol, int itmax)
+{
+ DECLARE_CCTK_PARAMETERS;
+ int ntotal = n1 * n2 * n3 * nvar, ii, j, it;
+ double *F, dmax, normres;
+ derivs u, dv;
+
+ F = dvector (0, ntotal - 1);
+ allocate_derivs (&dv, ntotal);
+ allocate_derivs (&u, ntotal);
+
+/* TestRelax(nvar, n1, n2, n3, v, dv.d0); */
+ for (j = 0; j < ntotal; j++)
+ v.d0[j] = 0;
+
+ it = 0;
+ dmax = 1;
+ while (dmax > tol && it < itmax)
+ {
+ if (it == 0)
+ {
+ F_of_v (nvar, n1, n2, n3, v, F, u);
+ dmax = -1;
+ for (j = 0; j < ntotal; j++)
+ {
+ if (fabs (F[j]) > dmax)
+ dmax = fabs (F[j]);
+ dv.d0[j] = 0;
+ }
+ }
+ printf ("Newton: it=%d \t |F|=%e \n", it, dmax);
+ ii =
+ bicgstab (nvar, n1, n2, n3, v, dv, 1, 100, dmax * 1.e-3, &normres);
+ for (j = 0; j < ntotal; j++)
+ v.d0[j] -= dv.d0[j];
+ F_of_v (nvar, n1, n2, n3, v, F, u);
+ dmax = -1;
+ for (j = 0; j < ntotal; j++)
+ {
+ if (fabs (F[j]) > dmax)
+ {
+ dmax = fabs (F[j]);
+ }
+ }
+ it += 1;
+ }
+ printf ("Newton: it=%d \t |F|=%e \n", it, dmax);
+
+ free_dvector (F, 0, ntotal - 1);
+ free_derivs (&dv, ntotal);
+ free_derivs (&u, ntotal);
+}
+
+// -------------------------------------------------------------------
diff --git a/src/TP_utilities.c b/src/TP_utilities.c
new file mode 100644
index 0000000..64edbfd
--- /dev/null
+++ b/src/TP_utilities.c
@@ -0,0 +1,951 @@
+// TwoPunctures: File "utilities.c"
+
+#include <math.h>
+#include <stdio.h>
+#include <stddef.h>
+#include <stdlib.h>
+#include "TP_utilities.h"
+
+//-----------------------------------------------------------------------------
+void
+nrerror (char error_text[])
+/* Numerical Recipes standard error handler */
+{
+ fprintf (stderr, "Numerical Recipes run-time error...\n");
+ fprintf (stderr, "%s\n", error_text);
+ fprintf (stderr, "...now exiting to system...\n");
+ exit (1);
+}
+
+//-----------------------------------------------------------------------------
+int *
+ivector (long nl, long nh)
+/* allocate an int vector with subscript range v[nl..nh] */
+{
+ int *v;
+
+ v = (int *) malloc ((size_t) ((nh - nl + 1 + NR_END) * sizeof (int)));
+ if (!v)
+ nrerror ("allocation failure in ivector()");
+ return v - nl + NR_END;
+}
+
+//-----------------------------------------------------------------------------
+double *
+dvector (long nl, long nh)
+/* allocate a double vector with subscript range v[nl..nh] */
+{
+ double *v;
+
+ v = (double *) malloc ((size_t) ((nh - nl + 1 + NR_END) * sizeof (double)));
+ if (!v)
+ nrerror ("allocation failure in dvector()");
+ return v - nl + NR_END;
+}
+
+//-----------------------------------------------------------------------------
+int **
+imatrix (long nrl, long nrh, long ncl, long nch)
+/* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */
+{
+ long i, nrow = nrh - nrl + 1, ncol = nch - ncl + 1;
+ int **m;
+
+ /* allocate pointers to rows */
+ m = (int **) malloc ((size_t) ((nrow + NR_END) * sizeof (int *)));
+ if (!m)
+ nrerror ("allocation failure 1 in matrix()");
+ m += NR_END;
+ m -= nrl;
+
+
+ /* allocate rows and set pointers to them */
+ m[nrl] = (int *) malloc ((size_t) ((nrow * ncol + NR_END) * sizeof (int)));
+ if (!m[nrl])
+ nrerror ("allocation failure 2 in matrix()");
+ m[nrl] += NR_END;
+ m[nrl] -= ncl;
+
+ for (i = nrl + 1; i <= nrh; i++)
+ m[i] = m[i - 1] + ncol;
+
+ /* return pointer to array of pointers to rows */
+ return m;
+}
+
+//-----------------------------------------------------------------------------
+double **
+dmatrix (long nrl, long nrh, long ncl, long nch)
+/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
+{
+ long i, nrow = nrh - nrl + 1, ncol = nch - ncl + 1;
+ double **m;
+
+ /* allocate pointers to rows */
+ m = (double **) malloc ((size_t) ((nrow + NR_END) * sizeof (double *)));
+ if (!m)
+ nrerror ("allocation failure 1 in matrix()");
+ m += NR_END;
+ m -= nrl;
+
+ /* allocate rows and set pointers to them */
+ m[nrl] =
+ (double *) malloc ((size_t) ((nrow * ncol + NR_END) * sizeof (double)));
+ if (!m[nrl])
+ nrerror ("allocation failure 2 in matrix()");
+ m[nrl] += NR_END;
+ m[nrl] -= ncl;
+
+ for (i = nrl + 1; i <= nrh; i++)
+ m[i] = m[i - 1] + ncol;
+
+ /* return pointer to array of pointers to rows */
+ return m;
+}
+
+//-----------------------------------------------------------------------------
+double ***
+d3tensor (long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
+/* allocate a double 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
+{
+ long i, j, nrow = nrh - nrl + 1, ncol = nch - ncl + 1, ndep = ndh - ndl + 1;
+ double ***t;
+
+ /* allocate pointers to pointers to rows */
+ t = (double ***) malloc ((size_t) ((nrow + NR_END) * sizeof (double **)));
+ if (!t)
+ nrerror ("allocation failure 1 in f3tensor()");
+ t += NR_END;
+ t -= nrl;
+
+ /* allocate pointers to rows and set pointers to them */
+ t[nrl] =
+ (double **)
+ malloc ((size_t) ((nrow * ncol + NR_END) * sizeof (double *)));
+ if (!t[nrl])
+ nrerror ("allocation failure 2 in f3tensor()");
+ t[nrl] += NR_END;
+ t[nrl] -= ncl;
+
+ /* allocate rows and set pointers to them */
+ t[nrl][ncl] =
+ (double *)
+ malloc ((size_t) ((nrow * ncol * ndep + NR_END) * sizeof (double)));
+ if (!t[nrl][ncl])
+ nrerror ("allocation failure 3 in f3tensor()");
+ t[nrl][ncl] += NR_END;
+ t[nrl][ncl] -= ndl;
+
+ for (j = ncl + 1; j <= nch; j++)
+ t[nrl][j] = t[nrl][j - 1] + ndep;
+ for (i = nrl + 1; i <= nrh; i++)
+ {
+ t[i] = t[i - 1] + ncol;
+ t[i][ncl] = t[i - 1][ncl] + ncol * ndep;
+ for (j = ncl + 1; j <= nch; j++)
+ t[i][j] = t[i][j - 1] + ndep;
+ }
+
+ /* return pointer to array of pointers to rows */
+ return t;
+}
+
+//-----------------------------------------------------------------------------
+void
+free_ivector (int *v, long nl, long nh)
+/* free an int vector allocated with ivector() */
+{
+ free ((FREE_ARG) (v + nl - NR_END));
+}
+
+//-----------------------------------------------------------------------------
+void
+free_dvector (double *v, long nl, long nh)
+/* free a double vector allocated with dvector() */
+{
+ free ((FREE_ARG) (v + nl - NR_END));
+}
+
+//-----------------------------------------------------------------------------
+void
+free_imatrix (int **m, long nrl, long nrh, long ncl, long nch)
+/* free an int matrix allocated by imatrix() */
+{
+ free ((FREE_ARG) (m[nrl] + ncl - NR_END));
+ free ((FREE_ARG) (m + nrl - NR_END));
+}
+
+//-----------------------------------------------------------------------------
+void
+free_dmatrix (double **m, long nrl, long nrh, long ncl, long nch)
+/* free a double matrix allocated by dmatrix() */
+{
+ free ((FREE_ARG) (m[nrl] + ncl - NR_END));
+ free ((FREE_ARG) (m + nrl - NR_END));
+}
+
+//-----------------------------------------------------------------------------
+void
+free_d3tensor (double ***t, long nrl, long nrh, long ncl, long nch,
+ long ndl, long ndh)
+/* free a double f3tensor allocated by f3tensor() */
+{
+ free ((FREE_ARG) (t[nrl][ncl] + ndl - NR_END));
+ free ((FREE_ARG) (t[nrl] + ncl - NR_END));
+ free ((FREE_ARG) (t + nrl - NR_END));
+}
+
+//-----------------------------------------------------------------------------
+int
+minimum2 (int i, int j)
+{
+ int result = i;
+ if (j < result)
+ result = j;
+ return result;
+}
+
+//-----------------------------------------------------------------------------
+int
+minimum3 (int i, int j, int k)
+{
+ int result = i;
+ if (j < result)
+ result = j;
+ if (k < result)
+ result = k;
+ return result;
+}
+
+//-----------------------------------------------------------------------------
+int
+maximum2 (int i, int j)
+{
+ int result = i;
+ if (j > result)
+ result = j;
+ return result;
+}
+
+//-----------------------------------------------------------------------------
+int
+maximum3 (int i, int j, int k)
+{
+ int result = i;
+ if (j > result)
+ result = j;
+ if (k > result)
+ result = k;
+ return result;
+}
+
+//-----------------------------------------------------------------------------
+int
+pow_int (int mantisse, int exponent)
+{
+ int i, result = 1;
+
+ for (i = 1; i <= exponent; i++)
+ result *= mantisse;
+
+ return result;
+}
+
+//-----------------------------------------------------------------------------
+double
+atanh (double x)
+{
+ return 0.5 * log ((1 + x) / (1 - x));
+}
+
+//-----------------------------------------------------------------------------
+double
+asinh (double x)
+{
+ return log (x + sqrt (1 + x * x));
+}
+
+//-----------------------------------------------------------------------------
+double
+acosh (double x)
+{
+ return log (x + sqrt (x * x - 1));
+}
+
+//-----------------------------------------------------------------------------
+dcomplex
+Cadd (dcomplex a, dcomplex b)
+{
+ dcomplex c;
+ c.r = a.r + b.r;
+ c.i = a.i + b.i;
+ return c;
+}
+
+//-----------------------------------------------------------------------------
+dcomplex
+Csub (dcomplex a, dcomplex b)
+{
+ dcomplex c;
+ c.r = a.r - b.r;
+ c.i = a.i - b.i;
+ return c;
+}
+
+//-----------------------------------------------------------------------------
+dcomplex
+Cmul (dcomplex a, dcomplex b)
+{
+ dcomplex c;
+ c.r = a.r * b.r - a.i * b.i;
+ c.i = a.i * b.r + a.r * b.i;
+ return c;
+}
+
+//-----------------------------------------------------------------------------
+dcomplex
+RCmul (double x, dcomplex a)
+{
+ dcomplex c;
+ c.r = x * a.r;
+ c.i = x * a.i;
+ return c;
+}
+
+//-----------------------------------------------------------------------------
+dcomplex
+Cdiv (dcomplex a, dcomplex b)
+{
+ dcomplex c;
+ double r, den;
+ if (fabs (b.r) >= fabs (b.i))
+ {
+ r = b.i / b.r;
+ den = b.r + r * b.i;
+ c.r = (a.r + r * a.i) / den;
+ c.i = (a.i - r * a.r) / den;
+ }
+ else
+ {
+ r = b.r / b.i;
+ den = b.i + r * b.r;
+ c.r = (a.r * r + a.i) / den;
+ c.i = (a.i * r - a.r) / den;
+ }
+ return c;
+}
+
+//-----------------------------------------------------------------------------
+dcomplex
+Complex (double re, double im)
+{
+ dcomplex c;
+ c.r = re;
+ c.i = im;
+ return c;
+}
+
+//-----------------------------------------------------------------------------
+dcomplex
+Conjg (dcomplex z)
+{
+ dcomplex c;
+ c.r = z.r;
+ c.i = -z.i;
+ return c;
+}
+
+//-----------------------------------------------------------------------------
+double
+Cabs (dcomplex z)
+{
+ double x, y, ans, temp;
+ x = fabs (z.r);
+ y = fabs (z.i);
+ if (x == 0.0)
+ ans = y;
+ else if (y == 0.0)
+ ans = x;
+ else if (x > y)
+ {
+ temp = y / x;
+ ans = x * sqrt (1.0 + temp * temp);
+ }
+ else
+ {
+ temp = x / y;
+ ans = y * sqrt (1.0 + temp * temp);
+ }
+ return ans;
+}
+
+//-----------------------------------------------------------------------------
+dcomplex
+Csqrt (dcomplex z)
+{
+ dcomplex c;
+ double x, y, w, r;
+ if ((z.r == 0.0) && (z.i == 0.0))
+ {
+ c.r = 0.0;
+ c.i = 0.0;
+ return c;
+ }
+ else
+ {
+ x = fabs (z.r);
+ y = fabs (z.i);
+ if (x >= y)
+ {
+ r = y / x;
+ w = sqrt (x) * sqrt (0.5 * (1.0 + sqrt (1.0 + r * r)));
+ }
+ else
+ {
+ r = x / y;
+ w = sqrt (y) * sqrt (0.5 * (r + sqrt (1.0 + r * r)));
+ }
+ if (z.r >= 0.0)
+ {
+ c.r = w;
+ c.i = z.i / (2.0 * w);
+ }
+ else
+ {
+ c.i = (z.i >= 0) ? w : -w;
+ c.r = z.i / (2.0 * c.i);
+ }
+ return c;
+ }
+}
+
+//-----------------------------------------------------------------------------
+dcomplex
+Cexp (dcomplex z)
+{
+ dcomplex c;
+ double exp_r = exp (z.r);
+
+ c.r = exp_r * cos (z.i);
+ c.i = exp_r * sin (z.i);
+
+ return c;
+}
+
+//-----------------------------------------------------------------------------
+dcomplex
+Clog (dcomplex z)
+{
+ dcomplex c;
+
+ c.r = 0.5 * log (z.r * z.r + z.i * z.i);
+ c.i = atan2 (z.i, z.r);
+
+ return c;
+}
+
+//-----------------------------------------------------------------------------
+dcomplex
+Csin (dcomplex z)
+{
+ dcomplex c;
+
+ c.r = sin (z.r) * cosh (z.i);
+ c.i = cos (z.r) * sinh (z.i);
+
+ return c;
+} //-----------------------------------------------------------------------------
+
+dcomplex
+Ccos (dcomplex z)
+{
+ dcomplex c;
+
+ c.r = cos (z.r) * cosh (z.i);
+ c.i = -sin (z.r) * sinh (z.i);
+
+ return c;
+}
+
+//-----------------------------------------------------------------------------
+dcomplex
+Ctan (dcomplex z)
+{
+ return Cdiv (Csin (z), Ccos (z));
+}
+
+//-----------------------------------------------------------------------------
+dcomplex
+Ccot (dcomplex z)
+{
+ return Cdiv (Ccos (z), Csin (z));
+}
+
+//-----------------------------------------------------------------------------
+dcomplex
+Csinh (dcomplex z)
+{
+ dcomplex c;
+
+ c.r = sinh (z.r) * cos (z.i);
+ c.i = cosh (z.r) * sin (z.i);
+
+ return c;
+} //-----------------------------------------------------------------------------
+
+dcomplex
+Ccosh (dcomplex z)
+{
+ dcomplex c;
+
+ c.r = cosh (z.r) * cos (z.i);
+ c.i = sinh (z.r) * sin (z.i);
+
+ return c;
+}
+
+//-----------------------------------------------------------------------------
+dcomplex
+Ctanh (dcomplex z)
+{
+ return Cdiv (Csinh (z), Ccosh (z));
+}
+
+//-----------------------------------------------------------------------------
+dcomplex
+Ccoth (dcomplex z)
+{
+ return Cdiv (Ccosh (z), Csinh (z));
+}
+
+//-----------------------------------------------------------------------------
+void
+chebft_Zeros (double u[], int n, int inv)
+{
+ int k, j, isignum;
+ double fac, sum, Pion, *c;
+
+ c = dvector (0, n);
+ Pion = Pi / n;
+ if (inv == 0)
+ {
+ fac = 2.0 / n;
+ isignum = 1;
+ for (j = 0; j < n; j++)
+ {
+ sum = 0.0;
+ for (k = 0; k < n; k++)
+ sum += u[k] * cos (Pion * j * (k + 0.5));
+ c[j] = fac * sum * isignum;
+ isignum = -isignum;
+ }
+ }
+ else
+ {
+ for (j = 0; j < n; j++)
+ {
+ sum = -0.5 * u[0];
+ isignum = 1;
+ for (k = 0; k < n; k++)
+ {
+ sum += u[k] * cos (Pion * (j + 0.5) * k) * isignum;
+ isignum = -isignum;
+ }
+ c[j] = sum;
+ }
+ }
+ for (j = 0; j < n; j++)
+ u[j] = c[j];
+ free_dvector (c, 0, n);
+}
+
+// -------------------------------------------------------------------------------
+
+void
+chebft_Extremes (double u[], int n, int inv)
+{
+ int k, j, isignum, N = n - 1;
+ double fac, sum, PioN, *c;
+
+ c = dvector (0, N);
+ PioN = Pi / N;
+ if (inv == 0)
+ {
+ fac = 2.0 / N;
+ isignum = 1;
+ for (j = 0; j < n; j++)
+ {
+ sum = 0.5 * (u[0] + u[N] * isignum);
+ for (k = 1; k < N; k++)
+ sum += u[k] * cos (PioN * j * k);
+ c[j] = fac * sum * isignum;
+ isignum = -isignum;
+ }
+ c[N] = 0.5 * c[N];
+ }
+ else
+ {
+ for (j = 0; j < n; j++)
+ {
+ sum = -0.5 * u[0];
+ isignum = 1;
+ for (k = 0; k < n; k++)
+ {
+ sum += u[k] * cos (PioN * j * k) * isignum;
+ isignum = -isignum;
+ }
+ c[j] = sum;
+ }
+ }
+ for (j = 0; j < n; j++)
+ u[j] = c[j];
+ free_dvector (c, 0, N);
+}
+
+// -------------------------------------------------------------------------------
+
+void
+chder (double *c, double *cder, int n)
+{
+ int j;
+
+ cder[n] = 0.0;
+ cder[n - 1] = 0.0;
+ for (j = n - 2; j >= 0; j--)
+ cder[j] = cder[j + 2] + 2 * (j + 1) * c[j + 1];
+}
+
+// -------------------------------------------------------------------------------
+double
+chebev (double a, double b, double c[], int m, double x)
+{
+ double d = 0.0, dd = 0.0, sv, y, y2;
+ int j;
+
+ y2 = 2.0 * (y = (2.0 * x - a - b) / (b - a));
+ for (j = m - 1; j >= 1; j--)
+ {
+ sv = d;
+ d = y2 * d - dd + c[j];
+ dd = sv;
+ }
+ return y * d - dd + 0.5 * c[0];
+}
+
+// -------------------------------------------------------------------------------
+void
+fourft (double *u, int N, int inv)
+{
+ int l, k, iy, M;
+ double x, x1, fac, Pi_fac, *a, *b;
+
+ M = N / 2;
+ a = dvector (0, M);
+ b = dvector (1, M); // Actually: b=vector(1,M-1) but this is problematic if M=1
+ fac = 1. / M;
+ Pi_fac = Pi * fac;
+ if (inv == 0)
+ {
+ for (l = 0; l <= M; l++)
+ {
+ a[l] = 0;
+ if (l > 0 && l < M)
+ b[l] = 0;
+ x1 = Pi_fac * l;
+ for (k = 0; k < N; k++)
+ {
+ x = x1 * k;
+ a[l] += fac * u[k] * cos (x);
+ if (l > 0 && l < M)
+ b[l] += fac * u[k] * sin (x);
+ }
+ }
+ u[0] = a[0];
+ u[M] = a[M];
+ for (l = 1; l < M; l++)
+ {
+ u[l] = a[l];
+ u[l + M] = b[l];
+ }
+ }
+ else
+ {
+ a[0] = u[0];
+ a[M] = u[M];
+ for (l = 1; l < M; l++)
+ {
+ a[l] = u[l];
+ b[l] = u[M + l];
+ }
+ iy = 1;
+ for (k = 0; k < N; k++)
+ {
+ u[k] = 0.5 * (a[0] + a[M] * iy);
+ x1 = Pi_fac * k;
+ for (l = 1; l < M; l++)
+ {
+ x = x1 * l;
+ u[k] += a[l] * cos (x) + b[l] * sin (x);
+ }
+ iy = -iy;
+ }
+ }
+ free_dvector (a, 0, M);
+ free_dvector (b, 1, M);
+}
+
+/* -----------------------------------------*/
+void
+fourder (double u[], double du[], int N)
+{
+ int l, M, lpM;
+
+ M = N / 2;
+ du[0] = 0.;
+ du[M] = 0.;
+ for (l = 1; l < M; l++)
+ {
+ lpM = l + M;
+ du[l] = u[lpM] * l;
+ du[lpM] = -u[l] * l;
+ }
+}
+
+/* -----------------------------------------*/
+void
+fourder2 (double u[], double d2u[], int N)
+{
+ int l, l2, M, lpM;
+
+ d2u[0] = 0.;
+ M = N / 2;
+ for (l = 1; l <= M; l++)
+ {
+ l2 = l * l;
+ lpM = l + M;
+ d2u[l] = -u[l] * l2;
+ if (l < M)
+ d2u[lpM] = -u[lpM] * l2;
+ }
+}
+
+/* ----------------------------------------- */
+double
+fourev (double *u, int N, double x)
+{
+ int l, M = N / 2;
+ double xl, result;
+
+ result = 0.5 * (u[0] + u[M] * cos (x * M));
+ for (l = 1; l < M; l++)
+ {
+ xl = x * l;
+ result += u[l] * cos (xl) + u[M + l] * sin (xl);
+ }
+ return result;
+}
+
+// -------------------------------------------------------------------------------
+void
+ludcmp (double **a, int n, int *indx, double *d)
+{ // Version of 'ludcmp' of the numerical recipes for
+ // matrices a[0:n-1][0:n-1]
+ int i, imax, j, k;
+ double big, dum, sum, temp;
+ double *vv;
+
+ vv = dvector (0, n - 1);
+ *d = 1.0;
+ for (i = 0; i < n; i++)
+ {
+ big = 0.0;
+ for (j = 0; j < n; j++)
+ if ((temp = fabs (a[i][j])) > big)
+ big = temp;
+ if (big == 0.0)
+ nrerror ("Singular matrix in routine ludcmp");
+ vv[i] = 1.0 / big;
+ }
+ for (j = 0; j < n; j++)
+ {
+ for (i = 0; i < j; i++)
+ {
+ sum = a[i][j];
+ for (k = 0; k < i; k++)
+ sum -= a[i][k] * a[k][j];
+ a[i][j] = sum;
+ }
+ big = 0.0;
+ for (i = j; i < n; i++)
+ {
+ sum = a[i][j];
+ for (k = 0; k < j; k++)
+ sum -= a[i][k] * a[k][j];
+ a[i][j] = sum;
+ if ((dum = vv[i] * fabs (sum)) >= big)
+ {
+ big = dum;
+ imax = i;
+ }
+ }
+ if (j != imax)
+ {
+ for (k = 0; k < n; k++)
+ {
+ dum = a[imax][k];
+ a[imax][k] = a[j][k];
+ a[j][k] = dum;
+ }
+ *d = -(*d);
+ vv[imax] = vv[j];
+ }
+ indx[j] = imax;
+ if (a[j][j] == 0.0)
+ a[j][j] = TINY;
+ if (j != n)
+ {
+ dum = 1.0 / (a[j][j]);
+ for (i = j + 1; i < n; i++)
+ a[i][j] *= dum;
+ }
+ }
+ free_dvector (vv, 0, n - 1);
+}
+
+// -------------------------------------------------------------------------------
+void
+lubksb (double **a, int n, int *indx, double b[])
+{ // Version of 'lubksb' of the numerical recipes for
+ // matrices a[0:n-1][0:n-1] and vectors b[0:n-1]
+
+ int i, ii = 0, ip, j;
+ double sum;
+
+ for (i = 0; i < n; i++)
+ {
+ ip = indx[i];
+ sum = b[ip];
+ b[ip] = b[i];
+ if (ii)
+ for (j = ii; j <= i - 1; j++)
+ sum -= a[i][j] * b[j];
+ else if (sum)
+ ii = i;
+ b[i] = sum;
+ }
+ for (i = n - 1; i >= 0; i--)
+ {
+ sum = b[i];
+ for (j = i + 1; j <= n; j++)
+ sum -= a[i][j] * b[j];
+ b[i] = sum / a[i][i];
+ }
+}
+
+// -----------------------------------------------------------------------------------
+void
+tridag (double a[], double b[], double c[], double r[], double u[], int n)
+{ // Version of 'tridag' of the numerical recipes for
+ // vectors a, b, c, r, u with indices in the range [0:n-1]
+ int j;
+ double bet, *gam;
+
+ gam = dvector (0, n - 1);
+ if (b[0] == 0.0)
+ nrerror ("Error 1 in tridag");
+ u[0] = r[0] / (bet = b[0]);
+ for (j = 1; j < n; j++)
+ {
+ gam[j] = c[j - 1] / bet;
+ bet = b[j] - a[j] * gam[j];
+ if (bet == 0.0)
+ nrerror ("Error 2 in tridag");
+ u[j] = (r[j] - a[j] * u[j - 1]) / bet;
+ }
+ for (j = (n - 2); j >= 0; j--)
+ u[j] -= gam[j + 1] * u[j + 1];
+ free_dvector (gam, 0, n - 1);
+}
+
+// -----------------------------------------------------------------------------------
+double
+norm1 (double *v, int n)
+{
+ int i;
+ double result = -1;
+
+ for (i = 0; i < n; i++)
+ if (fabs (v[i]) > result)
+ result = fabs (v[i]);
+
+ return result;
+}
+
+// -----------------------------------------------------------------------------------
+double
+norm2 (double *v, int n)
+{
+ int i;
+ double result = 0;
+
+ for (i = 0; i < n; i++)
+ result += v[i] * v[i];
+
+ return sqrt (result);
+}
+
+// -----------------------------------------------------------------------------------
+double
+scalarproduct (double *v, double *w, int n)
+{
+ int i;
+ double result = 0;
+
+ for (i = 0; i < n; i++)
+ result += v[i] * w[i];
+
+ return result;
+}
+
+// -----------------------------------------------------------------------------------
+double
+plgndr (int l, int m, double x)
+{
+ void nrerror (char error_text[]);
+ double fact, pll, pmm, pmmp1, somx2;
+ int i, ll;
+
+ if (m < 0 || m > l || fabs (x) > 1.0)
+ nrerror ("Bad arguments in routine plgndr");
+ pmm = 1.0;
+ if (m > 0)
+ {
+ somx2 = sqrt ((1.0 - x) * (1.0 + x));
+ fact = 1.0;
+ for (i = 1; i <= m; i++)
+ {
+ pmm *= -fact * somx2;
+ fact += 2.0;
+ }
+ }
+ if (l == m)
+ return pmm;
+ else
+ {
+ pmmp1 = x * (2 * m + 1) * pmm;
+ if (l == (m + 1))
+ return pmmp1;
+ else
+ {
+ for (ll = m + 2; ll <= l; ll++)
+ {
+ pll = (x * (2 * ll - 1) * pmmp1 - (ll + m - 1) * pmm) / (ll - m);
+ pmm = pmmp1;
+ pmmp1 = pll;
+ }
+ return pll;
+ }
+ }
+}
+
+// -----------------------------------------------------------------------------------
diff --git a/src/TP_utilities.h b/src/TP_utilities.h
new file mode 100644
index 0000000..37db3f0
--- /dev/null
+++ b/src/TP_utilities.h
@@ -0,0 +1,91 @@
+// TwoPunctures: File "utilities.h"
+
+#define Pi 3.14159265358979323846264338328
+#define Pih 1.57079632679489661923132169164 // Pi/2
+#define Piq 0.78539816339744830961566084582 // Pi/4
+
+#define TINY 1.0e-20
+#define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}
+#define NR_END 1
+#define FREE_ARG char*
+
+typedef struct DCOMPLEX
+{
+ double r, i;
+} dcomplex;
+
+#define nrerror TP_nrerror
+#define ivector TP_ivector
+#define dvector TP_dvector
+#define imatrix TP_imatrix
+#define dmatrix TP_dmatrix
+#define d3tensor TP_d3tensor
+#define free_ivector TP_free_ivector
+#define free_dvector TP_free_dvector
+#define free_imatrix TP_free_imatrix
+#define free_dmatrix TP_free_dmatrix
+#define free_d3tensor TP_free_d3tensor
+
+void nrerror (char error_text[]);
+int *ivector (long nl, long nh);
+double *dvector (long nl, long nh);
+int **imatrix (long nrl, long nrh, long ncl, long nch);
+double **dmatrix (long nrl, long nrh, long ncl, long nch);
+double ***d3tensor (long nrl, long nrh, long ncl, long nch, long ndl,
+ long ndh);
+void free_ivector (int *v, long nl, long nh);
+void free_dvector (double *v, long nl, long nh);
+void free_imatrix (int **m, long nrl, long nrh, long ncl, long nch);
+void free_dmatrix (double **m, long nrl, long nrh, long ncl, long nch);
+void free_d3tensor (double ***t, long nrl, long nrh, long ncl, long nch,
+ long ndl, long ndh);
+
+int minimum2 (int i, int j);
+int minimum3 (int i, int j, int k);
+int maximum2 (int i, int j);
+int maximum3 (int i, int j, int k);
+int pow_int (int mantisse, int exponent);
+double atanh (double x);
+double asinh (double x);
+double acosh (double x);
+
+dcomplex Cadd (dcomplex a, dcomplex b);
+dcomplex Csub (dcomplex a, dcomplex b);
+dcomplex Cmul (dcomplex a, dcomplex b);
+dcomplex RCmul (double x, dcomplex a);
+dcomplex Cdiv (dcomplex a, dcomplex b);
+dcomplex Complex (double re, double im);
+dcomplex Conjg (dcomplex z);
+double Cabs (dcomplex z);
+
+dcomplex Csqrt (dcomplex z);
+dcomplex Cexp (dcomplex z);
+dcomplex Clog (dcomplex z);
+dcomplex Csin (dcomplex z);
+dcomplex Ccos (dcomplex z);
+dcomplex Ctan (dcomplex z);
+dcomplex Ccot (dcomplex z);
+dcomplex Csinh (dcomplex z);
+dcomplex Ccosh (dcomplex z);
+dcomplex Ctanh (dcomplex z);
+dcomplex Ccoth (dcomplex z);
+
+void chebft_Zeros (double u[], int n, int inv);
+void chebft_Extremes (double u[], int n, int inv);
+void chder (double *c, double *cder, int n);
+double chebev (double a, double b, double c[], int m, double x);
+void fourft (double *u, int N, int inv);
+void fourder (double u[], double du[], int N);
+void fourder2 (double u[], double d2u[], int N);
+double fourev (double *u, int N, double x);
+
+
+void ludcmp (double **a, int n, int *indx, double *d);
+void lubksb (double **a, int n, int *indx, double b[]);
+void tridag (double a[], double b[], double c[], double r[], double u[],
+ int n);
+double norm1 (double *v, int n);
+double norm2 (double *v, int n);
+double scalarproduct (double *v, double *w, int n);
+
+double plgndr (int l, int m, double x);
diff --git a/src/TwoPunctures.c b/src/TwoPunctures.c
new file mode 100644
index 0000000..a43f22b
--- /dev/null
+++ b/src/TwoPunctures.c
@@ -0,0 +1,128 @@
+// TwoPunctures: File "TwoPunctures.c"
+
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <ctype.h>
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+#include "TP_utilities.h"
+#include "TwoPunctures.h"
+
+static inline double pow2 (const double x)
+{
+ return x*x;
+}
+
+// -------------------------------------------------------------------
+void
+TwoPunctures (CCTK_ARGUMENTS)
+{
+ DECLARE_CCTK_ARGUMENTS;
+ DECLARE_CCTK_PARAMETERS;
+
+ int nvar = 1, n1 = npoints_A, n2 = npoints_B, n3 = npoints_phi;
+
+ int i, j, k, ntotal = n1 * n2 * n3 * nvar;
+ double *F;
+ derivs u, v;
+
+ F = dvector (0, ntotal - 1);
+ allocate_derivs (&u, ntotal);
+ allocate_derivs (&v, ntotal);
+
+ CCTK_INFO ("Beginning elliptic solving");
+ Newton (nvar, n1, n2, n3, v, 1.e-10, 5);
+
+ F_of_v (nvar, n1, n2, n3, v, F, u);
+
+ CCTK_INFO ("Interpolating result");
+ if (CCTK_EQUALS(metric_type, "static conformal")) {
+ if (CCTK_EQUALS(conformal_storage, "factor")) {
+ *conformal_state = 1;
+ } else if (CCTK_EQUALS(conformal_storage, "factor+derivs")) {
+ *conformal_state = 2;
+ } else if (CCTK_EQUALS(conformal_storage, "factor+derivs+2nd derivs")) {
+ *conformal_state = 3;
+ }
+ } else {
+ *conformal_state = 0;
+ }
+
+ for (k = 0; k < cctk_lsh[2]; ++k)
+ {
+ for (j = 0; j < cctk_lsh[1]; ++j)
+ {
+ for (i = 0; i < cctk_lsh[0]; ++i)
+ {
+
+ const int ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
+
+ const double r_plus
+ = sqrt(pow2(x[ind] - par_b) + pow2(y[ind]) + pow2(z[ind]));
+ const double r_minus
+ = sqrt(pow2(x[ind] + par_b) + pow2(y[ind]) + pow2(z[ind]));
+
+ const double U = PunctIntPolAtArbitPosition
+ (0, nvar, n1, n2, n3, v, x[ind], y[ind], z[ind]);
+ const double psi1 = 1 + 0.5 * par_m_plus / r_plus
+ + 0.5 * par_m_minus / r_minus + U;
+ const double psi4 = pow2(pow2(psi1));
+ const double psim2 = 1.0/pow2(psi1);
+
+ double Aij[3][3];
+ BY_Aijofxyz (x[ind], y[ind], z[ind], Aij);
+
+ switch (*conformal_state) {
+ case 0:
+ gxx[ind] = psi4;
+ gxy[ind] = 0;
+ gxz[ind] = 0;
+ gyy[ind] = psi4;
+ gyz[ind] = 0;
+ gzz[ind] = psi4;
+ break;
+
+ case 3:
+ /* not yet supported */
+ assert (0);
+ /* fall through */
+
+ case 2:
+ /* not yet supported */
+ assert (0);
+ /* fall through */
+
+ case 1:
+ psi[ind] = psi1;
+
+ gxx[ind] = 1;
+ gxy[ind] = 0;
+ gxz[ind] = 0;
+ gyy[ind] = 1;
+ gyz[ind] = 0;
+ gzz[ind] = 1;
+ break;
+
+ default:
+ assert(0);
+ }
+
+ kxx[ind] = psim2 * Aij[0][0];
+ kxy[ind] = psim2 * Aij[0][1];
+ kxz[ind] = psim2 * Aij[0][2];
+ kyy[ind] = psim2 * Aij[1][1];
+ kyz[ind] = psim2 * Aij[1][2];
+ kzz[ind] = psim2 * Aij[2][2];
+
+ }
+ }
+ }
+
+ free_dvector (F, 0, ntotal - 1);
+ free_derivs (&u, ntotal);
+ free_derivs (&v, ntotal);
+}
diff --git a/src/TwoPunctures.h b/src/TwoPunctures.h
new file mode 100644
index 0000000..ab0d841
--- /dev/null
+++ b/src/TwoPunctures.h
@@ -0,0 +1,95 @@
+// TwoPunctures: File "TwoPunctures.h"
+
+#define StencilSize 19
+#define N_PlaneRelax 1
+#define NRELAX 200
+#define Step_Relax 1
+
+typedef struct DERIVS
+{
+ double *d0, *d1, *d2, *d3, *d11, *d12, *d13, *d22, *d23, *d33;
+} derivs;
+
+/*
+Files of "TwoPunctures":
+ TwoPunctures.c
+ FuncAndJacobian.c
+ CoordTransf.c
+ Equations.c
+ Newton.c
+ utilities.c (see utilities.h)
+**************************
+*/
+
+// Routines in "TwoPunctures.c"
+double TestSolution (double A, double B, double X, double R, double phi);
+void TestVector_w (double *par, int nvar, int n1, int n2, int n3, double *w);
+
+// Routines in "FuncAndJacobian.c"
+int Index (int ivar, int i, int j, int k, int nvar, int n1, int n2, int n3);
+void allocate_derivs (derivs * v, int n);
+void free_derivs (derivs * v, int n);
+void Derivatives_AB3 (int nvar, int n1, int n2, int n3, derivs v);
+void F_of_v (int nvar, int n1, int n2, int n3, derivs v,
+ double *F, derivs u);
+void J_times_dv (int nvar, int n1, int n2, int n3, derivs dv,
+ double *Jdv, derivs u);
+void JFD_times_dv (int i, int j, int k, int nvar, int n1,
+ int n2, int n3, derivs dv, derivs u, double *values);
+void SetMatrix_JFD (int nvar, int n1, int n2, int n3,
+ derivs u, int *ncols, int **cols, double **Matrix);
+double PunctEvalAtArbitPosition (double *v, double A, double B, double phi,
+ int n1, int n2, int n3);
+void calculate_derivs (int i, int j, int k, int ivar, int nvar, int n1,
+ int n2, int n3, derivs v, derivs vv);
+double interpol (double a, double b, double c, derivs v);
+double PunctIntPolAtArbitPosition (int ivar, int nvar, int n1,
+ int n2, int n3, derivs v, double x,
+ double y, double z);
+
+// Routines in "CoordTransf.c"
+void AB_To_XR (int nvar, double A, double B, double *X,
+ double *R, derivs U);
+void C_To_c (int nvar, double X, double R, double *x,
+ double *r, derivs U);
+void rx3_To_xyz (int nvar, double x, double r, double phi, double *y,
+ double *z, derivs U);
+
+// Routines in "Equations.c"
+double BY_KKofxyz (double x, double y, double z);
+void BY_Aijofxyz (double x, double y, double z, double Aij[3][3]);
+void NonLinEquations (double A, double B, double X, double R,
+ double x, double r, double phi,
+ double y, double z, derivs U, double *values);
+void LinEquations (double A, double B, double X, double R,
+ double x, double r, double phi,
+ double y, double z, derivs dU, derivs U, double *values);
+
+// Routines in "Newton.c"
+void resid (double *res, int ntotal, double *dv, double *rhs,
+ int *ncols, int **cols, double **JFD);
+void LineRelax_al (double *dv, int j, int k, int nvar, int n1, int n2, int n3,
+ double *rhs, int *ncols, int **cols, double **JFD);
+void LineRelax_be (double *dv, int i, int k, int nvar, int n1, int n2, int n3,
+ double *rhs, int *ncols, int **cols, double **JFD);
+void relax (double *dv, int nvar, int n1, int n2, int n3, double *rhs,
+ int *ncols, int **cols, double **JFD);
+void TestRelax (int nvar, int n1, int n2, int n3, derivs v,
+ double *dv);
+int bicgstab (int nvar, int n1, int n2, int n3, derivs v,
+ derivs dv, int output, int itmax, double tol, double *normres);
+void Newton (int nvar, int n1, int n2, int n3, derivs v,
+ double tol, int itmax);
+
+
+/*
+ 27: -1.325691774825335e-03
+ 37: -1.325691778944117e-03
+ 47: -1.325691778942711e-03
+
+ 17: -1.510625972641537e-03
+ 21: -1.511443006977708e-03
+ 27: -1.511440785153687e-03
+ 37: -1.511440809549005e-03
+ 39: -1.511440809597588e-03
+ */
diff --git a/src/make.code.defn b/src/make.code.defn
new file mode 100644
index 0000000..726aec8
--- /dev/null
+++ b/src/make.code.defn
@@ -0,0 +1,8 @@
+# Main make.code.defn file for thorn TwoPunctures
+# $Header$
+
+# Source files in this directory
+SRCS = CoordTransf.c Equations.c FuncAndJacobian.c Newton.c TwoPunctures.c TP_utilities.c
+
+# Subdirectories containing source files
+SUBDIRS =
diff --git a/test/twopunctures.par b/test/twopunctures.par
new file mode 100644
index 0000000..09c5676
--- /dev/null
+++ b/test/twopunctures.par
@@ -0,0 +1,53 @@
+# $Header$
+
+ActiveThorns = "Boundary CartGrid3D Time CoordBase SymBase PUGH PUGHReduce PUGHSlab IOASCII IOBasic IOUtil ADMBase StaticConformal TwoPunctures"
+
+Cactus::cctk_itlast = 0
+
+Time::dtfac = 0.25
+
+Grid::type = "byrange"
+Grid::domain = "full"
+Grid::xmin = -4
+Grid::xmax = 4
+Grid::ymin = -4
+Grid::ymax = 4
+Grid::zmin = -4
+Grid::zmax = 4
+Driver::global_nx = 35
+Driver::global_ny = 35
+Driver::global_nz = 35
+
+ADMBase::initial_data = "twopunctures"
+
+TwoPunctures::npoints_A = 20
+TwoPunctures::npoints_B = 20
+TwoPunctures::npoints_phi = 12
+
+TwoPunctures::par_b = 1.5
+
+TwoPunctures::par_m_plus = 1.5
+TwoPunctures::par_P_plus[1] = 2.0
+TwoPunctures::par_S_plus[1] = 0.5
+TwoPunctures::par_S_plus[2] = 0.5
+
+TwoPunctures::par_m_minus = 1.0
+TwoPunctures::par_P_minus[1] = -2.0
+TwoPunctures::par_S_minus[0] = -1.0
+TwoPunctures::par_S_minus[2] = 1.0
+
+ADMBase::metric_type = "static conformal"
+#StaticConformal::conformal_storage = "factor+derivs+2nd derivs"
+StaticConformal::conformal_storage = "factor"
+
+ADMBase::lapse_evolution_method = "static"
+ADMBase::initial_lapse = "one"
+
+IO::out_dir = $parfile
+IO::out_fileinfo = "axis labels"
+
+IOBasic::outScalar_every = 1
+IOBasic::outScalar_vars = "staticconformal::psi admbase::metric admbase::curv"
+
+IOASCII::out1D_every = 1
+IOASCII::out1D_vars = "staticconformal::psi admbase::metric admbase::curv"
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-06 16:01:21 +0000