[Commits] [svn:einsteintoolkit] GRHydro/trunk/ (Rev. 500)
rhaas at tapir.caltech.edu
rhaas at tapir.caltech.edu
Wed Mar 27 20:47:15 CDT 2013
User: rhaas
Date: 2013/03/27 08:47 PM
Added:
/trunk/src/
GRHydro_Con2PrimM_ptee.c
Modified:
/trunk/
interface.ccl
/trunk/src/
GRHydro_Con2PrimM.F90, make.code.defn
Log:
Finish the implementation of con2prim using the entropy equation.
* Tested successfully pointwisely. It still needs to be tested on evolution.
From: Bruno Coutinho Mundim <bcmsma at astro.rit.edu>
File Changes:
Directory: /trunk/src/
======================
File [modified]: GRHydro_Con2PrimM.F90
Delta lines: +34 -10
===================================================================
--- trunk/src/GRHydro_Con2PrimM.F90 2013-03-28 01:47:09 UTC (rev 499)
+++ trunk/src/GRHydro_Con2PrimM.F90 2013-03-28 01:47:14 UTC (rev 500)
@@ -334,7 +334,11 @@
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 @@
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 @@
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
File [added]: GRHydro_Con2PrimM_ptee.c
Delta lines: +592 -0
===================================================================
--- trunk/src/GRHydro_Con2PrimM_ptee.c (rev 0)
+++ trunk/src/GRHydro_Con2PrimM_ptee.c 2013-03-28 01:47:14 UTC (rev 500)
@@ -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 at astro.rit.edu .
+
+***********************************************************************************/
+
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdarg.h>
+#include <string.h>
+#include <math.h>
+#include <float.h>
+#include <complex.h>
+
+#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
File [modified]: make.code.defn
Delta lines: +1 -0
===================================================================
--- trunk/src/make.code.defn 2013-03-28 01:47:09 UTC (rev 499)
+++ trunk/src/make.code.defn 2013-03-28 01:47:14 UTC (rev 500)
@@ -46,6 +46,7 @@
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 \
Directory: /trunk/
==================
File [modified]: interface.ccl
Delta lines: +22 -0
===================================================================
--- trunk/interface.ccl 2013-03-28 01:47:09 UTC (rev 499)
+++ trunk/interface.ccl 2013-03-28 01:47:14 UTC (rev 500)
@@ -84,6 +84,27 @@
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 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
More information about the Commits
mailing list