|
HARM
harm and utilities
|
|
HARM
harm and utilities
|
Jon's U->P inversion. More...
Go to the source code of this file.
Macros | |
| #define | DEBUGINDEX 0 |
| 0 = no debug indexes set 1 = set global (file scope) index for debugging OPENMPMARK: This is not thread safe!, so only set to 1 if not using more than one OpenMP thread. More... | |
| #define | MAXNUMGUESSCHANGES (1000) |
| #define | FRACSs0 (0.5) |
| #define | FACTORSHIFT (2.0) |
| #define | MINRESIDEPRIME(Wp) ( Bsq*0.5+Qdotnp+3.0/5.0*Wp+(Bsq*Qtsq-QdotBsq)/(2.0*(Bsq+D+Wp)*(Bsq+D+Wp)) ) |
| residual for Eprime (old normal residual) when {u}^2=0, which is minimum value before problems More... | |
| #define | ALF 1.0e-4 |
| #define | N_CONV_TYPES 3 |
| #define | N_NITER_BINS (MAX_NEWT_ITER + 3) |
| #define | NBINS 200 |
Functions | |
| static int | coldgrhd (PFTYPE *lpflag, FTYPE Qtsq, FTYPE D, FTYPE *Wp) |
| re-look at this GODMARK Solution to the cold GRHD equation P^2 = ...[Wp] More... | |
| static int | vsqgtr1 (FTYPE W, FTYPE Bsq, FTYPE QdotBsq, FTYPE Qtsq) |
| returns v^2>=1 does NOT depend on EOS W^3 (W+2B^2) - S^2 (2W+B^2) W^2 (P^2-B^4) GODMARK: vsqgtr1 based upon energy equation information? What about P^2 equation inversion? More... | |
| static FTYPE | Eprime_Wp (FTYPE Wp, FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra) |
| static int | Psq_Wp (FTYPE Wp, FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra, FTYPE *resid, FTYPE *norm) |
| new version without poles returns actually W^2 P^2 Note that this does NOT use Qdotn or Qdotnp (from energy equation) so works for entropy evolution too More... | |
| static FTYPE | dPsqdWp_Wp (FTYPE Wp, FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra) |
| new version without poles result = d(W^2 P^2)/dWp = Note that this does NOT use Qdotn or Qdotnp (from energy equation) so works for entropy evolution too More... | |
| static FTYPE | gammasq_calc_W (FTYPE W, FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra) |
| returns ^2 Note that this does NOT use Qdotn or Qdotnp (from energy equation) so works for entropy evolution too More... | |
| static void | validate_Wp (FTYPE Wpold, FTYPE *Wpnew, FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra) |
| validate solution: More... | |
| static int | validate_vsq (FTYPE vsqold, FTYPE *vsqnew, FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra) |
| retvsq=0 : no problem retvsq=1 : >=1.0 retvsq=2 : < too negative retvsq=3 : >VSQ_TOO_BIG (forced failure even if converged) More... | |
| static void | check_utsq_fail (FTYPE Wp, FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra) |
| check on v^2>1 failure More... | |
| static FTYPE | Ss_W_vsq (int whicheos, FTYPE *EOSextra, FTYPE W, FTYPE D, FTYPE vsq) |
| EOS-related things for entropy evolution. More... | |
| static FTYPE | dSsdW_calc_vsq (int whicheos, FTYPE *EOSextra, FTYPE W, FTYPE D, FTYPE vsq) |
| static FTYPE | dSsdvsq_calc (int whicheos, FTYPE *EOSextra, FTYPE W, FTYPE D, FTYPE vsq) |
| static int | find_root_1D_gen (int showmessages, PFTYPE *lpflag, int eomtype, FTYPE x0, FTYPE *xnew, FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra, struct of_newtonstats *newtonstats) |
| Find Root wrappers: More... | |
| static int | find_root_1D_gen_scn (int showmessages, PFTYPE *lpflag, int eomtype, FTYPE x0, FTYPE *xnew, FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra, struct of_newtonstats *newtonstats) |
| static int | find_root_1D_gen_Eprime (int showmessages, PFTYPE *lpflag, int eomtype, FTYPE x0, FTYPE *xnew, FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra, struct of_newtonstats *newtonstats) |
| static int | find_root_1D_gen_Psq (int showmessages, PFTYPE *lpflag, int eomtype, FTYPE x0, FTYPE *xnew, FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra, struct of_newtonstats *newtonstats) |
| static int | find_root_3D_gen_Palpha (int showmessages, PFTYPE *lpflag, int eomtype, FTYPE x0, FTYPE *xnew, FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra, struct of_newtonstats *newtonstats) |
| below not done yet! More... | |
| static int | find_root_1D_gen_Sc (int showmessages, PFTYPE *lpflag, int eomtype, FTYPE x0, FTYPE *xnew, FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra, struct of_newtonstats *newtonstats) |
| static void | func_1d_orig_scn (FTYPE x[], FTYPE dx[], FTYPE resid[], FTYPE norm[], FTYPE(*jac)[NEWT_DIM], FTYPE *f, FTYPE *df, int n, FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra) |
| static FTYPE | res_sq_1d_orig_scn (FTYPE[], FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra) |
| static int | newt_repeatcheck (int n, FTYPE errx, FTYPE errx_old, FTYPE errx_oldest, FTYPE *dx, FTYPE *dx_old, FTYPE *x, FTYPE *x_old, FTYPE *x_older, FTYPE *x_olderer) |
| Newton checks: More... | |
| static int | newt_cyclecheck (int n, FTYPE trueerror, FTYPE errx, FTYPE errx_old, FTYPE errx_oldest, FTYPE *dx, FTYPE *dx_old, FTYPE *x, FTYPE *x_old, FTYPE *x_older, FTYPE *x_olderer, int n_iter, FTYPE(*xlist)[MAX_NEWT_ITER+EXTRA_NEWT_ITER_ULTRAREL+EXTRA_NEWT_ITER], FTYPE(*dxlist)[MAX_NEWT_ITER+EXTRA_NEWT_ITER_ULTRAREL+EXTRA_NEWT_ITER], FTYPE *errxlist, FTYPE *trueerrorlist) |
| see if repeating Newton cycle with no change in errx, so indicates can't do anything (even damping) to reduce error More... | |
| static int | newt_errorcheck (int n, FTYPE NEWT_TOL_VAR, FTYPE NEWT_TOL_ULTRAREL_VAR, FTYPE errx, FTYPE trueerrortotal, FTYPE x0, FTYPE dx0, FTYPE *x, FTYPE *dx, FTYPE *wglobal) |
| see if meet tolerance More... | |
| static int | newt_extracheck (int n, int EXTRA_NEWT_ITER_VAR, int EXTRA_NEWT_ITER_ULTRAREL_VAR, FTYPE errx, FTYPE x0, FTYPE dx0, FTYPE *x, FTYPE *dx, FTYPE *wglobal) |
| static int | compute_setup_quantities (FTYPE *prim, FTYPE *U, struct of_geom *ptrgeom, FTYPE *Qtcon, FTYPE *Bcon, FTYPE *Bcov, FTYPE *Bsq, FTYPE *QdotB, FTYPE *QdotBsq, FTYPE *Qtsq, FTYPE *Qdotn, FTYPE *Qdotnp, FTYPE *D, FTYPE *Sc, int whicheos, FTYPE *EOSextra) |
| Setup inversion: More... | |
| static int | compute_Wp (int showmessages, PFTYPE *lpflag, int eomtype, FTYPE *prim, struct of_geom *ptrgeom, FTYPE *W_last, FTYPE *Wp_last, FTYPE *Wpabs, FTYPE *etaabs, FTYPE *gamma, FTYPE *gammasq, FTYPE *rho0, FTYPE *u, FTYPE *p, FTYPE *wmrho0, FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra) |
| compute Wp from primitives also get normalizations wglobal More... | |
| static int | set_guess_Wp (int showmessages, PFTYPE *lpflag, int eomtype, FTYPE *prim, struct of_geom *ptrgeom, FTYPE *W_last, FTYPE *Wp_last, FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra) |
| SETUP ITERATIVE METHODS (good for GRMHD or entropy evolution or cold GRMHD) Doesn't use Sc since initial guess for W can be found from last primitives alone. More... | |
| static int | check_Wp (PFTYPE *lpflag, int eomtype, FTYPE *prim, FTYPE *U, struct of_geom *ptrgeom, FTYPE Wp_last, FTYPE Wp, int retval, FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra) |
| Check if solution was found. More... | |
| static int | Wp2prim (int showmessages, PFTYPE *lpflag, int eomtype, FTYPE *prim, FTYPE *pressure, FTYPE *U, struct of_geom *ptrgeom, FTYPE Wp, FTYPE *Qtcon, FTYPE *Bcon, FTYPE *Bcov, int retval, FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra) |
| static int | verify_Wlast (FTYPE u, FTYPE p, struct of_geom *ptrgeom, FTYPE *W_last, FTYPE *Wp_last, FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra) |
| Good for GRMHD or entropy GRMHD or cold GRMHD. More... | |
| static int | forcefree_inversion (int cleaned, struct of_geom *ptrgeom, FTYPE W, FTYPE *Qtcon, FTYPE Bsq, FTYPE *Bcon, FTYPE *Bcov, FTYPE Qtsq, FTYPE *U, FTYPE *prim, FTYPE *gammaret) |
| other inversions: More... | |
| int | Utoprim_jon_nonrelcompat_inputnorestmass (int showmessages, int eomtype, FTYPE *EOSextra, FTYPE *U, struct of_geom *ptrgeom, PFTYPE *lpflag, FTYPE *prim, FTYPE *pressure, struct of_newtonstats *newtonstats) |
| The ONLY global function: More... | |
| static FTYPE | gamma_calc_W (FTYPE W, FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra) |
| returns Note that this does NOT use Qdotn or Qdotnp (from energy equation) so works for entropy evolution too More... | |
| static int | Eprime_Wp_unopt (FTYPE Wp, FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra, FTYPE *resid, FTYPE *norm) |
| 0 = -E' + E'[Wp] for hot GRMHD (can be used for cold GRMHD too) while this has poles near bad roots, the E Wp for ultra relativistic case, so easy to find root More... | |
| static FTYPE | dEprimedWp_unopt (FTYPE Wp, FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra) |
| dE'/dWp [ Wp ] It is ok that there exists a term (1-dpdW) since even with ideal gas EOS in non-rel limit (1-dpdWp) still order unity as above, while poles near bad roots, linear dependence makes easier to find roots. More... | |
| static FTYPE | Eprime_Wp_new1 (FTYPE Wp, FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra) |
| this version has no poles assumes X2!=0, in which case previous old version gives inf anyways 0 = (B^2+W)^2 (-E' + E'[Wp]) for hot GRMHD (can be used for cold GRMHD too) More... | |
| static FTYPE | dEprimedWp_new1 (FTYPE Wp, FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra) |
| this version has no poles d((B^2+W)^2 E')/dWp [ Wp ] It is ok that there exists a term (1-dpdW) since even with ideal gas EOS in non-rel limit (1-dpdWp) still order unity More... | |
| static int | Sc_Wp_unopt (FTYPE Wp, FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra, FTYPE *resid, FTYPE *norm) |
| 0 = Wp(-Sc + Sc[Wp]) for entropy GRMHD Noticed that original residual behaves like log(Wp) near the root, so that the derivative shoots too far. More... | |
| static FTYPE | dScdWp_unopt (FTYPE Wp, FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra) |
| dSc/dWp [ Wp ] for entropy GRMHD More... | |
| static void | validate_x_2d (FTYPE x[NEWT_DIM], FTYPE x0[NEWT_DIM], FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra) |
| validate functions here change where Newton is, so ok to change arbitrarily since next Newton step will be in well-defined place (although may get stuck) More... | |
| static void | validate_x_1d (FTYPE x[1], FTYPE x0[1], FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra) |
| 1D methods using W' More... | |
| static void | pick_validate_x (int eomtype, void(**ptr_validate_x)(FTYPE x[NEWT_DIM], FTYPE x0[NEWT_DIM], FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra)) |
| set pointer to the correct validate_x function More... | |
| static void | func_Eprime_unopt (FTYPE x[], FTYPE dx[], FTYPE resid[], FTYPE norm[], FTYPE(*jac)[NEWT_DIM], FTYPE *f, FTYPE *df, int n, FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra) |
| unoptimized vresion of func_Eprime More... | |
| static FTYPE | res_sq_Eprime_unopt (FTYPE x[], FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra) |
| res_sq_Eprime() (or any pure residual function) is only called by line searching method More... | |
| static void | get_kinetics_part1 (FTYPE Wp, FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra, FTYPE *W, FTYPE *X, FTYPE *X2, FTYPE *X3, FTYPE *utsq, FTYPE *gamma, FTYPE *gammasq, FTYPE *rho0, FTYPE *wmrho0) |
| get first part of kinetic stuff for func() or resid_sq() for optimized version of func() More... | |
| static void | get_kinetics_part2 (FTYPE Wp, FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra, FTYPE W, FTYPE X, FTYPE X2, FTYPE X3, FTYPE utsq, FTYPE gamma, FTYPE gammasq, FTYPE rho0, FTYPE wmrho0, FTYPE *dvsq, FTYPE *dwmrho0dW, FTYPE *drho0dW, FTYPE *dwmrho0dvsq, FTYPE *drho0dvsq) |
| get kinetic stuff for func() or resid_sq() for optimized version of func() More... | |
| static void | func_Eprime_opt (FTYPE x[], FTYPE dx[], FTYPE resid[], FTYPE norm[], FTYPE(*jac)[NEWT_DIM], FTYPE *f, FTYPE *df, int n, FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra) |
| optimized version of func_Eprime() optimized in sense that minimizes extra kinetic calculations AND EOS lookups More... | |
| static FTYPE | res_sq_Eprime_opt (FTYPE x[], FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra) |
| res_sq_Eprime() (or any pure residual function) is only called by line searching method More... | |
| static void | func_Psq (FTYPE x[], FTYPE dx[], FTYPE resid[], FTYPE norm[], FTYPE(*jac)[NEWT_DIM], FTYPE *f, FTYPE *df, int n, FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra) |
| not that this is well-optimized. Psq_Wp() and dPsqdWp_Wp() don't use subfunctions like "_unopt" versions for hot and entropy inversions. Also no EOS functions to call. More... | |
| static FTYPE | res_sq_Psq (FTYPE x[], FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra) |
| static void | func_Sc_unopt (FTYPE x[], FTYPE dx[], FTYPE resid[], FTYPE norm[], FTYPE(*jac)[NEWT_DIM], FTYPE *f, FTYPE *df, int n, FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra) |
| static FTYPE | res_sq_Sc_unopt (FTYPE x[], FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra) |
| static void | func_Sc_opt (FTYPE x[], FTYPE dx[], FTYPE resid[], FTYPE norm[], FTYPE(*jac)[NEWT_DIM], FTYPE *f, FTYPE *df, int n, FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra) |
| optimized version of func_Sc() very similar to optimized version of func_Eprime() but different "getall_forinversion()" call and different residual and derivative More... | |
| static FTYPE | res_sq_Sc_opt (FTYPE x[], FTYPE *wglobal, FTYPE Bsq, FTYPE QdotB, FTYPE QdotBsq, FTYPE Qtsq, FTYPE Qdotn, FTYPE Qdotnp, FTYPE D, FTYPE Sc, int whicheos, FTYPE *EOSextra) |
| optimized version of res_sq_Sc() just calls func_Sc_opt() for now since more optimized than original res_sq_Sc() still. More... | |
Jon's U->P inversion.
OPENMPNOTE: This function should only have local variables in order to be parallelized, otherwise global shared variables can change while computing. OPENMPNOTE: Should avoid use of static variables inside functions, for example when first time called to initialize things. Only done right now in bin_newt_data() that is not used.
Definition in file utoprim_jon.c.
| #define ALF 1.0e-4 |
Definition at line 4968 of file utoprim_jon.c.
| #define DEBUGINDEX 0 |
0 = no debug indexes set 1 = set global (file scope) index for debugging OPENMPMARK: This is not thread safe!, so only set to 1 if not using more than one OpenMP thread.
Definition at line 17 of file utoprim_jon.c.
| #define FACTORSHIFT (2.0) |
| #define FRACSs0 (0.5) |
| #define MAXNUMGUESSCHANGES (1000) |
| #define MINRESIDEPRIME | ( | Wp | ) | ( Bsq*0.5+Qdotnp+3.0/5.0*Wp+(Bsq*Qtsq-QdotBsq)/(2.0*(Bsq+D+Wp)*(Bsq+D+Wp)) ) |
residual for Eprime (old normal residual) when {u}^2=0, which is minimum value before problems
Definition at line 2099 of file utoprim_jon.c.
| #define N_CONV_TYPES 3 |
Definition at line 5110 of file utoprim_jon.c.
| #define N_NITER_BINS (MAX_NEWT_ITER + 3) |
Definition at line 5111 of file utoprim_jon.c.
| #define NBINS 200 |
Definition at line 5112 of file utoprim_jon.c.
|
static |
check on v^2>1 failure
Definition at line 1404 of file utoprim_jon.c.
|
static |
Check if solution was found.
wglobal[1]; // normalize to old densities
Definition at line 1299 of file utoprim_jon.c.
re-look at this GODMARK Solution to the cold GRHD equation P^2 = ...[Wp]
Definition at line 1977 of file utoprim_jon.c.
|
static |
Setup inversion:
returns global quanities Bsq,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc returns other local quantities
Definition at line 789 of file utoprim_jon.c.
|
static |
compute Wp from primitives also get normalizations wglobal
Definition at line 958 of file utoprim_jon.c.
|
static |
this version has no poles d((B^2+W)^2 E')/dWp [ Wp ] It is ok that there exists a term (1-dpdW) since even with ideal gas EOS in non-rel limit (1-dpdWp) still order unity
Definition at line 2221 of file utoprim_jon.c.
|
static |
dE'/dWp [ Wp ] It is ok that there exists a term (1-dpdW) since even with ideal gas EOS in non-rel limit (1-dpdWp) still order unity as above, while poles near bad roots, linear dependence makes easier to find roots.
OPTMARK: Note that I already tried removing if-else and has no effect (05/20/09)
Definition at line 2145 of file utoprim_jon.c.
|
static |
new version without poles result = d(W^2 P^2)/dWp = Note that this does NOT use Qdotn or Qdotnp (from energy equation) so works for entropy evolution too
Definition at line 2076 of file utoprim_jon.c.
|
static |
dSc/dWp [ Wp ] for entropy GRMHD
Definition at line 2342 of file utoprim_jon.c.
|
static |
|
static |
this version has no poles assumes X2!=0, in which case previous old version gives inf anyways 0 = (B^2+W)^2 (-E' + E'[Wp]) for hot GRMHD (can be used for cold GRMHD too)
Definition at line 2182 of file utoprim_jon.c.
|
static |
0 = -E' + E'[Wp] for hot GRMHD (can be used for cold GRMHD too) while this has poles near bad roots, the E Wp for ultra relativistic case, so easy to find root
Definition at line 2103 of file utoprim_jon.c.
|
static |
Find Root wrappers:
Definition at line 3817 of file utoprim_jon.c.
|
static |
below not done yet!
Definition at line 3903 of file utoprim_jon.c.
|
static |
other inversions:
appears to give qualitatively similar, but still qualitatively different results than original phys.ffde.c code. Gives more oscillatory results for torus DISKFIELD force-free field (i.e. set rho=u=0 but otherwise start with init.fishmon.c type setup)
Definition at line 1721 of file utoprim_jon.c.
|
static |
optimized version of func_Eprime() optimized in sense that minimizes extra kinetic calculations AND EOS lookups
Definition at line 3125 of file utoprim_jon.c.
|
static |
unoptimized vresion of func_Eprime
Definition at line 3031 of file utoprim_jon.c.
|
static |
not that this is well-optimized. Psq_Wp() and dPsqdWp_Wp() don't use subfunctions like "_unopt" versions for hot and entropy inversions. Also no EOS functions to call.
Definition at line 3279 of file utoprim_jon.c.
|
static |
optimized version of func_Sc() very similar to optimized version of func_Eprime() but different "getall_forinversion()" call and different residual and derivative
Definition at line 3357 of file utoprim_jon.c.
|
static |
returns Note that this does NOT use Qdotn or Qdotnp (from energy equation) so works for entropy evolution too
Definition at line 1911 of file utoprim_jon.c.
|
static |
returns ^2 Note that this does NOT use Qdotn or Qdotnp (from energy equation) so works for entropy evolution too
Definition at line 1890 of file utoprim_jon.c.
|
static |
get first part of kinetic stuff for func() or resid_sq() for optimized version of func()
Definition at line 3075 of file utoprim_jon.c.
|
static |
get kinetic stuff for func() or resid_sq() for optimized version of func()
Definition at line 3104 of file utoprim_jon.c.
|
static |
see if repeating Newton cycle with no change in errx, so indicates can't do anything (even damping) to reduce error
Definition at line 4910 of file utoprim_jon.c.
|
static |
see if meet tolerance
Definition at line 4863 of file utoprim_jon.c.
|
static |
Definition at line 4944 of file utoprim_jon.c.
|
static |
Newton checks:
see if repeating Newton cycle with no change in errx, so indicates can't do anything (even damping) to reduce error
Definition at line 4886 of file utoprim_jon.c.
|
static |
set pointer to the correct validate_x function
Definition at line 2578 of file utoprim_jon.c.
|
static |
new version without poles returns actually W^2 P^2 Note that this does NOT use Qdotn or Qdotnp (from energy equation) so works for entropy evolution too
Definition at line 2055 of file utoprim_jon.c.
|
static |
res_sq_Eprime() (or any pure residual function) is only called by line searching method
Definition at line 3259 of file utoprim_jon.c.
|
static |
res_sq_Eprime() (or any pure residual function) is only called by line searching method
Definition at line 3057 of file utoprim_jon.c.
|
static |
Definition at line 3301 of file utoprim_jon.c.
|
static |
optimized version of res_sq_Sc() just calls func_Sc_opt() for now since more optimized than original res_sq_Sc() still.
Could optimize further, but don't require only residual unless doing line searching that don't currently do.
Definition at line 3476 of file utoprim_jon.c.
|
static |
0 = Wp(-Sc + Sc[Wp]) for entropy GRMHD Noticed that original residual behaves like log(Wp) near the root, so that the derivative shoots too far.
Multiplying by Wp solves that problem.
Definition at line 2307 of file utoprim_jon.c.
|
static |
SETUP ITERATIVE METHODS (good for GRMHD or entropy evolution or cold GRMHD) Doesn't use Sc since initial guess for W can be found from last primitives alone.
Definition at line 1046 of file utoprim_jon.c.
EOS-related things for entropy evolution.
| int Utoprim_jon_nonrelcompat_inputnorestmass | ( | int | showmessages, |
| int | eomtype, | ||
| FTYPE * | EOSextra, | ||
| FTYPE * | U, | ||
| struct of_geom * | ptrgeom, | ||
| PFTYPE * | lpflag, | ||
| FTYPE * | prim, | ||
| FTYPE * | pressure, | ||
| struct of_newtonstats * | newtonstats | ||
| ) |
The ONLY global function:
Definition at line 142 of file utoprim_jon.c.
|
static |
retvsq=0 : no problem retvsq=1 : >=1.0 retvsq=2 : < too negative retvsq=3 : >VSQ_TOO_BIG (forced failure even if converged)
Definition at line 2492 of file utoprim_jon.c.
|
static |
validate solution:
Wpold is previous Newton iteration's value *Wpnew is updated Newton iterations value.
Definition at line 2453 of file utoprim_jon.c.
|
static |
1D methods using W'
Definition at line 2569 of file utoprim_jon.c.
|
static |
validate functions here change where Newton is, so ok to change arbitrarily since next Newton step will be in well-defined place (although may get stuck)
Definition at line 2541 of file utoprim_jon.c.
|
static |
Good for GRMHD or entropy GRMHD or cold GRMHD.
Make sure that W is large enough so that v^2 < 1 :
GODMARK: apparently this is necessary for 1D method, otherwise it can't find solution (2D method seems to be ok)
Definition at line 1229 of file utoprim_jon.c.
returns v^2>=1 does NOT depend on EOS W^3 (W+2B^2) - S^2 (2W+B^2) W^2 (P^2-B^4) GODMARK: vsqgtr1 based upon energy equation information? What about P^2 equation inversion?
Definition at line 1842 of file utoprim_jon.c.
|
static |
1.8.3.1
