diff options
| author | allen <allen@5301f0c2-dbc4-4cee-b2f5-8d7afba4d129> | 1999-11-01 14:56:51 +0000 |
|---|---|---|
| committer | allen <allen@5301f0c2-dbc4-4cee-b2f5-8d7afba4d129> | 1999-11-01 14:56:51 +0000 |
| commit | 661c56c847f4334feb53fd8d94551b2ab1efdb9a (patch) | |
| tree | 3dbdeebf2c0b8f0f0fe30218e9a58fc0cd1638a1 | |
| parent | d8cf32e113a96b24744d231013e3d0b81bf55c30 (diff) | |
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.einsteintoolkit.org/cactus/EinsteinAnalysis/Extract/trunk@3 5301f0c2-dbc4-4cee-b2f5-8d7afba4d129
43 files changed, 4428 insertions, 0 deletions
diff --git a/COPYRIGHT b/COPYRIGHT new file mode 100644 index 0000000..45b8425 --- /dev/null +++ b/COPYRIGHT @@ -0,0 +1,377 @@ + ******************************** + Welcome to the Cactus Code + Version 4.0 + ******************************** + +All Cactus source code is (c) Copyright by the Authors listed in the +respective README and source files that you will find scattered thru +the directory structure of all the thorns. This is a truly +collaborative development project that benefits from input of many +people worldwide and we try our best to give credit to every +contribution. The cactus code web page at + http://www.cactuscode.org +contains an extensive list of authors to all the thorns that are +part of the Cactus Code project. + +The flesh of the Cactus Code and the design of the project is (c) +Copyright 1997-99 by Gabrielle Allen, Tom Goodale, Joan Masso and Paul +Walker from the Numerical Relativity group of the Max-Planck-Institut +fuer Gravitationphysik / Albert Einstein Institute +(http://www.aei-potsdam.mpg.de). + +The Cactus Code is distributed under the GNU general public license +version 2 which is available below. All thorns distributed with Cactus +follow this license unless it is explicitly stated otherwise in the +corresponding README files. + +Please, take a look at the documentation in the main doc directory or +check the cactus web page for instructions on how to run the code, +report bugs, give feedback, etc. + +Thanks for using cactus. We hope you find it useful. Love it or hate +it, please let us know what you think. +- The Cactus Team <cactus@cactuscode.org> + +************************************************************************** + + 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. + + <one line to give the program's name and a brief idea of what it does.> + Copyright (C) 19yy <name of author> + + 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. + + <signature of Ty Coon>, 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. + + @@ -0,0 +1,11 @@ +Cactus Code Thorn Extract +Authors : Gabrielle Allen +CVS info : $Header$ +-------------------------------------------------------------------------- + +Purpose of the thorn: + +Calculate gauge invariant waveforms, using the Moncrief formulism, +for cases when the spacetime approximates spherical symmetry + + diff --git a/interface.ccl b/interface.ccl new file mode 100644 index 0000000..6b24006 --- /dev/null +++ b/interface.ccl @@ -0,0 +1,12 @@ +# Interface definition for thorn Extract +# $Header$ + +implements: extract +inherits: grid einstein +private: + +REAL temps TYPE=GF DIM=3 +{ + temp3d + g00 +}
\ No newline at end of file diff --git a/par/Extract.par b/par/Extract.par new file mode 100644 index 0000000..3af0e13 --- /dev/null +++ b/par/Extract.par @@ -0,0 +1,54 @@ +################################################################# +# Extract.par: Extracts from Schwarzschild Black Hole +################################################################# + +ActiveThorns = "pugh cartgrid3d time einstein interp extract idanalyticbh" + +driver::global_nx = 30 +driver::global_ny = 30 +driver::global_nz = 30 + +grid::domain = "octant" +grid::type = "bySpacing" +grid::dxyz = 0.3 + +time::dtfac = 0.2 + +cactus::cctk_itlast = 10 + +einstein::initial_data = "onebh" + +### Extraction parameters + +interp::intorder = 2 + +extract::Extract_num_detectors = 2 +extract::Extract_itout = 1 + +extract::Extract_Nt = 150 +extract::Extract_Np = 150 +extract::Extract_origin_x = 0. +extract::Extract_origin_y = 0. +extract::Extract_origin_z = 0. + +extract::Extract_all_modes = "no" +extract::Extract_l_mode = 2 +extract::Extract_m_mode = 0 +extract::Extract_detector1 = 5.0 +extract::Extract_detector2 = 6.0 + + +### Sample Output (on Origin) +# +# +# Qeven_R1_20 +# +# 0.000000000000000 5.3832644668692326E-05 0.000000000000000 +# 5.9999999999999998E-02 5.3847884014513530E-05 0.000000000000000 +# 0.1200000000000000 5.3893621456197952E-05 0.000000000000000 +# +# mass_R1 +# +# 0.000000000000000 0.9999777739365948 +# 5.9999999999999998E-02 0.9999698763645316 +# 0.1200000000000000 0.9999461867852573 diff --git a/param.ccl b/param.ccl new file mode 100644 index 0000000..428c58a --- /dev/null +++ b/param.ccl @@ -0,0 +1,154 @@ +# Parameter definitions for thorn Extract +# $Header$ + + +shares: grid + +USES KEYWORD domain "" +{ +} + +shares: io + +USES STRING outdir "" +{ +} + +shares: einstein + +USES LOGICAL use_conformal "" +{ +} + +private: + +KEYWORD timecoord "Which time coordinate to use" +{ + "proper" :: "" + "coordinate" :: "" + "both" :: "" +} "both" + + +LOGICAL all_modes "Extract: all l,m modes up to l" +{ +} "yes" + + +LOGICAL Cauchy "Do Cauchy data extraction at given timestep" +{ +} "no" + +LOGICAL verbose "Say what is happening" +{ +} "no" + +LOGICAL doADMmass "Calculate ADM mass at extraction radii" +{ +} "no" + + +LOGICAL do_momentum "Calculate momentum at extraction radii" +{ +} "no" + +LOGICAL do_spin "Calculate spin at extraction radii" +{ +} "no" + +INT itout "How often to extract, in iterations" +{ + 0:* :: "" +} 1 + +INT l_mode "l mode" +{ +0:* :: "" +} 2 + +INT m_mode "m mode (ignore if extracting all modes" +{ +0:* :: "" +} 0 + +INT Nt "Number of theta divisions" +{ +0:* :: "" +} 100 + +INT Np "Number of phi divisions" +{ +0:* :: "" +} 100 + +INT Cauchy_timestep "Timestep for Cauchy data extraction" +{ +0:* :: "" +} 0 + +INT num_detectors "Number of detectors" +{ +0:* :: "" +} 0 + +REAL detector1 "Coordinate radius of detector 1" +{ +0:* :: "" +} 5.0 +REAL detector2 "Coordinate radius of detector 2" +{ +0:* :: "" +} 5.0 +REAL detector3 "Coordinate radius of detector 3" +{ +0:* :: "" +} 5.0 +REAL detector4 "Coordinate radius of detector 4" +{ +0:* :: "" +} 5.0 +REAL detector5 "Coordinate radius of detector 5" +{ +0:* :: "" +} 5.0 +REAL detector6 "Coordinate radius of detector 6" +{ +0:* :: "" +} 5.0 +REAL detector7 "Coordinate radius of detector 7" +{ +0:* :: "" +} 5.0 +REAL detector8 "Coordinate radius of detector 8" +{ +0:* :: "" +} 5.0 +REAL detector9 "Coordinate radius of detector 9" +{ +0:* :: "" +} 5.0 + +REAL origin_x "x-origin to extract about" +{ +*:* :: "" +} 0.0 +REAL origin_y "y-origin to extract about" +{ +*:* :: "" +} 0.0 +REAL origin_z "z-origin to extract about" +{ +*:* :: "" +} 0.0 + +REAL Cauchy_r1 "First radius for Cauchy data extraction" +{ +*:* :: "" +} 1.0 + + +REAL Cauchy_dr "Gridspacing for Cauchy data extraction" +{ +*:* :: "" +} 0.2 + diff --git a/schedule.ccl b/schedule.ccl new file mode 100644 index 0000000..bfb2cbb --- /dev/null +++ b/schedule.ccl @@ -0,0 +1,8 @@ +# Schedule definitions for thorn Extract +# $Header$ + +schedule Extract at POSTSTEP +{ + LANG: Fortran + STOR: temps +} "Extract waveforms" diff --git a/src/ADMmass_integrand3D.F b/src/ADMmass_integrand3D.F new file mode 100644 index 0000000..06b1198 --- /dev/null +++ b/src/ADMmass_integrand3D.F @@ -0,0 +1,170 @@ +#include "cctk.h" + +c ======================================================================== + + SUBROUTINE ADMmass_integrand3D(origin,Dx,Dy,Dz,x,y,z,gxx,gxy,gxz, + & gyy,gyz,gzz,ADMmass_int,Psi,Psi_power) + +c ------------------------------------------------------------------------ +c +c Estimates the ADM mass at a given radius using Equation (7) from +c O Murchadha and York, Phys Rev D, 10, 1974 page 2345 +c +c ------------------------------------------------------------------------ + + IMPLICIT NONE + +c Input variables + INTEGER,INTENT(IN) :: + & Psi_power + CCTK_REAL,INTENT(IN) :: + & Dx,Dy,Dz,origin(3) + CCTK_REAL,DIMENSION(:),INTENT(IN) :: + & x,y,z + CCTK_REAL,DIMENSION(:,:,:),INTENT(IN) :: + & gxx,gxy,gxz,gyy,gyz,gzz,Psi + +c Output variables + CCTK_REAL,DIMENSION(:,:,:),INTENT(OUT) :: + & ADMmass_int + +c Local variables, here only + INTEGER :: + & i,j,k,ip + CCTK_REAL,PARAMETER :: + & half = 0.5D0 + CCTK_REAL :: + & rad,ux,uy,uz,det,dxx,dxy,dxz,dyy,dyz,dzz,uxx,uxy,uxz,uyy, + & uyz,uzz,term1,term2,term3,dxx_y,dxx_z,dxy_x,dxy_y,dxy_z, + & dyy_x,dxz_x,dxz_y,dxz_z,dyz_x,dzz_x,dyy_z,dyz_y,dyz_z,dzz_y, + & Pi,idet,p,pip,pim,pjp,pjm,pkp,pkm + +c ------------------------------------------------------------------------ + + Pi = ACOS(-1D0) + +c Because other codes evolve Psi**4 + SELECT CASE (Psi_power) + + CASE (1) + ip = 4 + + CASE (4) + ip = 1 + + CASE DEFAULT + WRITE(*,*) "This value of Psi_power is not supported" + + END SELECT + + + DO k = 2, SIZE(z)-1 + DO j = 2, SIZE(y)-1 + DO i = 2, SIZE(x)-1 + + rad = SQRT((x(i)-origin(1))**2 + & +(y(j)-origin(2))**2 + & +(z(k)-origin(3))**2) + + IF (rad.NE.0) THEN + + ux = (x(i)-origin(1))/rad + uy = (y(j)-origin(2))/rad + uz = (z(k)-origin(3))/rad + +c Abbreviations for metric coefficients +c ------------------------------------- + p = psi(i,j,k)**ip + + dxx = p*gxx(i,j,k); dxy = p*gxy(i,j,k) + dxz = p*gxz(i,j,k); dyy = p*gyy(i,j,k) + dyz = p*gyz(i,j,k); dzz = p*gzz(i,j,k) + +c Determinant of 3-metric +c ----------------------- + det = (dxx*dyy*dzz + 2.0D0*dxy*dxz*dyz + & - (dxx*dyz**2 + dyy*dxz**2 + dzz*dxy**2)) + + idet = 1.0/det + +c Inverse 3-metric +c ---------------- + uxx = idet*(dyy*dzz - dyz**2) + uyy = idet*(dxx*dzz - dxz**2) + uzz = idet*(dxx*dyy - dxy**2) + uxy = idet*(dxz*dyz - dzz*dxy) + uxz = idet*(dxy*dyz - dyy*dxz) + uyz = idet*(dxy*dxz - dxx*dyz) + +c Derivatives of the 3-metric +c --------------------------- + pip = psi(i+1,j,k)**ip + pim = psi(i-1,j,k)**ip + pjp = psi(i,j+1,k)**ip + pjm = psi(i,j-1,k)**ip + pkp = psi(i,j,k+1)**ip + pkm = psi(i,j,k-1)**ip + + dxx_y = half/Dy*(pjp*gxx(i,j+1,k)-pjm*gxx(i,j-1,k)) + dxx_z = half/Dz*(pkp*gxx(i,j,k+1)-pkm*gxx(i,j,k-1)) + dxy_x = half/Dx*(pip*gxy(i+1,j,k)-pim*gxy(i-1,j,k)) + dxy_y = half/Dy*(pjp*gxy(i,j+1,k)-pjm*gxy(i,j-1,k)) + dxy_z = half/Dz*(pkp*gxy(i,j,k+1)-pkm*gxy(i,j,k-1)) + dyy_x = half/Dx*(pip*gyy(i+1,j,k)-pim*gyy(i-1,j,k)) + dyy_z = half/Dz*(pkp*gyy(i,j,k+1)-pkm*gyy(i,j,k-1)) + dxz_x = half/Dx*(pip*gxz(i+1,j,k)-pim*gxz(i-1,j,k)) + dxz_y = half/Dy*(pjp*gxz(i,j+1,k)-pjm*gxz(i,j-1,k)) + dxz_z = half/Dz*(pkp*gxz(i,j,k+1)-pkm*gxz(i,j,k-1)) + dyz_x = half/Dx*(pip*gyz(i+1,j,k)-pim*gyz(i-1,j,k)) + dyz_y = half/Dy*(pjp*gyz(i,j+1,k)-pjm*gyz(i,j-1,k)) + dyz_z = half/Dz*(pkp*gyz(i,j,k+1)-pkm*gyz(i,j,k-1)) + dzz_x = half/Dx*(pip*gzz(i+1,j,k)-pim*gzz(i-1,j,k)) + dzz_y = half/Dy*(pjp*gzz(i,j+1,k)-pjm*gzz(i,j-1,k)) + + term1 = uxy*(dxx_y-dxy_x)+uxz*(dxx_z-dxz_x) + & +uyy*(dxy_y-dyy_x)+uyz*(dxy_z-dyz_x) + & +uyz*(dxz_y-dyz_x)+uzz*(dxz_z-dzz_x) + + term2 = uyz*(dyy_z-dyz_y)+uxy*(dyy_x-dxy_y) + & +uzz*(dyz_z-dzz_y)+uxz*(dyz_x-dxz_y) + & +uxz*(dxy_z-dxz_y)+uxx*(dxy_x-dxx_y) + + term3 = uxz*(dzz_x-dxz_z)+uyz*(dzz_y-dyz_z) + & +uxx*(dxz_x-dxx_z)+uxy*(dxz_y-dxy_z) + & +uxy*(dyz_x-dxy_z)+uyy*(dyz_y-dyy_z) + + ADMmass_int(i,j,k) = 1.0D0/16.0D0/Pi*(ux*term1+ + & uy*term2+uz*term3)*SQRT(det)*rad**2 + + ELSE + + ADMmass_int(i,j,k) = 0.0D0 + + ENDIF + + ENDDO + ENDDO + ENDDO + + +c This is needed when the grid is an octant, but it does not hurt +c if it is not + + DO k = 2, size(z)-1 + DO j = 2, size(y)-1 + ADMmass_int(1,j,k) = ADMmass_int(2,j,k) + ENDDO + ENDDO + DO k = 2, size(z)-1 + DO i = 1, size(x)-1 + ADMmass_int(i,1,k) = ADMmass_int(i,2,k) + ENDDO + ENDDO + DO j = 1, size(y)-1 + DO i = 1, size(x)-1 + ADMmass_int(i,j,1) = ADMmass_int(i,j,2) + ENDDO + ENDDO + + + END SUBROUTINE ADMmass_integrand3D diff --git a/src/ADMmass_integrand3D_int.F b/src/ADMmass_integrand3D_int.F new file mode 100644 index 0000000..c103f55 --- /dev/null +++ b/src/ADMmass_integrand3D_int.F @@ -0,0 +1,27 @@ + +#include "cctk.h" + + MODULE ADMmass_integrand3D_int + +c ------------------------------------------------------------------ + + INTERFACE + + SUBROUTINE ADMmass_integrand3D(origin,Dx,Dy,Dz,x,y,z,gxx,gxy, + & gxz,gyy,gyz,gzz,ADMmass_int,Psi,Psi_power) + IMPLICIT NONE + INTEGER,INTENT(IN) :: + & Psi_power + CCTK_REAL,INTENT(IN) :: + & Dx,Dy,Dz,origin(3) + CCTK_REAL,DIMENSION(:),INTENT(IN) :: + & x,y,z + CCTK_REAL,DIMENSION(:,:,:),INTENT(IN) :: + & gxx,gxy,gxz,gyy,gyz,gzz,Psi + CCTK_REAL,DIMENSION(:,:,:),INTENT(OUT) :: + & ADMmass_int + END SUBROUTINE ADMmass_integrand3D + + END INTERFACE + + END MODULE ADMmass_integrand3D_int diff --git a/src/D2_extract.F b/src/D2_extract.F new file mode 100644 index 0000000..3b24ae0 --- /dev/null +++ b/src/D2_extract.F @@ -0,0 +1,918 @@ + +#include "cctk.h" + +c ======================================================================== + + SUBROUTINE D2_extract(eta,igrid,Nt,Np,theta,phi,all_modes,l,m,g00,g11, + & g12,g13,g22,g23,g33,dg22,dg23,dg33, + & do_ADMmass,ADMmass_int1,ADMmass_int2, + & do_momentum,momentum_int1,momentum_int2,momentum_int3, + & do_spin,spin_int1,spin_int2,spin_int3, + & mass,rsch,Qodd,Qeven,ADMmass,momentum,spin,dtaudt) + +c ------------------------------------------------------------------------ +c +c Extract the odd and even parity gauge invariant variables +c at a given isotropic radius (eta) from the center of a +c slightly perturbed Schwarzschild spacetime. +c +c The polar coordinates and metric variables must be given on +c an equally spaced grid, which is offset from the axis. +c +c Notice that for octant symmetry we only need calculate the real +c parts of the gauge invariant variables for the even-l,m modes. +c This subroutine is optimised (!!) for octant symmetry. +c +c ------------------------------------------------------------------------ + + IMPLICIT NONE + +c Input variables + INTEGER,INTENT(IN) :: + & igrid,l,m + CCTK_INT,INTENT(IN) :: + & Nt,Np,all_modes,do_momentum,do_spin + INTEGER,INTENT(IN) :: + & do_ADMmass(2) + CCTK_REAL,INTENT(IN) :: + & eta,theta(:),phi(:) + CCTK_REAL,DIMENSION (:,:) :: + & g00,g11,g12,g13,g22,g23,g33,dg22,dg23,dg33, + & ADMmass_int1,ADMmass_int2, + & momentum_int1,momentum_int2,momentum_int3, + & spin_int1,spin_int2,spin_int3 + +c Output variables + CCTK_REAL,INTENT(OUT) :: + & dtaudt,ADMmass(2),mass,rsch,Qodd(:,2:,0:),Qeven(:,2:,0:), + & momentum(3),spin(3) + +c Local Variables + LOGICAL, PARAMETER :: + & debug = .FALSE. + INTEGER :: + & i,j,il,im,lmin,lmax,mmin,mmax,lstep,mstep + CCTK_REAL,PARAMETER :: + & zero = 0.0D0, + & half = 0.5D0, + & one = 1.0D0, + & two = 2.0D0, + & four = 4.0D0, + & six = 6.0D0, + & eight = 8.0D0 + CCTK_REAL :: + & Dt,Dp,st,ist,drschdri,dridrsch, + & g22comb,g23comb,g33comb,S,rsch2,Lambda, + & Pi, + & fac_h1,fac_H2,fac_K,fac_G,fac_c1,fac_c2, + & fac_dG,fac_dK,fac_dc2,fac1,fac2,fac3,max_g11, + & sph_g11,sph_g22,sph_g33,sph_dg22,sph_dg33,error, + & sphere_int + CCTK_REAL,DIMENSION(2) :: + & Y,Y1,Y2,Y3,Y4,h1,H2,K,G,c1,c2,dG,dK,dc2 + CCTK_REAL,DIMENSION(Nt+1,Np+1) :: + & temp,h1i,H2i,Gi,Ki,c1i,c2i,dGi,dKi,dc2i + + +c ------------------------------------------------------------------------ +c +c 0. Initial things +c +c ------------------------------------------------------------------------ + + Pi = two*ASIN(one) + + SELECT CASE (igrid) + CASE (0) ! full grid + Dt = Pi/DBLE(Nt) + Dp = two*Pi/DBLE(Np) + lstep = 1 + mstep = 1 + CASE (1) ! octant grid + Dt = half*Pi/DBLE(Nt) + Dp = half*Pi/DBLE(Np) + lstep = 2 + mstep = 2 + CASE (2) ! cartoon + Dt = Pi/DBLE(Nt) + Dp = 2D0*Pi + lstep = 2 + mstep = 2 + CASE (3) ! bitant + Dt = half*Pi/DBLE(Nt) + Dp = two*Pi/DBLE(Np) + lstep = 1 + mstep = 1 + CASE DEFAULT + WRITE (*,*) "Grid incorrectly set in D2_extract" + STOP + END SELECT + +c Calculate ADM mass +c ------------------ + IF (do_ADMmass(1) == 1) THEN + +c Calculate ADM mass using standard formula + DO j = 2, Np+1 + DO i = 2, Nt+1 + temp(i,j) = ADMmass_int1(i-1,j-1)*DSIN(theta(i-1)) + ENDDO + ENDDO + ADMmass(1) = sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp) + + ENDIF + IF (do_ADMmass(2) == 1) THEN + +c Calculate ADM mass using conformal formula + DO j = 2, Np+1 + DO i = 2, Nt+1 + temp(i,j) = ADMmass_int2(i-1,j-1)*DSIN(theta(i-1)) + ENDDO + ENDDO + ADMmass(2) = sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp) + + END IF + +c Calculate momentum +c ------------------ + IF (do_momentum == 1) THEN + + DO j = 2, Np+1 + DO i = 2, Nt+1 + temp(i,j) = momentum_int1(i-1,j-1)*DSIN(theta(i-1)) + ENDDO + ENDDO + momentum(1) = sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp) + + DO j = 2, Np+1 + DO i = 2, Nt+1 + temp(i,j) = momentum_int2(i-1,j-1)*DSIN(theta(i-1)) + ENDDO + ENDDO + momentum(2) = sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp) + + DO j = 2, Np+1 + DO i = 2, Nt+1 + temp(i,j) = momentum_int3(i-1,j-1)*DSIN(theta(i-1)) + ENDDO + ENDDO + momentum(3) = sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp) + + ENDIF + +c Calculate spin +c -------------- + IF (do_spin == 1) THEN + + DO j = 2, Np+1 + DO i = 2, Nt+1 + temp(i,j) = spin_int1(i-1,j-1)*DSIN(theta(i-1)) + ENDDO + ENDDO + spin(1) = sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp) + + DO j = 2, Np+1 + DO i = 2, Nt+1 + temp(i,j) = spin_int2(i-1,j-1)*DSIN(theta(i-1)) + ENDDO + ENDDO + spin(2) = sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp) + + DO j = 2, Np+1 + DO i = 2, Nt+1 + temp(i,j) = spin_int3(i-1,j-1)*DSIN(theta(i-1)) + ENDDO + ENDDO + spin(3) = sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp) + + ENDIF + +c ------------------------------------------------------------------ +c +c 1. Find Spherical Parts of Metric Components +c +c ------------------------------------------------------------------ + +c ... g00 + DO j = 2, Np+1 + DO i = 2, Nt+1 + temp(i,j) = g00(i-1,j-1)*DSIN(theta(i-1)) + ENDDO + ENDDO + dtaudt = sqrt(one/(four*Pi)*sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp)) + +c ... g11 + DO j = 2, Np+1 + DO i = 2, Nt+1 + temp(i,j) = g11(i-1,j-1)*DSIN(theta(i-1)) + ENDDO + ENDDO + sph_g11 = one/(four*Pi)*sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp) + +c ... g22 + DO j = 2, Np+1 + DO i = 2, Nt+1 + temp(i,j) = g22(i-1,j-1)*DSIN(theta(i-1)) + ENDDO + ENDDO + sph_g22 = one/(four*Pi)*sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp) + +c ... dg22 + DO j = 2, Np+1 + DO i = 2, Nt+1 + temp(i,j) = dg22(i-1,j-1)*DSIN(theta(i-1)) + ENDDO + ENDDO + sph_dg22 = one/(four*Pi)*sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp) + + +c ... g33/sin(theta)**2 + DO j = 2, Np+1 + DO i = 2, Nt+1 + temp(i,j) = g33(i-1,j-1)/DSIN(theta(i-1)) + ENDDO + ENDDO + sph_g33 = one/(four*Pi)*sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp) + + +c Calculate error in orthonormality of spherical harmonic +c DO j = 2, Np+1 +c DO i = 2, Nt+1 +c CALL spher_harm_combs(theta(i),phi(j),l,m,Y,Y1,Y2,Y3,Y4) +c temp(i,j) = (Y(1)*Y(1)+Y(2)*Y(2))*sin(theta(i-1)) +c ENDDO +c ENDDO +c error = sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp) +c print *,"Error is",one-error + + + +c ------------------------------------------------------------------ +c +c 1. Compute the Schwarzschild radius +c +c ------------------------------------------------------------------ + + DO j = 2, Np+1 + DO i = 2, Nt+1 +c temp(i,j) = g22(i-1,j-1)*SIN(theta(i-1)) + temp(i,j) = (g22(i-1,j-1)+g33(i-1,j-1)/SIN(theta(i-1))**2) + & *SIN(theta(i-1))/(two) + ENDDO + ENDDO + + rsch2= one/(four*Pi)*sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp) + + rsch = DSQRT(rsch2) + + +c ------------------------------------------------------------------ +c +c 2. Compute the derivative of the schwarzschild radius with +c respect to the isotropic radius eta. +c +c ------------------------------------------------------------------ + + DO j = 2, Np+1 + DO i = 2, Nt+1 +c temp(i,j) = (dg22(i-1,j-1))*SIN(theta(i-1)) + temp(i,j) = (dg22(i-1,j-1)+dg33(i-1,j-1)/SIN(theta(i-1))**2) + & *SIN(theta(i-1))/(two) + ENDDO + ENDDO + + drschdri = one/(eight*Pi*rsch)* + & sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp) + + dridrsch = one/drschdri + + +c ------------------------------------------------------------------ +c +c 3. Calculate the Schwarzschild mass and S factor +c +c ------------------------------------------------------------------ + + S = drschdri**2/sph_g11 + mass = rsch*(one-S)/two + + + +c ------------------------------------------------------------------ +c +c 4. Subtract spherical part from metric +c Do not know how important this is +c +c ------------------------------------------------------------------ + + DO i=1,Nt + DO j=1,Np + st = SIN(theta(i)) + g11(i,j) = g11(i,j) - sph_g11 + g22(i,j) = g22(i,j) - sph_g22 + g33(i,j) = g33(i,j) - sph_g22*st**2 + dg22(i,j) = dg22(i,j) - sph_dg22 + dg33(i,j) = dg33(i,j) - sph_dg22*st**2 + ENDDO + ENDDO + + + IF (all_modes == 0) THEN + lmin = l ; lmax = l + ELSE + lmin = 2 ; lmax = l + ENDIF + + loop_l: DO il = lmin,lmax,lstep + + IF (all_modes == 0) THEN + mmin = m ; mmax = m + ELSE + mmin = 0 ; mmax = il + END IF + + loop_m: DO im = mmin,mmax,mstep + +c ------------------------------------------------------------------ +c +c 5. Calculate all Regge-Wheeler variables +c +c ------------------------------------------------------------------ + +c Factors independent of angular coordinates +c ------------------------------------------ + fac_h1 = dridrsch/DBLE(il*(il+1)) + fac_H2 = S*dridrsch**2 + fac_G = one/(rsch**2*DBLE(il*(il+1)*(il-1)*(il+2))) + fac_K = half/rsch**2 + fac_c1 = dridrsch/DBLE(il*(il+1)) + fac_c2 = two/DBLE(il*(il+1)*(il-1)*(il+2)) + fac_dG = fac_G + fac_dK = fac_K + fac_dc2 = fac_c2 + + +c ------------------------------------------------------------------ +c +c 5a. Real parts of the Regge-Wheeler variables +c +c ------------------------------------------------------------------ + +c Calculate integrands + + loop_phi1: DO j = 1, Np + loop_theta1: DO i = 1, Nt + + st = SIN(theta(i)) + ist = one/st + + CALL spher_harm_combs(theta(i),phi(j),il,im,Y,Y1,Y2,Y3 + & ,Y4) + + h1i(i+1,j+1) = st*g12(i,j)*Y1(1)+ist*g13(i,j)*Y2(1) + H2i(i+1,j+1) = st*g11(i,j)*Y(1) + Gi(i+1,j+1) = (st*g22(i,j)-ist*g33(i,j))*Y3(1) + & +four*ist*g23(i,j)*Y4(1) + Ki(i+1,j+1) = (st*g22(i,j)+ist*g33(i,j))*Y(1) + c1i(i+1,j+1) = g13(i,j)*Y1(1)-g12(i,j)*Y2(1) + c2i(i+1,j+1) = (g22(i,j)-ist**2*g33(i,j))*Y4(1) + & -g23(i,j)*Y3(1) + + g22comb = dridrsch*dg22(i,j)-two/rsch*g22(i,j) + g23comb = dridrsch*dg23(i,j)-two/rsch*g23(i,j) + g33comb = dridrsch*dg33(i,j)-two/rsch*g33(i,j) + + dGi(i+1,j+1) = (st*g22comb-ist*g33comb)*Y3(1) + & +four*ist*g23comb*Y4(1) + dKi(i+1,j+1) = (st*g22comb+ist*g33comb)*Y(1) + dc2i(i+1,j+1) = (dg22(i,j)-ist**2*dg33(i,j))*Y4(1) + & -dg23(i,j)*Y3(1) + + END DO loop_theta1 + END DO loop_phi1 + +c Integrations over the 2-sphere + + h1(1) = fac_h1*sphere_int(h1i,igrid,Nt+1,Np+1,Dt,Dp) + H2(1) = fac_H2*sphere_int(H2i,igrid,Nt+1,Np+1,Dt,Dp) + G(1) = fac_G*sphere_int(Gi,igrid,Nt+1,Np+1,Dt,Dp) + K(1) = fac_K*sphere_int(Ki,igrid,Nt+1,Np+1,Dt,Dp) + & +DBLE(il*(il+1))*half*G(1) + c1(1) = fac_c1*sphere_int(c1i,igrid,Nt+1,Np+1,Dt,Dp) + c2(1) = fac_c2*sphere_int(c2i,igrid,Nt+1,Np+1,Dt,Dp) + dG(1) = fac_dG*sphere_int(dGi,igrid,Nt+1,Np+1,Dt,Dp) + dK(1) = fac_dK*sphere_int(dKi,igrid,Nt+1,Np+1,Dt,Dp) + & +DBLE(il*(il+1))*half*dG(1) + dc2(1) = fac_dc2*sphere_int(dc2i,igrid,Nt+1,Np+1,Dt,Dp) + + +c ------------------------------------------------------------------ +c +c 5b. Imaginary parts of the Regge-Wheeler variables +c +c ------------------------------------------------------------------ + + complex_parts: IF (igrid /= 1 .AND. igrid /=2) THEN + +c Calculate integrands + + loop_phi2: DO j = 1, Np + + loop_theta2: DO i = 1, Nt + + st = SIN(theta(i)) + ist = one/st + + CALL spher_harm_combs(theta(i),phi(j),il,im,Y,Y1,Y2 + & ,Y3,Y4) + + h1i(i+1,j+1) =-st*g12(i,j)*Y1(2)-ist*g13(i,j)*Y2(2) + H2i(i+1,j+1) = -st*g11(i,j)*Y(2) + Gi(i+1,j+1) = -(st*g22(i,j)-ist*g33(i,j))*Y3(2) + & -four*ist*g23(i,j)*Y4(2) + Ki(i+1,j+1) = -(st*g22(i,j)+ist*g33(i,j))*Y(2) + c1i(i+1,j+1) = -g13(i,j)*Y1(2)+g12(i,j)*Y2(2) + c2i(i+1,j+1) = -(g22(i,j)-ist**2*g33(i,j))*Y4(2) + & +g23(i,j)*Y3(2) + + g22comb = dridrsch*dg22(i,j)-two/rsch*g22(i,j) + g23comb = dridrsch*dg23(i,j)-two/rsch*g23(i,j) + g33comb = dridrsch*dg33(i,j)-two/rsch*g33(i,j) + + dGi(i+1,j+1) = -(st*g22comb-ist*g33comb)*Y3(2) + & -four*ist*g23comb*Y4(2) + dKi(i+1,j+1) = -(st*g22comb+ist*g33comb)*Y(2) + dc2i(i+1,j+1) = -(dg22(i,j)-ist**2*dg33(i,j))*Y4(2) + & +dg23(i,j)*Y3(2) + + END DO loop_theta2 + END DO loop_phi2 + + +c Integrate over 2-sphere + + h1(2) = fac_h1*sphere_int(h1i,igrid,Nt+1,Np+1,Dt,Dp) + H2(2) = fac_H2*sphere_int(H2i,igrid,Nt+1,Np+1,Dt,Dp) + G(2) = fac_G*sphere_int(Gi,igrid,Nt+1,Np+1,Dt,Dp) + K(2) = fac_K*sphere_int(Ki,igrid,Nt+1,Np+1,Dt,Dp) + & +DBLE(il*(il+1))*half*G(2) + c1(2) = fac_c1*sphere_int(c1i,igrid,Nt+1,Np+1,Dt,Dp) + c2(2) = fac_c2*sphere_int(c2i,igrid,Nt+1,Np+1,Dt,Dp) + dG(2) = fac_dG*sphere_int(dGi,igrid,Nt+1,Np+1,Dt,Dp) + dK(2) = fac_dK*sphere_int(dKi,igrid,Nt+1,Np+1,Dt,Dp) + & +DBLE(il*(il+1))*half*dG(2) + dc2(2) = fac_dc2*sphere_int(dc2i,igrid,Nt+1,Np+1,Dt,Dp) + + ELSE + + h1(2) = zero + H2(2) = zero + G(2) = zero + K(2) = zero + c1(2) = zero + c2(2) = zero + dG(2) = zero + dK(2) = zero + dc2(2) = zero + + + END IF complex_parts + + + +c ------------------------------------------------------------------ +c +c 6. Construct gauge invariant variables +c +c ------------------------------------------------------------------ + + Qodd(:,il,im) = SQRT(two*DBLE((il+2)*(il+1)*il*(il-1)))*S/ + & rsch*(c1+half*(dc2-two/rsch*c2)) + + Lambda = DBLE((il-1)*(il+2))+3.0D0*(one-S) + + Qeven(:,il,im) = one/Lambda*SQRT(two*DBLE((il-1)*(il+2))/ + & DBLE(il*(il+1)))*(DBLE(il*(il+1))*S*(rsch**2*dG- + & two*h1)+two*rsch*S*(H2-rsch*dK) + & +Lambda*rsch*K) + + + END DO loop_m + + END DO loop_l + + + +c ------------------------------------------------------------------ +c +c 7. Write some diagnostics if required (for last l,m to be used) +c +c ------------------------------------------------------------------ + + IF (debug) THEN + WRITE(101,*) eta,h1 + WRITE(102,*) eta,H2 + WRITE(103,*) eta,G + WRITE(104,*) eta,K + WRITE(105,*) eta,c1 + WRITE(106,*) eta,c2 + WRITE(107,*) eta,dG + WRITE(108,*) eta,dK + WRITE(109,*) eta,dc2 + WRITE(201,*) eta,mass + WRITE(202,*) eta,rsch + WRITE(203,*) eta,drschdri + WRITE(301,*) eta,Qeven(:,lmax,mmax) + WRITE(302,*) eta,Qodd(:,lmax,mmax) + ENDIF + + + END SUBROUTINE D2_extract + +c =============================================================== +c +c +c +c +c +c =============================================================== + + SUBROUTINE simpsons(n,Dx,f,intf) + +c Simpsons rule to evaluate 1D integral equation + +c Variables In: +c n number of points [must be odd] +c Dx step size +c f(n) integrand at each abscissa + +c Variables Out: +c intf evaluated integral + +c Variables Local: +c i index + +c --------------------------------------------------------------- + + IMPLICIT NONE + +c Input variables + INTEGER,INTENT(IN) :: + & n + CCTK_REAL,INTENT(IN) :: + & Dx,f(n) + +c Output variables + CCTK_REAL,INTENT(OUT) :: + & intf + +c Local variables + INTEGER :: + & i + CCTK_REAL,PARAMETER :: + & two = 2.0D0, + & three = 3.0D0, + & four = 4.0D0 + +c --------------------------------------------------------------- + + IF (MOD(n,2) == 0) THEN + WRITE(*,*) "Call to -simpsons- with even number of points" + STOP + END IF + + intf = f(1) + four*f(n-1) + f(n) + + DO i = 2, n-3,2 + + intf = intf + four*f(i) + two*f(i+1) + + ENDDO + + intf = intf*Dx/three + + END SUBROUTINE simpsons + +c =============================================================== +c +c +c +c +c +c =============================================================== + + FUNCTION sphere_int(temp,igrid,Nt,Np,Dt,Dp) + +c --------------------------------------------------------------- +c +c Integration the function held in temp over the sphere. +c If the full grid is being used then Simpsons rule is used, +c if an octant is used then the symmetry means Simpsons rule +c becomes a simple sum. +c +c --------------------------------------------------------------- + + IMPLICIT NONE + +c Input variables + INTEGER,INTENT(IN) :: + & Nt,Np,igrid + CCTK_REAL,INTENT(IN) :: + & temp(Nt,Np),Dt,Dp + +c Output variables + CCTK_REAL :: + & sphere_int + +c Local variables + INTEGER :: + & i,j,Np_start + CCTK_REAL :: + & ft(Nt),fp(Np) + CCTK_REAL,PARAMETER :: + & two = 2.0D0, + & four = 4.0D0 + +c --------------------------------------------------------------- + + DO j = 2, Np + +c Integrate with respect to theta + + DO i = 2, Nt + ft(i) = temp(i,j) + ENDDO + ft(1) = ft(2) ! from symmetry across axis + + SELECT CASE (igrid) + CASE (0) + CALL simpsons(Nt,Dt,ft,fp(j)) + CASE (1) + fp(j)=two*Dt*SUM(ft(2:Nt)) + CASE (2) + CALL simpsons(Nt,Dt,ft,fp(j)) + CASE (3) + fp(j)=two*Dt*SUM(ft(2:Nt)) + END SELECT + + ENDDO + +c Integrate with respect to phi + + fp(1) = fp(2) ! from symmetry across axis + + + SELECT CASE (igrid) + CASE (0) + CALL simpsons(Np,Dp,fp,sphere_int) + CASE (1) + sphere_int=four*Dp*SUM(fp(2:Np)) + CASE (2) + sphere_int=Dp*fp(1) + CASE (3) + CALL simpsons(Np,Dp,fp,sphere_int) + END SELECT + + + END FUNCTION sphere_int + +c ================================================================== +c +c +c +c +c +c ================================================================== + + SUBROUTINE spherical_harmonic(l,m,theta,phi,Ylm) + +c ------------------------------------------------------------------ +c +c Calculate the (l,m) spherical harmonic at given angular +c coordinates. This number is in general complex, and +c +c Ylm = Ylm(1) + i Ylm(2) +c +c where +c +c a ( 2 l + 1 (l-|m|)! ) i m phi +c Ylm = (-1) SQRT( ------- ------- ) P_l|m| (cos(theta)) e +c ( 4 Pi (l+|m|)! ) +c +c and where +c +c a = m/2 (sign(m)+1) +c +c ------------------------------------------------------------------ + + IMPLICIT NONE + +c Input variables + INTEGER :: + & l,m + CCTK_REAL :: + & theta,phi + +c Output variables + CCTK_REAL :: + & Ylm(2) + +c Local variables + INTEGER :: + & i + CCTK_REAL,PARAMETER :: + & one = 1.0D0, + & two = 2.0D0, + & four = 4.0D0 + CCTK_REAL :: + & a,Pi,fac,plgndr + +c ------------------------------------------------------------------ + + Pi = ACOS(-one) + + fac = one + DO i = l-ABS(m)+1,l+ABS(m) + fac = fac*DBLE(i) + ENDDO + fac = one/fac + +c a = (-one)**((m*ISIGN(m,1)/ABS(m)+m)/2)*SQRT(DBLE(2*l+1) +c & /four/Pi*fac)*plgndr(l,ABS(m),cos(theta)) + + a = (-one)**MAX(m,0)*SQRT(DBLE(2*l+1) + & /four/Pi*fac)*plgndr(l,ABS(m),cos(theta)) + Ylm(1) = a*COS(DBLE(m)*phi) + Ylm(2) = a*SIN(DBLE(m)*phi) + + END SUBROUTINE spherical_harmonic + +c ================================================================== +c +c +c +c +c +c ================================================================== + + FUNCTION plgndr(l,m,x) + +c ------------------------------------------------------------------ +c +c From Numerical Recipes +c +c Calculates the associated Legendre polynomial Plm(x). +c Here m and l are integers satisfying 0 <= m <= l, +c while x lies in the range -1 <= x <= 1 +c +c ------------------------------------------------------------------ + +c Input variables + INTEGER,INTENT(IN) :: + & l,m + CCTK_REAL,INTENT(IN) :: + & x + +c Output variables + CCTK_REAL :: + & plgndr + +c Local Variables + INTEGER :: + & i,ll + CCTK_REAL,PARAMETER :: + & one = 1.0D0, + & two = 2.0D0 + CCTK_REAL :: + & pmm,somx2,fact,pmmp1,pll + +c ------------------------------------------------------------------ + + pmm = one + + IF (m.GT.0) THEN + somx2=SQRT((one-x)*(one+x)) + fact=one + DO i=1,m + pmm = -pmm*fact*somx2 + fact = fact+two + END DO + ENDIF + + IF (l.EQ.m) THEN + plgndr = pmm + ELSE + pmmp1 = x*(two*m+one)*pmm + IF (l.EQ.m+1) THEN + plgndr=pmmp1 + ELSE + DO ll=m+2,l + pll = (x*DBLE(2*ll-1)*pmmp1- + & DBLE(ll+m-1)*pmm)/DBLE(ll-m) + pmm = pmmp1 + pmmp1 = pll + END DO + plgndr = pll + ENDIF + ENDIF + + END FUNCTION plgndr + +c ================================================================== +c +c +c +c +c +c ================================================================== + + SUBROUTINE spher_harm_combs(theta,phi,l,m,Y,Y1,Y2,Y3,Y4) +c ------------------------------------------------------------------ +c +c Calculates the various combinations of spherical harmonics needed +c for the extraction (all are complex): +c +c Y = Ylm +c Y1 = Ylm,theta +c Y2 = Ylm,phi +c Y3 = Ylm,theta,theta-cot theta Ylm,theta-Ylm,phi,phi/sin^2 theta +c Y4 = Ylm,theta,phi-cot theta Ylm,phi +c +c The local variables Yplus is the spherical harmonic at (l+1,m) +c +c ------------------------------------------------------------------ + + IMPLICIT NONE + +c Input variables + INTEGER :: + & l,m + CCTK_REAL :: + & theta,phi + +c Output variables + CCTK_REAL,DIMENSION(2) :: + & Y,Y1,Y2,Y3,Y4 + +c Local variables + INTEGER :: + & i + CCTK_REAL :: + & Yplus(2),rl,rm,cot_theta,Pi + CCTK_REAL,PARAMETER :: + & half = 0.5D0, + & one = 1.0D0, + & two = 2.0D0 + +c ------------------------------------------------------------------ + + Pi = DACOS(-one) + + rl = DBLE(l) + rm = DBLE(m) + + cot_theta = DCOS(theta)/DSIN(theta) + + CALL spherical_harmonic(l+1,m,theta,phi,Yplus) + + +c Find Y ... + + CALL spherical_harmonic(l,m,theta,phi,Y) + + +c Find Y1 ... + + DO i = 1,2 + Y1(i) = -(rl+one)*cot_theta*Y(i)+Yplus(i)/DSIN(theta)* + & DSQRT(((rl+one)**2-rm**2)*(rl+half)/(rl+one+half)) + ENDDO + + +c Find Y2 ... + + Y2(1) = -rm*Y(2) + Y2(2) = rm*Y(1) + + +c Find Y3 ... + + DO i = 1,2 + Y3(i) = -two*cot_theta*Y1(i)+(two*rm*rm/(DSIN(theta)**2) + & -rl*(rl+one))*Y(i) + ENDDO + + +c Find Y4 ... + + Y4(1) = rm*(cot_theta*Y(2)-Y1(2)) + Y4(2) = rm*(Y1(1)-cot_theta*Y(1)) + + + END SUBROUTINE spher_harm_combs + +c ================================================================== + + + + diff --git a/src/D2_extract_int.F b/src/D2_extract_int.F new file mode 100644 index 0000000..43f6906 --- /dev/null +++ b/src/D2_extract_int.F @@ -0,0 +1,42 @@ + +#include "cctk.h" + + MODULE D2_extract_int + +c ------------------------------------------------------------------ + + INTERFACE + + SUBROUTINE D2_extract(eta,igrid,Nt,Np,theta,phi,all_modes,l,m,g00,g11, + & g12,g13,g22,g23,g33,dg22,dg23,dg33, + & do_ADMmass,ADMmass_int1,ADMmass_int2, + & do_momentum,momentum_int1,momentum_int2,momentum_int3, + & do_spin,spin_int1,spin_int2,spin_int3, + & mass,rsch,Qodd,Qeven,ADMmass,momentum,spin,dtaudt) + + IMPLICIT NONE + INTEGER,INTENT(IN) :: + & igrid,l,m + CCTK_INT,INTENT(IN) :: + & Nt,Np,all_modes,do_momentum,do_spin + INTEGER,INTENT(IN) :: + & do_ADMmass(2) + CCTK_REAL,INTENT(IN) :: + & eta + CCTK_REAL,INTENT(IN),DIMENSION (:) :: + & theta,phi + CCTK_REAL,INTENT(IN),DIMENSION (:,:) :: + & g00,g11,g12,g13,g22,g23,g33,dg22,dg23,dg33, + & ADMmass_int1, ADMmass_int2, + & momentum_int1,momentum_int2,momentum_int3, + & spin_int1,spin_int2,spin_int3 + CCTK_REAL,INTENT(OUT) :: + & ADMmass(2),mass,rsch,Qodd(:,:,:),Qeven(:,:,:),dtaudt, + & momentum(3),spin(3) + END SUBROUTINE + + END INTERFACE + + END MODULE D2_extract_int + + diff --git a/src/D3_extract.F b/src/D3_extract.F new file mode 100644 index 0000000..f5f6850 --- /dev/null +++ b/src/D3_extract.F @@ -0,0 +1,221 @@ + +#include "cctk.h" + +c ================================================================== + + SUBROUTINE D3_extract(cctkGH,do_ADMmass,do_momentum,do_spin,igrid, + & origin,myproc,Nt,Np,all_modes, + & l,m,x,y,z,Dx,Dy,Dz,Psi_power,Psi, + & g00,gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz, + & eta,ADMmass,momentum,spin,mass,rsch,Qodd,Qeven,Extract_temp3d,dtaudt) + +c ------------------------------------------------------------------ +c +c Extract odd and even parity gauge invariant variables on a +c 2-surface defined by constant isotropic radius ETA from ORIGIN, +c also returns the corresponding Schwarzschild radius RSCH. +c +c Works for a full grid (IGRID=0) or an octant (IGRID=1), given +c non-conformal (PSI=1,Gij) or conformal (PSI^PSI_POWER,Gij). +c +c IMPORTANT : Nt and Np must both be even numbers (so that Simpsons +c rule works later) +c +c Variables in: +c ------------ +c igrid Grid type (0=full,1=octant) +c [Could extract in an octant from a full grid, but why +c would you want to] +c origin The position of the origin of extraction in the x,y,z +c coordinate system +c myproc Cactus variable, must be zero if not using Cactus +c Nt,Np Number of theta,phi grid DIVISIONS to use [even numbers] +c all_modes if true the extract all l,m modes up to l +c l,m Spherical harmonic to extract if ALL_MODES is false. +c otherwise m is ignored and l is the maximum l mode +c x,y,z Cartesian coordinates in the 3-space, hopefully +c isotropic about ORIGIN +c Dx,Dy,Dz Grid spacings +c Psi_power The power of Psi which is being passed in, usually one +c but may be four +c Psi The conformal factor, set it to one if non-conformal +c data is passed in +c gij The 3-metric components in cartesian coordinates +c eta The coordinate radius from ORIGIN for extraction, must +c obviously give a 2-surface which is contained in the +c grid +c +c Variables out: +c ------------- +c ADMmass Estimate of ADM mass +c momentum Estimate of momentum +c spin Estimate of spin +c mass The extracted mass +c rsch The extracted Schwarzschild radius +c Qodd The extracted odd parity gauge invariant variable +c Qeven The extracted even parity gauge invariant variable +c +c ------------------------------------------------------------------ + + USE D3_to_D2_int + USE cartesian_to_spherical_int + USE unphysical_to_physical_int + USE D2_extract_int + +c ------------------------------------------------------------------ + + IMPLICIT NONE + +c Input variables + + CCTK_POINTER :: cctkGH + + INTEGER,INTENT(IN) :: + & igrid,l,m,Psi_power,myproc + CCTK_INT,INTENT(IN) :: + & Nt,Np,all_modes,do_momentum,do_spin + INTEGER,INTENT(IN) :: + & do_ADMmass(2) + CCTK_REAL,INTENT(IN) :: + & origin(3),Dx,Dy,Dz,eta + CCTK_REAL,INTENT(IN),DIMENSION(:,:,:) :: + & Psi,g00,gxx,gxy,gxz,gyy,gyz,gzz, + & hxx,hxy,hxz,hyy,hyz,hzz + CCTK_REAL,INTENT(IN),DIMENSION(:) :: + & x,y,z + CCTK_REAL,INTENT(INOUT),DIMENSION(:,:,:) :: + & Extract_temp3d + +c Output variables + CCTK_REAL,INTENT(OUT) :: + & ADMmass(2),mass,rsch,Qodd(:,2:,0:),Qeven(:,2:,0:),dtaudt, + & momentum(3),spin(3) + +c Local variables passed on + CCTK_REAL :: + & theta(Nt),phi(Np) + CCTK_REAL,DIMENSION(Nt,Np) :: + & Psis,g00s,gxxs,gxys,gxzs,gyys,gyzs,gzzs,dPsis,dgxxs,dgxys,dgxzs, + & dgyys,dgyzs,dgzzs,grr,grt,grp,gtt,gtp,gpp,dgtt,dgtp,dgpp, + & ADMmass_int1,ADMmass_int2, + & momentum_int1,momentum_int2,momentum_int3, + & spin_int1,spin_int2,spin_int3 + +c Local variables here only + INTEGER :: + & i,j + CCTK_REAL,PARAMETER :: + & half = 0.5D0, + & one = 1.0D0, + & two = 2.0D0 + CCTK_REAL :: + & Pi,Dt,Dp + +c ------------------------------------------------------------------ + + +c ------------------------------------------------------------------ +c +c 1. Specify polar coordinates on chosen 2-surface +c +c ------------------------------------------------------------------ + + Pi = two*ASIN(one) + +c Set polar gridspacing + + SELECT CASE (igrid) + CASE (0) ! full grid + Dt = Pi/DBLE(Nt) + Dp = two*Pi/DBLE(Np) + CASE (1) ! octant grid + Dt = half*Pi/DBLE(Nt) + Dp = half*Pi/DBLE(Np) + CASE (2) ! cartoon + Dt = Pi/DBLE(Nt) + Dp = 2D0*Pi + CASE (3) ! bitant + Dt = half*Pi/DBLE(Nt) + Dp = two*Pi/DBLE(Np) + CASE DEFAULT + WRITE (*,*) "Grid incorrectly set in D3_extract" + STOP + END SELECT + + +c Set coordinates + + DO i = 1, Nt + theta(i) = (DBLE(i) - half)* Dt + ENDDO + DO j = 1, Np + phi(j) = (DBLE(j) - half)* Dp + ENDDO + + IF (igrid == 2) phi(1) = 0D0 + + +c ------------------------------------------------------------------ +c +c 2. Project quantities onto the 2-surface +c +c ------------------------------------------------------------------ + + CALL D3_to_D2(cctkGH,do_ADMmass,do_momentum,do_spin, + & Psi_power,origin,myproc,Dx,Dy,Dz,Psi, + & g00,gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz, + & x,y,z,eta,Nt,Np,theta,phi,Psis,g00s,gxxs,gxys,gxzs, + & gyys,gyzs,gzzs,dPsis,dgxxs,dgxys,dgxzs,dgyys,dgyzs,dgzzs, + & ADMmass_int1,ADMmass_int2,momentum_int1,momentum_int2, + & momentum_int3,spin_int1,spin_int2,spin_int3,Extract_temp3d) + +c ------------------------------------------------------------------ +c Now this is done the grid on processor 0 has correct values and the +c grid on all other processors have junk. This means that processor +c 0 needs to calculate the wave forms and the other processors have +c to sit and do nothing. So we shall put this in a conditional +c to only happen on processor 0. +c PW Jan 31 98 + + if (myproc .eq. 0) then +c ------------------------------------------------------------------ +c +c 3. Transform tensor components to polar coordinates on 2-surface +c +c ------------------------------------------------------------------ + + CALL cartesian_to_spherical(theta,phi,eta,gxxs,gxys,gxzs,gyys, + & gyzs,gzzs,grr,grt,grp,gtt,gtp,gpp,dgxxs,dgxys,dgxzs,dgyys, + & dgyzs,dgzzs,dgtt,dgtp,dgpp) + + +c ------------------------------------------------------------------ +c +c 4. Calculate physical quantities on 2-surface +c +c ------------------------------------------------------------------ + + CALL unphysical_to_physical(grr,grt,grp,gtt,gtp,gpp,dgtt,dgtp, + & dgpp,Psis,dPsis,Psi_power) + + +c ------------------------------------------------------------------ +c +c 5. Extract gauge invariant quantities on 2-surface +c +c ------------------------------------------------------------------ + + CALL D2_extract(eta,igrid,Nt,Np,theta,phi,all_modes,l,m,g00s,grr, + & grt,grp,gtt,gtp,gpp,dgtt,dgtp,dgpp, + & do_ADMmass,ADMmass_int1,ADMmass_int2, + & do_momentum,momentum_int1,momentum_int2,momentum_int3, + & do_spin,spin_int1,spin_int2,spin_int3, + & mass,rsch,Qodd,Qeven,ADMmass,momentum,spin,dtaudt) + +c End of myproc = 0 conditional + endif + + END SUBROUTINE D3_extract + + + diff --git a/src/D3_extract_int.F b/src/D3_extract_int.F new file mode 100644 index 0000000..3bb617a --- /dev/null +++ b/src/D3_extract_int.F @@ -0,0 +1,44 @@ + +#include "cctk.h" + + MODULE D3_extract_int + +c ------------------------------------------------------------------ + + INTERFACE + + SUBROUTINE D3_extract(cctkGH,do_ADMmass,do_momentum,do_spin,igrid, + & origin,myproc, + & Nt,Np,all_modes,l,m,x,y,z,Dx,Dy,Dz,Psi_power,Psi, + & g00,gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz, + & eta,ADMmass,momentum,spin,mass,rsch,Qodd,Qeven, + & Extract_temp3d,dtaudt) + + IMPLICIT NONE + + CCTK_POINTER :: cctkGH + + INTEGER,INTENT(IN) :: + & igrid,l,m,Psi_power,myproc + CCTK_INT,INTENT(IN) :: + & Nt,Np,all_modes,do_momentum,do_spin + INTEGER,INTENT(IN) :: + & do_ADMmass(2) + CCTK_REAL,INTENT(IN) :: + & origin(3),Dx,Dy,Dz,eta + CCTK_REAL,INTENT(IN),DIMENSION(:) :: + & x,y,z + CCTK_REAL,INTENT(IN),DIMENSION(:,:,:) :: + & Psi,g00,gxx,gxy,gxz,gyy,gyz,gzz, + & hxx,hxy,hxz,hyy,hyz,hzz + CCTK_REAL,INTENT(INOUT),DIMENSION(:,:,:) :: + & Extract_Temp3d + CCTK_REAL,INTENT(OUT) :: + & ADMmass(2),mass,rsch,Qodd(:,:,:),Qeven(:,:,:),dtaudt, + & momentum(3),spin(3) + + END SUBROUTINE + + END INTERFACE + + END MODULE D3_extract_int diff --git a/src/D3_to_D2.F b/src/D3_to_D2.F new file mode 100644 index 0000000..884f8b8 --- /dev/null +++ b/src/D3_to_D2.F @@ -0,0 +1,398 @@ + +#include "cctk.h" + +c ======================================================================== + + SUBROUTINE D3_to_D2(cctkGH,do_ADMmass,do_momentum,do_spin, + & Psi_power,origin,myproc,Dx,Dy,Dz,Psi, + & g00,gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz, + & x,y,z,eta,Nt,Np,theta,phi, + & Psis,g00s,gxxs,gxys,gxzs,gyys,gyzs,gzzs,dPsis,dgxxs,dgxys,dgxzs, + & dgyys,dgyzs,dgzzs,ADMmass_int1,ADMmass_int2, + & momentum_int1,momentum_int2,momentum_int3, + & spin_int1,spin_int2,spin_int3,Extract_temp3d) + +c ------------------------------------------------------------------------ +c +c Project the 3-metric and its 1st radial derivatives onto a +c 2-sphere. +c +c ------------------------------------------------------------------------ + + USE met_rad_der_int + USE ADMmass_integrand3D_int + USE momentum_integrand3D_int + USE spin_integrand3D_int + + IMPLICIT NONE + +c Input variables + + CCTK_POINTER :: cctkGH + + INTEGER,INTENT(IN) :: + & myproc,Psi_power + CCTK_INT, INTENT(IN) :: + & Nt,Np,do_momentum,do_spin + CCTK_REAL,INTENT(IN) :: + & origin(3),Dx,Dy,Dz,eta + CCTK_REAL,INTENT(IN),DIMENSION(:) :: + & theta,phi,x,y,z + CCTK_REAL,INTENT(IN),DIMENSION(:,:,:) :: + & Psi,g00,gxx,gxy,gxz,gyy,gyz,gzz, + & hxx,hxy,hxz,hyy,hyz,hzz + CCTK_REAL,INTENT(INOUT),DIMENSION(:,:,:) :: + & Extract_temp3d + LOGICAL :: + & do_ADMmass(2) + +c Output variables + + CCTK_REAL,INTENT(OUT),DIMENSION(:,:) :: + & Psis,g00s,gxxs,gxys,gxzs,gyys, + & gyzs,gzzs,dPsis,dgxxs,dgxys,dgxzs,dgyys,dgyzs,dgzzs, + & ADMmass_int1,ADMmass_int2, + & momentum_int1,momentum_int2,momentum_int3, + & spin_int1,spin_int2,spin_int3 + + +c Local variables, passed on + + LOGICAL :: + & err_flag + INTEGER :: + & iorder,npoints,Nx,Ny,Nz + INTEGER,DIMENSION(Nt,Np) :: + & ib,jb,kb + CCTK_REAL,DIMENSION(Nt,Np) :: + & xs,ys,zs,ux,uy,uz,xb,yb,zb + +c Local variables, here only + + INTEGER :: + & i,j,handle,ierror + CCTK_REAL :: + & xmin,xmax,ymin,ymax,zmin,zmax + +c ------------------------------------------------------------------------ + +c Initial Stuff +c ------------- + Nx = SIZE(x) + Ny = SIZE(y) + Nz = SIZE(z) + + +c Compute Cartesian coordinates on the surface +c -------------------------------------------- + DO j = 1, Np + DO i = 1, Nt + + xs(i,j) = origin(1)+eta*SIN(theta(i))*COS(phi(j)) + ys(i,j) = origin(2)+eta*SIN(theta(i))*SIN(phi(j)) + zs(i,j) = origin(3)+eta*COS(theta(i)) + + ENDDO + ENDDO + + +c Only do interpolation on one processor +c -------------------------------------- + SELECT CASE (myproc) + + CASE (0) + npoints = Np*Nt + + CASE DEFAULT + npoints = 0 + + END SELECT + + +c Choose the interpolator +c ----------------------- + call CCTK_GetInterpHandle (handle, "simple") + +c Find the origin of the spatial coordinates +c ------------------------------------------ + call CCTK_CoordRange(ierror,cctkGH,xmin,xmax,"x") + call CCTK_CoordRange(ierror,cctkGH,ymin,ymax,"y") + call CCTK_CoordRange(ierror,cctkGH,zmin,zmax,"z") + + +c Project un-physical metric and conformal factor onto sphere +c ------------------------------------------------------------ + call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 8, 8, + $ nx, ny, nz, + $ xs, ys, zs, + $ CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, + $ CCTK_VARIABLE_REAL, + $ xmin, ymin, zmin, + $ dx, dy, dz, + $ Psi,g00,gxx,gxy,gxz,gyy,gyz,gzz, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ Psis,g00s,gxxs,gxys,gxzs,gyys,gyzs,gzzs, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL + $ ) + + +c Calculate radial derivatives and project onto sphere +c ---------------------------------------------------- + CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,Psi,Extract_temp3d) + CALL CCTK_SyncGroup(cctkGH,"extract::temps") + call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1, + $ nx, ny, nz, + $ xs, ys, zs, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ xmin, ymin, zmin, + $ dx, dy, dz, + $ Extract_temp3d, + $ CCTK_VARIABLE_REAL, + $ dPsis, + $ CCTK_VARIABLE_REAL + $ ) + + CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gxx,Extract_temp3d) + CALL CCTK_SyncGroup(cctkGH,"extract::temps") + call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1, + $ nx, ny, nz, + $ xs, ys, zs, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ xmin, ymin, zmin, + $ dx, dy, dz, + $ Extract_temp3d, + $ CCTK_VARIABLE_REAL, + $ dgxxs, + $ CCTK_VARIABLE_REAL + $ ) + + CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gxy,Extract_temp3d) + CALL CCTK_SyncGroup(cctkGH,"extract::temps") + call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1, + $ nx, ny, nz, + $ xs, ys, zs, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ xmin, ymin, zmin, + $ dx, dy, dz, + $ Extract_temp3d, + $ CCTK_VARIABLE_REAL, + $ dgxys, + $ CCTK_VARIABLE_REAL + $ ) + + CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gxz,Extract_temp3d) + CALL CCTK_SyncGroup(cctkGH,"extract::temps") + call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1, + $ nx, ny, nz, + $ xs, ys, zs, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ xmin, ymin, zmin, + $ dx, dy, dz, + $ Extract_temp3d, + $ CCTK_VARIABLE_REAL, + $ dgxzs, + $ CCTK_VARIABLE_REAL + $ ) + + CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gyy,Extract_temp3d) + CALL CCTK_SyncGroup(cctkGH,"extract::temps") + call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1, + $ nx, ny, nz, + $ xs, ys, zs, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ xmin, ymin, zmin, + $ dx, dy, dz, + $ Extract_temp3d, + $ CCTK_VARIABLE_REAL, + $ dgyys, + $ CCTK_VARIABLE_REAL + $ ) + + CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gyz,Extract_temp3d) + CALL CCTK_SyncGroup(cctkGH,"extract::temps") + call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1, + $ nx, ny, nz, + $ xs, ys, zs, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ xmin, ymin, zmin, + $ dx, dy, dz, + $ Extract_temp3d, + $ CCTK_VARIABLE_REAL, + $ dgyzs, + $ CCTK_VARIABLE_REAL + $ ) + + + CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gzz,Extract_temp3d) + CALL CCTK_SyncGroup(cctkGH,"extract::temps") + call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1, + $ nx, ny, nz, + $ xs, ys, zs, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ xmin, ymin, zmin, + $ dx, dy, dz, + $ Extract_temp3d, + $ CCTK_VARIABLE_REAL, + $ dgzzs, + $ CCTK_VARIABLE_REAL + $ ) + + +c Calculate integrands for ADM masses +c ----------------------------------- +c Standard equation + IF (do_ADMmass(1)) THEN + CALL ADMmass_integrand3D(origin,Dx,Dy,Dz,x,y,z,gxx,gxy, + & gxz,gyy,gyz,gzz,Extract_temp3d,Psi,Psi_power) + CALL CCTK_SyncGroup(cctkGH,"extract::temps") + call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1, + $ nx, ny, nz, + $ xs, ys, zs, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ xmin, ymin, zmin, + $ dx, dy, dz, + $ Extract_temp3d, + $ CCTK_VARIABLE_REAL, + $ ADMmass_int1, + $ CCTK_VARIABLE_REAL + $ ) + + END IF + +c Conformal equation + IF (do_ADMmass(2)) THEN + ADMmass_int2 = -eta**2/2D0/3.1416D0*dPsis + ENDIF + +c Calculate integrands for momentum +c --------------------------------- + IF (do_momentum==1) THEN + + CALL momentum_integrand3D(origin,Dx,Dy,Dz,x,y,z, + & gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz, + & Extract_temp3d,Psi,Psi_power) + CALL CCTK_SyncGroup(cctkGH,"extract::temps") + call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1, + $ nx, ny, nz, + $ xs, ys, zs, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ xmin, ymin, zmin, + $ dx, dy, dz, + $ Extract_temp3d, + $ CCTK_VARIABLE_REAL, + $ momentum_int1, + $ CCTK_VARIABLE_REAL + $ ) + + CALL momentum_integrand3D(origin,Dx,Dy,Dz,x,y,z, + & gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz, + & Extract_temp3d,Psi,Psi_power) + CALL CCTK_SyncGroup(cctkGH,"extract::temps") + call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1, + $ nx, ny, nz, + $ xs, ys, zs, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ xmin, ymin, zmin, + $ dx, dy, dz, + $ Extract_temp3d, + $ CCTK_VARIABLE_REAL, + $ momentum_int2, + $ CCTK_VARIABLE_REAL + $ ) + + + CALL momentum_integrand3D(origin,Dx,Dy,Dz,x,y,z, + & gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz, + & Extract_temp3d,Psi,Psi_power) + CALL CCTK_SyncGroup(cctkGH,"extract::temps") + call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1, + $ nx, ny, nz, + $ xs, ys, zs, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ xmin, ymin, zmin, + $ dx, dy, dz, + $ Extract_temp3d, + $ CCTK_VARIABLE_REAL, + $ momentum_int3, + $ CCTK_VARIABLE_REAL + $ ) + + END IF + +c Calculate integrands for spin +c ----------------------------- + IF (do_momentum==1) THEN + + CALL spin_integrand3D(origin,Dx,Dy,Dz,x,y,z, + & gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz, + & Extract_temp3d,Psi,Psi_power) + CALL CCTK_SyncGroup(cctkGH,"extract::temps") + call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1, + $ nx, ny, nz, + $ xs, ys, zs, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ xmin, ymin, zmin, + $ dx, dy, dz, + $ Extract_temp3d, + $ CCTK_VARIABLE_REAL, + $ spin_int1, + $ CCTK_VARIABLE_REAL + $ ) + + CALL spin_integrand3D(origin,Dx,Dy,Dz,x,y,z, + & gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz, + & Extract_temp3d,Psi,Psi_power) + CALL CCTK_SyncGroup(cctkGH,"extract::temps") + call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1, + $ nx, ny, nz, + $ xs, ys, zs, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ xmin, ymin, zmin, + $ dx, dy, dz, + $ Extract_temp3d, + $ CCTK_VARIABLE_REAL, + $ spin_int2, + $ CCTK_VARIABLE_REAL + $ ) + + + CALL spin_integrand3D(origin,Dx,Dy,Dz,x,y,z, + & gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz, + & Extract_temp3d,Psi,Psi_power) + CALL CCTK_SyncGroup(cctkGH,"extract::temps") + call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1, + $ nx, ny, nz, + $ xs, ys, zs, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ xmin, ymin, zmin, + $ dx, dy, dz, + $ Extract_temp3d, + $ CCTK_VARIABLE_REAL, + $ spin_int3, + $ CCTK_VARIABLE_REAL + $ ) + + END IF + + END SUBROUTINE D3_to_D2 + + + + + + + + + + + + + + + + + diff --git a/src/D3_to_D2_int.F b/src/D3_to_D2_int.F new file mode 100644 index 0000000..297e22e --- /dev/null +++ b/src/D3_to_D2_int.F @@ -0,0 +1,45 @@ + +#include "cctk.h" + + MODULE D3_to_D2_int + +c ------------------------------------------------------------------ + + INTERFACE + + SUBROUTINE D3_to_D2(cctkGH,do_ADMmass,do_momentum,do_spin, + & Psi_power,origin,myproc,Dx,Dy,Dz,Psi, + & g00,gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz, + & x,y,z,eta,Nt,Np,theta,phi,Psis,g00s,gxxs,gxys, + & gxzs,gyys,gyzs,gzzs,dPsis,dgxxs,dgxys,dgxzs,dgyys,dgyzs, + & dgzzs,ADMmass_int1,ADMmass_int2, + & momentum_int1,momentum_int2,momentum_int3, + & spin_int1,spin_int2,spin_int3,Extract_temp3d) + + IMPLICIT NONE + CCTK_POINTER :: cctkGH + INTEGER,INTENT(IN) :: + & myproc,Psi_power + CCTK_INT, INTENT(IN) :: + & Nt,Np,do_momentum,do_spin + INTEGER :: + & do_ADMmass(2) + CCTK_REAL,INTENT(IN) :: + & origin(3),Dx,Dy,Dz,eta + CCTK_REAL,INTENT(IN),DIMENSION(:,:,:) :: + & Psi,g00,gxx,gxy,gxz,gyy,gyz,gzz, + & hxx,hxy,hxz,hyy,hyz,hzz + CCTK_REAL,INTENT(IN),DIMENSION(:) :: + & x,y,z,theta,phi + CCTK_REAL,INTENT(OUT),DIMENSION(:,:) :: + & Psis,g00s,gxxs,gxys,gxzs,gyys,gyzs,gzzs,dPsis,dgxxs,dgxys,dgxzs, + & dgyys,dgyzs,dgzzs,ADMmass_int1,ADMmass_int2, + & momentum_int1,momentum_int2,momentum_int3, + & spin_int1,spin_int2,spin_int3 + CCTK_REAL,INTENT(INOUT),DIMENSION(:,:,:) :: + & Extract_temp3d + END SUBROUTINE + + END INTERFACE + + END MODULE D3_to_D2_int diff --git a/src/Extract.F b/src/Extract.F new file mode 100644 index 0000000..1aa8534 --- /dev/null +++ b/src/Extract.F @@ -0,0 +1,814 @@ +c/*@@ +c @file Extract.F +c @date 28th October 1997 +c @author Gab Allen +c @desc Waveform extraction +c +c @enddesc +c@@*/ + +#include "cctk.h" +#include "cctk_parameters.h" +#include "cctk_arguments.h" + +c/*@@ +c @routine Extract +c @date 28th October 1997 +c @author Gab Allen +c @desc Entry point for waveform extraction routines +c +c @enddesc +c @calls D3_extract +c @calledby No idea +c @history +c +c @endhistory +c@@*/ + + + SUBROUTINE Extract(CCTK_FARGUMENTS) + + USE D3_extract_int + + IMPLICIT NONE + + DECLARE_CCTK_FARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + +c Non-Cactus input variables for D3_extract + + INTEGER :: + & igrid,lmode,mmode,Psi_power=1 + INTEGER :: + & do_ADMmass(2) + CCTK_REAL :: + & orig(3),radius + CCTK_REAL :: + & x_1d(cctk_lsh(1)),y_1d(cctk_lsh(2)),z_1d(cctk_lsh(3)) + + +c Output variables from D3_extract + + CCTK_REAL :: + & dtaudt,ADMmass(2),mass,rsch,momentum(3),spin(3) + CCTK_REAL,ALLOCATABLE :: + & Qodd(:,:,:), Qeven(:,:,:) + + +c Local variables + + INTEGER :: + & ix,iy,iz,fn1,fn2,lmin,lmax,mmin,mmax,lstep,mstep, + & il,im,it,ndet,idet,nx_parfile,nchar,iconv,ioutput + INTEGER,SAVE :: openfiles = 1 + INTEGER,PARAMETER :: + & max_detectors = 9,number_timecoords = 2 + CCTK_REAL,PARAMETER :: + & two = 2.0D0 + CCTK_REAL :: + & time,r1,r2,dr,r_max,xmin,xmax,ymin,ymax,zmin,zmax,gf_min,gf_max + CCTK_REAL :: + & detector(max_detectors) + CHARACTER*100 :: + & out1,out2,massfile,rschfile,ADMmassfile1,ADMmassfile2, + & momentumfile1,momentumfile2,momentumfile3,spinfile1,spinfile2,spinfile3 + CHARACTER*3 :: + & timestring + + LOGICAL :: use_proper,use_coordinate,do_output + LOGICAL,SAVE :: firsttime + CCTK_REAL,DIMENSION(10),SAVE :: proper_time + CCTK_REAL,SAVE :: last_time + INTEGER, save :: open_file_level(10) = 0 + + INTEGER myproc,ierr + CCTK_REAL dx,dy,dz + character*200 filestr + character*80 infoline + +c ------------------------------------------------------------------ + +c Get the output directory sorted out + call CCTK_FortranString(nchar,outdir,filestr) + if (firsttime) then + call CCTK_mkdir(ierr,filestr); + firsttime = .FALSE. + end if + +c ------------------------------------------------------------------ + + myproc = CCTK_MyProc(cctkGH) + dx = cctk_delta_space(1) + dy = cctk_delta_space(2) + dz = cctk_delta_space(3) + +c ------------------------------------------------------------------ +c +c 0. Check to see if Extract should have been called +c +c ------------------------------------------------------------------ + + IF (MODULO(int(cctk_iteration),int(itout)) .NE. 0) RETURN + + if (verbose == 1) then + write(infoline,'(A24,G12.7)') + & 'Calling Extract at time ',cctk_time + call CCTK_INFO(infoline) + end if + +c ------------------------------------------------------------------ +c +c 1. Initial Stuff +c +c ------------------------------------------------------------------ + + +c Get the number of polar divisions + + IF (MOD(int(Nt),2) == 1 .OR. MOD(int(Np),2) == 1 .OR. + & Nt < 0 .OR. Np < 0) THEN + call CCTK_WARN(0,"Error in Nt or Np in Extract") + END IF + + +c Get the origin of spherical symmetry + + orig(1) = origin_x + orig(2) = origin_y + orig(3) = origin_z + + +c Set the value of igrid + + IF (CCTK_EQUALS(domain,"octant")) THEN + igrid = 1 + ELSEIF (CCTK_EQUALS(domain,"bitant")) THEN + igrid = 3 + ELSE + igrid = 0 + END IF + +#ifdef THORN_CARTOON_2D + IF (contains("cartoon_active","yes") .NE. 0) THEN + igrid = 2 + Np = 1 + END IF +#endif + +c Create 1D coordinate arrays + + do ix = 1, cctk_lsh(1) + x_1d(ix) = x(ix,1,1) + enddo + do iy = 1, cctk_lsh(2) + y_1d(iy) = y(1,iy,1) + enddo + do iz = 1, cctk_lsh(3) + z_1d(iz) = z(1,1,iz) + enddo + +c See if the ADM mass should be calculated and output + IF (doADMmass == 1) THEN + do_ADMmass(1) = 1 + IF (use_conformal == 1) THEN + do_ADMmass(2) = 1 + ELSE + do_ADMmass(2) = 0 + ENDIF + ELSE + do_ADMmass(:) = 0 + END IF + +c Set array containing gtt ... for the moment I am lazy and +c just make it the lapse and do not include the shift! + g00 = alp + +c What kind of timecoordinate do I want to use + if (CCTK_EQUALS(timecoord,"both")) then + use_proper = .true. + use_coordinate = .true. + if (cctk_iteration == 0) then + proper_time = 0 + last_time = cctk_time + end if + elseif (CCTK_EQUALS(timecoord,"proper")) then + use_proper = .true. + use_coordinate = .false. + if (cctk_iteration == 0) then + proper_time = 0 + last_time = cctk_time + end if + elseif (CCTK_EQUALS(timecoord,"coord")) then + use_proper = .false. + use_coordinate = .true. + else + use_proper = .false. + use_coordinate = .true. + endif + +c ------------------------------------------------------------------ +c +c 2. Sort out which modes to extract +c +c ------------------------------------------------------------------ + + lmode = l_mode + + IF (all_modes == 1) THEN + lmin = 2 + lmax = lmode +c Get m_mode in this case too, or the T3E has a junk value (PW) + mmode = 0 + ALLOCATE(Qodd(2,2:lmode,0:lmode),Qeven(2,2:lmode,0:lmode)) + ELSE + lmin = lmode + lmax = lmode + mmode = m_mode + ALLOCATE(Qodd(2,2:lmode,0:mmode),Qeven(2,2:lmode,0:mmode)) + ENDIF + + +c Check do not have bad values for the modes + + IF (lmode < 2 .OR. mmode > lmode .OR. mmode < 0) THEN + call CCTK_WARN(0,"Error in lmode,mmode in Extract") + END IF + IF (lmode > 9) THEN + call CCTK_WARN(4,"Filenames will probably be crazy in Extract") + END IF + +c If using an octant, do not do odd modes + + IF (igrid == 1 .OR. igrid == 2) THEN + lstep = 2 ; mstep = 2 + ELSE + lstep = 1 ; mstep = 1 + END IF + +c If we have cartoon, then jump past all the m-modes in loop + IF (igrid == 2) mstep = lmode+1 + + +c ------------------------------------------------------------------ +c +c 3. Find maximum radius of extraction on grid +c +c ------------------------------------------------------------------ + + call CCTK_CoordRange(ierr,cctkGH,xmin,xmax,"x") + call CCTK_CoordRange(ierr,cctkGH,ymin,ymax,"y") + call CCTK_CoordRange(ierr,cctkGH,zmin,zmax,"z") + + IF (igrid == 2) THEN + r_max = MIN(xmax-two*Dx-orig(1), + & zmax-two*Dz-orig(3)) + ELSE IF (igrid == 1) THEN + r_max = MIN(zmax-two*Dx-orig(1), + & ymax-two*Dy-orig(2), + & zmax-two*Dz-orig(3)) + ELSE IF (igrid == 3) THEN + r_max = MIN(xmax-two*Dx-orig(1), + & ymax-two*Dy-orig(2), + & zmax-two*Dz-orig(3), + & DABS(xmin+two*Dx-orig(1)), + & DABS(ymin+two*Dy-orig(2))) + ELSE + r_max = MIN(xmax-two*Dx-orig(1), + & ymax-two*Dy-orig(2), + & zmax-two*Dz-orig(3), + & DABS(xmin+two*Dx-orig(1)), + & DABS(ymin+two*Dy-orig(2)), + & DABS(zmin+two*Dz-orig(3))) + END IF + + + + +c ------------------------------------------------------------------ +c +c 4. Extract Cauchy initial data for a linear wave equation +c +c ------------------------------------------------------------------ + + test_Cauchy: IF (Cauchy == 1) THEN + + it = Cauchy_timestep + + +c check_timestep: IF (cctk_iteration == it) THEN + + check_timestep: IF (open_file_level(1) .eq. 0) THEN + +c Open output files + write (*,*) "Open output file on level ", 1 + open_file_level(1)=1; + + test_myproc1: IF (CCTK_MyProc(cctkGH) == 0 ) THEN + + out1 = filestr(1:nchar)//"/rsch_ini.rl" + OPEN(UNIT= 77,file = out1,STATUS="unknown") + out1 = filestr(1:nchar)//"/mass_ini.rl" + OPEN(UNIT= 78,file = out1,STATUS="unknown") + IF (do_ADMmass(1) == 1) THEN + out1 = filestr(1:nchar)//"/ADMmass_ini.rl" + OPEN(UNIT= 79,file = out1,STATUS="unknown") + ENDIF + IF (do_ADMmass(2) == 1) THEN + out1 = filestr(1:nchar)//"/ADMmassc_ini.rl" + OPEN(UNIT= 80,file = out1,STATUS="unknown") + ENDIF + + IF (do_momentum == 1) THEN + out1 = filestr(1:nchar)//"/momentum_x_ini.rl" + OPEN(UNIT=781,file = out1,STATUS="unknown") + out1 = filestr(1:nchar)//"/momentum_y_ini.rl" + OPEN(UNIT=782,file = out1,STATUS="unknown") + out1 = filestr(1:nchar)//"/momentum_z_ini.rl" + OPEN(UNIT=783,file = out1,STATUS="unknown") + ENDIF + IF (do_spin == 1) THEN + out1 = filestr(1:nchar)//"/spin_x_ini.rl" + OPEN(UNIT=784,file = out1,STATUS="unknown") + out1 = filestr(1:nchar)//"/spin_y_ini.rl" + OPEN(UNIT=785,file = out1,STATUS="unknown") + out1 = filestr(1:nchar)//"/spin_z_ini.rl" + OPEN(UNIT=786,file = out1,STATUS="unknown") + ENDIF + + loop_l1: DO il = lmin,lmax,lstep + + IF (all_modes == 0) THEN + mmin = mmode ; mmax = mmode + ELSE + mmin = 0 ; mmax = il + END IF + + loop_m1: DO im = mmin,mmax,lstep + + + fn1 = il*10+im + out1 = filestr(1:nchar)//"/Qeven_ini_"// + & CHAR(il+48)//CHAR(im+48)//".rl" + OPEN(UNIT=fn1,FILE=out1,STATUS="unknown") + +c Only print odd-parity if full grid + IF (igrid == 0) THEN + fn2 = 100+il*10+im + out2 = filestr(1:nchar)//"/Qodd_ini_"// + & CHAR(il+48)//CHAR(im+48)//".rl" + OPEN(UNIT=fn2,FILE=out2,STATUS="unknown") + END IF + + END DO loop_m1 + + END DO loop_l1 + + END IF test_myproc1 + + +c Find range of extraction radii + + r1 = Cauchy_r1 + r2 = r_max + dr = Cauchy_dr + + + +c Do extraction at each radius + + IF (verbose == 1) THEN + WRITE(*,*) + WRITE(*,*) "Extracting Cauchy initial data" + WRITE(*,*) "------------------------------" + WRITE(*,*) "From r =",r1 + WRITE(*,*) "To r =",r2 + WRITE(*,*) "In steps of ",dr + WRITE(*,*) + END IF + + radius = r1 + + extract_at_each_radius: DO WHILE (radius < r2) + + CALL D3_extract(cctkGH,do_ADMmass,do_momentum,do_spin, + & igrid,orig,myproc,Nt,Np,all_modes,lmode, + & mmode,x_1d,y_1d,z_1d,Dx,Dy,Dz,Psi_power,Psi,g00, + & gxx,gxy,gxz,gyy,gyz,gzz,kxx,kxy,kxz,kyy,kyz,kzz, + & radius,ADMmass,momentum,spin,mass,rsch, + & Qodd,Qeven,temp3d,dtaudt) + + IF (verbose == 1) THEN + WRITE(*,*) "Extracted at r =",radius + WRITE(*,*) "Schwarzschild radius =",rsch + WRITE(*,*) "Schwarzschild mass =",mass + if (do_ADMmass(1) == 1) then + WRITE(*,*) "ADM mass =",ADMmass(1) + end if + if (do_ADMmass(2) == 1) then + WRITE(*,*) "ADM mass (conformal) =",ADMmass(2) + end if + if (do_momentum == 1) then + WRITE(*,*) "momentum = (",momentum(1),",", + & momentum(2),",",momentum(3),")" + end if + if (do_spin == 1) then + WRITE(*,*) "spin = (",spin(1),",", + & spin(2),",",spin(3),")" + end if + ENDIF + + +c Write to file at each radius + + test_myproc2: IF (CCTK_MyProc(cctkGH) == 0) THEN + + WRITE( 77,*) radius,rsch + WRITE( 78,*) radius,mass + IF (do_ADMmass(1) == 1) WRITE( 79,*) radius,ADMmass(1) + IF (do_ADMmass(2) == 1) WRITE( 80,*) radius,ADMmass(2) + IF (do_momentum == 1) then + WRITE(781,*) radius,momentum(1) + WRITE(782,*) radius,momentum(2) + WRITE(783,*) radius,momentum(3) + end if + + IF (do_spin == 1) THEN + WRITE(784,*) radius,spin(1) + WRITE(785,*) radius,spin(2) + WRITE(786,*) radius,spin(3) + END IF + + loop_l2: DO il = lmin,lmax,lstep + + IF (all_modes == 0) THEN + mmin = mmode ; mmax = mmode + ELSE + mmin = 0 ; mmax = il + END IF + + loop_m2: DO im = mmin,mmax,mstep + + fn1 = il*10+im + WRITE(fn1,*) rsch,Qeven(:,il,im) + +c Only print odd-parity if full grid + IF (igrid == 0) THEN + fn2 = 100+il*10+im + WRITE(fn2,*) rsch,Qeven(:,il,im) + END IF + + END DO loop_m2 + + END DO loop_l2 + + END IF test_myproc2 + + radius = radius+dr + + END DO extract_at_each_radius + + +c Now close all the files + + test_myproc3: IF (myproc == 0) THEN + + CLOSE( 77) ! Schwarzschild radius + CLOSE( 78) ! Schwarzschild mass + IF (do_ADMmass(1) == 1) CLOSE( 79) ! ADM mass + IF (do_ADMmass(2) == 1) CLOSE( 80) ! ADM mass + + IF (do_momentum == 1) THEN + CLOSE(781) ! momentum + CLOSE(782) ! momentum + CLOSE(783) ! momentum + end if + + IF (do_spin == 1) THEN + CLOSE(784) ! spin + CLOSE(785) ! spin + CLOSE(786) ! spin + END IF + + loop_l3: DO il = lmin,lmax,lstep + + IF (all_modes == 0) THEN + mmin = mmode ; mmax = mmode + ELSE + mmin = 0 ; mmax = il + END IF + + loop_m3: DO im = mmin,mmax,mstep + + fn1 = il*10+im + CLOSE(fn1) ! Qeven_ini_lm.dat + + IF (igrid == 0) THEN + fn2 = 100+il*10+im + CLOSE(fn2) ! Qodd_ini_lm.dat + END IF + + END DO loop_m3 + + END DO loop_l3 + + END IF test_myproc3 + + END IF check_timestep + + END IF test_Cauchy + + +c ------------------------------------------------------------------ +c +c 5. Extract waveforms at detectors +c +c ------------------------------------------------------------------ + +c Cannot use the conformal equation for ADM mass now + do_ADMmass(2) = 0 + + ndet = num_detectors + + test_detectors: IF (ndet > 0) THEN + + IF (ndet > 9) THEN + call CCTK_WARN(0,"Too many detectors in Extract") + END IF + + detector(1) = detector1 + detector(2) = detector2 + detector(3) = detector3 + detector(4) = detector4 + detector(5) = detector5 + detector(6) = detector6 + detector(7) = detector7 + detector(8) = detector8 + detector(9) = detector9 + + DO idet = ndet+1, max_detectors + detector(idet) = 0.0D0 + END DO + + DO idet = 1, ndet + IF (detector(idet) > r_max) THEN + IF (openfiles == 1 .OR. cctk_iteration == it) THEN + call CCTK_WARN(8,"Removing detectors outside coordinate range") +c IF (iconv == 0) STOP + END IF + ndet =idet-1 + GOTO 404 + ENDIF + END DO + 404 CONTINUE + + IF (verbose == 1) THEN + IF (openfiles == 1) THEN + WRITE(*,*) + WRITE(*,*) "Extracting at detectors" + WRITE(*,*) "-----------------------" + WRITE(*,*) "Using ",ndet," detectors at" + DO idet=1,ndet + WRITE(*,*) "r = ",detector(idet) + ENDDO + END IF + END IF + + detector_loop: DO idet = 1,ndet + + radius = detector(idet) + + IF (verbose == 1) THEN + WRITE(*,*) "Calling extract at radius ... ",radius + END IF + + CALL D3_extract(cctkGH,do_ADMmass,do_momentum,do_spin, + & igrid,orig,myproc,Nt, + & Np,all_modes,lmode,mmode,x_1d,y_1d,z_1d, + & Dx,Dy,Dz,Psi_power,Psi,g00, + & gxx,gxy,gxz,gyy,gyz,gzz,kxx,kxy,kxz,kyy,kyz,kzz, + & radius,ADMmass,momentum,spin,mass,rsch, + & Qodd,Qeven,temp3d,dtaudt) + + test_verbose: + & IF (verbose == 1) THEN + + WRITE(*,*) + WRITE(*,*) "Detector ",idet," ..." + WRITE(*,*) "Schwarzschild radius =",rsch + WRITE(*,*) "Schwarzschild mass =",mass + + loop_l4: DO il = lmin,lmax,lstep + + IF (all_modes == 0) THEN + mmin = mmode ; mmax = mmode + ELSE + mmin = 0 ; mmax = il + END IF + + loop_m4: DO im = mmin,mmax,mstep + + WRITE(*,*) filestr(1:nchar)//"/Qeven_"// + & CHAR(il+48)//CHAR(im+48),Qeven(:,il,im) + WRITE(*,*) filestr(1:nchar)//"/Qodd_"// + & CHAR(il+48)//CHAR(im+48),Qodd(:,il,im) + + END DO loop_m4 + + END DO loop_l4 + + END IF test_verbose + +c Output to files + + test_myproc4: IF (myproc == 0) THEN + + DO ioutput=1,number_timecoords + + do_output = .false. + timestring = "ERR" + time = 0D0 + IF (ioutput == 1 .AND. use_coordinate) THEN + timestring = ".tl" + do_output = .true. + time = cctk_time + ELSEIF (ioutput == 2 .AND. use_proper) THEN + timestring = ".ul" + do_output = .true. + time = proper_time(idet) + (cctk_time-last_time)*dtaudt + proper_time(idet) = time + if (idet == ndet) last_time = cctk_time + ENDIF + + IF (do_output) THEN + +c Output extracted radius and mass + rschfile = filestr(1:nchar)//"/rsch_"//"R"// + & CHAR(idet+48)//timestring + massfile = filestr(1:nchar)//"/mass_"//"R"// + & CHAR(idet+48)//timestring + ADMmassfile1 = filestr(1:nchar)//"/ADMmass_"//"R"// + & CHAR(idet+48)//timestring + ADMmassfile2 = filestr(1:nchar)//"/ADMmassc_"//"R"// + & CHAR(idet+48)//timestring + momentumfile1 = filestr(1:nchar)//"/momentum_x_"//"R"// + & CHAR(idet+48)//timestring + momentumfile2 = filestr(1:nchar)//"/momentum_y_"//"R"// + & CHAR(idet+48)//timestring + momentumfile3 = filestr(1:nchar)//"/momentum_z_"//"R"// + & CHAR(idet+48)//timestring + spinfile1 = filestr(1:nchar)//"/spin_x_"//"R"// + & CHAR(idet+48)//timestring + spinfile2 = filestr(1:nchar)//"/spin_y_"//"R"// + & CHAR(idet+48)//timestring + spinfile3 = filestr(1:nchar)//"/spin_z_"//"R"// + & CHAR(idet+48)//timestring + + IF (openfiles == 1) THEN + OPEN(UNIT= 77,FILE=rschfile,STATUS="unknown") + OPEN(UNIT= 78,FILE=massfile,STATUS="unknown") + IF (do_ADMmass(1) == 1) + & OPEN(UNIT= 79,FILE=ADMmassfile1,STATUS="unknown") + IF (do_ADMmass(2) == 1) + & OPEN(UNIT= 80,FILE=ADMmassfile2,STATUS="unknown") + IF (do_momentum == 1) THEN + OPEN(UNIT=781,FILE=momentumfile1,STATUS="unknown") + OPEN(UNIT=782,FILE=momentumfile2,STATUS="unknown") + OPEN(UNIT=783,FILE=momentumfile3,STATUS="unknown") + END IF + + IF (do_spin == 1) THEN + OPEN(UNIT=784,FILE=spinfile1,STATUS="unknown") + OPEN(UNIT=785,FILE=spinfile2,STATUS="unknown") + OPEN(UNIT=786,FILE=spinfile3,STATUS="unknown") + END IF + ELSE + OPEN(UNIT= 77,FILE=rschfile,STATUS="old", + & POSITION="append") + OPEN(UNIT= 78,FILE=massfile,STATUS="old", + & POSITION="append") + IF (do_ADMmass(1) == 1) + & OPEN(UNIT= 79,FILE=ADMmassfile1,STATUS="old", + & POSITION="append") + IF (do_ADMmass(2) == 1) + & OPEN(UNIT= 80,FILE=ADMmassfile2,STATUS="old", + & POSITION="append") + IF (do_momentum == 1) THEN + OPEN(UNIT=781,FILE=momentumfile1,STATUS="old", + & POSITION="append") + OPEN(UNIT=782,FILE=momentumfile2,STATUS="old", + & POSITION="append") + OPEN(UNIT=783,FILE=momentumfile3,STATUS="old", + & POSITION="append") + END IF + + IF (do_spin == 1)THEN + OPEN(UNIT=784,FILE=spinfile1,STATUS="old", + & POSITION="append") + OPEN(UNIT=785,FILE=spinfile2,STATUS="old", + & POSITION="append") + OPEN(UNIT=786,FILE=spinfile3,STATUS="old", + & POSITION="append") + END IF + + ENDIF + WRITE( 77,21) time,rsch + WRITE( 78,21) time,mass + IF (do_ADMmass(1) == 1) WRITE( 79,21) time,ADMmass(1) + IF (do_ADMmass(2) == 1) WRITE( 80,21) time,ADMmass(2) + IF (do_momentum == 1) THEN + WRITE(781,*) cctk_time,momentum(1) + WRITE(782,*) cctk_time,momentum(2) + WRITE(783,*) cctk_time,momentum(3) + END IF + IF (do_spin == 1) THEN + WRITE(784,*) cctk_time,spin(1) + WRITE(785,*) cctk_time,spin(2) + WRITE(786,*) cctk_time,spin(3) + END IF + + CLOSE( 77) + CLOSE( 78) + IF (do_ADMmass(1) == 1) CLOSE( 79) + IF (do_ADMmass(2) == 1) CLOSE( 80) + + IF (do_momentum == 1) THEN + CLOSE(781) + CLOSE(782) + CLOSE(783) + END IF + IF (do_spin == 1) THEN + CLOSE(782) + CLOSE(783) + CLOSE(784) + END IF + + +c Output gauge invariant variables + + loop_l5: DO il = lmin,lmax,lstep + + IF (all_modes == 0) THEN + mmin = mmode ; mmax = mmode + ELSE + mmin = 0 ; mmax = il + END IF + + loop_m5: DO im = mmin,mmax,mstep + + fn1 = 11 + fn2 = 12 + + out1 = filestr(1:nchar)//"/Qeven_"//"R" + & //CHAR(idet+48)//"_"//CHAR(il+48)// + & CHAR(im+48)//timestring + out2 = filestr(1:nchar)//"/Qodd_"//"R" + & //CHAR(idet+48)//"_"//CHAR(il+48)// + & CHAR(im+48)//timestring + +c Write even parity waveforms + IF (openfiles == 1) THEN + OPEN(UNIT=fn1,FILE=out1,STATUS="unknown") + ELSE + OPEN(UNIT=fn1,FILE=out1,STATUS="old", + & POSITION="append") + ENDIF + + WRITE(fn1,20) time,Qeven(:,il,im) + + CLOSE(fn1) + +c Only write odd parity waveforms if full grid + IF (igrid == 0) THEN + IF (openfiles == 1) THEN + OPEN(UNIT=fn2,FILE=out2,STATUS="unknown") + ELSE + OPEN(UNIT=fn2,FILE=out2,STATUS="old", + & POSITION="append") + ENDIF + + WRITE(fn2,20) time,Qodd(:,il,im) + + CLOSE(fn2) + END IF + + END DO loop_m5 + + END DO loop_l5 + + ENDIF + END DO + END IF test_myproc4 + + END DO detector_loop + + END IF test_detectors + + DEALLOCATE(Qodd,Qeven) + + openfiles = 0 + + 20 format(e24.15,e24.15,e24.15) + 21 format(e24.15,e24.15) + + END SUBROUTINE Extract diff --git a/src/cartesian_to_spherical.F b/src/cartesian_to_spherical.F new file mode 100644 index 0000000..d55f84f --- /dev/null +++ b/src/cartesian_to_spherical.F @@ -0,0 +1,116 @@ + +#include "cctk.h" + +c ================================================================== + + SUBROUTINE cartesian_to_spherical(theta,phi,r,gxxs,gxys,gxzs, + & gyys,gyzs, gzzs,grrs,grts,grps,gtts,gtps,gpps,dgxxs,dgxys, + & dgxzs,dgyys,dgyzs,dgzzs,dgtts,dgtps,dgpps) + +c ------------------------------------------------------------------ +c +c Convert components of a symmetric 2nd rank 3D tensor on a 2-sphere +c from Cartesian to spherical polar coordinates +c +c (x,y,z) ----> (r,t,p) +c +c ( gxx gxy gxz ) ( grr grt grp ) +c ( . gyy gyz ) ----> ( . gpp gtp ) +c ( . . gzz ) ( . . gpp ) +c +c Also convert radial derivatives of the tensor components +c +c ------------------------------------------------------------------ + + IMPLICIT NONE + +c Input variables + CCTK_REAL,INTENT(IN) :: + & r + CCTK_REAL,INTENT(IN),DIMENSION (:) :: + & theta,phi + CCTK_REAL,INTENT(IN),DIMENSION (:,:) :: + & gxxs,gxys,gxzs,gyys,gyzs,gzzs,dgxxs,dgxys,dgxzs,dgyys,dgyzs, + & dgzzs + +c Output variables + CCTK_REAL,INTENT(OUT),DIMENSION (:,:) :: + & grrs,grts,grps,gtts,gtps, + & gpps,dgtts,dgtps,dgpps + +c Local variables + INTEGER :: + & it,ip,Nt,Np + CCTK_REAL :: + & ct,st,cp,sp,ct2,st2,cp2,sp2,r2,gxx,gxy,gxz,gyy,gyz,gzz,dgxx, + & dgxy,dgxz,dgyy,dgyz,dgzz + CCTK_REAL,PARAMETER :: + & two = 2.0D0 + +c ------------------------------------------------------------------ + + Nt = SIZE(theta) + Np = SIZE(phi) + + r2 = r*r + + DO it = 1, Nt + + ct = COS(theta(it)) ; ct2 = ct*ct + st = SIN(theta(it)) ; st2 = st*st + + DO ip = 1, Np + + cp = COS(phi(ip)) ; cp2 = cp*cp + sp = SIN(phi(ip)) ; sp2 = sp*sp + + gxx = gxxs(it,ip) ; dgxx = dgxxs(it,ip) + gxy = gxys(it,ip) ; dgxy = dgxys(it,ip) + gxz = gxzs(it,ip) ; dgxz = dgxzs(it,ip) + gyy = gyys(it,ip) ; dgyy = dgyys(it,ip) + gyz = gyzs(it,ip) ; dgyz = dgyzs(it,ip) + gzz = gzzs(it,ip) ; dgzz = dgzzs(it,ip) + + grrs(it,ip) = st2*cp2*gxx+st2*sp2*gyy+ct2*gzz + & +two*st2*cp*sp*gxy+two*st*cp*ct*gxz+two*st*ct*sp*gyz + + grts(it,ip) = r*(st*cp2*ct*gxx+two*st*ct*sp*cp*gxy + & +cp*(ct2-st2)*gxz+st*sp2*ct*gyy+sp*(ct2-st2)*gyz-ct*st*gzz) + + grps(it,ip) = r*st*(-st*sp*cp*gxx-st*(sp2-cp2)*gxy-sp*ct*gxz + & +st*sp*cp*gyy+ct*cp*gyz) + + gtts(it,ip) = r2*(ct2*cp2*gxx+two*ct2*sp*cp*gxy + & -two*st*ct*cp*gxz+ct2*sp2*gyy-two*st*sp*ct*gyz+st2*gzz) + + gtps(it,ip) = r2*st*(-cp*sp*ct*gxx-ct*(sp2-cp2)*gxy + & +st*sp*gxz+cp*sp*ct*gyy-st*cp*gyz) + + gpps(it,ip) = r2*st2*(sp2*gxx-two*cp*sp*gxy+cp2*gyy) + + dgtts(it,ip) = two/r*gtts(it,ip)+r2*(ct2*cp2*dgxx + & +two*ct2*sp*cp*dgxy-two*st*ct*cp*dgxz+ct2*sp2*dgyy + & -two*st*sp*ct*dgyz+st2*dgzz) + + dgtps(it,ip) = two/r*gtps(it,ip)+r2*st*(-cp*sp*ct*dgxx + & -ct*(sp2-cp2)*dgxy+st*sp*dgxz+cp*sp*ct*dgyy-st*cp*dgyz) + + dgpps(it,ip) = two/r*gpps(it,ip) +r2*st2*(sp2*dgxx + & -two*cp*sp*dgxy+cp2*dgyy) + + END DO + + END DO + + + END SUBROUTINE cartesian_to_spherical + + + + + + + + + + diff --git a/src/cartesian_to_spherical_int.F b/src/cartesian_to_spherical_int.F new file mode 100644 index 0000000..a4db3b2 --- /dev/null +++ b/src/cartesian_to_spherical_int.F @@ -0,0 +1,24 @@ + +#include "cctk.h" + + MODULE cartesian_to_spherical_int + +c ------------------------------------------------------------------ + + INTERFACE + + SUBROUTINE cartesian_to_spherical(theta,phi,rad,gxxs,gxys,gxzs, + & gyys,gyzs,gzzs,grrs,grts,grps,gtts,gtps,gpps,dgxxs,dgxys, + & dgxzs,dgyys,dgyzs,dgzzs,dgtts,dgtps,dgpps) + IMPLICIT NONE + CCTK_REAL :: + & rad,theta(:),phi(:) + CCTK_REAL, DIMENSION (:,:) :: + & gxxs,gxys,gxzs,gyys,gyzs,gzzs,dgxxs,dgxys,dgxzs,dgyys, + & dgyzs,dgzzs,grrs,grts,grps,gtts,gtps,gpps,dgtts,dgtps, + & dgpps + END SUBROUTINE + + END INTERFACE + + END MODULE cartesian_to_spherical_int diff --git a/src/energy.f b/src/energy.f new file mode 100644 index 0000000..19ef5e2 --- /dev/null +++ b/src/energy.f @@ -0,0 +1,205 @@ + PROGRAM energy + +c ----------------------------------------------------------------- +c +c Calculate energy flux from even parity variable (only real part) +c +c Now copes with having unequally spaced timesteps, it works out +c the smallest timestep, and uses this to interpolate to a new +c array with equally spaced timesteps. +c +c ----------------------------------------------------------------- + + IMPLICIT NONE + + INTEGER,PARAMETER :: nmax=5000 + INTEGER :: i,itotal,it1,it2,Nt + CHARACTER*50 filename + REAL*8 t(nmax),Qplus(nmax),dEdt(nmax),E,ddQ(nmax),Qnew(nmax) + REAL*8 Dtnew,fac,time,t1,t2,dQdt_l,dQdt_r,dtmin,dt + + WRITE(6,*) + WRITE(6,*) "Starting PROGRAM energy ..." + WRITE(6,*) "---------------------------" + + WRITE(6,131) + 131 FORMAT("Filename containing Q ? ",$) + READ(5,*) filename + WRITE(6,*) "Filename is ... ",filename + + OPEN(UNIT=12,FILE=filename,STATUS="old") + + DO i=1,nmax + READ(12,*,ERR=101) t(i),Qplus(i) + itotal = i + END DO + WRITE(6,*) "Didn't read all of data file" + STOP + + 101 CLOSE(12) + +c Give information about the data + WRITE(6,*) "We have data from t=",t(1)," to ",t(itotal) + +c Find minimum time step + dtmin = 100000.0 + do i=2,itotal + dt = t(i)-t(i-1) + if (dt<dtmin) dtmin = dt + end do + WRITE(6,*) "Minimum timestep is ",dtmin + +c Initial time for energy flux + 201 WRITE(6,132) + 132 FORMAT(" Initial time ? ",$) + READ(5,*) t1 + IF (t1 < t(1)) GOTO 201 + +c Calculate initial timeslice + do i=1,itotal + if (t(i)<t1) it1=i + end do + WRITE(6,*) "Initial timeslice is ...",it1 + +c Final time for energy flux + WRITE(6,133) + 133 FORMAT(" Final time ? ",$) + 301 READ(5,*) t2 + IF (t2 > t(itotal)) GOTO 301 + +c Calculate final timeslice + do i=1,itotal + if (t(i)<t2) it2=i + end do + WRITE(6,*) "Final timeslice is ...",it2 + + + Nt = (t2-t1)/dtmin + WRITE(6,*) "Interpolating to ",Nt," timesteps" + + +c Do a spline + Dtnew = Dtmin + dQdt_l = ( Qplus(2) - Qplus(1) ) / + & ( t(2) - t(1) ) + dQdt_r = (Qplus(itotal) - Qplus(itotal-1) ) / + & (t(itotal) - t(itotal-1)) + CALL SPLINE(t,Qplus,itotal,dQdt_l,dQdt_r,ddQ) + DO i = 1,Nt + time = t1+DBLE(i-1)*Dtnew + CALL SPLINT(t,Qplus,ddQ,itotal,time,Qnew(i)) + ENDDO + +c This is the normalisation used + fac = 1.0D0/32.0D0/ACOS(-1.0D0) + +c Calculate energy flux + dEdt(1) = fac*((Qnew(2)-Qnew(1))/Dtnew)**2 + DO i = 2,Nt-1 + dEdt(i) = fac*((Qnew(i+1)-Qnew(i-1))*0.5D0/Dtnew)**2 + END DO + dEdt(Nt) = fac*((Qnew(Nt)-Qnew(Nt-1))/Dtnew)**2 + +c Write results to file + OPEN(UNIT=101,FILE="dEdt",STATUS="unknown") + DO i = 1,Nt + time = t1+DBLE(i-1)*Dtnew + WRITE(101,*) time,dEdt(i) + END DO + CLOSE(101) + +c Use Simpsons rule to calculate integral + IF ( MOD(Nt,2) /= 0) THEN + WRITE(6,*) "Using ",Nt," timeslices. Removing one" + WRITE(6,*) "to give Simpsons rule an odd number ..." + END IF + + E = dEdt(1) + 4.0D0*dEdt(Nt-1)+dEdt(Nt) + DO i=2,Nt-3,2 + E = E + 4.0D0*dEdt(i) + 2.0D0*dEdt(i+1) + ENDDO + E = E*Dtnew/3.0D0 + +c Write out total energy + WRITE(6,*) + WRITE(6,*) "Total energy in this mode is : ",E + WRITE(6,*) "****************************" + + END PROGRAM energy + + + + + SUBROUTINE spline(x,y,n,yp1,ypn,y2) + + INTEGER, PARAMETER :: nmax=5000 + + INTEGER n,i,k + REAL*8 yp1,ypn,x(nmax),y(nmax),y2(nmax) + REAL*8 p,qn,sig,un,u(nmax) + + IF (yp1.GT..99d30) THEN + y2(1)=0.d0 + u(1)=0.d0 + ELSE + y2(1)=-0.5d0 + u(1)=(3.d0/(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1) + ENDIF + + DO i=2,n-1 + sig = (x(i)-x(i-1))/(x(i+1)-x(i-1)) + p = sig*y2(i-1)+2.d0 + y2(i) = (sig-1.d0)/p + u(i) = (6.d0*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1)) + & /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p + ENDDO + + IF (ypn.GT..99d30) then + qn=0.d0 + un=0.d0 + ELSE + qn=0.5d0 + un=(3.d0/(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1))) + ENDIF + + y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.d0) + + DO k=n-1,1,-1 + y2(k)=y2(k)*y2(k+1)+u(k) + ENDDO + + RETURN + END + + + + SUBROUTINE splint(xa,ya,y2a,n,x,y) + + INTEGER n + REAL*8 x,y,xa(n),y2a(n),ya(n) + INTEGER k,khi,klo + REAL*8 a,b,h + + klo=1 + khi=n +1 if (khi-klo.gt.1) then + k=(khi+klo)/2 + if(xa(k).gt.x)then + khi=k + else + klo=k + endif + goto 1 + endif + h=xa(khi)-xa(klo) + if (h.eq.0.d0) pause 'bad xa input in splint' + a=(xa(khi)-x)/h + b=(x-xa(klo))/h + y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h** + *2)/6.d0 + + RETURN + END + + + diff --git a/src/make.code.defn b/src/make.code.defn new file mode 100644 index 0000000..a0b7189 --- /dev/null +++ b/src/make.code.defn @@ -0,0 +1,27 @@ +# Main make.code.defn file for thorn Extract +# $Header$ + +# Source files in this directory +SRCS = D2_extract_int.F \ + D3_extract_int.F \ + D3_to_D2_int.F \ + cartesian_to_spherical_int.F \ + unphysical_to_physical_int.F \ + ADMmass_integrand3D_int.F \ + momentum_integrand3D_int.F\ + spin_integrand3D_int.F\ + met_rad_der_int.F \ + Extract.F \ + D3_extract.F \ + cartesian_to_spherical.F \ + unphysical_to_physical.F \ + D2_extract.F \ + met_rad_der.F \ + D3_to_D2.F \ + ADMmass_integrand3D.F \ + momentum_integrand3D.F \ + spin_integrand3D.F + +# Subdirectories containing source files +SUBDIRS = + diff --git a/src/make.code.deps b/src/make.code.deps new file mode 100644 index 0000000..126c2e9 --- /dev/null +++ b/src/make.code.deps @@ -0,0 +1,2 @@ +Extract.o: D3_extract_int.o +ADMmass_integrand3D.o: ADMmass_integrand_int3D.o diff --git a/src/met_rad_der.F b/src/met_rad_der.F new file mode 100644 index 0000000..31821f0 --- /dev/null +++ b/src/met_rad_der.F @@ -0,0 +1,91 @@ + +#include "cctk.h" + +c ======================================================================== + + SUBROUTINE met_rad_der(origin,Dx,Dy,Dz,x,y,z,g,dg) + +c ------------------------------------------------------------------------ +c +c Calculate isotropic radial derivative of unphysical metric, +c component or conformal factor at each point on the 3D grid, +c using derivatives in the Cartesian directions. +c +c ------------------------------------------------------------------------ + + IMPLICIT NONE + +c Input variables + CCTK_REAL,INTENT(IN) :: + & Dx,Dy,Dz,origin(3) + CCTK_REAL,DIMENSION(:),INTENT(IN) :: + & x,y,z + CCTK_REAL,DIMENSION(:,:,:),INTENT(IN) :: + & g + +c Output variables + CCTK_REAL,DIMENSION(:,:,:),INTENT(OUT) :: + & dg + +c Local variables, here only + INTEGER :: + & i,j,k + CCTK_REAL,PARAMETER :: + & zero = 0.0D0, + & half = 0.5D0 + CCTK_REAL :: + & rad + +c ------------------------------------------------------------------------ + + DO k = 2, SIZE(z)-1 + DO j = 2, SIZE(y)-1 + DO i = 2, SIZE(x)-1 + + rad = SQRT((x(i)-origin(1))**2 + & +(y(j)-origin(2))**2 + & +(z(k)-origin(3))**2) + + IF (rad.NE.0) THEN + + dg(i,j,k) = half/rad* + & ((x(i)-origin(1))/Dx*(g(i+1,j,k)-g(i-1,j,k)) + & +(y(j)-origin(2))/Dy*(g(i,j+1,k)-g(i,j-1,k)) + & +(z(k)-origin(3))/Dz*(g(i,j,k+1)-g(i,j,k-1))) + + ELSE + + dg(i,j,k) = zero + + ENDIF + + ENDDO + ENDDO + ENDDO + +c This is needed when the grid is an octant, but it does not hurt +c if it is not + + DO k = 2, size(z)-1 + DO j = 2, size(y)-1 + dg(1,j,k) = dg(2,j,k) + ENDDO + ENDDO + DO k = 2, size(z)-1 + DO i = 1, size(x)-1 + dg(i,1,k) = dg(i,2,k) + ENDDO + ENDDO + DO j = 1, size(y)-1 + DO i = 1, size(x)-1 + dg(i,j,1) = dg(i,j,2) + ENDDO + ENDDO + + END SUBROUTINE met_rad_der + + + + + + diff --git a/src/met_rad_der_int.F b/src/met_rad_der_int.F new file mode 100644 index 0000000..d7b2f32 --- /dev/null +++ b/src/met_rad_der_int.F @@ -0,0 +1,24 @@ + +#include "cctk.h" + + MODULE met_rad_der_int + +c ------------------------------------------------------------------ + + INTERFACE + + SUBROUTINE met_rad_der(origin,Dx,Dy,Dz,x,y,z,g,dg) + IMPLICIT NONE + CCTK_REAL,INTENT(IN) :: + & Dx,Dy,Dz,origin(3) + CCTK_REAL,DIMENSION(:),INTENT(IN) :: + & x,y,z + CCTK_REAL,DIMENSION(:,:,:),INTENT(IN) :: + & g + CCTK_REAL,DIMENSION(:,:,:),INTENT(OUT) :: + & dg + END SUBROUTINE met_rad_der + + END INTERFACE + + END MODULE met_rad_der_int diff --git a/src/momentum_integrand3D.F b/src/momentum_integrand3D.F new file mode 100644 index 0000000..011fa68 --- /dev/null +++ b/src/momentum_integrand3D.F @@ -0,0 +1,188 @@ + +#include "cctk.h" + +c ======================================================================== + + SUBROUTINE momentum_integrand3D(origin,Dx,Dy,Dz,x,y,z,gxx,gxy,gxz, + & gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz, + & momentum_int,Psi,Psi_power) + +c ------------------------------------------------------------------------ +c +c Estimates the momentum at a given radius using Equation (11.2.14) from +c Walds General Relativity +c +c ----------------------------------------------------------------------- + + IMPLICIT NONE + +c Input variables + INTEGER,INTENT(IN) :: + & Psi_power + CCTK_REAL,INTENT(IN) :: + & Dx,Dy,Dz,origin(3) + CCTK_REAL,DIMENSION(:),INTENT(IN) :: + & x,y,z + CCTK_REAL,DIMENSION(:,:,:),INTENT(IN) :: + & gxx,gxy,gxz,gyy,gyz,gzz,Psi, + & hxx,hxy,hxz,hyy,hyz,hzz + + +c Output variables + CCTK_REAL,DIMENSION(:,:,:),INTENT(OUT) :: + & momentum_int + +c Local variables, here only + INTEGER :: + & i,j,k,ip,count + + CCTK_REAL,PARAMETER :: + & half = 0.5D0 + CCTK_REAL :: + & rad,ux,uy,uz,det,dxx,dxy,dxz,dyy,dyz,dzz,uxx,uxy,uxz,uyy, + & uyz,uzz, + & Pi,idet,p,term1, + & tracek + + data count / 1 / + save count +c ------------------------------------------------------------------------ + + Pi = ACOS(-1D0) + +c Because other codes evolve Psi**4 + SELECT CASE (Psi_power) + + CASE (1) + ip = 4 + + CASE (4) + ip = 1 + + CASE DEFAULT + WRITE(*,*) "This value of Psi_power is not supported" + + END SELECT + + + DO k = 2, SIZE(z)-1 + DO j = 2, SIZE(y)-1 + DO i = 2, SIZE(x)-1 + + rad = SQRT((x(i)-origin(1))**2 + & +(y(j)-origin(2))**2 + & +(z(k)-origin(3))**2) + + IF (rad.NE.0) THEN + + ux = (x(i)-origin(1))/rad + uy = (y(j)-origin(2))/rad + uz = (z(k)-origin(3))/rad + +c Abbreviations for metric coefficients +c ------------------------------------- + p = psi(i,j,k)**ip + + dxx = p*gxx(i,j,k); dxy = p*gxy(i,j,k) + dxz = p*gxz(i,j,k); dyy = p*gyy(i,j,k) + dyz = p*gyz(i,j,k); dzz = p*gzz(i,j,k) + +c Determinant of 3-metric +c ----------------------- + det = (dxx*dyy*dzz + 2.0D0*dxy*dxz*dyz + & - (dxx*dyz**2 + dyy*dxz**2 + dzz*dxy**2)) + + idet = 1.0/det + +c Inverse 3-metric +c ---------------- + uxx = idet*(dyy*dzz - dyz**2) + uyy = idet*(dxx*dzz - dxz**2) + uzz = idet*(dxx*dyy - dxy**2) + uxy = idet*(dxz*dyz - dzz*dxy) + uxz = idet*(dxy*dyz - dyy*dxz) + uyz = idet*(dxy*dxz - dxx*dyz) + + +c Trace of extrinsic curvature +c ---------------------------- + + tracek = uxx*hxx(i,j,k) + uyy*hyy(i,j,k) + + & uzz*hzz(i,j,k) + + + +c Integrands +c ---------- + + if (count.eq.1) then + term1 = ux*hxx(i,j,k)+uy*hxy(i,j,k)+uz*hxz(i,j,k) + +c momentum_int(i,j,k) = 1.0D0/8.0D0/Pi* +c & (term1 - tracek*ux)*rad**2 + momentum_int(i,j,k) = 1.0D0/8.0D0/Pi* + & (term1 - tracek*(dxx*ux+dxy*uy+dxz*uz))*rad**2 + + else if (count.eq.2) then + term1 = ux*hxy(i,j,k)+uy*hyy(i,j,k)+uz*hyz(i,j,k) + +c momentum_int(i,j,k) = 1.0D0/8.0D0/Pi* +c & (term1 - tracek*uy)*rad**2 + momentum_int(i,j,k) = 1.0D0/8.0D0/Pi* + & (term1 - tracek*(dxy*ux+dyy*uy+dyz*uz))*rad**2 + + else if (count.eq.3) then + term1 = ux*hxz(i,j,k)+uy*hyz(i,j,k)+uz*hzz(i,j,k) + +c momentum_int(i,j,k) = 1.0D0/8.0D0/Pi* +c & (term1 - tracek*uz)*rad**2 + momentum_int(i,j,k) = 1.0D0/8.0D0/Pi* + & (term1 - tracek*(dxz*ux+dyz*uy+dzz*uz))*rad**2 + end if + + + ELSE + + momentum_int(i,j,k) = 0.0D0 + + ENDIF + + ENDDO + ENDDO + ENDDO + + +c This is needed when the grid is an octant, but it does not hurt +c if it is not + + DO k = 2, size(z)-1 + DO j = 2, size(y)-1 + momentum_int(1,j,k) = momentum_int(2,j,k) + ENDDO + ENDDO + DO k = 2, size(z)-1 + DO i = 1, size(x)-1 + momentum_int(i,1,k) = momentum_int(i,2,k) + ENDDO + ENDDO + DO j = 1, size(y)-1 + DO i = 1, size(x)-1 + momentum_int(i,j,1) = momentum_int(i,j,2) + ENDDO + ENDDO + + +c setting counter +c --------------- + + if (count.eq.3) then + count = 1 + else + count = count + 1 + end if + + +c end +c --- + + END SUBROUTINE momentum_integrand3D diff --git a/src/momentum_integrand3D_int.F b/src/momentum_integrand3D_int.F new file mode 100644 index 0000000..e97174f --- /dev/null +++ b/src/momentum_integrand3D_int.F @@ -0,0 +1,30 @@ + +#include "cctk.h" + + MODULE momentum_integrand3D_int + +c ------------------------------------------------------------------ + + INTERFACE + + SUBROUTINE momentum_integrand3D(origin,Dx,Dy,Dz,x,y,z,gxx, + & gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,momentum_int, + & Psi,Psi_power) + + IMPLICIT NONE + INTEGER,INTENT(IN) :: + & Psi_power + CCTK_REAL,INTENT(IN) :: + & Dx,Dy,Dz,origin(3) + CCTK_REAL,DIMENSION(:),INTENT(IN) :: + & x,y,z + CCTK_REAL,DIMENSION(:,:,:),INTENT(IN) :: + & gxx,gxy,gxz,gyy,gyz,gzz,Psi, + & hxx,hxy,hxz,hyy,hyz,hzz + CCTK_REAL,DIMENSION(:,:,:),INTENT(OUT) :: + & momentum_int + END SUBROUTINE momentum_integrand3D + + END INTERFACE + + END MODULE momentum_integrand3D_int diff --git a/src/spin_integrand3D.F b/src/spin_integrand3D.F new file mode 100644 index 0000000..3f1afc3 --- /dev/null +++ b/src/spin_integrand3D.F @@ -0,0 +1,193 @@ + +#include "cctk.h" + +c ======================================================================== + + SUBROUTINE spin_integrand3D(origin,Dx,Dy,Dz,x,y,z,gxx,gxy,gxz, + & gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz, + & spin_int,Psi,Psi_power) + +c ------------------------------------------------------------------------ +c +c Estimates the spin at a given radius using Equation (12) from +c <Time-asymmetric initial data for black holes and black-hole +c collisions>, Bowen and York, Phys. Rev. D, 21 2047(1980) +c +c ----------------------------------------------------------------------- + + IMPLICIT NONE + +c Input variables + INTEGER,INTENT(IN) :: + & Psi_power + CCTK_REAL,INTENT(IN) :: + & Dx,Dy,Dz,origin(3) + CCTK_REAL,DIMENSION(:),INTENT(IN) :: + & x,y,z + CCTK_REAL,DIMENSION(:,:,:),INTENT(IN) :: + & gxx,gxy,gxz,gyy,gyz,gzz,Psi, + & hxx,hxy,hxz,hyy,hyz,hzz + + +c Output variables + CCTK_REAL,DIMENSION(:,:,:),INTENT(OUT) :: + & spin_int + +c Local variables, here only + INTEGER :: + & i,j,k,ip,count + + CCTK_REAL,PARAMETER :: + & half = 0.5D0 + CCTK_REAL :: + & rad,ux,uy,uz,det,dxx,dxy,dxz,dyy,dyz,dzz,uxx,uxy,uxz,uyy, + & uyz,uzz, + & Pi,idet,p, + & term1,term2 + + data count / 1 / + save count +c ------------------------------------------------------------------------ + + Pi = ACOS(-1D0) + + SELECT CASE (Psi_power) + + CASE (1) + ip = 4 + + CASE (4) + ip = 1 + + CASE DEFAULT + WRITE(*,*) "This value of Psi_power is not supported" + + END SELECT + + + DO k = 2, SIZE(z)-1 + DO j = 2, SIZE(y)-1 + DO i = 2, SIZE(x)-1 + + rad = SQRT((x(i)-origin(1))**2 + & +(y(j)-origin(2))**2 + & +(z(k)-origin(3))**2) + + IF (rad.NE.0) THEN + + ux = (x(i)-origin(1))/rad + uy = (y(j)-origin(2))/rad + uz = (z(k)-origin(3))/rad + +c Abbreviations for metric coefficients +c ------------------------------------- + p = psi(i,j,k)**ip + + dxx = p*gxx(i,j,k); dxy = p*gxy(i,j,k) + dxz = p*gxz(i,j,k); dyy = p*gyy(i,j,k) + dyz = p*gyz(i,j,k); dzz = p*gzz(i,j,k) + +c Determinant of 3-metric +c ----------------------- + det = (dxx*dyy*dzz + 2.0D0*dxy*dxz*dyz + & - (dxx*dyz**2 + dyy*dxz**2 + dzz*dxy**2)) + + idet = 1.0/det + +c Inverse 3-metric +c ---------------- + uxx = idet*(dyy*dzz - dyz**2) + uyy = idet*(dxx*dzz - dxz**2) + uzz = idet*(dxx*dyy - dxy**2) + uxy = idet*(dxz*dyz - dzz*dxy) + uxz = idet*(dxy*dyz - dyy*dxz) + uyz = idet*(dxy*dxz - dxx*dyz) + + +c Integrands +c ---------- + + if (count.eq.1) then + term1 = uxz*(ux*hxx(i,j,k)+uy*hxy(i,j,k)+uz*hxz(i,j,k)) + & +uyz*(ux*hxy(i,j,k)+uy*hyy(i,j,k)) + & +uzz*(ux*hxz(i,j,k)+uy*hyz(i,j,k)+uz*hzz(i,j,k)) + + term2 = uxy*(ux*hxx(i,j,k)+uy*hxy(i,j,k)+uz*hxz(i,j,k)) + & +uyy*(ux*hxy(i,j,k)+uy*hyy(i,j,k)+uz*hyz(i,j,k)) + & +uyz*(ux*hxz(i,j,k)+uz*hzz(i,j,k)) + + spin_int(i,j,k) = 1.0D0/8.0D0/Pi* + & (uy*term1-uz*term2)*rad**3 + + else if (count.eq.2) then + term1 = uxx*(ux*hxx(i,j,k)+uy*hxy(i,j,k)+uz*hxz(i,j,k)) + & +uxy*(ux*hxy(i,j,k)+uy*hyy(i,j,k)+uz*hyz(i,j,k)) + & +uxz*(uy*hyz(i,j,k)+uz*hzz(i,j,k)) + + term2 = uxz*(ux*hxx(i,j,k)+uy*hxy(i,j,k)) + & +uyz*(ux*hxy(i,j,k)+uy*hyy(i,j,k)+uz*hyz(i,j,k)) + & +uzz*(ux*hxz(i,j,k)+uy*hyz(i,j,k)+uz*hzz(i,j,k)) + + spin_int(i,j,k) = 1.0D0/8.0D0/Pi* + & (uz*term1-ux*term2)*rad**3 + + else if (count.eq.3) then + term1 = uxy*(ux*hxx(i,j,k)+uz*hxz(i,j,k)) + & +uyy*(ux*hxy(i,j,k)+uy*hyy(i,j,k)+uz*hyz(i,j,k)) + & +uyz*(ux*hxz(i,j,k)+uy*hyz(i,j,k)+uz*hzz(i,j,k)) + + term2 = uxx*(ux*hxx(i,j,k)+uy*hxy(i,j,k)+uz*hxz(i,j,k)) + & +uxy*(uy*hyy(i,j,k)+uz*hyz(i,j,k)) + & +uxz*(ux*hxz(i,j,k)+uy*hyz(i,j,k)+uz*hzz(i,j,k)) + + spin_int(i,j,k) = 1.0D0/8.0D0/Pi* + & (ux*term1-uy*term2)*rad**3 + + end if + + + ELSE + + spin_int(i,j,k) = 0.0D0 + + ENDIF + + ENDDO + ENDDO + ENDDO + + +c This is needed when the grid is an octant, but it does not hurt +c if it is not + + DO k = 2, size(z)-1 + DO j = 2, size(y)-1 + spin_int(1,j,k) = spin_int(2,j,k) + ENDDO + ENDDO + DO k = 2, size(z)-1 + DO i = 1, size(x)-1 + spin_int(i,1,k) = spin_int(i,2,k) + ENDDO + ENDDO + DO j = 1, size(y)-1 + DO i = 1, size(x)-1 + spin_int(i,j,1) = spin_int(i,j,2) + ENDDO + ENDDO + + +c setting counter +c --------------- + + if (count.eq.3) then + count = 1 + else + count = count + 1 + end if + + +c end +c --- + + END SUBROUTINE spin_integrand3D diff --git a/src/spin_integrand3D_int.F b/src/spin_integrand3D_int.F new file mode 100644 index 0000000..cb2a1ca --- /dev/null +++ b/src/spin_integrand3D_int.F @@ -0,0 +1,30 @@ + +#include "cctk.h" + + MODULE spin_integrand3D_int + +c ------------------------------------------------------------------ + + INTERFACE + + SUBROUTINE spin_integrand3D(origin,Dx,Dy,Dz,x,y,z,gxx, + & gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,spin_int, + & Psi,Psi_power) + + IMPLICIT NONE + INTEGER,INTENT(IN) :: + & Psi_power + CCTK_REAL,INTENT(IN) :: + & Dx,Dy,Dz,origin(3) + CCTK_REAL,DIMENSION(:),INTENT(IN) :: + & x,y,z + CCTK_REAL,DIMENSION(:,:,:),INTENT(IN) :: + & gxx,gxy,gxz,gyy,gyz,gzz,Psi, + & hxx,hxy,hxz,hyy,hyz,hzz + CCTK_REAL,DIMENSION(:,:,:),INTENT(OUT) :: + & spin_int + END SUBROUTINE spin_integrand3D + + END INTERFACE + + END MODULE spin_integrand3D_int diff --git a/src/unphysical_to_physical.F b/src/unphysical_to_physical.F new file mode 100644 index 0000000..74b4ad8 --- /dev/null +++ b/src/unphysical_to_physical.F @@ -0,0 +1,56 @@ + +#include "cctk.h" + +c ================================================================== + + SUBROUTINE unphysical_to_physical(grr,grt,grp,gtt,gtp,gpp,dgtt, + & dgtp,dgpp,Psis,dPsis,Psi_power) + +c ------------------------------------------------------------------ +c +c Convert unphysical metric components and radial (eta) derivatives +c on the 2-sphere into physical quantities, using the conformal factor +c on the sphere +c +c ------------------------------------------------------------------ + + IMPLICIT NONE + +c Input variables + INTEGER :: + & Psi_power + CCTK_REAL,INTENT(INOUT),DIMENSION (:,:) :: + & grr,grt,grp,gtt,gtp,gpp,dgtt,dgtp,dgpp + CCTK_REAL,INTENT(IN),DIMENSION (:,:) :: + & Psis,dPsis + +c Output variables +c WE ARE CHANGING THE INPUT VARIABLES !!!! + +c ------------------------------------------------------------------ + + dgtt = Psis**4*dgtt + 4.0d0*Psis**3*dPsis*gtt + dgtp = Psis**4*dgtp + 4.0d0*Psis**3*dPsis*gtp + dgpp = Psis**4*dgpp + 4.0d0*Psis**3*dPsis*gpp + + grr = Psis**4*grr + grt = Psis**4*grt + grp = Psis**4*grp + gtt = Psis**4*gtt + gtp = Psis**4*gtp + gpp = Psis**4*gpp + + + END SUBROUTINE unphysical_to_physical + +c ================================================================== + + + + + + + + + + diff --git a/src/unphysical_to_physical_int.F b/src/unphysical_to_physical_int.F new file mode 100644 index 0000000..b53fa0b --- /dev/null +++ b/src/unphysical_to_physical_int.F @@ -0,0 +1,21 @@ + +#include "cctk.h" + + MODULE unphysical_to_physical_int + +c ------------------------------------------------------------------ + + INTERFACE + + SUBROUTINE unphysical_to_physical(grr,grt,grp,gtt,gtp,gpp,dgtt, + & dgtp,dgpp,Psis,dPsis,Psi_power) + IMPLICIT NONE + INTEGER :: + & Psi_power + CCTK_REAL, DIMENSION (:,:) :: + & grr,grt,grp,gtt,gtp,gpp,dgtt,dgtp,dgpp,Psis,dPsis + END SUBROUTINE + + END INTERFACE + + END MODULE unphysical_to_physical_int diff --git a/test/test_extract/Qeven_R1_20.tl b/test/test_extract/Qeven_R1_20.tl new file mode 100644 index 0000000..92cc299 --- /dev/null +++ b/test/test_extract/Qeven_R1_20.tl @@ -0,0 +1,6 @@ + 0.000000000000000E+00 -0.501417076884790E-05 0.000000000000000E+00 + 0.600000000000000E-01 -0.500918802832626E-05 0.000000000000000E+00 + 0.120000000000000E+00 -0.499425240358436E-05 0.000000000000000E+00 + 0.180000000000000E+00 -0.496940984704886E-05 0.000000000000000E+00 + 0.240000000000000E+00 -0.493474798703217E-05 0.000000000000000E+00 + 0.300000000000000E+00 -0.489039494629224E-05 0.000000000000000E+00 diff --git a/test/test_extract/Qeven_R1_22.tl b/test/test_extract/Qeven_R1_22.tl new file mode 100644 index 0000000..2b1b2ba --- /dev/null +++ b/test/test_extract/Qeven_R1_22.tl @@ -0,0 +1,6 @@ + 0.000000000000000E+00 0.115920453067371E-15 0.000000000000000E+00 + 0.600000000000000E-01 0.102781388986612E-15 0.000000000000000E+00 + 0.120000000000000E+00 0.118592170804785E-15 0.000000000000000E+00 + 0.180000000000000E+00 0.954262370907658E-16 0.000000000000000E+00 + 0.240000000000000E+00 0.107387444437799E-15 0.000000000000000E+00 + 0.300000000000000E+00 0.100602776639667E-15 0.000000000000000E+00 diff --git a/test/test_extract/Qeven_R1_40.tl b/test/test_extract/Qeven_R1_40.tl new file mode 100644 index 0000000..a45f933 --- /dev/null +++ b/test/test_extract/Qeven_R1_40.tl @@ -0,0 +1,6 @@ + 0.000000000000000E+00 0.560853851605455E-02 0.000000000000000E+00 + 0.600000000000000E-01 0.560617197989429E-02 0.000000000000000E+00 + 0.120000000000000E+00 0.559908697023419E-02 0.000000000000000E+00 + 0.180000000000000E+00 0.558735649093210E-02 0.000000000000000E+00 + 0.240000000000000E+00 0.557114337412382E-02 0.000000000000000E+00 + 0.300000000000000E+00 0.555070111549437E-02 0.000000000000000E+00 diff --git a/test/test_extract/Qeven_R1_42.tl b/test/test_extract/Qeven_R1_42.tl new file mode 100644 index 0000000..330cb7c --- /dev/null +++ b/test/test_extract/Qeven_R1_42.tl @@ -0,0 +1,6 @@ + 0.000000000000000E+00 -0.557956412077577E-16 0.000000000000000E+00 + 0.600000000000000E-01 -0.574565961468893E-16 0.000000000000000E+00 + 0.120000000000000E+00 -0.509845022367271E-16 0.000000000000000E+00 + 0.180000000000000E+00 -0.626531478147805E-16 0.000000000000000E+00 + 0.240000000000000E+00 -0.486937393077706E-16 0.000000000000000E+00 + 0.300000000000000E+00 -0.585556606036106E-16 0.000000000000000E+00 diff --git a/test/test_extract/Qeven_R1_44.tl b/test/test_extract/Qeven_R1_44.tl new file mode 100644 index 0000000..834c739 --- /dev/null +++ b/test/test_extract/Qeven_R1_44.tl @@ -0,0 +1,6 @@ + 0.000000000000000E+00 0.335000108029692E-02 0.000000000000000E+00 + 0.600000000000000E-01 0.334858768116714E-02 0.000000000000000E+00 + 0.120000000000000E+00 0.334435621377972E-02 0.000000000000000E+00 + 0.180000000000000E+00 0.333735032646192E-02 0.000000000000000E+00 + 0.240000000000000E+00 0.332766736797455E-02 0.000000000000000E+00 + 0.300000000000000E+00 0.331545888621867E-02 0.000000000000000E+00 diff --git a/test/test_extract/Qeven_R2_20.tl b/test/test_extract/Qeven_R2_20.tl new file mode 100644 index 0000000..ba564fc --- /dev/null +++ b/test/test_extract/Qeven_R2_20.tl @@ -0,0 +1,6 @@ + 0.000000000000000E+00 0.172079966298455E-05 0.000000000000000E+00 + 0.600000000000000E-01 0.172286369448982E-05 0.000000000000000E+00 + 0.120000000000000E+00 0.172904539451029E-05 0.000000000000000E+00 + 0.180000000000000E+00 0.173933695449188E-05 0.000000000000000E+00 + 0.240000000000000E+00 0.175375502737258E-05 0.000000000000000E+00 + 0.300000000000000E+00 0.177233802642148E-05 0.000000000000000E+00 diff --git a/test/test_extract/Qeven_R2_22.tl b/test/test_extract/Qeven_R2_22.tl new file mode 100644 index 0000000..ed0538b --- /dev/null +++ b/test/test_extract/Qeven_R2_22.tl @@ -0,0 +1,6 @@ + 0.000000000000000E+00 -0.287250898579437E-16 0.000000000000000E+00 + 0.600000000000000E-01 -0.554035706647425E-16 0.000000000000000E+00 + 0.120000000000000E+00 -0.269165513003967E-16 0.000000000000000E+00 + 0.180000000000000E+00 -0.366840077302228E-16 0.000000000000000E+00 + 0.240000000000000E+00 -0.345853744116724E-16 0.000000000000000E+00 + 0.300000000000000E+00 -0.211775808648429E-16 0.000000000000000E+00 diff --git a/test/test_extract/Qeven_R2_40.tl b/test/test_extract/Qeven_R2_40.tl new file mode 100644 index 0000000..5a3445e --- /dev/null +++ b/test/test_extract/Qeven_R2_40.tl @@ -0,0 +1,6 @@ + 0.000000000000000E+00 0.378538965351061E-02 0.000000000000000E+00 + 0.600000000000000E-01 0.378388666783014E-02 0.000000000000000E+00 + 0.120000000000000E+00 0.377937493624180E-02 0.000000000000000E+00 + 0.180000000000000E+00 0.377189764739652E-02 0.000000000000000E+00 + 0.240000000000000E+00 0.376159488954885E-02 0.000000000000000E+00 + 0.300000000000000E+00 0.374870133328616E-02 0.000000000000000E+00 diff --git a/test/test_extract/Qeven_R2_42.tl b/test/test_extract/Qeven_R2_42.tl new file mode 100644 index 0000000..b5f1d86 --- /dev/null +++ b/test/test_extract/Qeven_R2_42.tl @@ -0,0 +1,6 @@ + 0.000000000000000E+00 0.321146551025695E-16 0.000000000000000E+00 + 0.600000000000000E-01 0.112692196817111E-16 0.000000000000000E+00 + 0.120000000000000E+00 0.218806619482187E-16 0.000000000000000E+00 + 0.180000000000000E+00 0.299201220713587E-16 0.000000000000000E+00 + 0.240000000000000E+00 0.412911813148731E-16 0.000000000000000E+00 + 0.300000000000000E+00 0.250949838171749E-16 0.000000000000000E+00 diff --git a/test/test_extract/Qeven_R2_44.tl b/test/test_extract/Qeven_R2_44.tl new file mode 100644 index 0000000..7fba470 --- /dev/null +++ b/test/test_extract/Qeven_R2_44.tl @@ -0,0 +1,6 @@ + 0.000000000000000E+00 0.226136999516570E-02 0.000000000000000E+00 + 0.600000000000000E-01 0.226047170343381E-02 0.000000000000000E+00 + 0.120000000000000E+00 0.225777517507269E-02 0.000000000000000E+00 + 0.180000000000000E+00 0.225330622988509E-02 0.000000000000000E+00 + 0.240000000000000E+00 0.224714859635179E-02 0.000000000000000E+00 + 0.300000000000000E+00 0.223944252755130E-02 0.000000000000000E+00 diff --git a/test/test_extract/mass_R1.tl b/test/test_extract/mass_R1.tl new file mode 100644 index 0000000..b990cad --- /dev/null +++ b/test/test_extract/mass_R1.tl @@ -0,0 +1,6 @@ + 0.000000000000000E+00 0.999971751592259E+00 + 0.600000000000000E-01 0.999936346286986E+00 + 0.120000000000000E+00 0.999830166604527E+00 + 0.180000000000000E+00 0.999653227395764E+00 + 0.240000000000000E+00 0.999405559305214E+00 + 0.300000000000000E+00 0.999087209145958E+00 diff --git a/test/test_extract/mass_R2.tl b/test/test_extract/mass_R2.tl new file mode 100644 index 0000000..d402068 --- /dev/null +++ b/test/test_extract/mass_R2.tl @@ -0,0 +1,6 @@ + 0.000000000000000E+00 0.999974184437137E+00 + 0.600000000000000E-01 0.999949423545773E+00 + 0.120000000000000E+00 0.999875155689259E+00 + 0.180000000000000E+00 0.999751386553946E+00 + 0.240000000000000E+00 0.999578150001170E+00 + 0.300000000000000E+00 0.999355512752157E+00 diff --git a/test/test_extract/rsch_R1.tl b/test/test_extract/rsch_R1.tl new file mode 100644 index 0000000..99b86b2 --- /dev/null +++ b/test/test_extract/rsch_R1.tl @@ -0,0 +1,6 @@ + 0.000000000000000E+00 0.359984303230721E+01 + 0.600000000000000E-01 0.359970447161759E+01 + 0.120000000000000E+00 0.359928878117140E+01 + 0.180000000000000E+00 0.359859590654562E+01 + 0.240000000000000E+00 0.359762576012367E+01 + 0.300000000000000E+00 0.359637822181613E+01 diff --git a/test/test_extract/rsch_R2.tl b/test/test_extract/rsch_R2.tl new file mode 100644 index 0000000..413344a --- /dev/null +++ b/test/test_extract/rsch_R2.tl @@ -0,0 +1,6 @@ + 0.000000000000000E+00 0.408324884253190E+01 + 0.600000000000000E-01 0.408314103935444E+01 + 0.120000000000000E+00 0.408281762485777E+01 + 0.180000000000000E+00 0.408227856993223E+01 + 0.240000000000000E+00 0.408152382854459E+01 + 0.300000000000000E+00 0.408055333770619E+01 diff --git a/test/test_extract_par b/test/test_extract_par new file mode 100644 index 0000000..d930492 --- /dev/null +++ b/test/test_extract_par @@ -0,0 +1,42 @@ +# DESC "Extract waveforms from BH" + + +ActiveThorns = "cartgrid3d time idanalyticBH pugh Extract Einstein interp ioutil" + +driver::global_nx = 15 +driver::global_ny = 15 +driver::global_nz = 15 + +grid::type = "bySpacing" +grid::domain = "octant" +grid::dxyz = 0.3 + +time::dtfac = 0.2 + +cactus::cctk_itlast = 5 +IO::out_every = 1 + +einstein::initial_data = "schwarzschild" +idanalyticbh::mass = 1 + +io::outdir = "./test_Extract" + +### Extraction parameters + +interp::order = 2 + +extract::Extract_num_detectors = 2 +extract::Extract_itout = 1 + +extract::Extract_Nt = 100 +extract::Extract_Np = 100 +extract::Extract_origin_x = 0. +extract::Extract_origin_y = 0. +extract::Extract_origin_z = 0. + +extract::Extract_all_modes = "yes" +extract::Extract_l_mode = 4 +extract::Extract_m_mode = 0 + +extract::Extract_detector1 = 2.5 +extract::Extract_detector2 = 3.0 |