diff options
Diffstat (limited to 'src/patch/coords.cc')
-rw-r--r-- | src/patch/coords.cc | 282 |
1 files changed, 0 insertions, 282 deletions
diff --git a/src/patch/coords.cc b/src/patch/coords.cc index 578ab3d..5560936 100644 --- a/src/patch/coords.cc +++ b/src/patch/coords.cc @@ -2,17 +2,6 @@ // $Id$ // -// wr_coord::wr_coord -// wr_coord::pars::operator== -// wr_coord::pars::printme -// -// wr_coord::r_of_wr -//* r_of_wr_callback -// wr_coord::wr_of_r -//* wr_of_r_with_args - like wr_coord::wr_of_r() but with explicit args r0,a,b,c -// wr_coord::dwr_dr -// wr_coord::d2wr_dr -// // coords::{mu,nu,phi}_of_{{mu,nu,phi}} // coords::partial_{mu,nu,phi}_wrt_{mu,nu,phi}_at_const_{mu,nu,phi} // coords::partial2_{mu,nu,phi}_wrt_{{mu,nu,phi}}_at_const_{{mu,nu,phi}} @@ -36,287 +25,16 @@ #include <jt/stdc.h> #include <jt/util.h> #include <jt/util++.hh> -#include <jt/fortran.h> #include "policy.hh" #include "fp.hh" #include "misc.hh" #include "coords.hh" -#ifdef FP_IS_FLOAT - #include <jt/sfmm.h> -#endif -#ifdef FP_IS_DOUBLE - #include <jt/dfmm.h> -#endif - using jtutil::arctan_xy; using jtutil::pow2; using jtutil::hypot3; -// prototypes for functions local to this file -namespace { -fp mixed_210_wr_callback(const fp *pr); -fp wr_of_r_with_args(fp r0, fp a, fp b, fp c, fp r); - }; - -//***************************************************************************** -//***************************************************************************** -// wr_coord - wr(r) nonuniform gridding coordinate -//***************************************************************************** -//***************************************************************************** - -// -// This function constructs a wr_coord object. -// -wr_coord::wr_coord(const pars &pars_in) - : r0_(pars_in.r0), - a_(pars_in.a), b_(pars_in.b), c_(pars_in.c) -{ -} - -//***************************************************************************** - -// -// This function tests whether or not two struct pars objects -// are element-for-element equal. It returns true for equal, false for -// not-equal. -// -bool wr_coord::pars::operator==(const struct pars &rhs) -{ -return (r0 == rhs.r0) - && (a == rhs.a) - && (b == rhs.b) - && (c == rhs.c); -} - -//************************************* - -// -// This function tests whether or not a struct pars object is -// element-for-element equal to the corresponding "stuff" (members or -// access functions) of a wr_coord object. It returns true for equal, -// false for not-equal. -// -bool wr_coord::pars::operator==(const wr_coord &rhs) -{ -return (r0 == rhs.r0()) - && (a == rhs.a()) - && (b == rhs.b()) - && (c == rhs.c()); -} - -//***************************************************************************** - -// -// This function prints a struct_pars object in parameter-file format. -// -void wr_coord::pars::printme() const -{ -printf("\n"); -printf("wr_coord.r0\t%g\n", r0); -printf("wr_coord.a\t%g\n", a); -printf("wr_coord.b\t%g\n", b); -printf("wr_coord.c\t%g\n", c); -} - -//***************************************************************************** -//***************************************************************************** -//***************************************************************************** - -// -// This function computes the mixed-210 nonlinear grid warping function -// r = r(wr) by numerically inverting the wr = wr(r) function. We -// use Brent's zeroin() function, as given in -// G. E. Forsythe, M. A. Malcolm, and C. B. Moler -// "Computer Methods for Mathematical Computations" -// Prentice-Hall, Englewood Cliffs, 1977 -// with a bounding interval computed by successive doubling or halving -// in r, starting from the inner grid boundary r0. -// -// Arguments: -// {r0,a,b,c} = (in) The mixed-210 parameters. -// wr_star = (in) The mixed-210 wr coordinate. -// -// Result: -// The function returns the r coordinate. -// -// Bugs: -// - This function is somewhat inefficient. However, in practice it's -// only called during initialization to set up the coefficient arrays, -// so the inefficiency doesn't make much difference. -// - We implicitly assume that the Fortran zeroin() function uses <fp> -// precision; disaster is likely if this isn't true. -// - -// static variables for this function -static fp static_r0, static_a, static_b, static_c; -static fp static_wr_star; - -// the function itself -fp wr_coord::r_of_wr(fp wr_star) const -{ -// compute bounding interval by successive doubling or halving in r -// we assume that r=r0 corresponds to wr=0 -fp r_lo, r_hi; -if (wr_star >= 0.0) - then { - // search outwards from r0 by successive doubling in r - fp r; - for (r = 2.0 * r0() ; wr_of_r(r) <= wr_star ; r *= 2.0) - { - } - r_lo = 0.5 * r; - r_hi = r; - } - else { - // search inwards from r0 by successive *halving* in r - fp r; - for (r = 0.5 * r0() ; wr_of_r(r) >= wr_star ; r *= 0.5) - { - } - r_lo = r; - r_hi = 2.0 * r; - } - -// firewall check -if (! ((wr_of_r(r_lo) <= wr_star) && (wr_of_r(r_hi) >= wr_star)) ) - then error_exit(PANIC_EXIT, -"***** wr_coord::r_of_wr:\n" -" \"bounding interval\" [r_lo,r_hi] doesn't bound a solution!\n" -" (this should never happen!)\n" -" r0=%g a=%g b=%g c=%g\n" -" wr_star=%g\n" -" r_lo=%g wr_of_r(r0,a,b,c, r_lo)=%g\n" -" r_hi=%g wr_of_r(r0,a,b,c, r_hi)=%g\n" -, - double(r0()), double(a()), double(b()), double(c()), - double(wr_star), - double(r_lo), double(wr_of_r(r_lo)), - double(r_hi), double(wr_of_r(r_hi))); /*NOTREACHED*/ - -// solve wr(r) - wr_star = 0 for r -// using static variables to pass the parameters to the callback fn -static_r0 = r0(); -static_a = a(); -static_b = b(); -static_c = c(); -static_wr_star = wr_star; -fp tol = 0.0; // solve eqn as accurately as possible -return fmm::FORTRAN_NAME(zeroin)(&r_lo, &r_hi, &mixed_210_wr_callback, &tol); -} - -//***************************************************************************** - -// -// This is a callback function called by mixed_210_r_of_wr() (above) -// via Brent's zeroin() function. This function uses a C/Fortran-compatible -// interface. -// -// Arguments (implicit, as global variables): -// static_{r0,a,b,c} = (in) The {r0,a,b,c} mixed-210 parameters. -// static_wr_star = (in) The wr_star value from mixed_210_r_of_wr() . -// -// Arguments (explicit): -// pr --> (in) The r coordinate. -// -// Results: -// This function returns wr(r) - wr_star . -// -namespace { -fp mixed_210_wr_callback(const fp *pr) -{ -fp wr = wr_of_r_with_args(static_r0, static_a, static_b, static_c, *pr); -return wr - static_wr_star; -} - }; - -//***************************************************************************** - -// -// This function computes the mixed-210 grid warping function wr(r) , -// as per p.363 of my sssf notes. -// -// Arguments: -// r = (in) The r coordinate. -// -// Result: -// The function returns the mixed-210 wr coordinate. -// -fp wr_coord::wr_of_r(fp r) const -{ -return wr_of_r_with_args(r0(), a(), b(), c(), r); -} - -//***************************************************************************** - -// -// This function is the same as wr_coord::wr_of_r() (above), except -// that it's a non-member function with explicit arguments r0, a, b, and c. -// -// Arguments: -// {r0,a,b,c} = (in) The mixed-210 parameters. -// r = (in) The r coordinate. -// -// Result: -// The function returns the mixed-210 wr coordinate. -// -namespace { -fp wr_of_r_with_args(fp r0, fp a, fp b, fp c, fp r) -{ -fp x = r/r0; -fp a_term = (a == 0.0) ? 0.0 : (r0/a) * (1.0 - 1.0/x); -fp b_term = (b == 0.0) ? 0.0 : (r0/b) * log(x); -fp c_term = (c == 0.0) ? 0.0 : (r0/c) * (x - 1.0); - -return a_term + b_term + c_term; -} - }; - -//***************************************************************************** - -// -// This function computes the derivative of the mixed-210 grid warping -// function wr(r) , as per p.363 of my sssf notes. -// -// Arguments: -// r = (in) The r coordinate. -// -// Result: -// The function returns the derivative $d wr / d r$. -// -fp wr_coord::dwr_dr(fp r) const -{ -fp x = r/r0(); -fp a_term = (a() == 0.0) ? 0.0 : 1.0 / (a() * x*x); -fp b_term = (b() == 0.0) ? 0.0 : 1.0 / (b() * x); -fp c_term = (c() == 0.0) ? 0.0 : 1.0 / c(); - -return a_term + b_term + c_term; -} - -//***************************************************************************** - -// -// This function computes the second derivative of the mixed-210 grid warping -// function wr(r) , $d^2 wr / d r^2$. -// -// Arguments: -// {r0,a,b,c} = (in) The mixed-210 parameters. -// r = (in) The mixed-210 r coordinate. -// -// Result: -// The function returns the second derivative. -// -fp wr_coord::d2wr_dr2(fp r) const -{ -fp x = r/r0(); -fp a_term = (a() == 0.0) ? 0.0 : -2.0 / (a() * x*x*x * r0()); -fp b_term = (b() == 0.0) ? 0.0 : -1.0 / (b() * x*x * r0()); - -return( a_term + b_term ); -} - //***************************************************************************** //***************************************************************************** // coords - angular coordinate conversions and other misc coordinates stuff |