Commits January 2011

commits@lists.einsteintoolkit.org

7 participants
57 discussions

[svn:einsteintoolkit] Paper_EinsteinToolkit_2010/ (Rev. 6)
by schnetter＠cct.lsu.edu 17 Jan '11

17 Jan '11

User: eschnett Date: 2011/01/17 01:35 PM Modified: / ET.tex Log: Correct latex error. Put author list into comment. Remove all other bibliography files File Changes: Directory: / ============ File [modified]: ET.tex Delta lines: +15 -4 =================================================================== --- ET.tex 2011-01-17 16:49:52 UTC (rev 5) +++ ET.tex 2011-01-17 19:35:08 UTC (rev 6) @@ -54,9 +54,20 @@ \title[]{The Einstein Toolkit: A Community Computational Infrastructure for Relativistic Astrophysics} -\author{\textbf{check, sort}} +\author{\textbf{TBD}} +% POTENTIAL AUTHORS: M. Alcubierre, G. Allen, M. Ansorg, L. Baiotti, +% E. Bentivegna, S. Bernuzzi, T. Bode, S. Brandt, G. Corvino, R. De +% Pietri, P. Diener, H. Dimmelmeier, Y. El Khamra, T. Font, J. +% Frieben, B. Giacomazzo, T. Goodale, C. Gundlach, R. Haas, I. Hawke, +% S. Hawley, I. Hinder, T. Kellermann, G. Lanfermann, F. L\"{o}ffler, +% J. Masso, M. Miller, P. Montero, C. Ott, R. Paruchuri, D. Pollney, +% T. Radke, L. Rezzolla, D. Rideout, M. Ripeanu, E. Schnetter, E. +% Seidel, E. L. Seidel, J. Shalf, N. Stergioulas, R. Takahashi, J. +% Thornburg, A. Tonita, P. Walker + + \thanks{{}$^{(1)}$ Department of Computer Science, Louisiana State University, Baton Rouge, LA 70803, USA} @@ -361,7 +372,7 @@ \section{Examples} \todo{Frank} -\pararaph{Kerr-Schild} +\paragraph{Kerr-Schild} \paragraph{BH-binary} @@ -377,9 +388,9 @@ \section*{Acknowledgments} \todo{Frank, ADDME} + \bibliographystyle{amsplain-url} +\bibliography{manifest/einsteintoolkit} -%\bibliography{manifest/einsteintoolkit.bib,local_bibtex/CCT_CS,local_bibtex/references,local_bibtex/gw_references} -\bibliography{manifest/einsteintoolkit.bib,local_bibtex/CCT_CS,local_bibtex/references} \end{document}

1 0

[svn:einsteintoolkit] Paper_EinsteinToolkit_2010/ (Rev. 5)
by knarf＠cct.lsu.edu 17 Jan '11

17 Jan '11

User: knarf Date: 2011/01/17 10:49 AM Modified: / ET.tex Log: assignements, restructuring File Changes: Directory: / ============ File [modified]: ET.tex Delta lines: +24 -16 =================================================================== --- ET.tex 2011-01-17 15:55:53 UTC (rev 4) +++ ET.tex 2011-01-17 16:49:52 UTC (rev 5) @@ -54,12 +54,7 @@ \title[]{The Einstein Toolkit: A Community Computational Infrastructure for Relativistic Astrophysics} -\author{M. Alcubierre, G. Allen$^{(1,2,3)}$, M. Ansorg, L. Baiotti, E. Bentivegna$^{(1,2)}$, S. Bernuzzi, T. Bode, S. Brandt, -G. Corvino, R. De Pietri, P. Diener, H. Dimmelmeier, Y. El Khamra, T. Font, J. Frieben, -B. Giacomazzo, T. Goodale, C. Gundlach, R. Haas, I. Hawke, S. Hawley, I. Hinder, T. Kellermann, -G. Lanfermann, F. L\"{o}ffler$^{(3)}$, J. Masso, M. Miller, P. Montero, C. Ott, - R. Paruchuri$^{(3)}$, D. Pollney, T. Radke, L. Rezzolla, D. Rideout, M. Ripeanu, E. Schnetter$^{(2,3)}$, -E. Seidel, E. L. Seidel, J. Shalf, N. Stergioulas, R. Takahashi, J. Thornburg, A. Tonita, P. Walker\textbf{check, sort}} +\author{\textbf{check, sort}} \thanks{{}$^{(1)}$ Department of Computer Science, @@ -241,6 +236,7 @@ \section{Core Technologies} \subsection{Cactus Framework} +\todo{Gab} The \textbf{Cactus Framework}~\cite{CS_cactus_web,CS_Goodale02a,CS_cactususersguide} is @@ -260,6 +256,7 @@ \subsection{Carpet Mesh Refinement} +\todo{Erik} While Cactus is distributed with a structured-mesh unigrid MPI parallel driver (\codename{PUGH}\footnote{Recent results with a @@ -281,7 +278,8 @@ source and is openly developed, with the main development located at LSU (co-PI Schnetter) and contributions from AEI and others. -\subsection{CactusEinstein} The Cactus Framework was developed by the +\subsection{CactusEinstein} \todo{remove extra subsection for this, +move further up} The Cactus Framework was developed by the numerical relativity community, and although it is a general component framework that supports different application domains its core user group has remained from numerical relativity. The Cactus Team have @@ -295,17 +293,17 @@ The Einstein Toolkit collects the widely used parts of CactusEinstein, combined with contributions from the community. -\subsection{Whisky Code} - \subsection{Simulation Factory} +\todo{Erik} \subsection{Kranc} +\todo{Ian} - \section{Components} \subsection{Base Modules} +\todo{Frank} \paragraph{Vacuum Spacetimes} @@ -313,6 +311,7 @@ \subsection{Initial Data} +\todo{Josh} \paragraph{Gravitational Waves} @@ -322,6 +321,7 @@ \subsection{Evolution Methods} +\todo{Christian} \paragraph{Spacetime} @@ -329,6 +329,7 @@ \subsection{Analysis} +\todo{Tanja} \paragraph{Masses} @@ -338,6 +339,7 @@ \subsection{Relativity Tools} +\todo{Peter} \paragraph{Black Hole Excision} @@ -345,8 +347,10 @@ \paragraph{Object tracking} -\subsection{Numerical Methods} +\subsection{Infrastructure and Numerical Methods} +\todo{Erik} + \paragraph{Domains and Coordinates} \paragraph{Adaptive Mesh Refinement} @@ -355,19 +359,23 @@ \section{Examples} +\todo{Frank} -kerr-schild, BH-binary, TOV +\pararaph{Kerr-Schild} -\section{Results} +\paragraph{BH-binary} -Some success stories ... state of art of what Toolkit can do now, scaling results. +\paragraph{TOV+EOS} +\paragraph{Collapse} + + \section{Future Work} +\todo{Frank} - \section*{Acknowledgments} -\todo{ADDME} +\todo{Frank, ADDME} \bibliographystyle{amsplain-url}

1 0

[svn:einsteintoolkit] Paper_EinsteinToolkit_2010/ (Rev. 4)
by knarf＠cct.lsu.edu 17 Jan '11

17 Jan '11

User: knarf Date: 2011/01/17 09:55 AM Modified: / ET.tex Log: change structure a bit File Changes: Directory: / ============ File [modified]: ET.tex Delta lines: +41 -24 =================================================================== --- ET.tex 2010-12-14 21:49:21 UTC (rev 3) +++ ET.tex 2011-01-17 15:55:53 UTC (rev 4) @@ -59,7 +59,7 @@ B. Giacomazzo, T. Goodale, C. Gundlach, R. Haas, I. Hawke, S. Hawley, I. Hinder, T. Kellermann, G. Lanfermann, F. L\"{o}ffler$^{(3)}$, J. Masso, M. Miller, P. Montero, C. Ott, R. Paruchuri$^{(3)}$, D. Pollney, T. Radke, L. Rezzolla, D. Rideout, M. Ripeanu, E. Schnetter$^{(2,3)}$, -E. Seidel, E. L. Seidel, J. Shalf, N. Stergioulas, R. Takahashi, J. Thornburg, A. Tonita, P. Walker} +E. Seidel, E. L. Seidel, J. Shalf, N. Stergioulas, R. Takahashi, J. Thornburg, A. Tonita, P. Walker\textbf{check, sort}} \thanks{{}$^{(1)}$ Department of Computer Science, @@ -84,9 +84,13 @@ \section{Introduction} -Scientific progress in the field of numerical relativity has always been closely tied with the availability and ease-of-use of enabling software and computational infrastructure. This document describes the Einstein Toolkit. +Scientific progress in the field of numerical relativity has always +been closely tied with the availability and ease-of-use of enabling +software and computational infrastructure. This document describes +the Einstein Toolkit, which aims at providing such an infrastructure, +developed openly and available freely. -This is a particularly exciting time for numerical relativity and +Now is a particularly exciting time for numerical relativity and relativistic astrophysics. Recent computational breakthroughs~\cite{Pretorius:2005gq, Campanelli:2005dd, Baker:2005vv} have transformed the field with our ability to @@ -137,8 +141,8 @@ collaborators, including co-PI Faber, e.g., \cite{Etienne:2007jg,Liu:2008xy}), LSU (Lehner et al., \cite{Anderson:2007kz,Anderson:2008zp}), Caltech-Cornell (Duez et -al.~\cite{Duez:2008rb}), and at the AEI~(Rezzolla and collaborators, -\cite{Baiotti:2008ra,loeffler:06}). The first simulations in 3D GRHD +al.~\cite{Duez:2008rb}), and at the AEI and Sissa~(Rezzolla and collaborators, +\cite{Baiotti:2008ra,Loeffler06a}). The first simulations in 3D GRHD of massive star collapse to protoneutron stars have been carried out by co-PI Ott and collaborators \cite{ott:07prl,ott:07cqg} and by Shibata~et~al.~\cite{shibata:053D}. The collapse of rotating, @@ -223,11 +227,9 @@ \subsection{Academic and Social} \begin{itemize} - -\item A primary concern for research groups is securing reliable funding to maintain graduate students and postdoctoral researchers. - +\item A primary concern for research groups is securing reliable funding + to maintain graduate students and postdoctoral researchers. \item - \end{itemize} \section{Design and Strategy for the Einstein Toolkit} @@ -257,7 +259,7 @@ application domains. -\subsection{Carpet} +\subsection{Carpet Mesh Refinement} While Cactus is distributed with a structured-mesh unigrid MPI parallel driver (\codename{PUGH}\footnote{Recent results with a @@ -279,10 +281,20 @@ source and is openly developed, with the main development located at LSU (co-PI Schnetter) and contributions from AEI and others. +\subsection{CactusEinstein} The Cactus Framework was developed by the +numerical relativity community, and although it is a general component +framework that supports different application domains its core user +group has remained from numerical relativity. The Cactus Team have +traditionally developed and supported a set of core modules for numerical +relativity, as part of the {\tt CactusEinstein} arrangement. Over the +last few years however, the relevance of many of the modules has declined, +and more and more of the basic infrastructure for numerical relativity +has been provided by open modules provided and distributed by research +groups in the community. +The Einstein Toolkit collects the widely used parts of CactusEinstein, +combined with contributions from the community. -\subsection{CactusEinstein} The Cactus Framework was developed by the numerical relativity community, and although it is a general component framework that supports different application domains its core user group has remained from numerical relativity. The Cactus Team have traditionally developed and supported a set of core modules for numerical relativity, as part of the {\tt CactusEinstein} arrangement. Over the last few years however, the relevance of many of the modules has declined, and more and more of the basic infrastructure for numerical relativity has been provided by open modules provided and distributed by research groups in the community. - \subsection{Whisky Code} \subsection{Simulation Factory} @@ -291,27 +303,31 @@ +\section{Components} +\subsection{Base Modules} +\paragraph{Vacuum Spacetimes} -\section{Components} +\paragraph{Relativistic Matter Spacetimes} -\subsection{Vacuum Spacetimes} -\subsection{Relativistic Matter Spacetimes} - -\subsection{Gauge Conditions} - \subsection{Initial Data} \paragraph{Gravitational Waves} -\paragraph{Black Hole} +\paragraph{Black Holes} -\paragraph{Perturbed Black Hole} - \paragraph{Neutron Stars} + +\subsection{Evolution Methods} + +\paragraph{Spacetime} + +\paragraph{Relativistic Matter} + + \subsection{Analysis} \paragraph{Masses} @@ -320,13 +336,15 @@ \paragraph{Gravitational Waves} + \subsection{Relativity Tools} \paragraph{Black Hole Excision} - \paragraph{Vacuum-Matter Coupling} +\paragraph{Object tracking} + \subsection{Numerical Methods} \paragraph{Domains and Coordinates} @@ -336,10 +354,9 @@ \paragraph{Boundary Conditions} - \section{Examples} -kerr-schild, binary, TOV +kerr-schild, BH-binary, TOV \section{Results}

1 0

[svn:einsteintoolkit] manifest/trunk/ (Rev. 61)
by knarf＠cct.lsu.edu 17 Jan '11

17 Jan '11

User: knarf Date: 2011/01/17 09:33 AM Modified: /trunk/ einsteintoolkit.th Log: add the PunctureTracker to the ET File Changes: Directory: /trunk/ ================== File [modified]: einsteintoolkit.th Delta lines: +1 -0 =================================================================== --- trunk/einsteintoolkit.th 2010-12-28 17:13:56 UTC (rev 60) +++ trunk/einsteintoolkit.th 2011-01-17 15:33:57 UTC (rev 61) @@ -232,6 +232,7 @@ !URL = http://svn.aei.mpg.de/numrel/$1/$2/trunk !CHECKOUT = AEIThorns/AEILocalInterp +AEIThorns/PunctureTracker # Various thorns from LSU

1 0

[svn:einsteintoolkit] TwoPunctures/trunk/src/ (Rev. 119)
by knarf＠cct.lsu.edu 10 Jan '11

10 Jan '11

User: knarf Date: 2011/01/10 01:25 PM Modified: /trunk/src/ TwoPunctures.c Log: remove compiler warnings (declare functions, rename x1,y1,z1 (because of bessel functions in math libraries), initialize variables) File Changes: Directory: /trunk/src/ ====================== File [modified]: TwoPunctures.c Delta lines: +18 -17 =================================================================== --- trunk/src/TwoPunctures.c 2011-01-10 19:20:14 UTC (rev 118) +++ trunk/src/TwoPunctures.c 2011-01-10 19:25:56 UTC (rev 119) @@ -193,6 +193,7 @@ } /* -------------------------------------------------------------------*/ +void TwoPunctures(CCTK_ARGUMENTS); void TwoPunctures (CCTK_ARGUMENTS) { @@ -402,10 +403,10 @@ const int ind = CCTK_GFINDEX3D (cctkGH, i, j, k); - CCTK_REAL x1, y1, z1; - x1 = x[ind] - center_offset[0]; - y1 = y[ind] - center_offset[1]; - z1 = z[ind] - center_offset[2]; + CCTK_REAL xx, yy, zz; + xx = x[ind] - center_offset[0]; + yy = y[ind] - center_offset[1]; + zz = z[ind] - center_offset[2]; /* We implement swapping the x and z coordinates as follows. The bulk of the code that performs the actual calculations @@ -415,24 +416,24 @@ main loop-- we swap everything back. */ if (swap_xz) { /* Swap the x and z coordinates */ - SWAP (x1, z1); + SWAP (xx, zz); } CCTK_REAL r_plus - = sqrt(pow(x1 - par_b, 2) + pow(y1, 2) + pow(z1, 2)); + = sqrt(pow(xx - par_b, 2) + pow(yy, 2) + pow(zz, 2)); CCTK_REAL r_minus - = sqrt(pow(x1 + par_b, 2) + pow(y1, 2) + pow(z1, 2)); + = sqrt(pow(xx + par_b, 2) + pow(yy, 2) + pow(zz, 2)); CCTK_REAL U; switch (gsm) { case GSM_Taylor_expansion: U = PunctTaylorExpandAtArbitPosition - (0, nvar, n1, n2, n3, v, x1, y1, z1); + (0, nvar, n1, n2, n3, v, xx, yy, zz); break; case GSM_evaluation: U = PunctIntPolAtArbitPosition - (0, nvar, n1, n2, n3, v, x1, y1, z1); + (0, nvar, n1, n2, n3, v, xx, yy, zz); break; default: assert (0); @@ -463,9 +464,9 @@ CCTK_REAL static_psi = 1; CCTK_REAL Aij[3][3]; - BY_Aijofxyz (x1, y1, z1, Aij); + BY_Aijofxyz (xx, yy, zz, Aij); - CCTK_REAL old_alp; + CCTK_REAL old_alp=1.0; if (multiply_old_lapse) old_alp = alp[ind]; @@ -480,9 +481,9 @@ pyy = pyz = pzz = 0.0; /* first puncture */ - xp = x1 - par_b; - yp = y1; - zp = z1; + xp = xx - par_b; + yp = yy; + zp = zz; rp = sqrt (xp*xp + yp*yp + zp*zp); rp = pow (pow (rp, 4) + pow (TP_epsilon, 4), 0.25); if (rp < TP_Tiny) @@ -511,9 +512,9 @@ pzz += zp*zp*s5 + s3; /* second puncture */ - xp = x1 + par_b; - yp = y1; - zp = z1; + xp = xx + par_b; + yp = yy; + zp = zz; rp = sqrt (xp*xp + yp*yp + zp*zp); rp = pow (pow (rp, 4) + pow (TP_epsilon, 4), 0.25); if (rp < TP_Tiny)

1 0

[svn:einsteintoolkit] TwoPunctures/trunk/src/ (Rev. 118)
by knarf＠cct.lsu.edu 10 Jan '11

10 Jan '11

User: knarf Date: 2011/01/10 01:20 PM Modified: /trunk/src/ TwoPunctures.c Log: patch from Bernard Kelly: Output the individual puncture ADM mass post-solve in TwoPunctures File Changes: Directory: /trunk/src/ ====================== File [modified]: TwoPunctures.c Delta lines: +12 -1 =================================================================== --- trunk/src/TwoPunctures.c 2010-11-21 03:29:55 UTC (rev 117) +++ trunk/src/TwoPunctures.c 2011-01-10 19:20:14 UTC (rev 118) @@ -219,6 +219,7 @@ CCTK_REAL admMass; if (! F) { + CCTK_REAL Mp_adm, Mm_adm, up, um; /* Solve only when called for the first time */ F = dvector (0, ntotal - 1); allocate_derivs (&u, ntotal); @@ -255,7 +256,7 @@ target ADM masses target_M_plus and target_M_minus and with initial guesses given by par_m_plus and par_m_minus. */ if(!(give_bare_mass)) { - CCTK_REAL tmp, Mp_adm, Mm_adm, Mp_adm_err, Mm_adm_err, up, um; + CCTK_REAL tmp, Mp_adm_err, Mm_adm_err; char valbuf[100]; CCTK_REAL M_p = target_M_plus; @@ -316,6 +317,16 @@ "The two puncture masses are mp=%.17g and mm=%.17g", (double) *mp, (double) *mm); + up = PunctIntPolAtArbitPosition(0, nvar, n1, n2, n3, v, par_b, 0., 0.); + um = PunctIntPolAtArbitPosition(0, nvar, n1, n2, n3, v,-par_b, 0., 0.); + + /* Calculate the ADM masses from the current bare mass guess */ + Mp_adm = (1 + up) * *mp + *mp * *mm / (4. * par_b); + Mm_adm = (1 + um) * *mm + *mp * *mm / (4. * par_b); + + CCTK_VInfo (CCTK_THORNSTRING, "Puncture 1 ADM mass is %g", (double)Mp_adm); + CCTK_VInfo (CCTK_THORNSTRING, "Puncture 2 ADM mass is %g", (double)Mm_adm); + /* print out ADM mass, eq.: \Delta M_ADM=2*r*u=4*b*V for A=1,B=0,phi=0 */ admMass = (*mp + *mm - 4*par_b*PunctEvalAtArbitPosition(v.d0, 0, 1, 0, 0, nvar, n1, n2, n3));

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 210)
by knarf＠cct.lsu.edu 07 Jan '11

07 Jan '11

User: knarf Date: 2011/01/07 09:59 PM Modified: /trunk/src/ GRHydro_Prim2Con.F90 Log: correct variable declaration File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Prim2Con.F90 Delta lines: +1 -1 =================================================================== --- trunk/src/GRHydro_Prim2Con.F90 2011-01-08 01:37:53 UTC (rev 209) +++ trunk/src/GRHydro_Prim2Con.F90 2011-01-08 03:59:58 UTC (rev 210) @@ -46,7 +46,7 @@ integer :: i, j, k CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr - CCTK_REAL :: xtemp + CCTK_REAL :: xtemp(1) if(evolve_temper.ne.1) then !$OMP PARALLEL DO PRIVATE(i, j, avg_detl, avg_detr,&

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 209)
by cott＠tapir.caltech.edu 07 Jan '11

07 Jan '11

User: cott Date: 2011/01/07 07:37 PM Modified: /trunk/src/ GRHydro_Reconstruct.F90 Log: * fix a very nasty bug in PPM reconstruction that sneeked in with the very recent PPM changes. Old version of the code were not affected. File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Reconstruct.F90 Delta lines: +3 -3 =================================================================== --- trunk/src/GRHydro_Reconstruct.F90 2011-01-06 01:13:02 UTC (rev 208) +++ trunk/src/GRHydro_Reconstruct.F90 2011-01-08 01:37:53 UTC (rev 209) @@ -361,7 +361,7 @@ if(clean_divergence.ne.0) then !$OMP PARALLEL DO PRIVATE(i, j) do k = GRHydro_stencil, nz - GRHydro_stencil + 1 - do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 call SimplePPM_1dM(GRHydro_eos_handle,0,ny,CCTK_DELTA_SPACE(2),& rho(j,:,k),vely(j,:,k),velz(j,:,k),velx(j,:,k),& Bvecy(j,:,k),Bvecz(j,:,k),Bvecx(j,:,k),psidc(j,:,k),eps(j,:,k),press(j,:,k),& @@ -389,7 +389,7 @@ else !clean_divergence !$OMP PARALLEL DO PRIVATE(i, j) do k = GRHydro_stencil, nz - GRHydro_stencil + 1 - do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 call SimplePPM_1dM(GRHydro_eos_handle,0,ny,CCTK_DELTA_SPACE(2),& rho(j,:,k),vely(j,:,k),velz(j,:,k),velx(j,:,k),& Bvecy(j,:,k),Bvecz(j,:,k),Bvecx(j,:,k),dum(j,:,k),eps(j,:,k),press(j,:,k),& @@ -418,7 +418,7 @@ else !evolve_mhd !$OMP PARALLEL DO PRIVATE(i, j) do k = GRHydro_stencil, nz - GRHydro_stencil + 1 - do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + do j = GRHydro_stencil, nx - GRHydro_stencil + 1 call SimplePPM_1d(GRHydro_eos_handle,0,ny,CCTK_DELTA_SPACE(2),& rho(j,:,k),vely(j,:,k),velz(j,:,k),velx(j,:,k),eps(j,:,k),& press(j,:,k),rhominus(j,:,k),velyminus(j,:,k),velzminus(j,:,k),&

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 208)
by bcmsma＠astro.rit.edu 05 Jan '11

05 Jan '11

User: bmundim Date: 2011/01/05 07:13 PM Modified: /trunk/src/ GRHydro_Con2PrimM.F90, GRHydro_Con2PrimM_polytype_pt.c, GRHydro_Con2PrimM_pt.c Log: RIT MHD development: Fix race conditions arising in the Con2Prim routines. Now both GNU and Intel compilers agree on the results for the magnetic tests. File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Con2PrimM.F90 Delta lines: +7 -5 =================================================================== --- trunk/src/GRHydro_Con2PrimM.F90 2011-01-04 14:40:06 UTC (rev 207) +++ trunk/src/GRHydro_Con2PrimM.F90 2011-01-06 01:13:02 UTC (rev 208) @@ -94,7 +94,7 @@ !$OMP PARALLEL DO PRIVATE(i,j,itracer,& !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, epsnegative, & - !$OMP b2,local_gam,xrho,xpress,local_K,local_pgam,sc) + !$OMP b2,xrho,xeps,xpress,xtemp,local_K,local_pgam,sc,keyerr,anyerr) do k = 1, nz do j = 1, ny do i = 1, nx @@ -142,6 +142,7 @@ scon(i,j,k,:) = 0.d0 vel(i,j,k,:) = 0.d0 w_lorentz(i,j,k) = 1.d0 + xtemp = 0.0d0 call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr) @@ -225,13 +226,14 @@ call CCTK_WARN(GRHydro_NaN_verbose+2, 'Specific internal energy just went below 0, trying polytype.') !$OMP END CRITICAL - xrho=10.0d0 + xrho=1.0d0; xtemp=0.0d0; xeps=1.0d0 call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& - 1.d0,1.0d0,xtemp,xye,xpress,keyerr,anyerr) - local_K = xpress + xrho,xeps,xtemp,xye,xpress,keyerr,anyerr) + local_K = xpress; + xrho=10.0d0; xeps=1.0d0 call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& - xrho,1.0d0,xtemp,xye,xpress,keyerr,anyerr) + xrho,xeps,xtemp,xye,xpress,keyerr,anyerr) local_pgam=log(xpress/local_K)/log(xrho) sc = local_K*dens(i,j,k) File [modified]: GRHydro_Con2PrimM_polytype_pt.c Delta lines: +83 -166 =================================================================== --- trunk/src/GRHydro_Con2PrimM_polytype_pt.c 2011-01-04 14:40:06 UTC (rev 207) +++ trunk/src/GRHydro_Con2PrimM_polytype_pt.c 2011-01-06 01:13:02 UTC (rev 208) @@ -71,43 +71,33 @@ ***************************************************/ /* Local Globals */ -static CCTK_REAL Bsq,QdotBsq,Qtsq,Qdotn,D,half_Bsq,rho_gm1; -static CCTK_REAL W_for_gnr2, rho_for_gnr2, W_for_gnr2_old, rho_for_gnr2_old, drho_dW; -static CCTK_REAL t2,t4,t7,t24,two_Bsq,t300,t400,s200,s100; +struct LocGlob { + CCTK_REAL Bsq,QdotBsq,Qtsq,D; + CCTK_REAL W_for_gnr2, rho_for_gnr2, W_for_gnr2_old, rho_for_gnr2_old; + CCTK_REAL t4,t7,t24,two_Bsq,t300,t400,s200,s100; +} ; // Declarations: -static CCTK_REAL vsq_calc(CCTK_REAL W); - -static CCTK_INT oned_newton_raphson( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos, +static CCTK_INT oned_newton_raphson( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos, struct LocGlob *lgp, void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][1], - CCTK_REAL *, CCTK_REAL *, CCTK_INT, CCTK_REAL) ); + CCTK_REAL *, CCTK_REAL *, CCTK_INT, CCTK_REAL, + struct LocGlob *lgp) ); static void func_W( CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][1], - CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos); -static CCTK_REAL x1_of_x0(CCTK_REAL x0 ) ; + CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos, + struct LocGlob * lgp); -static CCTK_INT gnr2( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos, +static CCTK_INT gnr2( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos, struct LocGlob *lgp, void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], - CCTK_REAL [][1],CCTK_REAL *,CCTK_REAL *,CCTK_INT, CCTK_REAL) ); + CCTK_REAL [][1],CCTK_REAL *,CCTK_REAL *,CCTK_INT, CCTK_REAL, + struct LocGlob *lgp) ); + static void func_rho(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], - CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos); + CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos, + struct LocGlob *lgp); -static CCTK_REAL eos_info(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL *dpdw, CCTK_REAL *dpdvsq, CCTK_REAL gammaeos); -/* pressure as a function of rho0 and u */ -static CCTK_REAL pressure_rho0_u(CCTK_REAL rho0, CCTK_REAL u, CCTK_REAL gammaeos) -{ - return((gammaeos - 1.)*u) ; -} - -/* Pressure as a function of rho0 and w = rho0 + u + p */ -static CCTK_REAL pressure_rho0_w(CCTK_REAL rho0, CCTK_REAL w, CCTK_REAL gammaeos) -{ - return((gammaeos-1.)*(w - rho0)/gammaeos) ; -} - - void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_Polytype_pt) ( CCTK_INT *handle, CCTK_REAL *gamma_eos, CCTK_REAL *dens_in, @@ -220,8 +210,12 @@ CCTK_REAL detg = (*det); CCTK_REAL sqrt_detg = sqrt(detg); CCTK_REAL inv_sqrt_detg = 1./sqrt_detg; + CCTK_REAL t2, rho_gm1; CCTK_INT i,j, i_increase,ntries ; + struct LocGlob lg; + + /* Assume ok initially: */ *retval = 0.; *epsnegative = 0; @@ -275,7 +269,7 @@ // Calculate various scalars (Q.B, Q^2, etc) from the conserved variables: - Bsq = + lg.Bsq = (*gxx) * (*Bx) * (*Bx) + (*gyy) * (*By) * (*By) + (*gzz) * (*Bz) * (*Bz) + @@ -285,25 +279,23 @@ (*gyz) * (*By) * (*Bz) ); QdotB = (sx * (*Bx) + sy * (*By) + sz * (*Bz)) ; - QdotBsq = QdotB*QdotB ; + lg.QdotBsq = QdotB*QdotB ; - Qtsq = (usx * sx + usy * sy + usz * sz) ; + lg.Qtsq = (usx * sx + usy * sy + usz * sz) ; - D = dens; + lg.D = dens; - half_Bsq = 0.5*Bsq; + t2 = lg.D*lg.D; + lg.t4 = lg.QdotBsq*t2; + lg.t7 = lg.Bsq*lg.Bsq; + lg.t24 = 1/t2; + lg.two_Bsq = lg.Bsq + lg.Bsq; + lg.t300 = lg.QdotBsq*lg.Bsq*t2; + lg.t400 = lg.Qtsq*t2; - t2 = D*D; - t4 = QdotBsq*t2; - t7 = Bsq*Bsq; - t24 = 1/t2; - two_Bsq = Bsq + Bsq; - t300 = QdotBsq*Bsq*t2; - t400 = Qtsq*t2; - - s200 = D*gammaeos*(*sc_in); + lg.s200 = lg.D*gammaeos*(*sc_in); CCTK_REAL g_o_gm1 = (gammaeos/(gammaeos-1.)); - s100 = g_o_gm1*(*sc_in); + lg.s100 = g_o_gm1*(*sc_in); /* calculate W from last timestep and use for guess */ vsq = @@ -328,7 +320,7 @@ // Always calculate rho from D and gamma so that using D in EOS remains consistent // i.e. you don't get positive values for dP/d(vsq) . - rho0 = D / gamma ; + rho0 = lg.D / gamma ; CCTK_REAL gm1 = gammaeos-1.; rho_gm1 = pow(rho0,gm1); p = (*sc_in) * rho_gm1 / gamma; @@ -340,20 +332,20 @@ // Make sure that W is large enough so that v^2 < 1 : i_increase = 0; - while( (( W_last*W_last*W_last * ( W_last + 2.*Bsq ) - - QdotBsq*(2.*W_last + Bsq) ) <= W_last*W_last*(Qtsq-Bsq*Bsq)) + while( (( W_last*W_last*W_last * ( W_last + 2.*lg.Bsq ) + - lg.QdotBsq*(2.*W_last + lg.Bsq) ) <= W_last*W_last*(lg.Qtsq-lg.Bsq*lg.Bsq)) && (i_increase < 10) ) { W_last *= 10.; i_increase++; } - W_for_gnr2 = W_for_gnr2_old = W_last; - rho_for_gnr2 = rho_for_gnr2_old = rho0; + lg.W_for_gnr2 = lg.W_for_gnr2_old = W_last; + lg.rho_for_gnr2 = lg.rho_for_gnr2_old = rho0; // Calculate W: x_1d[0] = W_last; - *retval = 1.0*oned_newton_raphson( x_1d, 1, gammaeos, func_W ) ; + *retval = 1.0*oned_newton_raphson( x_1d, 1, gammaeos, &lg, func_W ) ; W = x_1d[0]; @@ -369,15 +361,15 @@ } } - rho_g = x_rho[0] = rho_for_gnr2; + rho_g = x_rho[0] = lg.rho_for_gnr2; ntries = 0; - while ( (*retval = gnr2( x_rho, 1, gammaeos, func_rho)) && ( ntries++ < 10 ) ) { + while ( (*retval = gnr2( x_rho, 1, gammaeos, &lg, func_rho)) && ( ntries++ < 10 ) ) { rho_g *= 10.; x_rho[0] = rho_g; } - rho_for_gnr2 = x_rho[0]; + lg.rho_for_gnr2 = x_rho[0]; if( (*retval != 0) ) { *retval = 10; @@ -385,10 +377,10 @@ } // Calculate v^2 : - rho0 = rho_for_gnr2; + rho0 = lg.rho_for_gnr2; rho_gm1 = pow(rho0,gm1); - utsq = (D-rho0)*(D+rho0)/(rho0*rho0); + utsq = (lg.D-rho0)*(lg.D+rho0)/(rho0*rho0); gamma_sq = 1.+utsq; gamma = sqrt(gamma_sq); @@ -421,9 +413,9 @@ *w_lorentz = gamma; *pressure = p ; - g_o_WBsq = 1./(W+Bsq); + g_o_WBsq = 1./(W+lg.Bsq); QdB_o_W = QdotB / W; - *bsq = Bsq * (1.-vsq) + QdB_o_W*QdB_o_W; + *bsq = lg.Bsq * (1.-vsq) + QdB_o_W*QdB_o_W; *velx = g_o_WBsq * ( usx + QdB_o_W*(*Bx) ) ; *vely = g_o_WBsq * ( usy + QdB_o_W*(*By) ) ; @@ -449,49 +441,8 @@ } -/**********************************************************************/ -/**************************************************************************** - vsq_calc(): - - -- evaluate v^2 (spatial, normalized velocity) from - W = \gamma^2 w - -****************************************************************************/ -static CCTK_REAL vsq_calc(CCTK_REAL W) -{ - CCTK_REAL Wsq,Xsq,Bsq_W; - - Wsq = W*W ; - Bsq_W = (Bsq + W); - Xsq = Bsq_W * Bsq_W; - - return( ( Wsq * Qtsq + QdotBsq * (Bsq_W + W)) / (Wsq*Xsq) ); -} - - /******************************************************************** - x1_of_x0(): - - -- calculates v^2 from W with some physical bounds checking; - -- asumes x0 is already physical - -- makes v^2 physical if not; - -*********************************************************************/ - -static CCTK_REAL x1_of_x0(CCTK_REAL x0 ) -{ - CCTK_REAL x1,vsq; - CCTK_REAL dv = 1.e-15; - - vsq = fabs(vsq_calc(x0)) ; // guaranteed to be positive - - return( ( vsq > 1. ) ? (1.0 - dv) : vsq ); - -} - -/******************************************************************** - validate_x(): -- makes sure that x[0,1] have physical values, based upon @@ -524,10 +475,10 @@ -- inspired in part by Num. Rec.'s routine newt(); *****************************************************************/ -static CCTK_INT oned_newton_raphson( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos, +static CCTK_INT oned_newton_raphson( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos, struct LocGlob *lgp, void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][1], CCTK_REAL *, - CCTK_REAL *, CCTK_INT, CCTK_REAL) ) + CCTK_REAL *, CCTK_INT, CCTK_REAL, struct LocGlob *) ) { CCTK_REAL f, df, dx[1], x_old[1], resid[1], jac[1][1]; @@ -554,7 +505,7 @@ keep_iterating = 1; while( keep_iterating ) { - (*funcd) (x, dx, resid, jac, &f, &df, n, gammaeos); /* returns with new dx, f, df */ + (*funcd) (x, dx, resid, jac, &f, &df, n, gammaeos, lgp); /* returns with new dx, f, df */ /* Save old values before calculating the new: */ errx = 0.; @@ -568,8 +519,8 @@ // //METHOD specific: // i_increase = 0; -// while( (( x[0]*x[0]*x[0] * ( x[0] + 2.*Bsq ) - -// QdotBsq*(2.*x[0] + Bsq) ) <= x[0]*x[0]*(Qtsq-Bsq*Bsq)) +// while( (( x[0]*x[0]*x[0] * ( x[0] + 2.*lgp->Bsq ) - +// lgp->QdotBsq*(2.*x[0] + lgp->Bsq) ) <= x[0]*x[0]*(lgp->Qtsq-lgp->Bsq*lgp->Bsq)) // && (i_increase < 10) ) { // x[0] -= (1.*i_increase) * dx[0] / 10. ; // i_increase++; @@ -615,7 +566,7 @@ if( (!finite(f)) || (!finite(df)) || (!finite(x[0])) ) { #if(DEBUG_CON2PRIMM) fprintf(stderr,"\ngnr not finite, f,df,x_o,x,W_o,W,rho_o,rho = %26.20e %26.20e %26.20e %26.20e %26.20e %26.20e %26.20e %26.20e \n", - f,df,x[0],x_old[0],W_for_gnr2_old,W_for_gnr2,rho_for_gnr2_old,rho_for_gnr2); fflush(stderr); + f,df,x[0],x_old[0],lgp->W_for_gnr2_old,lgp->W_for_gnr2,lgp->rho_for_gnr2_old,lgp->rho_for_gnr2); fflush(stderr); #endif return(2); } @@ -652,52 +603,51 @@ n = dimension of x[]; *********************************************************************************/ static void func_W(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], - CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos) + CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos, + struct LocGlob *lgp) { CCTK_INT retval, ntries; CCTK_REAL W, x_rho[1], rho, rho_g ; - CCTK_REAL t15, t200,t1000 ; + CCTK_REAL t15, t200,t1000, drho_dW, rho_gm1 ; W = x[0]; - W_for_gnr2_old = W_for_gnr2; - W_for_gnr2 = W; + lgp->W_for_gnr2_old = lgp->W_for_gnr2; + lgp->W_for_gnr2 = W; // get rho from NR: - rho_g = x_rho[0] = rho_for_gnr2; + rho_g = x_rho[0] = lgp->rho_for_gnr2; ntries = 0; - while ( (retval = gnr2( x_rho, 1, gammaeos, func_rho)) && ( ntries++ < 10 ) ) { + while ( (retval = gnr2( x_rho, 1, gammaeos, lgp, func_rho)) && ( ntries++ < 10 ) ) { rho_g *= 10.; x_rho[0] = rho_g; } #if(DEBUG_CON2PRIMM) if( x_rho[0] <= 0. ) { - fprintf(stderr,"gnr2 neg rho = %d ,rho_n,rho,rho_o,W,W_o = %26.20e %26.20e %26.20e %26.20e %26.20e \n", retval, x_rho[0], rho_for_gnr2, rho_for_gnr2_old, x[0], W_for_gnr2_old); + fprintf(stderr,"gnr2 neg rho = %d ,rho_n,rho,rho_o,W,W_o = %26.20e %26.20e %26.20e %26.20e %26.20e \n", retval, x_rho[0], lgp->rho_for_gnr2, lgp->rho_for_gnr2_old, x[0], lgp->W_for_gnr2_old); fflush(stderr); } if( retval ) { - fprintf(stderr,"gnr2 retval = %d ,rho_n,rho,rho_o,W,W_o = %26.20e %26.20e %26.20e %26.20e %26.20e \n", retval, x_rho[0], rho_for_gnr2, rho_for_gnr2_old, x[0], W_for_gnr2_old); + fprintf(stderr,"gnr2 retval = %d ,rho_n,rho,rho_o,W,W_o = %26.20e %26.20e %26.20e %26.20e %26.20e \n", retval, x_rho[0], lgp->rho_for_gnr2, lgp->rho_for_gnr2_old, x[0], lgp->W_for_gnr2_old); fflush(stderr); } #endif - rho_for_gnr2_old = rho_for_gnr2; - rho = rho_for_gnr2 = x_rho[0]; + lgp->rho_for_gnr2_old = lgp->rho_for_gnr2; + rho = lgp->rho_for_gnr2 = x_rho[0]; CCTK_REAL gm1 = gammaeos-1.; rho_gm1 = pow(rho,gm1); - drho_dW = -rho*rho/( -rho_gm1*s200 + W*rho); + drho_dW = -rho*rho/( -rho_gm1*lgp->s200 + W*rho); - // t6 = rho*rho; - // t100 = (D-rho)*(D+rho); // t2 - t6; - t15 = -(D-rho)*(D+rho); // t6-t2 - t200 = W + two_Bsq; + t15 = -(lgp->D-rho)*(lgp->D+rho); // t6-t2 + t200 = W + lgp->two_Bsq; t1000 = rho*drho_dW; - resid[0] = (t300+(t4+t4+(t400+t15*(t7+(t200)*W))*W)*W)*t24; - jac[0][0] = 2*(t4+(t400+t15*t7+(3.0*t15*Bsq+t7*t1000+(t15+t15+t1000*(t200))*W)*W)*W)*t24; + resid[0] = (lgp->t300+(lgp->t4+lgp->t4+(lgp->t400+t15*(lgp->t7+(t200)*W))*W)*W)*lgp->t24; + jac[0][0] = 2*(lgp->t4+(lgp->t400+t15*lgp->t7+(3.0*t15*lgp->Bsq+lgp->t7*t1000+(t15+t15+t1000*(t200))*W)*W)*W)*lgp->t24; dx[0] = -resid[0]/jac[0][0]; @@ -705,11 +655,11 @@ *f = -0.5*(*df); -// fprintf(stdout,"QdotBsq = %28.18e ; \n",QdotBsq ); +// fprintf(stdout,"QdotBsq = %28.18e ; \n",lgp->QdotBsq ); // fprintf(stdout,"Sc = %28.18e ; \n",Sc ); -// fprintf(stdout,"Bsq = %28.18e ; \n",Bsq ); -// fprintf(stdout,"Qtsq = %28.18e ; \n",Qtsq ); -// fprintf(stdout,"Dc = %28.18e ; \n",D ); +// fprintf(stdout,"Bsq = %28.18e ; \n",lgp->Bsq ); +// fprintf(stdout,"Qtsq = %28.18e ; \n",lgp->Qtsq ); +// fprintf(stdout,"Dc = %28.18e ; \n",lgp->D ); // fprintf(stdout,"drhodW = %28.18e ; \n",drho_dW ); // fprintf(stdout,"W = %28.18e ; \n",W ); // fprintf(stdout,"rho = %28.18e ; \n",rho ); @@ -730,9 +680,9 @@ -- used to calculate rho from W *****************************************************************/ -static CCTK_INT gnr2( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos, +static CCTK_INT gnr2( CCTK_REAL x[], CCTK_INT n, CCTK_REAL gammaeos, struct LocGlob *lgp, void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], - CCTK_REAL [][1],CCTK_REAL *,CCTK_REAL *,CCTK_INT,CCTK_REAL) ) + CCTK_REAL [][1],CCTK_REAL *,CCTK_REAL *,CCTK_INT,CCTK_REAL, struct LocGlob *lgp) ) { CCTK_REAL f, df, dx[1], x_old[1], resid[1], jac[1][1]; @@ -757,7 +707,7 @@ keep_iterating = 1; while( keep_iterating ) { - (*funcd) (x, dx, resid, jac, &f, &df, n,gammaeos); /* returns with new dx, f, df */ + (*funcd) (x, dx, resid, jac, &f, &df, n,gammaeos, lgp); /* returns with new dx, f, df */ /* Save old values before calculating the new: */ //-fast errx = 0.; @@ -801,7 +751,7 @@ if( (!finite(f)) || (!finite(df)) || (!finite(x[0])) ) { #if(DEBUG_CON2PRIMM) fprintf(stderr,"\ngnr2 not finite, f,df,x_o,x,W_o,W,rho_o,rho = %26.20e %26.20e %26.20e %26.20e %26.20e %26.20e %26.20e %26.20e \n", - f,df,x[0],x_old[0],W_for_gnr2_old,W_for_gnr2,rho_for_gnr2_old,rho_for_gnr2); fflush(stderr); + f,df,x[0],x_old[0],lgp->W_for_gnr2_old,lgp->W_for_gnr2,lgp->rho_for_gnr2_old,lgp->rho_for_gnr2); fflush(stderr); #endif return(2); } @@ -839,7 +789,8 @@ *********************************************************************************/ // for the isentropic version: eq. (27) static void func_rho(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], - CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos) + CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, CCTK_INT n, CCTK_REAL gammaeos, + struct LocGlob *lgp) { CCTK_REAL A, B, C, rho, W, B0; @@ -848,13 +799,13 @@ CCTK_REAL gm1 = gammaeos-1.; rho = x[0]; - W = W_for_gnr2; + W = lgp->W_for_gnr2; - rho_gm1 = t40 = pow(rho,gm1); + t40 = pow(rho,gm1); - resid[0] = (rho*W+(-t40*s100-D)*D); + resid[0] = (rho*W+(-t40*lgp->s100-lgp->D)*lgp->D); t14 = t40/rho; // rho^(g-2) - jac[0][0] = -t14*s200 + W; + jac[0][0] = -t14*lgp->s200 + W; // drho_dW = -rho/jac[0][0]; dx[0] = -resid[0]/jac[0][0]; @@ -862,7 +813,7 @@ *f = -0.5*(*df); // fprintf(stdout,"deriv3 %g %g %g %g %g \n",rho,W,resid[0],jac[0][0],dx[0]); -// fprintf(stdout,"Dc := %28.18e ; \n",D); +// fprintf(stdout,"Dc := %28.18e ; \n",lgp->D); // fprintf(stdout,"Sc := %28.18e ; \n",Sc); // fprintf(stdout,"rho := %28.18e ; \n",rho); // fprintf(stdout,"W := %28.18e ; \n",W); @@ -874,40 +825,6 @@ } -/********************************************************************** - ********************************************************************** - - The following routines specify the equation of state. All routines - above here should be indpendent of EOS. If the user wishes - to use another equation of state, the below functions must be replaced - by equivalent routines based upon the new EOS. - - ********************************************************************** -**********************************************************************/ - -/**********************************************************************/ -/********************************************************************** - eos_info(): - - -- returns with all the EOS-related values needed; - **********************************************************************/ -static CCTK_REAL eos_info(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL *dpdw, CCTK_REAL *dpdvsq, CCTK_REAL gammaeos) -{ - register double ftmp,gtmp; - - ftmp = 1. - vsq; - gtmp = sqrt(ftmp); - - CCTK_REAL gam_m1_o_gam = (gammaeos-1.)/gammaeos; - - *dpdw = gam_m1_o_gam * ftmp ; - *dpdvsq = gam_m1_o_gam * ( 0.5 * D/gtmp - W ) ; - - return( gam_m1_o_gam * ( W * ftmp - D * gtmp ) ); // p - -} - - /****************************************************************************** END ******************************************************************************/ File [modified]: GRHydro_Con2PrimM_pt.c Delta lines: +45 -40 =================================================================== --- trunk/src/GRHydro_Con2PrimM_pt.c 2011-01-04 14:40:06 UTC (rev 207) +++ trunk/src/GRHydro_Con2PrimM_pt.c 2011-01-06 01:13:02 UTC (rev 208) @@ -71,23 +71,26 @@ ***************************************************/ /* Local Globals */ -static CCTK_REAL Bsq,QdotBsq,Qtsq,Qdotn,D,half_Bsq ; +struct LocGlob { + CCTK_REAL Bsq,QdotBsq,Qtsq,Qdotn,D,half_Bsq ; +} ; // Declarations: -static CCTK_REAL vsq_calc(CCTK_REAL W); +static CCTK_REAL vsq_calc(CCTK_REAL W, struct LocGlob *lgp); -static CCTK_INT twod_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, +static CCTK_INT twod_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, struct LocGlob *lgp, void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][2], - CCTK_REAL *, CCTK_REAL *, CCTK_REAL) ); + CCTK_REAL *, CCTK_REAL *, CCTK_REAL, struct LocGlob *) ); static void func_vsq( CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][2], - CCTK_REAL *f, CCTK_REAL *df, CCTK_REAL gammaeos); -static CCTK_REAL x1_of_x0(CCTK_REAL x0 ) ; + CCTK_REAL *f, CCTK_REAL *df, CCTK_REAL gammaeos, struct LocGlob *lgp); +static CCTK_REAL x1_of_x0(CCTK_REAL x0, struct LocGlob *lgp ) ; + // EOS STUFF: -static CCTK_REAL eos_info(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL *dpdw, CCTK_REAL *dpdvsq, CCTK_REAL gammaeos); +static CCTK_REAL eos_info(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL *dpdw, CCTK_REAL *dpdvsq, CCTK_REAL gammaeos, struct LocGlob *lgp); /* pressure as a function of rho0 and u */ static CCTK_REAL pressure_rho0_u(CCTK_REAL rho0, CCTK_REAL u, CCTK_REAL gammaeos) { @@ -213,6 +216,8 @@ CCTK_REAL sqrt_detg = sqrt(detg); CCTK_REAL inv_sqrt_detg = 1./sqrt_detg; CCTK_INT i,j, i_increase ; + + struct LocGlob lg; gammaeos = *gamma_eos; @@ -268,7 +273,7 @@ // Calculate various scalars (Q.B, Q^2, etc) from the conserved variables: - Bsq = + lg.Bsq = (*gxx) * (*Bx) * (*Bx) + (*gyy) * (*By) * (*By) + (*gzz) * (*Bz) * (*Bz) + @@ -278,15 +283,15 @@ (*gyz) * (*By) * (*Bz) ); QdotB = (sx * (*Bx) + sy * (*By) + sz * (*Bz)) ; - QdotBsq = QdotB*QdotB ; + lg.QdotBsq = QdotB*QdotB ; - Qdotn = -(tau + dens) ; + lg.Qdotn = -(tau + dens) ; - Qtsq = (usx * sx + usy * sy + usz * sz) ; + lg.Qtsq = (usx * sx + usy * sy + usz * sz) ; - D = dens; + lg.D = dens; - half_Bsq = 0.5*Bsq; + lg.half_Bsq = 0.5*lg.Bsq; /* calculate W from last timestep and use for guess */ vsq = @@ -311,7 +316,7 @@ // Always calculate rho from D and gamma so that using D in EOS remains consistent // i.e. you don't get positive values for dP/d(vsq) . - rho0 = D / gamma ; + rho0 = lg.D / gamma ; u = (*epsilon) * rho0; p = pressure_rho0_u(rho0,u,gammaeos) ; // EOS w = rho0 + u + p ; @@ -321,8 +326,8 @@ // Make sure that W is large enough so that v^2 < 1 : i_increase = 0; - while( (( W_last*W_last*W_last * ( W_last + 2.*Bsq ) - - QdotBsq*(2.*W_last + Bsq) ) <= W_last*W_last*(Qtsq-Bsq*Bsq)) + while( (( W_last*W_last*W_last * ( W_last + 2.*lg.Bsq ) + - lg.QdotBsq*(2.*W_last + lg.Bsq) ) <= W_last*W_last*(lg.Qtsq-lg.Bsq*lg.Bsq)) && (i_increase < 10) ) { W_last *= 10.; i_increase++; @@ -330,8 +335,8 @@ // Calculate W and vsq: x_2d[0] = fabs( W_last ); - x_2d[1] = x1_of_x0( W_last ) ; - *retval = 1.0*twod_newton_raphson( x_2d, gammaeos, func_vsq ) ; + x_2d[1] = x1_of_x0( W_last, &lg ) ; + *retval = 1.0*twod_newton_raphson( x_2d, gammaeos, &lg, func_vsq ) ; W = x_2d[0]; vsq = x_2d[1]; @@ -357,7 +362,7 @@ // Recover the primitive variables from the scalars and conserved variables: gtmp = sqrt(1. - vsq); gamma = 1./gtmp ; - rho0 = D * gtmp; + rho0 = lg.D * gtmp; w = W * (1. - vsq) ; p = pressure_rho0_w(rho0,w,gammaeos) ; // EOS @@ -375,9 +380,9 @@ *w_lorentz = gamma; *pressure = p ; - g_o_WBsq = 1./(W+Bsq); + g_o_WBsq = 1./(W+lg.Bsq); QdB_o_W = QdotB / W; - *bsq = Bsq * (1.-vsq) + QdB_o_W*QdB_o_W; + *bsq = lg.Bsq * (1.-vsq) + QdB_o_W*QdB_o_W; *velx = g_o_WBsq * ( usx + QdB_o_W*(*Bx) ) ; *vely = g_o_WBsq * ( usy + QdB_o_W*(*By) ) ; @@ -411,15 +416,15 @@ W = \gamma^2 w ****************************************************************************/ -static CCTK_REAL vsq_calc(CCTK_REAL W) +static CCTK_REAL vsq_calc(CCTK_REAL W, struct LocGlob *lgp) { CCTK_REAL Wsq,Xsq,Bsq_W; Wsq = W*W ; - Bsq_W = (Bsq + W); + Bsq_W = (lgp->Bsq + W); Xsq = Bsq_W * Bsq_W; - return( ( Wsq * Qtsq + QdotBsq * (Bsq_W + W)) / (Wsq*Xsq) ); + return( ( Wsq * lgp->Qtsq + lgp->QdotBsq * (Bsq_W + W)) / (Wsq*Xsq) ); } @@ -433,12 +438,12 @@ *********************************************************************/ -static CCTK_REAL x1_of_x0(CCTK_REAL x0 ) +static CCTK_REAL x1_of_x0(CCTK_REAL x0, struct LocGlob *lgp ) { CCTK_REAL x1,vsq; CCTK_REAL dv = 1.e-15; - vsq = fabs(vsq_calc(x0)) ; // guaranteed to be positive + vsq = fabs(vsq_calc(x0,lgp)) ; // guaranteed to be positive return( ( vsq > 1. ) ? (1.0 - dv) : vsq ); @@ -478,10 +483,10 @@ -- inspired in part by Num. Rec.'s routine newt(); *****************************************************************/ -static CCTK_INT twod_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, +static CCTK_INT twod_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, struct LocGlob *lgp, void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][2], CCTK_REAL *, - CCTK_REAL *, CCTK_REAL) ) + CCTK_REAL *, CCTK_REAL, struct LocGlob *) ) { CCTK_REAL f, df, dx[2], x_old[2]; CCTK_REAL resid[2], jac[2][2]; @@ -507,7 +512,7 @@ keep_iterating = 1; while( keep_iterating ) { - (*funcd) (x, dx, resid, jac, &f, &df, gammaeos); /* returns with new dx, f, df */ + (*funcd) (x, dx, resid, jac, &f, &df, gammaeos, lgp); /* returns with new dx, f, df */ /* Save old values before calculating the new: */ @@ -599,7 +604,7 @@ *********************************************************************************/ static void func_vsq(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], - CCTK_REAL jac[][2], CCTK_REAL *f, CCTK_REAL *df, CCTK_REAL gammaeos) + CCTK_REAL jac[][2], CCTK_REAL *f, CCTK_REAL *df, CCTK_REAL gammaeos, struct LocGlob *lgp) { @@ -612,22 +617,22 @@ Wsq = W*W; - p_tmp = eos_info(W, vsq, &dPdW, &dPdvsq, gammaeos); + p_tmp = eos_info(W, vsq, &dPdW, &dPdvsq, gammaeos, lgp); // These expressions were calculated using Mathematica, but made into efficient // code using Maple. Since we know the analytic form of the equations, we can // explicitly calculate the Newton-Raphson step: - t2 = -half_Bsq+dPdvsq; - t3 = Bsq+W; + t2 = -lgp->half_Bsq+dPdvsq; + t3 = lgp->Bsq+W; t4 = t3*t3; t23 = 1/W; - t16 = QdotBsq*t23*t23; - t11 = Qtsq-vsq*t4+t16*(Bsq+W+W); - t18 = -Qdotn-half_Bsq*(1.0+vsq)+0.5*t16-W+p_tmp; + t16 = lgp->QdotBsq*t23*t23; + t11 = lgp->Qtsq-vsq*t4+t16*(lgp->Bsq+W+W); + t18 = -lgp->Qdotn-lgp->half_Bsq*(1.0+vsq)+0.5*t16-W+p_tmp; t24 = t16*t23; t25 = -1.0+dPdW-t24; - t35 = t25*t3+(Bsq-2.0*dPdvsq)*(t16+vsq*W)*t23; + t35 = t25*t3+(lgp->Bsq-2.0*dPdvsq)*(t16+vsq*W)*t23; // t21 = 1/t3; // t36 = 1/t35; t21 = 1/(t3*t35); @@ -666,7 +671,7 @@ -- returns with all the EOS-related values needed; **********************************************************************/ -static CCTK_REAL eos_info(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL *dpdw, CCTK_REAL *dpdvsq, CCTK_REAL gammaeos) +static CCTK_REAL eos_info(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL *dpdw, CCTK_REAL *dpdvsq, CCTK_REAL gammaeos, struct LocGlob *lgp) { register double ftmp,gtmp; @@ -676,9 +681,9 @@ double gam_m1_o_gam = ((gammaeos-1.)/gammaeos); *dpdw = gam_m1_o_gam * ftmp ; - *dpdvsq = gam_m1_o_gam * ( 0.5 * D/gtmp - W ) ; + *dpdvsq = gam_m1_o_gam * ( 0.5 * lgp->D/gtmp - W ) ; - return( gam_m1_o_gam * ( W * ftmp - D * gtmp ) ); // p + return( gam_m1_o_gam * ( W * ftmp - lgp->D * gtmp ) ); // p }

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 207)
by cott＠tapir.caltech.edu 04 Jan '11

04 Jan '11

User: cott Date: 2011/01/04 08:40 AM Modified: /trunk/src/ GRHydro_Prim2Con.F90 Log: * Update error handling in _hot routine of Prim2Con. File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Prim2Con.F90 Delta lines: +25 -3 =================================================================== --- trunk/src/GRHydro_Prim2Con.F90 2011-01-02 15:54:48 UTC (rev 206) +++ trunk/src/GRHydro_Prim2Con.F90 2011-01-04 14:40:06 UTC (rev 207) @@ -231,9 +231,11 @@ ! begin EOS Omni vars integer :: n,keytemp,anyerr,keyerr(1) + real*8 :: temp0(1) n = 1;keytemp = 0;anyerr = 0;keyerr(1) = 0 ! end EOS Omni vars + temp0 = temp w = 1.d0 / sqrt(1.d0 - (gxx*dvelx*dvelx + gyy*dvely*dvely + gzz & *dvelz*dvelz + 2*gxy*dvelx*dvely + 2*gxz*dvelx *dvelz + 2*gyz& *dvely*dvelz)) @@ -262,19 +264,39 @@ call CCTK_WARN(1,warnline) endif else - call CCTK_WARN(1,"EOS error in prim2con_hot") + ! This is a way of recovering even on finer refinement levels: + ! Use the average temperature at the interface instead of the + ! reconstructed specific internal energy. + call CCTK_WARN(1,"EOS error in prim2con_hot: NOW using averaged temp!") write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z call CCTK_WARN(1,warnline) - write(warnline,"(1P10E15.6)") drho,deps,temp,ye + write(warnline,"(1P10E15.6)") drho,deps,temp0,ye call CCTK_WARN(1,warnline) write(warnline,"(A7,i8)") "code: ",keyerr(1) call CCTK_WARN(1,warnline) write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel call CCTK_WARN(1,warnline) - call CCTK_WARN(0,"Aborting!!!") + keytemp=1 + temp = temp0 + call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& + drho,deps,temp,ye,dpress,keyerr,anyerr) + keytemp=0 + if(anyerr.ne.0) then + call CCTK_WARN(1,"EOS error in prim2con_hot") + write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z + call CCTK_WARN(1,warnline) + write(warnline,"(1P10E15.6)") drho,deps,temp,ye + call CCTK_WARN(1,warnline) + write(warnline,"(A7,i8)") "code: ",keyerr(1) + call CCTK_WARN(1,warnline) + write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel + call CCTK_WARN(1,warnline) + call CCTK_WARN(0,"Aborting!!!") + endif endif endif + vlowx = gxx*dvelx + gxy*dvely + gxz*dvelz vlowy = gxy*dvelx + gyy*dvely + gyz*dvelz vlowz = gxz*dvelx + gyz*dvely + gzz*dvelz

1 0

2026

2025

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

Commits January 2011