From 6d2aca99e8b9ff712ee47ebe23af38486e1b7d76 Mon Sep 17 00:00:00 2001 From: tradke Date: Wed, 14 Feb 2001 14:40:43 +0000 Subject: This commit was generated by cvs2svn to compensate for changes in r2, which included commits to RCS files with non-trunk default branches. git-svn-id: http://svn.cactuscode.org/arrangements/CactusPUGH/PUGHInterp/trunk@3 1c20744c-e24a-42ec-9533-f5004cb800e5 --- COPYRIGHT | 352 ++++++++++++++ README | 18 + doc/documentation.tex | 66 +++ interface.ccl | 4 + par/interp.par | 6 + param.ccl | 9 + schedule.ccl | 7 + src/Interpolate.c | 427 ++++++++++++++++ src/Operator.c | 1293 +++++++++++++++++++++++++++++++++++++++++++++++++ src/Startup.c | 253 ++++++++++ src/make.code.defn | 5 + src/pughInterpGH.h | 71 +++ 12 files changed, 2511 insertions(+) create mode 100644 COPYRIGHT create mode 100644 README create mode 100644 doc/documentation.tex create mode 100644 interface.ccl create mode 100644 par/interp.par create mode 100644 param.ccl create mode 100644 schedule.ccl create mode 100644 src/Interpolate.c create mode 100644 src/Operator.c create mode 100644 src/Startup.c create mode 100644 src/make.code.defn create mode 100644 src/pughInterpGH.h diff --git a/COPYRIGHT b/COPYRIGHT new file mode 100644 index 0000000..57b2896 --- /dev/null +++ b/COPYRIGHT @@ -0,0 +1,352 @@ + + +This thorn is (c) Copyright the authors listed in the sources files. This +thorn is distributed under the GNU GPL with the specific exception that +the GNU GPL does not migrate to code linking to or using infrastructure +from this thorn. + +- The Cactus Team + +************************************************************************** + + GNU GENERAL PUBLIC LICENSE + Version 2, June 1991 + + Copyright (C) 1989, 1991 Free Software Foundation, Inc. + 675 Mass Ave, Cambridge, MA 02139, USA + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + Preamble + + The licenses for most software are designed to take away your +freedom to share and change it. By contrast, the GNU General Public +License is intended to guarantee your freedom to share and change free +software--to make sure the software is free for all its users. This +General Public License applies to most of the Free Software +Foundation's software and to any other program whose authors commit to +using it. (Some other Free Software Foundation software is covered by +the GNU Library General Public License instead.) You can apply it to +your programs, too. + + When we speak of free software, we are referring to freedom, not +price. Our General Public Licenses are designed to make sure that you +have the freedom to distribute copies of free software (and charge for +this service if you wish), that you receive source code or can get it +if you want it, that you can change the software or use pieces of it +in new free programs; and that you know you can do these things. + + To protect your rights, we need to make restrictions that forbid +anyone to deny you these rights or to ask you to surrender the rights. +These restrictions translate to certain responsibilities for you if you +distribute copies of the software, or if you modify it. + + For example, if you distribute copies of such a program, whether +gratis or for a fee, you must give the recipients all the rights that +you have. You must make sure that they, too, receive or can get the +source code. And you must show them these terms so they know their +rights. + + We protect your rights with two steps: (1) copyright the software, and +(2) offer you this license which gives you legal permission to copy, +distribute and/or modify the software. + + Also, for each author's protection and ours, we want to make certain +that everyone understands that there is no warranty for this free +software. If the software is modified by someone else and passed on, we +want its recipients to know that what they have is not the original, so +that any problems introduced by others will not reflect on the original +authors' reputations. + + Finally, any free program is threatened constantly by software +patents. We wish to avoid the danger that redistributors of a free +program will individually obtain patent licenses, in effect making the +program proprietary. To prevent this, we have made it clear that any +patent must be licensed for everyone's free use or not licensed at all. + + The precise terms and conditions for copying, distribution and +modification follow. + + GNU GENERAL PUBLIC LICENSE + TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION + + 0. This License applies to any program or other work which contains +a notice placed by the copyright holder saying it may be distributed +under the terms of this General Public License. The "Program", below, +refers to any such program or work, and a "work based on the Program" +means either the Program or any derivative work under copyright law: +that is to say, a work containing the Program or a portion of it, +either verbatim or with modifications and/or translated into another +language. (Hereinafter, translation is included without limitation in +the term "modification".) Each licensee is addressed as "you". + +Activities other than copying, distribution and modification are not +covered by this License; they are outside its scope. The act of +running the Program is not restricted, and the output from the Program +is covered only if its contents constitute a work based on the +Program (independent of having been made by running the Program). +Whether that is true depends on what the Program does. + + 1. You may copy and distribute verbatim copies of the Program's +source code as you receive it, in any medium, provided that you +conspicuously and appropriately publish on each copy an appropriate +copyright notice and disclaimer of warranty; keep intact all the +notices that refer to this License and to the absence of any warranty; +and give any other recipients of the Program a copy of this License +along with the Program. + +You may charge a fee for the physical act of transferring a copy, and +you may at your option offer warranty protection in exchange for a fee. + + 2. You may modify your copy or copies of the Program or any portion +of it, thus forming a work based on the Program, and copy and +distribute such modifications or work under the terms of Section 1 +above, provided that you also meet all of these conditions: + + a) You must cause the modified files to carry prominent notices + stating that you changed the files and the date of any change. + + b) You must cause any work that you distribute or publish, that in + whole or in part contains or is derived from the Program or any + part thereof, to be licensed as a whole at no charge to all third + parties under the terms of this License. + + c) If the modified program normally reads commands interactively + when run, you must cause it, when started running for such + interactive use in the most ordinary way, to print or display an + announcement including an appropriate copyright notice and a + notice that there is no warranty (or else, saying that you provide + a warranty) and that users may redistribute the program under + these conditions, and telling the user how to view a copy of this + License. (Exception: if the Program itself is interactive but + does not normally print such an announcement, your work based on + the Program is not required to print an announcement.) + +These requirements apply to the modified work as a whole. If +identifiable sections of that work are not derived from the Program, +and can be reasonably considered independent and separate works in +themselves, then this License, and its terms, do not apply to those +sections when you distribute them as separate works. But when you +distribute the same sections as part of a whole which is a work based +on the Program, the distribution of the whole must be on the terms of +this License, whose permissions for other licensees extend to the +entire whole, and thus to each and every part regardless of who wrote it. + +Thus, it is not the intent of this section to claim rights or contest +your rights to work written entirely by you; rather, the intent is to +exercise the right to control the distribution of derivative or +collective works based on the Program. + +In addition, mere aggregation of another work not based on the Program +with the Program (or with a work based on the Program) on a volume of +a storage or distribution medium does not bring the other work under +the scope of this License. + + 3. You may copy and distribute the Program (or a work based on it, +under Section 2) in object code or executable form under the terms of +Sections 1 and 2 above provided that you also do one of the following: + + a) Accompany it with the complete corresponding machine-readable + source code, which must be distributed under the terms of Sections + 1 and 2 above on a medium customarily used for software interchange; or, + + b) Accompany it with a written offer, valid for at least three + years, to give any third party, for a charge no more than your + cost of physically performing source distribution, a complete + machine-readable copy of the corresponding source code, to be + distributed under the terms of Sections 1 and 2 above on a medium + customarily used for software interchange; or, + + c) Accompany it with the information you received as to the offer + to distribute corresponding source code. (This alternative is + allowed only for noncommercial distribution and only if you + received the program in object code or executable form with such + an offer, in accord with Subsection b above.) + +The source code for a work means the preferred form of the work for +making modifications to it. For an executable work, complete source +code means all the source code for all modules it contains, plus any +associated interface definition files, plus the scripts used to +control compilation and installation of the executable. However, as a +special exception, the source code distributed need not include +anything that is normally distributed (in either source or binary +form) with the major components (compiler, kernel, and so on) of the +operating system on which the executable runs, unless that component +itself accompanies the executable. + +If distribution of executable or object code is made by offering +access to copy from a designated place, then offering equivalent +access to copy the source code from the same place counts as +distribution of the source code, even though third parties are not +compelled to copy the source along with the object code. + + 4. You may not copy, modify, sublicense, or distribute the Program +except as expressly provided under this License. Any attempt +otherwise to copy, modify, sublicense or distribute the Program is +void, and will automatically terminate your rights under this License. +However, parties who have received copies, or rights, from you under +this License will not have their licenses terminated so long as such +parties remain in full compliance. + + 5. You are not required to accept this License, since you have not +signed it. However, nothing else grants you permission to modify or +distribute the Program or its derivative works. These actions are +prohibited by law if you do not accept this License. Therefore, by +modifying or distributing the Program (or any work based on the +Program), you indicate your acceptance of this License to do so, and +all its terms and conditions for copying, distributing or modifying +the Program or works based on it. + + 6. Each time you redistribute the Program (or any work based on the +Program), the recipient automatically receives a license from the +original licensor to copy, distribute or modify the Program subject to +these terms and conditions. You may not impose any further +restrictions on the recipients' exercise of the rights granted herein. +You are not responsible for enforcing compliance by third parties to +this License. + + 7. If, as a consequence of a court judgment or allegation of patent +infringement or for any other reason (not limited to patent issues), +conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot +distribute so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you +may not distribute the Program at all. For example, if a patent +license would not permit royalty-free redistribution of the Program by +all those who receive copies directly or indirectly through you, then +the only way you could satisfy both it and this License would be to +refrain entirely from distribution of the Program. + +If any portion of this section is held invalid or unenforceable under +any particular circumstance, the balance of the section is intended to +apply and the section as a whole is intended to apply in other +circumstances. + +It is not the purpose of this section to induce you to infringe any +patents or other property right claims or to contest validity of any +such claims; this section has the sole purpose of protecting the +integrity of the free software distribution system, which is +implemented by public license practices. Many people have made +generous contributions to the wide range of software distributed +through that system in reliance on consistent application of that +system; it is up to the author/donor to decide if he or she is willing +to distribute software through any other system and a licensee cannot +impose that choice. + +This section is intended to make thoroughly clear what is believed to +be a consequence of the rest of this License. + + 8. If the distribution and/or use of the Program is restricted in +certain countries either by patents or by copyrighted interfaces, the +original copyright holder who places the Program under this License +may add an explicit geographical distribution limitation excluding +those countries, so that distribution is permitted only in or among +countries not thus excluded. In such case, this License incorporates +the limitation as if written in the body of this License. + + 9. The Free Software Foundation may publish revised and/or new versions +of the General Public License from time to time. Such new versions will +be similar in spirit to the present version, but may differ in detail to +address new problems or concerns. + +Each version is given a distinguishing version number. If the Program +specifies a version number of this License which applies to it and "any +later version", you have the option of following the terms and conditions +either of that version or of any later version published by the Free +Software Foundation. If the Program does not specify a version number of +this License, you may choose any version ever published by the Free Software +Foundation. + + 10. If you wish to incorporate parts of the Program into other free +programs whose distribution conditions are different, write to the author +to ask for permission. For software which is copyrighted by the Free +Software Foundation, write to the Free Software Foundation; we sometimes +make exceptions for this. Our decision will be guided by the two goals +of preserving the free status of all derivatives of our free software and +of promoting the sharing and reuse of software generally. + + NO WARRANTY + + 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY +FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN +OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES +PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED +OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF +MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS +TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE +PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, +REPAIR OR CORRECTION. + + 12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING +WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR +REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, +INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING +OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED +TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY +YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER +PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE +POSSIBILITY OF SUCH DAMAGES. + + END OF TERMS AND CONDITIONS + + Appendix: How to Apply These Terms to Your New Programs + + If you develop a new program, and you want it to be of the greatest +possible use to the public, the best way to achieve this is to make it +free software which everyone can redistribute and change under these terms. + + To do so, attach the following notices to the program. It is safest +to attach them to the start of each source file to most effectively +convey the exclusion of warranty; and each file should have at least +the "copyright" line and a pointer to where the full notice is found. + + + Copyright (C) 19yy + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. + +Also add information on how to contact you by electronic and paper mail. + +If the program is interactive, make it output a short notice like this +when it starts in an interactive mode: + + Gnomovision version 69, Copyright (C) 19yy name of author + Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'. + This is free software, and you are welcome to redistribute it + under certain conditions; type `show c' for details. + +The hypothetical commands `show w' and `show c' should show the appropriate +parts of the General Public License. Of course, the commands you use may +be called something other than `show w' and `show c'; they could even be +mouse-clicks or menu items--whatever suits your program. + +You should also get your employer (if you work as a programmer) or your +school, if any, to sign a "copyright disclaimer" for the program, if +necessary. Here is a sample; alter the names: + + Yoyodyne, Inc., hereby disclaims all copyright interest in the program + `Gnomovision' (which makes passes at compilers) written by James Hacker. + + , 1 April 1989 + Ty Coon, President of Vice + +This General Public License does not permit incorporating your program into +proprietary programs. If your program is a subroutine library, you may +consider it more useful to permit linking proprietary applications with the +library. If this is what you want to do, use the GNU Library General +Public License instead of this License. + + diff --git a/README b/README new file mode 100644 index 0000000..b2b21d1 --- /dev/null +++ b/README @@ -0,0 +1,18 @@ +Cactus Code Thorn PUGHInterp +Authors : Paul Walker and Thomas Radke, some optimization ideas from Erik Schnetter +CVS info : $Header$ +-------------------------------------------------------------------------- + +1. Purpose of the thorn + + This thorn provides interpolation of arrays at arbitrary points + on a regular, uniform cartesian grid. + + +2. Additional Information + + More detailed information on thorn PUGHInterp is given in the + doc/documentation.tex file. + + For information on how to invoke interpolation operators via the + general flesh interpolation API please refer to the flesh documentation. diff --git a/doc/documentation.tex b/doc/documentation.tex new file mode 100644 index 0000000..b59a7a3 --- /dev/null +++ b/doc/documentation.tex @@ -0,0 +1,66 @@ +%version $Header$ +\documentclass{article} +\begin{document} + +\title{PUGHInterp} +\author{Paul Walker, Thomas Radke, Erik Schnetter} +\date{1997-2001} +\maketitle + +\abstract{Thorn PUGHInterp provides interpolation of arrays + at arbitrary points.} + +\section{Purpose} +Thorn PUGHInterp implements interpolation operators working on regular, +uniform cartesian grids. They can be applied to interpolate both CCTK grid +variables (grid functions and arrays distributed over all processors), and +processor-local arrays at arbitrary points.\\ +The operators for distributed/local arrays are invoked via the flesh +interpolation API routines {\tt CCTK\_InterpGV()} and {\tt CCTK\_InterpLocal()} +resp. They are registered with the flesh under the name {\it ''uniform +cartesian''} prepended by the interpolation order (eg. {\it ''second-order +uniform cartesian''}). Currently there is first, second, and third-order +interpolation implemented.\\ + +\section{Implementation Notes} +The interpolation operators registered for different orders are mapped +via wrappers (in {\tt Startup.c}) onto a single routine (in {\tt Operator.c}) +just passing the order as an additional argument.\\ +The routine for distributed arrays will then map all points to interpolate +at to the processors which own those points, and communicate the points' +coordinates and corresponding input arrays ({\tt MPI\_Alltoall()} is used +for this).\\ +Then the interpolation takes place in parallel on every processor, calling +a core interpolation routine (located in {\tt Interpolate.c}). This one +takes a list of input arrays and points and interpolates these to a +list of output arrays (one output value per interpolation point).\\ +Again, for distributed arrays, the interpolation results for remote points +are sent back to the requesting processors.\\[2ex] +% +Current limitations of the core interpolation routine's implementation are: +% +\begin{itemize} + \item arrays up to three ({\bf MAXDIM}) dimensions only can be handled + \item interpolation orders up to three ({\bf MAXORDER}) only are supported + \item coordinates must be given as {\bf CCTK\_REAL} types + \item input and output array types must be the same + (no type casting of interpolation results supported) +\end{itemize} +% +Despite of these limitations, the code it was programmed almost generic +in that it can easily be extended to support higher-dimensional arrays +or more interpolation orders.\\ +Please see the NOTES in this source file for details. + +\section{Comments} +% +For more information on how to invoke interpolation operators please refer +to the flesh documentation. +% +% Automatically created from the ccl files +% Do not worry for now. +\include{interface} +\include{param} +\include{schedule} + +\end{document} diff --git a/interface.ccl b/interface.ccl new file mode 100644 index 0000000..565e2d8 --- /dev/null +++ b/interface.ccl @@ -0,0 +1,4 @@ +# Interface definition for thorn PUGHInterp +# $Header$ + +implements: PUGHInterp diff --git a/par/interp.par b/par/interp.par new file mode 100644 index 0000000..f1ff599 --- /dev/null +++ b/par/interp.par @@ -0,0 +1,6 @@ +# interp.par +# + +ActiveThorns = "PUGH PUGHInterp CartGrid3D" + +PUGHInterp::demo = "yes" diff --git a/param.ccl b/param.ccl new file mode 100644 index 0000000..9f1c756 --- /dev/null +++ b/param.ccl @@ -0,0 +1,9 @@ +# Parameter definitions for thorn PUGHInterp +# $Header$ + +############################################################################# +### import CartGrid3D parameters +############################################################################# +shares: grid + +USES KEYWORD domain diff --git a/schedule.ccl b/schedule.ccl new file mode 100644 index 0000000..5e5a7bc --- /dev/null +++ b/schedule.ccl @@ -0,0 +1,7 @@ +# Schedule definitions for thorn PUGHInterp +# $Header$ + +schedule PUGHInterp_Startup at STARTUP after Driver_Startup +{ + LANG:C +} "PUGHInterp startup routine" diff --git a/src/Interpolate.c b/src/Interpolate.c new file mode 100644 index 0000000..392cb18 --- /dev/null +++ b/src/Interpolate.c @@ -0,0 +1,427 @@ + /*@@ + @file Interpolate.c + @date Wed 17 Jan 2001 + @author Thomas Radke + @desc + Interpolation of arrays to arbitrary points + + This interpolator is based on the Cactus 3.x Fortran version + written by Paul Walker. It also contains some nice optimization + features from Erik Schnetter. + @enddesc + @history + @date Wed 17 Jan 2001 + @author Thomas Radke + @hdesc Translation from Fortran to C + @endhistory + @version $Id$ + @@*/ + +#include + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "pughInterpGH.h" + +/* the rcs ID and its dummy function to use it */ +static char *rcsid = "$Id$"; +CCTK_FILEVERSION(CactusPUGH_PUGHInterp_Interpolate_c) + + +/* the highest order of interpolation we support so far */ +#define MAXORDER 3 + +/* the highest dimension for variables we can deal with (so far) */ +#define MAXDIM 3 + +/* macro to do the interpolation of in_array[n] to out_array[n] at point[] */ +#define INTERPOLATE(cctk_type, cctk_subtype, subpart, in_array, out_array, \ + order, point, dims, n, coeff) \ + { \ + int ii, jj, kk; \ + cctk_type *fi; \ + /* for some reason (probably loop unrolling & pipelining) the compiler \ + produces faster code if arrays are used instead of scalars */ \ + cctk_subtype interp_result, fj[1], fk[1]; \ + \ + \ + interp_result = 0; \ + \ + /* NOTE: support >3D arrays by adding more loops */ \ + for (kk = 0; kk <= order; kk++) \ + { \ + fk[0] = 0; \ + for (jj = 0; jj <= order; jj++) \ + { \ + /* NOTE: for >3D arrays adapt the index calculation here */ \ + fi = (cctk_type *) in_array + \ + point[0] + dims[0]*(point[1]+jj + dims[1]*(point[2]+kk)); \ + fj[0] = 0; \ + for (ii = 0; ii <= order; ii++) \ + { \ + fj[0] += fi[ii]subpart * coeff[0][ii]; \ + } \ + fk[0] += fj[0] * coeff[1][jj]; \ + } \ + interp_result += fk[0] * coeff[2][kk]; \ + } \ + \ + /* assign the result */ \ + ((cctk_type *) out_array)[n]subpart = interp_result; \ + } + +/* this is needed by some preprocessors to pass into INTERPOLATE + as a dummy macro */ +#define NOTHING + + +/*@@ + @routine PUGHInterp_Interpolate + @date Wed 17 Jan 2001 + @author Thomas Radke + @desc + This routine interpolates a set of input arrays + to a set of output arrays (one-to-one) at arbitrary points\ which are given by their coordinates and the underlying + regular, uniform grid. + + Current limitations of this implementation are: + + - arrays up to three (MAXDIM) dimensions only can be handled + - interpolation orders up to three (MAXORDER) only are supported + - coordinates must be given as CCTK_REAL types + - input and output array types must be the same + (no type casting of interpolation results supported) + + Despite of these limitations, the code was programmed almost + generic in that it can easily be extended to support higher- + dimensional arrays or more interpolation orders. + Please see the NOTES in this source file for details. + @enddesc + + @var num_points + @vdesc number of points to interpolate at + @vtype int + @vio in + @endvar + @var num_dims + @vdesc dimensionality of the input arrays + @vtype int + @vio in + @endvar + @var num_arrays + @vdesc number of input/output arrays + @vtype int + @vio in + @endvar + @var dims + @vdesc dimensions of the input arrays + @vtype int[ num_dims ] + @vio in + @endvar + @var coord + @vdesc coordinates to interpolate at + @vtype CCTK_REAL coord[ num_dims * num_points ] + @vio in + @endvar + @var origin + @vdesc origin of the underlying grid + @vtype CCTK_REAL origin[ num_dims ] + @vio in + @endvar + @var delta + @vdesc deltas of the underlying grid + @vtype CCTK_REAL delta[ num_dims ] + @vio in + @endvar + @var in_types + @vdesc CCTK variable types of input arrays + @vtype int in_types[ num_arrays ] + @vio in + @endvar + @var in_arrays + @vdesc list of input arrays + @vtype void *in_arrays[ num_arrays ] + @vio in + @endvar + @var out_types + @vdesc CCTK variable types of output arrays + @vtype int out_types[ num_arrays ] + @vio in + @endvar + @var out_arrays + @vdesc list of output arrays + @vtype void *out_arrays[ num_arrays ] + @vio out + @endvar + + @returntype int + @returndesc + 0 - successful interpolation + -1 - in case of any errors + @endreturndesc +@@*/ +int PUGHInterp_Interpolate (int order, + int num_points, + int num_dims, + int num_arrays, + int dims[ /* num_dims */ ], + CCTK_REAL coord[ /* num_dims * num_points */ ], + CCTK_REAL origin[ /* num_dims */ ], + CCTK_REAL delta[ /* num_dims */ ], + int in_types[ /* num_arrays */ ], + void *in_arrays[ /* num_arrays */ ], + int out_types[ /* num_arrays */ ], + void *out_arrays[ /* num_arrays */ ]) +{ + int retval; + int i, a, n, out_of_bounds, shift; + int max_dims[MAXDIM], point[MAXDIM]; + CCTK_REAL below[MAXDIM]; + CCTK_REAL offset[MAXDIM]; + CCTK_REAL coeff[MAXDIM][MAXORDER + 1]; + + + /* verify parameters and check against our restrictions */ + retval = -1; + if (num_dims < 1) + { + CCTK_WARN (1, "Number of dimensions must be positive"); + } + else if (num_dims > MAXDIM) + { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Interpolation of %d-dimensional arrays not implemented", + num_dims); + } + else if (order < 1) + { + CCTK_WARN (1, "Inperpolation order must be positive"); + } + else if (order > MAXORDER) + { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Interpolation order %d not implemented", order); + } + else if (num_points < 0) + { + CCTK_WARN (1, "Negative number of points given"); + } + else + { + retval = 0; + } + + /* immediately return in case of errors */ + if (retval) + { + return (retval); + } + + /* also immediately return if there's nothing to do */ + if (num_points == 0) + { + return (retval); + } + + /* duplicate the dims[] vector into one with MAXDIM-1 elements + (with the remaining elements zeroed out) + so that we can use nested loops over MAXDIM dimensions later on */ + memset (max_dims, 0, sizeof (max_dims)); + memcpy (max_dims, dims, (num_dims - 1) * sizeof (int)); + + /* zero out the coefficients and set the elements with index 'order' to one + so that we can use nested loops over MAXDIM dimensions later on */ + memset (coeff, 0, sizeof (coeff)); + for (i = num_dims; i < MAXDIM; i++) + { + coeff[i][0] = 1; + } + + /* zero out the iterator */ + memset (point, 0, sizeof (point)); + + /* loop over all points to interpolate at */ + for (n = 0; n < num_points; n++) + { + /* reset the out-of-bounds flag */ + out_of_bounds = 0; + + /* loop over all dimensions */ + for (i = 0; i < num_dims; i++) + { + /* nearest grid point for stencil */ + point[i] = (coord[num_dims*n + i] - origin[i]) / delta[i]; + + /* if beyond upper bound shift the grid point to the left */ + shift = point[i] + order - (dims[i] - 1); + if (shift > 0) + { + point[i] -= shift; + } + + /* physical coordinate of that grid point */ + below[i] = origin[i] + point[i] * delta[i]; + + /* offset from that grid point, in fractions of grid points */ + offset[i] = (coord[num_dims*n + i] - below[i]) / delta[i]; + + /* test bounds */ + out_of_bounds |= point[i] < 0 || point[i]+order >= dims[i]; + } + + /* check bounds */ + if (out_of_bounds) + { + char *msg; + + + /* put all information into a single message string for output */ + msg = (char *) malloc (100 + num_dims*(10 + 4*20)); + sprintf (msg, "Interpolation stencil out of bounds at grid point [%d", + point[0]); + for (i = 1; i < num_dims; i++) + { + sprintf (msg, "%s, %d", msg, point[i]); + } + sprintf (msg, "%s]\nrange would be min/max [%f / %f", msg, + (double) below[0], (double) (below[0] + offset[0])); + for (i = 1; i < num_dims; i++) + { + sprintf (msg, "%s, %f / %f", msg, + (double) below[i], (double) (below[i] + offset[i])); + } + sprintf (msg, "%s]\ngrid is min/max [%f / %f", msg, + (double) origin[0], (double) (origin[0] + dims[0]*delta[0])); + for (i = 1; i < num_dims; i++) + { + sprintf (msg, "%s, %f / %f", msg, + (double) origin[i], (double) (origin[i] + dims[i]*delta[i])); + } + sprintf (msg, "%s]", msg); + CCTK_WARN (1, msg); + free (msg); + + continue; + } + + /* compute the interpolation coefficients according to the order + Thanks to Erik for formulating these so nicely. */ + /* NOTE: support higher interpolation orders by adding the + appropriate coefficients in another else branch */ + if (order == 1) + { + /* first order (linear) 1D interpolation */ + for (i = 0; i < num_dims; i++) + { + coeff[i][0] = 1 - offset[i]; + coeff[i][1] = offset[i]; + } + } + else if (order == 2) + { + /* second order (quadratic) 1D interpolation */ + for (i = 0; i < num_dims; i++) + { + coeff[i][0] = (1-offset[i]) * (2-offset[i]) / ( 2 * 1 ); + coeff[i][1] = ( -offset[i]) * (2-offset[i]) / ( 1 * (-1)); + coeff[i][2] = ( -offset[i]) * (1-offset[i]) / ((-1) * (-2)); + } + } + else if (order == 3) + { + /* third order (cubic) 1D interpolation */ + for (i = 0; i < num_dims; i++) + { + coeff[i][0] = (1-offset[i]) * (2-offset[i]) * (3-offset[i]) / + ( 3 * 2 * 1 ); + coeff[i][1] = ( -offset[i]) * (2-offset[i]) * (3-offset[i]) / + ( 2 * 1 * (-1)); + coeff[i][2] = ( -offset[i]) * (1-offset[i]) * (3-offset[i]) / + ( 1 * (-1) * (-2)); + coeff[i][3] = ( -offset[i]) * (1-offset[i]) * (2-offset[i]) / + ((-1) * (-2) * (-3)); + } + } + else + { + CCTK_WARN (0, "Implementation error"); + } + + /* now loop over all arrays to interpolate at the current point */ + for (a = 0; a < num_arrays; a++) + { + /* sorry, for now input and output arrays must be of same type */ + if (in_types[a] != out_types[a]) + { + CCTK_WARN (1, "Type casting of interpolation results not implemented"); + continue; + } + + /* now do the interpolation according to the array type + we support all kinds of CCTK_REAL* and CCTK_COMPLEX* types here */ + if (in_types[a] == CCTK_VARIABLE_REAL) + { + INTERPOLATE (CCTK_REAL, CCTK_REAL, NOTHING, in_arrays[a], + out_arrays[a], order, point, max_dims, n, coeff); + } + else if (in_types[a] == CCTK_VARIABLE_COMPLEX) + { + INTERPOLATE (CCTK_COMPLEX, CCTK_REAL, .Re, in_arrays[a], + out_arrays[a], order, point, max_dims, n, coeff); + INTERPOLATE (CCTK_COMPLEX, CCTK_REAL, .Im, in_arrays[a], + out_arrays[a], order, point, max_dims, n, coeff); + } +#ifdef CCTK_REAL4 + else if (in_types[a] == CCTK_VARIABLE_REAL4) + { + INTERPOLATE (CCTK_REAL4, CCTK_REAL4, NOTHING, in_arrays[a], + out_arrays[a], order, point, max_dims, n, coeff); + } + else if (in_types[a] == CCTK_VARIABLE_COMPLEX8) + { + INTERPOLATE (CCTK_COMPLEX8, CCTK_REAL4, .Re, in_arrays[a], + out_arrays[a], order, point, max_dims, n, coeff); + INTERPOLATE (CCTK_COMPLEX8, CCTK_REAL4, .Im, in_arrays[a], + out_arrays[a], order, point, max_dims, n, coeff); + } +#endif +#ifdef CCTK_REAL8 + else if (in_types[a] == CCTK_VARIABLE_REAL8) + { + INTERPOLATE (CCTK_REAL8, CCTK_REAL8, NOTHING, in_arrays[a], + out_arrays[a], order, point, max_dims, n, coeff); + } + else if (in_types[a] == CCTK_VARIABLE_COMPLEX16) + { + INTERPOLATE (CCTK_COMPLEX16, CCTK_REAL8, .Re, in_arrays[a], + out_arrays[a], order, point, max_dims, n, coeff); + INTERPOLATE (CCTK_COMPLEX16, CCTK_REAL8, .Im, in_arrays[a], + out_arrays[a], order, point, max_dims, n, coeff); + } +#endif +#ifdef CCTK_REAL16 + else if (in_types[a] == CCTK_VARIABLE_REAL16) + { + INTERPOLATE (CCTK_REAL16, CCTK_REAL16, NOTHING, in_arrays[a], + out_arrays[a], order, point, max_dims, n, coeff); + } + else if (in_types[a] == CCTK_VARIABLE_COMPLEX32) + { + INTERPOLATE (CCTK_COMPLEX32, CCTK_REAL16, .Re, in_arrays[a], + out_arrays[a], order, point, max_dims, n, coeff); + INTERPOLATE (CCTK_COMPLEX32, CCTK_REAL16, .Im, in_arrays[a], + out_arrays[a], order, point, max_dims, n, coeff); + } +#endif + else + { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Unsupported variable type %d", (int) in_types[a]); + } + } /* end of loop over all arrays */ + + } /* end of loop over all points to interpolate at */ + + /* we're done */ + return (retval); +} diff --git a/src/Operator.c b/src/Operator.c new file mode 100644 index 0000000..76a71f5 --- /dev/null +++ b/src/Operator.c @@ -0,0 +1,1293 @@ + /*@@ + @file Operator.c + @date Tue Apr 15 18:22:45 1997 + @author Paul Walker + @desc + Definition of interpolation operators for regular uniform grids. + @enddesc + + @history + @date Sun Jul 04 1999 + @author Thomas Radke + @hdesc conversion to Cactus 4.0 (copied from pughGetPoints.c) + @date Wed 31 Jan 2001 + @author Thomas Radke + @hdesc translation of fortran interpolators into C + @endhistory + @version $Id$ + @@*/ + +#include +#include +#include + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "CactusPUGH/PUGH/src/include/pugh.h" +#include "pughInterpGH.h" + +/* the rcs ID and its dummy function to use it */ +static char *rcsid = "$Header$"; +CCTK_FILEVERSION(CactusPUGH_PUGHInterp_Operator_c) + +/* uncomment this to get some debugging output */ +/* #define PUGHINTERP_DEBUG 1 */ + +/* macro do sort interpolation results from a single communication buffer + into their appropriate output arrays */ +#define SORT_TYPED_ARRAY(cctk_type, source, dest, points, index_array, offset)\ + { \ + int i; \ + cctk_type *src, *dst; \ + \ + \ + src = (cctk_type *) source; \ + dst = (cctk_type *) dest; \ + for (i = 0; i < points; i++) \ + { \ + dst[index_array[i + offset]] = *src++; \ + } \ + source = (char *) src; \ + } + + + +#ifdef CCTK_MPI +/* internal structure describing a handle for a single CCTK data type */ +typedef struct +{ + int vtypesize; /* variable type's size in bytes */ + MPI_Datatype mpi_type; /* corresponding MPI datatype */ + int num_arrays; /* number of in/out arrays */ + int array; /* iterator for num_arrays */ + void *send_buffer; /* communication send buffer for this type */ + void *recv_buffer; /* communication receive buffer for this type */ + char *buffer; /* work pointer for send_buffer */ +} type_desc_t; +#endif + + +/* prototypes of routines defined in this source file */ +static int PUGHInterp_CheckArguments (cGH *GH, + int num_dims, + int num_points, + int num_in_arrays, + int num_out_arrays, + int interp_coord_array_types[]); +#ifdef CCTK_MPI +static int GetLocalCoords (cGH *GH, + int num_points, + const char *coord_system, + pGExtras *extras, + CCTK_REAL *coords[], + int *num_local_points, + CCTK_REAL **local_coords); +#endif + + +/*@@ + @routine PUGHInterp_InterpGV + @date Sun Jul 04 1999 + @author Thomas Radke + @desc + The interpolation operator registered with the CCTK + under the name "regular uniform cartesian". + + Interpolates a list of CCTK variables (domain-decomposed + grid functions or arrays) to a list of output arrays + (one-to-one) at a given number of interpolation points + (indicated by their coordinates). The points are located + on a coordinate system which is assumed to be a uniform + cartesian. + @enddesc + + @var GH + @vdesc Pointer to CCTK grid hierarchy + @vtype cGH * + @vio in + @endvar + @var order + @vdesc interpolation order + @vtype int + @vio in + @endvar + @var coord_system + @vdesc name of coordinate system to use for interpolation + @vtype const char * + @vio in + @endvar + @var num_points + @vdesc number of points to be interpolated on this processor + @vtype int + @vio in + @endvar + @var num_in_array_indices + @vdesc number of input arrays (given by their indices) + to interpolate from + @vtype int + @vio in + @endvar + @var num_out_arrays + @vdesc number of output arrays to interpolate to + @vtype int + @vio in + @endvar + @var interp_coord_arrays + @vdesc coordinates of points to interpolate at + @vtype void *[num_dims] + @vio in + @endvar + @var interp_coord_array_types + @vdesc CCTK data type of coordinate arrays + @vtype int [num_dims] + @vio in + @endvar + @var in_arrays + @vdesc list of input arrays to interpolate on + @vtype void *[num_in_arrays] + @vio in + @endvar + @var in_array_types + @vdesc CCTK data types of input arrays + @vtype int [num_in_arrays] + @vio in + @endvar + @var out_arrays + @vdesc list of output arrays to interpolate to + @vtype void *[num_out_arrays] + @vio out + @endvar + @var out_array_types + @vdesc CCTK data types of output arrays + @vtype int [num_out_arrays] + @vio in + @endvar + + @returntype int + @returndesc + 0 - successful interpolation + -1 - in case of any errors + @endreturndesc +@@*/ +int PUGHInterp_InterpGV (cGH *GH, + int order, + const char *coord_system, + int num_points, + int num_in_array_indices, + int num_out_arrays, + void *interp_coord_arrays[], + int interp_coord_array_types[], + int in_array_indices[], + void *out_arrays[], + int out_array_types[]) +{ + CCTK_REAL **data, *interp_local_coords; + CCTK_REAL *origin, *delta; + int nprocs; + int timelevel; + int dim, num_dims, *dims; + int array; + int point; + int retval; + void **in_arrays; + pGH *pughGH; + pGExtras *extras; + cGroupDynamicData group_data; +#ifdef CCTK_MPI + int num_local_points; + int offset; + int proc; + int type; + void **local_out_arrays; + pughInterpGH *myGH; + int max_type; + type_desc_t *type_desc; +#endif + + + /* get dimensionality of the coordinate system */ + num_dims = CCTK_CoordSystemDim (coord_system); + if (num_dims <= 0) + { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Cannot get dimensions of coordinate system '%s'", + coord_system); + return (-1); + } + + /* check other arguments */ + retval = PUGHInterp_CheckArguments (GH, num_dims, num_points, + num_in_array_indices, num_out_arrays, + interp_coord_array_types); + if (retval <= 0) + { + return (retval); + } + + /* get dimensions, origin, and delta of the processor-local grid */ + /* NOTE: getting the dimensions should be a flesh routine as well + for now we get the dimensions of every coordinate and take the + i'th element - this is inconsistent !! */ + dims = (int *) malloc (2 * num_dims * sizeof (int)); + origin = (CCTK_REAL *) malloc (2 * num_dims * sizeof (CCTK_REAL)); + delta = origin + num_dims; + in_arrays = malloc (num_in_array_indices * sizeof (void *)); + for (dim = 0; dim < num_dims; dim++) + { + CCTK_GrouplshVI (GH, num_dims, dims + num_dims, + CCTK_CoordIndex (dim + 1, NULL, coord_system)); + dims[dim] = dims[dim + num_dims]; + + CCTK_CoordLocalRange (GH, &origin[dim], &delta[dim], dim + 1, NULL, + coord_system); + delta[dim] = (delta[dim] - origin[dim]) / dims[dim]; + } + + /* get extension handle for PUGH */ + pughGH = (pGH *) CCTK_GHExtension (GH, "PUGH"); + extras = NULL; + + /* check that the input arrays dimensions match the coordinate system + (but their dimensionality can be less) */ + for (array = 0; array < num_in_array_indices; array++) + { + if (CCTK_GroupDynamicData (GH, + CCTK_GroupIndexFromVarI(in_array_indices[array]), + &group_data) < 0) + { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Invalid input array index %d", + in_array_indices[array]); + retval = -1; + continue; + } + + if (group_data.dim > num_dims) + { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Input array variable with index %d has more dimensions " + "than coordinate system '%s'", + in_array_indices[array], coord_system); + retval = -1; + continue; + } + + if (memcmp (group_data.lsh, dims, group_data.dim * sizeof (int))) + { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Dimensions of input array variable with index %d " + "doesn't match with coordinate system '%s'", + in_array_indices[array], coord_system); + retval = -1; + } + + /* the extras pointer will refer to the first input array with the + highest dimension (need this information later on) */ + if (! extras || extras->dim < group_data.dim) + { + extras = ((pGA *) pughGH->variables[in_array_indices[array]][0])->extras; + } + + /* get the data pointer to the input array (use current timelevel) */ + timelevel = CCTK_NumTimeLevelsFromVarI (in_array_indices[array]) - 1; + if (timelevel > 0) + { + timelevel--; + } + in_arrays[array] = CCTK_VarDataPtrI(GH, timelevel, in_array_indices[array]); + } + if (retval < 0) + { + free (in_arrays); + free (origin); + free (dims); + return (retval); + } + + /* single-processor case is easy: no communication or buffering, just direct + interpolation of interp_coord_arrays from in_arrays into out_arrays */ + nprocs = CCTK_nProcs (GH); + if (nprocs == 1) + { + /* sort the individual interpolation coordinate arrays into a single one */ + interp_local_coords = (CCTK_REAL *) malloc (num_dims * num_points * + sizeof (CCTK_REAL)); + data = (CCTK_REAL **) interp_coord_arrays; + for (point = 0; point < num_points; point++) + { + for (dim = 0; dim < num_dims; dim++) + { + *interp_local_coords++ = data[dim][point]; + } + } + interp_local_coords -= num_dims * num_points; + + /* call the interpolator function */ + retval = PUGHInterp_Interpolate (order, + num_points, num_dims, num_out_arrays, + dims, interp_local_coords, + origin, delta, + out_array_types, in_arrays, + out_array_types, out_arrays); + + /* free allocated resources */ + free (interp_local_coords); + free (in_arrays); + free (origin); + free (dims); + + return (retval); + } + +#ifdef CCTK_MPI + /*** Here follows the multi-processor case: + All processors locate their points to interpolate at + and exchange the coordinates so that every processor gets + those points which it can process locally. + After interpolation the results have to be send back to the + requesting processors. + For both communications MPI_Alltoallv() is used. + + In order to minimize the total number of MPI_Alltoallv() calls + (which are quite expensive) we collect the interpolation results + for all output arrays of the same CCTK data type into a single + communication buffer. That means, after communication the data + needs to be resorted from the buffer into the output arrays. + ***/ + + /* first of all, set up a structure with information of the + CCTK data types we have to deal with */ + + /* get the maximum value of the output array CCTK data types + NOTE: we assume that CCTK data types are defined as consecutive + positive constants starting from zero */ + max_type = 0; + for (array = 0; array < num_out_arrays; array++) + { + if (max_type < out_array_types[array]) + { + max_type = out_array_types[array]; + } + } + + /* now allocate an array of structures for all potential types */ + type_desc = (type_desc_t *) calloc (max_type + 1, sizeof (type_desc_t)); + + /* count the number of arrays of same type + (the num_arrays element was already initialized to zero by calloc() */ + for (array = 0; array < num_out_arrays; array++) + { + type_desc[out_array_types[array]].num_arrays++; + } + + /* fill in the type description information */ + retval = 0; + for (type = 0; type <= max_type; type++) + { + if (type_desc[type].num_arrays > 0) + { + /* get the variable type size in bytes */ + type_desc[type].vtypesize = CCTK_VarTypeSize (type); + if (type_desc[type].vtypesize <= 0) + { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Invalid variable type %d passed, " + "arrays of such type will be skipped during interpolation", + type); + type_desc[type].num_arrays = 0; + continue; + } + + /* get the MPI data type to use for communicating such a CCTK data type */ + type_desc[type].mpi_type = PUGH_MPIDataType (pughGH, type); + if (type_desc[type].mpi_type == MPI_DATATYPE_NULL) + { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "No MPI data type defined for variable type %d, " + "arrays of such type will be skipped during interpolation", + type); + type_desc[type].num_arrays = 0; + continue; + } + + retval++; + } + } + + /* check that there's at least one array with a valid CCTK data type */ + if (retval <= 0) + { + free (in_arrays); + free (dims); + free (origin); + free (type_desc); + return (-1); + } + + /* map the requested points to interpolate at onto the processors + they belong to and gather the coordinates of all the points + that this processor owns + the number of processor-local points is returned in num_local_points, + their coordinates in interp_local_coords */ + retval = GetLocalCoords (GH, num_points, coord_system, extras, + (CCTK_REAL **) interp_coord_arrays, + &num_local_points, &interp_local_coords); + if (retval) + { + free (in_arrays); + free (origin); + free (dims); + free (type_desc); + return (retval); + } + + /* allocate contiguous communication buffers for each CCTK data type + holding the local interpolation results from all input arrays + of that type + If there are no points to process on this processor + set the buffer pointer to an invalid but non-NULL value + otherwise we might get trouble with NULL pointers in MPI_Alltoallv () */ + for (type = 0; type <= max_type; type++) + { + if (type_desc[type].num_arrays > 0 && num_local_points > 0) + { + type_desc[type].send_buffer = malloc (num_local_points * + type_desc[type].num_arrays * + type_desc[type].vtypesize); + type_desc[type].buffer = (char *) type_desc[type].send_buffer; + type_desc[type].array = 0; + } + else + { + /* dereferencing such an address should code crash on most systems */ + type_desc[type].send_buffer = (void *) 1; + } + } + + /* get extension handle for interp */ + myGH = (pughInterpGH *) CCTK_GHExtension (GH, "PUGHInterp"); + + /* allocate new out_arrays array for local interpolation results + from this processor */ + local_out_arrays = calloc (num_out_arrays, sizeof (void *)); + + /* now, in a loop over all processors, do the interpolation + and put the results in the communication buffer at the proper offset */ + for (proc = 0; proc < nprocs; proc++) + { + for (type = 0; type <= max_type; type++) + { + if (type_desc[type].num_arrays > 0) + { + for (array = 0; array < num_out_arrays; array++) + { + if (out_array_types[array] == type) + { + local_out_arrays[array] = type_desc[type].buffer; + type_desc[type].buffer += myGH->num_points_to[proc] * + type_desc[type].vtypesize; + } + } + } + } + + /* call the interpolation operator to process all points of all + output arrays for this processor */ + PUGHInterp_Interpolate (order, + myGH->num_points_to[proc], num_dims, num_out_arrays, + dims, interp_local_coords, origin, delta, + out_array_types, in_arrays, + out_array_types, local_out_arrays); + + /* have to add offset for this processor to coordinates array */ + interp_local_coords += myGH->num_points_to[proc] * num_dims; + + } /* end of loop over all processors */ + + /* don't need these anymore */ + if (num_local_points > 0) + { + interp_local_coords -= num_local_points * num_dims; + free (interp_local_coords); + } + free (local_out_arrays); + free (in_arrays); + free (origin); + free (dims); + + /* now send the interpolation results back to the processors they were + requested from, also receive my own results that were computed + by other processors + Since all the locally computed results are in a single contiguous buffer + we need to call MPI_Alltoall() only once for each CCTK data type. */ + for (type = 0; type <= max_type; type++) + { + /* skip unused types */ + if (type_desc[type].num_arrays <= 0) + { + continue; + } + + /* set up the communication (this is type-independent) */ + myGH->send_count[0] = type_desc[type].num_arrays * myGH->num_points_to[0]; + myGH->recv_count[0] = type_desc[type].num_arrays * myGH->num_points_from[0]; + myGH->send_displ[0] = myGH->recv_displ[0] = 0; + for (proc = 1; proc < pughGH->nprocs; proc++) + { + myGH->send_count[proc] = type_desc[type].num_arrays * + myGH->num_points_to[proc]; + myGH->recv_count[proc] = type_desc[type].num_arrays * + myGH->num_points_from[proc]; + myGH->send_displ[proc] = myGH->send_displ[proc-1] + + myGH->send_count[proc-1]; + myGH->recv_displ[proc] = myGH->recv_displ[proc-1] + + myGH->recv_count[proc-1]; + } + + /* allocate buffer for receiving my own requested points */ + /* avoid NULL pointers here because MPI_Alltoallv() doesn't like it */ + if (num_points > 0) + { + type_desc[type].recv_buffer = malloc (num_points * + type_desc[type].num_arrays * + type_desc[type].vtypesize); + } + else + { + /* access to such a fake address should crash the code on most systems */ + type_desc[type].recv_buffer = (void *) 1; + } + + /* now exchange the data for this CCTK data type */ + CACTUS_MPI_ERROR (MPI_Alltoallv (type_desc[type].send_buffer, + myGH->send_count, myGH->send_displ, + type_desc[type].mpi_type, + type_desc[type].recv_buffer, + myGH->recv_count, myGH->recv_displ, + type_desc[type].mpi_type, + pughGH->PUGH_COMM_WORLD)); + + /* now that the data is sent we don't need the buffer anymore */ + if (num_local_points > 0) + { + free (type_desc[type].send_buffer); + } + + /* no sort neccessary if there are no points */ + if (num_points <= 0) + { + continue; + } + + /* Go back from processor-sorted data to input-ordered data. + The creation of the indices array above makes this not so bad. */ + offset = 0; + type_desc[type].buffer = (char *) type_desc[type].recv_buffer; + for (proc = 0; proc < nprocs; proc++) + { + for (array = 0; array < num_out_arrays; array++) + { + if (out_array_types[array] != type) + { + continue; + } + + /* now do the sorting according to the CCTK data type */ + if (out_array_types[array] == CCTK_VARIABLE_CHAR) + { + SORT_TYPED_ARRAY (CCTK_CHAR, type_desc[type].buffer, + out_arrays[array], myGH->num_points_from[proc], + myGH->indices, offset); + } + else if (out_array_types[array] == CCTK_VARIABLE_INT) + { + SORT_TYPED_ARRAY (CCTK_INT, type_desc[type].buffer, + out_arrays[array], myGH->num_points_from[proc], + myGH->indices, offset); + } + else if (out_array_types[array] == CCTK_VARIABLE_REAL) + { + SORT_TYPED_ARRAY (CCTK_REAL, type_desc[type].buffer, + out_arrays[array], myGH->num_points_from[proc], + myGH->indices, offset); + } + else if (out_array_types[array] == CCTK_VARIABLE_COMPLEX) + { + SORT_TYPED_ARRAY (CCTK_COMPLEX, type_desc[type].buffer, + out_arrays[array], myGH->num_points_from[proc], + myGH->indices, offset); + } +#ifdef CCTK_REAL4 + else if (out_array_types[array] == CCTK_VARIABLE_REAL4) + { + SORT_TYPED_ARRAY (CCTK_REAL4, type_desc[type].buffer, + out_arrays[array], myGH->num_points_from[proc], + myGH->indices, offset); + } + else if (out_array_types[array] == CCTK_VARIABLE_COMPLEX8) + { + SORT_TYPED_ARRAY (CCTK_COMPLEX8, type_desc[type].buffer, + out_arrays[array], myGH->num_points_from[proc], + myGH->indices, offset); + } +#endif +#ifdef CCTK_REAL8 + else if (out_array_types[array] == CCTK_VARIABLE_REAL8) + { + SORT_TYPED_ARRAY (CCTK_REAL8, type_desc[type].buffer, + out_arrays[array], myGH->num_points_from[proc], + myGH->indices, offset); + } + else if (out_array_types[array] == CCTK_VARIABLE_COMPLEX16) + { + SORT_TYPED_ARRAY (CCTK_COMPLEX16, type_desc[type].buffer, + out_arrays[array], myGH->num_points_from[proc], + myGH->indices, offset); + } +#endif +#ifdef CCTK_REAL16 + else if (out_array_types[array] == CCTK_VARIABLE_REAL16) + { + SORT_TYPED_ARRAY (CCTK_REAL16, type_desc[type].buffer, + out_arrays[array], myGH->num_points_from[proc], + myGH->indices, offset); + } + else if (out_array_types[array] == CCTK_VARIABLE_COMPLEX32) + { + SORT_TYPED_ARRAY (CCTK_COMPLEX32, type_desc[type].buffer, + out_arrays[array], myGH->num_points_from[proc], + myGH->indices, offset); + } +#endif + else + { + CCTK_WARN (0, "Implementation error"); + } + + } /* end of loop over all output arrays */ + + /* advance the offset into the communication receive buffer */ + offset += myGH->num_points_from[proc]; + + } /* end of loop over all processors */ + + /* this communication receive buffer isn't needed anymore */ + free (type_desc[type].recv_buffer); + + } /* end of loop over all types */ + + /* free remaining resources allocated within this run */ + if (myGH->whichproc) + { + free (myGH->whichproc); + myGH->whichproc = NULL; + } + free (type_desc); +#endif /* MPI */ + + return (0); +} + + +/*@@ + @routine PUGHInterp_InterpLocal + @date Sun Jul 04 1999 + @author Thomas Radke + @desc + The interpolation operator registered with the CCTK + under the name "uniform cartesian". + + Interpolates a list of non-distributed, processor-local + input arrays to a list of output arrays (one-to-one) + at a given number of interpolation points (indicated by + their coordinates). The points are located on a coordinate + system which is assumed to be a uniform cartesian. + @enddesc + + @var GH + @vdesc Pointer to CCTK grid hierarchy + @vtype cGH * + @vio in + @endvar + @var order + @vdesc interpolation order + @vtype int + @vio in + @endvar + @var num_points + @vdesc number of points to be interpolated on this processor + @vtype int + @vio in + @endvar + @var num_dims + @vdesc dimensionality of the underlying grid + @vtype int + @vio in + @endvar + @var num_in_arrays + @vdesc number of input arrays to interpolate from + @vtype int + @vio in + @endvar + @var num_out_arrays + @vdesc number of output arrays to interpolate to + @vtype int + @vio in + @endvar + @var coord_dims + @vdesc dimensions of the underlying grid + @vtype int [num_dims] + @vio in + @endvar + @var coord_arrays + @vdesc list of grid coordinate arrays + @vtype void *[num_dims] + @vio in + @endvar + @var coord_array_types + @vdesc CCTK data type of coordinate arrays + @vtype int [num_dims] + @vio int + @endvar + @var interp_coord_arrays + @vdesc coordinates of points to interpolate at + @vtype void *[num_dims] + @vio in + @endvar + @var interp_coord_array_types + @vdesc CCTK data type of coordinate arrays to interpolate at + @vtype int [num_dims] + @vio out + @endvar + @var in_arrays + @vdesc list of input arrays to interpolate on + @vtype void *[num_in_arrays] + @vio in + @endvar + @var in_array_types + @vdesc CCTK data types of input arrays + @vtype int [num_in_arrays] + @vio in + @endvar + @var out_arrays + @vdesc list of output arrays to interpolate to + @vtype void *[num_out_arrays] + @vio out + @endvar + @var out_array_types + @vdesc CCTK data types of output arrays + @vtype int [num_out_arrays] + @vio in + @endvar + + @returntype int + @returndesc + 0 - successful interpolation + -1 - in case of any errors + @endreturndesc +@@*/ +int PUGHInterp_InterpLocal (cGH *GH, + int order, + int num_points, + int num_dims, + int num_in_arrays, + int num_out_arrays, + int coord_dims[], + void *coord_arrays[], + int coord_array_types[], + void *interp_coord_arrays[], + int interp_coord_array_types[], + void *in_arrays[], + int in_array_types[], + void *out_arrays[], + int out_array_types[]) +{ + int dim, point, retval; + CCTK_REAL *coords, *origin, *delta, **data; + + + /* check arguments */ + retval = PUGHInterp_CheckArguments (GH, num_dims, num_points, + num_in_arrays, num_out_arrays, + interp_coord_array_types); + if (retval <= 0) + { + return (retval); + } + for (dim = 0; dim < num_dims; dim++) + { + if (coord_array_types[dim] != CCTK_VARIABLE_REAL) + { + CCTK_WARN (1, "Coordinates should be of type CCTK_REAL"); + return (-1); + } + if (interp_coord_array_types[dim] != CCTK_VARIABLE_REAL) + { + CCTK_WARN (1, "Interpolation coordinates should be of type CCTK_REAL"); + return (-1); + } + } + + /* get the grid spacings - this assumes a cartesian grid */ + origin = (CCTK_REAL *) malloc (2 * num_dims * sizeof (CCTK_REAL)); + delta = origin + num_dims; + data = (CCTK_REAL **) coord_arrays; + for (dim = 0; dim < num_dims; dim++) + { + origin[dim] = data[dim][0]; + delta[dim] = data[dim][1] - data[dim][0]; + } + + /* sort the individual interpolation coordinate arrays into a single one */ + coords = (CCTK_REAL *) malloc (num_dims * num_points * sizeof (CCTK_REAL)); + data = (CCTK_REAL **) interp_coord_arrays; + for (point = 0; point < num_points; point++) + { + for (dim = 0; dim < num_dims; dim++) + { + *coords++ = data[dim][point]; + } + } + coords -= num_dims * num_points; + + /* call the interpolator function */ + retval = PUGHInterp_Interpolate (order, + num_points, num_dims, num_out_arrays, + coord_dims, coords, origin, delta, + in_array_types, in_arrays, + out_array_types, out_arrays); + + /* free allocated resources */ + free (coords); + free (origin); + + return (retval); +} + + +/**************************************************************************/ +/* local routines */ +/**************************************************************************/ + +#ifdef CCTK_MPI + +/*@@ + @routine GetLocalCoords + @date Sun Jul 04 1999 + @author Thomas Radke + @desc + Collect the coordinates of all points to be processed locally + into an array coords[num_dims][num_local_points]. + + This means for the single-processor case to sort + + inCoords1[num_points], inCoords2[npoints], ..., inCoords[npoints] + into coords[num_dims][num_points] + + where num_points == num_local_points. + + In the multiprocessor case all processors map their points' coordinates + to the processor that owns this point and exchange this information + via MPI_Alltoall (). + num_local_points is then the number of all processors' points to be + interpolated locally on this processor. + + This routine returns the number of points to be processed locally and + a pointer to the allocated array of their coordinates. + @enddesc + + @var GH + @vdesc Pointer to CCTK grid hierarchy + @vtype cGH + @vio in + @endvar + @var num_points + @vdesc number of points to be interpolated on this processor + @vtype int + @vio in + @endvar + @var isGlobal + @vdesc flag indicating that coordinates are global (and need to be + collected over all processors) + @vtype int + @vio in + @endvar + @var num_dims + @vdesc number of coordinate dimensions for each point + @vtype int + @vio in + @endvar + @var coords + @vdesc coordinates of each point to be interpolated on this processor + @vtype CCTK_REAL array of size num_dims + @vio in + @endvar + @var num_local_points + @vdesc number of points to be processed by this processor + @vtype int * + @vio out + @endvar + @var local_coords + @vdesc coordinates of each point to be processed by this processor + @vtype pointer to CCTK_REAL array of dims[num_dims][num_local_points] + @vio out + @endvar + +@@*/ +static int GetLocalCoords (cGH *GH, + int num_points, + const char *coord_system, + pGExtras *extras, + CCTK_REAL *coords[], + int *num_local_points, + CCTK_REAL **local_coords) +{ + DECLARE_CCTK_PARAMETERS + int dim, point; + int tmp, proc, nprocs; + pGH *pughGH; + pughInterpGH *myGH; + CCTK_REAL *range_min, *range_max; + CCTK_REAL *origin, *delta; + CCTK_REAL *proc_coords; +#define FUDGE 0.00001 + + + /* Do a fliperooni in octant mode */ + /* This flips the sign around the origin only */ + /* if the left side (cx0) is bigger than the position to interpolate at */ + /* BoxInBox needs this for interpolating onto grids */ + if (CCTK_Equals (domain, "octant")) + { + for (point = 0; point < num_points; point++) + { + tmp = 0; + for (dim = 0; dim < extras->dim; dim++) + { + if (coords[dim][point] <= GH->cctk_origin_space[dim]) + { + tmp = 1; /* mark flipping */ + break; + } + } +/*** FIXME *** really change input arrays here ??? */ + if (tmp) + { + for (dim = 0; dim < extras->dim; dim++) + { + coords[dim][point] = fabs (coords[dim][point]); + } + } + } + } + + /* get GH extension handles for PUGHInterp and PUGH */ + myGH = (pughInterpGH *) CCTK_GHExtension (GH, "PUGHInterp"); + pughGH = (pGH *) CCTK_GHExtension (GH, "PUGH"); + + /* This holds the proccessor for *each* of num_points points */ + if (num_points > 0) + { + myGH->whichproc = (int *) malloc (2 * num_points * sizeof (int)); + } + else + { + myGH->whichproc = NULL; + } + /* indices[] is used to make the sorting easier + when receiving the output data */ + myGH->indices = myGH->whichproc + num_points; + + /* initialize whichproc with invalid processor number -1 */ + for (point = 0; point < num_points; point++) + { + myGH->whichproc[point] = -1; + } + + /* initialize num_points_from to 0 for counting it up in the following loop */ + nprocs = CCTK_nProcs (GH); + memset (myGH->num_points_from, 0, nprocs * sizeof (CCTK_INT)); + + /* allocate the ranges for my local coordinates */ + range_min = (CCTK_REAL *) malloc (4 * extras->dim * sizeof (CCTK_REAL)); + range_max = range_min + 1*extras->dim; + origin = range_min + 2*extras->dim; + delta = range_min + 3*extras->dim; + + /* get the global origin and delta of the coordinate system */ + for (dim = 0; dim < extras->dim; dim++) + { + CCTK_CoordRange (GH, &origin[dim], &delta[dim], dim+1, NULL, coord_system); + delta[dim] = (delta[dim] - origin[dim]) / extras->nsize[dim]; + } + + /* locate the points to interpolate at */ + for (proc = 0; proc < nprocs; proc++) + { + for (dim = 0; dim < extras->dim; dim++) + { + /* compute the coordinate ranges */ + range_min[dim] = origin[dim] + (extras->lb[proc][dim] - FUDGE)*delta[dim]; + range_max[dim] = origin[dim] + (extras->ub[proc][dim] + FUDGE)*delta[dim]; + } + + /* and now which point will be processed by what processor */ + for (point = 0; point < num_points; point++) + { + /* skip points which have already been located */ + if (myGH->whichproc[point] >= 0) + { + continue; + } + + /* check whether the point belongs to this processor + (must be within min/max in all dimensions) */ + tmp = 0; + for (dim = 0; dim < extras->dim; dim++) + { + if (coords[dim][point] >= range_min[dim] && + coords[dim][point] <= range_max[dim]) + { + tmp++; + } + } + if (tmp == extras->dim) + { + myGH->whichproc[point] = proc; + myGH->num_points_from[proc]++; + } + } + } + /* don't need this anymore */ + free (range_min); + + /* make sure that all points could be mapped onto a processor */ + tmp = 0; + for (point = 0; point < num_points; point++) + { + if (myGH->whichproc[point] < 0) + { + int i; + char *msg = (char *) malloc (80 + extras->dim*20); + + + sprintf (msg, "Unable to locate point %d [%f", + point, coords[0][point]); + for (i = 1; i < extras->dim; i++) + { + sprintf (msg, "%s %f", msg, coords[i][point]); + } + sprintf (msg, "%s]", msg); + CCTK_WARN (1, msg); + free (msg); + tmp = 1; /* mark as error */ + } + } + if (tmp) + { + if (myGH->whichproc) + { + free (myGH->whichproc); + myGH->whichproc = NULL; + } + return (-1); + } + + /* Now we want to resolve the num_points_from[]. Currently this is + the form of ( in 2 proc mode ) + P1: Num from P1 NFP2 + P2: NFP1 NFP2 + + and this needs to become + P1: Num to P1 NTP2 + P2: NTP1 NTP1 + + Since NTP1 = NFP2 (and this works in more proc mode too) + this is an all-to-all communication. + */ + CACTUS_MPI_ERROR (MPI_Alltoall (myGH->num_points_from, 1, PUGH_MPI_INT, + myGH->num_points_to, 1, PUGH_MPI_INT, + pughGH->PUGH_COMM_WORLD)); + +#ifdef PUGHINTERP_DEBUG + for (proc = 0; proc < nprocs; proc++) + { + printf ("processor %d <-> %d From: %d To: %d\n", + CCTK_MyProc (GH), proc, myGH->num_points_from[proc], + myGH->num_points_to[proc]); + } +#endif + + /* Great. Now we know how many to expect from each processor, + and how many to send to each processor. So first we have + to send the locations to the processors which hold our data. + This means I send coords[dim][point] to whichproc[point]. + I have num_points_from[proc] to send to each processor. + */ + + /* This is backwards in the broadcast location; the number of points + we are getting is how many everyone else is sending to us, + eg, num_points_to, not how many we get back from everyone else, + eg, num_points_from. The number we are sending, of course, is + all of our locations, eg, num_points */ + *num_local_points = 0; + for (proc = 0; proc < nprocs; proc++) + { + *num_local_points += myGH->num_points_to[proc]; + } + +#ifdef PUGHINTERP_DEBUG + printf ("processor %d gets %d points in total\n", + CCTK_MyProc (GH), *num_local_points); +#endif + + /* allocate the local coordinates array (sorted in processor order) + and the resulting coordinates array that I have to process */ + proc_coords = (CCTK_REAL *) malloc (extras->dim * num_points * + sizeof (CCTK_REAL)); + *local_coords = (CCTK_REAL *) malloc (extras->dim * *num_local_points * + sizeof (CCTK_REAL)); + + /* now sort my own coordinates as tupels of [extras->dim] */ + tmp = 0; + for (proc = 0; proc < nprocs; proc++) + { + for (point = 0; point < num_points; point++) + { + if (myGH->whichproc[point] == proc) + { + for (dim = 0; dim < extras->dim; dim++) + { + *proc_coords++ = coords[dim][point]; + } + myGH->indices[tmp++] = point; + } + } + } + proc_coords -= tmp * extras->dim; + + /* So load up the send and recv stuff */ + /* Send extras->dim elements per data point */ + myGH->send_count[0] = extras->dim * myGH->num_points_from[0]; + myGH->recv_count[0] = extras->dim * myGH->num_points_to[0]; + myGH->send_displ[0] = myGH->recv_displ[0] = 0; + for (proc = 1; proc < nprocs; proc++) + { + myGH->send_count[proc] = extras->dim * myGH->num_points_from[proc]; + myGH->recv_count[proc] = extras->dim * myGH->num_points_to[proc]; + myGH->send_displ[proc] = myGH->send_displ[proc-1] + + myGH->send_count[proc-1]; + myGH->recv_displ[proc] = myGH->recv_displ[proc-1] + + myGH->recv_count[proc-1]; + } + + /* Great, and now exchange the coordinates and collect the ones + that I have to process in *local_coords[] */ + CACTUS_MPI_ERROR (MPI_Alltoallv (proc_coords, myGH->send_count, + myGH->send_displ, PUGH_MPI_REAL, + *local_coords, myGH->recv_count, + myGH->recv_displ, PUGH_MPI_REAL, + pughGH->PUGH_COMM_WORLD)); + + /* don't need this anymore */ + free (proc_coords); + + return (0); +} +#endif /* CCTK_MPI */ + + + /*@@ + @routine PUGHInterp_CheckArguments + @date Thu 25 Jan 2001 + @author Thomas Radke + @desc + Checks the interpolation arguments passed in via + the flesh's general interpolation calling interface + + This routine also verifies that the parameters meet + the limitations of PUGHInterp's interpolation operators. + @enddesc + + @var GH + @vdesc Pointer to CCTK grid hierarchy + @vtype cGH * + @vio in + @endvar + @var num_dims + @vdesc dimensionality of the underlying grid + @vtype int + @vio in + @endvar + @var num_points + @vdesc number of points to interpolate at + @vtype int + @vio in + @endvar + @var num_in_arrays + @vdesc number of passed input arrays + @vtype int + @vio in + @endvar + @var num_out_arrays + @vdesc number of passed input arrays + @vtype int + @vio in + @endvar + @var interp_coord_array_types + @vdesc types of passed coordinates to interpolate at + @vtype int [num_dims] + @vio in + @endvar + + @returntype int + @returndesc + +1 for success + 0 for success but nothing to do + -1 for failure (wrong parameters passed or limitations not met) + @endreturndesc +@@*/ +static int PUGHInterp_CheckArguments (cGH *GH, + int num_dims, + int num_points, + int num_in_arrays, + int num_out_arrays, + int interp_coord_array_types[]) +{ + int i; + + + /* check for invalid arguments */ + if (num_dims < 0 || num_points < 0 || num_in_arrays < 0 || num_out_arrays < 0) + { + return (-1); + } + + /* check if there's anything to do at all */ + /* NOTE: num_points can be 0 in a collective call */ + if (num_dims == 0 || (CCTK_nProcs (GH) == 1 && num_points == 0) || + num_in_arrays == 0 || num_out_arrays == 0) + { + return (0); + } + + /* for now we can only deal with coordinates of type CCTK_REAL */ + for (i = 0; i < num_dims; i++) + { + if (interp_coord_array_types[i] != CCTK_VARIABLE_REAL) + { + CCTK_WARN (1, "Interpolation coordinates must be of type CCTK_REAL"); + return (-1); + } + } + + /* PUGHInterp's interpolation operators compute one output array + per input array */ + if (num_in_arrays != num_out_arrays) + { + CCTK_WARN (1, "Number of input arrays must match number of output arrays"); + return (-1); + } + + return (1); +} diff --git a/src/Startup.c b/src/Startup.c new file mode 100644 index 0000000..729e712 --- /dev/null +++ b/src/Startup.c @@ -0,0 +1,253 @@ + /*@@ + @file Startup.c + @date Sun Jul 04 1999 + @author Thomas Radke + @desc + Startup routines for PUGHInterp + @enddesc + @history + @endhistory + @version $Id$ + @@*/ + +#include + +#include "cctk.h" +#include "cctk_Interp.h" +#include "pughInterpGH.h" + +/* the rcs ID and its dummy function to use it */ +static char *rcsid = "$Header$"; +CCTK_FILEVERSION(CactusPUGH_PUGHInterp_Startup_c) + + +/* prototypes of routines defined in this source file */ +void PUGHInterp_Startup (void); +#ifdef CCTK_MPI +static void *PUGHInterp_SetupGH (tFleshConfig *config, + int convergence_level, + cGH *GH); +#endif + + +#ifdef CCTK_MPI + /*@@ + @routine PUGHInterp_SetupGH + @date Sun Jul 04 1999 + @author Thomas Radke + @desc + Allocates the GH extension structure and - for the MPI case - + the count/displacement buffers used in MPI_Alltoall() + and MPI_Alltoallv() + @enddesc + @calls CCTK_nProcs + + @returntype void * + @returndesc + pointer to the allocated GH extension structure + @endreturndesc +@@*/ +static void *PUGHInterp_SetupGH (tFleshConfig *config, + int convergence_level, + cGH *GH) +{ + int nprocs; + pughInterpGH *myGH; + + + /* suppress compiler warnings about unused variables */ + config = config; + convergence_level = convergence_level; + + myGH = (pughInterpGH *) malloc (sizeof (pughInterpGH)); + + /* allocate once for all fields of same type */ + nprocs = CCTK_nProcs (GH); + myGH->send_count = (int *) malloc (4 * nprocs * sizeof (int)); + myGH->send_displ = myGH->send_count + nprocs; + myGH->recv_count = myGH->send_displ + nprocs; + myGH->recv_displ = myGH->recv_count + nprocs; + + myGH->num_points_from = (CCTK_INT *) malloc (2 * nprocs * sizeof (CCTK_INT)); + myGH->num_points_to = myGH->num_points_from + nprocs; + + return (myGH); +} +#endif /* CCTK_MPI */ + + + /*@@ + @routine PUGHInterp_InterpGV_NthOrder + @date Wed 14 Feb 2001 + @author Thomas Radke + @desc + Wrappers for the different interpolation operators + registered for first/second/third order interpolation. + These wrappers just call the common interpolation routine + passing all arguments plus the interpolation order. + @enddesc + + @returntype int + @returndesc + the return code of the common interpolation routine + @endreturndesc +@@*/ +static int PUGHInterp_InterpGV_1stOrder (cGH *GH, + const char *coord_system, + int num_points, + int num_in_array_indices, + int num_out_arrays, + void *interp_coord_arrays[], + int interp_coord_array_types[], + int in_array_indices[], + void *out_arrays[], + int out_array_types[]) +{ + return (PUGHInterp_InterpGV (GH, 1, coord_system, num_points, + num_in_array_indices, num_out_arrays, + interp_coord_arrays, interp_coord_array_types, + in_array_indices, out_arrays, out_array_types)); +} + + +static int PUGHInterp_InterpGV_2ndOrder (cGH *GH, + const char *coord_system, + int num_points, + int num_in_array_indices, + int num_out_arrays, + void *interp_coord_arrays[], + int interp_coord_array_types[], + int in_array_indices[], + void *out_arrays[], + int out_array_types[]) +{ + return (PUGHInterp_InterpGV (GH, 2, coord_system, num_points, + num_in_array_indices, num_out_arrays, + interp_coord_arrays, interp_coord_array_types, + in_array_indices, out_arrays, out_array_types)); +} + + +static int PUGHInterp_InterpGV_3rdOrder (cGH *GH, + const char *coord_system, + int num_points, + int num_in_array_indices, + int num_out_arrays, + void *interp_coord_arrays[], + int interp_coord_array_types[], + int in_array_indices[], + void *out_arrays[], + int out_array_types[]) +{ + return (PUGHInterp_InterpGV (GH, 3, coord_system, num_points, + num_in_array_indices, num_out_arrays, + interp_coord_arrays, interp_coord_array_types, + in_array_indices, out_arrays, out_array_types)); +} + + +static int PUGHInterp_InterpLocal_1stOrder (cGH *GH, + int num_points, + int num_dims, + int num_in_arrays, + int num_out_arrays, + int coord_dims[], + void *coord_arrays[], + int coord_array_types[], + void *interp_coord_arrays[], + int interp_coord_array_types[], + void *in_arrays[], + int in_array_types[], + void *out_arrays[], + int out_array_types[]) +{ + return (PUGHInterp_InterpLocal (GH, 1, num_points, num_dims, + num_in_arrays, num_out_arrays, + coord_dims, coord_arrays, coord_array_types, + interp_coord_arrays, interp_coord_array_types, + in_arrays, in_array_types, + out_arrays, out_array_types)); +} + + +static int PUGHInterp_InterpLocal_2ndOrder (cGH *GH, + int num_points, + int num_dims, + int num_in_arrays, + int num_out_arrays, + int coord_dims[], + void *coord_arrays[], + int coord_array_types[], + void *interp_coord_arrays[], + int interp_coord_array_types[], + void *in_arrays[], + int in_array_types[], + void *out_arrays[], + int out_array_types[]) +{ + return (PUGHInterp_InterpLocal (GH, 2, num_points, num_dims, + num_in_arrays, num_out_arrays, + coord_dims, coord_arrays, coord_array_types, + interp_coord_arrays, interp_coord_array_types, + in_arrays, in_array_types, + out_arrays, out_array_types)); +} + + +static int PUGHInterp_InterpLocal_3rdOrder (cGH *GH, + int num_points, + int num_dims, + int num_in_arrays, + int num_out_arrays, + int coord_dims[], + void *coord_arrays[], + int coord_array_types[], + void *interp_coord_arrays[], + int interp_coord_array_types[], + void *in_arrays[], + int in_array_types[], + void *out_arrays[], + int out_array_types[]) +{ + return (PUGHInterp_InterpLocal (GH, 3, num_points, num_dims, + num_in_arrays, num_out_arrays, + coord_dims, coord_arrays, coord_array_types, + interp_coord_arrays, interp_coord_array_types, + in_arrays, in_array_types, + out_arrays, out_array_types)); +} + + + /*@@ + @routine PUGHInterp_Startup + @date Sun Jul 04 1999 + @author Thomas Radke + @desc + The startup registration routine for PUGHInterp. + Registers the PUGHInterp GH extensions setup routine + and the interpolation operators with the flesh. + @enddesc + @calls CCTK_RegisterGHExtensionSetupGH + CCTK_InterpRegisterOperatorGV + CCTK_InterpRegisterOperatorLocal +@@*/ +void PUGHInterp_Startup (void) +{ +#ifdef CCTK_MPI + CCTK_RegisterGHExtensionSetupGH (CCTK_RegisterGHExtension ("PUGHInterp"), + PUGHInterp_SetupGH); +#endif + + CCTK_InterpRegisterOperatorGV (PUGHInterp_InterpGV_1stOrder, + "first-order uniform cartesian"); + CCTK_InterpRegisterOperatorLocal (PUGHInterp_InterpLocal_1stOrder, + "first-order uniform cartesian"); + CCTK_InterpRegisterOperatorGV (PUGHInterp_InterpGV_2ndOrder, + "second-order uniform cartesian"); + CCTK_InterpRegisterOperatorLocal (PUGHInterp_InterpLocal_2ndOrder, + "second-order uniform cartesian"); + CCTK_InterpRegisterOperatorGV (PUGHInterp_InterpGV_3rdOrder, + "third-order uniform cartesian"); + CCTK_InterpRegisterOperatorLocal (PUGHInterp_InterpLocal_3rdOrder, + "third-order uniform cartesian"); +} diff --git a/src/make.code.defn b/src/make.code.defn new file mode 100644 index 0000000..b8d41c1 --- /dev/null +++ b/src/make.code.defn @@ -0,0 +1,5 @@ +# Main make.code.defn file for thorn PUGHInterp +# $Header$ + +# Source files in this directory +SRCS = Startup.c Operator.c Interpolate.c diff --git a/src/pughInterpGH.h b/src/pughInterpGH.h new file mode 100644 index 0000000..e6eef76 --- /dev/null +++ b/src/pughInterpGH.h @@ -0,0 +1,71 @@ + /*@@ + @file pughInterpGH.h + @date Sun Jul 04 1999 + @author Thomas Radke + @desc + Definition of GH extensions of thorn PUGHInterp and + declaration of function prototypes + @enddesc + @history + @endhistory +@@*/ + + +typedef struct +{ + int *send_count; /* number of elements to be sent */ + int *send_displ; /* offset of each element to be sent */ + int *recv_count; /* number of elements to be received */ + int *recv_displ; /* offset of each element to be received */ + + CCTK_INT *num_points_to; /* number of coordinate points to be sent */ + CCTK_INT *num_points_from; /* number of coordinate points to be received */ + + int *whichproc; /* which processor does point i belong to */ + int *indices; /* indices to sort from processor-sorted into + point-sorted order */ +} pughInterpGH; + + +/* prototypes of interpolation operators to be registered */ +int PUGHInterp_InterpGV (cGH *GH, + int order, + const char *coord_system, + int num_points, + int num_in_array_indices, + int num_out_arrays, + void *interp_coord_arrays[], + int interp_coord_array_types[], + int in_array_indices[], + void *out_arrays[], + int out_array_types[]); + +int PUGHInterp_InterpLocal (cGH *GH, + int order, + int num_points, + int num_dims, + int num_in_arrays, + int num_out_arrays, + int coord_dims[], + void *coord_arrays[], + int coord_array_types[], + void *interp_coord_arrays[], + int interp_coord_array_types[], + void *in_arrays[], + int in_array_types[], + void *out_arrays[], + int out_array_types[]); + +/* prototypes of routines used internally */ +int PUGHInterp_Interpolate (int order, + int num_points, + int num_dims, + int num_arrays, + int dims[ /* num_dims */ ], + CCTK_REAL coord[ /* num_dims * num_points */ ], + CCTK_REAL origin[ /* num_dims */ ], + CCTK_REAL delta[ /* num_dims */ ], + int in_types[ /* num_arrays */ ], + void *in_arrays[ /* num_arrays */ ], + int out_types[ /* num_arrays */ ], + void *out_arrays[ /* num_arrays */ ]); -- cgit v1.2.3