From 748569b8a25ce2d5d3b774f18df0345a3fe68768 Mon Sep 17 00:00:00 2001 From: knarf Date: Tue, 5 Oct 2004 11:07:13 +0000 Subject: - add source terms provided by aliased function - small fix suggested by Marcus concerning the initialisation of the dvs git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/TwoPunctures/trunk@19 b2a53a04-0f4f-0410-87ed-f9f25ced00cf --- interface.ccl | 12 +++- param.ccl | 7 +++ schedule.ccl | 5 +- src/Equations.c | 6 +- src/FuncAndJacobian.c | 159 +++++++++++++++++++++++++++++++++++--------------- src/Newton.c | 35 ++++++----- src/ParamCheck.c | 32 ++++++++++ src/TwoPunctures.c | 4 +- src/TwoPunctures.h | 18 +++--- src/make.code.defn | 2 +- 10 files changed, 201 insertions(+), 79 deletions(-) create mode 100644 src/ParamCheck.c diff --git a/interface.ccl b/interface.ccl index a969173..3bb592b 100644 --- a/interface.ccl +++ b/interface.ccl @@ -5,6 +5,14 @@ IMPLEMENTS: TwoPunctures INHERITS: ADMBase StaticConformal grid - - REAL puncture_u TYPE=gf + +CCTK_INT FUNCTION Rho_ADM( \ + CCTK_POINTER_TO_CONST IN cctkGH, \ + CCTK_INT IN size, \ + CCTK_POINTER IN source, \ + CCTK_POINTER IN x, \ + CCTK_POINTER IN y, \ + CCTK_POINTER IN z \ + ) +USES FUNCTION Rho_ADM diff --git a/param.ccl b/param.ccl index 3703c97..82463f6 100644 --- a/param.ccl +++ b/param.ccl @@ -101,3 +101,10 @@ REAL par_S_minus[3] "spin of the m- puncture" { (*:*) :: "" } 0.0 + + +INT par_use_sources "Use sources?" +{ + 0:1 :: "0 for no (default), 1 for yes" +} 0 + diff --git a/schedule.ccl b/schedule.ccl index 797d4d6..37262fb 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -5,7 +5,10 @@ if (keep_u_around) { STORAGE: puncture_u } - +SCHEDULE TwoPunctures_ParamCheck AT PARAMCHECK +{ + LANG: C +} "Check parameters and thorn needs" SCHEDULE TwoPunctures IN ADMBase_InitialData { diff --git a/src/Equations.c b/src/Equations.c index e1de7b5..0ea280a 100644 --- a/src/Equations.c +++ b/src/Equations.c @@ -136,7 +136,8 @@ BY_Aijofxyz (double x, double y, double z, double Aij[3][3]) //******* Nonlinear Equations ********** //--------------------------------------------------------------- void -NonLinEquations (double A, double B, double X, double R, +NonLinEquations (CCTK_REAL rho_adm, + double A, double B, double X, double R, double x, double r, double phi, double y, double z, derivs U, double *values) { @@ -154,7 +155,8 @@ NonLinEquations (double A, double B, double X, double R, psi7 = psi * psi2 * psi4; values[0] = - U.d11[0] + U.d22[0] + U.d33[0] + 0.125 * BY_KKofxyz (x, y, z) / psi7; + U.d11[0] + U.d22[0] + U.d33[0] + 0.125 * BY_KKofxyz (x, y, z) / psi7 + + 2.0 * Pi * psi4*psi * rho_adm; } diff --git a/src/FuncAndJacobian.c b/src/FuncAndJacobian.c index 7c30da0..c2c972e 100644 --- a/src/FuncAndJacobian.c +++ b/src/FuncAndJacobian.c @@ -7,6 +7,7 @@ #include #include #include +#include "cctk.h" #include "cctk_Parameters.h" #include "TP_utilities.h" #include "TwoPunctures.h" @@ -188,19 +189,78 @@ 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) +F_of_v (CCTK_POINTER_TO_CONST cctkGH, + 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 "+/-" + DECLARE_CCTK_PARAMETERS; int i, j, k, ivar, indx; double al, be, A, B, X, R, x, r, phi, y, z, Am1, *values; derivs U; + CCTK_REAL *sources; values = dvector (0, nvar - 1); allocate_derivs (&U, nvar); + if (par_use_sources) + { + CCTK_REAL *s_x, *s_y, *s_z; + CCTK_INT i3D; + s_x =calloc(n1*n2*n3, sizeof(CCTK_REAL)); + s_y =calloc(n1*n2*n3, sizeof(CCTK_REAL)); + s_z =calloc(n1*n2*n3, sizeof(CCTK_REAL)); + sources=calloc(n1*n2*n3, sizeof(CCTK_REAL)); + for (i = 0; i < n1; i++) + for (j = 0; j < n2; j++) + for (k = 0; k < n3; k++) + { + i3D = i*n2*n3 + j*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, &(s_x[i3D]), &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, s_x[i3D], r, phi, &(s_y[i3D]), &(s_z[i3D]), U); + } + Rho_ADM(cctkGH, n1*n2*n3, sources, s_x, s_y, s_z); + free(s_z); + free(s_y); + free(s_x); + } + else + for (i = 0; i < n1; i++) + for (j = 0; j < n2; j++) + for (k = 0; k < n3; k++) + sources[i*n1*n2 + j*n2 + k]=0.0; + Derivatives_AB3 (nvar, n1, n2, n3, v); for (i = 0; i < n1; i++) @@ -210,55 +270,58 @@ F_of_v (int nvar, int n1, int n2, int n3, derivs v, double *F, 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 - } + 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 (sources[i*n2*n3 + j*n3 + k], + 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 + } } } } + if (par_use_sources) + free(sources); free_dvector (values, 0, nvar - 1); free_derivs (&U, nvar); } diff --git a/src/Newton.c b/src/Newton.c index 593789f..0dbb627 100644 --- a/src/Newton.c +++ b/src/Newton.c @@ -173,7 +173,8 @@ relax (double *dv, int nvar, int n1, int n2, int n3, // ----------------------------------------------------------------------------------- void -TestRelax (int nvar, int n1, int n2, int n3, derivs v, +TestRelax (CCTK_POINTER_TO_CONST cctkGH, + int nvar, int n1, int n2, int n3, derivs v, double *dv) { DECLARE_CCTK_PARAMETERS; @@ -190,7 +191,7 @@ TestRelax (int nvar, int n1, int n2, int n3, derivs v, cols = imatrix (0, ntotal - 1, 0, maxcol - 1); ncols = ivector (0, ntotal - 1); - F_of_v (nvar, n1, n2, n3, v, F, u); + F_of_v (cctkGH, nvar, n1, n2, n3, v, F, u); SetMatrix_JFD (nvar, n1, n2, n3, u, ncols, cols, JFD); @@ -225,8 +226,9 @@ TestRelax (int nvar, int n1, int n2, int n3, derivs v, // ----------------------------------------------------------------------------------- int -bicgstab (int nvar, int n1, int n2, int n3, derivs v, - derivs dv, int output, int itmax, double tol, double *normres) +bicgstab (CCTK_POINTER_TO_CONST cctkGH, + 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; @@ -246,7 +248,7 @@ bicgstab (int nvar, int n1, int n2, int n3, derivs v, cols = imatrix (0, ntotal - 1, 0, maxcol - 1); ncols = ivector (0, ntotal - 1); - F_of_v (nvar, n1, n2, n3, v, F, u); + F_of_v (cctkGH, nvar, n1, n2, n3, v, F, u); SetMatrix_JFD (nvar, n1, n2, n3, u, ncols, cols, JFD); /* temporary storage */ @@ -389,7 +391,8 @@ bicgstab (int nvar, int n1, int n2, int n3, derivs v, // ------------------------------------------------------------------- void -Newton (int nvar, int n1, int n2, int n3, +Newton (CCTK_POINTER_TO_CONST cctkGH, + int nvar, int n1, int n2, int n3, derivs v, double tol, int itmax) { DECLARE_CCTK_PARAMETERS; @@ -401,7 +404,7 @@ Newton (int nvar, int n1, int n2, int n3, allocate_derivs (&dv, ntotal); allocate_derivs (&u, ntotal); -/* TestRelax(nvar, n1, n2, n3, v, dv.d0); */ +/* TestRelax(nvar, n1, n2, n3, v, dv.d0); */ for (j = 0; j < ntotal; j++) v.d0[j] = 0; @@ -411,28 +414,28 @@ Newton (int nvar, int n1, int n2, int n3, { if (it == 0) { - F_of_v (nvar, n1, n2, n3, v, F, u); + F_of_v (cctkGH, 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; - } + if (fabs (F[j]) > dmax) + dmax = fabs (F[j]); } + for (j = 0; j < ntotal; j++) + dv.d0[j] = 0; printf ("Newton: it=%d \t |F|=%e\n", it, dmax); fflush(stdout); ii = - bicgstab (nvar, n1, n2, n3, v, dv, verbose, 100, dmax * 1.e-3, &normres); + bicgstab (cctkGH, + nvar, n1, n2, n3, v, dv, verbose, 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); + F_of_v (cctkGH, nvar, n1, n2, n3, v, F, u); dmax = -1; for (j = 0; j < ntotal; j++) { if (fabs (F[j]) > dmax) { - dmax = fabs (F[j]); + dmax = fabs (F[j]); } } it += 1; diff --git a/src/ParamCheck.c b/src/ParamCheck.c new file mode 100644 index 0000000..ca06f64 --- /dev/null +++ b/src/ParamCheck.c @@ -0,0 +1,32 @@ +// TwoPunctures: File "TwoPunctures.c" + +#include +#include +#include +#include +#include +#include +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" +#include "TP_utilities.h" +#include "TwoPunctures.h" + +// ------------------------------------------------------------------- +void +TwoPunctures_ParamCheck (CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + if (par_use_sources) + { + CCTK_INFO("Solving for BH-NS"); + if (CCTK_IsFunctionAliased ("Rho_ADM")) + CCTK_INFO("Aliased Functions found"); + else + CCTK_WARN(0, "I found no (aliased) function for matter sources, but " + "was said to use matter.\n"); + } + else + CCTK_INFO("not using sources (only BHs)"); +} diff --git a/src/TwoPunctures.c b/src/TwoPunctures.c index 1c792b6..e926b65 100644 --- a/src/TwoPunctures.c +++ b/src/TwoPunctures.c @@ -46,9 +46,9 @@ TwoPunctures (CCTK_ARGUMENTS) allocate_derivs (&v, ntotal); CCTK_INFO ("Solving puncture equation"); - Newton (nvar, n1, n2, n3, v, Newton_tol, Newton_maxit); + Newton (cctkGH, nvar, n1, n2, n3, v, Newton_tol, Newton_maxit); - F_of_v (nvar, n1, n2, n3, v, F, u); + F_of_v (cctkGH, nvar, n1, n2, n3, v, F, u); /* print out ADM mass, eq.: \Delta M_ADM=2*r*u=4*b*V for A=1,B=0,phi=0 */ admMass = (par_m_plus + par_m_minus diff --git a/src/TwoPunctures.h b/src/TwoPunctures.h index 9915181..19763e7 100644 --- a/src/TwoPunctures.h +++ b/src/TwoPunctures.h @@ -30,7 +30,8 @@ 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, +void F_of_v (CCTK_POINTER_TO_CONST cctkGH, + 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); @@ -61,7 +62,8 @@ void rx3_To_xyz (int nvar, double x, double r, double phi, double *y, // 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, +void NonLinEquations (CCTK_REAL rho_adm, + 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, @@ -77,12 +79,14 @@ 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, +void TestRelax (CCTK_POINTER_TO_CONST cctkGH, + int nvar, int n1, int n2, int n3, derivs v, double *dv); +int bicgstab (CCTK_POINTER_TO_CONST cctkGH, + 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); +void Newton (CCTK_POINTER_TO_CONST cctkGH, + int nvar, int n1, int n2, int n3, derivs v, + double tol, int itmax); /* diff --git a/src/make.code.defn b/src/make.code.defn index 726aec8..512024e 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -2,7 +2,7 @@ # $Header$ # Source files in this directory -SRCS = CoordTransf.c Equations.c FuncAndJacobian.c Newton.c TwoPunctures.c TP_utilities.c +SRCS = CoordTransf.c Equations.c FuncAndJacobian.c Newton.c TwoPunctures.c TP_utilities.c ParamCheck.c # Subdirectories containing source files SUBDIRS = -- cgit v1.2.3