From 894ff0213c8d34b0ed45d11693687841e8eed74b Mon Sep 17 00:00:00 2001 From: jthorn Date: Tue, 15 Apr 2003 17:26:40 +0000 Subject: Switch to using new version of Lagrange polynomial interpolation which has separate tensor-product and maximum-degree interpolation operators, with the former now the default. This should fix the problem Yosef Zlowcher (sp?) found with the interpolant not being continuous in multiple dimensions. git-svn-id: http://svn.cactuscode.org/arrangements/CactusBase/LocalInterp/trunk@145 df1f8a13-aa1d-4dd4-9681-27ded5b42416 --- .../InterpLocalUniform.c | 238 ++++++++++++++++----- .../InterpLocalUniform.h | 54 +++-- src/GeneralizedPolynomial-Uniform/README | 23 +- src/GeneralizedPolynomial-Uniform/startup.c | 18 +- src/GeneralizedPolynomial-Uniform/template.c | 24 +++ 5 files changed, 276 insertions(+), 81 deletions(-) (limited to 'src') diff --git a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c index d0b4c78..9e0bf9f 100644 --- a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c +++ b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c @@ -3,8 +3,9 @@ /* ** *****data structures and functions local to this file ***** ** - * LocalInterp_ILU_Lagrange - external entry point for Lagrange interpolator - * LocalInterp_ILU_Hermite - external entry point for Hermite interpolator + * LocalInterp_ILU_Lagrange_TP - driver for Lagrange interpolator (tensor prod) + * LocalInterp_ILU_Lagrange_MD - driver for Lagrange interpolator (max degree) + * LocalInterp_ILU_Hermite - driver point for Hermite interpolator ** ** LocalInterp_InterpLocalUniform - main driver routine ** @@ -56,7 +57,8 @@ #include "util_Table.h" #include "InterpLocalUniform.h" -#include "Lagrange/all_prototypes.h" +#include "Lagrange-tensor-product/all_prototypes.h" +#include "Lagrange-maximum-degree/all_prototypes.h" #include "Hermite/all_prototypes.h" /* the rcs ID and its dummy function to use it */ @@ -78,8 +80,11 @@ CCTK_FILEVERSION(CactusBase_LocalInterp_src_GeneralizedPolynomialUniform_InterpL /* this enum describes which interpolation operator we're using */ enum interp_operator { - interp_operator_Lagrange, /* "Lagrange polynomial interpolation" */ - interp_operator_Hermite, /* "Hermite polynomial interpolation" */ + interp_operator_Lagrange_TP, /* "Lagrange polynomial interpolation */ + /* (tensor product)" */ + interp_operator_Lagrange_MD, /* "Lagrange polynomial interpolation + /* (maximum degree)" */ + interp_operator_Hermite, /* "Hermite polynomial interpolation" */ N_INTERP_OPERATORS /* this must be the last entry in the enum */ }; @@ -185,6 +190,8 @@ typedef * legal subscripts are of course 0...X. But for these axes we actually * only use 1...X. (Having the array oversize like this slightly simplifies * the array subscripting.) + * + * In the comments on the initializers, "n.i." = "not implemented". */ static const p_interp_fn_t p_interp_fn_table[N_INTERP_OPERATORS] [MAX_N_DIMS+1] /* see above */ @@ -193,9 +200,69 @@ static const p_interp_fn_t p_interp_fn_table[N_INTERP_OPERATORS] [MAX_SMOOTHING+1] = { - /* in the comments on the initializers, "n.i." = "not implemented" */ { - /* interp_operator == interp_operator_Lagrange */ + /* interp_operator == interp_operator_Lagrange_TP */ + { + /* N_dims = 0 ==> unused */ + { + /* molecule_family = molecule_family_cube */ + { NULL_INTERP_FN_PTR, }, /* order=0 (unused) */ + { NULL_INTERP_FN_PTR, }, /* order=1 (unused) */ + { NULL_INTERP_FN_PTR, }, /* order=2 (unused) */ + { NULL_INTERP_FN_PTR, }, /* order=3 (unused) */ + { NULL_INTERP_FN_PTR, }, /* order=4 (unused) */ + { NULL_INTERP_FN_PTR, }, /* order=5 (unused) */ + { NULL_INTERP_FN_PTR, }, /* order=6 (unused) */ + }, + }, + + { + /* N_dims = 1 */ + { + /* molecule_family = molecule_family_cube */ + { NULL_INTERP_FN_PTR, }, /* order=0 (unused) */ + { LocalInterp_U_LagrTP_1d_cube10, }, /* order=1 */ + { LocalInterp_U_LagrTP_1d_cube20, }, /* order=2 */ + { LocalInterp_U_LagrTP_1d_cube30, }, /* order=3 */ + { LocalInterp_U_LagrTP_1d_cube40, }, /* order=4 */ + { LocalInterp_U_LagrTP_1d_cube50, }, /* order=5 */ + { LocalInterp_U_LagrTP_1d_cube60, }, /* order=6 */ + }, + }, + + { + /* N_dims = 2 */ + { + /* molecule_family = molecule_family_cube */ + { NULL_INTERP_FN_PTR, }, /* order=0 (unused) */ + { LocalInterp_U_LagrTP_2d_cube10, }, /* order=1 */ + { LocalInterp_U_LagrTP_2d_cube20, }, /* order=2 */ + { LocalInterp_U_LagrTP_2d_cube30, }, /* order=3 */ + { LocalInterp_U_LagrTP_2d_cube40, }, /* order=4 */ + { NULL_INTERP_FN_PTR, }, /* order=5 (n.i.) */ + { NULL_INTERP_FN_PTR, }, /* order=6 (n.i.) */ + }, + }, + + { + /* N_dims = 3 */ + { + /* molecule_family = molecule_family_cube */ + { NULL_INTERP_FN_PTR, }, /* order=0 (unused) */ + { LocalInterp_U_LagrTP_3d_cube10, }, /* order=1 */ + { LocalInterp_U_LagrTP_3d_cube20, }, /* order=2 */ + { LocalInterp_U_LagrTP_3d_cube30, }, /* order=3 */ + { LocalInterp_U_LagrTP_3d_cube40, }, /* order=4 */ + { NULL_INTERP_FN_PTR, }, /* order=5 (n.i.) */ + { NULL_INTERP_FN_PTR, }, /* order=6 (n.i.) */ + }, + }, + + /* end of interp_operator == interp_operator_Lagrange_TP */ + }, + + { + /* interp_operator == interp_operator_Lagrange_MD */ { /* N_dims = 0 ==> unused */ { @@ -215,12 +282,12 @@ static const p_interp_fn_t p_interp_fn_table[N_INTERP_OPERATORS] { /* molecule_family = molecule_family_cube */ { NULL_INTERP_FN_PTR, }, /* order=0 (unused) */ - { LocalInterp_ULagrange_1d_cube10, }, /* order=1 */ - { LocalInterp_ULagrange_1d_cube20, }, /* order=2 */ - { LocalInterp_ULagrange_1d_cube30, }, /* order=3 */ - { LocalInterp_ULagrange_1d_cube40, }, /* order=4 */ - { LocalInterp_ULagrange_1d_cube50, }, /* order=5 */ - { LocalInterp_ULagrange_1d_cube60, }, /* order=6 */ + { LocalInterp_U_LagrMD_1d_cube10, }, /* order=1 */ + { LocalInterp_U_LagrMD_1d_cube20, }, /* order=2 */ + { LocalInterp_U_LagrMD_1d_cube30, }, /* order=3 */ + { LocalInterp_U_LagrMD_1d_cube40, }, /* order=4 */ + { LocalInterp_U_LagrMD_1d_cube50, }, /* order=5 */ + { LocalInterp_U_LagrMD_1d_cube60, }, /* order=6 */ }, }, @@ -229,10 +296,10 @@ static const p_interp_fn_t p_interp_fn_table[N_INTERP_OPERATORS] { /* molecule_family = molecule_family_cube */ { NULL_INTERP_FN_PTR, }, /* order=0 (unused) */ - { LocalInterp_ULagrange_2d_cube10, }, /* order=1 */ - { LocalInterp_ULagrange_2d_cube20, }, /* order=2 */ - { LocalInterp_ULagrange_2d_cube30, }, /* order=3 */ - { LocalInterp_ULagrange_2d_cube40, }, /* order=4 */ + { LocalInterp_U_LagrMD_2d_cube10, }, /* order=1 */ + { LocalInterp_U_LagrMD_2d_cube20, }, /* order=2 */ + { LocalInterp_U_LagrMD_2d_cube30, }, /* order=3 */ + { LocalInterp_U_LagrMD_2d_cube40, }, /* order=4 */ { NULL_INTERP_FN_PTR, }, /* order=5 (n.i.) */ { NULL_INTERP_FN_PTR, }, /* order=6 (n.i.) */ }, @@ -243,16 +310,16 @@ static const p_interp_fn_t p_interp_fn_table[N_INTERP_OPERATORS] { /* molecule_family = molecule_family_cube */ { NULL_INTERP_FN_PTR, }, /* order=0 (unused) */ - { LocalInterp_ULagrange_3d_cube10, }, /* order=1 */ - { LocalInterp_ULagrange_3d_cube20, }, /* order=2 */ - { LocalInterp_ULagrange_3d_cube30, }, /* order=3 */ - { LocalInterp_ULagrange_3d_cube40, }, /* order=4 */ + { LocalInterp_U_LagrMD_3d_cube10, }, /* order=1 */ + { LocalInterp_U_LagrMD_3d_cube20, }, /* order=2 */ + { LocalInterp_U_LagrMD_3d_cube30, }, /* order=3 */ + { LocalInterp_U_LagrMD_3d_cube40, }, /* order=4 */ { NULL_INTERP_FN_PTR, }, /* order=5 (n.i.) */ { NULL_INTERP_FN_PTR, }, /* order=6 (n.i.) */ }, }, - /* end of interp_operator == interp_operator_Lagrange */ + /* end of interp_operator == interp_operator_Lagrange_MD */ }, { @@ -314,7 +381,7 @@ static const p_interp_fn_t p_interp_fn_table[N_INTERP_OPERATORS] }, /* end of interp_operator == interp_operator_Hermite */ - } + }, }; @@ -323,7 +390,7 @@ static const p_interp_fn_t p_interp_fn_table[N_INTERP_OPERATORS] /******************************************************************************/ /*@@ - @routine LocalInterp_ILU_Lagrange + @routine LocalInterp_ILU_Lagrange_TP @date 22 Oct 2001 @author Jonathan Thornburg @desc @@ -332,10 +399,85 @@ static const p_interp_fn_t p_interp_fn_table[N_INTERP_OPERATORS] It does Lagrange interpolation of a set of input arrays defined on a uniform N-dimensional tensor-product grid, to a set of - interpolation points. It can also do differentiation and/or - smoothing in the process of the interpolation. It can also - return information about the Jacobian of the interpolated - values with respect to the input arrays. + interpolation points. For multidimensional interpolation, it + uses the tensor-product basis for the interpolating function, + so the interpolant is guaranteed to be continuous and to pass + through the input data at the grid points (up to floating-point + roundoff errors). + + This interpolator can also + * do differentiation and/or smoothing in the process of the + interpolation + * return information about the Jacobian of the interpolated + values with respect to the input arrays + + This function does nothing except pass all its arguments + down to LocalInterp_InterpLocalUniform() with an extra + 2 arguments added to indicate the operator handle. See that + function for details on this function's arguments/results. + @enddesc + @@*/ +int LocalInterp_ILU_Lagrange_TP(int N_dims, + int param_table_handle, + /***** coordinate system *****/ + const CCTK_REAL coord_origin[], + const CCTK_REAL coord_delta[], + /***** interpolation points *****/ + int N_interp_points, + int interp_coords_type_code, + const void *const interp_coords[], + /***** input arrays *****/ + int N_input_arrays, + const CCTK_INT input_array_dims[], + const CCTK_INT input_array_type_codes[], + const void *const input_arrays[], + /***** output arrays *****/ + int N_output_arrays, + const CCTK_INT output_array_type_codes[], + void *const output_arrays[]) +{ +return LocalInterp_InterpLocalUniform(interp_operator_Lagrange_TP, + "Lagrange polynomial interpolation (tensor product)", + LAGRANGE_BNDRY_OFF_CNTR_TOL_DEF, + LAGRANGE_BNDRY_EXTRAP_TOL_DEF, + N_dims, + param_table_handle, + coord_origin, coord_delta, + N_interp_points, + interp_coords_type_code, + interp_coords, + N_input_arrays, + input_array_dims, + input_array_type_codes, + input_arrays, + N_output_arrays, + output_array_type_codes, + output_arrays); +} + +/******************************************************************************/ + +/*@@ + @routine LocalInterp_ILU_Lagrange_MD + @date 22 Oct 2001 + @author Jonathan Thornburg + @desc + This function is a top-level driver for + LocalInterp::CCTK_InterpLocalUniform(). + + It does Lagrange interpolation of a set of input arrays defined + on a uniform N-dimensional tensor-product grid, to a set of + interpolation points. For multidimensional interpolation, it + uses the maximum-degree basis for the interpolating function. + This means the interpolant is in general *not* continuous and + in general does *not* pass through the input data at the grid + points. + + This interpolator can also + * do differentiation and/or smoothing in the process of the + interpolation + * return information about the Jacobian of the interpolated + values with respect to the input arrays This function does nothing except pass all its arguments down to LocalInterp_InterpLocalUniform() with an extra @@ -343,27 +485,27 @@ static const p_interp_fn_t p_interp_fn_table[N_INTERP_OPERATORS] function for details on this function's arguments/results. @enddesc @@*/ -int LocalInterp_ILU_Lagrange(int N_dims, - int param_table_handle, - /***** coordinate system *****/ - const CCTK_REAL coord_origin[], - const CCTK_REAL coord_delta[], - /***** interpolation points *****/ - int N_interp_points, - int interp_coords_type_code, - const void *const interp_coords[], - /***** input arrays *****/ - int N_input_arrays, - const CCTK_INT input_array_dims[], - const CCTK_INT input_array_type_codes[], - const void *const input_arrays[], - /***** output arrays *****/ - int N_output_arrays, - const CCTK_INT output_array_type_codes[], - void *const output_arrays[]) +int LocalInterp_ILU_Lagrange_MD(int N_dims, + int param_table_handle, + /***** coordinate system *****/ + const CCTK_REAL coord_origin[], + const CCTK_REAL coord_delta[], + /***** interpolation points *****/ + int N_interp_points, + int interp_coords_type_code, + const void *const interp_coords[], + /***** input arrays *****/ + int N_input_arrays, + const CCTK_INT input_array_dims[], + const CCTK_INT input_array_type_codes[], + const void *const input_arrays[], + /***** output arrays *****/ + int N_output_arrays, + const CCTK_INT output_array_type_codes[], + void *const output_arrays[]) { -return LocalInterp_InterpLocalUniform(interp_operator_Lagrange, - "Lagrange polynomial interpolation", +return LocalInterp_InterpLocalUniform(interp_operator_Lagrange_MD, + "Lagrange polynomial interpolation (maximum degree)", LAGRANGE_BNDRY_OFF_CNTR_TOL_DEF, LAGRANGE_BNDRY_EXTRAP_TOL_DEF, N_dims, diff --git a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.h b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.h index 8d9c551..d4821c4 100644 --- a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.h +++ b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.h @@ -192,24 +192,42 @@ struct Jacobian_info */ /* functions in InterpLocalUniform.c */ -int LocalInterp_ILU_Lagrange(int N_dims, - int param_table_handle, - /***** coordinate system *****/ - const CCTK_REAL coord_origin[], - const CCTK_REAL coord_delta[], - /***** interpolation points *****/ - int N_interp_points, - int interp_coords_type_code, - const void *const interp_coords[], - /***** input arrays *****/ - int N_input_arrays, - const CCTK_INT input_array_dims[], - const CCTK_INT input_array_type_codes[], - const void *const input_arrays[], - /***** output arrays *****/ - int N_output_arrays, - const CCTK_INT output_array_type_codes[], - void *const output_arrays[]); +int LocalInterp_ILU_Lagrange_TP(int N_dims, + int param_table_handle, + /***** coordinate system *****/ + const CCTK_REAL coord_origin[], + const CCTK_REAL coord_delta[], + /***** interpolation points *****/ + int N_interp_points, + int interp_coords_type_code, + const void *const interp_coords[], + /***** input arrays *****/ + int N_input_arrays, + const CCTK_INT input_array_dims[], + const CCTK_INT input_array_type_codes[], + const void *const input_arrays[], + /***** output arrays *****/ + int N_output_arrays, + const CCTK_INT output_array_type_codes[], + void *const output_arrays[]); +int LocalInterp_ILU_Lagrange_MD(int N_dims, + int param_table_handle, + /***** coordinate system *****/ + const CCTK_REAL coord_origin[], + const CCTK_REAL coord_delta[], + /***** interpolation points *****/ + int N_interp_points, + int interp_coords_type_code, + const void *const interp_coords[], + /***** input arrays *****/ + int N_input_arrays, + const CCTK_INT input_array_dims[], + const CCTK_INT input_array_type_codes[], + const void *const input_arrays[], + /***** output arrays *****/ + int N_output_arrays, + const CCTK_INT output_array_type_codes[], + void *const output_arrays[]); int LocalInterp_ILU_Hermite(int N_dims, int param_table_handle, /***** coordinate system *****/ diff --git a/src/GeneralizedPolynomial-Uniform/README b/src/GeneralizedPolynomial-Uniform/README index b90cc26..15c6d27 100644 --- a/src/GeneralizedPolynomial-Uniform/README +++ b/src/GeneralizedPolynomial-Uniform/README @@ -3,8 +3,9 @@ $Header$ This directory contains the interpolator registered for CCTK_InterpLocalUniform() under the names - "Lagrange polynomial interpolation" -and + "Lagrange polynomial interpolation" # a synonym for the next one... + "Lagrange polynomial interpolation (tensor product)" + "Lagrange polynomial interpolation (maximum degree)" "Hermite polynomial interpolation" @@ -47,10 +48,13 @@ The source code files are as follows: test cases stored in a table in the code. The subdirectories are as follows: -* [123]d.coeffs/ are empty directory trees left-over from previous - versions of this interpolator (CVS doesn't have any easy way to - remove directories) -* Lagrange/ contains the code for the Lagrange polynomial interpolator. +* Lagrange/ and [123]d.coeffs/ are empty directory trees left-over + from previous versions of this interpolator (CVS doesn't have any + easy way to remove directories) +* Lagrange-tensor-product/ contains the code for the Lagrange polynomial + interpolator using tensor-product bases for multiple dimensions. +* Lagrange-maximum-degree/ contains the code for the Lagrange polynomial + interpolator using maximum-degree bases for multiple dimensions. * Hermite/ contains the code for the Hermite polynomial interpolator. * common/ contains low-level code common to both the Lagrange and the Hermite interpolators. @@ -63,9 +67,10 @@ to this interpolator, you need to * edit the makefile to create the appropriate directory or directories, and if necessary add any new targets * 'make' as appropriate to create the new coefficients etc - [note this can be rather computationally intensive -- the current - set of {1d,2d,3d} coefficients take {4 seconds, 15 seconds, 8 minutes} - to generate using Maple 7 on a xeon (1.7 GHz Pentium 4)] + [note this can be rather computationally intensive -- for example, + the set of {1d,2d,3d} Lagrange-tensor-product coefficients take about + {4 seconds, 40 seconds, 20 minutes} to generate using Maple 7 on + an AEI xeon (1.7 GHz Pentium 4 with 1GB memory)] * edit InterpLocalUniform.c to add the new case to the decoding switch statements * create an appropriate "script" file which defines the right macros, diff --git a/src/GeneralizedPolynomial-Uniform/startup.c b/src/GeneralizedPolynomial-Uniform/startup.c index 0c6be52..5b1ff01 100644 --- a/src/GeneralizedPolynomial-Uniform/startup.c +++ b/src/GeneralizedPolynomial-Uniform/startup.c @@ -33,16 +33,22 @@ void LocalInterp_GPU_Startup(void); @@*/ void LocalInterp_GPU_Startup(void) { -CCTK_InterpRegisterOpLocalUniform(LocalInterp_ILU_Lagrange, - "Lagrange polynomial interpolation", +CCTK_InterpRegisterOpLocalUniform(LocalInterp_ILU_Lagrange_TP, + "Lagrange polynomial interpolation (tensor product)", CCTK_THORNSTRING); - -/* old operator name for backwards compatability */ -CCTK_InterpRegisterOpLocalUniform(LocalInterp_ILU_Lagrange, - "generalized polynomial interpolation", +CCTK_InterpRegisterOpLocalUniform(LocalInterp_ILU_Lagrange_MD, + "Lagrange polynomial interpolation (maximum degree)", CCTK_THORNSTRING); CCTK_InterpRegisterOpLocalUniform(LocalInterp_ILU_Hermite, "Hermite polynomial interpolation", CCTK_THORNSTRING); + +/* synonym operator names for backwards compatability */ +CCTK_InterpRegisterOpLocalUniform(LocalInterp_ILU_Lagrange_TP, + "Lagrange polynomial interpolation", + CCTK_THORNSTRING); +CCTK_InterpRegisterOpLocalUniform(LocalInterp_ILU_Lagrange_TP, + "generalized polynomial interpolation", + CCTK_THORNSTRING); } diff --git a/src/GeneralizedPolynomial-Uniform/template.c b/src/GeneralizedPolynomial-Uniform/template.c index c12a308..454f733 100644 --- a/src/GeneralizedPolynomial-Uniform/template.c +++ b/src/GeneralizedPolynomial-Uniform/template.c @@ -863,6 +863,19 @@ int pt; } } +#ifdef PICKLE + #if (N_DIMS == 1) +printf("pt=%d interp coord %.16g\n", + pt, (double) interp_coords_fp[0]); + #endif + #if (N_DIMS == 3) +printf("pt=%d interp coords %.16g %.16g %.16g\n", + pt, (double) interp_coords_fp[0], + (double) interp_coords_fp[1], + (double) interp_coords_fp[2]); + #endif +#endif + /* * compute position of interpolation molecules with respect to @@ -934,6 +947,17 @@ int pt; const fp z = xyz_temp [Z_AXIS]; #endif +#ifdef PICKLE + #if (N_DIMS == 1) +printf("interp center_i = %d\n", center_i); +printf("interp grid-relative x = %.16g\n", (double) x); + #endif + #if (N_DIMS == 3) +printf("interp center_ijk = %d %d %d\n", center_i, center_j, center_k); +printf("interp grid-relative xyz = %.16g %.16g %.16g\n", + (double) x, (double) y, (double) z); + #endif +#endif /* * compute 1-d position of molecule center in input arrays -- cgit v1.2.3