diff options
| -rw-r--r-- | COPYING | 3 | ||||
| -rw-r--r-- | GPLv2 | 339 | ||||
| -rw-r--r-- | README | 14 | ||||
| -rw-r--r-- | doc/documentation.tex | 151 | ||||
| -rw-r--r-- | interface.ccl | 84 | ||||
| -rw-r--r-- | param.ccl | 103 | ||||
| -rw-r--r-- | schedule.ccl | 80 | ||||
| -rw-r--r-- | src/CheckParam.c | 15 | ||||
| -rw-r--r-- | src/Whisky_C2P2C.F90 | 131 | ||||
| -rw-r--r-- | src/Whisky_Con2Prim.F90 | 108 | ||||
| -rw-r--r-- | src/Whisky_Only_Atmo.F90 | 114 | ||||
| -rw-r--r-- | src/Whisky_ReadConformalData.F90 | 189 | ||||
| -rw-r--r-- | src/Whisky_ReconstructTest.F90 | 84 | ||||
| -rw-r--r-- | src/Whisky_ShockTube.F90 | 155 | ||||
| -rw-r--r-- | src/Whisky_SimpleWave.F90 | 238 | ||||
| -rw-r--r-- | src/Whisky_TOV.F90 | 324 | ||||
| -rw-r--r-- | src/make.code.defn | 16 | ||||
| -rw-r--r-- | src/make.configuration.defn | 4 | ||||
| -rw-r--r-- | src/make.configuration.deps | 8 | ||||
| -rw-r--r-- | src/tov_fcn.F | 101 | ||||
| -rw-r--r-- | src/tov_lsoda.F | 308 |
21 files changed, 2569 insertions, 0 deletions
@@ -0,0 +1,3 @@ +You can use this code under the GPL General Public License, version 2, or later + +The license text is included in the file GPLv2. @@ -0,0 +1,339 @@ + GNU GENERAL PUBLIC LICENSE + Version 2, June 1991 + + Copyright (C) 1989, 1991 Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 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 Lesser General Public License instead.) You can apply it to +your programs, too. + + When we speak of free software, we are referring to freedom, not +price. Our General Public Licenses are designed to make sure that you +have the freedom to distribute copies of free software (and charge for +this service if you wish), that you receive source code or can get it +if you want it, that you can change the software or use pieces of it +in new free programs; and that you know you can do these things. + + To protect your rights, we need to make restrictions that forbid +anyone to deny you these rights or to ask you to surrender the rights. +These restrictions translate to certain responsibilities for you if you +distribute copies of the software, or if you modify it. + + For example, if you distribute copies of such a program, whether +gratis or for a fee, you must give the recipients all the rights that +you have. You must make sure that they, too, receive or can get the +source code. And you must show them these terms so they know their +rights. + + We protect your rights with two steps: (1) copyright the software, and +(2) offer you this license which gives you legal permission to copy, +distribute and/or modify the software. + + Also, for each author's protection and ours, we want to make certain +that everyone understands that there is no warranty for this free +software. If the software is modified by someone else and passed on, we +want its recipients to know that what they have is not the original, so +that any problems introduced by others will not reflect on the original +authors' reputations. + + Finally, any free program is threatened constantly by software +patents. We wish to avoid the danger that redistributors of a free +program will individually obtain patent licenses, in effect making the +program proprietary. To prevent this, we have made it clear that any +patent must be licensed for everyone's free use or not licensed at all. + + The precise terms and conditions for copying, distribution and +modification follow. + + GNU GENERAL PUBLIC LICENSE + TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION + + 0. This License applies to any program or other work which contains +a notice placed by the copyright holder saying it may be distributed +under the terms of this General Public License. The "Program", below, +refers to any such program or work, and a "work based on the Program" +means either the Program or any derivative work under copyright law: +that is to say, a work containing the Program or a portion of it, +either verbatim or with modifications and/or translated into another +language. (Hereinafter, translation is included without limitation in +the term "modification".) Each licensee is addressed as "you". + +Activities other than copying, distribution and modification are not +covered by this License; they are outside its scope. The act of +running the Program is not restricted, and the output from the Program +is covered only if its contents constitute a work based on the +Program (independent of having been made by running the Program). +Whether that is true depends on what the Program does. + + 1. You may copy and distribute verbatim copies of the Program's +source code as you receive it, in any medium, provided that you +conspicuously and appropriately publish on each copy an appropriate +copyright notice and disclaimer of warranty; keep intact all the +notices that refer to this License and to the absence of any warranty; +and give any other recipients of the Program a copy of this License +along with the Program. + +You may charge a fee for the physical act of transferring a copy, and +you may at your option offer warranty protection in exchange for a fee. + + 2. You may modify your copy or copies of the Program or any portion +of it, thus forming a work based on the Program, and copy and +distribute such modifications or work under the terms of Section 1 +above, provided that you also meet all of these conditions: + + a) You must cause the modified files to carry prominent notices + stating that you changed the files and the date of any change. + + b) You must cause any work that you distribute or publish, that in + whole or in part contains or is derived from the Program or any + part thereof, to be licensed as a whole at no charge to all third + parties under the terms of this License. + + c) If the modified program normally reads commands interactively + when run, you must cause it, when started running for such + interactive use in the most ordinary way, to print or display an + announcement including an appropriate copyright notice and a + notice that there is no warranty (or else, saying that you provide + a warranty) and that users may redistribute the program under + these conditions, and telling the user how to view a copy of this + License. (Exception: if the Program itself is interactive but + does not normally print such an announcement, your work based on + the Program is not required to print an announcement.) + +These requirements apply to the modified work as a whole. If +identifiable sections of that work are not derived from the Program, +and can be reasonably considered independent and separate works in +themselves, then this License, and its terms, do not apply to those +sections when you distribute them as separate works. But when you +distribute the same sections as part of a whole which is a work based +on the Program, the distribution of the whole must be on the terms of +this License, whose permissions for other licensees extend to the +entire whole, and thus to each and every part regardless of who wrote it. + +Thus, it is not the intent of this section to claim rights or contest +your rights to work written entirely by you; rather, the intent is to +exercise the right to control the distribution of derivative or +collective works based on the Program. + +In addition, mere aggregation of another work not based on the Program +with the Program (or with a work based on the Program) on a volume of +a storage or distribution medium does not bring the other work under +the scope of this License. + + 3. You may copy and distribute the Program (or a work based on it, +under Section 2) in object code or executable form under the terms of +Sections 1 and 2 above provided that you also do one of the following: + + a) Accompany it with the complete corresponding machine-readable + source code, which must be distributed under the terms of Sections + 1 and 2 above on a medium customarily used for software interchange; or, + + b) Accompany it with a written offer, valid for at least three + years, to give any third party, for a charge no more than your + cost of physically performing source distribution, a complete + machine-readable copy of the corresponding source code, to be + distributed under the terms of Sections 1 and 2 above on a medium + customarily used for software interchange; or, + + c) Accompany it with the information you received as to the offer + to distribute corresponding source code. (This alternative is + allowed only for noncommercial distribution and only if you + received the program in object code or executable form with such + an offer, in accord with Subsection b above.) + +The source code for a work means the preferred form of the work for +making modifications to it. For an executable work, complete source +code means all the source code for all modules it contains, plus any +associated interface definition files, plus the scripts used to +control compilation and installation of the executable. However, as a +special exception, the source code distributed need not include +anything that is normally distributed (in either source or binary +form) with the major components (compiler, kernel, and so on) of the +operating system on which the executable runs, unless that component +itself accompanies the executable. + +If distribution of executable or object code is made by offering +access to copy from a designated place, then offering equivalent +access to copy the source code from the same place counts as +distribution of the source code, even though third parties are not +compelled to copy the source along with the object code. + + 4. You may not copy, modify, sublicense, or distribute the Program +except as expressly provided under this License. Any attempt +otherwise to copy, modify, sublicense or distribute the Program is +void, and will automatically terminate your rights under this License. +However, parties who have received copies, or rights, from you under +this License will not have their licenses terminated so long as such +parties remain in full compliance. + + 5. You are not required to accept this License, since you have not +signed it. However, nothing else grants you permission to modify or +distribute the Program or its derivative works. These actions are +prohibited by law if you do not accept this License. Therefore, by +modifying or distributing the Program (or any work based on the +Program), you indicate your acceptance of this License to do so, and +all its terms and conditions for copying, distributing or modifying +the Program or works based on it. + + 6. Each time you redistribute the Program (or any work based on the +Program), the recipient automatically receives a license from the +original licensor to copy, distribute or modify the Program subject to +these terms and conditions. You may not impose any further +restrictions on the recipients' exercise of the rights granted herein. +You are not responsible for enforcing compliance by third parties to +this License. + + 7. If, as a consequence of a court judgment or allegation of patent +infringement or for any other reason (not limited to patent issues), +conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot +distribute so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you +may not distribute the Program at all. For example, if a patent +license would not permit royalty-free redistribution of the Program by +all those who receive copies directly or indirectly through you, then +the only way you could satisfy both it and this License would be to +refrain entirely from distribution of the Program. + +If any portion of this section is held invalid or unenforceable under +any particular circumstance, the balance of the section is intended to +apply and the section as a whole is intended to apply in other +circumstances. + +It is not the purpose of this section to induce you to infringe any +patents or other property right claims or to contest validity of any +such claims; this section has the sole purpose of protecting the +integrity of the free software distribution system, which is +implemented by public license practices. Many people have made +generous contributions to the wide range of software distributed +through that system in reliance on consistent application of that +system; it is up to the author/donor to decide if he or she is willing +to distribute software through any other system and a licensee cannot +impose that choice. + +This section is intended to make thoroughly clear what is believed to +be a consequence of the rest of this License. + + 8. If the distribution and/or use of the Program is restricted in +certain countries either by patents or by copyrighted interfaces, the +original copyright holder who places the Program under this License +may add an explicit geographical distribution limitation excluding +those countries, so that distribution is permitted only in or among +countries not thus excluded. In such case, this License incorporates +the limitation as if written in the body of this License. + + 9. The Free Software Foundation may publish revised and/or new versions +of the General Public License from time to time. Such new versions will +be similar in spirit to the present version, but may differ in detail to +address new problems or concerns. + +Each version is given a distinguishing version number. If the Program +specifies a version number of this License which applies to it and "any +later version", you have the option of following the terms and conditions +either of that version or of any later version published by the Free +Software Foundation. If the Program does not specify a version number of +this License, you may choose any version ever published by the Free Software +Foundation. + + 10. If you wish to incorporate parts of the Program into other free +programs whose distribution conditions are different, write to the author +to ask for permission. For software which is copyrighted by the Free +Software Foundation, write to the Free Software Foundation; we sometimes +make exceptions for this. Our decision will be guided by the two goals +of preserving the free status of all derivatives of our free software and +of promoting the sharing and reuse of software generally. + + NO WARRANTY + + 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY +FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN +OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES +PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED +OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF +MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS +TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE +PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, +REPAIR OR CORRECTION. + + 12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING +WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR +REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, +INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING +OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED +TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY +YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER +PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE +POSSIBILITY OF SUCH DAMAGES. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Programs + + If you develop a new program, and you want it to be of the greatest +possible use to the public, the best way to achieve this is to make it +free software which everyone can redistribute and change under these terms. + + To do so, attach the following notices to the program. It is safest +to attach them to the start of each source file to most effectively +convey the exclusion of warranty; and each file should have at least +the "copyright" line and a pointer to where the full notice is found. + + <one line to give the program's name and a brief idea of what it does.> + Copyright (C) <year> <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., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +Also add information on how to contact you by electronic and paper mail. + +If the program is interactive, make it output a short notice like this +when it starts in an interactive mode: + + Gnomovision version 69, Copyright (C) year name of author + Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'. + This is free software, and you are welcome to redistribute it + under certain conditions; type `show c' for details. + +The hypothetical commands `show w' and `show c' should show the appropriate +parts of the General Public License. Of course, the commands you use may +be called something other than `show w' and `show c'; they could even be +mouse-clicks or menu items--whatever suits your program. + +You should also get your employer (if you work as a programmer) or your +school, if any, to sign a "copyright disclaimer" for the program, if +necessary. Here is a sample; alter the names: + + Yoyodyne, Inc., hereby disclaims all copyright interest in the program + `Gnomovision' (which makes passes at compilers) written by James Hacker. + + <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 Lesser General +Public License instead of this License. @@ -0,0 +1,14 @@ +Cactus Code Thorn Whisky_Init_Data +Author(s) : Luca Baiotti + Ian Hawke + Scott Hawley + Frank Löffler + Christian David Ott +Maintainer(s): Whisky developers team +Licence : GPLv2+ +-------------------------------------------------------------------------- + +1. Purpose + +Provides simple initial data (shock tubes, ...) for the Whisky code. + diff --git a/doc/documentation.tex b/doc/documentation.tex new file mode 100644 index 0000000..968e98a --- /dev/null +++ b/doc/documentation.tex @@ -0,0 +1,151 @@ +\documentclass{article} + +\begin{document} + +\title{Whisky\_Init\_Data} +\author{Luca Baiotti, Ian Hawke, Scott Hawley} +\date{24/07/2008} +\maketitle + +\abstract{Whisky\_Init\_Data - some initial data for {\tt Whisky}} + +\section{Introduction} +\label{sec:intro} + +This thorn generates some initial data for the Whisky code. There are more (and more physically +interesting) initial-data codes in other thorns. As with the Whisky code itself, please feel free to +add, alter or extend any part of this code. However please keep the documentation up to date (even, +or especially, if it's just to say what doesn't work). + +Currently this thorn contains a few tests that should really be test suites, some shock-tube +initial data, +%some (largely untested and unmaintained; please use the TOV initial-data codes in the +%dedicated thorns) TOV initial data solver routines, +a routine to set atmopshere everywhere on the +grid (for tests), a routine to read initial data from files (not very generic, tough) and a routine +to set up the simple-wave initial data . + + +\subsection{Tests} +\label{sec:tests} + +There are tests of the TVD reconstruction routine and of the routines +that convert between conservative and primitive variables. These all +just produce output to the screen or to {\tt fort.*} files. The +reconstruction test outputs the function to be reconstructed and the +boundary-extended values. The conservative-to-primitive test just +outputs the two sets of variables. If you haven't altered the code an if you set +\begin{verbatim} +eos_polytrope::eos_gamma = 2.0 +eos_polytrope::eos_k = 100.0 +\end{verbatim} +(which are the defaults), the output should be +% I checked that in 2008 the number were still right (but one gets them with eos_gamma=2.0 and +% eos_k ~ 0.1865) and committed the new ones for the default eos parameter values. + +\begin{verbatim} + primitive variables: + rho : 1.29047172182043 + velx : 9.902578465178671E-004 + vely : 9.902578465178671E-004 + velz : 9.902578465178671E-004 + eps : 0.374770481293314 + press : 166.531726481819 + w_lor : 1.00000147091915 +\end{verbatim} +The conservative to primitive to conservative test outputs the initial +and final data which should agree. + +\subsection{Shocktube tests} +\label{sec:shock} + +There are three possible shock-tube problems, referred to as {\tt Sod}, +{\tt Simple} and {\tt Blast}, with initial data +\begin{center} + \begin{tabular}[c]{|c|c|c|c|c|c|c|} + \hline Type & $\rho_{_L}$ & $v^i_{_L}$ & $\varepsilon_{_L}$ & $\rho_{_R}$ & $v^i_{_R}$ + & $\varepsilon_{_R}$ \\ \hline + Sod & 1 & 0 & 1.5 & 0.125 & 0 & 0.15 \\ + Simple & 10 & 0 & 20 & 1 & 0 & $10^{-6}$ \\ + Blast & 1 & 0 & 1500 & 1 & 0 & $1.5\cdot 10^{-2}$ \\ \hline + \end{tabular} +\end{center} +The shock shape can be planar (along each axis or along the main diagonal) or spherical and the +position of the plane or of the center of the sphere can be chosen though parameters. +If a diagonal shock is selected, the initial data is set to either the left or right +state depending on where the centre of the cell falls. Cleverer +routines that weight the initial data to avoid ``staircasing'' may be +added if there is demand. For more discussion on shock tubes +see~\cite{livrevsrrfd}. + + +%\subsection{TOV stars} +%\label{sec:tov} +% +%The Tolman-Oppenheimer-Volkoff solution is a spherically symmetric +%fluid ball matched to a Schwarzschild exterior. Typically an +%atmosphere is placed in the exterior to stop the equations of motion +%of the fluid being singular. Given an equation of state, the central +%density $\rho_c$ is specified. Then the solution is found by +%integrating the radial equations +%\begin{eqnarray} +% \label{eq:tov} +% \frac{\partial P(r)}{\partial r} & = & - \frac{(\rho + \rho \epsilon +% + P)(m + 4\pi r^3 P)}{r (r - 2m)} \\ +% \frac{\partial (\log \alpha (r))}{\partial r} & = & \frac{m + 4 \pi +% r^3 P}{r (r - 2m)} \\ +% \frac{\partial m(r)}{\partial r} & = & 4 \pi r^2 (\rho + \rho +% \epsilon) \\ +% \gamma_{rr}(r) & = & \left( 1-\frac{2m(r)}{r} \right)^{-1}, +%\end{eqnarray} +%where $m$ is the mass energy contained in a sphere radius $r$, +%$\gamma_{ij}$ the 3-metric, and $\alpha$ the lapse in standard +%Schwarzschild like coordinates. For more details see~\cite{hydro1}. +% +%The routines here, written by Scott Hawley, use the LSODA library to +%integrate the equations and then interpolate onto the Cartesian +%grid. +% +%This routine is untested and unmaintained: please use other thorns providing TOV initial data. They +%are located in dedicated thorns of this arrangement. + +\subsection{Only atmosphere} +\label{sec:only-atmo} + +For testing purposes, this routine sets all the points to the values of the atmosphere. + + +\subsection{Simple wave} +\label{sec:simple-wave} + +This routine stes initial data for a simple wave with sinusoidal initial function for the velocity, +as described in Anile, Miller, Motta, {\it Formation and damping of relativistic strong + shocks},Phys. Fluids {\bf 26}, 1450 (1983). + + + +\begin{thebibliography}{1} + +\bibitem{livrevsrrfd} +J.~M. Mart{\'{\i}} and E.~M{\"u}ller. +\newblock Numerical hydrodynamics in {S}pecial {R}elativity. +\newblock {\em Living Rev. Relativity}, {\bf 3}, 1999. +\newblock [Article in online journal], cited on 31/7/01, + http://www.livingreviews.org/Articles/Volume2/1999-3marti/index.html. + +\bibitem{hydro1} +J.~A. Font, M. Miller, W. Suen and M. Tobias. +\newblock Three Dimensional Numerical General Relativistic +Hydrodynamics I: Formulations, Methods, and Code Tests +\newblock {\em Phys. Rev.}, {\bf D61}, 044011, 2000. + +\end{thebibliography} + + +\include{interface} +\include{param} +\include{schedule} + +\end{document} + + diff --git a/interface.ccl b/interface.ccl new file mode 100644 index 0000000..48a5456 --- /dev/null +++ b/interface.ccl @@ -0,0 +1,84 @@ +# Interface definition for thorn Whisky_Init_Data +# $Header$ + +implements: whisky_init_data +inherits: whisky grid + +USES INCLUDE: EOS_Base.inc +USES INCLUDE: SpaceMask.h +USES INCLUDE: EOS_Base.h + + +SUBROUTINE SpatialDet(CCTK_REAL IN gxx, CCTK_REAL IN gxy, \ + CCTK_REAL IN gxz, CCTK_REAL IN gyy, \ + CCTK_REAL IN gyz, CCTK_REAL IN gzz, \ + CCTK_REAL OUT det) + + +SUBROUTINE Prim2ConPoly(CCTK_INT IN handle, \ + CCTK_REAL IN gxx, CCTK_REAL IN gxy, CCTK_REAL IN gxz, \ + CCTK_REAL IN gyy, CCTK_REAL IN gyz, CCTK_REAL IN gzz, \ + CCTK_REAL IN det, CCTK_REAL OUT dens, \ + CCTK_REAL OUT sx, CCTK_REAL OUT sy, \ + CCTK_REAL OUT sz, CCTK_REAL OUT tau, \ + CCTK_REAL IN rho, CCTK_REAL IN velx, \ + CCTK_REAL IN vely, \ + CCTK_REAL IN velz, CCTK_REAL OUT epsilon, \ + CCTK_REAL OUT press, CCTK_REAL OUT w_lorentz) + + +SUBROUTINE Prim2ConGen(CCTK_INT IN handle, \ + CCTK_REAL IN gxx, CCTK_REAL IN gxy, \ + CCTK_REAL IN gxz, CCTK_REAL IN gyy, \ + CCTK_REAL IN gyz, CCTK_REAL IN gzz, \ + CCTK_REAL IN det, CCTK_REAL OUT dens, \ + CCTK_REAL OUT sx, CCTK_REAL OUT sy, \ + CCTK_REAL OUT sz, CCTK_REAL OUT tau, \ + CCTK_REAL IN rho, CCTK_REAL IN velx, \ + CCTK_REAL IN vely, \ + CCTK_REAL IN velz, CCTK_REAL IN epsilon, \ + CCTK_REAL OUT press, CCTK_REAL OUT w_lorentz) + +SUBROUTINE Con2PrimPoly(CCTK_INT IN handle, CCTK_REAL OUT dens, \ + CCTK_REAL OUT sx, CCTK_REAL OUT sy, \ + CCTK_REAL OUT sz, CCTK_REAL OUT tau, \ + CCTK_REAL OUT rho, CCTK_REAL OUT velx, \ + CCTK_REAL OUT vely, CCTK_REAL OUT velz, \ + CCTK_REAL OUT epsilon, CCTK_REAL OUT press, \ + CCTK_REAL OUT w_lorentz, CCTK_REAL IN uxx, \ + CCTK_REAL IN uxy, CCTK_REAL IN uxz, CCTK_REAL IN uyy, \ + CCTK_REAL IN uyz, CCTK_REAL IN uzz, CCTK_REAL IN det, \ + CCTK_REAL IN x, CCTK_REAL IN y, CCTK_REAL IN z, \ + CCTK_REAL IN r, CCTK_REAL IN rho_min, \ + CCTK_INT IN whisky_reflevel, CCTK_REAL IN whisky_C2P_failed) + + +USES FUNCTION SpatialDet +USES FUNCTION Prim2ConPoly +USES FUNCTION Prim2ConGen +USES FUNCTION Con2PrimPoly + +protected: + +CCTK_REAL simple_wave_grid_functions TYPE=GF TIMELEVELS=1 tags='checkpoint="no"' +{ + simple_tmp + c_s +} "1D arrays for the simple-wave routine" + +CCTK_REAL simple_wave_scalars TYPE=scalar +{ + simple_rho_0 + simple_eps_0 +} "values at v=0" + +CCTK_REAL simple_wave_output TYPE=GF TIMELEVELS=1 tags='checkpoint="no"' +{ + simple_rho + simple_eps +# simple_entropy +} "output variables for the simple-wave routine" + +private: + +int whisky_init_data_reflevel type = SCALAR tags='checkpoint="no"' "Refinement level Whisky is working on right now" diff --git a/param.ccl b/param.ccl new file mode 100644 index 0000000..e782674 --- /dev/null +++ b/param.ccl @@ -0,0 +1,103 @@ +# Parameter definitions for thorn Whisky_Init_Data +# $Header$ + +shares:HydroBase + +USES CCTK_INT timelevels + +shares:ADMBase + +EXTENDS KEYWORD initial_data "" +{ +# "shocktube" :: "Shock tube initial data for Whisky" + "con2primtest" :: "Testing the con -> prim conversion" + "con2prim2con_test" :: "Testing the con -> prim -> con conversion" + "reconstruction_test" :: "Testing reconstruction" +} + +private: + +KEYWORD whisky_initial_data "Initial data for hydro variables only" +{ + "none" :: "Nothing set here" + "shocktube" :: "Shocktube type" + "only_atmo" :: "Set only a low atmosphere" + "read_conformal":: "After reading in initial alp, rho and gxx from h5 files, sets the other quantities" + "simple_wave" :: "Set initial data from Anile Miller Motta, Phys.Fluids. 26, 1450 (1983)" +} "none" + +KEYWORD shocktube_type "Diagonal or parallel shock?" +{ + "diagshock" :: "Diagonal across all axes" + "xshock" :: "Parallel to x axis" + "yshock" :: "Parallel to y axis" + "zshock" :: "Parallel to z axis" + "sphere" :: "spherically symmetric shock" +} "diagshock" + +KEYWORD shock_case "Simple, Sod's problem or other?" +{ + "Simple" :: "GRAstro_Hydro test case" + "Sod" :: "Sod's problem" + "Blast" :: "Strong blast wave" +} "Sod" + +REAL shock_xpos "Position of shock plane: x" +{ + *:* :: "Anything" +} 0.0 + +REAL shock_ypos "Position of shock plane: y" +{ + *:* :: "Anything" +} 0.0 + +REAL shock_zpos "Position of shock plane: z" +{ + *:* :: "Anything" +} 0.0 + +REAL shock_radius "Radius of sperical shock" +{ + 0.0:* :: "Anything positive" +} 1.0 + +BOOLEAN change_shock_direction "Change the shock direction" +{ +} "no" + +REAL simple_wave_constant_c_0 "The c_0 constant in Anile Miller Motta, Phys.Fluids. 26, 1450 (1983)" +{ + 0:1 :: "It is the sound speed where the fluid velocity is zero" +} 0.3 + +REAL simple_wave_v_max "The v_max constant in Anile Miller Motta, Phys.Fluids. 26, 1450 (1983)" +{ + 0:1 :: "It is the maximum velocity in the initial configuration (see p. 1457, bottom of first column)" +} 0.7 + +# For the "atmosphere" + +REAL atmosphere_vel[3] "Velocity of the atmosphere if non-trivial" +{ + *:* :: "Anything" +} 0.0 + +BOOLEAN attenuate_atmosphere "Attenuate the velocity in the atmosphere" +{ +} "no" + +shares:whisky + +USES real whisky_rho_central +USES real whisky_eps_min "" +USES real whisky_perc_ptol "" +USES real whisky_del_ptol "" +USES string whisky_eos_type "" +USES string whisky_eos_table "" +USES string tvd_limiter "" +USES real rho_abs_min +USES real rho_rel_min +USES REAL initial_rho_abs_min +USES REAL initial_rho_rel_min +USES REAL initial_atmosphere_factor diff --git a/schedule.ccl b/schedule.ccl new file mode 100644 index 0000000..36f1461 --- /dev/null +++ b/schedule.ccl @@ -0,0 +1,80 @@ +# Schedule definitions for thorn Whisky_Init_Data + +schedule Whisky_InitData_CheckParameters AT CCTK_PARAMCHECK +{ + LANG: C +} "Check parameters" + +if (CCTK_Equals(whisky_initial_data,"shocktube")) { + schedule Whisky_shocktube in HydroBase_Initial + { + LANG: Fortran + } "Shocktube initial data" +} + +if (CCTK_Equals(initial_data,"con2primtest")) { + STORAGE:whisky_init_data_reflevel + schedule Whisky_Init_Data_RefinementLevel IN HydroBase_Initial BEFORE Whisky_con2primtest + { + LANG: Fortran + } "Calculate current refinement level" + + schedule Whisky_con2primtest in HydroBase_Initial + { + LANG: Fortran + } "Testing the conservative to primitive solver" +} + +if (CCTK_Equals(initial_data,"con2prim2con_test")) { + STORAGE:whisky_init_data_reflevel + schedule Whisky_Init_Data_RefinementLevel IN HydroBase_Initial BEFORE c2p2c + { + LANG: Fortran + } "Calculate current refinement level" + + schedule c2p2c in HydroBase_Initial + { + LANG: Fortran + } "Testing conservative to primitive to conservative" +} + +if (CCTK_Equals(initial_data,"reconstruction_test")) { + schedule Whisky_reconstruction_test in HydroBase_Initial + { + LANG: Fortran + STORAGE: whisky_prim_bext + STORAGE: whisky_scalars + OPTIONS: global loop-local + } "Testing the reconstruction" +} + +if (CCTK_Equals(whisky_initial_data,"only_atmo")) { + schedule Whisky_Only_Atmo in HydroBase_Initial + { + LANG: Fortran + } "Only atmosphere as initial data" +} + +if (CCTK_Equals(whisky_initial_data,"read_conformal")) { + schedule Whisky_ReadConformalData in HydroBase_Initial + { + LANG: Fortran + } "Set the missing quantities, after reading in from file initial data from conformally-flat codes (Garching)" +} + +if (CCTK_Equals(whisky_initial_data,"simple_wave")) { + STORAGE: simple_wave_grid_functions + STORAGE: simple_wave_scalars + schedule Whisky_SimpleWave in HydroBase_Initial + { + LANG: Fortran + } "Set initial data from Anile Miller Motta, Phys.Fluids. 26, 1450 (1983)" + + STORAGE: simple_wave_output + schedule Whisky_SimpleWave_Analysis AT CCTK_ANALYSIS AFTER Whisky_Entropy + { + LANG: Fortran + } "Compute some output variables for the Simple Wave" + + +} diff --git a/src/CheckParam.c b/src/CheckParam.c new file mode 100644 index 0000000..68f446c --- /dev/null +++ b/src/CheckParam.c @@ -0,0 +1,15 @@ +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" + +void Whisky_InitData_CheckParameters(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + if (timelevels < 2) + { + CCTK_PARAMWARN("You have to set 'HydroBase::timelevels to at least 2"); + } +} + diff --git a/src/Whisky_C2P2C.F90 b/src/Whisky_C2P2C.F90 new file mode 100644 index 0000000..3c55917 --- /dev/null +++ b/src/Whisky_C2P2C.F90 @@ -0,0 +1,131 @@ + /*@@ + @file Whisky_C2P2C.F90 + @date Sat Jan 26 02:44:43 2002 + @author Luca Baiotti + @desc + A test of the conservative <--> primitive variable exchange + @enddesc + @@*/ + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" + + /*@@ + @routine c2p2c + @date Sat Jan 26 02:45:19 2002 + @author Luca Baiotti + @desc + Testing the conservative <--> primitive variable transformations. + The values before and after should match. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +subroutine c2p2c(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + integer didit,i,j,k,nx,ny,nz + CCTK_REAL det,uxx,uxy,uxz,uyy,uyz,uzz + CCTK_REAL dens_send,sx_send,sy_send,sz_send,tau_send + CCTK_REAL rho_send,velx_send,vely_send,velz_send,eps_send + CCTK_REAL press_send,w_lorentz_send,x_send,y_send,z_send,r_send + logical epsnegative + + call CCTK_WARN(1,"This test works only with Ideal_Fluid EoS") + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + + x_send = 0.0d0 + y_send = 0.0d0 + z_send = 0.0d0 + r_send = 0.0d0 + + gxx = 1.0d0 + gyy = 1.0d0 + gzz = 1.0d0 + gxy = 0.0d0 + gxz = 0.0d0 + gyz = 0.0d0 + + det = 1.0d0 + + uxx = 1.0d0 + uyy = 1.0d0 + uzz = 1.0d0 + uxy = 0.0d0 + uxz = 0.0d0 + uyz = 0.0d0 + + dens_send = 1.29047362d0 + sx_send = 0.166666658d0 + sy_send = 0.166666658d0 + sz_send = 0.166666658d0 + tau_send = 0.484123939d0 + + eps_send = 1.0d-6 + press_send = 6.666666666666667d-7 + w_lorentz_send = 1.0d0 + + epsnegative = .false. + + whisky_rho_min = 1.e-10 + + write(*,*) 'C2P2C test: initial values.' + write(*,*) ' conservative variables: ' + write(*,*) ' dens: ',dens_send + write(*,*) ' sx : ',sx_send + write(*,*) ' sy : ',sy_send + write(*,*) ' sz : ',sz_send + write(*,*) ' tau : ',tau_send + write(*,*) ' eps : ',eps_send + write(*,*) ' W : ',w_lorentz_send + + write(*,*) 'C2P2C test: getting the associated primitive variables.' + call Con2Prim_pt(whisky_eos_handle,dens_send,sx_send,sy_send,sz_send, & + tau_send,rho_send,velx_send,vely_send,velz_send, & + eps_send,press_send,w_lorentz_send, & + uxx,uxy,uxz,uyy,uyz,uzz,det,x_send,y_send,z_send,r_send,& + epsnegative,whisky_rho_min, whisky_init_data_reflevel,whisky_C2P_failed) + + write(*,*) 'C2P2C test: the primitive variables are' + write(*,*) ' primitive variables: ' + write(*,*) ' rho : ',rho_send + write(*,*) ' velx : ',velx_send + write(*,*) ' vely : ',vely_send + write(*,*) ' velz : ',velz_send + write(*,*) ' press : ',press_send + write(*,*) ' eps : ',eps_send + write(*,*) ' W : ',w_lorentz_send + + write(*,*) 'C2P2C test: converting back to conserved variables.' + call prim2con(whisky_eos_handle,gxx, gxy, gxz, gyy, gyz, gzz, det, & + dens_send, sx_send, sy_send, sz_send, tau_send, rho_send, & + velx_send, vely_send, velz_send, eps_send, press_send, w_lorentz_send) + + write(*,*) 'C2P2C test: the conserved variables are' + write(*,*) ' conservative variables: ' + write(*,*) ' dens: ',dens_send + write(*,*) ' sx : ',sx_send + write(*,*) ' sy : ',sy_send + write(*,*) ' sz : ',sz_send + write(*,*) ' tau : ',tau_send + write(*,*) ' eps : ',eps_send + write(*,*) ' W : ',w_lorentz_send + + STOP + + return + +end subroutine c2p2c diff --git a/src/Whisky_Con2Prim.F90 b/src/Whisky_Con2Prim.F90 new file mode 100644 index 0000000..37c2916 --- /dev/null +++ b/src/Whisky_Con2Prim.F90 @@ -0,0 +1,108 @@ + /*@@ + @file Whisky_Con2Prim.F90 + @date Sat Jan 26 02:49:32 2002 + @author Luca Baiotti + @desc + A test of the conservative to primitive exchange. + @enddesc + @@*/ + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" + + /*@@ + @routine Whisky_con2primtest + @date Sat Jan 26 02:49:58 2002 + @author Luca Baiotti + @desc + A test of the conservative to primitive variable solver. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +subroutine Whisky_Init_Data_RefinementLevel(CCTK_ARGUMENTS) + + implicit none + DECLARE_CCTK_ARGUMENTS + + whisky_init_data_reflevel = aint(log10(dble(cctk_levfac(1)))/log10(2.0d0)) + +end subroutine Whisky_Init_Data_RefinementLevel + + +subroutine Whisky_con2primtest(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + integer didit,i,j,k,nx,ny,nz + CCTK_REAL det,uxx,uxy,uxz,uyy,uyz,uzz + CCTK_REAL dens_send,sx_send,sy_send,sz_send,tau_send + CCTK_REAL rho_send,velx_send,vely_send,velz_send,eps_send + CCTK_REAL press_send,w_lorentz_send,x_send,y_send,z_send,r_send + logical epsnegative + + call CCTK_WARN(1,"For this test, remember to use a polytropic EoS and to set eos_gamma = 2.0 and eos_k = 100.0") + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + + x_send = 0.0d0 + y_send = 0.0d0 + z_send = 0.0d0 + r_send = 0.0d0 + + det = 1.0d0 + uxx = 1.0d0 + uyy = 1.0d0 + uzz = 1.0d0 + uxy = 0.0d0 + uxz = 0.0d0 + uyz = 0.0d0 + + dens_send = 1.29047362d0 + sx_send = 0.166666658d0 + sy_send = 0.166666658d0 + sz_send = 0.166666658d0 + tau_send = 0.484123939d0 + + rho_send = 1.0d0 + velx_send = 0.0d0 + vely_send = 0.0d0 + velz_send = 0.0d0 + eps_send = 1.0d-6 + press_send = 6.666666666666667d-7 + w_lorentz_send = 1.0d0 + + epsnegative = .false. + + write(*,*) 'Con2Prim test: converting to primitive variables' + call Con2Prim_pt(whisky_eos_handle,dens_send,sx_send,sy_send,sz_send,& + tau_send,rho_send,velx_send,vely_send,velz_send,& + eps_send,press_send,w_lorentz_send, & + uxx,uxy,uxz,uyy,uyz,uzz,det,x_send,y_send,z_send,r_send,& + epsnegative,whisky_rho_min, whisky_init_data_reflevel, whisky_C2P_failed) + write(*,*) 'Con2Prim test: the primitive variables are' + write(*,*) ' primitive variables: ' + write(*,*) ' rho : ',rho_send + write(*,*) ' velx : ',velx_send + write(*,*) ' vely : ',vely_send + write(*,*) ' velz : ',velz_send + write(*,*) ' eps : ',eps_send + write(*,*) ' press : ',press_send + write(*,*) ' w_lor : ',w_lorentz_send + + STOP + + return + +end subroutine Whisky_con2primtest diff --git a/src/Whisky_Only_Atmo.F90 b/src/Whisky_Only_Atmo.F90 new file mode 100644 index 0000000..7bdc79f --- /dev/null +++ b/src/Whisky_Only_Atmo.F90 @@ -0,0 +1,114 @@ + /*@@ + @file Whisky_Only_Atmo.F90 + @date Sat Jan 29 04:44:25 2002 + @author Luca Baiotti + @desc + Hydro initial data for vacuum. + @enddesc + @@*/ + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" + +#define velx(i,j,k) vel(i,j,k,1) +#define vely(i,j,k) vel(i,j,k,2) +#define velz(i,j,k) vel(i,j,k,3) +#define sx(i,j,k) scon(i,j,k,1) +#define sy(i,j,k) scon(i,j,k,2) +#define sz(i,j,k) scon(i,j,k,3) + + /*@@ + @routine Whisky_Only_Atmo + @date Sat Jan 29 04:53:43 2002 + @author Luca Baiotti + @desc + Hydro initial data for vacuum. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +subroutine Whisky_Only_Atmo(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + CCTK_INT i,j,k,nx,ny,nz + CCTK_REAL det, atmo + + call CCTK_INFO("Setting only atmosphere as initial data.") + + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + + if (rho_abs_min > 0.0) then + atmo = rho_abs_min + whisky_rho_min = rho_abs_min + else + call CCTK_WARN(0,"For initial data 'only_atmo' the parameter 'rho_abs_min' must be set") + end if + + if (initial_rho_abs_min > 0.0) then + atmo = initial_rho_abs_min + else if (initial_rho_rel_min > 0.0) then + call CCTK_WARN(3,"For initial data 'only_atmo' the parameter 'initial_rho_rel_min' is ignored") + end if + + if (initial_atmosphere_factor > 0.0) then + atmo = atmo * initial_atmosphere_factor + end if + + rho = atmo + if (rho_abs_min < initial_rho_abs_min) then + velx(:,:,:) = atmosphere_vel(1) + vely(:,:,:) = atmosphere_vel(2) + velz(:,:,:) = atmosphere_vel(3) + else + velx(:,:,:) = 0.d0 + vely(:,:,:) = 0.d0 + velz(:,:,:) = 0.d0 + end if + if (attenuate_atmosphere .ne. 0) then + velx(:,:,:) = velx(:,:,:) * (1.d0 - exp(-2.d0 * (r - 1.d0)**2)) + vely(:,:,:) = vely(:,:,:) * (1.d0 - exp(-2.d0 * (r - 1.d0)**2)) + velz(:,:,:) = velz(:,:,:) * (1.d0 - exp(-2.d0 * (r - 1.d0)**2)) + end if + eps = whisky_eps_min + + do i=1,nx + do j=1,ny + do k=1,nz + +!!$ if (x(i,j,k) < 0.0) then +!!$ velx(i,j,k) = -velx(i,j,k) +!!$ vely(i,j,k) = -vely(i,j,k) +!!$ velz(i,j,k) = -velz(i,j,k) +!!$ end if + + call SpatialDeterminant(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),& + gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),det) + call prim2con(whisky_eos_handle,gxx(i,j,k),gxy(i,j,k),& + gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),& + det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),& + eps(i,j,k),press(i,j,k),w_lorentz(i,j,k)) + enddo + enddo + enddo + + densrhs = 0.d0 + srhs = 0.d0 + taurhs = 0.d0 + + return + +end subroutine Whisky_Only_Atmo diff --git a/src/Whisky_ReadConformalData.F90 b/src/Whisky_ReadConformalData.F90 new file mode 100644 index 0000000..0d31976 --- /dev/null +++ b/src/Whisky_ReadConformalData.F90 @@ -0,0 +1,189 @@ + /*@@ + @file Whisky_ReadConformalData.F90 + @date Fri Jul 6 12:29:27 2007 + @author Luca Baiotti + @desc + Set the missing quantities, after reading in from file initial data from conformally-flat codes (Garching) + @enddesc + @@*/ + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#include "SpaceMask.h" + +#define velx(i,j,k) vel(i,j,k,1) +#define vely(i,j,k) vel(i,j,k,2) +#define velz(i,j,k) vel(i,j,k,3) +#define sx(i,j,k) scon(i,j,k,1) +#define sy(i,j,k) scon(i,j,k,2) +#define sz(i,j,k) scon(i,j,k,3) + + /*@@ + @routine Whisky_ReadConformalData + @date Fri Jul 6 12:31:18 2007 + @author Luca Baiotti + @desc + Set the missing quantities, after reading in from file initial data from conformally-flat codes (Garching) + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +subroutine Whisky_ReadConformalData(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + +!#include "EOS_Base.h" +#ifdef _EOS_BASE_INC_ +#undef _EOS_BASE_INC_ +#endif +#include "EOS_Base.inc" + + + CCTK_INT :: i,j,k,handle,ierr + CCTK_REAL :: eos_k, eos_gamma, rho_min, det + + + + ! only gxx has been read in; copy it into gyy and gzz as well + gyy = gxx + gzz = gxx + + ! set the other components to zero + gxy = 0.0 + gxz = 0.0 + gyz = 0.0 + + ! set the extrinsic curvature to zero + kxx = 0.0 + kxy = 0.0 + kxz = 0.0 + + kyy = 0.0 + kyz = 0.0 + kzz = 0.0 + + ! set the shift to zero (the lapse is read in) + betax = 0.0 + betay = 0.0 + betaz = 0.0 + + + ! atmosphere + + if (rho_abs_min > 0.0) then + rho_min = rho_abs_min + else + call CCTK_WARN(0,"For now, with ReadConformalData you need to set rho_abs_min") + end if + + if (initial_rho_abs_min > 0.0) then + rho_min = initial_rho_abs_min; + else if (initial_rho_rel_min > 0.0) then + call CCTK_WARN(1,"Setting initial_rho_rel_min is ignored with ReadConformalData, for the time being") + end if + + if (initial_atmosphere_factor > 0.0) rho_min = rho_min * initial_atmosphere_factor + + ierr = 0 + do i=1,cctk_lsh(1) + do j=1,cctk_lsh(2) + do k=1,cctk_lsh(3) + + ! write(*,*) i,j,k, rho(i,j,k) + ! check on the read-in rho + if (rho(i,j,k) > 0.0) then + write(*,*) i,j,k, rho(i,j,k) + ierr = ierr+1 + end if + + ! set the atmosphere + if (rho(i,j,k) <= 0.d0) rho(i,j,k) = rho_min + +! if (rho(i,j,k) > 1.d-9) write(*,*) i,j,k,rho(i,j,k) + end do + end do + end do + +! if (ierr == 0) call CCTK_WARN(0,"rho.h5 contains only zeroes: stopping the simulation") + +! write(*,*)"rho min", rho_min + + +!!$ where(rho < rho_min) +!!$ rho = rho_min +!!$ end where + + + ! set pressure and eps from rho, using the polytropic EoS + +! handle = EOS_Handle("2D_Polytrope") + +! if (handle < 0) call CCTK_WARN(0,"For this hack you need to compile with EOS_Polytrope") + + eos_k = EOS_Pressure(whisky_eos_handle, 1.0, 1.0) + eos_gamma = 1.0 + eos_k / EOS_SpecificIntEnergy(whisky_eos_handle,1.0,1.0) + + +! press = eos_k * rho**eos_gamma + + do i=1,cctk_lsh(1) + do j=1,cctk_lsh(2) + do k=1,cctk_lsh(3) + eps(i,j,k) = EOS_SpecificIntEnergy(whisky_eos_handle,rho(i,j,k),press(i,j,k)) + end do + end do + end do + + + ! set velocities to zero + + velx(:,:,:) = 0.0 + vely(:,:,:) = 0.0 + velz(:,:,:) = 0.0 + w_lorentz = 1.0 + + + ! Compute conserved variables + + do i=1,cctk_lsh(1) + do j=1,cctk_lsh(2) + do k=1,cctk_lsh(3) + + call SpatialDet(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),& + gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),det) + + if (CCTK_EQUALS(whisky_eos_type,"Polytype")) then + call Prim2ConPoly(whisky_eos_handle,gxx(i,j,k),gxy(i,j,k),& + gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),& + det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + tau(i,j,k),rho(i,j,k),& + velx(i,j,k),vely(i,j,k),velz(i,j,k),& + eps(i,j,k),press(i,j,k),w_lorentz(i,j,k)) + else + call Prim2ConGen(whisky_eos_handle,gxx(i,j,k),gxy(i,j,k),& + gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),& + det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + tau(i,j,k),rho(i,j,k),& + velx(i,j,k),vely(i,j,k),velz(i,j,k),& + eps(i,j,k),press(i,j,k),w_lorentz(i,j,k)) + end if + + end do + end do + end do + + + return + +end subroutine Whisky_ReadConformalData diff --git a/src/Whisky_ReconstructTest.F90 b/src/Whisky_ReconstructTest.F90 new file mode 100644 index 0000000..d19f786 --- /dev/null +++ b/src/Whisky_ReconstructTest.F90 @@ -0,0 +1,84 @@ + /*@@ + @file Whisky_ReconstructTest.F90 + @date Sat Jan 26 02:51:49 2002 + @author Luca Baiotti + @desc + A test of the reconstruction algorithm. + @enddesc + @@*/ + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" + + /*@@ + @routine Whisky_reconstruction_test + @date Sat Jan 26 02:52:17 2002 + @author Luca Baiotti + @desc + Tests the reconstruction method by giving smoothly varying initial + data with a shock along the x axis. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +subroutine Whisky_reconstruction_test(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + integer i, j, k, nx, ny, nz + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + + xoffset = 1 + yoffset = 0 + zoffset = 0 + + do k=1, cctk_lsh(3) + do j=1, cctk_lsh(2) + do i=1, cctk_lsh(1) + if (i < 8) then + rho(i,j,k) = 5+sin(real(i+j+k)/2.) + else + rho(i,j,k) = 1+sin(real(i+j+k)/2.) + end if + end do + end do + end do + + call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & + rho, rhoplus, rhominus) + + open(unit=11,file='original.dat',status='unknown') + do i=1, CCTK_LSH(1) + write(11,*) i, rho(i,4,7) + write(*,*) i, rho(i,4,7) + end do + close(11) + + open(unit=12,file='extensions.dat',status='unknown') + do i=2, CCTK_LSH(1)-1 + write(12,*) i-0.5, rhominus(i,4,7) + write(12,*) i+0.5, rhoplus(i,4,7) + end do + close(12) + + + write(*,*) 'test_reconstruction: done tvdreconstruct.' + write(*,*) 'data written in files original.dat and extension.dat' + + stop + + return + +end subroutine Whisky_reconstruction_test diff --git a/src/Whisky_ShockTube.F90 b/src/Whisky_ShockTube.F90 new file mode 100644 index 0000000..8d92dbb --- /dev/null +++ b/src/Whisky_ShockTube.F90 @@ -0,0 +1,155 @@ + /*@@ + @file Whisky_ShockTube.F90 + @date Sat Jan 26 02:53:25 2002 + @author Ian Hawke + @desc + Initial data of the shock tube type. + @enddesc + @@*/ + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" + +#define velx(i,j,k) vel(i,j,k,1) +#define vely(i,j,k) vel(i,j,k,2) +#define velz(i,j,k) vel(i,j,k,3) +#define sx(i,j,k) scon(i,j,k,1) +#define sy(i,j,k) scon(i,j,k,2) +#define sz(i,j,k) scon(i,j,k,3) + + /*@@ + @routine Whisky_shocktube + @date Sat Jan 26 02:53:49 2002 + @author Ian Hawke + @desc + Initial data for shock tubes. Either diagonal or parallel to + a coordinate axis. Either Sods problem or the standard shock tube. + @enddesc + @calls + @calledby + @history + Expansion and alteration of the test code from GRAstro_Hydro, + written by Mark Miller. + @endhistory + +@@*/ + +subroutine Whisky_shocktube(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + CCTK_INT :: i, j, k, nx, ny, nz + CCTK_REAL :: direction, det + CCTK_REAL :: rhol, rhor, velxl, velxr, velyl, velyr, & + velzl, velzr, epsl, epsr + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + + do i=1,nx + do j=1,ny + do k=1,nz + if (CCTK_EQUALS(shocktube_type,"diagshock")) then + direction = x(i,j,k) - shock_xpos + & + y(i,j,k) - shock_ypos + z(i,j,k) - shock_zpos + else if (CCTK_EQUALS(shocktube_type,"xshock")) then + direction = x(i,j,k) - shock_xpos + else if (CCTK_EQUALS(shocktube_type,"yshock")) then + direction = y(i,j,k) - shock_ypos + else if (CCTK_EQUALS(shocktube_type,"zshock")) then + direction = z(i,j,k) - shock_zpos + else if (CCTK_EQUALS(shocktube_type,"sphere")) then + direction = sqrt((x(i,j,k)-shock_xpos)**2+& + (y(i,j,k)-shock_ypos)**2+& + (z(i,j,k)-shock_zpos)**2)-shock_radius + end if + if (CCTK_EQUALS(shock_case,"Simple")) then + rhol = 10.d0 + rhor = 1.d0 + velxl = 0.d0 + velxr = 0.d0 + velyl = 0.d0 + velyr = 0.d0 + velzl = 0.d0 + velzr = 0.d0 + epsl = 2.d0 + epsr = 1.d-6 + else if (CCTK_EQUALS(shock_case,"Sod")) then + rhol = 1.d0 + rhor = 0.125d0 + velxl = 0.d0 + velxr = 0.d0 + velyl = 0.d0 + velyr = 0.d0 + velzl = 0.d0 + velzr = 0.d0 + epsl = 1.5d0 + epsr = 1.2d0 +!!$This line only for polytrope, k=1 +!!$ epsr = 0.375d0 + else if (CCTK_EQUALS(shock_case,"Blast")) then + rhol = 1.d0 + rhor = 1.d0 + velxl = 0.d0 + velxr = 0.d0 + velyl = 0.d0 + velyr = 0.d0 + velzl = 0.d0 + velzr = 0.d0 + epsl = 1500.d0 + epsr = 1.5d-2 + else + call CCTK_WARN(0,"Shock case not recognized") + end if + + if ( ((change_shock_direction==0).and.(direction .lt. 0.0d0)).or.& + ((change_shock_direction==1).and.(direction .gt. 0.0d0)) ) then + rho(i,j,k) = rhol + velx(i,j,k) = velxl + vely(i,j,k) = velyl + velz(i,j,k) = velzl + eps(i,j,k) = epsl + else + rho(i,j,k) = rhor + velx(i,j,k) = velxr + vely(i,j,k) = velyr + velz(i,j,k) = velzr + eps(i,j,k) = epsr + end if + + call SpatialDet(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),& + gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),det) + + if (CCTK_EQUALS(whisky_eos_type,"Polytype")) then + call Prim2ConPoly(whisky_eos_handle,gxx(i,j,k),gxy(i,j,k),& + gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),& + det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + tau(i,j,k),rho(i,j,k),& + velx(i,j,k),vely(i,j,k),velz(i,j,k),& + eps(i,j,k),press(i,j,k),w_lorentz(i,j,k)) + else + call Prim2ConGen(whisky_eos_handle,gxx(i,j,k),gxy(i,j,k),& + gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),& + det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + tau(i,j,k),rho(i,j,k),& + velx(i,j,k),vely(i,j,k),velz(i,j,k),& + eps(i,j,k),press(i,j,k),w_lorentz(i,j,k)) + end if + enddo + enddo + enddo + + densrhs = 0.d0 + srhs = 0.d0 + taurhs = 0.d0 + + return + +end subroutine Whisky_shocktube diff --git a/src/Whisky_SimpleWave.F90 b/src/Whisky_SimpleWave.F90 new file mode 100644 index 0000000..b840822 --- /dev/null +++ b/src/Whisky_SimpleWave.F90 @@ -0,0 +1,238 @@ + + /*@@ + @file Whisky_SimpleWave.F90 + @date Thu Aug 2 15:17:35 2007 + @author Luca Baiotti + @desc + Initial data for a simple wave with sinusoidal initial function for the velocity + See Anile, Miller, Motta, Formation and damping of relativistic strong shocks, + Phys. Fluids 26, 1450 (1983) + @enddesc + @@*/ + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" + +#define velx(i,j,k) vel(i,j,k,1) +#define vely(i,j,k) vel(i,j,k,2) +#define velz(i,j,k) vel(i,j,k,3) +#define sx(i,j,k) scon(i,j,k,1) +#define sy(i,j,k) scon(i,j,k,2) +#define sz(i,j,k) scon(i,j,k,3) + + /*@@ + @routine Whisky_SimpleWave + @date Thu Aug 2 15:20:28 2007 + @author Luca Baiotti + @desc + Initial data for a simple wave with sinusoidal initial function for the velocity + See Anile, Miller, Motta, Formation and damping of relativistic strong shocks, + Phys. Fluids 26, 1450 (1983) + @enddesc + @calls + @calledby + @history + @endhistory + +@@*/ + +subroutine Whisky_SimpleWave(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + +#ifdef _EOS_BASE_INC_ +#undef _EOS_BASE_INC_ +#endif +#include "EOS_Base.inc" + + CCTK_INT :: i, j, k, nx, ny, nz + CCTK_REAL :: dr, k1, k2, k3, k4, in_data, old_data, source_data, new_data, c_0, det, pi + + call CCTK_INFO("Setting initial data for a simple wave as Anile Miller Motta") + + call CCTK_WARN(1, "The simple-wave initial-data routine works only for unigrid and on one node.") + + pi = 4.d0 * atan(1.0) + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + + velx(:,:,:) = 0.0 + rho = 0.0 + eps = 0.0 + press = 0.0 + w_lorentz = 0.0 + simple_tmp = 0.0 + + i = 0 + do while (minval(velx(:,1,1)) >= 0.0) + i=i+1 + if (x(i,1,1) > 0.0) velx(i,1,1) = simple_wave_v_max * sin(pi * x(i,1,1)) + end do + + velx(i,1,1) = 0.0 ! set the first term that became negative to zero + + c_0 = simple_wave_constant_c_0 ! a parameter + + ! compute quantities at v=0 + + simple_eps_0 = 9.d0*c_0**2 /(4.d0 * (1.d0 - 3.d0 * c_0**2)) + + simple_rho_0 = simple_eps_0**3.d0 + +! press_0 = eps_0**4.d0 + +! do j = 1,i + do j = 1,CCTK_LSH(1) + + simple_tmp(j,1,1) = ( 1.d0 + c_0*sqrt(3.d0) )/( 1.d0 - c_0*sqrt(3.d0) ) & + * ( (1+velx(j,1,1))/(1-velx(j,1,1)) )**(1/(2.d0*sqrt(3.d0))) + + c_s(j,1,1) = ( simple_tmp(j,1,1) - 1.d0 ) / ( sqrt(3.d0)* ( simple_tmp(j,1,1) + 1.d0 ) ) + + eps(j,1,1) = 9.d0*c_s(j,1,1)**2 /(4.d0 * (1.d0 - 3.d0 * c_s(j,1,1)**2)) + + rho(j,1,1) = eps(j,1,1)**3.d0 +! write(*,*) j, x(j,1,1), rho(j,1,1) +! rho(j,1,1) = rho_abs_min * rho(j,1,1) + + press(j,1,1) = eps(j,1,1)**(4.d0) + + w_lorentz(j,1,1) = 1.d0/sqrt(1.d0-velx(j,1,1)**2) ! flat spacetime + + end do + + !arrays + gxx = 1.d0; gyy=1.d0; gzz=1.d0; gxy=0.d0; gxz=0.d0; gyz=0.d0 + + + +!!$ do i = 1,CCTK_LSH(1) +!!$ +!!$ write(*,*) i, x(i,1,1), velx(i,1,1), rho(i,1,1) +!!$ write(*,*) eps(i,1,1), press(i,1,1), w_lorentz(i,1,1) +!!$ +!!$ end do + + + + do i=1,nx + + ! atmosphere + + if ( (rho(i,1,1) < whisky_rho_min).OR.(velx(i,1,1) < 0) ) then + rho(i,1,1) = rho_abs_min +! rho(i,1,1) = 1.0 !the value of rho_min for the initial data + eps(i,1,1) = rho_abs_min**(1.d0/3.d0) + velx(i,1,1) = 0.d0 + w_lorentz(i,1,1) = 1.d0 + press(i,1,1) = EOS_Pressure(whisky_polytrope_handle, rho(i,1,1), 1.0d0) + ! polytrope only (initial data) + end if + +! write(*,*) 'p',i, x(i,1,1), rho(i,1,1)**(4.d0/3.d0)/press(i,1,1) + + call SpatialDet(gxx(i,1,1),gxy(i,1,1),gxz(i,1,1),& + gyy(i,1,1),gyz(i,1,1),gzz(i,1,1),det) + +! if (CCTK_EQUALS(whisky_eos_type,"Polytype")) then + ! always use polytype for initial data + call Prim2ConPoly(whisky_polytrope_handle,gxx(i,1,1),gxy(i,1,1),& + gxz(i,1,1),gyy(i,1,1),gyz(i,1,1),gzz(i,1,1),& + det, dens(i,1,1),sx(i,1,1),sy(i,1,1),sz(i,1,1),& + tau(i,1,1),rho(i,1,1),& + velx(i,1,1),vely(i,1,1),velz(i,1,1),& + eps(i,1,1),press(i,1,1),w_lorentz(i,1,1)) +!!$ else +!!$ call Prim2ConGen(whisky_eos_handle,gxx(i,1,1),gxy(i,1,1),& +!!$ gxz(i,1,1),gyy(i,1,1),gyz(i,1,1),gzz(i,1,1),& +!!$ det, dens(i,1,1),sx(i,1,1),sy(i,1,1),sz(i,1,1),& +!!$ tau(i,1,1),rho(i,1,1),& +!!$ velx(i,1,1),vely(i,1,1),velz(i,1,1),& +!!$ eps(i,1,1),press(i,1,1),w_lorentz(i,1,1)) +!!$ end if +!!$ + +! write(*,*) 'd',i, x(i,1,1), rho(i,1,1)**(4.d0/3.d0)/press(i,1,1) + + enddo + + + ! planar symmetry + do j=1,ny + do k=1,nz + + velx(:,j,k) = velx(:,1,1) + rho(:,j,k) = rho(:,1,1) + eps(:,j,k) = eps(:,1,1) + press(:,j,k) = press(:,1,1) + w_lorentz(:,j,k) = w_lorentz(:,1,1) + + sx(:,j,k) = sx(:,1,1) + dens(:,j,k) = dens(:,1,1) + tau(:,j,k) = tau(:,1,1) + + enddo + enddo + + + vely(:,:,:) = 0.d0 + velz(:,:,:) = 0.d0 + + densrhs = 0.d0 + srhs = 0.d0 + taurhs = 0.d0 + + simple_tmp = rho + + + +!!$ do i = 1,CCTK_LSH(1) +!!$ +!!$ do j = 1,CCTK_LSH(2) +!!$ +!!$ do k = 1,CCTK_LSH(3) +!!$ +!!$ write(*,*) i, x(i,j,k), velx(i,j,k) +!!$ write(*,*) eps(i,j,k), press(i,j,k), rho(i,j,k) +!!$ +!!$ end do +!!$ end do +!!$ end do + + +call CCTK_INFO("Finished initial data") + + return + +end subroutine Whisky_SimpleWave + + +subroutine Whisky_SimpleWave_Analysis(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + if (velx(CCTK_LSH(1),1,1) == 0.0) then + simple_eps_0 = eps(CCTK_LSH(1),1,1) + else + call CCTK_WARN(1,"The wave has reached the outer boundary: the computation of simple_eps is now wrong") + end if + + simple_rho = (rho - simple_rho_0)/simple_rho_0 + simple_eps = (eps - simple_eps_0)/simple_eps_0 + + + return + +end subroutine Whisky_SimpleWave_Analysis diff --git a/src/Whisky_TOV.F90 b/src/Whisky_TOV.F90 new file mode 100644 index 0000000..c62044d --- /dev/null +++ b/src/Whisky_TOV.F90 @@ -0,0 +1,324 @@ + /*@@ + @file Whisky_TOV.F90 + @date Wed Feb 13 02:53:25 2002 + @author Scott Hawley + @desc + Initial data of the TOV type. + @enddesc + @@*/ + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" + + /*@@ + @routine Whisky_TOV + @date Sat Jan 26 02:53:49 2002 + @author Scott Hawley + @desc + Initial data for TOV stars. + @enddesc + @calls + @calledby + @history + Alteration of the Whisky_shockwave written by Ian Hawke + @endhistory + +@@*/ + +subroutine Whisky_TOV(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + real*8 rho0 + common / com_tov / rho0 + + integer nr, n_succ + integer i,j,k,nx,ny,nz + integer CCTK_Equals + CCTK_REAL rmax, dr + CCTK_REAL xlb, xub + CCTK_REAL ylb, yub + CCTK_REAL zlb, zub + CCTK_REAL tol + parameter ( tol = 1e-10 ) + + CCTK_REAL, allocatable, dimension (:) :: gtt, grr + + CCTK_REAL det + integer rc + +#include "EOS_Base.inc" + +#ifdef HAVE_ODEPACK + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + +!----------------------------------------------------------- +! Set up the 1D spherical grid +!----------------------------------------------------------- +! decide maximum radius to integrate out to: +! Just use the maximum diagonal width of the grid +! BUG: this may not be the proper procedure call + call CCTK_CoordRange(cctkGH, xlb, xub, -1, "x", "cart3d") + call CCTK_CoordRange(cctkGH, ylb, yub, -1, "y", "cart3d") + call CCTK_CoordRange(cctkGH, zlb, zub, -1, "z", "cart3d") + rmax = sqrt( (xub-xlb)**2 + (yub-ylb)**2 + (zub-zlb)**2 ) + + nr = tov_nr + dr = rmax / (nr - 1) + do j = 1, nr + tov_r(j) = (j-1) * dr + enddo + + +!----------------------------------------------------------- +! Assign data at r=0 +!----------------------------------------------------------- + rho0 = whisky_rho_central + tov_rho(1) = rho0 + tov_alpha(1) = 1.0 + tov_m(1) = 0.0 + +!----------------------------------------------------------- +! Get the 1D TOV data +!----------------------------------------------------------- + call tov_lsoda (tov_rho, tov_p, tov_m, tov_alpha, & + tov_r, nr, tol , n_succ, rc) + +!----------------------------------------------------------- +! Handle any errors from integration routine +!----------------------------------------------------------- + if (rc .ne. 0) then + call CCTK_WARN(0,"Fatal error in TOV integration") + endif + if (n_succ .lt. nr) then + call CCTK_WARN(0,"TOV integration did not extend as far out as requested") + endif + nr = n_succ + +!----------------------------------------------------------- +! fill in additional physical info +!----------------------------------------------------------- + do j=1,nr +! BUG: the call looks something like this... + tov_rho(j) = EOS_RestMassDens(whisky_eos_handle,tov_p(j)) + end do + +!----------------------------------------------------------- +! Generate 3D data from 1D data +! BUG: These Spin1D functions do not yet exist! +!----------------------------------------------------------- + call Spin1D(tov_rho, tov_r, nr, rho, x, y, z, nx, ny, nz) + call Spin1D(tov_p, tov_r, nr, press, x, y, z, nx, ny, nz) + call Spin1D(tov_rho, tov_r, nr, eps, x, y, z, nx, ny, nz) + + allocate(grr(nr), gtt(nr)) + do i = 1, nr + grr(i) = tov_r(i) / (tov_r(i) - 2*tov_m(i)) + gtt(i) = - tov_alpha(i)**2 + enddo + call Spin1DMetric(tov_m, tov_r, nr, gxx, gxy, gxz, gyy, gyz, gzz, & + x, y, z, nx, ny, nz) + + call Spin1DMetric(grr, gtt, gxx, gxy, gxz, gyy, gyz, gzz, nx, ny, nz, & + x, y, z, r) + deallocate(gtt) + deallocate(grr) + + do i=1,nx + do j=1,ny + do k=1,nz + velx(i,j,k) = 0.0 + vely(i,j,k) = 0.0 + velz(i,j,k) = 0.0 + w_lorentz(i,j,k) = 1.0 + end do + end do + enddo + +!----------------------------------------------------------- +! Set additional variables +!----------------------------------------------------------- +! Conserved Variables + det = 1.0d0 + call prim2con(whisky_eos_handle,gxx, gxy, gxz, gyy, gyz, gzz, det, & + dens, sx, sy, sz, tau, rho, & + velx, vely, velz, eps, press, w_lorentz) + +! Copy data to other time steps + do i=1,nx + do j=1,ny + do k=1,nz + +! Primitive variables + rho_p(i,j,k) = rho(i,j,k) + press_p(i,j,k) = press(i,j,k) + eps_p(i,j,k) = eps(i,j,k) + velx_p(i,j,k) = velx(i,j,k) + vely_p(i,j,k) = vely(i,j,k) + velz_p(i,j,k) = velz(i,j,k) + w_lorentz_p(i,j,k) = w_lorentz(i,j,k) + +! Conserved variables + dens_p(i,j,k) = dens(i,j,k) + tau_p(i,j,k) = tau(i,j,k) + sx_p(i,j,k) = sx(i,j,k) + sy_p(i,j,k) = sy(i,j,k) + sz_p(i,j,k) = sz(i,j,k) + + end do + end do + end do + +! BUG: Why am I doing this? + densrhs = 0.d0 + sxrhs = 0.d0 + syrhs = 0.d0 + szrhs = 0.d0 + taurhs = 0.d0 + + return + +end subroutine Whisky_TOV + + + +! ======================================================================= + + subroutine Spin1D(rad_func, r, nr, cart_func, x, y, z, & + nx, ny, nz) +! ................................................................ +! Takes a "radial grid function" and creates a 3D grid function +! BUG/FEATURE: This doesnt do any fancy interpolation at all! +! + implicit none + + integer nr, nx, ny, nz + CCTK_REAL rad_func(nr), r(nr) + CCTK_REAL cart_func(nx,ny,nz), x(nx,ny,nz), y(nx,ny,nz), z(nx,ny,nz) + + integer i,j,k, ri + real rloc + + do k = 1, nz + do j = 1, ny + do i = 1, nx + +! Given a point in 3D cartesian space, substitute next-further-out data +! value in 1D spherical space + rloc = sqrt(x(i,j,k)**2 + y(i,j,k)**2 + z(i,j,k)**2) + ri = 1 +100 continue + if (r(ri) .lt. rloc) then + ri = ri+1 + endif + + cart_func(i,j,k) = rad_func(ri) + end do + end do + end do + + end subroutine Spin1D + + + +! ======================================================================= + + + subroutine Spin1DMetric(grr,gtt,gxx,gxy, & + gxz,gyy,gyz,gzz,nx,ny,nz,x,y,z,r) + +! ................................................................ +! this subroutine converts spherical metric components to cartesian +! +! This code was totally ripped off of Francisco Guzmans spheretocart +! code. Used with his persmission. + + implicit NONE + integer nx,ny,nz + integer i,j,k + + CCTK_REAL x(nx,ny,nz),y(nx,ny,nz),z(nx,ny,nz),r(nx,ny,nz) + CCTK_REAL grr(nx,ny,nz),gtt(nx,ny,nz) + CCTK_REAL gxx(nx,ny,nz),gxy(nx,ny,nz),gxz(nx,ny,nz) + CCTK_REAL gyy(nx,ny,nz),gyz(nx,ny,nz),gzz(nx,ny,nz) + + CCTK_REAL xp,yp,zp,rp + +! define derivatives drx = (dr/dx) + CCTK_REAL drx,dry,drz + CCTK_REAL dtx,dty,dtz + CCTK_REAL dpx,dpy,dpz + +! this equation for gxx gyy gzz assumes that the metric is diagonal +! and that grr = 1 at the origin, but should compute for all points +! including the origin and the line x=y=0. + + + do k = 1, nz + do j = 1, ny + do i = 1, nx + + xp = x(i,j,k) + yp = y(i,j,k) + zp = z(i,j,k) + rp = r(i,j,k) + + if (rp.eq.0) then + + gxx(i,j,k) = 1. + gyy(i,j,k) = 1. + gzz(i,j,k) = 1. + gxy(i,j,k) = 0. + gxz(i,j,k) = 0. + gyz(i,j,k) = 0. + + else + + gxx(i,j,k) = (xp/rp)**2 * grr(i,j,k) + & + (zp**2 + yp**2)/(rp**2) + gyy(i,j,k) = (yp/rp)**2 * grr(i,j,k) + & + (xp**2 + zp**2)/(rp**2) + gzz(i,j,k) = (zp/rp)**2 * grr(i,j,k) + & + (xp**2 + yp**2)/(rp**2) + + if ((xp.eq.0) .and. (yp.eq.0)) then + + gxy(i,j,k) = 0. + gyz(i,j,k) = 0. + gxz(i,j,k) = 0. + + else + + gxy(i,j,k) = (xp * yp/(rp**2) * grr(i,j,k) & + + (xp*yp*zp**2)/(rp**2*(xp**2+yp**2)) & + - (yp*xp)/(xp**2+yp**2)) + + gxz(i,j,k) = (xp * zp/(rp**2) & + * grr(i,j,k) - (xp * zp/(rp**2))) + + gyz(i,j,k) = (yp * zp/(rp**2) & + * grr(i,j,k) - (yp * zp/(rp**2))) + + endif + + endif + + enddo + enddo + enddo +#else + write(6,*) 'Whisky_TOV: Does nothing. Recompile with odepack' +#endif + + return + end + + + diff --git a/src/make.code.defn b/src/make.code.defn new file mode 100644 index 0000000..2aefd19 --- /dev/null +++ b/src/make.code.defn @@ -0,0 +1,16 @@ +# Main make.code.defn file for thorn Hydra_Init_Data +# $Header: /Whisky/Whisky_Dev/Whisky_Init_Data/src/make.code.defn,v 1.6 2004/04/28 17:30:18 hawke Exp $ + +# Source files in this directory +SRCS = Whisky_C2P2C.F90 \ + Whisky_Con2Prim.F90 \ + Whisky_ReconstructTest.F90 \ + Whisky_ShockTube.F90 \ + Whisky_Only_Atmo.F90 \ + Whisky_ReadConformalData.F90 \ + Whisky_SimpleWave.F90 \ + CheckParam.c + +# Subdirectories containing source files +SUBDIRS = + diff --git a/src/make.configuration.defn b/src/make.configuration.defn new file mode 100644 index 0000000..87e4b0e --- /dev/null +++ b/src/make.configuration.defn @@ -0,0 +1,4 @@ +# Main make.code.defn file for thorn Hydra_Init_Data +# $Header$ + + diff --git a/src/make.configuration.deps b/src/make.configuration.deps new file mode 100644 index 0000000..19d41f1 --- /dev/null +++ b/src/make.configuration.deps @@ -0,0 +1,8 @@ + +# the TOV solver uses lsoda, which comes with the odepack library +# Currently the Cactus configure script doesn't check for this +# Thus the user needs to set an environment variable, ODEPACK_DIR +# -Scott Hawley + +# Remove by IH. + diff --git a/src/tov_fcn.F b/src/tov_fcn.F new file mode 100644 index 0000000..0e87845 --- /dev/null +++ b/src/tov_fcn.F @@ -0,0 +1,101 @@ +#ifdef HAVE_ODEPACK +c------------------------------------------------------------------------ +c equation of state, converts p to rho +c currently using a polytrope +c------------------------------------------------------------------------ + real*8 function p2rho(p) + implicit none + real*8 p + + p2rho = p ** (3/4.0d0) + + return + end + +c------------------------------------------------------------------------ +c equation of state, converts rho to p +c currently using a polytrope +c------------------------------------------------------------------------ + real*8 function rho2p(rho) + implicit none + real*8 rho + + rho2p = rho ** (4/3.0d0) + + return + end + + +c----------------------------------------------------------- +c Integrates static equations of motion for spherically +c symmetric, general-relativistic TOV star: +c +c dp/dr = - (p + rho)*(m+4*pi*r**3*p)/ r / (r-2*m) +c dm/dr = 4*pi*r**2 * rho +c d alpha/dr = alpha*(m+4*pi*r**3*rho )/ r / (r-2*m) +c +c (see, e.g. Walds textbook pp126-127.) +c where g_tt = - alpha**2, and g_rr = 1/(1-2m/r) +c +c Note that (my rho) = (the-rest-of-Whiskys rho)*(1+epsilon) +c +c----------------------------------------------------------- +c neq = 3 +c----------------------------------------------------------- + + subroutine tov_fcn(neq,r,y,yprime) + implicit none + + real*8 rho0 + common / com_tov / rho0 + + real*8 p2rho + + integer neq + real*8 r, y(neq), yprime(neq) + + real*8 pi + parameter ( pi = 3.1415926 5358 9793 d0 ) + + real*8 p, m, alf, rho, fourpr2 + + p = y(1) + m = y(2) + alf = y(3) + + if( r .eq. 0.0d0 ) then + yprime(1) = 0.0d0 + yprime(2) = 0.0d0 + yprime(3) = 0.0d0 + else + fourpr2 = 4*pi*r*r +c BUG: should be using call to CactusEOS + rho = p2rho(p) + + yprime(1) = - (p + rho)*(m+fourpr2*r*p)/ r / (r-2*m) + + yprime(2) = fourpr2 * rho + + yprime(3) = alf*(m+fourpr2*r*rho)/ r / (r-2*m) + endif + + return + end + +c----------------------------------------------------------- +c Implements Jacobian (optional) ... +c Does absolutely nothing +c----------------------------------------------------------- + subroutine tov_jac + implicit none + + real*8 rho0 + common / com_tov / rho0 + + return + end +#else + subroutine tov_fcn() + return + end +#endif diff --git a/src/tov_lsoda.F b/src/tov_lsoda.F new file mode 100644 index 0000000..b7a1e9f --- /dev/null +++ b/src/tov_lsoda.F @@ -0,0 +1,308 @@ +#ifdef HAVE_ODEPACK +c=========================================================== +c History: sode.f +c +c Driver routine which integrates ODEs defining +c TOV star +c +c Based on bstar.f by Matthew Choptuik +c=========================================================== + subroutine tov_lsoda ( rho, p, m, alf, + & r, nr, tol , n_succ, ssrc) + + implicit none + + + integer vsrc, vsxynt + real*8 rho2p, p2rho + + integer nr + real*8 rho(nr),p(nr), m(nr), alf(nr) + real*8 r(nr) + real*8 tol + integer n_succ, ssrc + + character*12 codenm + parameter ( codenm = 'tov_lsoda' ) + + integer iargc, indlnb, i4arg + real*8 r8arg + + real*8 r8_never + parameter ( r8_never = -1.0d-60 ) + +c----------------------------------------------------------- +c Order of system. +c----------------------------------------------------------- + integer neq + parameter ( neq = 3 ) + +c----------------------------------------------------------- +c Storage for solution at requested output radii. +c----------------------------------------------------------- + integer maxout + parameter ( maxout = 10 000 ) + + real*8 y0(neq) + real*8 vxout(maxout), vyout(maxout,neq), + & work(maxout) + integer nxout, ixout, nxout_succ + + integer ieq + + logical ltrace + parameter ( ltrace = .false. ) + + integer maxdump + parameter ( maxdump = 50 ) + + logical lsodatrace + parameter ( lsodatrace = .false. ) + + +c----------------------------------------------------------- +c LSODA Variables. +c----------------------------------------------------------- + external tov_fcn, tov_jac + real*8 y(neq) + real*8 tbgn, tend + integer itol + real*8 rtol, atol + integer itask, istate, iopt + integer lrw + +c parameter ( lrw = 22 + neq * max(16 ,neq + 9) ) + parameter ( lrw = 22 + neq * (16) ) + real*8 rwork(lrw) + + integer liw + parameter ( liw = 20 + neq ) + integer iwork(liw) + integer jt + + integer i + +c----------------------------------------------------------- +c Common communication with routine 'tov_fcn' in 'tov_fcn.f' ... +c----------------------------------------------------------- +c rho0: Value of rho at r=0 +c----------------------------------------------------------- + real*8 rho0 + common / com_tov / rho0 + +c----------------------------------------------------------- +c Parse arguments. +c----------------------------------------------------------- + rho0 = rho(1) + + nxout = nr + do ixout = 1, nxout + vxout(ixout) = r(ixout) + enddo + +c----------------------------------------------------------- +c Initialize those "boundary" conditions which are +c fixed by regularity, or which are arbitrary. +c----------------------------------------------------------- + y0(1) = rho2p(rho0) + y0(2) = 0.0d0 + y0(3) = 1.0 + + +c----------------------------------------------------------- +c Echo command line arguments if "local tracing" +c enabled ... +c----------------------------------------------------------- + if( ltrace ) then + write(0,*) 'rho0: ', rho0 + write(0,*) 'p(1): ', y0(1) + write(0,*) 'm(1): ', y0(2) + write(0,*) 'alf(1): ', y0(3) + end if + +c----------------------------------------------------------- +c Get output radii +c----------------------------------------------------------- + if( ltrace ) then + if( nxout .le. maxdump ) then + call dvdump(vxout,nxout,codenm//': output radii',0) + else + call dvdmp1(vxout,work, + & int(1.0d0 * nxout / maxdump + 0.5d0), + & nxout,codenm//': selected output radii',0) + end if + write(0,*) + write(0,*) codenm//': Initial time: ', vxout(1) + write(0,*) + end if + +c----------------------------------------------------------- +c Set LSODA parameters ... +c----------------------------------------------------------- + itol = 1 + rtol = tol + atol = tol + itask = 1 + iopt = 0 + jt = 2 + +c----------------------------------------------------------- +c Initialize the solution ... +c----------------------------------------------------------- + do ieq = 1 , neq + y(ieq) = y0(ieq) + vyout(1,ieq) = y0(ieq) + end do + +c----------------------------------------------------------- +c Do the integration ... +c----------------------------------------------------------- + ssrc = 0 + do ixout = 2 , nxout + istate = 1 +c----------------------------------------------------------- +c Need these temporaries since lsoda overwrites +c tend ... +c----------------------------------------------------------- +100 continue + tbgn = vxout(ixout-1) + tend = vxout(ixout) + call lsoda(tov_fcn,neq,y,tbgn,tend, + & itol,rtol,atol,itask, + & istate,iopt,rwork,lrw,iwork,liw,tov_jac,jt) + if( lsodatrace ) then + write(0,1000) codenm, ixout, vxout(ixout), + & vxout(ixout+1) +1000 format(' ',a,': Step ',i4,' t = ',1pe10.3, + & ' .. ',1pe10.3) + write(0,*) codenm//': lsoda reurns ', istate + end if + +c if (istate .eq. -2) then +c write(0,*) 'istate = -2, trying again...' +c itol = 1 +c rtol = tol +c atol = tol +c itask = 1 +c iopt = 0 +c jt = 2 +c istate = 3 +c goto 100 +c endif + + if( istate .lt. -1 ) then + write(0,1500) codenm, istate, ixout, nxout, + & vxout(ixout-1), vxout(ixout) +1500 format(/' ',a,': Error return ',i2, + & ' from integrator LSODA.'/ + & ' At output time ',i5,' of ',i5/ + & ' Interval ',1pe11.3,' .. ', + & 1pe11.3/) + nxout_succ = ixout - 1 + ssrc = 1 + go to 500 + + end if + + do ieq = 1 , neq + vyout(ixout,ieq) = y(ieq) + end do + + end do + nxout_succ = nxout + +500 continue + +c----------------------------------------------------------- +c OUTPUT +c----------------------------------------------------------- + do ixout = 1 , nxout_succ + p(ixout) = vyout(ixout, 1) + m(ixout) = vyout(ixout, 2) + alf(ixout) = vyout(ixout, 3) + rho(ixout) = p2rho(p(ixout)) + end do + + n_succ = nxout_succ + + return + + end + + + + + +c=========================================================== +c Some double precision vector utility routines. +c +c Originally from ~matt/vutil/dveclib.f +c=========================================================== + +c----------------------------------------------------------- +c Dumps vector labelled with LABEL on UNIT. +c----------------------------------------------------------- + subroutine dvdump(v,n,label,unit) + implicit none + + real*8 v(*) + character*(*) label + integer i, n, st, unit + + if( n .lt. 1 ) go to 130 + write(unit,100) label + 100 format(/' <<< ',A,' >>>'/) + st = 1 + 110 continue + write(unit,120) ( v(i) , i = st , min(st+3,n)) + 120 FORMAT(' ',4(1PE19.10)) + st = st + 4 + if( st .le. n ) go to 110 + + 130 continue + + return + end + +c----------------------------------------------------------- +c Extension of dvdump which dumps every 'inc'th element. +c----------------------------------------------------------- + subroutine dvdmp1(v,w,inc,n,label,unit) + implicit none + + real*8 v(*), w(*) + character*(*) label + integer inc, n, unit + + call dvinj(v,w,inc,n) + call dvdump(w,1+(n-1)/inc,label,unit) + + return + + end + +c--------------------------------------------------------- +c Injects every incth element of v1 into v2 +c--------------------------------------------------------- + subroutine dvinj(v1,v2,inc,n) + implicit none + + real*8 v1(*), v2(*) + integer i, inc, j, n + + j = 1 + do i = 1 , n , inc + v2(j) = v1(i) + j = j + 1 + end do + + + return + end + +#else + subroutine tov_lsoda() + return + end + +#endif |