From d14686b8cb54a26334a7a66a05f5d7910203e966 Mon Sep 17 00:00:00 2001 From: rhaas Date: Thu, 28 Mar 2013 01:47:14 +0000 Subject: Finish the implementation of con2prim using the entropy equation. * Tested successfully pointwisely. It still needs to be tested on evolution. From: Bruno Coutinho Mundim git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@500 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45 --- interface.ccl | 22 ++ src/GRHydro_Con2PrimM.F90 | 44 +++- src/GRHydro_Con2PrimM_ptee.c | 592 +++++++++++++++++++++++++++++++++++++++++++ src/make.code.defn | 1 + 4 files changed, 649 insertions(+), 10 deletions(-) create mode 100644 src/GRHydro_Con2PrimM_ptee.c diff --git a/interface.ccl b/interface.ccl index 89cae27..db6615d 100644 --- a/interface.ccl +++ b/interface.ccl @@ -84,6 +84,27 @@ void FUNCTION Con2PrimGenM(CCTK_INT INOUT handle, CCTK_INT INOUT keytemp, CCTK_R CCTK_INT OUT epsnegative, \ CCTK_REAL OUT retval) +void FUNCTION Con2PrimGenMee(CCTK_INT INOUT handle, CCTK_INT INOUT keytemp, \ + CCTK_REAL INOUT prec, CCTK_REAL INOUT gamma_eos, \ + CCTK_REAL INOUT dens, \ + CCTK_REAL INOUT sx, CCTK_REAL INOUT sy, CCTK_REAL INOUT sz, \ + CCTK_REAL INOUT tau, \ + CCTK_REAL INOUT Bconsx, CCTK_REAL INOUT Bconsy, CCTK_REAL INOUT Bconsz, \ + CCTK_REAL INOUT entropycons, \ + CCTK_REAL INOUT y_e, CCTK_REAL INOUT temp, CCTK_REAL INOUT rho, \ + CCTK_REAL INOUT velx, CCTK_REAL INOUT vely, CCTK_REAL INOUT velz, \ + CCTK_REAL INOUT epsilon, CCTK_REAL INOUT pressure, \ + CCTK_REAL INOUT Bvecx, CCTK_REAL INOUT Bvecy, CCTK_REAL INOUT Bvecz, \ + CCTK_REAL INOUT Bvecsq, \ + CCTK_REAL INOUT w_lorentz, \ + CCTK_REAL INOUT gxx, CCTK_REAL INOUT gxy, CCTK_REAL INOUT gxz, \ + CCTK_REAL INOUT gyy, CCTK_REAL INOUT gyz, CCTK_REAL INOUT gzz, \ + CCTK_REAL INOUT uxx, CCTK_REAL INOUT uxy, CCTK_REAL INOUT uxz, \ + CCTK_REAL INOUT uyy, CCTK_REAL INOUT uyz, CCTK_REAL INOUT uzz, \ + CCTK_REAL INOUT det, \ + CCTK_INT OUT epsnegative, \ + CCTK_REAL OUT retval) + #void FUNCTION Con2PrimGenM(CCTK_INT INOUT handle, CCTK_INT INOUT keytemp, CCTK_REAL INOUT prec,CCTK_REAL INOUT gamma_eos, CCTK_REAL INOUT dens, \ # CCTK_REAL INOUT sx, CCTK_REAL INOUT sy, CCTK_REAL INOUT sz, \ # CCTK_REAL INOUT tau, CCTK_REAL INOUT Bconsx, CCTK_REAL INOUT Bconsy, CCTK_REAL INOUT Bconsz, \ @@ -207,6 +228,7 @@ PROVIDES FUNCTION UpperMet WITH UpperMetric LANGUAGE Fortran #PROVIDES FUNCTION Con2Prim WITH Con2Prim_pt LANGUAGE Fortran PROVIDES FUNCTION Con2PrimPoly WITH Con2Prim_ptPolytype LANGUAGE Fortran PROVIDES FUNCTION Con2PrimGenM WITH GRHydro_Con2PrimM_pt LANGUAGE Fortran +PROVIDES FUNCTION Con2PrimGenMee WITH GRHydro_Con2PrimM_ptee LANGUAGE Fortran PROVIDES FUNCTION Con2PrimGen WITH Con2Prim_pt LANGUAGE Fortran PROVIDES FUNCTION Con2PrimPolyM WITH GRHydro_Con2PrimM_Polytype_pt LANGUAGE Fortran PROVIDES FUNCTION Prim2ConGen WITH prim2con LANGUAGE Fortran diff --git a/src/GRHydro_Con2PrimM.F90 b/src/GRHydro_Con2PrimM.F90 index d59cc9c..0960b0d 100644 --- a/src/GRHydro_Con2PrimM.F90 +++ b/src/GRHydro_Con2PrimM.F90 @@ -334,7 +334,11 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS) Bvecz_tmp = Bprim(i,j,k,3) keytemp = 0 - + !Watch out for the values returned to b2. Here b2 is the Bprim^2 + !while inside the point-wise con2prim routines it is the square + !of the comoving B-field, b^{\mu} b_{\mu}. It is overwritten + !in this routine, but we may need to find a better notation + !avoid future confusions. call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, keytemp, & GRHydro_eos_rf_prec, local_gam(1), dens(i,j,k), & scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), tau(i,j,k), & @@ -347,7 +351,13 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS) epsnegative,GRHydro_C2P_failed(i,j,k)) if(evolve_entropy.ne.0) then - entropy(i,j,k) = entropycons(i,j,k)/dens(i,j,k)*rho(i,j,k) + if(GRHydro_C2P_failed(i,j,k).ne.0) then + !Use previous time step for rho: + entropy(i,j,k) = entropycons(i,j,k)/dens(i,j,k)*rho(i,j,k) + else + !Use the current correct value of rho returned by con2prim: + entropy(i,j,k) = entropycons(i,j,k)/dens(i,j,k)*rho_tmp + endif endif if(GRHydro_C2P_failed(i,j,k).ne.0) then @@ -430,14 +440,28 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS) Bvecy_tmp = Bprim(i,j,k,2) Bvecz_tmp = Bprim(i,j,k,3) - call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, & - dens(i,j,k),scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), sc, & - Bcons(i,j,k,1), Bcons(i,j,k,2), Bcons(i,j,k,3),rho_tmp,& - velx_tmp,vely_tmp,velz_tmp,eps_tmp,press_tmp,& - Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,w_lorentz_tmp,& - g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k), & - uxx,uxy,uxz,uyy,uyz,uzz,det, & - epsnegative,GRHydro_C2P_failed(i,j,k)) + if(evolve_entropy.ne.0) then + call GRHydro_Con2PrimM_ptee(GRHydro_eos_handle, keytemp, & + GRHydro_eos_rf_prec, local_gam(1), dens(i,j,k), & + scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), tau(i,j,k), & + Bcons(i,j,k,1),Bcons(i,j,k,2),Bcons(i,j,k,3), & + entropycons(i,j,k), xye(1), & + xtemp(1),rho_tmp,velx_tmp,vely_tmp,velz_tmp,& + eps_tmp,press_tmp,Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,& + w_lorentz_tmp,g11(i,j,k),g12(i,j,k),g13(i,j,k),& + g22(i,j,k),g23(i,j,k),g33(i,j,k), & + uxx,uxy,uxz,uyy,uyz,uzz,det, & + epsnegative,GRHydro_C2P_failed(i,j,k)) + else + call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, & + dens(i,j,k),scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), sc, & + Bcons(i,j,k,1), Bcons(i,j,k,2), Bcons(i,j,k,3),rho_tmp,& + velx_tmp,vely_tmp,velz_tmp,eps_tmp,press_tmp,& + Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,w_lorentz_tmp,& + g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k), & + uxx,uxy,uxz,uyy,uyz,uzz,det, & + epsnegative,GRHydro_C2P_failed(i,j,k)) + end if rho(i,j,k) = rho_tmp press(i,j,k) = press_tmp diff --git a/src/GRHydro_Con2PrimM_ptee.c b/src/GRHydro_Con2PrimM_ptee.c new file mode 100644 index 0000000..de58d4c --- /dev/null +++ b/src/GRHydro_Con2PrimM_ptee.c @@ -0,0 +1,592 @@ +/*********************************************************************************** + Copyright 2006 Scott C. Noble, Charles F. Gammie, Jonathan C. McKinney, + and Luca Del Zanna. + + PVS_GRMHD + + This file was derived from PVS_GRMHD. The authors of PVS_GRMHD include + Scott C. Noble, Charles F. Gammie, Jonathan C. McKinney, and Luca Del Zanna. + PVS_GRMHD is available under the GPL from: + http://rainman.astro.uiuc.edu/codelib/ + + You are morally obligated to cite the following paper in his/her + scientific literature that results from use of this file: + + [1] Noble, S. C., Gammie, C. F., McKinney, J. C., \& Del Zanna, L. \ 2006, + Astrophysical Journal, 641, 626. + + PVS_GRMHD 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. + + PVS_GRMHD 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 PVS_GRMHD; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + + If the user has any questions, please direct them to Scott C. Noble at + scn@astro.rit.edu . + +***********************************************************************************/ + + + +#include +#include +#include +#include +#include +#include +#include + +#include "cctk.h" +#include "cctk_Parameters.h" + +/* Set this to be 1 if you want debug output */ +#define DEBUG_CON2PRIMM (0) + + +/* Adiabatic index used for the state equation */ + +#define MAX_NEWT_ITER (30) /* Max. # of Newton-Raphson iterations for find_root_2D(); */ +#define NEWT_TOL (1.0e-10) /* Min. of tolerance allowed for Newton-Raphson iterations */ +#define MIN_NEWT_TOL (1.0e-10) /* Max. of tolerance allowed for Newton-Raphson iterations */ +#define EXTRA_NEWT_ITER (2) + +#define NEWT_TOL2 (1.0e-15) /* TOL of new 1D^*_{v^2} gnr2 method */ +#define MIN_NEWT_TOL2 (1.0e-10) /* TOL of new 1D^*_{v^2} gnr2 method */ + +#define W_TOO_BIG (1.e20) /* \gamma^2 (\rho_0 + u + p) is assumed + to always be smaller than this. This + is used to detect solver failures */ + +#define FAIL_VAL (1.e30) /* Generic value to which we set variables when a problem arises */ + +/************************************************** + The following functions assume a Gamma-law EOS: +***************************************************/ + +/* Local Globals */ +struct LocGlob { + CCTK_REAL Bsq, QdotBsq, Qtsq, Qdotn, D, half_Bsq, Sc, g_o_gm1, + W_for_gnr2, rho_for_gnr2, W_for_gnr2_old, rho_for_gnr2_old, drho_dW ; +} ; + +// Declarations: + + +void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_ptee) ( + CCTK_INT *handle, CCTK_INT *keytemp, CCTK_REAL *prec, + CCTK_REAL *gamma_eos, + CCTK_REAL *dens_in, + CCTK_REAL *sx_in, CCTK_REAL *sy_in, CCTK_REAL *sz_in, + CCTK_REAL *tau_in, + CCTK_REAL *Bconsx_in, CCTK_REAL *Bconsy_in, CCTK_REAL *Bconsz_in, + CCTK_REAL *entropycons_in, + CCTK_REAL *y_e_in, CCTK_REAL* temp_in, + CCTK_REAL *rho, + CCTK_REAL *velx, CCTK_REAL *vely, CCTK_REAL *velz, + CCTK_REAL *epsilon, CCTK_REAL *pressure, + CCTK_REAL *Bx, CCTK_REAL *By, CCTK_REAL *Bz, + CCTK_REAL *bsq, + CCTK_REAL *w_lorentz, + CCTK_REAL *gxx, CCTK_REAL *gxy, CCTK_REAL *gxz, + CCTK_REAL *gyy, CCTK_REAL *gyz, CCTK_REAL *gzz, + CCTK_REAL *uxx, CCTK_REAL *uxy, CCTK_REAL *uxz, + CCTK_REAL *uyy, CCTK_REAL *uyz, CCTK_REAL *uzz, + CCTK_REAL *det, + CCTK_INT *epsnegative, + CCTK_REAL *retval) ; + +static CCTK_INT general_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, + struct LocGlob *lgp, + void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], + CCTK_REAL [][1], CCTK_REAL *, + CCTK_REAL *, CCTK_REAL, struct LocGlob *) ); + +static void func_rho(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], + CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, + CCTK_REAL gammaeos, struct LocGlob *lgp); + +/**********************************************************************/ +/********************************************************************************** + + Con2PrimM_pt(): + ----------------------------- + + -- Attempts an inversion from GRMHD conserved variables to primitive variables assuming a guess. + + -- Uses the 2D method of Noble et al. (2006): + -- Solves for two independent variables (W,v^2) via a 2D + Newton-Raphson method + -- Can be used (in principle) with a general equation of state. + + -- Minimizes two residual functions using a homemade Newton-Raphson routine. + -- It is homemade so that it can catch exceptions and handle them correctly, plus it is + optimized for this particular problem. + + -- Note that the notation used herein is that of Noble et al. (2006) except for the argument + list. + + +INPUT: (using GRHydro variable defintions) + + s[x,y,z] = scons[0,1,2] = \alpha \sqrt(\gamma) T^0_i + dens, tau = as defined in GRHydro and are assumed to be densitized (i.e. with sqrt(\gamma)) + dens = D = \sqrt(\gamma) W \rho + tau = \alpha^2 \sqrt(\gamma) T^{00} - D + g[x,y,z][x,y,x] = spatial metric corresponding to \gamma + u[x,y,z][x,y,z] = inverse of the spatial metric, g[x,y,z][x,y,x] + det = sqrt(\gamma) + B[x,y,z] = Bvec[0,1,2] + bsq = b^\mu b_\mu + + epsnegative = (integer) + = 0 if rho and epsilon are positive + != 0 otherwise + + + -- (users should set B[x,y,z] = 0 for hydrodynamic runs) + + +OUTPUT: (using GRHydro variable defintions) + rho, eps = as defined in GRHydro, primitive variables + vel[x,y,z] = as defined in GRHydro, primitive variables + + +RETURN VALUE: of retval = (i*100 + j) where + i = 0 -> Newton-Raphson solver either was not called (yet or not used) + or returned successfully; + 1 -> Newton-Raphson solver did not converge to a solution with the + given tolerances; + 2 -> Newton-Raphson procedure encountered a numerical divergence + (occurrence of "nan" or "+/-inf" ; + + j = 0 -> success + 1 -> failure: some sort of failure in Newton-Raphson; + 2 -> failure: unphysical vsq = v^2 value at initial guess; + 3 -> failure: W<0 or W>W_TOO_BIG + 4 -> failure: v^2 > 1 + ( used to be 5 -> failure: rho,uu <= 0 but now sets epsnegative to non-zero ) + +**********************************************************************************/ + +void CCTK_FCALL CCTK_FNAME(GRHydro_Con2PrimM_ptee) ( + CCTK_INT *handle, CCTK_INT *keytemp, CCTK_REAL *prec, + CCTK_REAL *gamma_eos, + CCTK_REAL *dens_in, + CCTK_REAL *sx_in, CCTK_REAL *sy_in, CCTK_REAL *sz_in, + CCTK_REAL *tau_in, + CCTK_REAL *Bconsx_in, CCTK_REAL *Bconsy_in, CCTK_REAL *Bconsz_in, + CCTK_REAL *entropycons_in, + CCTK_REAL *y_e_in, CCTK_REAL* temp_in, + CCTK_REAL *rho, + CCTK_REAL *velx, CCTK_REAL *vely, CCTK_REAL *velz, + CCTK_REAL *epsilon, CCTK_REAL *pressure, + CCTK_REAL *Bx, CCTK_REAL *By, CCTK_REAL *Bz, + CCTK_REAL *bsq, + CCTK_REAL *w_lorentz, + CCTK_REAL *gxx, CCTK_REAL *gxy, CCTK_REAL *gxz, + CCTK_REAL *gyy, CCTK_REAL *gyz, CCTK_REAL *gzz, + CCTK_REAL *uxx, CCTK_REAL *uxy, CCTK_REAL *uxz, + CCTK_REAL *uyy, CCTK_REAL *uyz, CCTK_REAL *uzz, + CCTK_REAL *det, + CCTK_INT *epsnegative, + CCTK_REAL *retval) + +{ + CCTK_REAL x_1d[1]; + CCTK_REAL sx, sy, sz; + CCTK_REAL usx, usy, usz; + CCTK_REAL tau, dens, gammaeos; + CCTK_REAL QdotB; + CCTK_REAL rho0,u,p,w,gammasq,gamma,W,vsq; + CCTK_REAL g_o_WBsq, QdB_o_W; + CCTK_REAL detg = (*det); + CCTK_REAL sqrt_detg = sqrt(detg); + CCTK_REAL inv_sqrt_detg = 1./sqrt_detg; + CCTK_REAL rho_gm1; + CCTK_REAL utsq; + + DECLARE_CCTK_PARAMETERS; + + struct LocGlob lg; + + gammaeos = *gamma_eos; + + /* Assume ok initially: */ + // BCM: But let the driver function take care of its initialization + //*retval = 0.; + *epsnegative = 0; + +#if(DEBUG_CON2PRIMM) + fprintf(stdout," *dens = %26.16e \n", *dens_in ); + fprintf(stdout," *sx = %26.16e \n", *sx_in ); + fprintf(stdout," *sy = %26.16e \n", *sy_in ); + fprintf(stdout," *sz = %26.16e \n", *sz_in ); + fprintf(stdout," *tau = %26.16e \n", *tau_in ); + fprintf(stdout," *Bconsx = %26.16e \n", *Bconsx_in ); + fprintf(stdout," *Bconsy = %26.16e \n", *Bconsy_in ); + fprintf(stdout," *Bconsz = %26.16e \n", *Bconsz_in ); + fprintf(stdout," *entropycons = %26.16e \n", *entropycons_in ); + fprintf(stdout," *rho = %26.16e \n", *rho ); + fprintf(stdout," *velx = %26.16e \n", *velx ); + fprintf(stdout," *vely = %26.16e \n", *vely ); + fprintf(stdout," *velz = %26.16e \n", *velz ); + fprintf(stdout," *epsilon = %26.16e \n", *epsilon ); + fprintf(stdout," *pressure = %26.16e \n", *pressure ); + fprintf(stdout," *Bx = %26.16e \n", *Bx ); + fprintf(stdout," *By = %26.16e \n", *By ); + fprintf(stdout," *Bz = %26.16e \n", *Bz ); + fprintf(stdout," *bsq = %26.16e \n", *bsq ); + fprintf(stdout," *w_lorentz = %26.16e \n", *w_lorentz ); + fprintf(stdout," *gxx = %26.16e \n", *gxx ); + fprintf(stdout," *gxy = %26.16e \n", *gxy ); + fprintf(stdout," *gxz = %26.16e \n", *gxz ); + fprintf(stdout," *gyy = %26.16e \n", *gyy ); + fprintf(stdout," *gyz = %26.16e \n", *gyz ); + fprintf(stdout," *gzz = %26.16e \n", *gzz ); + fprintf(stdout," *uxx = %26.16e \n", *uxx ); + fprintf(stdout," *uxy = %26.16e \n", *uxy ); + fprintf(stdout," *uxz = %26.16e \n", *uxz ); + fprintf(stdout," *uyy = %26.16e \n", *uyy ); + fprintf(stdout," *uyz = %26.16e \n", *uyz ); + fprintf(stdout," *uzz = %26.16e \n", *uzz ); + fprintf(stdout," *det = %26.16e \n", *det ); + fprintf(stdout," *epsnegative = %10d \n", *epsnegative ); + fprintf(stdout," *retval = %26.16e \n", *retval ); + fflush(stdout); +#endif + + /* First undensitize all conserved variables : */ + sx = ( *sx_in) * inv_sqrt_detg; + sy = ( *sy_in) * inv_sqrt_detg; + sz = ( *sz_in) * inv_sqrt_detg; + tau = ( *tau_in) * inv_sqrt_detg; + dens = (*dens_in) * inv_sqrt_detg; + + usx = (*uxx)*sx + (*uxy)*sy + (*uxz)*sz; + usy = (*uxy)*sx + (*uyy)*sy + (*uyz)*sz; + usz = (*uxz)*sx + (*uyz)*sy + (*uzz)*sz; + + *Bx = (*Bconsx_in) * inv_sqrt_detg; + *By = (*Bconsy_in) * inv_sqrt_detg; + *Bz = (*Bconsz_in) * inv_sqrt_detg; + + lg.Sc = (*entropycons_in) * inv_sqrt_detg; + // Calculate various scalars (Q.B, Q^2, etc) from the conserved variables: + + lg.g_o_gm1 = gammaeos/(gammaeos-1.0); + + lg.Bsq = + (*gxx) * (*Bx) * (*Bx) + + (*gyy) * (*By) * (*By) + + (*gzz) * (*Bz) * (*Bz) + + 2*( + (*gxy) * (*Bx) * (*By) + + (*gxz) * (*Bx) * (*Bz) + + (*gyz) * (*By) * (*Bz) ); + + QdotB = (sx * (*Bx) + sy * (*By) + sz * (*Bz)) ; + lg.QdotBsq = QdotB*QdotB ; + + lg.Qdotn = -(tau + dens) ; + + lg.Qtsq = (usx * sx + usy * sy + usz * sz) ; + + lg.D = dens; + + lg.half_Bsq = 0.5*lg.Bsq; + + /* calculate W from last timestep and use for guess */ + vsq = + (*gxx) * (*velx) * (*velx) + + (*gyy) * (*vely) * (*vely) + + (*gzz) * (*velz) * (*velz) + + 2*( + (*gxy) * (*velx) * (*vely) + + (*gxz) * (*velx) * (*velz) + + (*gyz) * (*vely) * (*velz) ); + + if( (vsq < 0.) && (fabs(vsq) < 1.0e-13) ) { + vsq = fabs(vsq); + } + if(vsq < 0. || vsq > 1. ) { + *retval = 2.; + return; + } + + gammasq = 1. / (1. - vsq); + gamma = sqrt(gammasq); + + // 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 = lg.D / gamma ; + rho_gm1 = pow(rho0,(gammaeos-1.)); + p = lg.Sc * rho_gm1 / gamma; + u = p / (gammaeos-1.); + w = rho0 + u + p ; + +// W_last = w*gammasq ; + + // Calculate W and vsq: + x_1d[0] = rho0; +// *retval = 1.0*twod_newton_raphson( x_2d, gammaeos, &lg, func_vsq ) ; + *retval = general_newton_raphson( x_1d, gammaeos, &lg, func_rho ) ; + rho0 = x_1d[0]; + + /* Problem with solver, so return denoting error before doing anything further */ + if( ((*retval) != 0.) || (rho0 == FAIL_VAL) ) { + *retval = *retval*100.+1.; + return; + } + else{ + if( rho0 > W_TOO_BIG) { + *retval = 3.; + return; + } + } + + // Calculate v^2: + utsq = (lg.D-rho0)*(lg.D+rho0)/(rho0*rho0); + gammasq = 1.+utsq; + gamma = sqrt(gammasq); + + if( utsq < 0. ) { + *retval = 4.; + return; + } + + + // Recover the primitive variables from the scalars and conserved variables: + rho0 = lg.D / gamma; + rho_gm1 = pow(rho0,(gammaeos-1.)); + p = lg.Sc * rho_gm1 / gamma; + u = p / (gammaeos-1.); + w = rho0 + u + p ; + W = w * gammasq; + + // User may want to handle this case differently, e.g. do NOT return upon + // a negative rho/u, calculate v^i so that rho/u can be floored by other routine: + if( (rho0 <= 0.) || (u <= 0.) ) { + *epsnegative = 1; + return; + } + + *rho = rho0; + *epsilon = u / rho0; + *w_lorentz = gamma; + *pressure = p ; + + g_o_WBsq = 1./(W+lg.Bsq); + QdB_o_W = QdotB / W; + *bsq = lg.Bsq / gammasq + QdB_o_W*QdB_o_W; + + *velx = g_o_WBsq * ( usx + QdB_o_W*(*Bx) ) ; + *vely = g_o_WBsq * ( usy + QdB_o_W*(*By) ) ; + *velz = g_o_WBsq * ( usz + QdB_o_W*(*Bz) ) ; + + if (*rho <= rho_abs_min*(1.0+GRHydro_atmo_tolerance) ) { + *rho = rho_abs_min; + *velx = 0.0; + *vely = 0.0; + *velz = 0.0; + *w_lorentz = 1.0; + } + + +#if(DEBUG_CON2PRIMM) + fprintf(stdout,"rho = %26.16e \n",*rho ); + fprintf(stdout,"epsilon = %26.16e \n",*epsilon ); + fprintf(stdout,"pressure = %26.16e \n",*pressure ); + fprintf(stdout,"w_lorentz = %26.16e \n",*w_lorentz); + fprintf(stdout,"bsq = %26.16e \n",*bsq ); + fprintf(stdout,"velx = %26.16e \n",*velx ); + fprintf(stdout,"vely = %26.16e \n",*vely ); + fprintf(stdout,"velz = %26.16e \n",*velz ); + fprintf(stdout,"gammaeos = %26.16e \n",gammaeos ); + fflush(stdout); +#endif + + /* done! */ + return; + +} + +/**********************************************************************/ +/************************************************************ + + general_newton_raphson(): + + -- performs Newton-Rapshon method on an arbitrary system. + + -- inspired in part by Num. Rec.'s routine newt(); + + Arguements: + + -- x[] = set of independent variables to solve for; + -- n = number of independent variables and residuals; + -- funcd = name of function that calculates residuals, etc.; + +*****************************************************************/ +static CCTK_INT general_newton_raphson( CCTK_REAL x[], CCTK_REAL gammaeos, + struct LocGlob *lgp, + void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], + CCTK_REAL [][1], CCTK_REAL *, + CCTK_REAL *, CCTK_REAL, struct LocGlob *) ) +{ + CCTK_REAL f, df, dx[1], x_old[1], resid[1], + jac[1][1]; + CCTK_REAL errx, x_orig[1]; + CCTK_INT n_iter, i_extra, doing_extra; + CCTK_REAL W,W_old; + + int keep_iterating; + + + // Initialize various parameters and variables: + errx = 1. ; + df = f = 1.; + i_extra = doing_extra = 0; + x_old[0] = x_orig[0] = x[0] ; + + W = W_old = 0.; + + n_iter = 0; + + + /* Start the Newton-Raphson iterations : */ + keep_iterating = 1; + while( keep_iterating ) { + + (*funcd) (x, dx, resid, jac, &f, &df, gammaeos, lgp); /* returns with new dx, f, df */ + + /* Save old values before calculating the new: */ + errx = 0.; + x_old[0] = x[0] ; + /* don't use line search : */ + x[0] += dx[0] ; + + /****************************************/ + /* Calculate the convergence criterion */ + /****************************************/ + + /* For the new criterion, always look at error in "W" : */ + // METHOD specific: + errx = (x[0]==0.) ? fabs(dx[0]) : fabs(dx[0]/x[0]); + + + /****************************************/ + /* Make sure that the new x[] is physical : */ + /****************************************/ + x[0] = fabs(x[0]); + + + /*****************************************************************************/ + /* If we've reached the tolerance level, then just do a few extra iterations */ + /* before stopping */ + /*****************************************************************************/ + + if( (fabs(errx) <= NEWT_TOL) && (doing_extra == 0) && (EXTRA_NEWT_ITER > 0) ) { + doing_extra = 1; + } + + if( doing_extra == 1 ) i_extra++ ; + + if( ((fabs(errx) <= NEWT_TOL)&&(doing_extra == 0)) || + (i_extra > EXTRA_NEWT_ITER) || (n_iter >= (MAX_NEWT_ITER-1)) ) { + keep_iterating = 0; + } + + n_iter++; + + } // END of while(keep_iterating) + + + /* Check for bad untrapped divergences : */ + if( (!finite(f)) || (!finite(df)) || (!finite(x[0])) ) { + return(2); + } + + + if( fabs(errx) > MIN_NEWT_TOL){ + return(1); + } + if( (fabs(errx) <= MIN_NEWT_TOL) && (fabs(errx) > NEWT_TOL) ){ + return(0); + } + if( fabs(errx) <= NEWT_TOL ){ + return(0); + } + + return(0); + +} + +/**********************************************************************/ +/********************************************************************************* + func_rho(): + + -- residual/jacobian routine to calculate rho Qtsq equation with + the definition of W + W = ( 1 + GAMMA * K_atm * rho^(GAMMA-1)/(GAMMA-1) ) D^2 / rho + substituted in. + + Arguments: + x = current value of independent var's (on input & output); + dx = Newton-Raphson step (on output); + resid = residuals based on x (on output); + jac = Jacobian matrix based on x (on output); + f = resid.resid/2 (on output) + df = -2*f; (on output) + n = dimension of x[]; + *********************************************************************************/ +static void func_rho(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], + CCTK_REAL jac[][1], CCTK_REAL *f, CCTK_REAL *df, + CCTK_REAL gammaeos, struct LocGlob *lgp) +{ + + CCTK_REAL t1,t2; + CCTK_REAL t12; + CCTK_REAL t17; + CCTK_REAL t3, t100,rhosq,t200,rho,W,dWdrho,dvsqdrho,vsq; + + rho = x[0]; + rhosq = rho*rho; + t200 = 1./(lgp->D*lgp->D); + t100 = lgp->g_o_gm1*lgp->Sc*pow(rho,(gammaeos-1.)); + W = lgp->D*( lgp->D + t100 ) / rho ; + dWdrho = lgp->D * ( -lgp->D + t100*(gammaeos-2.) ) / rhosq; + t1 = W*W; + t2 = lgp->Bsq+W; + // t3 = pow(Bsq+W,2.0); + t3 = t2*t2; + vsq = (lgp->D-rho)*(lgp->D+rho)*t200; + dvsqdrho = -2*rho*t200; + resid[0] = t1*(lgp->Qtsq-vsq*t3)+lgp->QdotBsq*(t2+W); + t12 = lgp->Bsq*lgp->Bsq; + t17 = dWdrho*vsq; + jac[0][0] = 2*lgp->QdotBsq*dWdrho + +((lgp->Qtsq-vsq*t12)*2*dWdrho+(-6*t17*lgp->Bsq-dvsqdrho*t12 + +(-2*dvsqdrho*lgp->Bsq-4*t17-dvsqdrho*W)*W)*W)*W; + + dx[0] = -resid[0]/jac[0][0]; + *f = 0.5*resid[0]*resid[0]; + *df = -2. * (*f); + + return; + +} + + +/****************************************************************************** + END + ******************************************************************************/ + + +#undef DEBUG_CON2PRIMM diff --git a/src/make.code.defn b/src/make.code.defn index c16764e..7bdd9f6 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -46,6 +46,7 @@ SRCS = Utils.F90 \ GRHydro_Tmunu.F90 \ GRHydro_Con2PrimM.F90 \ GRHydro_Con2PrimM_pt.c \ + GRHydro_Con2PrimM_ptee.c \ GRHydro_Con2PrimM_polytype_pt.c \ GRHydro_EigenproblemM.F90 \ GRHydro_FluxM.F90 \ -- cgit v1.2.3