From 42a98ca518377949e8cfc902343a4b2d97a25564 Mon Sep 17 00:00:00 2001 From: schnetter Date: Mon, 15 Jul 2002 17:15:29 +0000 Subject: Added thorn TGRtensor. This is a utility thorn from the TGR code. It contains generic tensor operations, some useful constants, code for numeric derivatives, interfaces to Lapack and BLAS, and interfaces for Cactus routines. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinUtils/TGRtensor/trunk@2 b716e942-a2de-43ad-8f52-f3dfe468e4e7 --- COPYING | 340 ++++++++++++++++++++++++++++++++++++++ README | 8 + interface.ccl | 4 + param.ccl | 2 + schedule.ccl | 2 + src/cactus.F90 | 389 ++++++++++++++++++++++++++++++++++++++++++++ src/constants.F90 | 52 ++++++ src/conversion.F90 | 185 +++++++++++++++++++++ src/covderivs.F90 | 201 +++++++++++++++++++++++ src/covderivs2.F90 | 204 +++++++++++++++++++++++ src/depend.pl | 46 ++++++ src/derivs.F90 | 46 ++++++ src/derivs2.F90 | 40 +++++ src/hyperslab.F90 | 33 ++++ src/lapack.F90 | 51 ++++++ src/make.code.defn | 25 +++ src/make.code.deps | 11 ++ src/make.configuration.deps | 6 + src/matexp.F90 | 38 +++++ src/matinv.F90 | 54 ++++++ src/pointwise.F90 | 252 ++++++++++++++++++++++++++++ src/pointwise2.F90 | 148 +++++++++++++++++ src/rotation.F90 | 59 +++++++ src/tensor.F90 | 127 +++++++++++++++ src/tensor2.F90 | 107 ++++++++++++ 25 files changed, 2430 insertions(+) create mode 100644 COPYING create mode 100644 README create mode 100644 interface.ccl create mode 100644 param.ccl create mode 100644 schedule.ccl create mode 100644 src/cactus.F90 create mode 100644 src/constants.F90 create mode 100644 src/conversion.F90 create mode 100644 src/covderivs.F90 create mode 100644 src/covderivs2.F90 create mode 100644 src/depend.pl create mode 100644 src/derivs.F90 create mode 100644 src/derivs2.F90 create mode 100644 src/hyperslab.F90 create mode 100644 src/lapack.F90 create mode 100644 src/make.code.defn create mode 100644 src/make.code.deps create mode 100644 src/make.configuration.deps create mode 100644 src/matexp.F90 create mode 100644 src/matinv.F90 create mode 100644 src/pointwise.F90 create mode 100644 src/pointwise2.F90 create mode 100644 src/rotation.F90 create mode 100644 src/tensor.F90 create mode 100644 src/tensor2.F90 diff --git a/COPYING b/COPYING new file mode 100644 index 0000000..5b6e7c6 --- /dev/null +++ b/COPYING @@ -0,0 +1,340 @@ + GNU GENERAL PUBLIC LICENSE + Version 2, June 1991 + + Copyright (C) 1989, 1991 Free Software Foundation, Inc. + 59 Temple Place, Suite 330, Boston, MA 02111-1307 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 + + 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) + + 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 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) year 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..30b1a69 --- /dev/null +++ b/README @@ -0,0 +1,8 @@ +Cactus Code Thorn TGRtensor +Authors : Erik Schnetter +CVS info : $Header$ +-------------------------------------------------------------------------- + +Purpose of the thorn: + +It contains generic tensor operations and other generic stuff. diff --git a/interface.ccl b/interface.ccl new file mode 100644 index 0000000..5eb7d12 --- /dev/null +++ b/interface.ccl @@ -0,0 +1,4 @@ +# Interface definition for thorn TGRtensor +# $Header$ + +IMPLEMENTS: TGRtensor diff --git a/param.ccl b/param.ccl new file mode 100644 index 0000000..e5004ff --- /dev/null +++ b/param.ccl @@ -0,0 +1,2 @@ +# Parameter definitions for thorn TGRtensor +# $Header$ diff --git a/schedule.ccl b/schedule.ccl new file mode 100644 index 0000000..b6b05ea --- /dev/null +++ b/schedule.ccl @@ -0,0 +1,2 @@ +# Schedule definitions for thorn TGRtensor +# $Header$ diff --git a/src/cactus.F90 b/src/cactus.F90 new file mode 100644 index 0000000..5d0c278 --- /dev/null +++ b/src/cactus.F90 @@ -0,0 +1,389 @@ +! $Header$ + +#include "cctk.h" + +module cactus + implicit none + public + + interface + + ! from Cactus: + + subroutine cctk_barrier (ierr, cctkgh) + implicit none + integer ierr + CCTK_POINTER cctkgh + end subroutine cctk_barrier + + subroutine cctk_coordrange (ierr, cctkgh, lower, upper, dir, name, & + & systemname) + implicit none + integer ierr + CCTK_POINTER cctkgh + CCTK_REAL lower + CCTK_REAL upper + integer dir + character*(*) name + character*(*) systemname + end subroutine cctk_coordrange + + subroutine cctk_coordsystemhandle (handle, coordsystem) + implicit none + integer handle + character*(*) coordsystem + end subroutine cctk_coordsystemhandle + + function cctk_equals (arg1, arg2) + implicit none + integer cctk_equals + CCTK_POINTER arg1 + character*(*) arg2 + end function cctk_equals + + subroutine cctk_fortranstring (clength, cstring, fortstring) + implicit none + CCTK_INT clength ! intent(out) + CCTK_POINTER cstring + character*(*) fortstring + end subroutine cctk_fortranstring + + subroutine cctk_groupbboxgn (ierr, cctkgh, nbbox, bbox, groupname) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer nbbox + integer bbox(nbbox) + character*(*) groupname + end subroutine cctk_groupbboxgn + + subroutine cctk_groupgshgn (ierr, cctkgh, ngsh, gsh, groupname) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer ngsh + integer gsh(ngsh) + character*(*) groupname + end subroutine cctk_groupgshgn + + subroutine cctk_grouplbndgn (ierr, cctkgh, nlbnd, lbnd, groupname) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer nlbnd + integer lbnd(nlbnd) + character*(*) groupname + end subroutine cctk_grouplbndgn + + subroutine cctk_grouplshgn (ierr, cctkgh, nlsh, lsh, groupname) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer nlsh + integer lsh(nlsh) + character*(*) groupname + end subroutine cctk_grouplshgn + + subroutine cctk_groupnghostzonesgn & + (ierr, cctkgh, nnghostzones, nghostzones, groupname) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer nnghostzones + integer nghostzones(nnghostzones) + character*(*) groupname + end subroutine cctk_groupnghostzonesgn + + subroutine cctk_groupubndgn (ierr, cctkgh, nubnd, ubnd, groupname) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer nubnd + integer ubnd(nubnd) + character*(*) groupname + end subroutine cctk_groupubndgn + + subroutine cctk_info (thorn, message) + implicit none + character*(*) thorn + character*(*) message + end subroutine cctk_info + + subroutine cctk_interphandle (handle, interp) + implicit none + integer handle + character*(*) interp + end subroutine cctk_interphandle + + function cctk_isthornactive (name) + implicit none + integer cctk_isthornactive + character*(*) name + end function cctk_isthornactive + + function cctk_myproc (cctkgh) + implicit none + integer cctk_myproc + CCTK_POINTER cctkgh + end function cctk_myproc + + function cctk_nprocs (cctkgh) + implicit none + integer cctk_nprocs + CCTK_POINTER cctkgh + end function cctk_nprocs + + subroutine cctk_paramwarn (thorn, message) + implicit none + character*(*) thorn + character*(*) message + end subroutine cctk_paramwarn + + subroutine cctk_reductionarrayhandle (handle, reduction) + implicit none + integer handle + character*(*) reduction + end subroutine cctk_reductionarrayhandle + + subroutine cctk_registerbanner (ierr, banner) + implicit none + integer ierr + character*(*) banner + end subroutine cctk_registerbanner + + subroutine cctk_syncgroup (ierr, cctkgh, groupname) + implicit none + integer ierr + CCTK_POINTER cctkgh + character*(*) groupname + end subroutine cctk_syncgroup + + subroutine cctk_syncgroupi (ierr, cctkgh, group) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer group + end subroutine cctk_syncgroupi + + subroutine cctk_syncgroupwithvar (ierr, cctkgh, varname) + implicit none + integer ierr + CCTK_POINTER cctkgh + character*(*) varname + end subroutine cctk_syncgroupwithvar + + subroutine cctk_syncgroupwithvari (ierr, cctkgh, var) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer var + end subroutine cctk_syncgroupwithvari + + subroutine cctk_varindex (vi, vname) + implicit none + integer vi + character*(*) vname + end subroutine cctk_varindex + + subroutine cctk_vartypei (vartype, varindex) + implicit none + integer vartype + integer varindex + end subroutine cctk_vartypei + + subroutine cctk_warn (level, line, file, thorn, message) + implicit none + integer level + integer line + character*(*) file + character*(*) thorn + character*(*) message + end subroutine cctk_warn + + + + ! from thorn CactusBase/Boundary: + subroutine bndcopygn (ierr, cctkgh, sw, togroupname, fromgroupname) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer sw(3) + character*(*) togroupname + character*(*) fromgroupname + end subroutine bndcopygn + + subroutine bndflatgn (ierr, cctkgh, sw, groupname) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer sw(3) + character*(*) groupname + end subroutine bndflatgn + + subroutine bndradiativegn (ierr, cctkgh, sw, var0, v0, & + & togroupname, fromgroupname) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer sw(3) + CCTK_REAL var0 + CCTK_REAL v0 + character*(*) togroupname + character*(*) fromgroupname + end subroutine bndradiativegn + + subroutine bndradiativevn (ierr, cctkgh, sw, var0, v0, & + & tovarname, fromvarname) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer sw(3) + CCTK_REAL var0 + CCTK_REAL v0 + character*(*) tovarname + character*(*) fromvarname + end subroutine bndradiativevn + + subroutine bndrobinvn (ierr, cctkgh, sw, finf, npow, varname) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer sw(3) + CCTK_REAL finf + integer npow + character*(*) varname + end subroutine bndrobinvn + + subroutine bndscalargn (ierr, cctkgh, sw, var0, groupname) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer sw(3) + CCTK_REAL var0 + character*(*) groupname + end subroutine bndscalargn + + subroutine bndscalarvn (ierr, cctkgh, sw, var0, varname) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer sw(3) + CCTK_REAL var0 + character*(*) varname + end subroutine bndscalarvn + + + + ! from thorn CactusBase/CartGrid3D: + + subroutine cartsymgn (ierr, cctkgh, groupname) + implicit none + integer ierr + CCTK_POINTER cctkgh + character*(*) groupname + end subroutine cartsymgn + + subroutine setcartsymvn (ierr, cctkgh, sw, varname) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer sw(3) + character*(*) varname + end subroutine setcartsymvn + + + + ! from thorn AlphaThorns/Cart3D: + + subroutine cart3dsymvi (ierr, cctkgh, vi) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer vi + end subroutine cart3dsymvi + + subroutine cart3dsymvn (ierr, cctkgh, vn) + implicit none + integer ierr + CCTK_POINTER cctkgh + character*(*) vn + end subroutine cart3dsymvn + + subroutine cart3dsymgi (ierr, cctkgh, gi) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer gi + end subroutine cart3dsymgi + + subroutine cart3dsymgn (ierr, cctkgh, gn) + implicit none + integer ierr + CCTK_POINTER cctkgh + character*(*) gn + end subroutine cart3dsymgn + + subroutine cart3dsettensortypevi (ierr, cctkgh, nvars, vi, typename) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer nvars + integer vi(nvars) + character*(*) typename + end subroutine cart3dsettensortypevi + + subroutine cart3dsettensortypevn (ierr, cctkgh, vn, typename) + implicit none + integer ierr + CCTK_POINTER cctkgh + character*(*) vn + character*(*) typename + end subroutine cart3dsettensortypevn + + subroutine cart3dgetsymmetriesvi (ierr, cctkgh, sym, vi) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer sym(6) + integer vi + end subroutine cart3dgetsymmetriesvi + + subroutine cart3dgetsymmetriesvn (ierr, cctkgh, sym, vn) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer sym(6) + character*(*) vn + end subroutine cart3dgetsymmetriesvn + + subroutine cart3dgetsymmetryboundaries (ierr, cctkgh, symbnd) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer symbnd(6) + end subroutine cart3dgetsymmetryboundaries + + subroutine cart3dgetsymmetric (ierr, cctkgh, symmetric) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer symmetric(6) + end subroutine cart3dgetsymmetric + + subroutine cart3dgetstaggered (ierr, cctkgh, staggered) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer staggered(6) + end subroutine cart3dgetstaggered + + subroutine cart3dgetperiodic (ierr, cctkgh, periodic) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer periodic(6) + end subroutine cart3dgetperiodic + +end interface + +end module cactus diff --git a/src/constants.F90 b/src/constants.F90 new file mode 100644 index 0000000..d9b948d --- /dev/null +++ b/src/constants.F90 @@ -0,0 +1,52 @@ +! $Header$ + +#include "cctk.h" +#include "cctk_Parameters.h" + +module constants + implicit none + DECLARE_CCTK_PARAMETERS + private + + public delta2, delta3, delta4 + public eta4 + public epsilon2, epsilon3 + + public pi + + + + CCTK_REAL, parameter :: zero = 0 + + integer, parameter :: rk = kind(zero) + + + CCTK_REAL, parameter :: delta2(2,2) & + = reshape((/ 1,0, 0,1 /), (/2,2/)) + + CCTK_REAL, parameter :: delta3(3,3) & + = reshape((/ 1,0,0, 0,1,0, 0,0,1 /), (/3,3/)) + + CCTK_REAL, parameter :: delta4(0:3,0:3) & + = reshape((/ 1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1 /), (/4,4/)) + + + + CCTK_REAL, parameter :: eta4(0:3,0:3) & + = reshape((/ -1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1 /), (/4,4/)) + + + + CCTK_REAL, parameter :: epsilon2(2,2) & + = reshape ((/ 0,+1, -1,0 /), (/2,2/)) + + CCTK_REAL, parameter :: epsilon3(3,3,3) & + = reshape ((/ 0,0,0, 0,0,+1, 0,-1,0, & + & 0,0,-1, 0,0,0, +1,0,0, & + & 0,+1,0, -1,0,0, 0,0,0 /), (/3,3,3/)) + + + + CCTK_REAL, parameter :: pi = 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068_rk + +end module constants diff --git a/src/conversion.F90 b/src/conversion.F90 new file mode 100644 index 0000000..d6e30e9 --- /dev/null +++ b/src/conversion.F90 @@ -0,0 +1,185 @@ +! $Header$ + +#include "cctk.h" +#include "cctk_Parameters.h" + +module conversion + implicit none + DECLARE_CCTK_PARAMETERS + private + + public make_cylindrical2cartesian + public make_spherical2cartesian + public make_cartesian2spherical + +contains + + subroutine make_cylindrical2cartesian (xx, tt) + CCTK_REAL, intent(in) :: xx(3) + CCTK_REAL, intent(out) :: tt(3,3) + CCTK_REAL, parameter :: eps = 1.0d-20 + CCTK_REAL :: x,y,z + CCTK_REAL :: rho + + ! cylindrical coordinates: rho phi z + ! Cartesian coordinates: x y z + ! rho^2 = x^2 + y^2 + ! tan(phi) = y/x + + ! A[cart](p)_i = TT(p)_i^j A[cyl](p)_j + ! where p = p[cart] + + ! Transformation: + ! x = rho cos(phi) + ! y = rho sin(phi) + ! z = z + ! rho = sqrt(x^2+y^2) + ! phi = atan(y/x) + ! z = z + + ! Jacobian: + ! drho/dx = x/rho + ! drho/dy = y/rho + ! drho/dz = 0 + ! dphi/dx = -y/rho^2 + ! dphi/dy = x/rho^2 + ! dphi/dz = 0 + ! dz /dx = 0 + ! dz /dy = 0 + ! dz /dz = 1 + + x = xx(1) + y = xx(2) + z = xx(3) + + rho = sqrt(x**2+y**2) + + tt(1,1) = x/(rho+eps) + tt(2,1) = y/(rho+eps) + tt(3,1) = 0 + tt(1,2) = -y/(rho**2+eps) + tt(2,2) = x/(rho**2+eps) + tt(3,2) = 0 + tt(1,3) = 0 + tt(2,3) = 0 + tt(3,3) = 1 + + end subroutine make_cylindrical2cartesian + + subroutine make_spherical2cartesian (xx, tt) + CCTK_REAL, intent(in) :: xx(3) + CCTK_REAL, intent(out) :: tt(3,3) + CCTK_REAL, parameter :: eps = 1.0d-20 + CCTK_REAL :: x,y,z + CCTK_REAL :: rho,r + + ! Cartesian coordinates: x y z + ! spherical coordinates: r theta phi + ! r^2 = x^2 + y^2 + z^2 + ! cos(theta) = z/r + ! tan(phi) = y/x + + ! A[cart](p)_i = TT(p)_i^j A[spher](p)_j + ! where p = p[cart] + + ! Definition: + ! rho = r sin(theta) + ! z/rho = cos(theta) / sin(theta) + + ! Transformation: + ! x = r sin(theta) cos(phi) + ! y = r sin(theta) sin(phi) + ! z = r cos(theta) + ! r = sqrt(x^2+y^2+z^2) + ! theta = acos(z/r) + ! phi = atan(y/x) + + ! Jacobian: + ! dr /dx = x/r + ! dr /dy = y/r + ! dr /dz = z/r + ! dtheta/dx = xz/r^2rho ! 0 + ! dtheta/dy = yz/r^2rho ! 0 + ! dtheta/dz = -rho/r^2 ! 1/rho + ! dphi /dx = -y/rho^2 + ! dphi /dy = x/rho^2 + ! dphi /dz = 0 + + x = xx(1) + y = xx(2) + z = xx(3) + + rho = sqrt(x**2+y**2) + r = sqrt(x**2+y**2+z**2) + + tt(1,1) = x/(r+eps) + tt(2,1) = y/(r+eps) + tt(3,1) = z/(r+eps) + tt(1,2) = x*z/(r**2*rho+eps) + tt(2,2) = y*z/(r**2*rho+eps) + tt(3,2) = -rho**2/(r**2*rho+eps) + tt(1,3) = -y/(rho**2+eps) + tt(2,3) = x/(rho**2+eps) + tt(3,3) = 0 + + end subroutine make_spherical2cartesian + + subroutine make_cartesian2spherical (xx, tt) + CCTK_REAL, intent(in) :: xx(3) + CCTK_REAL, intent(out) :: tt(3,3) + CCTK_REAL, parameter :: eps = 1.0d-20 + CCTK_REAL :: x,y,z + CCTK_REAL :: rho,r + + ! Cartesian coordinates: x y z + ! spherical coordinates: r theta phi + ! r^2 = x^2 + y^2 + z^2 + ! cos(theta) = z/r + ! tan(phi) = y/x + + ! A[spher](p)_i = TT(p)_i^j A[cart](p)_j + ! where p = p[cart] + + ! Definition: + ! rho = r sin(theta) + ! z/rho = cos(theta) / sin(theta) + + ! Transformation: + ! x = r sin(theta) cos(phi) + ! y = r sin(theta) sin(phi) + ! z = r cos(theta) + ! r = sqrt(x^2+y^2+z^2) + ! theta = acos(z/r) + ! phi = atan(y/x) + + ! Jacobian: + ! dx/dr = sin(theta) cos(phi) = x/r + ! dx/dtheta = r cos(theta) cos(phi) = x z/rho + ! dx/dphi = - r sin(theta) sin(phi) = -y + ! dy/dr = sin(theta) sin(phi) = y/r + ! dy/dtheta = r cos(theta) sin(phi) = y z/rho + ! dy/dphi = r sin(theta) cos(phi) = x + ! dz/dr = cos(theta) = z/r + ! dz/dtheta = - r sin(theta) = -rho + ! dz/dphi = 0 = 0 + + x = xx(1) + y = xx(2) + z = xx(3) + + rho = sqrt(x**2+y**2) + r = sqrt(x**2+y**2+z**2) + + tt(1,1) = x/(r+eps) + tt(2,1) = x * z/(rho+eps) + tt(3,1) = -y + tt(1,2) = y/(r+eps) + tt(2,2) = y * z/(rho+eps) + tt(3,2) = x + tt(1,3) = z/(r+eps) + tt(2,3) = -rho + tt(3,3) = 0 + + end subroutine make_cartesian2spherical + +end module conversion diff --git a/src/covderivs.F90 b/src/covderivs.F90 new file mode 100644 index 0000000..50e3624 --- /dev/null +++ b/src/covderivs.F90 @@ -0,0 +1,201 @@ +! $Header$ + +#include "cctk.h" +#include "cctk_Parameters.h" + +module covderivs + implicit none + DECLARE_CCTK_PARAMETERS + private + + public calc_scalargrad + public calc_covectorgrad + public calc_vectorgrad + public calc_tensorgrad + public calc_scalargradgrad + public calc_covectorgradgrad + public calc_longitudinal + public calc_gradlongitudinal + +contains + + subroutine calc_scalargrad (f, df, gamma, gf) + CCTK_REAL, intent(in) :: f + CCTK_REAL, intent(in) :: df(3) + CCTK_REAL, intent(in) :: gamma(3,3,3) + CCTK_REAL, intent(out) :: gf(3) + integer :: i + ! f;i = f,i + do i=1,3 + gf(i) = df(i) + end do + end subroutine calc_scalargrad + + subroutine calc_covectorgrad (f, df, gamma, gf) + CCTK_REAL, intent(in) :: f(3) + CCTK_REAL, intent(in) :: df(3,3) + CCTK_REAL, intent(in) :: gamma(3,3,3) + CCTK_REAL, intent(out) :: gf(3,3) + integer :: i,j,k + ! f_i;j = f_i,j - Gamma^k_ij f_k + do i=1,3 + do j=1,3 + gf(i,j) = df(i,j) + do k=1,3 + gf(i,j) = gf(i,j) - gamma(k,i,j) * f(k) + end do + end do + end do + end subroutine calc_covectorgrad + + subroutine calc_vectorgrad (f, df, gamma, gf) + CCTK_REAL, intent(in) :: f(3) + CCTK_REAL, intent(in) :: df(3,3) + CCTK_REAL, intent(in) :: gamma(3,3,3) + CCTK_REAL, intent(out) :: gf(3,3) + integer :: i,j,k + ! f^i;j = f^i,j + Gamma^i_kj f^k + do i=1,3 + do j=1,3 + gf(i,j) = df(i,j) + do k=1,3 + gf(i,j) = gf(i,j) + gamma(i,k,j) * f(k) + end do + end do + end do + end subroutine calc_vectorgrad + + subroutine calc_tensorgrad (f, df, gamma, gf) + CCTK_REAL, intent(in) :: f(3,3) + CCTK_REAL, intent(in) :: df(3,3,3) + CCTK_REAL, intent(in) :: gamma(3,3,3) + CCTK_REAL, intent(out) :: gf(3,3,3) + integer :: i,j,k,l + ! f_ij;k = f_ij,k - Gamma^l_ik f_lj - Gamma^l_jk f_il + do i=1,3 + do j=1,3 + do k=1,3 + gf(i,j,k) = df(i,j,k) + do l=1,3 + gf(i,j,k) = gf(i,j,k) - gamma(l,i,k) * f(l,j) & + & - gamma(l,j,k) * f(i,l) + end do + end do + end do + end do + end subroutine calc_tensorgrad + + + + subroutine calc_scalargradgrad (f, df, ddf, gamma, ggf) + CCTK_REAL, intent(in) :: f + CCTK_REAL, intent(in) :: df(3) + CCTK_REAL, intent(in) :: ddf(3,3) + CCTK_REAL, intent(in) :: gamma(3,3,3) + CCTK_REAL, intent(out) :: ggf(3,3) + integer :: i,j,k + ! f;ij = f,ij - Gamma^k_ij f,k + do i=1,3 + do j=1,3 + ggf(i,j) = ddf(i,j) + do k=1,3 + ggf(i,j) = ggf(i,j) - gamma(k,i,j) * df(k) + end do + end do + end do + end subroutine calc_scalargradgrad + + subroutine calc_covectorgradgrad (f, df, ddf, gamma, dgamma, ggf) + CCTK_REAL, intent(in) :: f(3) + CCTK_REAL, intent(in) :: df(3,3) + CCTK_REAL, intent(in) :: ddf(3,3,3) + CCTK_REAL, intent(in) :: gamma(3,3,3), dgamma(3,3,3,3) + CCTK_REAL, intent(out) :: ggf(3,3,3) + CCTK_REAL :: gf(3,3) + integer :: i,j,k,l + ! f_i;j = f_i,j - Gamma^l_ij f_l + ! f_i;jk = f_i,jk - Gamma^l_ij,k f_l - Gamma^l_ij f_l,k + ! - Gamma^l_ik f_l;j - Gamma^l_jk f_i;l + call calc_covectorgrad (f, df, gamma, gf) + do i=1,3 + do j=1,3 + do k=1,3 + ggf(i,j,k) = ddf(i,j,k) + do l=1,3 + ggf(i,j,k) = ggf(i,j,k) & + - dgamma(l,i,j,k) * f(l) - gamma(l,i,j) * df(l,k) & + - gamma(l,i,k) * gf(l,j) - gamma(l,j,k) * gf(i,l) + end do + end do + end do + end do + end subroutine calc_covectorgradgrad + + subroutine calc_vectorgradgrad (f, df, ddf, gamma, dgamma, ggf) + CCTK_REAL, intent(in) :: f(3) + CCTK_REAL, intent(in) :: df(3,3) + CCTK_REAL, intent(in) :: ddf(3,3,3) + CCTK_REAL, intent(in) :: gamma(3,3,3), dgamma(3,3,3,3) + CCTK_REAL, intent(out) :: ggf(3,3,3) + CCTK_REAL :: gf(3,3) + integer :: i,j,k,l + ! f^i;j = f^i,j + Gamma^i_lj f^l + ! f^i;jk = f^i,jk + Gamma^i_lj,k f^l + Gamma^i_lj f^l,k + ! - Gamma^i_lk f^l;j - Gamma^l_jk f^i;l + call calc_vectorgrad (f, df, gamma, gf) + do i=1,3 + do j=1,3 + do k=1,3 + ggf(i,j,k) = ddf(i,j,k) + do l=1,3 + ggf(i,j,k) = ggf(i,j,k) & + + dgamma(i,l,j,k) * f(l) + gamma(i,l,j) * df(l,k) & + - gamma(i,l,k) * gf(l,j) - gamma(l,j,k) * gf(i,l) + end do + end do + end do + end do + end subroutine calc_vectorgradgrad + + + + subroutine calc_longitudinal (gf, gg, gu, lf) + CCTK_REAL, intent(in) :: gf(3,3) + CCTK_REAL, intent(in) :: gg(3,3), gu(3,3) + CCTK_REAL, intent(out) :: lf(3,3) + integer :: i,j,k,l + ! (Lf)_ij = D_i f_j + D_j f_i - 2/3 g_ij D^l f_l + do i=1,3 + do j=1,3 + lf(i,j) = gf(i,j) + gf(j,i) + do k=1,3 + do l=1,3 + lf(i,j) = lf(i,j) - (2.d0/3.d0) * gg(i,j) * gu(k,l) * gf(k,l) + end do + end do + end do + end do + end subroutine calc_longitudinal + + subroutine calc_gradlongitudinal (ggf, gg, gu, glf) + CCTK_REAL, intent(in) :: ggf(3,3,3) + CCTK_REAL, intent(in) :: gg(3,3), gu(3,3) + CCTK_REAL, intent(out) :: glf(3,3,3) + integer :: i,j,k,l,m + ! D_k (Lf)_ij = D_k (D_i f_j + D_j f_i - 2/3 g_ij D^l f_l) + do i=1,3 + do j=1,3 + do k=1,3 + glf(i,j,k) = ggf(j,i,k) + ggf(i,j,k) + do l=1,3 + do m=1,3 + glf(i,j,k) = glf(i,j,k) & + - (2.d0/3.d0) * gg(i,j) * gu(l,m) * ggf(l,m,k) + end do + end do + end do + end do + end do + end subroutine calc_gradlongitudinal + +end module covderivs diff --git a/src/covderivs2.F90 b/src/covderivs2.F90 new file mode 100644 index 0000000..15e0143 --- /dev/null +++ b/src/covderivs2.F90 @@ -0,0 +1,204 @@ +! $Header$ + +#include "cctk.h" +#include "cctk_Parameters.h" + +module covderivs2 + implicit none + DECLARE_CCTK_PARAMETERS + private + + public calc_2scalargrad + public calc_2covectorgrad + public calc_2vectorgrad + public calc_2tensorgrad + + public calc_2scalargradgrad + public calc_2covectorgradgrad + public calc_2vectorgradgrad + + public calc_2longitudinal + public calc_2gradlongitudinal + +contains + + subroutine calc_2scalargrad (f, df, gamma, gf) + CCTK_REAL, intent(in) :: f + CCTK_REAL, intent(in) :: df(2) + CCTK_REAL, intent(in) :: gamma(2,2,2) + CCTK_REAL, intent(out) :: gf(2) + integer :: i + ! f;i = f,i + do i=1,2 + gf(i) = df(i) + end do + end subroutine calc_2scalargrad + + subroutine calc_2covectorgrad (f, df, gamma, gf) + CCTK_REAL, intent(in) :: f(2) + CCTK_REAL, intent(in) :: df(2,2) + CCTK_REAL, intent(in) :: gamma(2,2,2) + CCTK_REAL, intent(out) :: gf(2,2) + integer :: i,j,k + ! f_i;j = f_i,j - Gamma^k_ij f_k + do i=1,2 + do j=1,2 + gf(i,j) = df(i,j) + do k=1,2 + gf(i,j) = gf(i,j) - gamma(k,i,j) * f(k) + end do + end do + end do + end subroutine calc_2covectorgrad + + subroutine calc_2vectorgrad (f, df, gamma, gf) + CCTK_REAL, intent(in) :: f(2) + CCTK_REAL, intent(in) :: df(2,2) + CCTK_REAL, intent(in) :: gamma(2,2,2) + CCTK_REAL, intent(out) :: gf(2,2) + integer :: i,j,k + ! f^i;j = f^i,j + Gamma^i_kj f^k + do i=1,2 + do j=1,2 + gf(i,j) = df(i,j) + do k=1,2 + gf(i,j) = gf(i,j) + gamma(i,k,j) * f(k) + end do + end do + end do + end subroutine calc_2vectorgrad + + subroutine calc_2tensorgrad (f, df, gamma, gf) + CCTK_REAL, intent(in) :: f(2,2) + CCTK_REAL, intent(in) :: df(2,2,2) + CCTK_REAL, intent(in) :: gamma(2,2,2) + CCTK_REAL, intent(out) :: gf(2,2,2) + integer :: i,j,k,l + ! f_ij;k = f_ij,k - Gamma^l_ik f_lj - Gamma^l_jk f_il + do i=1,2 + do j=1,2 + do k=1,2 + gf(i,j,k) = df(i,j,k) + do l=1,2 + gf(i,j,k) = gf(i,j,k) - gamma(l,i,k) * f(l,j) & + & - gamma(l,j,k) * f(i,l) + end do + end do + end do + end do + end subroutine calc_2tensorgrad + + + + subroutine calc_2scalargradgrad (f, df, ddf, gamma, ggf) + CCTK_REAL, intent(in) :: f + CCTK_REAL, intent(in) :: df(2) + CCTK_REAL, intent(in) :: ddf(2,2) + CCTK_REAL, intent(in) :: gamma(2,2,2) + CCTK_REAL, intent(out) :: ggf(2,2) + integer :: i,j,k + ! f;ij = f,ij - Gamma^k_ij f,k + do i=1,2 + do j=1,2 + ggf(i,j) = ddf(i,j) + do k=1,2 + ggf(i,j) = ggf(i,j) - gamma(k,i,j) * df(k) + end do + end do + end do + end subroutine calc_2scalargradgrad + + subroutine calc_2covectorgradgrad (f, df, ddf, gamma, dgamma, ggf) + CCTK_REAL, intent(in) :: f(2) + CCTK_REAL, intent(in) :: df(2,2) + CCTK_REAL, intent(in) :: ddf(2,2,2) + CCTK_REAL, intent(in) :: gamma(2,2,2), dgamma(2,2,2,2) + CCTK_REAL, intent(out) :: ggf(2,2,2) + CCTK_REAL :: gf(2,2) + integer :: i,j,k,l + ! f_i;j = f_i,j - Gamma^l_ij f_l + ! f_i;jk = f_i,jk - Gamma^l_ij,k f_l - Gamma^l_ij f_l,k + ! - Gamma^l_ik f_l;j - Gamma^l_jk f_i;l + call calc_2covectorgrad (f, df, gamma, gf) + do i=1,2 + do j=1,2 + do k=1,2 + ggf(i,j,k) = ddf(i,j,k) + do l=1,2 + ggf(i,j,k) = ggf(i,j,k) & + - dgamma(l,i,j,k) * f(l) - gamma(l,i,j) * df(l,k) & + - gamma(l,i,k) * gf(l,j) - gamma(l,j,k) * gf(i,l) + end do + end do + end do + end do + end subroutine calc_2covectorgradgrad + + subroutine calc_2vectorgradgrad (f, df, ddf, gamma, dgamma, ggf) + CCTK_REAL, intent(in) :: f(2) + CCTK_REAL, intent(in) :: df(2,2) + CCTK_REAL, intent(in) :: ddf(2,2,2) + CCTK_REAL, intent(in) :: gamma(2,2,2), dgamma(2,2,2,2) + CCTK_REAL, intent(out) :: ggf(2,2,2) + CCTK_REAL :: gf(2,2) + integer :: i,j,k,l + ! f^i;j = f^i,j + Gamma^i_lj f^l + ! f^i;jk = f^i,jk + Gamma^i_lj,k f^l + Gamma^i_lj f^l,k + ! - Gamma^i_lk f^l;j - Gamma^l_jk f^i;l + call calc_2vectorgrad (f, df, gamma, gf) + do i=1,2 + do j=1,2 + do k=1,2 + ggf(i,j,k) = ddf(i,j,k) + do l=1,2 + ggf(i,j,k) = ggf(i,j,k) & + + dgamma(i,l,j,k) * f(l) + gamma(i,l,j) * df(l,k) & + - gamma(i,l,k) * gf(l,j) - gamma(l,j,k) * gf(i,l) + end do + end do + end do + end do + end subroutine calc_2vectorgradgrad + + + + subroutine calc_2longitudinal (gf, gg, gu, lf) + CCTK_REAL, intent(in) :: gf(2,2) + CCTK_REAL, intent(in) :: gg(2,2), gu(2,2) + CCTK_REAL, intent(out) :: lf(2,2) + integer :: i,j,k,l + ! (Lf)_ij = D_i f_j + D_j f_i - 2/3 g_ij D^l f_l + do i=1,2 + do j=1,2 + lf(i,j) = gf(i,j) + gf(j,i) + do k=1,2 + do l=1,2 + lf(i,j) = lf(i,j) - (2.d0/3.d0) * gg(i,j) * gu(k,l) * gf(k,l) + end do + end do + end do + end do + end subroutine calc_2longitudinal + + subroutine calc_2gradlongitudinal (ggf, gg, gu, glf) + CCTK_REAL, intent(in) :: ggf(2,2,2) + CCTK_REAL, intent(in) :: gg(2,2), gu(2,2) + CCTK_REAL, intent(out) :: glf(2,2,2) + integer :: i,j,k,l,m + ! D_k (Lf)_ij = D_k (D_i f_j + D_j f_i - 2/3 g_ij D^l f_l) + do i=1,2 + do j=1,2 + do k=1,2 + glf(i,j,k) = ggf(j,i,k) + ggf(i,j,k) + do l=1,2 + do m=1,2 + glf(i,j,k) = glf(i,j,k) & + - (2.d0/3.d0) * gg(i,j) * gu(l,m) * ggf(l,m,k) + end do + end do + end do + end do + end do + end subroutine calc_2gradlongitudinal + +end module covderivs2 diff --git a/src/depend.pl b/src/depend.pl new file mode 100644 index 0000000..b41b8f6 --- /dev/null +++ b/src/depend.pl @@ -0,0 +1,46 @@ +#!/usr/bin/perl -w +# $Header$ + +# Create dependencies for Fortran 90 "use" and "include" statements + +$src = $ARGV[0]; +$srcfile = $ARGV[1]; +$dest = $ARGV[2]; +$srcdir = $ARGV[3]; +$moddir = $ARGV[4]; +$srcsuffix = $ARGV[5]; +$modsuffix = $ARGV[6]; +@otherdirs = @ARGV[7..$#ARGV]; + +print "# $src\n"; +print "$dest:"; + +open (FILE, $src); +while () { + if (/^\s*include\s*['"]([^'"]+)['"]/i) { + print " $srcdir$1"; + } elsif (/^\s*use\s+(\w+)/i) { + $found = 0; + if (! $found) { + if (-e "$srcdir$1$srcsuffix") { + $found = 1; + print " $moddir$1$modsuffix"; + } + } + if (! $found) { + foreach $dir (@otherdirs) { + if (-e "$dir$1$modsuffix") { + $found = 1; + print " $dir$1$modsuffix"; + last; + } + } + } + if (! $found) { + die "\nWhile tracing depencencies:\nFortran module $1 (referenced in file $srcfile) not found.\nAborting.\n\n"; + } + } +} +close (FILE); + +print "\n"; diff --git a/src/derivs.F90 b/src/derivs.F90 new file mode 100644 index 0000000..1d425c8 --- /dev/null +++ b/src/derivs.F90 @@ -0,0 +1,46 @@ +! $Header$ + +#ifndef TGR_INCLUDED + +#include "cctk.h" +#include "cctk_Parameters.h" + +module derivs + implicit none + private + public get_derivs + public get_derivs2 +contains +#endif + subroutine get_derivs (a, f, pos, off, dx) + CCTK_REAL, intent(in) :: a(*) + CCTK_REAL, intent(out) :: f(3) + integer, intent(in) :: pos, off(3) + CCTK_REAL, intent(in) :: dx(3) + integer :: i + do i=1,3 + f(i) = (a(pos+off(i)) - a(pos-off(i))) / (2*dx(i)) + end do + end subroutine get_derivs + subroutine get_derivs2 (a, f, pos, off, dx) + CCTK_REAL, intent(in) :: a(*) + CCTK_REAL, intent(out) :: f(3,3) + integer, intent(in) :: pos, off(3) + CCTK_REAL, intent(in) :: dx(3) + integer :: i + do i=1,3 + f(i,i) = (a(pos+off(i)) - 2*a(pos) + a(pos-off(i))) / dx(i)**2 + end do + f(1,2) = ( a(pos-off(1)-off(2)) - a(pos+off(1)-off(2)) & + & - a(pos-off(1)+off(2)) + a(pos+off(1)+off(2))) / (4*dx(1)*dx(2)) + f(2,1) = f(1,2) + f(1,3) = ( a(pos-off(1)-off(3)) - a(pos+off(1)-off(3)) & + & - a(pos-off(1)+off(3)) + a(pos+off(1)+off(3))) / (4*dx(1)*dx(3)) + f(3,1) = f(1,3) + f(2,3) = ( a(pos-off(2)-off(3)) - a(pos+off(2)-off(3)) & + & - a(pos-off(2)+off(3)) + a(pos+off(2)+off(3))) / (4*dx(2)*dx(3)) + f(3,2) = f(2,3) + end subroutine get_derivs2 +#ifndef TGR_INCLUDED +end module derivs +#endif diff --git a/src/derivs2.F90 b/src/derivs2.F90 new file mode 100644 index 0000000..4ca51bf --- /dev/null +++ b/src/derivs2.F90 @@ -0,0 +1,40 @@ +! $Header$ + +#ifndef TGR_INCLUDED + +#include "cctk.h" +#include "cctk_Parameters.h" + +module derivs2 + implicit none + private + public get_2derivs + public get_2derivs2 +contains +#endif + subroutine get_2derivs (a, f, pos, off, dx) + CCTK_REAL, intent(in) :: a(*) + CCTK_REAL, intent(out) :: f(2) + integer, intent(in) :: pos, off(2) + CCTK_REAL, intent(in) :: dx(2) + integer :: i + do i=1,2 + f(i) = (a(pos+off(i)) - a(pos-off(i))) / (2*dx(i)) + end do + end subroutine get_2derivs + subroutine get_2derivs2 (a, f, pos, off, dx) + CCTK_REAL, intent(in) :: a(*) + CCTK_REAL, intent(out) :: f(2,2) + integer, intent(in) :: pos, off(2) + CCTK_REAL, intent(in) :: dx(2) + integer :: i + do i=1,2 + f(i,i) = (a(pos+off(i)) - 2*a(pos) + a(pos-off(i))) / dx(i)**2 + end do + f(1,2) = ( a(pos-off(1)-off(2)) - a(pos+off(1)-off(2)) & + & - a(pos-off(1)+off(2)) + a(pos+off(1)+off(2))) / (4*dx(1)*dx(2)) + f(2,1) = f(1,2) + end subroutine get_2derivs2 +#ifndef TGR_INCLUDED +end module derivs2 +#endif diff --git a/src/hyperslab.F90 b/src/hyperslab.F90 new file mode 100644 index 0000000..e540cb2 --- /dev/null +++ b/src/hyperslab.F90 @@ -0,0 +1,33 @@ +! $Header$ + +#include "cctk.h" + +module hyperslab + implicit none + public + + interface + + subroutine hyperslab_fillhyperslab (ierr, cctkgh, target_proc, & + vindex, vtimelevel, & + hdim, global_startpoint, directions, lengths, downsample, & + nhdata, hdata, hsize) + implicit none + integer ierr + CCTK_POINTER cctkgh + integer target_proc + integer vindex + integer vtimelevel + integer hdim + integer global_startpoint(*) ! vdim + integer directions(*) ! vdim + integer lengths(hdim) + integer downsample(hdim) + integer nhdata + CCTK_REAL hdata(*) + integer hsize(hdim) + end subroutine hyperslab_fillhyperslab + + end interface + +end module hyperslab diff --git a/src/lapack.F90 b/src/lapack.F90 new file mode 100644 index 0000000..5bcd70b --- /dev/null +++ b/src/lapack.F90 @@ -0,0 +1,51 @@ +! $Header$ + +#include "cctk.h" + +module lapack + implicit none + public + + interface geev + subroutine sgeev (jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, & + & ldvr, work, lwork, info) + character jobvl, jobvr + integer info, lda, ldvl, ldvr, lwork, n + real a(lda,n), vl(ldvl,n), vr(ldvr,n), & + & wi(n), work(lwork), wr(n) + end subroutine sgeev + subroutine dgeev (jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, & + & ldvr, work, lwork, info) + character jobvl, jobvr + integer info, lda, ldvl, ldvr, lwork, n + double precision a(lda,n), vl(ldvl,n), vr(ldvr,n), & + & wi(n), work(lwork), wr(n) + end subroutine dgeev + end interface + + interface gesv + subroutine sgesv (n, nrhs, a, lda, ipiv, b, ldb, info) + implicit none + integer n + integer nrhs + integer lda + real a(lda,n) + integer ipiv(n) + integer ldb + real b(ldb,nrhs) + integer info + end subroutine sgesv + subroutine dgesv (n, nrhs, a, lda, ipiv, b, ldb, info) + implicit none + integer n + integer nrhs + integer lda + double precision a(lda,n) + integer ipiv(n) + integer ldb + double precision b(ldb,nrhs) + integer info + end subroutine dgesv + end interface + +end module lapack diff --git a/src/make.code.defn b/src/make.code.defn new file mode 100644 index 0000000..4b157e4 --- /dev/null +++ b/src/make.code.defn @@ -0,0 +1,25 @@ +# Main make.code.defn file for thorn TGRtensor -*-Makefile-*- +# $Header$ + +# Source files in this directory +SRCS = cactus.F90 \ + constants.F90 \ + conversion.F90 \ + covderivs.F90 \ + covderivs2.F90 \ + derivs.F90 \ + derivs2.F90 \ + hyperslab.F90 \ + hyperslab.c \ + lapack.F90 \ + matexp.F90 \ + matinv.F90 \ + rotation.F90 \ + tensor.F90 \ + tensor2.F90 \ + pointwise.F90 \ + pointwise2.F90 + +# Subdirectories containing source files +SUBDIRS = + diff --git a/src/make.code.deps b/src/make.code.deps new file mode 100644 index 0000000..e8ea1b9 --- /dev/null +++ b/src/make.code.deps @@ -0,0 +1,11 @@ +# Main make.code.deps file for thorn TGRtensor -*-Makefile-*- +# $Header$ + +USESTHORNS = + +# Automatically create dependencies for Fortran modules and Fortran includes +define F90_DEPENDENCIES +$(F_DEPEND) $(INC_DIRS:%=-I%) $(EXTRA_DEFINES:%=-D%) -DFCODE $< $(F_DEPEND_OUT) +$(DEPENDENCY_FIXER) +dir=`pwd`; $(CPP) $(INC_DIRS:%=-I%) $(EXTRA_DEFINES:%=-D%) -DFCODE $< | $(PERL) $(SRCDIR)/depend.pl - $< $(basename $(notdir $<)).F90.o $(SRCDIR)/ ./ .F90 .F90.o $(USESTHORNS:%=$$dir/../%/) >> $@ || { rm $@; exit 1; } +endef diff --git a/src/make.configuration.deps b/src/make.configuration.deps new file mode 100644 index 0000000..6a85f9c --- /dev/null +++ b/src/make.configuration.deps @@ -0,0 +1,6 @@ +# Main make.configuration.deps file for thorn TGRtensor -*-Makefile-*- +# $Header$ + +USESTHORNS = + +$(CCTK_LIBDIR)$(DIRSEP)libTGRtensor.a: $(USESTHORNS:%=$(CCTK_LIBDIR)$(DIRSEP)lib%.a) diff --git a/src/matexp.F90 b/src/matexp.F90 new file mode 100644 index 0000000..d2ecc77 --- /dev/null +++ b/src/matexp.F90 @@ -0,0 +1,38 @@ +! $Header$ + +#include "cctk.h" +#include "cctk_Parameters.h" + +module matexp + use constants + implicit none + DECLARE_CCTK_PARAMETERS + private + + public calc_exp3 + +contains + + subroutine calc_exp3 (h3, g3) + CCTK_REAL, intent(in) :: h3(3,3) + CCTK_REAL, intent(out) :: g3(3,3) + CCTK_REAL :: nfact + CCTK_REAL :: tmp(3,3) + integer :: n + integer :: i, j + + g3 = delta3 + + nfact = 1 + tmp = delta3 + do n = 1, 18 + nfact = nfact * n + tmp = matmul (h3, tmp) + + ! exp(x) = sum_n x^n / n! + g3 = g3 + tmp / nfact + end do + + end subroutine calc_exp3 + +end module matexp diff --git a/src/matinv.F90 b/src/matinv.F90 new file mode 100644 index 0000000..ab6c515 --- /dev/null +++ b/src/matinv.F90 @@ -0,0 +1,54 @@ +! $Header$ + +#include "cctk.h" +#include "cctk_Parameters.h" + +module matinv + use constants + use lapack + implicit none + DECLARE_CCTK_PARAMETERS + private + + public calc_inv3 + public calc_inv4 + +contains + + subroutine calc_inv3 (g3, gu3) + CCTK_REAL, intent(in) :: g3(3,3) + CCTK_REAL, intent(out) :: gu3(3,3) + CCTK_REAL :: tmp(3,3) + integer :: ipiv(3) + integer :: info + character :: msg*1000 + + tmp = g3 + gu3 = delta3 + call gesv (3, 3, tmp, 3, ipiv, gu3, 3, info) + + if (info /= 0) then + write (msg, '("Error in call to GESV, info=",i2)') info + call CCTK_WARN (0, trim(msg)) + end if + end subroutine calc_inv3 + + subroutine calc_inv4 (g4, gu4) + CCTK_REAL, intent(in) :: g4(0:3,0:3) + CCTK_REAL, intent(out) :: gu4(0:3,0:3) + CCTK_REAL :: tmp(0:3,0:3) + integer :: ipiv(0:3) + integer :: info + character :: msg*1000 + + tmp = g4 + gu4 = delta4 + call gesv (4, 4, tmp, 4, ipiv, gu4, 4, info) + + if (info /= 0) then + write (msg, '("Error in call to GESV, info=",i2)') info + call CCTK_WARN (0, trim(msg)) + end if + end subroutine calc_inv4 + +end module matinv diff --git a/src/pointwise.F90 b/src/pointwise.F90 new file mode 100644 index 0000000..827ad6e --- /dev/null +++ b/src/pointwise.F90 @@ -0,0 +1,252 @@ +! $Header$ + +#include "cctk.h" +#include "cctk_Parameters.h" + +module pointwise +!!$ use derivs + implicit none + DECLARE_CCTK_PARAMETERS + private + public calc_position + public calc_offsets + public get_scalar + public get_vector + public get_tensor + public get_connections + public set_scalar + public set_vector + public set_tensor + public set_connections + public get_scalarderivs + public get_vectorderivs + public get_tensorderivs + public get_scalarderivs2 + public get_vectorderivs2 + public get_tensorderivs2 +contains +#define TGR_INCLUDED +#include "derivs.F90" +#undef TGR_INCLUDED + subroutine calc_position (shape, i,j,k, pos) + integer, intent(in) :: shape(3) + integer, intent(in) :: i,j,k + integer, intent(out) :: pos + pos = 1 + i-1 + shape(1) * (j-1 + shape(2) * (k-1)) + end subroutine calc_position + subroutine calc_offsets (shape, off) + integer, intent(in) :: shape(3) + integer, intent(out) :: off(3) + off(1) = 1 + off(2) = shape(1) + off(3) = shape(1) * shape(2) + end subroutine calc_offsets + subroutine get_scalar (a, f, pos) + CCTK_REAL, intent(in) :: a(*) + CCTK_REAL, intent(out) :: f + integer, intent(in) :: pos + f = a(pos) + end subroutine get_scalar + subroutine get_vector (ax,ay,az, f, pos) + CCTK_REAL, intent(in) :: ax(*),ay(*),az(*) + CCTK_REAL, intent(out) :: f(3) + integer, intent(in) :: pos + f(1) = ax(pos) + f(2) = ay(pos) + f(3) = az(pos) + end subroutine get_vector + subroutine get_tensor (axx,axy,axz,ayy,ayz,azz, f, pos) + CCTK_REAL, intent(in) :: axx(*),axy(*),axz(*),ayy(*),ayz(*),azz(*) + CCTK_REAL, intent(out) :: f(3,3) + integer, intent(in) :: pos + f(1,1) = axx(pos) + f(1,2) = axy(pos) + f(1,3) = axz(pos) + f(2,1) = axy(pos) + f(2,2) = ayy(pos) + f(2,3) = ayz(pos) + f(3,1) = axz(pos) + f(3,2) = ayz(pos) + f(3,3) = azz(pos) + end subroutine get_tensor + subroutine get_connections (& + gammaxxx,gammaxxy,gammaxxz,gammaxyy,gammaxyz,gammaxzz, & + gammayxx,gammayxy,gammayxz,gammayyy,gammayyz,gammayzz, & + gammazxx,gammazxy,gammazxz,gammazyy,gammazyz,gammazzz, gamma, pos) + CCTK_REAL, intent(in) :: gammaxxx(*),gammaxxy(*),gammaxxz(*) + CCTK_REAL, intent(in) :: gammaxyy(*),gammaxyz(*),gammaxzz(*) + CCTK_REAL, intent(in) :: gammayxx(*),gammayxy(*),gammayxz(*) + CCTK_REAL, intent(in) :: gammayyy(*),gammayyz(*),gammayzz(*) + CCTK_REAL, intent(in) :: gammazxx(*),gammazxy(*),gammazxz(*) + CCTK_REAL, intent(in) :: gammazyy(*),gammazyz(*),gammazzz(*) + CCTK_REAL, intent(out) :: gamma(3,3,3) + integer, intent(in) :: pos + gamma(1,1,1) = gammaxxx(pos) + gamma(1,1,2) = gammaxxy(pos) + gamma(1,1,3) = gammaxxz(pos) + gamma(1,2,1) = gammaxxy(pos) + gamma(1,2,2) = gammaxyy(pos) + gamma(1,2,3) = gammaxyz(pos) + gamma(1,3,1) = gammaxxz(pos) + gamma(1,3,2) = gammaxyz(pos) + gamma(1,3,3) = gammaxzz(pos) + gamma(2,1,1) = gammayxx(pos) + gamma(2,1,2) = gammayxy(pos) + gamma(2,1,3) = gammayxz(pos) + gamma(2,2,1) = gammayxy(pos) + gamma(2,2,2) = gammayyy(pos) + gamma(2,2,3) = gammayyz(pos) + gamma(2,3,1) = gammayxz(pos) + gamma(2,3,2) = gammayyz(pos) + gamma(2,3,3) = gammayzz(pos) + gamma(3,1,1) = gammazxx(pos) + gamma(3,1,2) = gammazxy(pos) + gamma(3,1,3) = gammazxz(pos) + gamma(3,2,1) = gammazxy(pos) + gamma(3,2,2) = gammazyy(pos) + gamma(3,2,3) = gammazyz(pos) + gamma(3,3,1) = gammazxz(pos) + gamma(3,3,2) = gammazyz(pos) + gamma(3,3,3) = gammazzz(pos) + end subroutine get_connections + subroutine set_scalar (f, a, pos) + CCTK_REAL, intent(in) :: f + CCTK_REAL :: a(*) + integer, intent(in) :: pos + a(pos) = f + end subroutine set_scalar + subroutine set_vector (f, ax,ay,az, pos) + CCTK_REAL, intent(in) :: f(3) + CCTK_REAL :: ax(*),ay(*),az(*) + integer, intent(in) :: pos + ax(pos) = f(1) + ay(pos) = f(2) + az(pos) = f(3) + end subroutine set_vector + subroutine set_tensor (f, axx,axy,axz,ayy,ayz,azz, pos) + CCTK_REAL, intent(in) :: f(3,3) + CCTK_REAL :: axx(*),axy(*),axz(*),ayy(*),ayz(*),azz(*) + integer, intent(in) :: pos + axx(pos) = f(1,1) + axy(pos) = f(1,2) + axz(pos) = f(1,3) + ayy(pos) = f(2,2) + ayz(pos) = f(2,3) + azz(pos) = f(3,3) + end subroutine set_tensor + subroutine set_connections (gamma, & + gammaxxx,gammaxxy,gammaxxz,gammaxyy,gammaxyz,gammaxzz, & + gammayxx,gammayxy,gammayxz,gammayyy,gammayyz,gammayzz, & + gammazxx,gammazxy,gammazxz,gammazyy,gammazyz,gammazzz, pos) + CCTK_REAL, intent(in) :: gamma(3,3,3) + CCTK_REAL :: gammaxxx(*),gammaxxy(*),gammaxxz(*) + CCTK_REAL :: gammaxyy(*),gammaxyz(*),gammaxzz(*) + CCTK_REAL :: gammayxx(*),gammayxy(*),gammayxz(*) + CCTK_REAL :: gammayyy(*),gammayyz(*),gammayzz(*) + CCTK_REAL :: gammazxx(*),gammazxy(*),gammazxz(*) + CCTK_REAL :: gammazyy(*),gammazyz(*),gammazzz(*) + integer, intent(in) :: pos + gammaxxx(pos) = gamma(1,1,1) + gammaxxy(pos) = gamma(1,1,2) + gammaxxz(pos) = gamma(1,1,3) + gammaxyy(pos) = gamma(1,2,2) + gammaxyz(pos) = gamma(1,2,3) + gammaxzz(pos) = gamma(1,3,3) + gammayxx(pos) = gamma(2,1,1) + gammayxy(pos) = gamma(2,1,2) + gammayxz(pos) = gamma(2,1,3) + gammayyy(pos) = gamma(2,2,2) + gammayyz(pos) = gamma(2,2,3) + gammayzz(pos) = gamma(2,3,3) + gammazxx(pos) = gamma(3,1,1) + gammazxy(pos) = gamma(3,1,2) + gammazxz(pos) = gamma(3,1,3) + gammazyy(pos) = gamma(3,2,2) + gammazyz(pos) = gamma(3,2,3) + gammazzz(pos) = gamma(3,3,3) + end subroutine set_connections + subroutine get_scalarderivs (a, f, pos, off, dx) + CCTK_REAL, intent(in) :: a(*) + CCTK_REAL, intent(out) :: f(3) + integer, intent(in) :: pos, off(3) + CCTK_REAL, intent(in) :: dx(3) + call get_derivs (a, f, pos, off, dx) + end subroutine get_scalarderivs + subroutine get_vectorderivs (ax,ay,az, f, pos, off, dx) + CCTK_REAL, intent(in) :: ax(*),ay(*),az(*) + CCTK_REAL, intent(out) :: f(3,3) + integer, intent(in) :: pos, off(3) + CCTK_REAL, intent(in) :: dx(3) + CCTK_REAL :: fx(3),fy(3),fz(3) + call get_derivs (ax, fx, pos, off, dx) + call get_derivs (ay, fy, pos, off, dx) + call get_derivs (az, fz, pos, off, dx) + f(1,:) = fx + f(2,:) = fy + f(3,:) = fz + end subroutine get_vectorderivs + subroutine get_tensorderivs (axx,axy,axz,ayy,ayz,azz, f, pos, off, dx) + CCTK_REAL, intent(in) :: axx(*),axy(*),axz(*),ayy(*),ayz(*),azz(*) + CCTK_REAL, intent(out) :: f(3,3,3) + integer, intent(in) :: pos, off(3) + CCTK_REAL, intent(in) :: dx(3) + CCTK_REAL :: fxx(3),fxy(3),fxz(3),fyy(3),fyz(3),fzz(3) + call get_derivs (axx, fxx, pos, off, dx) + call get_derivs (axy, fxy, pos, off, dx) + call get_derivs (axz, fxz, pos, off, dx) + call get_derivs (ayy, fyy, pos, off, dx) + call get_derivs (ayz, fyz, pos, off, dx) + call get_derivs (azz, fzz, pos, off, dx) + f(1,1,:) = fxx + f(1,2,:) = fxy + f(1,3,:) = fxz + f(2,1,:) = fxy + f(2,2,:) = fyy + f(2,3,:) = fyz + f(3,1,:) = fxz + f(3,2,:) = fyz + f(3,3,:) = fzz + end subroutine get_tensorderivs + subroutine get_scalarderivs2 (a, f, pos, off, dx) + CCTK_REAL, intent(in) :: a(*) + CCTK_REAL, intent(out) :: f(3,3) + integer, intent(in) :: pos, off(3) + CCTK_REAL, intent(in) :: dx(3) + call get_derivs2 (a, f, pos, off, dx) + end subroutine get_scalarderivs2 + subroutine get_vectorderivs2 (ax,ay,az, f, pos, off, dx) + CCTK_REAL, intent(in) :: ax(*),ay(*),az(*) + CCTK_REAL, intent(out) :: f(3,3,3) + integer, intent(in) :: pos, off(3) + CCTK_REAL, intent(in) :: dx(3) + CCTK_REAL :: fx(3,3),fy(3,3),fz(3,3) + call get_derivs2 (ax, fx, pos, off, dx) + call get_derivs2 (ay, fy, pos, off, dx) + call get_derivs2 (az, fz, pos, off, dx) + f(1,:,:) = fx + f(2,:,:) = fy + f(3,:,:) = fz + end subroutine get_vectorderivs2 + subroutine get_tensorderivs2 (axx,axy,axz,ayy,ayz,azz, f, pos, off, dx) + CCTK_REAL, intent(in) :: axx(*),axy(*),axz(*),ayy(*),ayz(*),azz(*) + CCTK_REAL, intent(out) :: f(3,3,3,3) + integer, intent(in) :: pos, off(3) + CCTK_REAL, intent(in) :: dx(3) + CCTK_REAL :: fxx(3,3),fxy(3,3),fxz(3,3),fyy(3,3),fyz(3,3),fzz(3,3) + call get_derivs2 (axx, fxx, pos, off, dx) + call get_derivs2 (axy, fxy, pos, off, dx) + call get_derivs2 (axz, fxz, pos, off, dx) + call get_derivs2 (ayy, fyy, pos, off, dx) + call get_derivs2 (ayz, fyz, pos, off, dx) + call get_derivs2 (azz, fzz, pos, off, dx) + f(1,1,:,:) = fxx + f(1,2,:,:) = fxy + f(1,3,:,:) = fxz + f(2,1,:,:) = fxy + f(2,2,:,:) = fyy + f(2,3,:,:) = fyz + f(3,1,:,:) = fxz + f(3,2,:,:) = fyz + f(3,3,:,:) = fzz + end subroutine get_tensorderivs2 +end module pointwise diff --git a/src/pointwise2.F90 b/src/pointwise2.F90 new file mode 100644 index 0000000..1759b28 --- /dev/null +++ b/src/pointwise2.F90 @@ -0,0 +1,148 @@ +! $Header$ + +#include "cctk.h" +#include "cctk_Parameters.h" + +module pointwise2 +!!$ use derivs2 + implicit none + DECLARE_CCTK_PARAMETERS + private + public calc_2position + public calc_2offsets + public get_2scalar + public get_2vector + public get_2tensor + public set_2scalar + public set_2vector + public set_2tensor + public get_2scalarderivs + public get_2vectorderivs + public get_2tensorderivs + public get_2scalarderivs2 + public get_2vectorderivs2 + public get_2tensorderivs2 +contains +#define TGR_INCLUDED +#include "derivs2.F90" +#undef TGR_INCLUDED + subroutine calc_2position (shape, i,j, pos) + integer, intent(in) :: shape(2) + integer, intent(in) :: i,j + integer, intent(out) :: pos + pos = 1 + i-1 + shape(1) * (j-1) + end subroutine calc_2position + subroutine calc_2offsets (shape, off) + integer, intent(in) :: shape(2) + integer, intent(out) :: off(2) + off(1) = 1 + off(2) = shape(1) + end subroutine calc_2offsets + subroutine get_2scalar (a, f, pos) + CCTK_REAL, intent(in) :: a(*) + CCTK_REAL, intent(out) :: f + integer, intent(in) :: pos + f = a(pos) + end subroutine get_2scalar + subroutine get_2vector (ax,ay, f, pos) + CCTK_REAL, intent(in) :: ax(*),ay(*) + CCTK_REAL, intent(out) :: f(2) + integer, intent(in) :: pos + f(1) = ax(pos) + f(2) = ay(pos) + end subroutine get_2vector + subroutine get_2tensor (axx,axy,ayy, f, pos) + CCTK_REAL, intent(in) :: axx(*),axy(*),ayy(*) + CCTK_REAL, intent(out) :: f(2,2) + integer, intent(in) :: pos + f(1,1) = axx(pos) + f(1,2) = axy(pos) + f(2,1) = axy(pos) + f(2,2) = ayy(pos) + end subroutine get_2tensor + subroutine set_2scalar (f, a, pos) + CCTK_REAL, intent(in) :: f + CCTK_REAL :: a(*) + integer, intent(in) :: pos + a(pos) = f + end subroutine set_2scalar + subroutine set_2vector (f, ax,ay, pos) + CCTK_REAL, intent(in) :: f(2) + CCTK_REAL :: ax(*),ay(*) + integer, intent(in) :: pos + ax(pos) = f(1) + ay(pos) = f(2) + end subroutine set_2vector + subroutine set_2tensor (f, axx,axy,ayy, pos) + CCTK_REAL, intent(in) :: f(2,2) + CCTK_REAL :: axx(*),axy(*),ayy(*) + integer, intent(in) :: pos + axx(pos) = f(1,1) + axy(pos) = f(1,2) + ayy(pos) = f(2,2) + end subroutine set_2tensor + subroutine get_2scalarderivs (a, f, pos, off, dx) + CCTK_REAL, intent(in) :: a(*) + CCTK_REAL, intent(out) :: f(2) + integer, intent(in) :: pos, off(2) + CCTK_REAL, intent(in) :: dx(2) + call get_2derivs (a, f, pos, off, dx) + end subroutine get_2scalarderivs + subroutine get_2vectorderivs (ax,ay, f, pos, off, dx) + CCTK_REAL, intent(in) :: ax(*),ay(*) + CCTK_REAL, intent(out) :: f(2,2) + integer, intent(in) :: pos, off(2) + CCTK_REAL, intent(in) :: dx(2) + CCTK_REAL :: fx(2),fy(2) + call get_2derivs (ax, fx, pos, off, dx) + call get_2derivs (ay, fy, pos, off, dx) + f(1,:) = fx + f(2,:) = fy + end subroutine get_2vectorderivs + subroutine get_2tensorderivs (axx,axy,ayy, f, pos, off, dx) + CCTK_REAL, intent(in) :: axx(*),axy(*),ayy(*) + CCTK_REAL, intent(out) :: f(2,2,2) + integer, intent(in) :: pos, off(2) + CCTK_REAL, intent(in) :: dx(2) + CCTK_REAL :: fxx(2),fxy(2),fyy(2) + call get_2derivs (axx, fxx, pos, off, dx) + call get_2derivs (axy, fxy, pos, off, dx) + call get_2derivs (ayy, fyy, pos, off, dx) + f(1,1,:) = fxx + f(1,2,:) = fxy + f(2,1,:) = fxy + f(2,2,:) = fyy + end subroutine get_2tensorderivs + subroutine get_2scalarderivs2 (a, f, pos, off, dx) + CCTK_REAL, intent(in) :: a(*) + CCTK_REAL, intent(out) :: f(2,2) + integer, intent(in) :: pos, off(2) + CCTK_REAL, intent(in) :: dx(2) + call get_2derivs2 (a, f, pos, off, dx) + end subroutine get_2scalarderivs2 + subroutine get_2vectorderivs2 (ax,ay, f, pos, off, dx) + CCTK_REAL, intent(in) :: ax(*),ay(*) + CCTK_REAL, intent(out) :: f(2,2,2) + integer, intent(in) :: pos, off(2) + CCTK_REAL, intent(in) :: dx(2) + CCTK_REAL :: fx(2,2),fy(2,2) + call get_2derivs2 (ax, fx, pos, off, dx) + call get_2derivs2 (ay, fy, pos, off, dx) + f(1,:,:) = fx + f(2,:,:) = fy + end subroutine get_2vectorderivs2 + subroutine get_2tensorderivs2 (axx,axy,ayy, f, pos, off, dx) + CCTK_REAL, intent(in) :: axx(*),axy(*),ayy(*) + CCTK_REAL, intent(out) :: f(2,2,2,2) + integer, intent(in) :: pos, off(2) + CCTK_REAL, intent(in) :: dx(2) + CCTK_REAL :: fxx(2,2),fxy(2,2),fyy(2,2) + call get_2derivs2 (axx, fxx, pos, off, dx) + call get_2derivs2 (axy, fxy, pos, off, dx) + call get_2derivs2 (ayy, fyy, pos, off, dx) + f(1,1,:,:) = fxx + f(1,2,:,:) = fxy + f(2,1,:,:) = fxy + f(2,2,:,:) = fyy + end subroutine get_2tensorderivs2 +end module pointwise2 diff --git a/src/rotation.F90 b/src/rotation.F90 new file mode 100644 index 0000000..d76973e --- /dev/null +++ b/src/rotation.F90 @@ -0,0 +1,59 @@ +! $Header$ + +#include "cctk.h" +#include "cctk_Parameters.h" + +module rotation + use constants + implicit none + DECLARE_CCTK_PARAMETERS + private + + public make_euler + +contains + + subroutine make_euler1 (phi, tt) + CCTK_REAL, intent(in) :: phi + CCTK_REAL, intent(out) :: tt(0:3,0:3) + ! y^a = T^a_b x^b + tt = delta4 + tt(1,1) = cos(phi) + tt(1,2) = -sin(phi) + tt(2,1) = sin(phi) + tt(2,2) = cos(phi) + end subroutine make_euler1 + + subroutine make_euler2 (theta, tt) + CCTK_REAL, intent(in) :: theta + CCTK_REAL, intent(out) :: tt(0:3,0:3) + ! y^a = T^a_b x^b + tt = delta4 + tt(2,2) = cos(theta) + tt(2,3) = -sin(theta) + tt(3,2) = sin(theta) + tt(3,3) = cos(theta) + end subroutine make_euler2 + + subroutine make_euler3 (psi, tt) + CCTK_REAL, intent(in) :: psi + CCTK_REAL, intent(out) :: tt(0:3,0:3) + ! y^a = T^a_b x^b + tt = delta4 + tt(1,1) = cos(psi) + tt(1,2) = -sin(psi) + tt(2,1) = sin(psi) + tt(2,2) = cos(psi) + end subroutine make_euler3 + + subroutine make_euler (phi, theta, psi, tt) + CCTK_REAL, intent(in) :: phi, theta, psi + CCTK_REAL, intent(out) :: tt(0:3,0:3) + CCTK_REAL :: tt1(0:3,0:3), tt2(0:3,0:3), tt3(0:3,0:3) + call make_euler1 (phi, tt1) + call make_euler2 (theta, tt2) + call make_euler3 (psi, tt3) + tt = matmul(tt3, matmul(tt2, tt1)) + end subroutine make_euler + +end module rotation diff --git a/src/tensor.F90 b/src/tensor.F90 new file mode 100644 index 0000000..0fef831 --- /dev/null +++ b/src/tensor.F90 @@ -0,0 +1,127 @@ +! $Header$ + +#include "cctk.h" +#include "cctk_Parameters.h" + +module tensor + implicit none + DECLARE_CCTK_PARAMETERS + private + + public calc_trace + + public calc_det + public calc_detderiv + public calc_detdot + + public calc_inv + public calc_invderiv + public calc_invdot + +contains + + subroutine calc_trace (kk, gu, trk) + CCTK_REAL, intent(in) :: kk(3,3) + CCTK_REAL, intent(in) :: gu(3,3) + CCTK_REAL, intent(out) :: trk + integer :: i,j + trk = 0 + do i=1,3 + do j=1,3 + trk = trk + gu(i,j) * kk(i,j) + end do + end do + end subroutine calc_trace + + + + subroutine calc_det (gg, dtg) + CCTK_REAL, intent(in) :: gg(3,3) + CCTK_REAL, intent(out) :: dtg + dtg = gg(1,1) * gg(2,2) * gg(3,3) & + & + 2*gg(1,2) * gg(1,3) * gg(2,3) & + & - gg(1,1) * gg(2,3)**2 & + & - gg(2,2) * gg(1,3)**2 & + & - gg(3,3) * gg(1,2)**2 + end subroutine calc_det + + subroutine calc_pdet (gg, pgg, pdtg) + CCTK_REAL, intent(in) :: gg(3,3), pgg(3,3) + CCTK_REAL, intent(out) :: pdtg + pdtg = pgg(1,1) * gg(2,2) * gg(3,3) & + & + gg(1,1) * pgg(2,2) * gg(3,3) & + & + gg(1,1) * gg(2,2) * pgg(3,3) & + & + 2*pgg(1,2) * gg(1,3) * gg(2,3) & + & + 2* gg(1,2) * pgg(1,3) * gg(2,3) & + & + 2* gg(1,2) * gg(1,3) * pgg(2,3) & + & - pgg(1,1) * gg(2,3)**2 & + & - gg(1,1) * 2*gg(2,3) * pgg(2,3) & + & - pgg(2,2) * gg(1,3)**2 & + & - gg(2,2) * 2*gg(1,3) * pgg(1,3) & + & - pgg(3,3) * gg(1,2)**2 & + & - gg(3,3) * 2*gg(1,2) * pgg(1,2) + end subroutine calc_pdet + + subroutine calc_detderiv (gg, dgg, ddtg) + CCTK_REAL, intent(in) :: gg(3,3), dgg(3,3,3) + CCTK_REAL, intent(out) :: ddtg(3) + integer :: i + do i=1,3 + call calc_pdet (gg, dgg(:,:,i), ddtg(i)) + end do + end subroutine calc_detderiv + + subroutine calc_detdot (gg, gg_dot, dtg_dot) + CCTK_REAL, intent(in) :: gg(3,3), gg_dot(3,3) + CCTK_REAL, intent(out) :: dtg_dot + call calc_pdet (gg, gg_dot, dtg_dot) + end subroutine calc_detdot + + + + subroutine calc_inv (gg, dtg, gu) + CCTK_REAL, intent(in) :: gg(3,3), dtg + CCTK_REAL, intent(out) :: gu(3,3) + gu(1,1) = (gg(2,2) * gg(3,3) - gg(2,3) ** 2 ) / dtg + gu(2,2) = (gg(1,1) * gg(3,3) - gg(1,3) ** 2 ) / dtg + gu(3,3) = (gg(1,1) * gg(2,2) - gg(1,2) ** 2 ) / dtg + gu(1,2) = (gg(1,3) * gg(2,3) - gg(1,2) * gg(3,3)) / dtg + gu(1,3) = (gg(1,2) * gg(2,3) - gg(1,3) * gg(2,2)) / dtg + gu(2,3) = (gg(1,3) * gg(1,2) - gg(2,3) * gg(1,1)) / dtg + gu(2,1) = gu(1,2) + gu(3,1) = gu(1,3) + gu(3,2) = gu(2,3) + end subroutine calc_inv + + subroutine calc_pinv (gu, pgg, pgu) + CCTK_REAL, intent(in) :: gu(3,3), pgg(3,3) + CCTK_REAL, intent(out) :: pgu(3,3) + integer :: i,j,k,l + do i=1,3 + do j=1,3 + pgu(i,j) = 0 + do k=1,3 + do l=1,3 + pgu(i,j) = pgu(i,j) - gu(i,k) * gu(j,l) * pgg(k,l) + end do + end do + end do + end do + end subroutine calc_pinv + + subroutine calc_invderiv (gu, dgg, dgu) + CCTK_REAL, intent(in) :: gu(3,3), dgg(3,3,3) + CCTK_REAL, intent(out) :: dgu(3,3,3) + integer :: i + do i=1,3 + call calc_pinv (gu, dgg(:,:,i), dgu(:,:,i)) + end do + end subroutine calc_invderiv + + subroutine calc_invdot (gu, gg_dot, gu_dot) + CCTK_REAL, intent(in) :: gu(3,3), gg_dot(3,3) + CCTK_REAL, intent(out) :: gu_dot(3,3) + call calc_pinv (gu, gg_dot, gu_dot) + end subroutine calc_invdot + +end module tensor diff --git a/src/tensor2.F90 b/src/tensor2.F90 new file mode 100644 index 0000000..b02cd06 --- /dev/null +++ b/src/tensor2.F90 @@ -0,0 +1,107 @@ +! $Header$ + +#include "cctk.h" +#include "cctk_Parameters.h" + +module tensor2 + implicit none + DECLARE_CCTK_PARAMETERS + private + + public calc_2trace + + public calc_2det + public calc_2detderiv + public calc_2detdot + + public calc_2inv + public calc_2invderiv + public calc_2invdot + +contains + + subroutine calc_2trace (kk, gu, trk) + CCTK_REAL, intent(in) :: kk(2,2) + CCTK_REAL, intent(in) :: gu(2,2) + CCTK_REAL, intent(out) :: trk + integer :: i,j + trk = 0 + do i=1,2 + do j=1,2 + trk = trk + gu(i,j) * kk(i,j) + end do + end do + end subroutine calc_2trace + + + + subroutine calc_2det (gg, dtg) + CCTK_REAL, intent(in) :: gg(2,2) + CCTK_REAL, intent(out) :: dtg + dtg = gg(1,1) * gg(2,2) - gg(1,2)**2 + end subroutine calc_2det + + subroutine calc_2pdet (gg, pgg, pdtg) + CCTK_REAL, intent(in) :: gg(2,2), pgg(2,2) + CCTK_REAL, intent(out) :: pdtg + pdtg = pgg(1,1) * gg(2,2) + gg(1,1) * pgg(2,2) - 2 * gg(1,2) * pgg(1,2) + end subroutine calc_2pdet + + subroutine calc_2detderiv (gg, dgg, ddtg) + CCTK_REAL, intent(in) :: gg(2,2), dgg(2,2,2) + CCTK_REAL, intent(out) :: ddtg(2) + integer :: i + do i=1,2 + call calc_2pdet (gg, dgg(:,:,i), ddtg(i)) + end do + end subroutine calc_2detderiv + + subroutine calc_2detdot (gg, gg_dot, dtg_dot) + CCTK_REAL, intent(in) :: gg(2,2), gg_dot(2,2) + CCTK_REAL, intent(out) :: dtg_dot + call calc_2pdet (gg, gg_dot, dtg_dot) + end subroutine calc_2detdot + + + + subroutine calc_2inv (gg, dtg, gu) + CCTK_REAL, intent(in) :: gg(2,2), dtg + CCTK_REAL, intent(out) :: gu(2,2) + gu(1,1) = gg(2,2) / dtg + gu(2,2) = gg(1,1) / dtg + gu(1,2) = gg(1,2) / dtg + gu(2,1) = gu(1,2) + end subroutine calc_2inv + + subroutine calc_2pinv (gu, pgg, pgu) + CCTK_REAL, intent(in) :: gu(2,2), pgg(2,2) + CCTK_REAL, intent(out) :: pgu(2,2) + integer :: i,j,k,l + do i=1,2 + do j=1,2 + pgu(i,j) = 0 + do k=1,2 + do l=1,2 + pgu(i,j) = pgu(i,j) - gu(i,k) * gu(j,l) * pgg(k,l) + end do + end do + end do + end do + end subroutine calc_2pinv + + subroutine calc_2invderiv (gu, dgg, dgu) + CCTK_REAL, intent(in) :: gu(2,2), dgg(2,2,2) + CCTK_REAL, intent(out) :: dgu(2,2,2) + integer :: i + do i=1,2 + call calc_2pinv (gu, dgg(:,:,i), dgu(:,:,i)) + end do + end subroutine calc_2invderiv + + subroutine calc_2invdot (gu, gg_dot, gu_dot) + CCTK_REAL, intent(in) :: gu(2,2), gg_dot(2,2) + CCTK_REAL, intent(out) :: gu_dot(2,2) + call calc_2pinv (gu, gg_dot, gu_dot) + end subroutine calc_2invdot + +end module tensor2 -- cgit v1.2.3