diff options
-rw-r--r-- | src/gr/AHFinderDirect.hh | 17 | ||||
-rw-r--r-- | src/gr/Jacobian.hh | 30 | ||||
-rw-r--r-- | src/gr/auxiliary.maple | 6 | ||||
-rw-r--r-- | src/gr/cg.hh | 28 | ||||
-rw-r--r-- | src/gr/doit.maple | 4 | ||||
-rw-r--r-- | src/gr/driver.cc | 2 | ||||
-rw-r--r-- | src/gr/gfn.hh | 7 | ||||
-rw-r--r-- | src/gr/horizon.maple | 158 | ||||
-rw-r--r-- | src/gr/horizon_Jacobian.cc | 161 | ||||
-rw-r--r-- | src/gr/horizon_function.cc | 58 | ||||
-rw-r--r-- | src/gr/maple.log | 64 | ||||
-rw-r--r-- | src/gr/setup_gr_gfas.maple | 38 |
12 files changed, 496 insertions, 77 deletions
diff --git a/src/gr/AHFinderDirect.hh b/src/gr/AHFinderDirect.hh index 59dacb8..8bd0225 100644 --- a/src/gr/AHFinderDirect.hh +++ b/src/gr/AHFinderDirect.hh @@ -21,10 +21,11 @@ struct geometry_interpolator_info // // This structure holds all the information we need about the Cactus grid -// and gridfns outside the top-level driver. At present we use -// CCTK_InterpLocalUniform() -// to access the Cactus gridfns, and much of this information is specific -// to that API. +// and gridfns outside the top-level driver. At present we use the +// CCTK_InterpLocalUniform() local interpolator to access the Cactus +// gridfns. Much of the information in this structure is specific to +// that API, and (alas) will need changing when we eventually switch to +// a "global" multiprocessor/distributed interpolator. // struct cactus_grid_info { @@ -67,4 +68,10 @@ extern "C" // horizon_function.cc void horizon_function(patch_system& ps, const struct cactus_grid_info& cgi, - const struct geometry_interpolator_info& gii); + const struct geometry_interpolator_info& gii, + bool Jacobian_flag); + +// horizon_Jacobian.cc +class Jacobian; +void horizon_Jacobian(const patch_system& ps, + Jacobian& Jac); diff --git a/src/gr/Jacobian.hh b/src/gr/Jacobian.hh new file mode 100644 index 0000000..626d080 --- /dev/null +++ b/src/gr/Jacobian.hh @@ -0,0 +1,30 @@ +// Jacobian.hh -- data structures for the Jacobian matrix +// $Id$ + +// +// prerequisites: +// + +// +// A Jacobian object stores the (a) Jacobian matrix of H with respect to h. +// Jacobian is an abstract base class, from which we derive a separate +// concrete class for each type of (storage scheme for the) Jacobian matrix. +// +// A row/column index of the Jacobian (denoted II/JJ) is a linearization +// of the patch system's (patch,irho,isigma). +// + +class Jacobian + { +public: + // access a matrix element + virtual const fp& operator()(int II, int JJ) const; // rvalue + virtual fp& operator()(int II, int JJ) const; // lvalue + + // zero the whole matrix or a single row + virtual void zero(); + virtual void zero_row(int II); + + // add to an element: Jacobian(II,JJ) += delta + virtual void add_to_element(int II, int JJ, fp delta); + }; diff --git a/src/gr/auxiliary.maple b/src/gr/auxiliary.maple index fecfc48..99069e7 100644 --- a/src/gr/auxiliary.maple +++ b/src/gr/auxiliary.maple @@ -14,7 +14,7 @@ # algebraic functions of the basic ADM fields {g_dd,K_dd}, and optionally # also generates C code for them (the auxiliar variables). # -# Code files are written to the directory cg/ . +# Code files are written to the directory ../gr.cg/ . # # Inputs (Maple): # N @@ -73,7 +73,7 @@ assert_fnd_exists(g_uu, fnd); g_uu__fnd := linalg[inverse](g_dd); if (cg_flag) - then codegen2(g_uu__fnd, 'g_uu', "cg/inverse_metric.c"); + then codegen2(g_uu__fnd, 'g_uu', "../gr.cg/inverse_metric.c"); fi; NULL; @@ -132,7 +132,7 @@ K__fnd := simplify( if (cg_flag) then codegen2([K__fnd, K_uu__fnd], ['K', 'K_uu'], - "cg/extrinsic_curvature_trace_raise.c"); + "../gr.cg/extrinsic_curvature_trace_raise.c"); fi; NULL; diff --git a/src/gr/cg.hh b/src/gr/cg.hh index b6174eb..3095a5a 100644 --- a/src/gr/cg.hh +++ b/src/gr/cg.hh @@ -9,6 +9,8 @@ // gfn.hh // +//****************************************************************************** + // // The machine-generated code uses the following variables: // patch& p // current patch @@ -37,6 +39,8 @@ #define PARTIAL_SIGMA_SIGMA(ghosted_gridfn_name) \ p.partial_sigma_sigma(ghosted_gfns::gfn__ ## ghosted_gridfn_name, irho,isigma) +//****************************************************************************** + // // ghosted-grid gridfns // n.b. since we're always evaluating on the surface, r == h @@ -44,6 +48,8 @@ #define h p.ghosted_gridfn(ghosted_gfns::gfn__h, irho,isigma) #define r h +//****************************************************************************** + // // nominal-grid gridfns // @@ -81,9 +87,22 @@ #define H p.gridfn(nominal_gfns::gfn__H, irho,isigma) +#define partial_H_wrt_partial_d_h_1 \ + p.gridfn(nominal_gfns::gfn__partial_H_wrt_partial_d_h_1, irho,isigma) +#define partial_H_wrt_partial_d_h_2 \ + p.gridfn(nominal_gfns::gfn__partial_H_wrt_partial_d_h_2, irho,isigma) +#define partial_H_wrt_partial_dd_h_11 \ + p.gridfn(nominal_gfns::gfn__partial_H_wrt_partial_dd_h_11, irho,isigma) +#define partial_H_wrt_partial_dd_h_12 \ + p.gridfn(nominal_gfns::gfn__partial_H_wrt_partial_dd_h_12, irho,isigma) +#define partial_H_wrt_partial_dd_h_22 \ + p.gridfn(nominal_gfns::gfn__partial_H_wrt_partial_dd_h_22, irho,isigma) + +//****************************************************************************** + // -// pseudo-gridfns -// (used only at a single grid point, not actually stored as gridfns) +// pseudo-gridfns, i.e. temporaries used only at a single grid point +// (these are not actually stored as gridfns) // fp g_uu_11; fp g_uu_12; @@ -121,3 +140,8 @@ fp partial_d_g_uu_313; fp partial_d_g_uu_322; fp partial_d_g_uu_323; fp partial_d_g_uu_333; + +fp HA; +fp HB; +fp HC; +fp HD; diff --git a/src/gr/doit.maple b/src/gr/doit.maple index 8d5edf1..86cdebe 100644 --- a/src/gr/doit.maple +++ b/src/gr/doit.maple @@ -6,10 +6,10 @@ read "setup_gr_gfas.mm"; setup_gr_gfas(); read "auxiliary.mm"; -read "metric_derivs.mm"; +read "curvature.mm"; read "horizon.mm"; `saveit/level` := 10; auxiliary(true); -metric_derivs(true); +curvature(true); horizon(true); diff --git a/src/gr/driver.cc b/src/gr/driver.cc index fe64dce..e95c783 100644 --- a/src/gr/driver.cc +++ b/src/gr/driver.cc @@ -200,7 +200,7 @@ else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, // if (STRING_EQUAL(method, "horizon")) then { - horizon_function(ps, cgi, gii); + horizon_function(ps, cgi, gii, false); ps.print_gridfn_with_xyz(nominal_gfns::gfn__H, true, ghosted_gfns::gfn__h, "H.dat"); diff --git a/src/gr/gfn.hh b/src/gr/gfn.hh index 16af52f..24c56fe 100644 --- a/src/gr/gfn.hh +++ b/src/gr/gfn.hh @@ -48,6 +48,11 @@ enum { gfn__partial_d_g_dd_323, gfn__partial_d_g_dd_333, gfn__H, - max_gfn = gfn__H + gfn__partial_H_wrt_partial_d_h_1, + gfn__partial_H_wrt_partial_d_h_2, + gfn__partial_H_wrt_partial_dd_h_11, + gfn__partial_H_wrt_partial_dd_h_12, + gfn__partial_H_wrt_partial_dd_h_22, + max_gfn = gfn__partial_H_wrt_partial_dd_h_22 // no comma }; }; diff --git a/src/gr/horizon.maple b/src/gr/horizon.maple index 6950e56..0707af4 100644 --- a/src/gr/horizon.maple +++ b/src/gr/horizon.maple @@ -4,7 +4,8 @@ # horizon - overall driver ## non_unit_normal - compute s_d ## non_unit_normal_deriv - compute partial_d_s_d -## horizon_function - compute HA,HB,HC,HD,H = fn(s_d) +## horizon_function - compute H[ABCD],H = fn(s_d__fnd) = fn(h) +## horizon_Jacobian - compute Jacobian coeffs for H[ABCD], H # ################################################################################ @@ -27,6 +28,7 @@ proc(cg_flag::boolean) non_unit_normal(); non_unit_normal_deriv(); horizon_function(cg_flag); +horiozn_Jacobian(cg_flag); NULL; end proc; @@ -34,9 +36,14 @@ end proc; ################################################################################ # -# This function computes the non-unit outward-pointing normal (covector) -# to the horizon, s_d = n_u(h), using the equation on p.7.1 of my apparent -# horizon finding notes. This function does *NOT* generate any C code. +# This function computes the xyz components of the non-unit +# outward-pointing normal (covector) to the horizon, s_d = s_d(h), +# using the equation on p.7.1 of my apparent horizon finding notes. +# +# This function does *NOT* generate any C code. +# +# Outputs: +# s_d = s_d__fnd( x_xyz, X_ud, Diff(h,y_rs) ) # non_unit_normal := proc() @@ -64,10 +71,18 @@ end proc; ################################################################################ # -# This function computes the 1st derivatives of the non-unit outward-pointing +# This function computes the xyz derivatives of the non-unit outward-pointing # normal (covector) to the horizon, using the equation on p.7.2 of my -# apparent horizon finding notes. This function does *NOT* generate any -# C code. +# apparent horizon finding notes. +# +# This function does *NOT* generate any C code. +# +# Inputs: +# s_d = s_d__fnd( x_xyz, X_ud, Diff(h,y_rs) ) +# +# Outputs: +# partial_d_s_d = partial_d_s_d__fnd( x_xyz, X_ud, X_udd, +# Diff(h,y_rs), Diff(h,y_rs,y_rs) ) # non_unit_normal_deriv := proc() @@ -113,9 +128,32 @@ end proc; # as a function of 1st and 2nd angular partial derivatives of the # apparent horizon radius $h$, in the Maple gridfn array H__fnd, # using equation (15) in my 1996 apparent horizon finding paper, -# and using partial_d_g_uu[k,i,j] instead of Diff(g_uu[i,j], x_xyz[k]). +# but using partial_d_g_uu[k,i,j] instead of Diff(g_uu[i,j], x_xyz[k]). +# +# These equation gives H[ABCD] and H as functions of s_d, but here we +# use s_d__fnd, which we assume is already given in terms of angular +# derivatives of h. The result is that we compute H[ABCD] and H directly +# in terms of 1st and 2nd angular derivatives of h, without s_d ever +# appearing in our final results. +# # This function also optionally generates C code for this computation. # +# Inputs: +# s_d = s_d__fnd( x_xyz, X_ud, Diff(h,y_rs) ) +# partial_d_s_d = partial_d_s_d__fnd( x_xyz, X_ud, X_udd, +# Diff(h,y_rs), Diff(h,y_rs,y_rs) ) +# +# Outputs: +# HA = HA__fnd( x_xyz, X_ud, X_udd, g_uu, partial_d_g_uu, +# Diff(h,y_rs), Diff(h,y_rs,y_rs) ) +# HB = HB__fnd( x_xyz, X_ud, X_udd, +# g_uu, partial_d_g_uu, partial_d_ln_sqrt_g, +# Diff(h,y_rs), Diff(h,y_rs,y_rs) ) +# HC = HC__fnd( x_xyz, X_ud, K_uu, Diff(h,y_rs) ) +# HC = HC__fnd( x_xyz, X_ud, g_uu, Diff(h,y_rs) ) +# --------------------------------------------------------------- +# H = H__fnd( HA__fnd, HB__fnd, HC__fnd, HD__fnd ) +# horizon_function := proc(cg_flag::boolean) global @@ -163,7 +201,109 @@ HD__fnd := msum('g_uu[i,j]*s_d__fnd[i]*s_d__fnd[j]', 'i'=1..N, 'j'=1..N); H__fnd := HA__fnd/HD__fnd^(3/2) + HB__fnd/HD__fnd^(1/2) + HC__fnd/HD__fnd - K; if (cg_flag) - then codegen2(H__fnd, 'H', "cg/horizon_function.c"); + then codegen2([ HA__fnd, HB__fnd, HC__fnd, HD__fnd, H__fnd], + ['HA', 'HB', 'HC', 'HD', 'H' ], + "../gr.cg/horizon_function.c"); +fi; + +NULL; +end proc; + +################################################################################ + +# +# This function computes the Jacobian coefficients for the LHS of the +# apparent horizon equation with respect to 1st and 2nd angular partial +# derivatives of the apparent horizon radius $h$. These coefficients +# are (or at least should be :) the same as those in equation (A1) in +# my 1996 apparent horizon finding paper, bue here we compute them in +# a different manner: we have Maple directly differentiate H__fnd with +# respect to Diff(h,y_rs[u]) and Diff(h,y_rs[u],y_rs[v]). We use the +# Maple frontend() function to do this. +# +# This function also optionally generates C code for the Jacobian +# coefficients. +# +# Inputs: +# HA = HA__fnd( x_xyz, X_ud, X_udd, g_uu, partial_d_g_uu, +# Diff(h,y_rs), Diff(h,y_rs,y_rs) ) +# HB = HB__fnd( x_xyz, X_ud, X_udd, +# g_uu, partial_d_g_uu, partial_d_ln_sqrt_g, +# Diff(h,y_rs), Diff(h,y_rs,y_rs) ) +# HC = HC__fnd( x_xyz, X_ud, K_uu, Diff(h,y_rs) ) +# HC = HC__fnd( x_xyz, X_ud, g_uu, Diff(h,y_rs) ) +# H = H__fnd( HA__fnd, HB__fnd, HC__fnd, HD__fnd ) +# +# Outputs: +# partial_HA_wrt_partial_d_h partial_HA_wrt_partial_dd_h +# partial_HB_wrt_partial_d_h partial_HB_wrt_partial_dd_h +# partial_HC_wrt_partial_d_h +# partial_HD_wrt_partial_d_h +# --------------------------------------------------------------- +# partial_H_wrt_partial_d_h partial_H_wrt_partial_dd_h +# +horizon_Jacobian := +proc(cg_flag::boolean) +global + @include "../maple/coords.minc", + @include "../maple/gfa.minc", + @include "../gr/gr_gfas.minc"; +local u,v, + temp; + +printf("%a...\n", procname); + +assert_fnd_exists(g_uu); +assert_fnd_exists(K_uu); +assert_fnd_exists(partial_d_ln_sqrt_g); +assert_fnd_exists(partial_d_g_uu); +assert_fnd_exists(HA) +assert_fnd_exists(HB) +assert_fnd_exists(HC) +assert_fnd_exists(HD) + +# Jacobian coefficients of H[ABCD] and H wrt Diff(h,y_rs[u]) + for u from 1 to N_ang + do + partial_HA_wrt_partial_d_h__fnd[u] + := frontend('diff', [HA__fnd, Diff(h,y_rs[u])]); + partial_HB_wrt_partial_d_h__fnd[u] + := frontend('diff', [HB__fnd, Diff(h,y_rs[u])]); + partial_HC_wrt_partial_d_h__fnd[u] + := frontend('diff', [HC__fnd, Diff(h,y_rs[u])]); + partial_HD_wrt_partial_d_h__fnd[u] + := frontend('diff', [HD__fnd, Diff(h,y_rs[u])]); + + # equation (A1a) in my 1996 apparent horizon finding paper + temp := (3/2)*HA/HD^(5/2) + (1/2)*HB/HD^(3/2) + HC/HD^2; + partial_H_wrt_partial_d_h__fnd[u] + := + partial_HA_wrt_partial_d_h__fnd[u] / HD^(3/2) + + partial_HB_wrt_partial_d_h__fnd[u] / HD^(1/2) + + partial_HC_wrt_partial_d_h__fnd[u] / HD + - partial_HD_wrt_partial_d_h__fnd[u] * temp; + end do; + +# Jacobian coefficients of H[AB] and H wrt Diff(h,y_rs[u],y_rs[v]) + for u from 1 to N_ang + do + for v from u to N_ang + do + partial_HA_wrt_partial_dd_h__fnd[u,v] + := frontend('diff', [HA__fnd, Diff(h,y_rs[u],y_rs[v])]); + partial_HB_wrt_partial_dd_h__fnd[u,v] + := frontend('diff', [HB__fnd, Diff(h,y_rs[u],y_rs[v])]); + + # equation (A1b) in my 1996 apparent horizon finding paper + partial_H_wrt_partial_dd_h__fnd[u,v] + := + partial_HA_wrt_partial_dd_h__fnd[u,v] / HD^(3/2) + + partial_HB_wrt_partial_dd_h__fnd[u,v] / HD^(1/2) + end do; + +if (cg_flag) + then codegen2([partial_H_wrt_partial_d_h__fnd, + partial_H_wrt_partial_dd_h__fnd], + ['partial_H_wrt_partial_d_h', 'partial_H_wrt_partial_dd_h'], + "../gr.cg/horizon_Jacobian.c"); fi; NULL; diff --git a/src/gr/horizon_Jacobian.cc b/src/gr/horizon_Jacobian.cc new file mode 100644 index 0000000..40fada6 --- /dev/null +++ b/src/gr/horizon_Jacobian.cc @@ -0,0 +1,161 @@ +// horizon_Jacobian.cc -- evaluate Jacobian matrix of LHS function H(h) +// $Id$ +// +// <<<prototypes for functions local to this file>>> +// horizon_Jacobian - top-level driver +// + +#include <stdio.h> +#include <assert.h> +#include <math.h> +#include <vector> + +#include "util_Table.h" +#include "cctk.h" +#include "cctk_Arguments.h" + +#include "jt/stdc.h" +#include "jt/util.hh" +#include "jt/array.hh" +#include "jt/cpm_map.hh" +#include "jt/linear_map.hh" +using jtutil::error_exit; + +#include "../config.hh" + +#include "../util/coords.hh" +#include "../util/grid.hh" +#include "../util/fd_grid.hh" +#include "../util/patch.hh" +#include "../util/patch_edge.hh" +#include "../util/patch_interp.hh" +#include "../util/ghost_zone.hh" +#include "../util/patch_system.hh" + +#include "gfn.hh" +#include "Jacobian.hh" +#include "AHFinderDirect.hh" + +//****************************************************************************** + +// +// ***** prototypes for functions local to this file ***** +// + +namespace { + }; + +//****************************************************************************** +//****************************************************************************** +//****************************************************************************** + +// +// This function computes the Jacobian matrix of the LHS function H(h) +// from the Jacobian coefficient (angular) gridfns. The computation is +// done a Jacobian row at a time, using equation (25) of my 1996 apparent +// horizon finding paper. +// +// Inputs (angular gridfns, on ghosted grid): +// h # shape of trial surface +// partial_H_wrt_partial_d_h, # Jacobian coefficients +// partial_H_wrt_partial_dd_h, +// +// Outputs: +// The Jacobian matrix is stored in the Jacobian object Jac. +void horizon_Jacobian(const patch_system& ps, + Jacobian& Jac); +{ +CCTK_VInfo(CCTK_THORNSTRING, " horizon Jacobian"); + +Jac.zero(II); + + for (int xpn = 0 ; xpn < ps.N_patches() ; ++xpn) + { + patch& xp = ps.ith_patch(xpn); + + for (int x_irho = xp.min_irho() ; x_irho <= xp.max_irho() ; ++x_irho) + { + for (int x_isigma = xp.min_isigma() ; + x_isigma <= xp.max_isigma() ; + ++x_isigma) + { + // + // compute the Jacobian row for this grid point, i.e. + // partial H(this point x, Jacobian row II) + // --------------------------------------------- + // partial h(other points y, Jacobian column JJ) + // + // FIXME FIXME: we still have to take into account the + // position-dependence of the coefficients, + // cf the difference between J[3H(h)] and J[2H(h)]; + // here we're sort of computing the former, but + // Newton's method really wants the latter + // + + // Jacobian row index for this point + const int II = ps.gpn_of_patch_irho_isigma(xp,x_irho,x_isigma); + + // Jacobian coefficients for this point + const fp Jacobian_coeff_rho + = xp.gridfn(gfn__partial_H_wrt_partial_d_h_1, x_irho,x_isigma); + const fp Jacobian_coeff_sigma + = xp.gridfn(gfn__partial_H_wrt_partial_d_h_2, x_irho,x_isigma); + const fp Jacobian_coeff_rho_rho + = xp.gridfn(gfn__partial_H_wrt_partial_dd_h_11, x_irho,x_isigma); + const fp Jacobian_coeff_rho_sigma + = xp.gridfn(gfn__partial_H_wrt_partial_dd_h_12, x_irho,x_isigma); + const fp Jacobian_coeff_sigma_sigma + = xp.gridfn(gfn__partial_H_wrt_partial_dd_h_22, x_irho,x_isigma); + + // partial_rho, partial_rho_rho + for (int m_irho = xp.molecule_min_m() ; + m_irho <= xp.molecule_max_m() ; + ++m_irho) + { + const fp Jac_rho = Jacobian_coeff_rho * xp.partial_rho_coeff(m_rho); + const fp Jac_rho_rho = Jacobian_coeff_rho_rho + * xp.partial_rho_rho_coeff(m_rho); + const fp Jac_sum = Jac_rho + Jac_rho_rho; + if (xp.is_in_nominal_grid(x_irho+m_irho, x_isigma)) + then { + const int JJ + = ps.gpn_of_patch_irho_isigma(xp, x_irho+m_irho, + x_isigma); + Jac.add_to_element(II,JJ, Jac_sum); + } + else { + const patch_edge& xe = (m_irho < 0) + ? xp.min_rho_patch_edge() + : xp.max_rho_patch_edge() + const ghost_zone& gz = xp.ghost_zone_on_edge(xe); + const int x_iperp = xe.iperp_of_irho_isigma(x_irho, + x_isigma); + const int x_ipar = xe.ipar_of_irho_isigma(x_irho, x_isigma); + gz.compute_Jacobian(gfn__h, gfn__h); + const patch& yp = gz.Jacobian_y_patch(); + const patch_edge& ye = gz.Jacobian_y_edge(); + const int min_y_ipar_m = gz.Jacobian_min_y_ipar_m(); + const int max_y_ipar_m = gz.Jacobian_max_y_ipar_m(); + const int y_iperp = gz.Jacobian_y_iperp(x_iperp); + const int y_ipar_posn + = gz.Jacobian_y_ipar_posn(x_iperp, x_ipar); + for (int y_ipar_m = min_y_ipar_m ; + y_ipar_m <= max_y_ipar_m ; + ++y_ipar_m) + { + const int y_ipar = y_ipar_posn + y_ipar_m; + const int y_irho + = ye.irho_of_iperp_ipar(y_iperp, y_ipar); + const int y_isigma + = ye.isigma_of_iperp_ipar(y_iperp, y_ipar); + const int JJ + = ps.gpn_of_patch_irho_isigma(yp, y_irho, y_isigma); + const fp gz_Jac = gz.Jacobian(x_iperp,x_ipar, y_ipar_m); + Jac.add_to_element(II,JJ, Jac_sum*gz_Jac); + } + } + } + } + } + } +} diff --git a/src/gr/horizon_function.cc b/src/gr/horizon_function.cc index 6ea0a8d..745760b 100644 --- a/src/gr/horizon_function.cc +++ b/src/gr/horizon_function.cc @@ -57,7 +57,9 @@ void compute_H(patch_system& ps); //****************************************************************************** // -// This function computes the LHS function H(h). +// This function computes the LHS function H(h), and optionally also +// its Jacobian coefficients (from which the Jacobian matrix may be +// computed later). // // Inputs (angular gridfns, on ghosted grid): // ... defined on ghosted grid @@ -86,9 +88,14 @@ void compute_H(patch_system& ps); // partial_d_g_dd_{1,2,3}{11,12,13,22,23,33} # $\partial_k g_{ij}$ // H # $H = H(h)$ // +// Arguments: +// Jacobian_flag = true to compute the Jacobian coefficients, +// false to skip this. +// void horizon_function(patch_system& ps, const struct cactus_grid_info& cgi, - const struct geometry_interpolator_info& ii) + const struct geometry_interpolator_info& ii, + bool Jacobian_flag) { CCTK_VInfo(CCTK_THORNSTRING, " horizon function"); @@ -101,9 +108,9 @@ setup_xyz_posns(ps); // interpolate $g_{ij}$, $K_{ij}$ --> $g_{ij}$, $K_{ij}$, $\partial_k g_{ij}$ interpolate_geometry(ps, cgi, ii); -// compute remaining gridfns --> $H$ +// compute remaining gridfns --> $H$ and optionally Jacobian coefficients // by algebraic ops and angular finite differencing -compute_H(ps); +compute_H(ps, Jacobian_flag); } //****************************************************************************** @@ -302,8 +309,9 @@ void* const output_arrays[N_OUTPUT_ARRAYS] bool first_time = true; if (first_time) then { - // store derivative info in interpolator parameter table first_time = false; + + // store derivative info in interpolator parameter table CCTK_VInfo(CCTK_THORNSTRING, " setting up derivative info for interpolator"); @@ -404,11 +412,17 @@ if (status < 0) //****************************************************************************** // -// This function computes H(h) by a mixture of algebraic operations -// and (rho,sigma) finite differencing, on the angular grid. +// This function computes H(h), and optionally its Jacobian coefficients, +// (from which the Jacobian matrix may be computed later). This function +// uses a mixture of algebraic operations and (rho,sigma) finite differencing. +// The computation is done (entirely) on the nominal angular grid. +// +// Arguments: +// Jacobian_flag = true to compute the Jacobian coefficients, +// false to skip this. // namespace { -void compute_H(patch_system& ps) +void compute_H(patch_system& ps, bool Jacobian_flag) { CCTK_VInfo(CCTK_THORNSTRING, " computing H(h)"); @@ -460,22 +474,40 @@ CCTK_VInfo(CCTK_THORNSTRING, " computing H(h)"); // ... each cg/*.c file has a separate set of temp variables, // and so must be inside its own set of { } braces // + + // gridfn #defins #include "cg.hh" + { - #include "cg/inverse_metric.c" + // g_uu + #include "../gr.cg/inverse_metric.c" } + { - #include "cg/extrinsic_curvature_trace_raise.c" + // K, K_uu + #include "../gr.cg/extrinsic_curvature_trace_raise.c" } + { - #include "cg/inverse_metric_gradient.c" + // partial_d_g_uu + #include "../gr.cg/inverse_metric_gradient.c" } + { - #include "cg/metric_det_gradient.c" + // partial_d_ln_sqrt_g + #include "../gr.cg/metric_det_gradient.c" } + { - #include "cg/horizon_function.c" + // HA, HB, HC, HD, H + #include "../gr.cg/horizon_function.c" } + + if (Jacobian_flag) + then { + // partial_H_wrt_partial_d_h, partial_H_wrt_partial_dd_h + #include "../gr.cg/horizon_Jacobian.c" + } } } } diff --git a/src/gr/maple.log b/src/gr/maple.log index 7732b11..2e5ddd0 100644 --- a/src/gr/maple.log +++ b/src/gr/maple.log @@ -634,7 +634,8 @@ HA, HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd; assert_fnd_exists(g_dd); assert_fnd_exists(g_uu, fnd); g_uu__fnd := linalg[inverse](g_dd); - if cg_flag then codegen2(g_uu__fnd, 'g_uu', "cg/inverse_metric.c") + if cg_flag then + codegen2(g_uu__fnd, 'g_uu', "../gr.cg/inverse_metric.c") end if; NULL end proc @@ -661,7 +662,7 @@ HA, HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd; end do end do; if cg_flag then codegen2([K__fnd, K_uu__fnd], ['K', 'K_uu'], - "cg/extrinsic_curvature_trace_raise.c") + "../gr.cg/extrinsic_curvature_trace_raise.c") end if; NULL end proc @@ -692,7 +693,7 @@ HA, HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd; end do end do; if cg_flag then codegen2(partial_d_g_uu__fnd, 'partial_d_g_uu', - "cg/inverse_metric_gradient.c") + "../gr.cg/inverse_metric_gradient.c") end if; NULL end proc @@ -714,7 +715,7 @@ HA, HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd; 'g_uu[i, j]*partial_d_g_dd[k, i, j]', 'i' = 1 .. N, 'j' = 1 .. N)) end do; if cg_flag then codegen2(partial_d_ln_sqrt_g__fnd, - 'partial_d_ln_sqrt_g', "cg/metric_det_gradient.c") + 'partial_d_ln_sqrt_g', "../gr.cg/metric_det_gradient.c") end if; NULL end proc @@ -807,7 +808,8 @@ HA, HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd; H__fnd := HA__fnd/HD__fnd^(3/2) + HB__fnd/HD__fnd^(1/2) + HC__fnd/HD__fnd - K ; - if cg_flag then codegen2(H__fnd, 'H', "cg/horizon_function.c") end if; + if cg_flag then codegen2(H__fnd, 'H', "../gr.cg/horizon_function.c") + end if; NULL end proc @@ -817,7 +819,7 @@ end proc > auxiliary(true); inverse_metric... -codegen2(g_uu) --> "cg/inverse_metric.c" +codegen2(g_uu) --> "../gr.cg/inverse_metric.c" --> `codegen2/input` convert --> equation list --> `codegen2/eqnlist` @@ -826,7 +828,7 @@ codegen2(g_uu) --> "cg/inverse_metric.c" find temporary variables --> `codegen2/temps` convert Diff(expr,rho,sigma) --> PARTIAL_RHO_SIGMA(expr) etc -bytes used=1000208, alloc=917336, time=0.10 +bytes used=1000148, alloc=917336, time=0.12 --> `codegen2/fix_Diff` convert R_dd[2,3] --> R_dd_23 etc --> `codegen2/unindex` @@ -834,12 +836,12 @@ bytes used=1000208, alloc=917336, time=0.10 --> `codegen2/fix_rationals` writing C code extrinsic_curvature_trace_raise... -codegen2([K, K_uu]) --> "cg/extrinsic_curvature_trace_raise.c" +codegen2([K, K_uu]) --> "../gr.cg/extrinsic_curvature_trace_raise.c" --> `codegen2/input` convert --> equation list --> `codegen2/eqnlist` optimizing computation sequence -bytes used=2000708, alloc=1441528, time=0.15 +bytes used=2000796, alloc=1441528, time=0.22 --> `codegen2/optimize` find temporary variables --> `codegen2/temps` @@ -847,39 +849,39 @@ bytes used=2000708, alloc=1441528, time=0.15 --> `codegen2/fix_Diff` convert R_dd[2,3] --> R_dd_23 etc --> `codegen2/unindex` -bytes used=3000940, alloc=1572576, time=0.21 +bytes used=3001776, alloc=1572576, time=0.28 convert p/q --> RATIONAL(p/q) --> `codegen2/fix_rationals` writing C code > metric_derivs(true); inverse_metric_gradient... -bytes used=4001284, alloc=1638100, time=0.26 -codegen2(partial_d_g_uu) --> "cg/inverse_metric_gradient.c" +bytes used=4002196, alloc=1638100, time=0.36 +codegen2(partial_d_g_uu) --> "../gr.cg/inverse_metric_gradient.c" --> `codegen2/input` convert --> equation list --> `codegen2/eqnlist` optimizing computation sequence -bytes used=5002720, alloc=1638100, time=0.35 +bytes used=5002380, alloc=1638100, time=0.43 --> `codegen2/optimize` find temporary variables --> `codegen2/temps` convert Diff(expr,rho,sigma) --> PARTIAL_RHO_SIGMA(expr) etc --> `codegen2/fix_Diff` -bytes used=6002996, alloc=1703624, time=0.41 +bytes used=6004052, alloc=1703624, time=0.50 convert R_dd[2,3] --> R_dd_23 etc --> `codegen2/unindex` -bytes used=7003240, alloc=1703624, time=0.46 +bytes used=7004228, alloc=1703624, time=0.57 convert p/q --> RATIONAL(p/q) --> `codegen2/fix_rationals` writing C code -bytes used=8003492, alloc=1703624, time=0.49 +bytes used=8004548, alloc=1703624, time=0.63 metric_det_gradient... -codegen2(partial_d_ln_sqrt_g) --> "cg/metric_det_gradient.c" +codegen2(partial_d_ln_sqrt_g) --> "../gr.cg/metric_det_gradient.c" --> `codegen2/input` convert --> equation list --> `codegen2/eqnlist` optimizing computation sequence -bytes used=9003660, alloc=1703624, time=0.60 +bytes used=9004804, alloc=1703624, time=0.74 --> `codegen2/optimize` find temporary variables --> `codegen2/temps` @@ -894,40 +896,40 @@ codegen/C/expression: Unknown function: RATIONAL will be left as is. > horizon(true); non_unit_normal... non_unit_normal_deriv... -bytes used=10003824, alloc=1769148, time=0.66 +bytes used=10004956, alloc=1769148, time=0.85 horizon_function... -bytes used=11004072, alloc=1834672, time=0.73 -codegen2(H) --> "cg/horizon_function.c" +bytes used=11005112, alloc=1834672, time=0.93 +codegen2(H) --> "../gr.cg/horizon_function.c" --> `codegen2/input` convert --> equation list --> `codegen2/eqnlist` optimizing computation sequence -bytes used=12004376, alloc=1834672, time=0.79 -bytes used=13004892, alloc=2031244, time=0.87 -bytes used=14005160, alloc=2031244, time=0.96 +bytes used=12005688, alloc=1834672, time=1.01 +bytes used=13005864, alloc=2031244, time=1.12 +bytes used=14006148, alloc=2031244, time=1.21 --> `codegen2/optimize` find temporary variables --> `codegen2/temps` convert Diff(expr,rho,sigma) --> PARTIAL_RHO_SIGMA(expr) etc -bytes used=15005740, alloc=2096768, time=1.01 +bytes used=15006312, alloc=2096768, time=1.28 --> `codegen2/fix_Diff` convert R_dd[2,3] --> R_dd_23 etc -bytes used=16006176, alloc=2096768, time=1.07 +bytes used=16006476, alloc=2096768, time=1.35 --> `codegen2/unindex` -bytes used=17006396, alloc=2096768, time=1.13 +bytes used=17006632, alloc=2096768, time=1.41 convert p/q --> RATIONAL(p/q) -bytes used=18006808, alloc=2096768, time=1.18 +bytes used=18006900, alloc=2096768, time=1.48 --> `codegen2/fix_rationals` writing C code codegen/C/expression: Unknown function: PARTIAL_RHO will be left as is. codegen/C/expression: Unknown function: PARTIAL_SIGMA will be left as is. -bytes used=19007204, alloc=2096768, time=1.22 +bytes used=19007116, alloc=2096768, time=1.54 codegen/C/expression: Unknown function: PARTIAL_RHO_RHO will be left as is. codegen/C/expression: Unknown function: PARTIAL_RHO_SIGMA will be left as is. codegen/C/expression: Unknown function: PARTIAL_SIGMA_SIGMA will be left as is. -bytes used=20007392, alloc=2096768, time=1.38 +bytes used=20007280, alloc=2096768, time=1.66 > quit -bytes used=20239428, alloc=2096768, time=1.40 +bytes used=20238456, alloc=2096768, time=1.70 diff --git a/src/gr/setup_gr_gfas.maple b/src/gr/setup_gr_gfas.maple index ede82ca..dc079cb 100644 --- a/src/gr/setup_gr_gfas.maple +++ b/src/gr/setup_gr_gfas.maple @@ -22,7 +22,7 @@ make_gfa('K_dd', {inert}, [1..N, 1..N], symmetric); # index-raised and contracted metric and extrinsic curvature make_gfa('g_uu', {inert,fnd}, [1..N, 1..N], symmetric); make_gfa('K_uu', {inert,fnd}, [1..N, 1..N], symmetric); -make_gfa('K', {inert, fnd}, [], none); +make_gfa('K', {inert,fnd}, [], none); # xyz partial derivatives of metric # ... as far as Maple is concerned, these are indeed separate gridfns; @@ -37,22 +37,40 @@ make_gfa('partial_d_ln_sqrt_g', {inert,fnd}, [1..N], none); make_gfa('partial_d_g_uu', {inert,fnd}, [1..N, 1..N, 1..N], symmetric3_23); # radius of horizon -make_gfa('h', {inert, fnd}, [], none); +make_gfa('h', {inert,fnd}, [], none); # outward-pointing *non*-unit normal (covector) to horizon # and it's xyz-coordinate partial derivatives make_gfa('s_d', {inert,fnd}, [1..N], none); make_gfa('partial_d_s_d', {inert,fnd}, [1..N, 1..N], none); -# LHS of apparent horizon equation -make_gfa('H', {inert, fnd}, [], none); - # subexpressions for computing LHS of apparent horizon equation # ... these are defined by (15) in my 1996 apparent horizon finding paper -make_gfa('HA', {inert, fnd}, [], none); -make_gfa('HB', {inert, fnd}, [], none); -make_gfa('HC', {inert, fnd}, [], none); -make_gfa('HD', {inert, fnd}, [], none); +make_gfa('HA', {inert,fnd}, [], none); +make_gfa('HB', {inert,fnd}, [], none); +make_gfa('HC', {inert,fnd}, [], none); +make_gfa('HD', {inert,fnd}, [], none); + +# LHS of apparent horizon equation +make_gfa('H', {inert,fnd}, [], none); + +# Jacobian coefficients for H[ABCD] +# these are the partial derivatives of H[ABCD] +# ... wrt Diff(h,x_rs[x]) +make_gfa('partial_HA_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none); +make_gfa('partial_HB_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none); +make_gfa('partial_HC_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none); +make_gfa('partial_HD_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none); +# ... wrt Diff(h,x_rs[x],x_rs[y]) +make_gfa('partial_HA_wrt_partial_dd_h', {inert,fnd}, + [1..N_ang, 1..N_ang], symmetric); +make_gfa('partial_HB_wrt_partial_dd_h', {inert,fnd}, + [1..N_ang, 1..N_ang], symmetric); + +# Jacobian coefficients for H itself +make_gfa('partial_H_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none); +make_gfa('partial_H_wrt_partial_dd_h', {inert,fnd}, + [1..N_ang, 1..N_ang], symmetric); NULL; end proc; @@ -75,7 +93,7 @@ end proc; # var_seq = (in) (varargs) An expression sequence of the variables to # differentiate with respect to. # -`Diff/gridfn` := +`Diff/gridfn2` := proc(operand) # varargs option remember; # performance optimization global |