HARM
harm and utilities
 All Data Structures Files Functions Variables Typedefs Macros Pages
utoprim_jon.c
Go to the documentation of this file.
1 
10 #include "u2p_defs.h" // only includes #define's
11 #include "utoprim_jon.h" // only includes #define's
12 
13 
17 #define DEBUGINDEX 0
18 
19 #if(DEBUGINDEX)
20 //
21 static int ifileglobal,jfileglobal,kfileglobal,pfileglobal;
22 #endif
23 
24 
26 static int Utoprim_new_body(int showmessages, int eomtype, PFTYPE *lpflag, int whicheos, FTYPE *EOSextra, FTYPE U[NPR], struct of_geom *ptrgeom, FTYPE *prim, FTYPE *pressure, struct of_newtonstats *newtonstats);
27 
29 static void raise_g(FTYPE vcov[], FTYPE *gcon, FTYPE vcon[]);
30 static void lower_g(FTYPE vcon[], FTYPE *gcov, FTYPE vcov[]);
31 static void ncov_calc_fromlapse(FTYPE lapse,FTYPE ncov[]) ;
32 static void ncov_calc(FTYPE gcon[SYMMATRIXNDIM],FTYPE ncov[]) ;
33 static int coldgrhd(PFTYPE *lpflag, FTYPE Qtsq, FTYPE D, FTYPE *Wp);
34 static int vsqgtr1(FTYPE W,FTYPE Bsq,FTYPE QdotBsq, FTYPE Qtsq);
35 
37 static FTYPE vsq_calc(FTYPE W, FTYPE *wglobal,FTYPE Bsq,FTYPE QdotB,FTYPE QdotBsq,FTYPE Qtsq,FTYPE Qdotn,FTYPE Qdotnp,FTYPE D,FTYPE Sc, int whicheos, FTYPE *EOSextra);
38 static FTYPE dvsq_dW(FTYPE W, FTYPE *wglobal,FTYPE Bsq,FTYPE QdotB,FTYPE QdotBsq,FTYPE Qtsq,FTYPE Qdotn,FTYPE Qdotnp,FTYPE D,FTYPE Sc, int whicheos, FTYPE *EOSextra);
39 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);
40 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);
41 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);
42 static FTYPE utsq_calc(FTYPE W, FTYPE *wglobal,FTYPE Bsq,FTYPE QdotB,FTYPE QdotBsq,FTYPE Qtsq,FTYPE Qdotn,FTYPE Qdotnp,FTYPE D,FTYPE Sc, int whicheos, FTYPE *EOSextra);
43 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);
44 static FTYPE utsq_calc(FTYPE W, FTYPE *wglobal,FTYPE Bsq,FTYPE QdotB,FTYPE QdotBsq,FTYPE Qtsq,FTYPE Qdotn,FTYPE Qdotnp,FTYPE D,FTYPE Sc, int whicheos, FTYPE *EOSextra);
45 
47 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);
48 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);
49 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);
50 
52 static FTYPE pressure_W_vsq(int whicheos, FTYPE *EOSextra, FTYPE W, FTYPE D, FTYPE vsq) ;
53 static FTYPE pressure_Wp_utsq(int whicheos, FTYPE *EOSextra, FTYPE Wp, FTYPE D, FTYPE utsq);
54 static FTYPE dpdW_calc_vsq(int whicheos, FTYPE *EOSextra, FTYPE W, FTYPE D, FTYPE vsq);
55 static FTYPE dpdWp_calc_utsq(int whicheos, FTYPE *EOSextra, FTYPE Wp, FTYPE D, FTYPE utsq);
56 static FTYPE dpdvsq_calc2_Wp(int whicheos, FTYPE *EOSextra, FTYPE Wp, FTYPE D, FTYPE utsq);
57 static FTYPE dpdvsq_calc(int whicheos, FTYPE *EOSextra, FTYPE W, FTYPE D, FTYPE vsq);
58 static FTYPE dpdWp_calc_utsq(int whicheos, FTYPE *EOSextra, FTYPE Wp, FTYPE D, FTYPE utsq);
59 
61 static FTYPE Ss_W_vsq(int whicheos, FTYPE *EOSextra, FTYPE W, FTYPE D, FTYPE vsq) ;
62 static FTYPE Ss_Wp_utsq(int whicheos, FTYPE *EOSextra, FTYPE Wp, FTYPE D, FTYPE utsq);
63 static FTYPE dSsdW_calc_vsq(int whicheos, FTYPE *EOSextra, FTYPE W, FTYPE D, FTYPE vsq);
64 static FTYPE dSsdWp_calc_utsq(int whicheos, FTYPE *EOSextra, FTYPE Wp, FTYPE D, FTYPE utsq);
65 static FTYPE dSsdvsq_calc2_Wp(int whicheos, FTYPE *EOSextra, FTYPE Wp, FTYPE D, FTYPE utsq);
66 static FTYPE dSsdvsq_calc(int whicheos, FTYPE *EOSextra, FTYPE W, FTYPE D, FTYPE vsq);
67 static FTYPE dSsdWp_calc_utsq(int whicheos, FTYPE *EOSextra, FTYPE Wp, FTYPE D, FTYPE utsq);
68 
70 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);
71 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);
72 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);
73 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);
74 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);
75 static int find_root_2D_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);
76 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);
77 
79 static int general_newton_raphson(int showmessages, PFTYPE *lpflag, int eomtype, FTYPE x[], int n, int do_line_search,
80  void (*funcd) (FTYPE [], FTYPE [], FTYPE [], FTYPE [], FTYPE [][NEWT_DIM], FTYPE *, FTYPE *, int, FTYPE *wglobal,FTYPE Bsq,FTYPE QdotB,FTYPE QdotBsq,FTYPE Qtsq,FTYPE Qdotn,FTYPE Qdotnp,FTYPE D,FTYPE Sc, int whicheos, FTYPE *EOSextra),
81  FTYPE (*res_func) (FTYPE [], FTYPE *wglobal,FTYPE Bsq,FTYPE QdotB,FTYPE QdotBsq,FTYPE Qtsq,FTYPE Qdotn,FTYPE Qdotnp,FTYPE D,FTYPE Sc, int whicheos, FTYPE *EOSextra), 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);
82 
83 static void func_1d_orig(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);
84 static FTYPE res_sq_1d_orig( FTYPE [], FTYPE *wglobal,FTYPE Bsq,FTYPE QdotB,FTYPE QdotBsq,FTYPE Qtsq,FTYPE Qdotn,FTYPE Qdotnp,FTYPE D,FTYPE Sc, int whicheos, FTYPE *EOSextra);
85 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);
86 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);
87 static void func_1d_orig(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);
88 static FTYPE res_sq_1d_orig( FTYPE [], FTYPE *wglobal,FTYPE Bsq,FTYPE QdotB,FTYPE QdotBsq,FTYPE Qtsq,FTYPE Qdotn,FTYPE Qdotnp,FTYPE D,FTYPE Sc, int whicheos, FTYPE *EOSextra);
89 static void func_1d_orig(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);
90 static FTYPE res_sq_1d_orig( FTYPE [] , FTYPE *wglobal,FTYPE Bsq,FTYPE QdotB,FTYPE QdotBsq,FTYPE Qtsq,FTYPE Qdotn,FTYPE Qdotnp,FTYPE D,FTYPE Sc, int whicheos, FTYPE *EOSextra);
91 static void func_1d_orig(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);
92 static FTYPE res_sq_1d_orig( FTYPE [] , FTYPE *wglobal,FTYPE Bsq,FTYPE QdotB,FTYPE QdotBsq,FTYPE Qtsq,FTYPE Qdotn,FTYPE Qdotnp,FTYPE D,FTYPE Sc, int whicheos, FTYPE *EOSextra);
93 static void func_vsq( FTYPE [], FTYPE [], FTYPE [], FTYPE [], FTYPE [][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);
94 static void func_vsq2( FTYPE [], FTYPE [], FTYPE [], FTYPE [], FTYPE [][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);
95 static FTYPE res_sq_vsq( FTYPE [], FTYPE *wglobal,FTYPE Bsq,FTYPE QdotB,FTYPE QdotBsq,FTYPE Qtsq,FTYPE Qdotn,FTYPE Qdotnp,FTYPE D,FTYPE Sc, int whicheos, FTYPE *EOSextra);
96 static FTYPE res_sq_vsq2( FTYPE [] , FTYPE *wglobal,FTYPE Bsq,FTYPE QdotB,FTYPE QdotBsq,FTYPE Qtsq,FTYPE Qdotn,FTYPE Qdotnp,FTYPE D,FTYPE Sc, int whicheos, FTYPE *EOSextra);
97 static int x1_of_x0(FTYPE x0, FTYPE *x1, FTYPE *wglobal,FTYPE Bsq,FTYPE QdotB,FTYPE QdotBsq,FTYPE Qtsq,FTYPE Qdotn,FTYPE Qdotnp,FTYPE D,FTYPE Sc, int whicheos, FTYPE *EOSextra);
98 static void my_lnsrch(int eomtype, int n, FTYPE xold[], FTYPE fold, FTYPE g[], FTYPE p[], FTYPE x[],
99  FTYPE *f, FTYPE TOLX, FTYPE stpmax, int *check, FTYPE (*func) (FTYPE [], FTYPE *wglobal,FTYPE Bsq,FTYPE QdotB,FTYPE QdotBsq,FTYPE Qtsq,FTYPE Qdotn,FTYPE Qdotnp,FTYPE D,FTYPE Sc, int whicheos, FTYPE *EOSextra), FTYPE *wglobal,FTYPE Bsq,FTYPE QdotB,FTYPE QdotBsq,FTYPE Qtsq,FTYPE Qdotn,FTYPE Qdotnp,FTYPE D,FTYPE Sc, int whicheos, FTYPE *EOSextra);
100 
101 
103 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);
104 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);
105 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);
106 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);
107 static void bin_newt_data(FTYPE MAX_NEWT_ITER_VAR, FTYPE errx, int niters, int conv_type, int print_now ) ;
108 
109 static void validate_x(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),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);
110 static void check_on_inversion(int showmessages, FTYPE *prim, 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, struct of_newtonstats *newtonstats);
111 
112 
113 
115 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);
116 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);
117 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);
118 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);
119 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);
120 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);
121 
122 
124 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);
125 
126 
127 
128 
129 
130 #include "u2p_util.c" // only includes static functions
131 #include "utoprim_jon_eos.c" // only includes static functions
132 
133 
134 /******************************************************************
135 
136  Driver for new prim. var. solver. The driver just translates
137  between the two sets of definitions for U and P.
138 
139 ******************************************************************/
140 
142 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)
143 {
144 
145  FTYPE U_tmp[NPR], U_tmp2[NPR], prim_tmp[NPR];
146  int i, j, k,ret;
147  int pl,pl2;
148  FTYPE gcov[SYMMATRIXNDIM], gcon[SYMMATRIXNDIM], alpha;
149  const int ltrace = 0;
150  /* these variables need to be shared between the functions
151  Utoprim_1D, residual, and utsq */
152  int real_dimen_newton;
153  int whicheos;
154  int pliter;
155 
156  // showmessages=1; // FORCE
157 
158  // DEBUG:
159  // PLOOP(pliter,pl) dualfprintf(fail_file,"%d %d %d : pl=%d U=%21.15g\n",ptrgeom->i,ptrgeom->j,ptrgeom->k,pl,U[pl]);
160 
161 
162 #if(0)
163  U[0]= 100.03640347219;
164  U[1]=-2.21974020354085e-25;
165  U[2]=-2.50158167405965e-12;
166  U[3]= 0;
167  U[4]= 0;
168  U[5]= 0;
169  U[6]= 0;
170  U[7]= 0;
171 
172  prim[0]= 99.8798655540244;
173  prim[1]= 9.6414760309742e-26;
174  prim[2]=-9.98706197921322e-14;
175  prim[3]= 0;
176  prim[4]= 0;
177  prim[5]= 0;
178  prim[6]= 0;
179  prim[7]= 0;
180 #endif
181 
182 
183 
185  //
186  // EOS functions used during inversion
187  // GODMARK: WHICHEOS could be input instead of using global definition
188  //
190  pickeos_eomtype(WHICHEOS,eomtype,&whicheos);
191 
192 
193 
194 
195 #if(!OPTIMIZED)
196  // make sure nan and inf's are not in initial guess or U
197  PLOOPALLINVERT(pl){
198  if(!isfinite(prim[pl])){
199  dualfprintf(fail_file,"Guess for p is NAN\n");
200  PLOOPALLINVERT(pl2) dualfprintf(fail_file,"prim[%d]=%21.15g\n",pl2,prim[pl2]);
201  *lpflag= UTOPRIMFAILNANGUESS;
202  myexit(876);
203  }
204  if(!isfinite(U[pl])){
205  dualfprintf(fail_file,"Value of U is NAN\n");
206  PLOOPALLINVERT(pl2) dualfprintf(fail_file,"prim[%d]=%21.15g\n",pl2,U[pl2]);
207  *lpflag= UTOPRIMFAILNANGUESS;
208  myexit(877);
209  }
210  }
211 #endif
212 
213 
214 
216  //
217  // initialize iteration counter
218  //
220  newtonstats->lntries=0;
221 
222 
223 
225  //
226  // Setup dimensionality of method for a given eomtype and choice if inverter
227  //
229  if(eomtype==EOMGRMHD){
230  if(WHICHHOTINVERTER==1) real_dimen_newton=1;
231  else if(WHICHHOTINVERTER==2) real_dimen_newton=2;
232  else if(WHICHHOTINVERTER==3) real_dimen_newton=1;
233  else{
234  dualfprintf(fail_file,"No such hot inverter %d\n",WHICHHOTINVERTER);
235  myexit(89285323);
236  }
237  }
238  else if(eomtype==EOMCOLDGRMHD){
239  if(WHICHCOLDINVERTER==0) real_dimen_newton=1;
240  else if(WHICHCOLDINVERTER==1) real_dimen_newton=1;
241  else if(WHICHCOLDINVERTER==2) real_dimen_newton=1;
242  else if(WHICHCOLDINVERTER==3) real_dimen_newton=1;
243  else{
244  dualfprintf(fail_file,"No such cold inverter %d\n",WHICHCOLDINVERTER);
245  myexit(89285324);
246  }
247  }
248  else if(eomtype==EOMENTROPYGRMHD){
249  if(WHICHENTROPYINVERTER==0) real_dimen_newton=1;
250  else{
251  dualfprintf(fail_file,"No such entropy inverter %d\n",WHICHENTROPYINVERTER);
252  myexit(89285325);
253  }
254  }
255  else real_dimen_newton=0; // not using Newton's method
256 
257 
258 
259 #if( WHICHVEL != VELREL4 )
260  stderrfprintf("Utoprim_2d() Not implemented for WHICHVEL = %d \n", WHICHVEL );
261  return(1);
262 #endif
263 
264 
265 
267  //
268  // First update the primitive B-fields
269  //
271  for(i = BCON1; i <= BCON3; i++) prim[i] = U[i] ;
272 
273 
274 
275 
277  //
278  // Calculate the transform the CONSERVATIVE variables into the new system
279  //
281 
282  // Set the geometry variables
283  // alpha = 1.0/sqrt(-ptrgeom->gcon[GIND(0,0)]);
284  alpha = ptrgeom->alphalapse;
285 
286  for( i = RHO; i <= BCON3; i++ ) {
287  U_tmp[i] = alpha * U[i] ;
288  }
289 
290 #if(DOENTROPY!=DONOENTROPY)
291  i=ENTROPY;
292  U_tmp[i] = alpha * U[i] ;
293 #endif
294 
295 
296 #if(CRAZYDEBUG)
297  if(ptrgeom->i==0 && ptrgeom->j==31 && nstep==9 && steppart==2){
298  PLOOPALLINVERT(k) dualfprintf(fail_file,"U_tmp[%d]=%21.15g\n",k,U_tmp[k]);
299  }
300 #endif
301 
302 
304  //
305  // Calculate the transform the PRIMITIVE variables into the new system
306  //
308 
309  // default is to treat all quantities as scalars
310  // Note that need PALLLOOP here since while not setting all quantities, may use all quantities
311  PALLLOOP(i) prim_tmp[i] = prim[i];
312  // field
313  for( i = BCON1; i <= BCON3; i++ ) prim_tmp[i] = alpha*prim[i];
314 
315 
316 
318  //
319  // Perform the inversion from U -> P
320  //
322  ret = Utoprim_new_body(showmessages, eomtype, lpflag, whicheos,EOSextra, U_tmp, ptrgeom, prim_tmp,pressure,newtonstats);
323 
324 
325 
327  //
328  // Check conservative variable transformation
329  //
331 #if(!OPTIMIZED)
332  if( ltrace ) {
333  for(i = 0; i < NDIM; i++ ) {
334  for(j = 0; j < NDIM; j++ ) {
335  gcov[GIND(i,j)] = ptrgeom->gcov[GIND(i,j)];
336  gcon[GIND(i,j)] = ptrgeom->gcon[GIND(i,j)];
337  dualfprintf(fail_file,"gcov,gcon %d %d = %21.15g %21.15g \n", i, j, gcov[GIND(i,j)], gcon[GIND(i,j)]);
338  }
339  }
340  dualfprintf(fail_file,"gdet = %21.15g \n", ptrgeom->g);
341  // primtoU_g(ptrgeom, prim_tmp, gcov, gcon, U_tmp2 );
342  // PALLLOOP(i){
343  // dualfprintf(fail_file, "Utoprim_1d(): Utmp1[%d] = %21.15g , Utmp2[%d] = %21.15g , dUtmp[%d] = %21.15g \n",
344  // i, U_tmp[i], i, U_tmp2[i], i, fabs( (U_tmp[i]-U_tmp2[i]) / ( (U_tmp2[i]!=0.) ? U_tmp2[i] : 1. ) ) );
345  // }
346  }
347 #endif
348 
349 
350 
352  //
353  // Transform new primitive variables back
354  //
356  if( ret == 0 ) {
357  // only changed RHO through U3 (otherwise have to remove alpha from B^i)
358  for( i = 0; i <=U3; i++ ) prim[i] = prim_tmp[i];
359 #if(DOENTROPY!=DONOENTROPY)
360  // if entropy is also copy entropy primitive that is internal energy density
361  i=ENTROPY;
362  prim[i] = prim_tmp[i];
363 #endif
364 
365  }
366  else{
367  // ensure failure reported if ret!=0
368  if(IFUTOPRIMNOFAILORFIXED(*lpflag)){
369  // dualfprintf(fail_file,"Internal ret=%d failure\n",ret);
370  *lpflag= ret+UTOPRIMFAILCONVRET+1;// related to UTOPRIMFAILCONVRET
371  }
372  }
373 
374 
375  // PLOOP(pliter,pl) dualfprintf(fail_file,"internal pl=%d prim=%g\n",pl,prim[pl]);
376 
377 
379  //
380  // revert EOS for WHICHEOS,EOMTYPE
381  // not really necessary unless whicheos used after the below call
382  //
384  pickeos_eomtype(WHICHEOS,EOMTYPE,&whicheos); // see comments in eos.c for pickeos_eomtype_old()
385 
386  return( 0 ) ;
387 
388 }
389 
390 
391 /**********************************************************************************
392 
393 attempt an inversion from U to prim using the initial guess
394 prim.
395 
396  -- assumes that
397  / rho gamma \
398  U = | alpha T^t_\mu |
399  \ alpha B^i /
400  \ ... /
401  \ specific entropy /
402 
403 
404  / rho \
405  P = | uu |
406  | \tilde{u}^i |
407  \ alpha B^i /
408  \ ... /
409  \ uu from entropy (copied to uu) /
410 
411 
412  *lpflag: (i*100 + j) where
413  i = 0 -> Newton-Raphson solver either was not called (yet or not used) or returned successfully;
414  1 -> Newton-Raphson solver did not converge to a solution with the given tolerances;
415  2 -> Newton-Raphson procedure encountered a numerical divergence (occurrence of "nan" or "+/-inf" ;
416 
417  j = 0 -> success
418  1 -> failure: some sort of failure in Newton-Raphson;
419  2 -> failure: utsq<0 w/ initial p[] guess;
420  3 -> failure: W<0 or W>W_TOO_BIG
421  4 -> failure: utsq<0 or utsq > UTSQ_TOO_BIG with new W;
422  5 -> failure: rho,uu <= 0 ;
423 
424 **********************************************************************************/
425 
426 static int Utoprim_new_body(int showmessages, int eomtype, PFTYPE *lpflag, int whicheos, FTYPE *EOSextra, FTYPE U[NPR], struct of_geom *ptrgeom, FTYPE *prim, FTYPE *pressure, struct of_newtonstats *newtonstats)
427 {
428  FTYPE Qconp[NDIM],Qcovp[NDIM];
429  FTYPE Qtcon[4],Bcon[4],Bcov[4];
430  FTYPE primother[NPR];
431  FTYPE W_last,W,Wp_last,Wp0,Wp,Wp1,Wp2,Wp3,Wpother;
432  int i,j, retval0, retval, retval1,retval2, retval3 ;
433  int retvalother,retvalother2;
434  // below used to be static global in this file
435  FTYPE Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc ;
436  FTYPE wglobal[3]; // [0] is normalization for Wp and [1] is normalization for Wtest, and [2] is for rho0gammasqresidabs term
437  int Wdependentsolution;
438 
439 
440 #if(0) // old globals used for debugging -- can't be used during OpenMP computations
441  static FTYPE Qcovorig[NDIM],Qconorig[NDIM],Qdotnorig,Qsqorig,Qtsqorig;
442  static FTYPE Qtsqnew;
443 #endif
444 
445 
446 
447 #if(DEBUGINDEX)
448  // initialize debug indices
449  ifileglobal=ptrgeom->i;
450  jfileglobal=ptrgeom->j;
451  kfileglobal=ptrgeom->k;
452  pfileglobal=ptrgeom->p;
453 #endif
454 
455  const int ltrace = 0;
456  const int ltrace2 = 0;
457 
458 
459 
460 
461 
462 #if(!OPTIMIZED)
463  if( ltrace ) {
464  for(i=0;i<NPR;i++) dualfprintf(fail_file,"%d %21.15g %21.15g\n",i,prim[i],U[i]) ;
465 
466  }
467 #endif
468 
469 
470 
471 
473  //
474  // Assume no inversion failure initially
475  //
477  *lpflag= UTOPRIMNOFAIL;
478  retval=0;
479 
480 
482  //
483  // compute quantities necessary for inversion in a certain frame
484  //
486  compute_setup_quantities(prim, U, ptrgeom, Qtcon, Bcon, Bcov, &Bsq,&QdotB,&QdotBsq,&Qtsq,&Qdotn,&Qdotnp,&D,&Sc,whicheos,EOSextra);
487 
488 
489 
490 
491 #if(0)
492  // DEBUG:
493  // if(ptrgeom->i==1 && ptrgeom->j==18){
494  DLOOPA(i) dualfprintf(fail_file,"Qtcon[%d]=%21.15g Bcon[%d]=%21.15g Bcov[%d]=%21.15g\n",i,Qtcon[i],i,Bcon[i],i,Bcov[i]);
495  dualfprintf(fail_file,"Bsq=%21.15g QdotB=%21.15g QdotBsq=%21.15g Qtsq=%21.15g Qdotn=%21.15g Qdotnp=%21.15g D=%21.15g Sc=%21.15g\n",Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc);
496  // }
497 #endif
498 
499 
500 
501 
502 
504  //
505  // CHECK WHICH EOM TYPE and proceed accordingly
506  //
508 
509  if(eomtype==EOMFFDE || eomtype==EOMFFDE2){
510  Wdependentsolution=0;
511  if(eomtype==EOMFFDE){
512  FTYPE gammaret0;
513  W_last=0.0;
514  int cleaned=1;
515  if(forcefree_inversion(cleaned, ptrgeom, W_last, Qtcon, Bsq, Bcon, Bcov, Qtsq, U, prim, &gammaret0)) retval++;
516  }
517  else if(eomtype==EOMFFDE2){ // does not so good.
518  int cleaned=0;
519  FTYPE gammaret0;
520 
521  // get Wp, wglobal, etaabs, gamma, gammaabs
522  FTYPE gammasq,gamma,rho0,u,p,wmrho0;
523  FTYPE Wpabs, etaabs;
524  compute_Wp(showmessages, lpflag, eomtype, prim, ptrgeom, &W_last, &Wp_last, &Wpabs, &etaabs, &gamma, &gammasq, &rho0, &u, &p, &wmrho0, wglobal, Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra);
525 
526  if(forcefree_inversion(cleaned, ptrgeom, W_last, Qtcon, Bsq, Bcon, Bcov, Qtsq, U, prim, &gammaret0)) retval++;
527  }
528  else{ // doesn't do as well as EOMFFDE, so smallest instant value of gamma not only important thing, since eventually leads to larger gamma.
529 
530  // get Wp, wglobal, etaabs, gamma, gammaabs
531  FTYPE gammasq,gamma,rho0,u,p,wmrho0;
532  FTYPE Wpabs, etaabs;
533  compute_Wp(showmessages, lpflag, eomtype, prim, ptrgeom, &W_last, &Wp_last, &Wpabs, &etaabs, &gamma, &gammasq, &rho0, &u, &p, &wmrho0, wglobal, Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra);
534 
535  // then perform the analytic inversion with Lorentz factor limit
536  FTYPE gammaret1;
537  FTYPE prim1[NPR];
538  int pliter,pl;
539  PLOOP(pliter,pl) prim1[pl] = prim[pl];
540  int ret1=forcefree_inversion(0, ptrgeom, W_last, Qtcon, Bsq, Bcon, Bcov, Qtsq, U, prim1, &gammaret1);
541 
542  FTYPE gammaret2;
543  FTYPE prim2[NPR];
544  PLOOP(pliter,pl) prim2[pl] = prim[pl];
545  int ret2=forcefree_inversion(1, ptrgeom, W_last, Qtcon, Bsq, Bcon, Bcov, Qtsq, U, prim2, &gammaret2);
546  if(gammaret1>gammaret2){
547  PLOOP(pliter,pl) prim[pl] = prim2[pl];
548  if(ret2) retval++;
549  }
550  else{
551  PLOOP(pliter,pl) prim[pl] = prim1[pl];
552  if(ret1) retval++;
553  }
554  }
555  }
556  else{
557  Wdependentsolution=1;
558 
560  //
561  // SETUP ITERATIVE METHODS (good for GRMHD or entropy GRMHD or cold GRMHD)
562  //
564  retval+=set_guess_Wp(showmessages, lpflag, eomtype, prim, ptrgeom, &W_last, &Wp_last,wglobal, Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra);
565 
566 
568  //
569  //
570  // HOT GRMHD INVERTERS (obtain Wp)
571  //
573  if(eomtype==EOMGRMHD){
574 
575  // METHOD specific:
576 #if(WHICHHOTINVERTER==3)
577  retval+= find_root_1D_gen_Eprime(showmessages, lpflag, eomtype, Wp_last, &Wp, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra,newtonstats);
578 #elif(WHICHHOTINVERTER==2)
579  retval+= find_root_2D_gen(showmessages, lpflag, eomtype, Wp_last, &Wp, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra,newtonstats);
580 #elif(WHICHHOTINVERTER==1)
581  retval+= find_root_1D_gen(showmessages, lpflag, eomtype, Wp_last, &Wp, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra,newtonstats);
582 #else
583  dualfprintf(fail_file,"No such WHICHHOTINVERTER=%d\n",WHICHHOTINVERTER);
584  myexit(89285221);
585 #endif
586 
587 
588 #if(0)
589  retval = find_root_2D_gen(showmessages, lpflag, eomtype, Wp_last, &Wp, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra,newtonstats);
590  retvalother = find_root_1D_gen(showmessages, lpflag, eomtype, Wp_last, &Wpother, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra,newtonstats);
591  retvalother2=find_root_1D_gen_scn(showmessages, lpflag, eomtype, W_last, &W, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra,newtonstats);
592 
593  // Wp=Wpother;
594  if(retval!=retvalother || (fabs(Wp-Wpother)/(Wp+Wpother)>1E-13)){
595  dualfprintf(fail_file,"2D/1D: i=%d j=%d :: %d %d %d :: Wp_last=%21.15g Wp=%21.15g Wpother=%21.15g\n",ptrgeom->i,ptrgeom->j,retval,retvalother,retvalother2,Wp_last,Wp,Wpother);
596  dualfprintf(fail_file,"scn W_last=%g W=%g D=%g Wp_last=%g Wp=%g\n",W_last,W,D,W_last-D,W-D);
597  dualfprintf(fail_file,"rho0=%g u=%g p=%g gamma=%g gammasq=%g w=%g utsq=%g\n",rho0,u,p,gamma,gammasq,w,utsq);
598  }
599 #endif // end comparing SCN and Jon method
600 
601 
602 
603  } // end if hot GRMHD
604 
605 
607  //
608  //
609  // COLD GRMHD INVERTERS (obtain Wp)
610  //
612  else if(eomtype==EOMCOLDGRMHD){
613 
614 
615 
616 #if(WHICHCOLDINVERTER==0)
617  retval0 = find_root_1D_gen(showmessages, lpflag, eomtype, Wp_last, &Wp0, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra,newtonstats);
618 
619 
620  // if(IFUTOPRIMFAIL(*lpflag)){
621  // dualfprintf(fail_file,"E'[W,v^2] equation gave bad answer: Wp=%21.15g Qdotnp=%21.15g Qdotn=%21.15g D=%21.15g\n",Wp0,Qdotnp,Qdotn,D);
622  // *lpflag=0; // let E' try again
623  // }
624  // else{
625  // dualfprintf(fail_file,"E'[W',v^2] equation gave good answer: Wp=%21.15g Qdotnp=%21.15g\n",Wp0,Qdotnp);
626  // }
627 
628  retval+=retval0;
629  Wp=Wp0;
630 
631 
632 
633 
634 #elif(WHICHCOLDINVERTER==1)
635 
636 
637 
638  // Do inversion using E'[W'] equation
639  retval1 = find_root_1D_gen_Eprime(showmessages, lpflag, eomtype, Wp_last, &Wp1, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra,newtonstats);
640 
641 
642 
643 
644 
645  // if(IFUTOPRIMFAIL(*lpflag)){
646  // dualfprintf(fail_file,"E'[W'] equation gave bad answer: Wp=%21.15g Bsq=%21.15g Qdotnp=%21.15g Qdotn=%21.15g D=%21.15g\n",Wp,Bsq,Qdotnp,Qdotn,D);
647 
648  // for(i=0;i<4;i++) Qcovp[i] = U[QCOV0+i] ;
649  // raise_g(Qcovp,ptrgeom->gcon,Qconp) ;
650 
651  // for(i=0;i<4;i++) dualfprintf(fail_file,"%d Qcovp=%21.15g Qconp=%21.15g\n",i,Qcovp[i],Qconp[i]);
652 
653 
654  // *lpflag=0; // let momentum try again
655  // }
656  // else{
657  // dualfprintf(fail_file,"Wp=%21.15g Bsq=%21.15g QdotB=%21.15g QdotBsq=%21.15g Qtsq=%21.15g Qdotnp=%21.15g D=%21.15g\n",Wp,Bsq,QdotB,QdotBsq,Qtsq,Qdotnp,D,Sc);
658  // dualfprintf(fail_file,"E'[W'] equation gave good answer: Wp=%21.15g Qdotnp=%21.15g\n",Wp,Qdotnp);
659  // }
660 
661 
662  retval+=retval1;
663  Wp=Wp1;
664 
665 
666 
667 #elif(WHICHCOLDINVERTER==2)
668  // Do inversion using P^2[W'] equation
669 
670 
671 #if(CRAZYDEBUG&&DEBUGINDEX)
672  if(ifileglobal==0 && jfileglobal==63 && nstep==9 && steppart==2){
673  dualfprintf(fail_file,"Wp_last=%21.15g\n",Wp_last);
674  }
675 #endif
676 
677  retval2 = find_root_1D_gen_Psq(showmessages, lpflag, eomtype, Wp_last, &Wp2, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra,newtonstats);
678  // if(IFUTOPRIMFAIL(*lpflag)){
679  // dualfprintf(fail_file,"P^2[W'] equation gave bad answer: Wp=%21.15g Qdotnp=%21.15g Qdotn=%21.15g\n",Wp2,Qdotnp,Qdotn);
680  // myexit(0);
681  // }
682  // else{
683  // dualfprintf(fail_file,"P^2[W'] equation gave good answer: Wp=%21.15g Qdotnp=%21.15g\n",Wp2,Qdotnp);
684  // }
685 
686 
687  // Wp = (Wp>Wp2) ? Wp2 : Wp;
688  retval+=retval2;
689  Wp=Wp2;
690 
691 
692 
693 #elif(WHICHCOLDINVERTER==3)
694  // Do inversion using P^\alpha[W'] equation
695 
696 
697  retval3 = find_root_3D_gen_Palpha(showmessages, lpflag, eomtype, Wp_last, &Wp3, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra,newtonstats);
698 
699 
700 
701  // if(IFUTOPRIMFAIL(*lpflag)){
702  // dualfprintf(fail_file,"P^2[W'] equation gave bad answer: Wp=%21.15g Qdotnp=%21.15g Qdotn=%21.15g\n",Wp2,Qdotnp,Qdotn);
703  // myexit(0);
704  // }
705  // else{
706  // dualfprintf(fail_file,"P^2[W'] equation gave good answer: Wp=%21.15g Qdotnp=%21.15g\n",Wp2,Qdotnp);
707  // }
708 
709 
710  // Wp = (Wp>Wp3) ? Wp3 : Wp;
711  retval+=retval3;
712  Wp=Wp3;
713 
714 #else
715  dualfprintf(fail_file,"No such WHICHCOLDINVERTER=%d\n",WHICHCOLDINVERTER);
716  myexit(89285224);
717 #endif
718 
719 
720 
721 
722  } // end if eomtype==EOMCOLDGRMHD
723 
724 
725 
727  //
728  //
729  // ENTROPY GRMHD INVERTERS (obtain Wp)
730  //
732  else if(eomtype==EOMENTROPYGRMHD){
733 
734 
735 
736  // METHOD specific:
737 #if(WHICHENTROPYINVERTER==0)
738  retval+= find_root_1D_gen_Sc(showmessages, lpflag, eomtype, Wp_last, &Wp, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra,newtonstats);
739 #else
740  dualfprintf(fail_file,"No such WHICHENTROPYINVERTER=%d\n",WHICHENTROPYINVERTER);
741  myexit(89285225);
742 #endif
743 
744 
745 
746  } // end if-else structure over eomtype
747  }// end if doing non-analytic methods
748 
749 
750 
751 
752 
754  //
755  // Use root if no failure in finding root
756  //
758 
759  // check on result of inversion even if failure (i.e. even if retval>0)
760  // Can be done for Wdependentsolution==0, but ignore Wp
761  check_on_inversion(showmessages, prim, U, ptrgeom, Wp, Qtcon, Bcon, Bcov,retval, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra,newtonstats);
762 
763 
765  if(retval==0 && Wdependentsolution==1){
766  // Check if solution was found
767  retval+=check_Wp(lpflag, eomtype, prim, U, ptrgeom, Wp_last, Wp, retval, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra); // should add in case retval!=0 before this call
768  if(retval) return(retval);
769 
770 
771  // find solution
772  retval+=Wp2prim(showmessages, lpflag, eomtype,prim, pressure, U, ptrgeom, Wp, Qtcon, Bcon, Bcov, retval, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra);
773  if(retval) return(retval);
774  }
775 
776 
777  /* done! */
778  return(retval) ;
779 
780 }
781 
782 
783 
784 
785 
786 
789 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)
790 {
791  int i,j;
792  FTYPE Qcov[4],Qcon[4],Qcovp[4],Qconp[4],ncov[4],ncon[4],Qsq,Qtcov[4];
793  FTYPE alpha,alphasq,Dfactor;
794 
795 
796 
797  // for(i=0;i<NPR;i++) dualfprintf(fail_file,"U[%d]=%g\n",i,U[i]);
798  // exit(0);
799 
800  *D = U[RHO] ;
801 
802 #if(DOENTROPY!=DONOENTROPY)
803  *Sc = U[ENTROPY];
804 #endif
805 
806  // get B^i
807  for(i = BCON1; i <= BCON3; i++) prim[i] = U[i] ;
808  Bcon[0] = 0. ;
809  for(i=1;i<4;i++) Bcon[i] = U[BCON1+i-1] ;
810 
811  // get B_\mu
812  lower_vec(Bcon,ptrgeom,Bcov) ;
813 
814  // get B^2
815  *Bsq = 0. ;
816  for(i=1;i<4;i++) *Bsq += Bcon[i]*Bcov[i] ;
817 
818  // get \eta^\mu and \eta_\mu
819  ncov_calc_fromlapse(ptrgeom->alphalapse,ncov) ;
820  raise_vec(ncov,ptrgeom,ncon);
821 
822  // -P_\alpha (P_t has D term inside it)
823  // Q_\nu = \alpha \{T^t_\nu\} = -T.\eta
824 
825  // Q'_\nu = Q_\nu + \delta^t_\nu D
826  for(i=0;i<4;i++) Qcovp[i] = U[QCOV0+i] ;
827  raise_vec(Qcovp,ptrgeom,Qconp) ;
828 
829 
830 
831  // P^2, which can be written without use of T^t_t
832  // Qsq = 0. ;
833  //for(i=0;i<4;i++) Qsq += Qcovp[i]*Qconp[i] ;
834 
835  // P^2 = Qtilde^2 can be written without T^t_t
836  // Notice U !=Qcon, but only a problem for time term that is not used
837  // T^t_t not really needed (cancels in the end)
838  // DLOOP(j) Qtcov[j]=0.0;
839  // DSLOOP(j,i) Qtcov[j] = U[QCOV0+i]*(delta(i,j)+ncon[i]*ncov[j]) ;
840  // raise_vec(Qtcov,ptrgeom,Qtcon) ;
841 
842  // P^\alpha = Qtcon^\alpha
843  DLOOPA(j) Qtcon[j]=0.0;
844  SSLOOP(j,i) Qtcon[j] += U[QCOV0+i]*(ptrgeom->gcon[GIND(i,j)] + ncon[i]*ncon[j]) ;
845  // for(i=0;i<4;i++){
846  // for(j=1;j<4;j++){
847  // Qtcon[i] = U[QCOV0+j]*(ptrgeom->gcon[GIND(j,i)]+ncon[j]*ncon[i]);
848  // }
849  // }
850  lower_vec(Qtcon,ptrgeom,Qtcov);
851  *Qtsq = 0. ;
852  for(i=0;i<4;i++) *Qtsq += Qtcov[i]*Qtcon[i] ;
853 
854 #if(0)
855  DLOOPA(j) dualfprintf(fail_file,"Qtcon[%d]=%21.15g Qtcov[%d]=%21.15g\n",j,Qtcon[j],j,Qtcov[j]);
856 #endif
857 
858 
859  // Qcon/Qcov = -\eta.T = \alpha T^t_\nu :: where T^t_t *includes* rest-mass
860  // avoid using this Q if possible to avoid large errors when v or u is small
861  // put back in the D
862  for(i=0;i<4;i++){
863  if(i==0) Qcov[i]=Qcovp[i] - (*D);
864  else Qcov[i]=Qcovp[i];
865  }
866  raise_vec(Qcov,ptrgeom,Qcon) ;
867 
868 
869 
870 
871 
872  // dualfprintf(fail_file,"D=%g\n",D);
873  // for(i=0;i<4;i++) dualfprintf(fail_file,"Qcon[%d]=%g Qcov[%d]=%g\n",i,Qcon[i],i,Qcov[i]);
874  // exit(0);
875 
876 
877 
878  // P.B (P_t not used because B^t=0, so ok to use Qcov)
879  *QdotB = 0. ;
880  for(i=1;i<4;i++) *QdotB += Qcov[i]*Bcon[i] ;
881  *QdotBsq = (*QdotB)*(*QdotB) ;
882 
883  // Energy
884  // Qdotnp=-E'=-E+D automatically
885 #if(0)
886  // This does not resolve metric term to have gravity taken into account in non-rel limit
887  // Note that gravity will be de-resolved once $r\sim 10^{14}M$
888  *Qdotnp = Qconp[0]*ncov[0] + (*D)*(1.0-1.0/fabs(ncov[TT])) ; // -Qdotn-W = -Qdotnp-Wp
889 #else
890 
891  //
892  // Get 1-1/|\eta_t| in way that avoids catastrophic cancellation in non-rel gravity limit
893  //
895  // 1-1/|\eta_t| = 1-1/\alpha = (\alpha-1)/\alpha = (\alpha^2-1)/(\alpha(\alpha+1)) = (-gcovpert[TT] + \alpha^2 g_{ti} g^{ti})/(\alpha(\alpha+1))
896  //
897  // \beta^i \beta_i / \alpha^2 = g^{ti} g_{ti}
898  // betasqoalphasq = 0.0;
899  // SLOOPA(j) betasqoalphasq += (ptrgeom->gcov[GIND(TT,j)])*(ptrgeom->gcon[GIND(TT,j)]);
900 
901  alpha = ptrgeom->alphalapse;
902  alphasq = alpha*alpha;
903 
904  Dfactor = (-ptrgeom->gcovpert[TT] + alphasq*(ptrgeom->betasqoalphasq))/(alphasq+alpha);
905 
906  // now Compute Qdotnp
907  *Qdotnp = Qconp[0]*ncov[0] + (*D)*(Dfactor) ; // -Qdotn-W = -Qdotnp-Wp
908 #endif
909 
910 
911  *Qdotn = Qcon[0]*ncov[0] ;
912 
913 
914 #if(0)
915  // deprecated DEBUG:
916  // used to check on old vs. new method
917  *Qtsqnew=*Qtsq;
918  // dualfprintf(fail_file,"1: Qtsq=%21.15g\n",Qtsq);
919  for(i=0;i<4;i++) Qcovorig[i] = U[QCOV0+i] ;
920  Qcovorig[TT] -=D;
921  raise_vec(Qcovorig,ptrgeom,Qconorig) ;
922  Qdotnorig = Qconorig[0]*ncov[0] ;
923  *Qsqorig = 0. ;
924  for(i=0;i<4;i++) *Qsqorig += Qcovorig[i]*Qconorig[i] ;
925  *Qtsqorig = *Qsqorig + Qdotnorig*Qdotnorig ;
926  // dualfprintf(fail_file,"2: Qtsq=%21.15g\n",*Qtsq);
927  *Qtsq=*Qtsqorig;
928 #endif
929 
930 
931  // P^2 (since D already removed from E, then cancelling terms are small)
932  // so no problem with catastrophic cancellation since involves energy term
933  // *Qtsq = *Qsq + (*Qdotnp)*(*Qdotnp) ; // P^2 in eta frame
934 
935 
936  // dualfprintf(fail_file,"Qdotn=%g Qtsq=%g\n",*Qdotn,*Qtsq);
937  //exit(0);
938 
939 
940 #if(CRAZYDEBUG)
941  // debug Qsq
942  // deprecated DEBUG:
943  if(ptrgeom->i==0 && ptrgeom->j==31 && nstep==9 && steppart==2){
944  Qsq = 0. ;
945  for(i=0;i<4;i++) *Qsq += Qcov[i]*Qcon[i] ;
946  for(i=0;i<4;i++) dualfprintf(fail_file,"Qcon[%d]=%21.15g Qcov[%d]=%21.15g Qsq=%21.15g Qtsq=%21.15g\n",i,Qcon[i],i,Qcov[i],Qsq,Qtsq);
947  }
948 #endif
949 
950 
951  return(0);
952 
953 }
954 
955 
958 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)
959 {
960  FTYPE utsq;
961  FTYPE Ss;
962  FTYPE w;
963  int i,j;
964  int jj;
965  FTYPE Bcon[NDIM],Bcov[NDIM],bsq;
966  FTYPE nuabs;
967 
968  /* calculate W from last timestep and use
969  for guess */
970  utsq = 0. ;
971  for(i=1;i<4;i++)
972  for(j=1;j<4;j++) utsq += ptrgeom->gcov[GIND(i,j)]*prim[UTCON1+i-1]*prim[UTCON1+j-1] ;
973 
974 
975  if( (utsq < 0.) && (fabs(utsq) < MAXNEGUTSQ) ) {
976  if(debugfail>=2) dualfprintf(fail_file,"utsq was %21.15g now 0.0 @ i=%d j=%d k=%d\n",utsq,ptrgeom->i,ptrgeom->j,ptrgeom->k);
977  utsq = 0.0;
978  }
979  if(utsq >= 0. && utsq < UTSQ_TOO_BIG) {
980  // then good, if nan then fails with next else
981  }
982  else{
983  if( debugfail>=2 ){
984  dualfprintf(fail_file,"Utoprim_new(): utsq < 0 or utsq>UTSQ_TOO_BIG in utoprim_jon attempt. Resetting to utsq=0.0 from utsq = %21.15g @ i=%d j=%d k=%d\n", utsq,ptrgeom->i,ptrgeom->j,ptrgeom->k) ;
985  for(i=1;i<4;i++)
986  for(j=1;j<4;j++) dualfprintf(fail_file,"i=%d j=%d :: gcov=%21.15g primi=%21.15g primj=%21.15g\n",i,j,ptrgeom->gcov[GIND(i,j)],prim[UTCON1+i-1],prim[UTCON1+j-1]) ;
987  }
988 
989  // then set some kind of guess since guess may just be bad and may still be able to find a solution
990  utsq = 0.0;
991 
992  // *lpflag= UTOPRIMFAILCONVGUESSUTSQ; // guess failure actually
993  // check on v^2>1 failure
994  //check_utsq_fail(-1E30, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra);
995  //return(1) ;
996  }
997 
998  *gammasq = 1. + utsq ;
999  *gamma = sqrt(*gammasq);
1000 
1001  // Always calculate rho from D and gamma so that using D in EOS remains consistent
1002  // i.e. you don't get positive values for dP/d(vsq) .
1003  *rho0 = (D) / (*gamma) ;
1004  *u = prim[UU] ;
1005  if(whicheos==KAZFULL){
1006  *u = MAX(0.0,*u); // near degeneracy, allow somewhat "hot" solution so guess doesn't give immediately bad Wp and utsq(Wp)
1007  }
1008  *p = pressure_rho0_u(whicheos,EOSextra,*rho0,*u) ; // p(rho0,u) just used for guess, while p(rho0,\chi) used to get solution since assume don't know \chi yet. Even if had initial \chi, wouldn't be final \chi anyways.
1009  *wmrho0=(*u)+(*p);
1010  w = (*rho0) + (*wmrho0) ;
1011 
1012 
1013 
1014 
1015  // W'=W-D (D removed from W)
1016  *Wp_last = (D)*utsq/(1.0+(*gamma)) + ((*u) + (*p))*(*gammasq);
1017  *W_last = *Wp_last + (D);
1018 
1019 
1020  // need b^2 to normalize W since error in W limited by b^2
1021  // GODMARK: Computing b^2 this way is too expensive, just use B^2 since just used to make an error estimate
1022  // Misses u.B term
1023  // bsq_calc(prim, ptrgeom, &bsq);
1024  Bcon[TT]=0.0; SLOOPA(jj) Bcon[jj]=prim[B1-1+jj];
1025  lower_vec(Bcon,ptrgeom,Bcov);
1026  bsq = 0.0; SLOOPA(jj) bsq += Bcon[jj]*Bcov[jj];
1027  bsq *= ((ptrgeom->alphalapse * ptrgeom->alphalapse))/(*gammasq);
1028 
1029  // Get normalizations
1030  *Wpabs = fabs((D)*utsq)/(1.0+fabs(*gamma)) + (fabs(*u)+fabs(*p))*fabs(*gammasq) + SQRTMINNUMREPRESENT;
1031  nuabs = fabs(*u) + fabs(*p) + fabs(bsq);
1032  *etaabs = fabs(*rho0) + nuabs;
1033  wglobal[2]=GAMMASQCHECKRESID * fabs(*rho0);
1034 
1035 
1036  return(0);
1037 }
1038 
1039 
1046 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)
1047 {
1048  FTYPE Ss;
1049  FTYPE gammasq,gamma,rho0,u,p,wmrho0;
1050  FTYPE Wpabs, etaabs;
1051  int i,j;
1052  FTYPE Wpcoldhd;
1053  int numattemptstofixguess;
1054  int jj;
1055 
1056 
1057  // get Wp, wglobal, etaabs, gamma, gammaabs
1058  compute_Wp(showmessages, lpflag, eomtype, prim, ptrgeom, W_last, Wp_last, &Wpabs, &etaabs, &gamma, &gammasq, &rho0, &u, &p, &wmrho0, wglobal, Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra);
1059 
1060 
1061 
1062 
1063  // dualfprintf(fail_file,"Wlast=%g Wplast=%g\n",*W_last,*Wp_last);
1064 
1065  // *wglobal=w+bsq; // used to validate normalized version of W
1066 
1068  //SUPERGODMARK: what if in nonrel case rho0 is identically zero? looks like rel. case but it's not (need rhorel for non-rel case)
1069  if( JONGENNEWTSTUFF && gammasq * etaabs >= (wglobal[2]) ){
1070  // RIGHT
1071 
1072 
1074  // if( JONGENNEWTSTUFF && (w+bsq)*gammasq >= GAMMASQCHECKRESID ) { // WRONG
1076 
1077  int gotcold=0;
1078  if(coldgrhd(lpflag, Qtsq, D, &Wpcoldhd)){
1079  Wpcoldhd=*Wp_last; // if failed to get coldGRHD Wp, then just use Wp_last
1080  gotcold=1;
1081  }
1082 
1083  // *wglobal=MAX(MAX(Wpcoldhd,*Wp_last),Wpabs); // used to validate normalized version of W -- accurate for cold GRHD
1084  *wglobal= MAX(MAX(Wpcoldhd,*Wp_last),Wpabs); // used to validate normalized version of W -- accurate for cold GRHD
1085  //*wglobal=Wpabs;
1086 
1087  if(gotcold) *Wp_last=MAX(2.0*Wpcoldhd,*Wp_last); // use cold HD to upgrade Wp if needed. Don't upgrade from cold if just set to Wp_last from guess.
1088 
1089  // fudge
1090  // *Wp_last = *W_last;
1091  }
1092  else{
1093  *wglobal=Wpabs; // non-relativistically correct
1094  }
1095 
1096  // include Wp-like term with rest-mass so becomes W-like term
1097  wglobal[1]=fabs(rho0)+wglobal[0]+SQRTMINNUMREPRESENT;
1098 
1099 
1100  // if((nstep==4)&&(ptrgeom->i==0)&&(ptrgeom->j==47)){
1101  // dualfprintf(fail_file,"utsq=%g p=%g w=%g W_last=%g Wp_last=%g\n",utsq,p,w,*W_last,*Wp_last);
1102  // exit(0);
1103  // }
1104 
1105 #if(CRAZYDEBUG&&DEBUGINDEX)
1106  if(ifileglobal==0 && jfileglobal==63 && nstep==9 && steppart==2){
1107  dualfprintf(fail_file,"pre verify\n");
1108  dualfprintf(fail_file,"rho0=%21.15g u=%21.15g p=%21.15g\n",rho0,u,p);
1109 
1110  dualfprintf(fail_file,"W_last=%21.15g Wp_last=%21.15g gammasq=%21.15g w=%21.15g\n",*W_last,*Wp_last,gammasq,w);
1111  }
1112 #endif
1113 
1114 
1115 
1116  verify_Wlast(u, p, ptrgeom, W_last, Wp_last, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra); // W_last and Wp_last and rest of variables are already pointers
1117 
1118 #if(CRAZYDEBUG&&DEBUGINDEX)
1119  if(ifileglobal==0 && jfileglobal==63 && nstep==9 && steppart==2){
1120  dualfprintf(fail_file,"post verify\n");
1121 
1122  dualfprintf(fail_file,"W_last=%21.15g Wp_last=%21.15g gammasq=%21.15g w=%21.15g\n",*W_last,*Wp_last,gammasq,w);
1123  }
1124 #endif
1125 
1126 
1127 
1128 
1129 
1130  // dualfprintf(fail_file,"MID: Wlast=%g Wplast=%g\n",*W_last,*Wp_last);
1131 
1132 
1133 #define MAXNUMGUESSCHANGES (1000)
1134 
1135 
1136  // sometimes above gives invalid guess (Wp=0 or utsq<0) so fix
1137  numattemptstofixguess=0;
1138  FTYPE utsq,Ss0;
1139  while(1){
1140  // under all eomtype's, ensure utsq reasonable for guess so Newton's method starts at reasonable value around which to work from
1141  // check utsq from this guess for W
1142  utsq=utsq_calc(*W_last, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra); // check for precision problems
1143 
1144  if(eomtype==EOMENTROPYGRMHD){
1145  // for Entropy method, must also ensure Ss is reasonable
1146  Ss=Ss_Wp_utsq(whicheos,EOSextra,*Wp_last,D,utsq);
1147  // see what initial entropy is. Even if Ss is not a nan or inf, can be so small that Newton stepping will be wrong, so keep increasing Wp until actually hotter than older entropy. This ensures start hotter than required and not bumping-up against limit of method.
1148  // Noticed that entropy inversion will think it succeeds (low errx and dx) but it really fails very catastrophically! (conserved quantities don't match) for this reason. Typically would give very low internal energy result, when should have been much higher. Doesn't do that anymore.
1149  // But, once that very low entropy occurs, then noticed a secondary bad effect. The updated conserved quantity would then be VERY large, leading to a VERY large entropy and huge internal energy. While this is avoided with no internal energy drop-outs, still a concern that very low entropy regions can cause problematic fluxes and give high entropy. Probably need larger floor on entropy than SMALL in EOS.
1150  Ss0=compute_specificentropy_wmrho0(whicheos,EOSextra,rho0,wmrho0);
1151  }
1152  else{
1153  // just some dummy value so passes below conditional without affect
1154  Ss=SMALL;
1155  Ss0=0.0;
1156  }
1157 
1158 #define FRACSs0 (0.5) // want to have initial entropy not too far below previous entropy value
1159  if(utsq>=0.0 && utsq==utsq && isfinite(utsq) && (Ss==Ss && isfinite(Ss) && Ss>=Ss0-FRACSs0*fabs(Ss0))){
1160  // if utsq=nan or inf, will fail to reach here
1161  // if Ss=nan or inf, will fail to reach here
1162  if(numattemptstofixguess>0) if(showmessages && debugfail>=3) dualfprintf(fail_file,"GOOD Initial guess #%d/%d [i=%d j=%d k=%d] for W=%21.15g Wp=%21.15g Wp/D=%21.15g gives bad utsq=%21.15g Ss=%21.15g D=%21.15g u=%21.15g p=%21.15g gamma=%21.15g Ss0=%21.15g\n",numattemptstofixguess,MAXNUMGUESSCHANGES,ptrgeom->i,ptrgeom->j,ptrgeom->k,*W_last,*Wp_last,*Wp_last/D,utsq,Ss,D,u,p,gamma,Ss0);
1163  break;
1164  }
1165  else{
1166 #if(PRODUCTION==0)
1167  if(showmessages && debugfail>=3) dualfprintf(fail_file,"Initial guess #%d/%d [i=%d j=%d k=%d] for W=%21.15g Wp=%21.15g Wp/D=%21.15g gives bad utsq=%21.15g Ss=%21.15g D=%21.15g u=%21.15g p=%21.15g gamma=%21.15g Ss0=%21.15g\n",numattemptstofixguess,MAXNUMGUESSCHANGES,ptrgeom->i,ptrgeom->j,ptrgeom->k,*W_last,*Wp_last,*Wp_last/D,utsq,Ss,D,u,p,gamma,Ss0);
1168 #endif
1169 
1170 
1171  PFTYPE coldpflag;
1172  FTYPE coldWp;
1173  coldgrhd(&coldpflag, Qtsq, D, &coldWp);
1174  if(*Wp_last<coldWp) *Wp_last=coldWp;
1175 
1176 
1177  if(numattemptstofixguess>0){
1178  // then cold fix was not enough:
1179  //#define FACTORSHIFT (10.0) // bit much
1180 #define FACTORSHIFT (2.0)
1181  *Wp_last = MAX(MAX(fabs(*Wp_last)*FACTORSHIFT,NUMEPSILON*100.0*fabs(D)),KINDASMALL);
1182  }
1183 
1184  // set new W
1185  *W_last = *Wp_last + (D);
1186  }
1187 
1188 
1189  if(numattemptstofixguess>MAXNUMGUESSCHANGES){
1190  *lpflag= UTOPRIMFAILCONVGUESSUTSQ; // guess failure actually
1191  // check on v^2>1 failure
1192  return(1) ;
1193  }
1194 
1195 
1196  numattemptstofixguess++;
1197 
1198  }
1199 
1200 
1201  // dualfprintf(fail_file,"NEW: Wlast=%g Wplast=%g\n",*W_last,*Wp_last);
1202 
1203 
1204  // DEBUG:
1205  // dualfprintf(fail_file,"DEBUG: %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g\n",utsq,rho0,u,p,bsq,*Wp_last,wglobal[2]);
1206  //int pl; PALLLOOP(pl) dualfprintf(fail_file,"DEBUG2: prim[%d]=%21.15g\n",pl,prim[pl]);
1207  // int jjm,kkm; DLOOP(jjm,kkm) dualfprintf(fail_file,"DEBUG3: gcov=%21.15g gcon=%21.15g\n",ptrgeom->gcov[GIND(jjm,kkm)],ptrgeom->gcon[GIND(jjm,kkm)]);
1208  //for(pl=FIRSTEOSGLOBAL;pl<=LASTEOSGLOBAL;pl++) dualfprintf(fail_file,"DEBUG3: EOSextra[%d]=%21.15g\n",pl,EOSextra[pl]);
1209 
1210 
1211 
1212 
1213 
1214  return(0);
1215 
1216 
1217 }
1218 
1219 
1229 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)
1230 {
1231  int i_increase;
1232 
1233 
1234  // FTYPE Wpnewsmall=-3.0*MAX(MAX(MAX(MAX(MAX(fabs(D),fabs(Qdotn)),fabs(Bsq)),fabs(Qdotnp)),fabs(wglobal[0])),fabs(wglobal[1]));
1235  // allow this constraint on W_last and Wp_last so starts hotter than could be.
1236 
1237  if(*W_last<-(D) || *Wp_last<0){
1238  *W_last = -(D);
1239  *Wp_last = 0;
1240 #if(DEBUGINDEX)
1241  if(ifileglobal==0 && jfileglobal==63 && nstep==9 && steppart==2){
1242  dualfprintf(fail_file,"W_last=%21.15g Wp_last=%21.15g\n",*W_last,*Wp_last);
1243  }
1244 #endif
1245  }
1246 
1247 
1248  if(1){
1249  i_increase = 0;
1250  while( vsqgtr1(*W_last,Bsq,QdotBsq,Qtsq) && (i_increase < 10) ) {
1251 
1252  // dualfprintf(fail_file,"bumped up %d %d %g %g\n",ptrgeom->i,ptrgeom->j,*W_last,*Wp_last);
1253 
1254 
1255  *Wp_last*= 10.;
1256  *W_last = *Wp_last+(D);
1257 
1258  // dualfprintf(fail_file,"bumped up %d %d\n",ptrgeom->i,ptrgeom->j);
1259 
1260 
1261  i_increase++;
1262 #if(!OPTIMIZED)
1263  dualfprintf(fail_file,"badval : W = %21.15g, i_increase = %d \n", *W_last, i_increase);
1264 #endif
1265  }
1266 #if(!OPTIMIZED)
1267  if( i_increase >= 10 ) {
1268  dualfprintf(fail_file,"i_increase is too large, i_increase = %d , W = %21.15g \n", i_increase, *W_last);
1269  }
1270 #endif
1271 
1272 
1273 #if(!OPTIMIZED)
1274  if( ltrace ) {
1275  dualfprintf(fail_file,"u = %21.15g, p = %21.15g, Bsq = %21.15g, Qtsq = %21.15g \n",u,p,Bsq,Qtsq);
1276  dualfprintf(fail_file,"Bcon[0-3] = %21.15g %21.15g %21.15g %21.15g \n", Bcon[0],Bcon[1],Bcon[2],Bcon[3]);
1277  dualfprintf(fail_file,"Bcov[0-3] = %21.15g %21.15g %21.15g %21.15g \n", Bcov[0],Bcov[1],Bcov[2],Bcov[3]);
1278  dualfprintf(fail_file,"Qcon[0-3] = %21.15g %21.15g %21.15g %21.15g \n", Qcon[0],Qcon[1],Qcon[2],Qcon[3]);
1279  dualfprintf(fail_file,"Qcov[0-3] = %21.15g %21.15g %21.15g %21.15g \n", Qcov[0],Qcov[1],Qcov[2],Qcov[3]);
1280  dualfprintf(fail_file,"call find_root\n") ;
1281 
1282  }
1283 #endif
1284  }
1285 
1286 
1287  return(0);
1288 
1289 }
1290 
1291 
1292 
1293 
1299 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)
1300 {
1301 
1302  FTYPE Wtest;
1303 
1304 
1305 
1306 
1307  /* Problem with solver, so return denoting error before doing anything further */
1308  if( (retval != 0) || (Wp == FAIL_VAL) ) {
1309  if( debugfail>=2 ) {
1310  dualfprintf(fail_file, "Failed to find a prim. var. solution!! Wp_last=%21.15g Bsq=%21.15g QdotBsq=%21.15g Qdotn=%21.15g D=%21.15g Qtsq=%21.15g \n",Wp_last,Bsq,QdotBsq,Qdotn,D,Qtsq);
1311  dualfprintf(fail_file, "Utoprim_new_body(): bad newt failure, nstep,steppart :: t,i,j,k, p[0-%d], U[0-%d] = %ld %d :: %21.15g %d %d %d ", NPR, NPR, nstep,steppart,t, ptrgeom->i, ptrgeom->j,ptrgeom->k );
1312  int pliter,pl;
1313  PINVERTLOOP(pliter,pl){
1314  dualfprintf(fail_file, "%21.15g ", prim[pl]);
1315  }
1316  PINVERTLOOP(pliter,pl){
1317  dualfprintf(fail_file, "%21.15g ", U[pl]);
1318  }
1319  dualfprintf(fail_file, "\n");
1320  }
1321  *lpflag= retval+UTOPRIMFAILCONVRET+1;// related to UTOPRIMFAILCONVRET
1322  if(debugfail>=1) dualfprintf(fail_file,"Got here Wp_last=%21.15g Wp=%21.15g retval=%d\n",Wp_last,Wp,retval);
1323  return(retval);
1324  }
1325  else{
1326  // Note: Wp comparable to D in scale/dimension
1327  // Note: Wp/wglobal[1] comparable to \gamma^2
1328  Wtest=Wp;
1329 
1330  // GODMARK
1331  // if(Wtest<=0. && fabs(Wtest)>1E-4){
1332  // Complex result if W+D<0 or Wp-D+D<0 or Wp<0
1333  FTYPE Wpnewsmall=-3.0*MAX(MAX(MAX(MAX(MAX(fabs(D),fabs(Qdotn)),fabs(Bsq)),fabs(Qdotnp)),fabs(wglobal[0])),fabs(wglobal[1]));
1334  // if(Wtest<=-(D)){
1335  if(Wtest<Wpnewsmall){
1336  if(debugfail>=1){
1337  dualfprintf(fail_file,"Wtest1 failed :: Wp=%21.15g D=%21.15g wglobal=%21.15g\n",Wp,(D),wglobal[1]);
1338  dualfprintf(fail_file,"nstep=%ld steppart=%d :: t=%21.15g :: i=%d j=%d k=%d p=%d\n",nstep,steppart,t,ptrgeom->i,ptrgeom->j,ptrgeom->k,ptrgeom->p);
1339  }
1340  *lpflag= retval+UTOPRIMFAILCONVRET+1;// related to UTOPRIMFAILCONVRET
1341  // reset Wp
1342  // Wp = -(D)+fabs((D))*0.01;
1343  Wp = Wpnewsmall;
1344  //
1345  return(retval) ;
1346  }
1347  else if(Wtest < Wpnewsmall || Wtest/wglobal[1] > GAMMASQ_TOO_BIG) {
1348  // dualfprintf(fail_file,"Wtest2 failed Wp=%21.15g\n",Wp);
1349 
1350 
1351 
1352  //if(Wtest <= 0.) {
1353  if( debugfail>=2 ) {
1354  dualfprintf(fail_file,"Wtest2 failure %21.15g %21.15g %21.15g %21.15g\n",Wtest,D,wglobal[1],GAMMASQ_TOO_BIG) ;
1355  dualfprintf(fail_file, "Utoprim_new_body(): Wtest<Wpnewsmall or Wtest=toobig failure, t,i,j,k, p[0-7], U[0-7] = %21.15g %d %d %d ", t, ptrgeom->i, ptrgeom->j,ptrgeom->k );
1356  int pliter,pl;
1357  PINVERTLOOP(pliter,pl){
1358  dualfprintf(fail_file, "%21.15g ", prim[pl]);
1359  }
1360  PINVERTLOOP(pliter,pl){
1361  if(isfinite(U[pl])) dualfprintf(fail_file, "%21.15g ", U[pl]);
1362  else dualfprintf(fail_file, "nan ");
1363  }
1364  dualfprintf(fail_file, "\n");
1365  }
1366 
1367  *lpflag= UTOPRIMFAILCONVW;
1368  return(retval) ;
1369  }
1370  }
1371 
1372 
1373 #if(!OPTIMIZED)
1374  if( ltrace ) {
1375  dualfprintf(fail_file,"(Wp,Wp_last,Bsq,Qtsq,QdotB,gammasq,Qdotn) %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g\n",
1376  Wp,Wp_last,
1377  Bsq,Qtsq,QdotB,gammasq,Qdotn) ;
1378  dualfprintf(fail_file,"done find_root\n") ;
1379  }
1380  if( ltrace2 ) {
1381  dualfprintf(fail_file, "\n <--------- %21.15g %21.15g %21.15g %21.15g %21.15g \n", Bsq,QdotBsq,Qdotn,D,Qtsq);
1382 
1383  }
1384 #endif
1385 
1386  /*
1387  if( nstep == 3 ) {
1388  if( ((ptrgeom->i==89) && (ptrgeom->j==88||ptrgeom->j==111)) || (ptrgeom->i==92&&(ptrgeom->j==86||ptrgeom->j==113))
1389  || (ptrgeom->i==92&&(ptrgeom->j==86||ptrgeom->j==113) ) || (ptrgeom->i==94&&(ptrgeom->j==84||ptrgeom->j==115) )
1390  || (ptrgeom->i==105&&(ptrgeom->j==84||ptrgeom->j==115) ) || (ptrgeom->i==107&&(ptrgeom->j==86||ptrgeom->j==113) )
1391  || (ptrgeom->i==110&&(ptrgeom->j==88||ptrgeom->j==111) ) ) {
1392  dualfprintf(fail_file, "\n <--------- %21.15g %21.15g %21.15g %21.15g %21.15g \n", Bsq,QdotBsq,Qdotn,D,Qtsq);
1393  }
1394  }
1395  */
1396 
1397  return(0);
1398 
1399 
1400 }
1401 
1402 
1404 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)
1405 {
1406 
1407 #if(DEBUGINDEX)
1408  if(steppart==1 && nstep==32 && ifileglobal==0 && jfileglobal==19){
1409  dualfprintf(fail_file,"Qtsq=%21.15g Bsq=%21.15g D=%21.15g QdotB=%21.15g Wp=%21.15g\n",Qtsq,Bsq,D,QdotB,Wp);
1410  myexit(12245);
1411  }
1412 #endif
1413 
1414 }
1415 
1416 
1417 static void check_on_inversion(int showmessages, FTYPE *prim, 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, struct of_newtonstats *newtonstats)
1418 {
1419  int checki;
1420 
1421 #if(CHECKONINVERSION)
1422  checki=0;
1423  strcpy(newtonstats->invpropertytext[checki++],"Qdotnp");
1424  strcpy(newtonstats->invpropertytext[checki++],"Qtsq");
1425  // strcpy(newtonstats->invpropertytext[checki++],"Qtsqorig");
1426  strcpy(newtonstats->invpropertytext[checki++],"Bsq");
1427  strcpy(newtonstats->invpropertytext[checki++],"DD");
1428  strcpy(newtonstats->invpropertytext[checki++],"Sc");
1429  strcpy(newtonstats->invpropertytext[checki++],"QdotB");
1430  strcpy(newtonstats->invpropertytext[checki++],"WWp");
1431  strcpy(newtonstats->invpropertytext[checki++],"Qtcon1");
1432  strcpy(newtonstats->invpropertytext[checki++],"Qtcon2");
1433  strcpy(newtonstats->invpropertytext[checki++],"Qtcon3");
1434  strcpy(newtonstats->invpropertytext[checki++],"Bcon1");
1435  strcpy(newtonstats->invpropertytext[checki++],"Bcon2");
1436  strcpy(newtonstats->invpropertytext[checki++],"Bcon3");
1437  // then store quantities required to check up on inversion in mathematica
1438  checki=0;
1439  newtonstats->invproperty[checki++]=Qdotnp;
1440  newtonstats->invproperty[checki++]=Qtsq;
1441  // newtonstats->invproperty[checki++]=Qtsqorig;
1442  // newtonstats->invproperty[checki++]=Qtsqnew;
1443  newtonstats->invproperty[checki++]=Bsq;
1444  newtonstats->invproperty[checki++]=D;
1445  newtonstats->invproperty[checki++]=Sc;
1446  newtonstats->invproperty[checki++]=QdotB;
1447  newtonstats->invproperty[checki++]=Wp;
1448  newtonstats->invproperty[checki++]=Qtcon[1];
1449  newtonstats->invproperty[checki++]=Qtcon[2];
1450  newtonstats->invproperty[checki++]=Qtcon[3];
1451  newtonstats->invproperty[checki++]=Bcon[1];
1452  newtonstats->invproperty[checki++]=Bcon[2];
1453  newtonstats->invproperty[checki++]=Bcon[3];
1454  if(checki>NUMINVPROPERTY){
1455  dualfprintf(fail_file,"Not enough memory for invproperty: checki=%d NUMINVPROPERTY=%d\n",checki,NUMINVPROPERTY);
1456  }
1457 #endif
1458 }
1459 
1460 
1461 
1462 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)
1463 {
1464  FTYPE W,vsq;
1465  FTYPE utsq;
1466  int i;
1467  FTYPE rho0;
1468  FTYPE gtmp,gamma,gammasq,w,wmrho0,p,u,tmpdiff;
1469  FTYPE wmrho0_compute_vsq(FTYPE Wp, FTYPE D, FTYPE vsq, FTYPE gamma,FTYPE gammasq);
1470  FTYPE wmrho0_compute_utsq(FTYPE Wp, FTYPE D, FTYPE utsq, FTYPE gamma,FTYPE gammasq);
1471  int retvsq;
1472  FTYPE pr0[NPR];
1473  int pl,pl2;
1474  int foundnan;
1475 
1476 #if(DEBUGINDEX)
1477  // if(ifileglobal==1 && jfileglobal==10 && nstep==137 && steppart==3){
1478  //Wp = -0.513243850397847;
1479  // Wp = -0.4043038500443759;
1480  //Wp = -0.0010560020881852734;
1481  //Wp = 0.001055708263718976;
1482  //dualfprintf(fail_file,"Testing\n");
1483  //}
1484 #endif
1485 
1486 #if(!OPTIMIZED)
1487  if(!isfinite(Wp)){
1488  PLOOPALLINVERT(pl) prim[pl]=pr0[pl]; // since might contaminate solution
1489  *lpflag= UTOPRIMFAILCONVW;
1490  return(1);
1491  }
1492 #endif
1493 
1494 
1495 
1496 
1497 
1499  //
1500  // Invert (trivially) field components
1501  //
1503  for(i = BCON1; i <= BCON3; i++) prim[i] = U[i] ;
1504 
1506  //
1507  // set backup of primitive
1508  //
1510  PLOOPALLINVERT(pl) pr0[pl]=prim[pl];
1511 
1512 
1514  //
1515  // Set dependent quantity W[W']
1516  //
1518  W=Wp+D;
1519 
1520 
1521 
1522 
1524  //
1525  // Compute:
1526  // 1) Either v^2 or \tilde{u}^2
1527  // 2) \gamma
1528  // 3) \chi = u+p
1529  // 4) \rho_0
1530  // 5) u
1531  //
1533 
1534 #if(PRIMFROMVSQ==1)
1535  // Calculate vsq
1536  vsq = vsq_calc(W, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra) ;
1537 
1538  retvsq=validate_vsq(vsq,&vsq, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra);
1539 
1540 
1541  if( retvsq==1 || retvsq==2) {
1542  //(vsq >= 1.0) || (vsq<0.0) ) {
1543  if( debugfail>=2 ) {
1544  dualfprintf(fail_file,"vsq failure: vsq = %21.15g , W = %21.15g \n",vsq, W) ;
1545  dualfprintf(fail_file, "Utoprim_new_body(): utsq==bad failure, t,i,j,k, p[0-7], U[0-7] = %21.15g %d %d %d ", t, ptrgeom->i, ptrgeom->j,ptrgeom->k );
1546  PINVERTLOOP(pliter,pl){
1547  dualfprintf(fail_file, "%21.15g ", prim[pl]);
1548  }
1549  PINVERTLOOP(pliter,pl){
1550  dualfprintf(fail_file, "%21.15g ", U[pl]);
1551  }
1552  dualfprintf(fail_file, "\n");
1553  }
1554 
1555  // check on v^2>1 failure
1556  check_utsq_fail(Wp, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra);
1557 
1558  *lpflag= UTOPRIMFAILCONVUTSQVERYBAD;
1559  return(retval) ;
1560  }
1561  else if(retvsq==3){
1562  // if( vsq > VSQ_TOO_BIG ) {
1563 
1564  // check on v^2>VSQ_TOO_BIG failure (not really a failure of convergence, but of confidence in the solution)
1565  check_utsq_fail(Wp, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra);
1566 
1567  *lpflag= UTOPRIMFAILCONVUTSQ;
1568  return(retval) ;
1569  }
1570 
1571  gtmp = sqrt(1. - vsq);
1572  gamma = 1./gtmp ;
1573  gammasq=gamma*gamma;
1574  // can compute utsq or gamma from Wp alone instead of using vsq
1575  utsq = gammasq*vsq;
1576 
1577  // GODMARK: can compute u+p directly from utsq, Wp, and D.
1578  wmrho0=wmrho0_compute_vsq(Wp,D,vsq,gamma,gammasq);
1579 
1580 
1581 #elif(PRIMFROMVSQ==0)
1582 
1583  utsq=utsq_calc(W, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra);
1584 
1585  gammasq=1.0+utsq;
1586  gamma=sqrt(gammasq);
1587 
1588  wmrho0=wmrho0_compute_utsq(Wp,D,utsq,gamma,gammasq);
1589 
1590 
1591 #endif
1592 
1593  rho0 = D / gamma;
1594 
1595  // w = W * (1. - vsq) ;
1596  // not used, but compute correctly anyways
1597  // w = \rho_0 + u + p
1598  w = rho0+wmrho0;
1599 
1600  p = pressure_wmrho0(whicheos,EOSextra,rho0, wmrho0) ;
1601  // u = w - (rho0 + p) ;
1602  u = wmrho0-p;
1603 
1604 
1605  *pressure=p; // return pressure since may be non-trivial to obtain from p(rho0,u) or inconsistent with p(rho0,\chi) if computed otherwise
1606 
1607  // if((nstep==4)&&(ptrgeom->i==0)&&(ptrgeom->j==47)){
1608  // dualfprintf(fail_file,"Wp=%g W=%g D=%g vsq=%g gamma=%g wmrho0=%g p=%g u=%g\n",Wp,W,D,vsq,gamma,wmrho0,p,u);
1609  // exit(0);
1610  // }
1611 
1612  // dualfprintf(fail_file,"Wp=%g W=%g D=%g vsq=%g gamma=%g wmrho0=%g p=%g u=%g\n",Wp,W,D,vsq,gamma,wmrho0,p,u);
1613 
1614  // dualfprintf(fail_file,"Wp=%g W=%g D=%g utsq=%g gamma=%g wmrho0=%g p=%g u=%g\n",Wp,W,D,utsq,gamma,wmrho0,p,u);
1615 
1616 
1617  // GODMARK
1618  // fix checks for cold case
1619  if( (rho0 <= 0.)
1620  || ((p < 0.)&&(eomtype!=EOMCOLDGRMHD))
1621  ) {
1622  if(showmessages && debugfail>=2 ) {
1623  tmpdiff = w - rho0;
1624  dualfprintf(fail_file,
1625  "rho or uu < 0 failure: rho,w,(w-rho),p,u = %21.15g %21.15g %21.15g %21.15g %21.15g gamma,utsq = %21.15g %21.15g : nstep=%ld steppart=%d :: t=%21.15g :: i=%d j=%d k=%d p=%d\n",rho0,w,tmpdiff,p,u,gamma, utsq,nstep,steppart,t,ptrgeom->i,ptrgeom->j,ptrgeom->k,ptrgeom->p);
1626  }
1627  // Below changed to pressure since internal energy can be negative due to arbitrary energy per baryon offset
1628  if((rho0<=0.)&&(p>=0.)) *lpflag= UTOPRIMFAILRHONEG;
1629  if((rho0>0.)&&(p<0.)) *lpflag= UTOPRIMFAILUNEG;
1630  if((rho0<=0.)&&(p<0.)) *lpflag= UTOPRIMFAILRHOUNEG;
1631  if(UTOPRIMFAILRETURNTYPE==UTOPRIMRETURNNOTADJUSTED) return(retval) ; // else let assign -- used to check how bad failure is.
1632  }
1633 
1634 
1635 
1637  //
1638  // If didn't exit so far, then assume valid inversion and so assign primitives (i.e. prim[])
1639  // Assign:
1640  // 1) rho_0
1641  // 2) u (either from energy or entropy evolution)
1642  // 3) \tilde{u}^i[W']
1643  // 4) [Already inverted field]
1644  //
1646 
1647  prim[RHO] = rho0 ;
1648 
1649 
1650 
1651  if(eomtype==EOMCOLDGRMHD){
1652  prim[UU] = 0;
1653  }
1654  else{
1655  prim[UU] = u ;
1656  }
1657 
1658 
1659 #if(DOENTROPY!=DONOENTROPY)
1660  // put result in both UU and ENTROPY for primitive
1661  prim[ENTROPY]=prim[UU];
1662 #endif
1663 
1664 
1665 
1666  // for(i=1;i<4;i++) Qtcon[i] = Qcon[i] + ncon[i] * Qdotn; // Qtcon already defined in good way
1667 
1668  // GODMARK: gamma should be computed from W and conserved quantities directly as \gamma[W] (see notes)
1669  for(i=1;i<4;i++) prim[UTCON1+i-1] = gamma/(W+Bsq) * ( Qtcon[i] + QdotB*Bcon[i]/W ) ;
1670 
1671 
1672  // dualfprintf(fail_file,"%d %d : %g %g %g\n",ptrgeom->i,ptrgeom->j,vsq,W,Wp);
1673  //if(fabs(prim[U1])>1E-30)
1674  // for(i=1;i<4;i++) dualfprintf(fail_file,"u[%d]= %g %g\n",i,Qtcon[i],prim[UTCON1+i-1]);
1675  // dualfprintf(fail_file,"retval=%d\n",retval);
1676 
1677 
1678 
1679 
1681  //
1682  // make sure nan and inf's don't pass through
1683  //
1685  foundnan=0;
1686  PLOOPALLINVERT(pl){
1687  if(!isfinite(prim[pl])){
1688  foundnan=1;
1689  }
1690  }
1691 
1692  if(foundnan){
1693  *lpflag= UTOPRIMFAILNANRESULT;
1694 
1695  if(debugfail>=2) PLOOPALLINVERT(pl2) dualfprintf(fail_file,"prim is nan: prim[%d]=%21.15g -> pr0[%d]=%21.15g\n",pl2,prim[pl2],pl2,pr0[pl2]);
1696 
1697  // remove nan since might contaminate solution or cause super-slowdown
1698  PLOOPALLINVERT(pl) prim[pl]=pr0[pl];
1699 
1700  return(1);
1701  }
1702 
1703 
1704 #if(!OPTIMIZED)
1705  if( ltrace ) {
1706  dualfprintf(fail_file," rho final = %21.15g , u final = %21.15g \n", rho0, u);
1707  }
1708 #endif
1709 
1710 
1711 
1712  return(0);
1713 
1714 
1715 }
1716 
1717 
1718 
1719 
1721 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)
1722 {
1723  // force-free variables
1724  FTYPE Qtconclean[NDIM],Qtcovclean[NDIM],Qtsqclean;
1725  FTYPE alphasq,gammamax,vcon[NDIM],vcov[NDIM],QtdotB,vconff[NDIM],vcovff[NDIM];
1726  int i;
1727  FTYPE vsq,gamma;
1728  FTYPE realBsq;
1729 
1730  // note that doing this means B.B!=Bsq
1731  Bsq+=SMALL; // in case Bsq=0 identically
1732 
1733 
1734  alphasq=GAMMAMAX*GAMMAMAX/(GAMMAMAX*GAMMAMAX-1.0);
1735 
1736 
1737  prim[RHO]=prim[UU]=0.0;
1738 
1739 
1740  // first 3-velocity
1741  vcon[TT]=0.0; // Qtcon[TT]=0
1742  for(i=1;i<4;i++) vcon[i] = Qtcon[i]/Bsq ;
1743 
1744  // first project out parallel velocity from Qtcon
1745  // note that \tilde{Q}.B = Q.B
1746  QtdotB=0.0;
1747  for(i=1;i<4;i++) QtdotB += Qtcon[i]*Bcov[i];
1748 
1749  Qtconclean[TT]=0.0;
1750  if(cleaned==1){
1751  for(i=1;i<4;i++) Qtconclean[i] = Qtcon[i] - QtdotB*Bcon[i]/Bsq;
1752 
1753  // define new cleaned 3-velocity
1754  vconff[TT]=0.0;
1755  for(i=1;i<4;i++) vconff[i] = Qtconclean[i]/Bsq;
1756  }
1757  else{
1758  for(i=1;i<4;i++) Qtconclean[i] = Qtcon[i];
1759 
1760  // define new non-cleaned 3-velocity
1761  vconff[TT]=0.0;
1762  for(i=1;i<4;i++) vconff[i] = (Qtconclean[i] + 0.0*QtdotB*Bcon[i]/fabs(W) )/(fabs(W) + fabs(Bsq));
1763 
1764  }
1765 
1766 
1767  // v_\alpha cleaned
1768  lower_vec(vconff,ptrgeom,vcovff) ;
1769 
1770  // v^2 cleaned
1771  vsq=0.0;
1772  for(i=1;i<4;i++) vsq+=vconff[i]*vcovff[i]; // vcon[TT]=0
1773 
1774  // dualfprintf(fail_file,"vsq=%21.15g Bsq=%21.15g alphasq=%21.15g Qtsq=%21.15g\n",vsq,Bsq,alphasq,Qtsq);
1775 
1776  // now limit the 3-velocity to be physical (no change in direction, so QtdotB=0 will be preserved)
1777  if( (vsq<0.0) || (vsq>1.0/alphasq) ){
1778  // then need to limit Lorentz factor
1779  // dissipation is expected to occur without affecting force-free equations
1780 
1781  // define cleaned Qtsq
1782  lower_vec(Qtconclean,ptrgeom,Qtcovclean) ;
1783  Qtsqclean=0.0;
1784  for(i=1;i<4;i++) Qtsqclean+=Qtconclean[i]*Qtcovclean[i];
1785 
1786  // define new vconff (overwrites above) that is cleaned AND has dissipation
1787  // v^\alpha = P^\alpha/(B^2+\sqrt{\alpha^2 P^2})
1788  vconff[TT]=0.0;
1789  // for(i=1;i<4;i++) vconff[i] = Qtconclean[i]/(Bsq+sqrt(alphasq*Qtsqclean)) ;
1790  //realBsq=sqrt(Bsq+sqrt(alphasq*Qtsqclean));
1791  // Note that E^2 = P^2/B^2, and then use same limiting procedure as in phys.ffde.c
1792  realBsq=sqrt(Qtsqclean*alphasq);
1793 
1794  if(cleaned==1){
1795  for(i=1;i<4;i++) vconff[i] = Qtconclean[i]/realBsq ;
1796  }
1797  else{
1798  for(i=1;i<4;i++) vconff[i] = (Qtconclean[i] + 0.0*QtdotB*Bcon[i]/fabs(W) )/(fabs(W) + fabs(realBsq));
1799  }
1800 
1801  // v_\alpha cleaned and w/ dissipation
1802  lower_vec(vconff,ptrgeom,vcovff) ;
1803 
1804  // v^2 cleaned and w/ dissipation
1805  vsq=0.0;
1806  for(i=1;i<4;i++) vsq+=vconff[i]*vcovff[i]; // vcon[TT]=0
1807 
1808  }
1809 
1810  // notice that vsq<1 enforced now, so should be ok
1811  if( (vsq<-VSQ_TOO_BIG) || (vsq>1.0/alphasq + 100.0*NUMEPSILON) ){
1812  dualfprintf(fail_file,"Failed to limit vcon\n");
1813  if(vsq<-VSQ_TOO_BIG) vsq=-VSQ_TOO_BIG;
1814  else return(UTOPRIMFAILCONVRET+1);
1815  }
1816 
1817  // gamma cleaned and w/ dissipation if necessary
1818  // clearly limits accuracy of gamma -- unavoidable when using B^\alpha instead of b^\alpha as fundamental variable
1819  gamma=1.0/sqrt(1.0-vsq);
1820 
1821  // set relative 4-velocity
1822  for(i=1;i<4;i++) prim[UTCON1+i-1] = gamma*vconff[i];
1823 
1824  // set field (not really used)
1825  for(i = BCON1; i <= BCON3; i++) prim[i] = U[i] ;
1826 
1827 
1828  // DEBUG:
1829  // PLOOPALLINVERT(i) dualfprintf(fail_file,"vsq=%21.15g prim[%d]=%21.15g\n",vsq,i,prim[i]);
1830 
1831  *gammaret=gamma;
1832 
1833  return(0);
1834 
1835 }
1836 
1837 
1842 static int vsqgtr1(FTYPE W,FTYPE Bsq,FTYPE QdotBsq, FTYPE Qtsq)
1843 {
1844 
1845  return(W*W*W * ( W + 2.*Bsq ) - QdotBsq*(2.*W + Bsq) <= W*W*(Qtsq-Bsq*Bsq) );
1846 
1847 }
1848 
1849 
1850 
1851 
1852 
1853 
1855 //
1856 // GENERAL FUNCTIONS USED BY INVERSION
1857 //
1859 
1860 
1863 static FTYPE utsq_calc(FTYPE W, FTYPE *wglobal,FTYPE Bsq,FTYPE QdotB,FTYPE QdotBsq,FTYPE Qtsq,FTYPE Qdotn,FTYPE Qdotnp,FTYPE D,FTYPE Sc, int whicheos, FTYPE *EOSextra)
1864 {
1865  FTYPE utsqtop,utsqbottom,utsq;
1866  FTYPE S,SoW;
1867 
1868  S = QdotB;
1869  SoW = S/W;
1870 
1871  utsqtop = Qtsq + (Bsq+2.0*W)*SoW*SoW;
1872 
1873  utsqbottom = (Bsq*Bsq-Qtsq) + 2.0*Bsq*W + W*W - SoW*SoW*(Bsq+2.0*W);
1874 
1875  utsq = utsqtop/utsqbottom;
1876 
1877 #if(DEBUGINDEX)
1878  // if(ifileglobal>90 && CRAZYNEWCHECK){
1879  // if(CRAZYNEWCHECK&&0){
1880  // dualfprintf(fail_file,"utsqtop=%21.15g utsqbottom=%21.15g Qtsq=%21.15g\n",utsqtop,utsqbottom,Qtsq);
1881  // }
1882 #endif
1883 
1884  return(utsq);
1885 
1886 }
1887 
1890 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)
1891 {
1892  FTYPE gammatop,gammabottom,gammasq;
1893  FTYPE S,SoW;
1894 
1895  S = QdotB;
1896  SoW = S/W;
1897 
1898  gammatop = Bsq+W;
1899 
1900  gammabottom = W*W+2.0*Bsq*W + (Bsq*Bsq-Qtsq)-2.0*SoW*SoW*W - Bsq*SoW*SoW;
1901 
1902  gammasq = gammatop*gammatop/gammabottom;
1903 
1904  return(gammasq);
1905 
1906 
1907 }
1908 
1911 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)
1912 {
1913  FTYPE gammasq,gamma;
1914 
1915  gammasq = gammasq_calc_W(W, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra);
1916 
1917  // should I check sign of gammasq? GODMARK
1918 
1919  gamma=sqrt(gammasq);
1920 
1921  return(gamma);
1922 
1923 }
1924 
1925 
1926 
1927 
1931 static FTYPE vsq_calc(FTYPE W, FTYPE *wglobal,FTYPE Bsq,FTYPE QdotB,FTYPE QdotBsq,FTYPE Qtsq,FTYPE Qdotn,FTYPE Qdotnp,FTYPE D,FTYPE Sc, int whicheos, FTYPE *EOSextra)
1932 {
1933  FTYPE Wsq,Xsq,Ssq;
1934  FTYPE SoW;
1935 
1936  Wsq = W*W ;
1937  Xsq = (Bsq + W) * (Bsq + W);
1938  // Ssq = QdotBsq / Bsq;
1939 
1940  SoW = QdotB/W;
1941 
1942  //return( Ssq * ( 1./Wsq - 1./Xsq ) + Qtsq / Xsq );
1943  // return( ( Wsq * Qtsq + QdotBsq * (Bsq + 2.*W)) / (Wsq*Xsq) );
1944 
1945  return( (SoW*SoW*(2.0*W+Bsq) + Qtsq)/Xsq );
1946 
1947 }
1948 
1952 static FTYPE dvsq_dW(FTYPE W, FTYPE *wglobal,FTYPE Bsq,FTYPE QdotB,FTYPE QdotBsq,FTYPE Qtsq,FTYPE Qdotn,FTYPE Qdotnp,FTYPE D,FTYPE Sc, int whicheos, FTYPE *EOSextra)
1953 {
1954  FTYPE W3,X3,Ssq,Wsq,X;
1955 
1956  X = Bsq + W;
1957  Wsq = W*W;
1958  W3 = Wsq*W ;
1959  X3 = X*X*X;
1960 
1961  // if(fabs(Bsq)==0.0) Ssq=0.0;
1962  // else Ssq = QdotBsq / Bsq;
1963 
1964  //return( -2.*( Ssq * ( 1./W3 - 1./X3 ) + Qtsq / X3 ) );
1965  // return( -2.*( W3*Qtsq + QdotBsq * ( 3*W*X + Bsq*Bsq ) ) / ( W3 * X3 ) );
1966 
1967  return( -2.0/X3 * ( Qtsq + QdotBsq * (3.0*W*X + Bsq*Bsq)/W3 ) ); // RIGHT (avoids catastrophic cancellation with W^3 term in numerator)
1968 
1969  // return( -2.*( Qtsq/X3 + QdotBsq/Bsq * (1.0/W3 - 1.0/X3) ) ); // RIGHT (said was WRONG!)
1970 
1971 
1972 }
1973 
1974 
1977 static int coldgrhd(PFTYPE *lpflag, FTYPE Qtsq, FTYPE D, FTYPE *Wp)
1978 {
1979  FTYPE det;
1980 
1981  det = D*D+Qtsq;
1982 
1983  if(det<0.0){
1984  // dualfprintf(fail_file,"coldgrhd: no solution: det=%21.15g\n",det);
1985  *lpflag=UTOPRIMFAILCONVUTSQ;
1986  return(1);
1987  }
1988  else{
1989  if(D>0.0){
1990  *Wp = Qtsq/(D+sqrt(det));
1991  }
1992  else{
1993  *Wp = Qtsq/(D-sqrt(det));
1994  }
1995  return(0);
1996  }
1997 
1998  return(0); // shouldn't get here -- only used if coldgrhd() used
1999 }
2000 
2001 
2005 static FTYPE Psq_Wp_old(FTYPE Wp, FTYPE *wglobal,FTYPE Bsq,FTYPE QdotB,FTYPE QdotBsq,FTYPE Qtsq,FTYPE Qdotn,FTYPE Qdotnp,FTYPE D,FTYPE Sc, int whicheos, FTYPE *EOSextra)
2006 {
2007  FTYPE result;
2008  FTYPE S,W,W2,DoW,SoW;
2009  FTYPE Y;
2010 
2011  W = Wp+D;
2012  Y = Bsq+2.0*W;
2013  W2 = W*W;
2014  S = QdotB;
2015  DoW = D/W;
2016  SoW = S/W;
2017 
2018 #if(CRAZYDEBUG&&DEBUGINDEX)
2019  if(ifileglobal==0 && jfileglobal==63 && nstep==9 && steppart==2){
2020  dualfprintf(fail_file,"Qtsq=%21.15g Bsq=%21.15g D=%21.15g S=%21.15g\n",Qtsq,Bsq,D,S);
2021  }
2022 #endif
2023 
2024  // result = -Qtsq + (Bsq+2.0*D)*Bsq + 2.0*(Bsq+D)*Wp + Wp*Wp - (Bsq*DoW*DoW + SoW*SoW)*(Bsq+2.0*W);
2025 
2026  result = -Qtsq + Wp*(D+W)*(1.0 + Bsq*Y/W2) - Y*SoW*SoW;
2027 
2028  return(result);
2029 
2030 }
2031 
2032 
2036 static FTYPE dPsqdWp_Wp_old(FTYPE Wp, FTYPE *wglobal,FTYPE Bsq,FTYPE QdotB,FTYPE QdotBsq,FTYPE Qtsq,FTYPE Qdotn,FTYPE Qdotnp,FTYPE D,FTYPE Sc, int whicheos, FTYPE *EOSextra)
2037 {
2038  FTYPE result;
2039  FTYPE S,W,DoW,SoW;
2040 
2041  W = Wp+D;
2042  S = QdotB;
2043  DoW = D/W;
2044  SoW = S/W;
2045 
2046  result = 2.0*(Bsq+W)*(Bsq/W*DoW*DoW + SoW*SoW/W + 1.0);
2047 
2048  return(result);
2049 
2050 }
2051 
2055 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)
2056 {
2057  FTYPE S,W,W2,DoW,SoW;
2058  FTYPE Y;
2059 
2060  W = Wp+D;
2061  Y = Bsq+2.0*W;
2062  W2 = W*W;
2063  S = QdotB;
2064 
2065  *resid = -Qtsq*W2 + Wp*(D+W)*(W2 + Bsq*Y) - Y*S*S;
2066  *norm = fabs(Qtsq*W2) + fabs(Wp*(D+W)*(W2)) + fabs(Wp*(D+W)*(Bsq*Y)) + fabs(Y*S*S);
2067 
2068  return(0);
2069 
2070 }
2071 
2072 
2076 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)
2077 {
2078  FTYPE result;
2079  FTYPE S,W,DoW,SoW;
2080 
2081  W = Wp+D;
2082  S = QdotB;
2083 
2084  result = 2.0*(D*( (Bsq+D)*(Bsq+D) - Qtsq) - (S*S) + ((Bsq + D)*(Bsq + 5.0*D) - Qtsq)*Wp + 3.0*(Bsq + 2.0*D)*(Wp*Wp) + 2.0*(Wp*Wp*Wp));
2085 
2086  return(result);
2087 
2088 }
2089 
2090 
2091 
2092 
2093 
2094 
2095 
2096 
2097 
2099 #define MINRESIDEPRIME(Wp) ( Bsq*0.5+Qdotnp+3.0/5.0*Wp+(Bsq*Qtsq-QdotBsq)/(2.0*(Bsq+D+Wp)*(Bsq+D+Wp)) )
2100 
2103 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)
2104 {
2105  FTYPE X,X2,W;
2106  FTYPE utsq;
2107  FTYPE p_tmp;
2108 
2109 
2110  W = Wp+D;
2111  X = Bsq + W + SMALL;
2112  X2 = X*X;
2113 
2114 
2115  utsq=utsq_calc(W, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra); // check for precision problems
2116 
2117 
2118  // PERFORMANCEMARK: Tested code with and without these type of if/else (all of them in utoprim_jon.c) and all rest used by func_Eprime_unopt() and found only changed performance from 7% of CPU to 6.90% CPU. Not interesting change
2119 
2120  if(utsq<UTSQNEGLIMIT || utsq!=utsq){ // also checks if nan
2121  *resid=-VERYBIG; // force residual to be large and negative so indicates a problem
2122  *norm=-VERYBIG;
2123  }
2124  else{
2125  p_tmp=pressure_Wp_utsq(whicheos,EOSextra,Wp,D,utsq);
2126 
2127  // Note that -E'=Qdotnp
2128  *resid = Qdotnp + Wp - p_tmp + 0.5*Bsq + 0.5*(Bsq*Qtsq-QdotBsq)/X2;
2129  *norm = fabs(Qdotnp) + fabs(Wp) + fabs(p_tmp) + fabs(0.5*Bsq) + fabs(0.5*(Bsq*Qtsq)/X2) + fabs(0.5*(QdotBsq)/X2);
2130 
2131  if(p_tmp!=p_tmp || *resid!=*resid) *resid=-VERYBIG; // indicates still problem
2132  }
2133 
2134  return(0);
2135 
2136 }
2137 
2138 
2139 
2140 
2145 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)
2146 {
2147  FTYPE result;
2148  FTYPE X,X3,W;
2149  FTYPE utsq;
2150  FTYPE dvsq,dp1,dp2,dpdWp;
2151 
2152 
2153  W = Wp+D;
2154  X = Bsq + W + SMALL;
2155  X3 = X*X*X;
2156 
2157  utsq=utsq_calc(W, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra); // check for precision problems
2158  if(utsq<UTSQNEGLIMIT || utsq!=utsq){ // also checks for nan
2159  result=0.0; // indicates a problem
2160  }
2161  else{
2162  dvsq = dvsq_dW(W, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra); // no precision problem
2163  dp1 = dpdWp_calc_utsq(whicheos,EOSextra,Wp,D,utsq);
2164  dp2 = dpdvsq_calc2_Wp(whicheos,EOSextra, Wp, D, utsq );
2165  dpdWp = dp1 + dp2*dvsq; // dp/dW = dp/dWp
2166 
2167  result = (1.0 - dpdWp) - (Bsq*Qtsq - QdotBsq)/X3;
2168 
2169  // check for nan (but not inf)
2170  if(dvsq!=dvsq || dp1!=dp1 || dp2!=dp2 || dpdWp!=dpdWp || result!=result) result=0.0; // indicates still problem
2171 
2172  }
2173 
2174  return(result);
2175 
2176 }
2177 
2178 
2182 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)
2183 {
2184  FTYPE result;
2185  FTYPE X,X2,W;
2186  FTYPE utsq;
2187  FTYPE p_tmp;
2188 
2189 
2190  W = Wp+D;
2191  X = Bsq + W + SMALL;
2192  X2 = X*X;
2193 
2194  utsq=utsq_calc(W, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra); // check for precision problems
2195  if(utsq<UTSQNEGLIMIT || utsq!=utsq){
2196  result=0.0; // indicates problem
2197  }
2198  else{
2199  p_tmp=pressure_Wp_utsq(whicheos,EOSextra,Wp,D,utsq);
2200 
2201  result = X2*(Qdotnp + Wp - p_tmp + 0.5*Bsq) + 0.5*(Bsq*Qtsq-QdotBsq);
2202 
2203  if(p_tmp!=p_tmp || result!=result) result=0.0; // indicates still problem
2204  }
2205 
2206 #if(DEBUGINDEX)
2207  // if(ifileglobal>90 && CRAZYNEWCHECK){
2208  // if(CRAZYNEWCHECK&&0){
2209  // dualfprintf(fail_file,"Eprime: Wp=%21.15g W=%21.15g X2=%21.15g p_tmp=%21.15g utsq=%21.15g D=%21.15g result=%21.15g\n",Wp,W,X2,p_tmp,utsq,D,result);
2210  // }
2211 #endif
2212 
2213  return(result);
2214 
2215 }
2216 
2217 
2221 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)
2222 {
2223  FTYPE result;
2224  FTYPE X,X3,W;
2225  FTYPE utsq;
2226  FTYPE dvsq,dp1,dp2,dpdWp;
2227  FTYPE p_tmp;
2228 
2229 
2230  W = Wp+D;
2231  X = Bsq + W + SMALL;
2232  X3 = X*X*X;
2233 
2234  utsq=utsq_calc(W, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra); // check for precision problems
2235  if(utsq<UTSQNEGLIMIT || utsq!=utsq){
2236  result=0.0; // indicates problem
2237  }
2238  else{
2239 
2240  dvsq = dvsq_dW(W, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra); // no precision problem
2241  dp1 = dpdWp_calc_utsq(whicheos,EOSextra,Wp,D,utsq);
2242  dp2 = dpdvsq_calc2_Wp(whicheos,EOSextra, Wp, D, utsq );
2243  dpdWp = dp1 + dp2*dvsq; // dp/dW = dp/dWp
2244 
2245  // just copied from Eprime_Wp() to avoid recomputing utsq_calc()
2246  p_tmp=pressure_Wp_utsq(whicheos,EOSextra,Wp,D,utsq);
2247 
2248  result = X*(2.0*(Qdotnp + Wp - p_tmp + 0.5*Bsq) + (1.0 - dpdWp)*X );
2249 
2250  if(dvsq!=dvsq || dp1!=dp1 || dp2!=dp2 || dpdWp!=dpdWp || p_tmp!=p_tmp || result!=result) result=0.0; // indicates still problem
2251 
2252 
2253  }
2254 
2255 #if(DEBUGINDEX)
2256  // if(ifileglobal>90 && CRAZYNEWCHECK){
2257  // if(CRAZYNEWCHECK&&0){
2258  // dualfprintf(fail_file,"dEprime: Wp=%21.15g W=%21.15g X3=%21.15g p_tmp=%21.15g utsq=%21.15g dp1=%21.15g dp2=%21.15g dpdWp=%21.15g result=%21.15g\n",Wp,W,X3,p_tmp,utsq,dp1,dp2,dpdWp,result);
2259  // }
2260 #endif
2261 
2262  return(result);
2263 
2264 }
2265 
2266 
2267 
2268 
2269 
2270 
2271 
2272 
2273 
2275 static FTYPE Sc_Wp_old(FTYPE Wp, FTYPE *wglobal,FTYPE Bsq,FTYPE QdotB,FTYPE QdotBsq,FTYPE Qtsq,FTYPE Qdotn,FTYPE Qdotnp,FTYPE D,FTYPE Sc, int whicheos, FTYPE *EOSextra)
2276 {
2277  FTYPE result;
2278  FTYPE W;
2279  FTYPE utsq;
2280  FTYPE Ss_tmp;
2281 
2282 
2283  W = Wp+D;
2284  utsq=utsq_calc(W, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra); // check for precision problems
2285 
2286 
2287  if(utsq<UTSQNEGLIMIT || utsq!=utsq){ // also checks if nan
2288  result=-VERYBIG; // force residual to be large and negative so indicates a problem
2289  }
2290  else{
2291  Ss_tmp=Ss_Wp_utsq(whicheos,EOSextra,Wp,D,utsq);
2292 
2293  result = -Sc + D*Ss_tmp;
2294 
2295  if(Ss_tmp!=Ss_tmp || result!=result) result=-VERYBIG; // indicates still problem
2296  }
2297 
2298  // DEBUG:
2299  // dualfprintf(fail_file,"Sc_Wp_old: Wp=%21.15g W=%21.15g D=%21.15g utsq=%21.15g :: Ss_tmp=%21.15g result=%21.15g Sc=%21.15g\n",Wp,W,D,utsq,Ss_tmp,result,Sc);
2300 
2301  return(result);
2302 
2303 }
2304 
2307 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)
2308 {
2309  FTYPE W;
2310  FTYPE utsq;
2311  FTYPE Ss_tmp;
2312 
2313 
2314  W = Wp+D;
2315  utsq=utsq_calc(W, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra); // check for precision problems
2316 
2317 
2318  if(utsq<UTSQNEGLIMIT || utsq!=utsq){ // also checks if nan
2319  *resid=-VERYBIG; // force residual to be large and negative so indicates a problem
2320  *norm=-VERYBIG; // force residual to be large and negative so indicates a problem
2321  }
2322  else{
2323  Ss_tmp=Ss_Wp_utsq(whicheos,EOSextra,Wp,D,utsq);
2324 
2325  *resid = Wp*(-Sc + D*Ss_tmp);
2326  *norm = fabs(Wp)*fabs(Sc) + fabs(Wp)*fabs(D*Ss_tmp);
2327 
2328  if(Ss_tmp!=Ss_tmp || *resid!=*resid) *resid=-VERYBIG; // indicates still problem
2329  }
2330 
2331  // DEBUG:
2332  // dualfprintf(fail_file,"Sc_Wp_unopt: Wp=%21.15g W=%21.15g D=%21.15g utsq=%21.15g :: Ss_tmp=%21.15g *resid=%21.15g Sc=%21.15g\n",Wp,W,D,utsq,Ss_tmp,*resid,Sc);
2333 
2334  return(0);
2335 
2336 }
2337 
2338 
2339 
2340 
2342 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)
2343 {
2344  FTYPE result;
2345  FTYPE W;
2346  FTYPE utsq;
2347  FTYPE dvsq,dSs1,dSs2,dSsdWp;
2348  FTYPE Ss;
2349 
2350  W = Wp+D;
2351  utsq=utsq_calc(W, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra); // check for precision problems
2352  if(utsq<UTSQNEGLIMIT || utsq!=utsq){ // also checks for nan
2353  result=0.0; // indicates a problem
2354  }
2355  else{
2356  dvsq = dvsq_dW(W, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra); // no precision problem
2357  dSs1 = dSsdWp_calc_utsq(whicheos,EOSextra,Wp,D,utsq);
2358  dSs2 = dSsdvsq_calc2_Wp(whicheos,EOSextra, Wp, D, utsq );
2359  dSsdWp = dSs1 + dSs2*dvsq; // dSs/dW = dSs/dWp
2360  Ss=Ss_Wp_utsq(whicheos,EOSextra,Wp,D,utsq);
2361 
2362  // old result for old residual
2363  // result = D*dSsdWp;
2364 
2365  // new result for new residual
2366  result = -Sc + D*(dSsdWp*Wp + Ss);
2367 
2368  // check for nan (but not inf)
2369  if(dvsq!=dvsq || dSs1!=dSs1 || dSs2!=dSs2 || dSsdWp!=dSsdWp || result!=result) result=0.0; // indicates still problem
2370 
2371  }
2372 
2373  // DEBUG:
2374  // dualfprintf(fail_file,"dScdWp: Wp=%21.15g W=%21.15g D=%21.15g utsq=%21.15g :: dvsq=%21.15g dSs1=%21.15g dSs2=%21.15g dSsdWp=%21.15g result=%21.15g\n",Wp,W,D,utsq,dvsq,dSs1,dSs2,dSsdWp,result);
2375 
2376 
2377  return(result);
2378 
2379 }
2380 
2381 
2382 
2383 
2384 
2385 
2386 
2387 
2388 
2389 
2390 
2391 
2392 
2393 /********************************************************************
2394 
2395  x1_of_x0(): Only Used For 2D Inversion
2396 
2397  -- calculates x1 from x0 depending on value of vartype;
2398 
2399  x0 x1
2400  vartype = 2 : Wp, gamma
2401  vartype = 3 : Wp, vsq
2402  vartype = 4 : Wp, utsq
2403 
2404  -- asumes x0 is already physical
2405 
2406 *********************************************************************/
2407 
2408 static int x1_of_x0(FTYPE x0, FTYPE *x1, FTYPE *wglobal,FTYPE Bsq,FTYPE QdotB,FTYPE QdotBsq,FTYPE Qtsq,FTYPE Qdotn,FTYPE Qdotnp,FTYPE D,FTYPE Sc, int whicheos, FTYPE *EOSextra)
2409 {
2410  FTYPE vsq;
2411  // FTYPE dv = 1.e-15;
2412  FTYPE W;
2413  int retval;
2414 
2415 
2416 
2417  W = x0+D;
2418 
2419 #if(CHANGEDTOOLDER)
2420  vsq = fabs(vsq_calc(W, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra)) ; // guaranteed to be positive (>=0)
2421 #else
2422  vsq = vsq_calc(W, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra) ; // ok if not positive
2423 #endif
2424 
2425  retval=validate_vsq(vsq,&vsq, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra);
2426 
2427  // return( ( vsq > 1. ) ? VSQ_TOO_BIG : vsq );
2428 
2429  *x1=vsq;
2430  return(retval); // indicates failure of v^2 or not
2431 
2432 }
2433 
2434 /********************************************************************
2435 
2436  validate_x():
2437 
2438  -- makes sure that x[0,1] have physical values, based upon their definitions:
2439 
2440  // added wglobal to properly normalize -- global variable set before find_root_2d()
2441 
2442  // avoided truncating bottom value of x[0] and x[1] since then can't converge properly to slightly unphysical values or temporarily reach out into unphysical space and come back in.
2443 
2444  // avoid truncating x[1] for v^2<0 since all equations are ok for v^2<0 and can help convergence to nearly physical or even physical solutions
2445 
2446  // must avoid x[0] too large and v^2>1, however, since that produces nan's.
2447 
2448  *********************************************************************/
2449 
2450 
2453 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)
2454 {
2455  FTYPE Wp;
2456  FTYPE dv = NUMEPSILON;
2457  PFTYPE vlpflag;
2458 
2459  /* Always take the absolute value of Wpnew[0] and check to see if it's too big: */
2460  // Wpnew[0] = fabs(Wpold);
2461  // below allows cases when can be negative internal energy density. Primarily allowed by wglobal[0].
2462  FTYPE Wpnewsmall=-3.0*MAX(MAX(MAX(MAX(MAX(fabs(D),fabs(Qdotn)),fabs(Bsq)),fabs(Qdotnp)),fabs(wglobal[0])),fabs(wglobal[1]));
2463  // dualfprintf(fail_file,"Wpnewsmall=%Lg : %Lg %Lg %Lg %Lg %Lg %Lg\n",Wpnewsmall,D,Qtsq,Qdotn,*wglobal,Bsq,Qdotnp);
2464  if(Wpnew[0]<Wpnewsmall){
2465  // if(Wpnew[0]<-D && D>0.0){
2466  // Wpnew[0] = 10.0*Qtsq/(2.0*D+Wp0);
2467  // Wpnew[0]=10.0*D+dv;
2468 
2469  // reset above HD solution
2470  coldgrhd(&vlpflag, Qtsq, D, &Wp); // ignore vlpflag result
2471  Wpnew[0] = Wp*10.0;
2472 
2473 #if(CRAZYDEBUG)
2474  dualfprintf(fail_file,"Wpnew[0] changed from %21.15g to %21.15g (D=%21.15g Qtsq=%21.15g)\n",Wpold,Wpnew[0],D,Qtsq);
2475 #endif
2476 
2477  }
2478 
2479  // revert to previous Wp if "new" Wp very bad
2480  // can't go back to previous Wp otherwise won't go anywhere
2481  // Wpnew[0]/global = gamma^2
2482  Wpnew[0] = (fabs(Wpnew[0]/wglobal[1]) > GAMMASQ_TOO_BIG) ? 0.5*(Wpold+Wpnew[0]) : Wpnew[0];
2483 
2484 
2485 }
2486 
2487 
2492 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)
2493 {
2494  int retval;
2495 
2496 
2497  retval=0; // assume no problem
2498 
2499  if(vsqnew[0]>1.0){
2500  retval=1;
2501  vsqnew[0]=VSQ_TOO_BIG;
2502  return(retval);
2503  }
2504 
2505 
2506  if(vsqnew[0]>VSQ_TOO_BIG){
2507  retval=3;
2508  vsqnew[0]=VSQ_TOO_BIG;
2509  return(retval);
2510  }
2511 
2512 #if(VSQNEGCHECK==1)
2513  // unnecessary limitation on Newton's method
2514  if( (vsqnew[0] < 0.) && (fabs(vsqnew[0]) < MAXNEGVSQ) ) {
2515  vsqnew[0] = 0.0;
2516  retval=2;
2517  return(retval);
2518  }
2519 #elif(VSQNEGCHECK==2)
2520  // best can do, but need to indicate that did it since Newton's method will fail to compute residual and derivative correctly
2521  if(vsqnew[0] <=-VSQ_TOO_BIG) {
2522  vsqnew[0] = -VSQ_TOO_BIG;
2523  retval=2;
2524  return(retval);
2525  }
2526 #elif(VSQNEGCHECK==3)
2527  if(vsqnew[0]<0.0) {
2528  vsqnew[0]=0.0;
2529  retval=2;
2530  return(retval);
2531  }
2532 #endif
2533 
2534 
2535  return(retval);
2536 
2537 }
2538 
2539 
2541 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)
2542 {
2543 
2544  validate_Wp(x0[0], &x[0], wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra);
2545 
2546  validate_vsq(x0[1], &x[1], wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra);
2547 
2548 }
2549 
2550 
2552 static void validate_x_1d_old(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)
2553 {
2554  FTYPE dv = NUMEPSILON;
2555 
2556  /* Always take the absolute value of x[0] and check to see if it's too big: */
2557  // x[0] = fabs(x[0]);
2558  FTYPE Wpnewsmall=-3.0*MAX(MAX(MAX(MAX(MAX(fabs(D),fabs(Qdotn)),fabs(Bsq)),fabs(Qdotnp)),fabs(wglobal[0])),fabs(wglobal[1]));
2559  // dualfprintf(fail_file,"Wpnewsmall=%Lg : %Lg %Lg %Lg %Lg %Lg %Lg\n",Wpnewsmall,D,Qtsq,Qdotn,*wglobal,Bsq,Qdotnp);
2560  if(x[0]<Wpnewsmall) x[0]=D+dv;
2561  // x[0]/global = gamma^2
2562  x[0] = (x[0]/wglobal[1] > GAMMASQ_TOO_BIG) ? x0[0] : x[0];
2563 
2564 }
2565 
2566 
2567 
2569 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)
2570 {
2571 
2572  validate_Wp(x0[0], &x[0], wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra);
2573 
2574 }
2575 
2576 
2578 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))
2579 {
2580 
2581  if(
2582  (WHICHHOTINVERTER==1 && (eomtype==EOMGRMHD) ) ||
2583  (WHICHHOTINVERTER==3 && (eomtype==EOMGRMHD) ) ||
2584  ( (WHICHCOLDINVERTER==0 || WHICHCOLDINVERTER==1) && eomtype==EOMCOLDGRMHD) ||
2585  (WHICHCOLDINVERTER==2 && eomtype==EOMCOLDGRMHD)
2586  || (eomtype==EOMENTROPYGRMHD)
2587  ){
2588  *ptr_validate_x = &validate_x_1d;
2589  }
2590  else if(WHICHHOTINVERTER==2 && (eomtype==EOMGRMHD)){
2591  *ptr_validate_x = &validate_x_2d;
2592  }
2593  else {
2594  dualfprintf(fail_file,"pick_validate_x: Shouldn't reach here\n");
2595  myexit(999);
2596  }
2597 
2598 
2599 }
2600 
2601 static void validate_x(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),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)
2602 {
2603  (*ptr_validate_x)(x,x0,wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra);
2604 
2605 }
2606 
2607 
2608 
2610 //
2611 //
2613 //
2614 //
2616 
2617 
2618 
2619 
2620 /****************************************************
2621 *****************************************************/
2622 
2623 static void func_1d_orig(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)
2624 {
2625  FTYPE vsq,vsqfake,gamma,gamma_sq,rho0,w,W,Wp,drdW,dpdW,Wsq,p_tmp, dvsq,dp1,dp2,a_term,ap_term ;
2626  FTYPE X,X2,X3;
2627  // FTYPE dv = 1.0e-10;
2628  const int ltrace = 0;
2629  int retvalvalidatevsq;
2630 
2631  Wp = x[0];
2632  W = Wp+D;
2633  Wsq = W*W;
2634  X = Bsq+W;
2635  X2 = X*X;
2636  X3 = X2*X;
2637 
2638  //dualfprintf(fail_file,"call utsq in residual\n") ;
2639  vsq = vsq_calc(W, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra) ;
2640  //dualfprintf(fail_file,"done w/ utsq in residual\n") ;
2641 
2642  // GODMARK DANGEROUS
2643  // vsqfake = ( vsq < -1.e-10 ) ? dv : fabs(vsq) ;
2644  // vsqfake = ( vsq > 1. ) ? (1.-dv) : vsq ;
2645 
2646  retvalvalidatevsq=validate_vsq(vsq,&vsqfake, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra);
2647 
2648 #if(0)
2649  // no longer use below methods
2650 #if(OPTIMIZED)
2651  vsqfake = ( vsq < -MAXNEGVSQ ) ? 0.0 : fabs(vsq) ;
2652  vsqfake = ( vsq > 1. ) ? VSQ_TOO_BIG : vsq ;
2653 #elif(0)
2654  if(vsq<0.0){
2655  if(vsq>-MAXNEGVSQ) vsqfake=0.0;
2656  else{
2657  dualfprintf(fail_file,"func_1d_orig: vsq=%21.15g\n",vsq);
2658  myexit(1);
2659  }
2660  }
2661  if(vsq>VSQ_TOO_BIG){
2662  if(vsq<=1.0) vsqfake=VSQ_TOO_BIG;
2663  else{
2664  dualfprintf(fail_file,"func_1d_orig: vsq=%21.15g\n",vsq);
2665  myexit(1);
2666  }
2667  }
2668 #endif
2669 #endif
2670 
2671  if(retvalvalidatevsq==0){
2672 
2673  //rho0 = D * sqrt(1. - vsq) ;
2674  //dualfprintf(fail_file,"rho0,D,gamma: %21.15g %21.15g %21.15g\n",rho0,D,gamma) ;
2675  // w = W * (1. - vsq);
2676 
2677  p_tmp = pressure_Wp_vsq(whicheos,EOSextra,Wp,D,vsqfake); // needs to be correct
2678 
2679 
2680 
2681  // Same as previous residual, however, jacobian is calculated using full differentiation w.r.t. W
2682  //
2683  // dp/dW' = dp/dW + dP/dv^2 dv^2/dW
2684 
2685 
2686  dvsq = dvsq_dW(W, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra); // no precision problem
2687  dp1 = dpdWp_calc_vsq(whicheos,EOSextra, Wp, D, vsq ); // vsq can be unphysical
2688  dp2 = dpdvsq_calc_Wp(whicheos,EOSextra, Wp, D, vsq ); // vsq must be physical
2689  dpdW = dp1 + dp2*dvsq; // dp/dW = dp/dWp
2690 
2691 
2692 #if(0) // Scott
2693 
2694  resid[0] =
2695  + Wp
2696  + 0.5 * Bsq * ( 1. + vsq ) // any vsq
2697  - 0.5*QdotBsq/Wsq
2698  + Qdotnp
2699  - p_tmp;
2700 
2701  norm[0] =
2702  + fabs(Wp)
2703  + fabs(0.5 * Bsq * ( 1. + vsq ))
2704  + fabs(0.5*QdotBsq/Wsq)
2705  + fabs(Qdotnp)
2706  + fabs(p_tmp);
2707 
2708  jac[0][0] = drdW = 1. - dpdW + QdotBsq/(Wsq*W) + 0.5*Bsq*dvsq;
2709 
2710 #elif(0)
2711  // resid[0] =
2712  // + W*Wsq
2713  // + 0.5 * Wsq * Bsq * ( 1. + vsq )
2714  // - 0.5*QdotBsq
2715  // - Qdotn*Wsq
2716  // - p_tmp*Wsq;
2717  //
2718  // jac[0][0] = drdW = W * ( 3.*W + Bsq * (1. + vsq + 0.5*W*dvsq) - 2.*Qdotn - 2.*p_tmp - W*dpdW );
2719 
2720 #elif(1) // JON:
2721 
2722  resid[0] = Qdotnp + Wp - p_tmp + 0.5*Bsq + (Bsq*Qtsq - QdotBsq)/X2;
2723  norm[0] = fabs(Qdotnp) + fabs(Wp) + fabs(p_tmp) + fabs(0.5*Bsq) + fabs((Bsq*Qtsq)/X2) + fabs((QdotBsq)/X2);
2724 
2725  jac[0][0] = drdW = 1. - dpdW + (Bsq*Qtsq - QdotBsq)/X3*(-2.0);
2726 
2727 #endif
2728 
2729 
2730  dx[0] = -resid[0]/drdW;
2731 
2732  *f = 0.5*resid[0]*resid[0];
2733  *df = -2. * (*f);
2734 
2735 #if(!OPTIMIZED)
2736  if( ltrace ) {
2737  a_term = 0.5*Bsq*(1.+vsq) + Qdotnp - p_tmp;
2738  ap_term = dvsq * ( 0.5*Bsq - dp2 ) - dp1;
2739 
2740  dualfprintf(fail_file,"func_1d_orig(): x = %21.15g, dx = %21.15g, resid = %21.15g, drdW = %21.15g \n", x[0], dx[0], resid[0], drdW);
2741  dualfprintf(fail_file,"func_1d_orig(): dvsq = %21.15g, dp = %21.15g , dp1 = %21.15g , dp2 = %21.15g \n", dvsq, dpdW,dp1,dp2);
2742  dualfprintf(fail_file,"func_1d_orig(): W = %21.15g , vsq = %21.15g , a = %21.15g, a_prime = %21.15g \n", x[0], vsq, a_term, ap_term);
2743 
2744  }
2745 #endif
2746  }
2747  else{ // indicate failure for this 1-D method since depends on v^2 or \tilde{u}^2 to be well-defined as a function of W'
2748  resid[0]=-VERYBIG;
2749  norm[0]=-VERYBIG;
2750  jac[0][0]=0.0;
2751  dx[0]=0.0;
2752  *f=-VERYBIG;
2753  *df=-VERYBIG;
2754  }
2755 
2756 }
2757 
2758 
2759 
2760 
2761 
2762 /****************************************************
2763  Routine for line searching for the 1D method .
2764 *****************************************************/
2765 
2766 static FTYPE res_sq_1d_orig(FTYPE x[], FTYPE *wglobal,FTYPE Bsq,FTYPE QdotB,FTYPE QdotBsq,FTYPE Qtsq,FTYPE Qdotn,FTYPE Qdotnp,FTYPE D,FTYPE Sc, int whicheos, FTYPE *EOSextra)
2767 {
2768  FTYPE vsq,W,Wp,Wsq,p_tmp,resid[1],norm[1];
2769  FTYPE vsqfake;
2770  FTYPE X,X2;
2771  // FTYPE dv = 1.0e-10;
2772  const int ltrace = 0;
2773  int retvalvalidatevsq;
2774 
2775 
2776  Wp = x[0];
2777  W = Wp+D;
2778  Wsq = W*W;
2779  X = Bsq+W;
2780  X2 = X*X;
2781 
2782  vsq = vsq_calc(W, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra) ;
2783 
2784 
2785  retvalvalidatevsq=validate_vsq(vsq,&vsqfake, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra);
2786 
2787  if(retvalvalidatevsq==0){
2788 
2789 #if(0)
2790  // GODMARK DANGEROUS
2791 #if(OPTIMIZED)
2792  // scn's version
2793  // vsqfake = ( vsq < -dv ) ? dv : fabs(vsq) ;
2794  //vsqfake = ( vsq > 1. ) ? (1.-dv) : vsq ;
2795  // jon's version
2796  vsqfake = ( vsq < 0.0 ) ? 0.0 : vsq ;
2797  vsqfake = ( vsq > VSQ_TOO_BIG ) ? VSQ_TOO_BIG : vsq ; // GODMARK: must let iterate if vsq>1. Can't just force vsq or else Newton's method won't work.
2798 #else
2799  if(vsq<0.0){
2800  if(vsq>-MAXNEGVSQ) vsqfake=0.0;
2801  else{
2802  dualfprintf(fail_file,"func_1d_orig: vsq=%21.15g\n",vsq);
2803  myexit(1);
2804  }
2805  }
2806  if(vsq>VSQ_TOO_BIG){
2807  if(vsq<=1.0) vsqfake=VSQ_TOO_BIG;
2808  else{
2809  dualfprintf(fail_file,"func_1d_orig: vsq=%21.15g\n",vsq);
2810  myexit(1);
2811  }
2812  }
2813 #endif
2814 #endif
2815 
2816  p_tmp = pressure_Wp_vsq(whicheos,EOSextra,Wp,D,vsqfake); // vsq<=1 for pressure
2817 
2818 
2819 
2820 #if(0) // Scott
2821  resid[0] =
2822  + Wp
2823  + 0.5 * Bsq * ( 1. + vsq ) // residual must use unphysical vsq
2824  - 0.5*QdotBsq/Wsq
2825  + Qdotnp
2826  - p_tmp;
2827 
2828  norm[0] =
2829  + fabs(Wp)
2830  + fabs(0.5 * Bsq * ( 1. + vsq ))
2831  + fabs(0.5*QdotBsq/Wsq)
2832  + fabs(Qdotnp)
2833  + fabs(p_tmp);
2834 
2835 #elif(1) // JON:
2836 
2837  resid[0] = Qdotnp + Wp - p_tmp + 0.5*Bsq + (Bsq*Qtsq - QdotBsq)/X2;
2838  norm[0] = fabs(Qdotnp) + fabs(Wp) + fabs(p_tmp) + fabs(0.5*Bsq) + fabs((Bsq*Qtsq)/X2) + fabs((QdotBsq)/X2);
2839 
2840 #endif
2841  }
2842  else{
2843  resid[0]=-VERYBIG;
2844  norm[0]=-VERYBIG;
2845  }
2846 
2847 
2848  return( 0.5*resid[0]*resid[0] );
2849 
2850 
2851 }
2852 
2853 
2854 
2855 
2856 
2857 
2859 //
2860 // SCN versions of 1D method
2861 //
2863 
2864 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)
2865 {
2866  FTYPE vsq,vsqfake,gamma,gamma_sq,rho0,w,W,drdW,dpdW,Wsq,p_tmp, dvsq,dp1,dp2,a_term,ap_term ;
2867  // FTYPE dv = 1.0e-10;
2868  const int ltrace = 0;
2869 
2870 
2871  W = x[0];
2872  Wsq = W*W;
2873 
2874  //dualfprintf(fail_file,"call utsq in residual\n") ;
2875  vsq = vsq_calc(W, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra) ;
2876  //dualfprintf(fail_file,"done w/ utsq in residual\n") ;
2877 
2878  // GODMARK DANGEROUS
2879  // vsq = ( vsq < -1.e-10 ) ? dv : fabs(vsq) ;
2880  // vsq = ( vsq > 1. ) ? (1.-dv) : vsq ;
2881 
2882 #if(OPTIMIZED)
2883  vsqfake = ( vsq < -MAXNEGVSQ ) ? 0.0 : fabs(vsq) ;
2884  vsqfake = ( vsq > 1. ) ? VSQ_TOO_BIG : vsq ;
2885 #else
2886  if(vsq<0.0){
2887  if(vsq>-MAXNEGVSQ) vsqfake=0.0;
2888  else{
2889  dualfprintf(fail_file,"func_1d_orig: vsq=%21.15g\n",vsq);
2890  myexit(1);
2891  }
2892  }
2893  if(vsq>VSQ_TOO_BIG){
2894  if(vsq<=1.0) vsqfake=VSQ_TOO_BIG;
2895  else{
2896  dualfprintf(fail_file,"func_1d_orig: vsq=%21.15g\n",vsq);
2897  myexit(1);
2898  }
2899  }
2900 #endif
2901 
2902  //rho0 = D * sqrt(1. - vsq) ;
2903  //dualfprintf(fail_file,"rho0,D,gamma: %21.15g %21.15g %21.15g\n",rho0,D,gamma) ;
2904  // w = W * (1. - vsq);
2905 
2906  p_tmp = pressure_W_vsq_scn(whicheos,EOSextra, W,D,vsqfake);
2907 
2908 
2909 
2910  // Same as previous residual, however, jacobian is calculated using full differentiation w.r.t. W
2911 
2912  dvsq = dvsq_dW(W, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra);
2913  dp1 = dpdW_calc_vsq_scn(whicheos,EOSextra, W, D, vsq ); // use real vsq
2914  dp2 = dpdvsq_calc_scn(whicheos,EOSextra, W, D, vsqfake );
2915  dpdW = dp1 + dp2*dvsq;
2916 
2917  resid[0] =
2918  + W
2919  + 0.5 * Bsq * ( 1. + vsq ) // use real vsq
2920  - 0.5*QdotBsq/Wsq
2921  + Qdotn
2922  - p_tmp;
2923 
2924  norm[0] =
2925  + fabs(W )
2926  + fabs(0.5 * Bsq * ( 1. + vsq ))
2927  + fabs(0.5*QdotBsq/Wsq)
2928  + fabs(Qdotn)
2929  + fabs(p_tmp);
2930 
2931  jac[0][0] = drdW = 1. - dpdW + QdotBsq/(Wsq*W) + 0.5*Bsq*dvsq;
2932 
2933  // resid[0] =
2934  // + W*Wsq
2935  // + 0.5 * Wsq * Bsq * ( 1. + vsq )
2936  // - 0.5*QdotBsq
2937  // - Qdotn*Wsq
2938  // - p_tmp*Wsq;
2939  //
2940  // jac[0][0] = drdW = W * ( 3.*W + Bsq * (1. + vsq + 0.5*W*dvsq) - 2.*Qdotn - 2.*p_tmp - W*dpdW );
2941 
2942  dx[0] = -resid[0]/drdW;
2943 
2944  *f = 0.5*resid[0]*resid[0];
2945  *df = -2. * (*f);
2946 
2947 #if(!OPTIMIZED)
2948  if( ltrace ) {
2949  a_term = 0.5*Bsq*(1.+vsq) + Qdotn - p_tmp;
2950  ap_term = dvsq * ( 0.5*Bsq - dp2 ) - dp1;
2951 
2952  dualfprintf(fail_file,"func_1d_orig(): x = %21.15g, dx = %21.15g, resid = %21.15g, drdW = %21.15g \n", x[0], dx[0], resid[0], drdW);
2953  dualfprintf(fail_file,"func_1d_orig(): dvsq = %21.15g, dp = %21.15g , dp1 = %21.15g , dp2 = %21.15g \n", dvsq, dpdW,dp1,dp2);
2954  dualfprintf(fail_file,"func_1d_orig(): W = %21.15g , vsq = %21.15g , a = %21.15g, a_prime = %21.15g \n", x[0], vsq, a_term, ap_term);
2955 
2956  }
2957 #endif
2958 
2959 }
2960 
2961 /****************************************************
2962  Routine for line searching for the 1D method .
2963 *****************************************************/
2964 
2965 static FTYPE res_sq_1d_orig_scn(FTYPE x[], FTYPE *wglobal,FTYPE Bsq,FTYPE QdotB,FTYPE QdotBsq,FTYPE Qtsq,FTYPE Qdotn,FTYPE Qdotnp,FTYPE D,FTYPE Sc, int whicheos, FTYPE *EOSextra)
2966 {
2967  FTYPE vsq,vsqfake,W,Wsq,p_tmp,resid[1],norm[1];
2968  // FTYPE dv = 1.0e-10;
2969  const int ltrace = 0;
2970 
2971 
2972  W = x[0];
2973  Wsq = W*W;
2974 
2975  vsq = vsq_calc(W, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra) ;
2976 
2977  // GODMARK DANGEROUS
2978 #if(OPTIMIZED)
2979  // scn's version
2980  // vsq = ( vsq < -dv ) ? dv : fabs(vsq) ;
2981  //vsq = ( vsq > 1. ) ? (1.-dv) : vsq ;
2982  // jon's version
2983  vsqfake = ( vsq < 0.0 ) ? 0.0 : vsq ;
2984  vsqfake = ( vsq > VSQ_TOO_BIG ) ? VSQ_TOO_BIG : vsq ;
2985 #else
2986  if(vsq<0.0){
2987  if(vsq>-MAXNEGVSQ) vsqfake=0.0;
2988  else{
2989  dualfprintf(fail_file,"func_1d_orig: vsq=%21.15g\n",vsq);
2990  myexit(1);
2991  }
2992  }
2993  if(vsq>VSQ_TOO_BIG){
2994  if(vsq<=1.0) vsqfake=VSQ_TOO_BIG;
2995  else{
2996  dualfprintf(fail_file,"func_1d_orig: vsq=%21.15g\n",vsq);
2997  myexit(1);
2998  }
2999  }
3000 #endif
3001 
3002  p_tmp = pressure_W_vsq_scn(whicheos,EOSextra,W,D,vsqfake);
3003 
3004 
3005  resid[0] =
3006  + W
3007  + 0.5 * Bsq * ( 1. + vsq )
3008  - 0.5*QdotBsq/Wsq
3009  + Qdotn
3010  - p_tmp;
3011 
3012  norm[0] =
3013  + fabs(W)
3014  + fabs(0.5 * Bsq * ( 1. + vsq ))
3015  + fabs(0.5*QdotBsq/Wsq)
3016  + fabs(Qdotn)
3017  + fabs(p_tmp);
3018 
3019  return( 0.5*resid[0]*resid[0] );
3020 
3021 
3022 }
3023 
3024 
3025 
3026 
3027 
3028 
3029 
3031 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)
3032 {
3033  FTYPE drdW;
3034  // FTYPE dv = 1.0e-10;
3035  const int ltrace = 0;
3036  FTYPE Wp;
3037 
3038 
3039  Wp = x[0];
3040 
3041  Eprime_Wp_unopt(Wp, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra, &resid[0], &norm[0]);
3042  jac[0][0] = drdW = dEprimedWp_unopt(Wp, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra);
3043 
3044 
3045  dx[0] = -resid[0]/drdW;
3046  // if((fabs(dx[0])/(Wp+D)>0.1){ // don't trust such a big jump
3047  // dx[0] = sign(dx[0])*0.1;
3048  // }
3049  *f = 0.5*resid[0]*resid[0];
3050  *df = -2. * (*f);
3051 
3052 
3053 }
3054 
3055 
3057 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)
3058 {
3059  FTYPE Wp;
3060  FTYPE resid,norm;
3061 
3062  Wp = x[0];
3063  Eprime_Wp_unopt(Wp, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra, &resid, &norm);
3064 
3065  return( 0.5*resid*resid );
3066 
3067 
3068 }
3069 
3070 
3071 
3072 
3075 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)
3076 {
3077 
3079  //
3080  // compute common quantities
3081  //
3083 
3084 
3085  *W = Wp+D;
3086  *X = Bsq + *W + SMALL;
3087  *X2 = (*X) *(*X);
3088  *X3 = (*X2) *(*X);
3089 
3090  *utsq=utsq_calc(*W, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra); // check for precision problems
3091 
3092  // get \gamma^2, \gamma, \chi, and \rho_0 (no EOS used):
3093  *gammasq=1.0+ (*utsq);
3094  *gamma=sqrt(*gammasq);
3095  *wmrho0=wmrho0_compute_utsq(Wp, D, *utsq, *gamma, *gammasq);
3096  *rho0=D/(*gamma);
3097 
3098 
3099 
3100 }
3101 
3104 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)
3105 {
3106 
3107  // total derivative dv^2/dW'
3108  *dvsq = dvsq_dW(W, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra); // no precision problem -- and no EOS calls inside -- purely kinetic
3109 
3110  *dwmrho0dW = 1.0/gammasq; // holding utsq fixed
3111  *drho0dW = 0.0; // because \rho=D/\gamma and holding utsq fixed
3112 
3113  *dwmrho0dvsq = (D*(gamma*0.5-1.0) - Wp); // holding Wp fixed
3114  *drho0dvsq = -D*gamma*0.5; // because \rho=D/\gamma and holding Wp fixed
3115 
3116 
3117 }
3118 
3119 
3120 
3121 
3122 
3125 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)
3126 {
3127  // FTYPE dv = 1.0e-10;
3128  const int ltrace = 0;
3129 
3130 
3132  //
3133  // Get some kinetic quantities
3134  //
3136  FTYPE Wp;
3137  FTYPE X,X2,X3,W;
3138  FTYPE utsq;
3139  FTYPE gamma,gammasq;
3140  FTYPE rho0,wmrho0;
3141  Wp = x[0];
3142  get_kinetics_part1(Wp, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc, whicheos, EOSextra, &W, &X, &X2, &X3, &utsq, &gamma, &gammasq, &rho0, &wmrho0);
3143 
3144 
3146  //
3147  // Now can get thermodynamical quantities and results from them
3148  //
3150  // PERFORMANCEMARK: Tested code with and without these type of if/else (all of them in utoprim_jon.c) and all rest used by func_Eprime() and found only changed performance from 7% of CPU to 6.90% CPU. Not interesting change
3151  FTYPE Eprime,dEprimedWp;
3152  if(gammasq<0.0 || utsq<UTSQNEGLIMIT || utsq!=utsq){ // also checks if nan
3153  Eprime=-VERYBIG; // force residual to be large and negative so indicates a problem
3154  dEprimedWp=0.0; // indicates a problem
3155 
3156  resid[0]=Eprime;
3157  norm[0]=Eprime;
3158  jac[0][0]=dEprimedWp;
3159  dx[0] = VERYBIG;
3160  *f = 0.5*(Eprime*Eprime);
3161  *df = -2. * (*f);
3162 
3163  }
3164  else{
3166  //
3167  // Now can get pressure and other thermodynamical quantities
3168  // Optimize for any EOS table lookup by computing P(\rho_0,\chi), dP/d\rho_0 and dP/d\chi now.
3169  // Otherwise, Eprime_Wp() and dEprimeWp() (as previously coded) repeated these calls and made calls for each EOS term while now can make a pipelined call and obtain all 3 EOS values at once.
3170  //
3172  FTYPE pofchi,dpdrho,dpdchi;
3173  getall_forinversion(whicheos,EOMGRMHD,CHIDIFF,EOSextra,rho0,wmrho0,&pofchi,&dpdrho,&dpdchi);
3174 
3176  //
3177  // get Eprime
3178  // Note that -E'=Qdotnp
3179  //
3181  Eprime = Qdotnp + Wp - pofchi + 0.5*Bsq + 0.5*(Bsq*Qtsq-QdotBsq)/X2;
3182  FTYPE normEprime = fabs(Qdotnp) + fabs(Wp) + fabs(pofchi) + fabs(0.5*Bsq) + fabs(0.5*(Bsq*Qtsq)/X2) + fabs(0.5*(QdotBsq)/X2);
3183 
3184 
3185  if(pofchi!=pofchi || Eprime!=Eprime) Eprime=-VERYBIG; // indicates still problem
3186 
3188  //
3189  // Get derivative-related kinetics
3190  //
3192  FTYPE dvsq,dwmrho0dW,drho0dW;
3193  FTYPE dwmrho0dvsq,drho0dvsq;
3194  get_kinetics_part2(Wp, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc, whicheos, EOSextra, W, X, X2, X3, utsq, gamma, gammasq, rho0, wmrho0, &dvsq,&dwmrho0dW,&drho0dW,&dwmrho0dvsq,&drho0dvsq);
3195 
3197  //
3198  // get mixed kinetic-thermal derivatives
3199  //
3201  FTYPE dpdW,dpdvsq,dpdWp;
3202  // dvsq = dv^2/dW' total derivative
3203  dpdW = drho0dW *dpdrho + dwmrho0dW *dpdchi; // dP/dW' holding utsq fixed
3204  dpdvsq = drho0dvsq *dpdrho + dwmrho0dvsq *dpdchi; // dP/dvsq holding W' fixed where utsq and vsq are same fixed quantity
3205  dpdWp = dpdW + dpdvsq*dvsq; // dp/dW = dp/dWp [total derivative]
3206 
3208  //
3209  // get dEprimedWp
3210  //
3212  dEprimedWp = (1.0 - dpdWp) - (Bsq*Qtsq - QdotBsq)/X3;
3213 
3214 
3215  // check for nan (but not inf)
3216  if(dvsq!=dvsq || dpdW!=dpdW || dpdvsq!=dpdvsq || dpdWp!=dpdWp || dEprimedWp!=dEprimedWp) dEprimedWp=0.0; // indicates still problem
3217 
3218 
3219 
3220 
3221 
3222 
3224  //
3225  // residual, jacobian, and step
3226  //
3228  resid[0]=Eprime;
3229  norm[0]=normEprime;
3230  jac[0][0]=dEprimedWp;
3231  dx[0] = -Eprime/dEprimedWp;
3232  *f = 0.5*(Eprime*Eprime);
3233  *df = -2. * (*f);
3234 
3235 
3236  // DEBUG:
3237  // if(x[0]+dx[0]<SMALL){
3238  // dualfprintf(fail_file,"EPRIMEPROBLEM: %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g :: %21.15g %21.15g %21.15g :: %21.15g %21.15g %21.15g %21.15g %21.15g :: %21.15g %21.15g %21.15g %21.15g %21.15g :: %21.15g %21.15g\n"
3239  // ,Wp, Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc, W, X, X2, X3, utsq, gamma, gammasq, rho0, wmrho0
3240  // ,pofchi,dpdrho,dpdchi
3241  // ,dvsq,dwmrho0dW,drho0dW,dwmrho0dvsq,drho0dvsq
3242  // ,dpdW,dpdvsq,dpdWp,Eprime,dEprimedWp
3243  // ,x[0],dx[0]);
3244  // int extraloop;
3245  // for(extraloop=FIRSTEOSGLOBAL;extraloop<=LASTEOSGLOBAL;extraloop++) dualfprintf(fail_file,"EPRIMEPROBELMEXTRA: %21.15g\n",EOSextra[extraloop]);
3246  // }
3247 
3248 
3249  }
3250 
3251 
3252 
3253 
3254 
3255 }
3256 
3257 
3259 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)
3260 {
3261  FTYPE jac[NEWT_DIM][NEWT_DIM]; // dummy variable, return value not used
3262  FTYPE f,df; // dummy, return value not used
3263  int n=1; // true for this residual
3264  FTYPE resid[NEWT_DIM];
3265  FTYPE norm[NEWT_DIM];
3266 
3267  // just call func_Eprime_opt() since more optimized than before -- could optimize further but don't often or at all use this function
3268  func_Eprime_opt(x, dx, resid, norm, jac, &f, &df, n, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc, whicheos, EOSextra);
3269 
3270  return( 0.5*resid[0]*resid[0] );
3271 
3272 
3273 }
3274 
3275 
3276 
3277 
3279 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)
3280 {
3281  FTYPE drdW,Wp;
3282  // FTYPE dv = 1.0e-10;
3283  const int ltrace = 0;
3284 
3285 
3286  Wp = x[0];
3287 
3288 
3289  Psq_Wp(Wp, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra, &resid[0], &norm[0]);
3290  jac[0][0] = drdW = dPsqdWp_Wp(Wp, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra);
3291 
3292 
3293  dx[0] = -resid[0]/drdW;
3294  *f = 0.5*resid[0]*resid[0];
3295  *df = -2. * (*f);
3296 
3297 
3298 }
3299 
3300 
3301 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)
3302 {
3303  FTYPE Wp;
3304  FTYPE resid,norm;
3305 
3306  Wp = x[0];
3307  Psq_Wp(Wp, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra, &resid, &norm);
3308 
3309  return( 0.5*resid*resid );
3310 
3311 
3312 }
3313 
3314 
3315 
3316 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)
3317 {
3318  FTYPE drdW;
3319  // FTYPE dv = 1.0e-10;
3320  const int ltrace = 0;
3321  FTYPE Wp;
3322 
3323 
3324  Wp = x[0];
3325 
3326  Sc_Wp_unopt(Wp, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra, &resid[0], &norm[0]);
3327  jac[0][0] = drdW = dScdWp_unopt(Wp, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra);
3328 
3329 
3330  dx[0] = -resid[0]/drdW;
3331  // if((fabs(dx[0])/(Wp+D)>0.1){ // don't trust such a big jump
3332  // dx[0] = sign(dx[0])*0.1;
3333  // }
3334  *f = 0.5*resid[0]*resid[0];
3335  *df = -2. * (*f);
3336 
3337 
3338 }
3339 
3340 
3341 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)
3342 {
3343  FTYPE Wp;
3344  FTYPE resid,norm;
3345 
3346  Wp = x[0];
3347  Sc_Wp_unopt(Wp, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra, &resid, &norm);
3348 
3349  return( 0.5*resid*resid );
3350 
3351 
3352 }
3353 
3354 
3357 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)
3358 {
3359 
3360 
3362  //
3363  // Get some kinetic quantities
3364  //
3366  FTYPE Wp;
3367  FTYPE X,X2,X3,W;
3368  FTYPE utsq;
3369  FTYPE gamma,gammasq;
3370  FTYPE rho0,wmrho0;
3371  Wp = x[0];
3372  get_kinetics_part1(Wp, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc, whicheos, EOSextra, &W, &X, &X2, &X3, &utsq, &gamma, &gammasq, &rho0, &wmrho0);
3373 
3374 
3376  //
3377  // Now can get thermodynamical quantities and results from them
3378  //
3380  // PERFORMANCEMARK: Tested code with and without these type of if/else (all of them in utoprim_jon.c) and all rest used by func_Eprime() and found only changed performance from 7% of CPU to 6.90% CPU. Not interesting change
3381  FTYPE Scprime,dScprimedWp;
3382  if(gammasq<0.0 || utsq<UTSQNEGLIMIT || utsq!=utsq){ // also checks if nan
3383  Scprime=-VERYBIG; // force residual to be large and negative so indicates a problem
3384  dScprimedWp=0.0; // indicates a problem
3385 
3386  resid[0]=Scprime;
3387  norm[0]=Scprime;
3388  jac[0][0]=dScprimedWp;
3389  dx[0] = -Scprime/dScprimedWp;
3390  *f = 0.5*(Scprime*Scprime);
3391  *df = -2. * (*f);
3392 
3393  }
3394  else{
3396  //
3397  // Now can get pressure and other thermodynamical quantities
3398  // Optimize for any EOS table lookup by computing P(\rho_0,\chi), dP/d\rho_0 and dP/d\chi now.
3399  // Otherwise, Eprime_Wp() and dEprimeWp() (as previously coded) repeated these calls and made calls for each EOS term while now can make a pipelined call and obtain all 3 EOS values at once.
3400  //
3402  FTYPE Ssofchi,dSsdrho,dSsdchi;
3403  getall_forinversion(whicheos,EOMENTROPYGRMHD,CHIDIFF,EOSextra,rho0,wmrho0,&Ssofchi,&dSsdrho,&dSsdchi);
3404 
3406  //
3407  // get Scprime
3408  //
3410  // old residual
3411  Scprime = -Sc + D*Ssofchi;
3412  FTYPE Scprimenorm = fabs(Sc) + fabs(D*Ssofchi);
3413 
3414  // new residual (problem near Wp=0 of course that could diverge to Wp=0 and be wrong solution. Was using this "new" residual to help reach solution faster, but in general doesn't help so much and hurts since sometimes with tabulated EOS reaches Wp=0 and results in bad inversion.)
3415  // Also, for tabulated EOS, new residual more complicated and error doesn't drop. Even though old residual has many iterations like new residual, at least old residual has very small error. This is probably because derivative is not exactly consistent with function in the tabulated EOS yet both are used in the derivative expression.
3417 
3418 
3419  if(Ssofchi!=Ssofchi || Scprime!=Scprime) Scprime=-VERYBIG; // indicates still problem
3420 
3422  //
3423  // Get derivative-related kinetics
3424  //
3426  FTYPE dvsq,dwmrho0dW,drho0dW;
3427  FTYPE dwmrho0dvsq,drho0dvsq;
3428  get_kinetics_part2(Wp, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc, whicheos, EOSextra, W, X, X2, X3, utsq, gamma, gammasq, rho0, wmrho0, &dvsq,&dwmrho0dW,&drho0dW,&dwmrho0dvsq,&drho0dvsq);
3429 
3431  //
3432  // get mixed kinetic-thermal derivatives
3433  //
3435  FTYPE dSsdW,dSsdvsq,dSsdWp;
3436  dSsdW = drho0dW *dSsdrho + dwmrho0dW *dSsdchi; // dSs/dW' holding utsq fixed
3437  dSsdvsq = drho0dvsq *dSsdrho + dwmrho0dvsq *dSsdchi;
3438  dSsdWp = dSsdW + dSsdvsq*dvsq; // dSs/dW = dSs/dWp [total derivative]
3439 
3441  //
3442  // get dScprimeWp
3443  //
3445  // old result for old residual
3446  dScprimedWp = D*dSsdWp;
3447 
3448  // new result for new residual
3450 
3451  // check for nan (but not inf)
3452  if(dvsq!=dvsq || dSsdW!=dSsdW || dSsdvsq!=dSsdvsq || dSsdWp!=dSsdWp || dScprimedWp!=dScprimedWp) dScprimedWp=0.0; // indicates still problem
3453 
3454 
3456  //
3457  // residual, jacobian, and step
3458  //
3460  resid[0]=Scprime;
3461  norm[0]=Scprimenorm;
3462  jac[0][0]=dScprimedWp;
3463  dx[0] = -Scprime/dScprimedWp;
3464  *f = 0.5*(Scprime*Scprime);
3465  *df = -2. * (*f);
3466 
3467  }
3468 
3469 
3470 
3471 }
3472 
3473 
3476 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)
3477 {
3478  FTYPE jac[NEWT_DIM][NEWT_DIM]; // dummy variable, return value not used
3479  FTYPE f,df; // dummy, return value not used
3480  int n=1; // true for this residual
3481  FTYPE resid[NEWT_DIM];
3482  FTYPE norm[NEWT_DIM];
3483 
3484  // just call func_Eprime_opt() since more optimized than before -- could optimize further but don't often or at all use this function
3485  func_Sc_opt(x, dx, resid, norm, jac, &f, &df, n, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc, whicheos, EOSextra);
3486 
3487  return( 0.5*resid[0]*resid[0] );
3488 
3489 
3490 }
3491 
3492 
3493 
3494 
3496 //
3498 //
3500 
3501 
3502 
3503 
3504 /*********************************************************
3505 returns residual of vsq-1
3506 **********************************************************/
3507 
3508 static void func_vsq(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)
3509 {
3510 
3511 
3512  FTYPE Wp,W, vsq, Wsq, p_tmp, dPdvsq, dPdW, temp, detJ,tmp2,tmp3;
3513  const int ltrace = 0;
3514 
3515 
3516  Wp = x[0];
3517  W = Wp+D;
3518  vsq = x[1];
3519 
3520  Wsq = W*W;
3521 
3522  p_tmp = pressure_Wp_vsq(whicheos,EOSextra, Wp, D, vsq );
3523  dPdW = dpdWp_calc_vsq(whicheos,EOSextra, Wp, D, vsq );
3524  dPdvsq = dpdvsq_calc_Wp(whicheos,EOSextra, Wp, D, vsq );
3525 
3526 
3527 #if(CRAZYDEBUG&&DEBUGINDEX)
3528  if(ifileglobal==0 && jfileglobal==31 && nstep==9 && steppart==2){
3529  dualfprintf(fail_file,"p_tmp=%21.15g dPdW=%21.15g dPdvsq=%21.15g\n",p_tmp,dPdW,dPdvsq);
3530  }
3531 #endif
3532 
3533 
3534  /* These expressions were calculated using Mathematica */
3535  /* Since we know the analytic form of the equations, we can explicitly
3536  calculate the Newton-Raphson step: */
3537 
3538  dx[0] = (-Bsq/2. + dPdvsq)*(Qtsq - vsq*((Bsq+W)*(Bsq+W)) +
3539  (QdotBsq*(Bsq + 2*W))/Wsq) +
3540  ((Bsq+W)*(Bsq+W))*(-Qdotnp - (Bsq*(1 + vsq))/2. + QdotBsq/(2.*Wsq) -Wp + p_tmp);
3541 
3542  dx[1] = -((-1 + dPdW - QdotBsq/(Wsq*W))*
3543  (Qtsq - vsq*((Bsq+W)*(Bsq+W)) + (QdotBsq*(Bsq + 2*W))/Wsq)) -
3544  2*(vsq + QdotBsq/(Wsq*W))*(Bsq + W)*
3545  (-Qdotnp - (Bsq*(1 + vsq))/2. + QdotBsq/(2.*Wsq) - Wp + p_tmp);
3546 
3547  detJ = (Bsq + W)*((-1 + dPdW - QdotBsq/(Wsq*W))*(Bsq + W) +
3548  ((Bsq - 2*dPdvsq)*(QdotBsq + vsq*(Wsq*W)))/(Wsq*W));
3549 
3550  dx[0] /= -(detJ) ;
3551  dx[1] /= -(detJ) ;
3552 
3553  jac[0][0] = -2*(vsq + QdotBsq/(Wsq*W))*(Bsq + W);
3554  jac[0][1] = -(Bsq + W)*(Bsq + W);
3555  jac[1][0] = -1 - QdotBsq/(Wsq*W) + dPdW;
3556  jac[1][1] = -Bsq/2. + dPdvsq;
3557  resid[0] = Qtsq - vsq*(Bsq + W)*(Bsq + W) + (QdotBsq*(Bsq + 2*W))/Wsq;
3558  resid[1] = -Qdotnp - (Bsq*(1 + vsq))/2. + QdotBsq/(2.*Wsq) - Wp + p_tmp;
3559 
3560  norm[0] = fabs(Qtsq) + fabs(vsq*(Bsq + W)*(Bsq + W)) + fabs((QdotBsq*(Bsq + 2*W))/Wsq);
3561  norm[1] = fabs(Qdotnp) + fabs((Bsq*(1 + vsq))/2.) + fabs(QdotBsq/(2.*Wsq)) + fabs(Wp) + fabs(p_tmp);
3562 
3563  *df = -resid[0]*resid[0] - resid[1]*resid[1];
3564 
3565  *f = -0.5 * ( *df );
3566 
3567 
3568 #if(!OPTIMIZED)
3569  if( ltrace ) {
3570  dualfprintf(fail_file,"func_vsq(): x[0] = %21.15g , x[1] = %21.15g , dx[0] = %21.15g , dx[1] = %21.15g , f = %21.15g \n",
3571  x[0], x[1],dx[0], dx[1], *f);
3572 
3573  }
3574 #endif
3575 
3576 }
3577 
3578 
3579 /**********************************************************/
3580 
3581 static void func_vsq2(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)
3582 {
3583 
3584 
3585  FTYPE W, Wp,vsq, Wsq, p_tmp, dPdvsq, dPdW, temp, detJ,tmp2,tmp3;
3586  FTYPE t11;
3587  FTYPE t16;
3588  FTYPE t18;
3589  FTYPE t2;
3590  FTYPE t21;
3591  FTYPE t23;
3592  FTYPE t24;
3593  FTYPE t25;
3594  FTYPE t3;
3595  FTYPE t35;
3596  FTYPE t36;
3597  FTYPE t4;
3598  FTYPE t40;
3599  FTYPE t9;
3600 
3601  const int ltrace = 0;
3602 
3603 
3604  Wp = x[0];
3605  W = Wp+D;
3606  vsq = x[1];
3607 
3608  Wsq = W*W;
3609 
3610  p_tmp = pressure_Wp_vsq(whicheos,EOSextra, Wp, D, vsq );
3611  dPdW = dpdWp_calc_vsq(whicheos,EOSextra, Wp, D, vsq );
3612  dPdvsq = dpdvsq_calc_Wp(whicheos,EOSextra, Wp, D, vsq );
3613 
3614 #if(CRAZYDEBUG&&DEBUGINDEX)
3615  if(ifileglobal==0 && jfileglobal==31 && nstep==9 && steppart==2){
3616  dualfprintf(fail_file,"W=%21.15g D=%21.15g vsq=%21.15g\n",W,D,vsq);
3617  dualfprintf(fail_file,"p_tmp=%21.15g dPdW=%21.15g dPdvsq=%21.15g\n",p_tmp,dPdW,dPdvsq);
3618  }
3619 #endif
3620 
3621 #if(0)
3622  dualfprintf(fail_file,"p_tmp=%21.15g dPdW=%21.15g dPdvsq=%21.15g\n",p_tmp,dPdW,dPdvsq);
3623 #endif
3624 
3625  /* These expressions were calculated using Mathematica, but made into efficient code using Maple */
3626  /* Since we know the analytic form of the equations, we can explicitly
3627  calculate the Newton-Raphson step: */
3628 
3629  t2 = -0.5*Bsq+dPdvsq;
3630  t3 = Bsq+W;
3631  t4 = t3*t3;
3632  t9 = 1/Wsq;
3633  t11 = Qtsq-vsq*t4+QdotBsq*(Bsq+2.0*W)*t9;
3634  FTYPE t11norm = fabs(Qtsq)+fabs(vsq*t4)+fabs(QdotBsq*(Bsq+2.0*W)*t9);
3635  t16 = QdotBsq*t9;
3636  t18 = -Qdotnp-0.5*Bsq*(1.0+vsq)+0.5*t16-Wp+p_tmp;
3637  FTYPE t18norm = fabs(Qdotnp)+fabs(0.5*Bsq*(1.0+vsq))+fabs(0.5*t16)+fabs(Wp)+fabs(p_tmp);
3638  //t18 = -Qdotn-0.5*Bsq*(1.0+vsq)+0.5*t16-W+p_tmp;
3639  t21 = 1/t3;
3640  t23 = 1/W;
3641  t24 = t16*t23;
3642  t25 = -1.0+dPdW-t24;
3643  t35 = t25*t3+(Bsq-2.0*dPdvsq)*(QdotBsq+vsq*Wsq*W)*t9*t23;
3644  t36 = 1/t35;
3645  dx[0] = -(t2*t11+t4*t18)*t21*t36;
3646 
3647 #if(0)
3648  dualfprintf(fail_file,"Qdotn=%21.15g Qdotnp=%21.15g Wp=%21.15g p_tmp=%21.15g\n",Qdotn,Qdotnp,Wp,p_tmp);
3649  dualfprintf(fail_file,"t2=%21.15g t3=%21.15g t4=%21.15g t9=%21.15g t11=%21.15g t16=%21.15g t18=%21.15g t21=%21.15g t23=%21.15g t24=%21.15g t25=%21.15g t35=%21.15g t36=%21.15g dx[0]=%21.15g\n",t2,t3,t4,t9,t11,t16,t18,t21,t23,t24,t25,t35,t36,dx[0]);
3650 #endif
3651 
3652  t40 = (vsq+t24)*t3;
3653  dx[1] = -(-t25*t11-2.0*t40*t18)*t21*t36;
3654  detJ = t3*t35;
3655  jac[0][0] = -2.0*t40;
3656  jac[0][1] = -t4;
3657  jac[1][0] = t25;
3658  jac[1][1] = t2;
3659  resid[0] = t11;
3660  resid[1] = t18;
3661  norm[0] = t11norm;
3662  norm[1] = t18norm;
3663 
3664  *df = -resid[0]*resid[0] - resid[1]*resid[1];
3665 
3666  *f = -0.5 * ( *df );
3667 
3668 #if(CRAZYDEBUG&&DEBUGINDEX)
3669  if(ifileglobal==0 && jfileglobal==31 && nstep==9 && steppart==2){
3670  dualfprintf(fail_file,"Qdotn=%21.15g QdotBsq=%21.15g Qtsq=%21.15g Bsq=%21.15g vsq=%21.15g W=%21.15g p_tmp=%21.15g -Qdotnp-Wp=%21.15g\n",Qdotn,QdotBsq,Qtsq,Bsq,vsq,W,p_tmp,-Qdotnp-Wp);
3671 
3672  dualfprintf(fail_file,"t2=%21.15g\nt3=%21.15g\nt4=%21.15g\nt9=%21.15g\nt11=%21.15g\nt16=%21.15g\nt18=%21.15g\nt21=%21.15g\nt23=%21.15g\nt24=%21.15g\nt25=%21.15g\nt35=%21.15g\nt36=%21.15g\ndx[0]=%21.15g\nt40=%21.15g\ndx[1]=%21.15g\ndetJ=%21.15g\njac[0][0]=%21.15g\njac[0][1]=%21.15g\njac[1][0]=%21.15g\njac[1][1]=%21.15g\nresid[0]=%21.15g\nresid[1]=%21.15g\n",t2,t3,t4,t9,t11,t16,t18,t21,t23,t24,t25,t35,t36,dx[0],t40,dx[1],detJ,jac[0][0],jac[0][1],jac[1][0],jac[1][1],resid[0],resid[1]);
3673 
3674  dualfprintf(fail_file,"*df=%21.15g *f=%21.15g Wp=%21.15g W=%21.15g\n",*df,*f,Wp,W);
3675  }
3676 #endif
3677 
3678 
3679 #if(!OPTIMIZED)
3680  if( ltrace ) {
3681  dualfprintf(fail_file,"func_vsq2(): x[0] = %21.15g , x[1] = %21.15g , dx[0] = %21.15g , dx[1] = %21.15g , f = %21.15g \n",
3682  x[0], x[1],dx[0], dx[1], *f);
3683 
3684  }
3685 #endif
3686 
3687 }
3688 
3689 
3690 /*********************************************************
3691 **********************************************************/
3692 static FTYPE res_sq_vsq(FTYPE x[], FTYPE *wglobal,FTYPE Bsq,FTYPE QdotB,FTYPE QdotBsq,FTYPE Qtsq,FTYPE Qdotn,FTYPE Qdotnp,FTYPE D,FTYPE Sc, int whicheos, FTYPE *EOSextra)
3693 {
3694 
3695 
3696  FTYPE W, Wp, vsq, Wsq, p_tmp, dPdvsq, dPdW, temp, detJ,tmp2,tmp3;
3697  FTYPE resid[2],norm[2];
3698  const int ltrace = 0;
3699 
3700 
3701  Wp = x[0];
3702  W = Wp+D;
3703  vsq = x[1];
3704 
3705  Wsq = W*W;
3706 
3707  p_tmp = pressure_Wp_vsq(whicheos,EOSextra, Wp, D, vsq );
3708  // dPdW = dpdW_calc_vsq(whicheos,EOSextra, Wp, D, vsq );
3709  // dPdvsq = dpdvsq_calc_Wp(whicheos,EOSextra, Wp, D, vsq );
3710 
3711  /* These expressions were calculated using Mathematica */
3712  /* Since we know the analytic form of the equations, we can explicitly
3713  calculate the Newton-Raphson step: */
3714 
3715  resid[0] = Qtsq - vsq*(Bsq + W)*(Bsq + W) + (QdotBsq*(Bsq + 2*W))/Wsq;
3716  resid[1] = -Qdotnp - (Bsq*(1 + vsq))/2. + QdotBsq/(2.*Wsq) - Wp + p_tmp;
3717 
3718  norm[0] = fabs(Qtsq) + fabs(vsq*(Bsq + W)*(Bsq + W)) + fabs((QdotBsq*(Bsq + 2*W))/Wsq);
3719  norm[1] = fabs(Qdotnp) + fabs((Bsq*(1 + vsq))/2.) + fabs(QdotBsq/(2.*Wsq)) + fabs(Wp) + fabs(p_tmp);
3720 
3721 
3722  tmp3 = 0.5 * ( resid[0]*resid[0] + resid[1]*resid[1] ) ;
3723 
3724 #if(!OPTIMIZED)
3725  if( ltrace ) {
3726  dualfprintf(fail_file,"res_sq_vsq(): W,vsq,resid0,resid1,f = %21.15g %21.15g %21.15g %21.15g %21.15g \n",
3727  W,vsq,resid[0],resid[1],tmp3);
3728  dualfprintf(fail_file,"res_sq_vsq(): Qtsq, Bsq, QdotBsq, Qdotn = %21.15g %21.15g %21.15g %21.15g \n",
3729  Qtsq, Bsq, QdotBsq, Qdotn);
3730  }
3731 #endif
3732 
3733  return( tmp3 );
3734 }
3735 
3736 
3737 /*********************************************************
3738 **********************************************************/
3739 static FTYPE res_sq_vsq2(FTYPE x[], FTYPE *wglobal,FTYPE Bsq,FTYPE QdotB,FTYPE QdotBsq,FTYPE Qtsq,FTYPE Qdotn,FTYPE Qdotnp,FTYPE D,FTYPE Sc, int whicheos, FTYPE *EOSextra)
3740 {
3741 
3742  FTYPE W, Wp, vsq, Wsq, p_tmp, dPdvsq, dPdW, temp, detJ,tmp2,tmp3;
3743  FTYPE resid[2],norm[2];
3744  FTYPE t3, t4 ;
3745  FTYPE t9 ;
3746  FTYPE t11;
3747  FTYPE t16;
3748  FTYPE t18;
3749 
3750  const int ltrace = 0;
3751 
3752 
3753  Wp = x[0];
3754  W = Wp+D;
3755  vsq = x[1];
3756 
3757  Wsq = W*W;
3758 
3759  p_tmp = pressure_Wp_vsq(whicheos,EOSextra, Wp, D, vsq );
3760 
3761  /* These expressions were calculated using Mathematica and made into efficient code using Maple*/
3762  /* Since we know the analytic form of the equations, we can explicitly
3763  calculate the Newton-Raphson step: */
3764 
3765  t3 = Bsq+W;
3766  t4 = t3*t3;//
3767  t9 = 1/Wsq; //
3768  t11 = Qtsq-vsq*t4+QdotBsq*(Bsq+2.0*W)*t9; //
3769  FTYPE t11norm = fabs(Qtsq)+fabs(vsq*t4)+fabs(QdotBsq*(Bsq+2.0*W)*t9); //
3770  t16 = QdotBsq*t9;//
3771  t18 = -Qdotnp-0.5*Bsq*(1.0+vsq)+0.5*t16-Wp+p_tmp;//
3772  FTYPE t18norm = fabs(Qdotnp)+fabs(0.5*Bsq*(1.0+vsq))+fabs(0.5*t16)+fabs(Wp)+fabs(p_tmp);//
3773  resid[0] = t11;
3774  resid[1] = t18;
3775  norm[0] = t11norm;
3776  norm[1] = t18norm;
3777 
3778  tmp3 = 0.5 * ( resid[0]*resid[0] + resid[1]*resid[1] ) ;
3779 
3780 #if(!OPTIMIZED)
3781  if( ltrace ) {
3782  dualfprintf(fail_file,"res_sq_vsq(): W,vsq,resid0,resid1,f = %21.15g %21.15g %21.15g %21.15g %21.15g \n",
3783  W,vsq,resid[0],resid[1],tmp3);
3784  dualfprintf(fail_file,"res_sq_vsq(): Qtsq, Bsq, QdotBsq, Qdotn = %21.15g %21.15g %21.15g %21.15g \n",
3785  Qtsq, Bsq, QdotBsq, Qdotn);
3786  }
3787 #endif
3788 
3789 #if(CRAZYDEBUG&&DEBUGINDEX)
3790  if(ifileglobal==0 && jfileglobal==31 && nstep==9 && steppart==2){
3791  dualfprintf(fail_file,"res: W=%21.15g Wsq=%21.15g p_tmp=%21.15g t3=%21.15g t4=%21.15g t9=%21.15g t11=%21.15g t16=%21.15g t18=%21.15g resid[0]=%21.15g resid[1]=%21.15g tmp3=%21.15g\n",W,Wsq,p_tmp,t3,t4,t9,t11,t16,t18,resid[0],resid[1],tmp3);
3792  }
3793 #endif
3794 
3795 
3796  return( tmp3 );
3797 }
3798 
3799 
3800 
3801 
3802 
3803 
3804 
3805 
3806 
3807 
3808 /****************************************************
3809 
3810 
3811 GENERAL PROCEDURES BELOW (except validate_x() and must choose function and residual)
3812 
3813 
3814 *****************************************************/
3815 
3816 
3817 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)
3818 {
3819  FTYPE x_1d[1];
3820  int retval;
3821 
3822 
3823  x_1d[0] = x0;
3824 
3825  if( (retval=general_newton_raphson(showmessages, lpflag, eomtype, x_1d, 1, USE_LINE_SEARCH, func_1d_orig, res_sq_1d_orig, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra,newtonstats) ) ) {
3826 #if(!OPTIMIZED)
3827  if( debugfail>=2 ) dualfprintf(fail_file, "GNR failed, x_1d[0] = %21.15g \n", x_1d[0] );
3828 #endif
3829  }
3830 
3831  *xnew = x_1d[0];
3832 
3833  return(retval);
3834 }
3835 
3836 
3837 
3838 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)
3839 {
3840  FTYPE x_1d[1];
3841  int retval;
3842 
3843 
3844  x_1d[0] = x0;
3845 
3846  if( (retval=general_newton_raphson(showmessages, lpflag, eomtype, x_1d, 1, USE_LINE_SEARCH, func_1d_orig_scn, res_sq_1d_orig_scn, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra,newtonstats) ) ) {
3847 #if(!OPTIMIZED)
3848  if( debugfail>=2 ) dualfprintf(fail_file, "GNR failed, x_1d[0] = %21.15g \n", x_1d[0] );
3849 #endif
3850  }
3851 
3852  *xnew = x_1d[0];
3853 
3854  return(retval);
3855 }
3856 
3857 
3858 
3859 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)
3860 {
3861  FTYPE x_1d[1];
3862  int retval;
3863 
3864 
3865  x_1d[0] = x0;
3866 
3867  // Note u sing func_Eprime_opt() and res_sq_Eprime_opt() : optimized versions.
3868  if( (retval=general_newton_raphson(showmessages, lpflag, eomtype, x_1d, 1, USE_LINE_SEARCH, func_Eprime_opt, res_sq_Eprime_opt, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra,newtonstats) ) ) {
3869 #if(!OPTIMIZED)
3870  if( debugfail>=2 ) dualfprintf(fail_file, "GNR failed, x_1d[0] = %21.15g \n", x_1d[0] );
3871 #endif
3872  }
3873 
3874  *xnew = x_1d[0];
3875 
3876  return(retval);
3877 }
3878 
3879 
3880 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)
3881 {
3882  FTYPE x_1d[1];
3883  int retval;
3884 
3885 
3886  x_1d[0] = x0;
3887 
3888  if( (retval=general_newton_raphson(showmessages, lpflag, eomtype, x_1d, 1, USE_LINE_SEARCH, func_Psq, res_sq_Psq, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra,newtonstats) ) ) {
3889 #if(!OPTIMIZED)
3890  if( debugfail>=2 ) dualfprintf(fail_file, "GNR failed, x_1d[0] = %21.15g \n", x_1d[0] );
3891 #endif
3892  }
3893 
3894  *xnew = x_1d[0];
3895 
3896  return(retval);
3897 }
3898 
3899 
3900 
3901 
3903 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)
3904 {
3905  FTYPE x_3d[3];
3906  int retval;
3907 
3908 
3909  x_3d[0] = x0;
3910  x_3d[1] = x0;
3911  x_3d[2] = x0;
3912 
3913  if( (retval=general_newton_raphson(showmessages, lpflag, eomtype, x_3d, 3, USE_LINE_SEARCH, func_Psq, res_sq_Psq, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra,newtonstats) ) ) {
3914 #if(!OPTIMIZED)
3915  if( debugfail>=2 ) dualfprintf(fail_file, "GNR failed, x_3d[0] = %21.15g x_3d[1] = %21.15g x_3d[2] = %21.15g \n", x_3d[0], x_3d[1], x_3d[2] );
3916 #endif
3917  }
3918 
3919  *xnew = x_3d[0];
3920  *xnew = x_3d[1];
3921  *xnew = x_3d[2];
3922 
3923  return(retval);
3924 }
3925 
3926 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)
3927 {
3928  FTYPE x_1d[1];
3929  int retval;
3930 
3931 
3932  x_1d[0] = x0;
3933 
3934  // Note, using optimized versions of func_Sc() and res_sq_Sc()
3935  if( (retval=general_newton_raphson(showmessages, lpflag, eomtype, x_1d, 1, USE_LINE_SEARCH, func_Sc_opt, res_sq_Sc_opt, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra,newtonstats) ) ) {
3936 #if(!OPTIMIZED)
3937  if( debugfail>=2 ) dualfprintf(fail_file, "GNR failed, x_1d[0] = %21.15g \n", x_1d[0] );
3938 #endif
3939  }
3940 
3941  *xnew = x_1d[0];
3942 
3943  return(retval);
3944 }
3945 
3946 
3947 
3948 
3949 /********************************************************************
3950 
3951  find_root_2D_gen():
3952 
3953  -- performs a 2D Newton-Raphson method for solving the \tilde{Q}^2
3954  and Q.n equations for two variables in general;
3955  (usually W and something else like utsq, vsq, or gamma)
3956 
3957  residual vector = { Qtsq eq. , Q.n eq. }
3958  indep. var. vector, x = { W, vsq/utsq/gamma/? }
3959 
3960  vartype = 2 : W, gamma
3961  vartype = 3 : W, vsq
3962  vartype = 4 : W, utsq
3963 
3964  -- Uses general_newton_raphson();
3965 
3966 *********************************************************************/
3967 
3968 
3969 static int find_root_2D_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)
3970 {
3971 
3972  FTYPE x[2], x_orig[2];
3973  int ntries = 0;
3974  int retval, n = 2;
3975  int it;
3976 
3977  const int ltrace = 0;
3978 
3979  /* Set presets: */
3980  x_orig[0] = x[0] = fabs(x0) ;
3981  x1_of_x0( x0, &x[1], wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra) ;
3982  x_orig[1] = x[1];
3983 
3984 #if(CRAZYDEBUG&&DEBUGINDEX)
3985  if(ifileglobal==0 && jfileglobal==31 && nstep==9 && steppart==2){
3986  dualfprintf(fail_file,"x[0]=%21.15g W=%21.15g x[1]=%21.15g\n",x[0],x[0]+D,x[1]);
3987  }
3988 #endif
3989 
3990 #if(!OPTIMIZED)
3991  if( ltrace ) {
3992  dualfprintf(fail_file, "find_root_2D_gen(): x[0] = %21.15g , x[1] = %21.15g \n", x[0], x[1]);
3993 
3994  }
3995 #endif
3996 
3997 #if(0)
3998  for(it=0;it<n;it++) dualfprintf(fail_file,"x[%d]=%21.15g xnew[%d]=%21.15g\n",it,x[it],it,xnew[it]);
3999 #endif
4000 
4001  retval = general_newton_raphson(showmessages, lpflag, eomtype, x, n, USE_LINE_SEARCH, func_vsq2, res_sq_vsq2, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra,newtonstats) ;
4002  ntries++;
4003 
4004  while( (retval==1) && (ntries <= MAX_NEWT_RETRIES ) ) {
4005  // x[0] = x_orig[0] * (1. + ( (1.*rand())/(1.*RAND_MAX) - 0.5 ) );
4006  x[0] = x_orig[0];
4007  for( it=0; it<ntries; it++) x[0] *= 10.0;
4008  x1_of_x0( x[0],&x[1], wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra) ;
4009  // retval = general_newton_raphson(showmessages, lpflag, eomtype, x, n, USE_LINE_SEARCH, func_vsq, res_sq_vsq, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra,newtonstats) ;
4010  retval = general_newton_raphson(showmessages, lpflag, eomtype, x, n, USE_LINE_SEARCH, func_vsq2, res_sq_vsq2, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra,newtonstats) ;
4011 
4012 #if(!OPTIMIZED)
4013  if( ltrace ) {
4014  dualfprintf(fail_file, "find_root_2D_gen(): ntries, x[0,1] = %4i %21.15g %21.15g \n", ntries, x[0], x[1]);
4015 
4016  }
4017 #endif
4018  ntries++;
4019  }
4020 
4021 #if( MAX_NEWT_RETRIES > 0 )
4022  if( (ntries > MAX_NEWT_RETRIES) && (retval==1) ) {
4023  if( debugfail >= 2 ) {
4024  dualfprintf(fail_file, "find_root_2D_gen(): Bad exit value from general_newton_raphson() !! \n");
4025  dualfprintf(fail_file, "find_root_2D_gen(): ntries = %d , x[0] = %21.15g , x[1] = %21.15g \n", ntries, x[0], x[1]);
4026  }
4027  }
4028 #endif
4029 
4030  *xnew = x[0];
4031  return(retval);
4032 
4033 }
4034 
4035 
4036 
4037 
4038 /************************************************************
4039 
4040  general_newton_raphson():
4041 
4042  -- performs Newton-Rapshon method on an arbitrary system.
4043 
4044  -- inspired in part by Num. Rec.'s routine newt();
4045 
4046 *****************************************************************/
4047 static int general_newton_raphson(int showmessages, PFTYPE *lpflag, int eomtype, FTYPE x[], int n, int do_line_search,
4048  void (*funcd) (FTYPE [], FTYPE [], FTYPE [], FTYPE [], FTYPE [][NEWT_DIM], FTYPE *, FTYPE *, int, FTYPE *wglobal,FTYPE Bsq,FTYPE QdotB,FTYPE QdotBsq,FTYPE Qtsq,FTYPE Qdotn,FTYPE Qdotnp,FTYPE D,FTYPE Sc, int whicheos, FTYPE *EOSextra),
4049  FTYPE (*res_func) (FTYPE [], FTYPE *wglobal,FTYPE Bsq,FTYPE QdotB,FTYPE QdotBsq,FTYPE Qtsq,FTYPE Qdotn,FTYPE Qdotnp,FTYPE D,FTYPE Sc, int whicheos, FTYPE *EOSextra), 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)
4050 {
4051  FTYPE f, f_old, df, df_old, dx[NEWT_DIM], dx_old[NEWT_DIM], x_old[NEWT_DIM], x_older[NEWT_DIM], x_olderer[NEWT_DIM], resid[NEWT_DIM], norm[NEWT_DIM], jac[NEWT_DIM][NEWT_DIM];
4052  FTYPE errx, errx_old, errx_oldest, x_orig[NEWT_DIM];
4053  int n_iter, id, jd, i_extra, doing_extra;
4054  FTYPE randtmp, tmp;
4055  FTYPE dW,dvsq,vsq_old,vsq,W,W_old;
4056  FTYPE resid_norm, resid_check, grad_check;
4057 
4058  FTYPE res_func_val, res_func_old, res_func_new;
4059  FTYPE dn[NEWT_DIM], del_f[NEWT_DIM];
4060 
4061  FTYPE trueerror[NEWT_DIM],trueerrortotal;
4062 
4063  // Jon's modification
4064  FTYPE DAMPFACTOR[NEWT_DIM];
4065 
4066  int keep_iterating, i_increase, retval2,retval = 0;
4067  const int ltrace = 0;
4068  const int ltrace2 = 1;
4069  int it;
4070  int foundnan;
4071  int diddamp=0;
4072  int didcycle=0;
4073  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);
4074  int newtoniter;
4075 
4076 
4078  //
4079  // use inputs from newtonstats that sets up iterations
4080  // see global.structs.h: of_newtonstats{}
4081  //
4083  FTYPE NEWT_TOL_VAR;
4084  if(newtonstats->tryconv>0.0 && newtonstats->tryconv>NEWT_TOL) NEWT_TOL_VAR=newtonstats->tryconv;
4085  else NEWT_TOL_VAR=NEWT_TOL;
4086 
4087  FTYPE NEWT_TOL_ULTRAREL_VAR;
4088  if(newtonstats->tryconvultrarel>0.0 && newtonstats->tryconvultrarel>NEWT_TOL_ULTRAREL) NEWT_TOL_ULTRAREL_VAR=newtonstats->tryconvultrarel;
4089  else NEWT_TOL_ULTRAREL_VAR=NEWT_TOL_ULTRAREL;
4090 
4091  FTYPE MIN_NEWT_TOL_VAR;
4092  if(newtonstats->mintryconv>0.0 && newtonstats->mintryconv>MIN_NEWT_TOL) MIN_NEWT_TOL_VAR=newtonstats->mintryconv;
4093  else MIN_NEWT_TOL_VAR=MIN_NEWT_TOL;
4094 
4095  FTYPE MAX_NEWT_ITER_VAR;
4096  if(newtonstats->maxiter>=0 && newtonstats->maxiter<MAX_NEWT_ITER) MAX_NEWT_ITER_VAR=newtonstats->maxiter;
4097  else MAX_NEWT_ITER_VAR=MAX_NEWT_ITER;
4098 
4099  int EXTRA_NEWT_ITER_VAR;
4100  if(newtonstats->extra_newt_iter>=0 && newtonstats->extra_newt_iter<EXTRA_NEWT_ITER) EXTRA_NEWT_ITER_VAR=newtonstats->extra_newt_iter;
4101  else EXTRA_NEWT_ITER_VAR=EXTRA_NEWT_ITER;
4102 
4103  int EXTRA_NEWT_ITER_ULTRAREL_VAR;
4104  if(newtonstats->extra_newt_iter_ultrarel>=0 && newtonstats->extra_newt_iter_ultrarel<EXTRA_NEWT_ITER_ULTRAREL) EXTRA_NEWT_ITER_ULTRAREL_VAR=newtonstats->extra_newt_iter_ultrarel;
4105  else EXTRA_NEWT_ITER_ULTRAREL_VAR=EXTRA_NEWT_ITER_ULTRAREL;
4106 
4107 
4108  // setup error list
4109  FTYPE xlist[NEWT_DIM][MAX_NEWT_ITER+EXTRA_NEWT_ITER_ULTRAREL+EXTRA_NEWT_ITER],errxlist[MAX_NEWT_ITER+EXTRA_NEWT_ITER_ULTRAREL+EXTRA_NEWT_ITER];
4110  FTYPE dxlist[NEWT_DIM][MAX_NEWT_ITER+EXTRA_NEWT_ITER_ULTRAREL+EXTRA_NEWT_ITER],trueerrorlist[MAX_NEWT_ITER+EXTRA_NEWT_ITER_ULTRAREL+EXTRA_NEWT_ITER];
4111  int listi;
4112  for(it=0;it<n;it++){
4113  for(listi=0;listi<MAX_NEWT_ITER+EXTRA_NEWT_ITER_ULTRAREL+EXTRA_NEWT_ITER;listi++){
4114  xlist[it][listi]=errxlist[listi]=BIG;
4115  dxlist[it][listi]=trueerrorlist[listi]=BIG;
4116  }
4117  }
4118 
4119 
4120 
4121  // pick version of validate_x depending upon scheme
4122  pick_validate_x(eomtype, &ptr_validate_x);
4123 
4124 
4125  retval = 0;
4126 
4127 
4128  errx = 1. ;
4129  errx_old = 2.;
4130  df = df_old = f = f_old = 1.;
4131  i_extra = doing_extra = 0;
4132  for( id = 0; id < n ; id++) x_olderer[id]=x_older[id]=x_old[id] = x_orig[id] = x[id] ; // save initial guess as orig and old
4133  vsq_old = vsq = W = W_old = 0.;
4134 
4135  for( id = 0; id < n ; id++) DAMPFACTOR[id]=1.0; // Jon's mod: start out fast
4136  for( id = 0; id < n ; id++) dx_old[id]=VERYBIG;// Assume bad error at first (notice not normalized to x[id] since x[id] might be 0
4137 
4138 
4139 
4140 
4142  //
4143  /* Start the Newton-Raphson iterations : */
4144  //
4146  n_iter = 0;
4147  keep_iterating = 1;
4148  newtoniter=0;
4149  while( keep_iterating ) {
4150  (newtonstats->nstroke)++;
4151  (newtonstats->lntries)++;
4152 
4153  if((int)(newtonstats->lntries) > ITERDAMPSTART && diddamp==0 && didcycle==0){
4154  for( id = 0; id < n ; id++) DAMPFACTOR[id]=0.5; // Jon's mod: try to catch bad closed loop cycles that can occur, e.g., when cubic type behavior
4155  }
4156 
4157 
4158 
4159 #if(CRAZYDEBUG&&DEBUGINDEX)
4160  // DEBUG:
4161  if(ifileglobal==0 && jfileglobal==63 && nstep==9 && steppart==2){
4162  for(it=0;it<n;it++) dualfprintf(fail_file,"before funcd: x[%d]=%21.15g dx[%d]=%2.15g\n",it,x[it],it,dx[it]);
4163  }
4164 #endif
4165 
4166 
4167  /* returns with new dx, f, df */
4168  (*funcd) (x, dx, resid, norm, jac, &f, &df, n, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra);
4169 
4170  // get true error
4171  for(it=0;it<n;it++) trueerror[it]=fabs(resid[it]/(SMALL+fabs(norm[it])));
4172  trueerrortotal=0.0;
4173  for(it=0;it<n;it++) trueerrortotal = MAX(trueerrortotal,(1.0/(FTYPE)(n))*(fabs(resid[it])/(SMALL+fabs(norm[it]))));
4174 
4175  // store in list
4176  errxlist[n_iter] = errx;
4177  trueerrorlist[n_iter] = trueerrortotal;
4178  for(it=0;it<n;it++){
4179  dxlist[it][n_iter] = fabs(dx[it]);
4180  xlist[it][n_iter] = x[it];
4181  }
4182 
4183 
4184  // below now used currently
4185  id=0; if(resid[id] <=-VERYBIG || jac[0][id]==0.0){
4186  // get error and dx
4187  newtoniter++;
4188  }
4189 
4190 
4191 
4192  // DEBUG:
4193 #if(0)
4194  // if(nstep>=26070){
4195  // if(1||myid==5 && nstep==1 && steppart==0 && ifileglobal==19 && jfileglobal==15){
4196  for(it=0;it<n;it++) dualfprintf(fail_file,"lntries=%d after funcd: x[%d]=%26.20g dx[%d]=%26.20g f=%26.20g df=%26.20g errx=%26.20g trueerror=%26.20g (%26.20g %26.20g) diddamp=%d dampfactor=%26.20g didcycle=%d : %g %g %g %g %g %g %g %g %g\n",(int)(newtonstats->lntries),it,x[it],it,dx[it],f,df,errx,resid[it]/norm[it],NEWT_TOL_VAR,NEWT_TOL_ULTRAREL_VAR,diddamp,DAMPFACTOR[it],didcycle,wglobal[2],Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc);
4197  // }
4198  // }
4199 #endif
4200 
4201 #if(CRAZYDEBUG&&DEBUGINDEX)
4202  // DEBUG:
4203  if(ifileglobal==0 && jfileglobal==63 && nstep==9 && steppart==2){
4204  for(it=0;it<n;it++) dualfprintf(fail_file,"after funcd: x[%d]=%21.15g dx[%d]=%2.15g\n",it,x[it],it,dx[it]);
4205  }
4206 #endif
4207 
4208 
4209 #if(!OPTIMIZED)
4210  /* Check for bad untrapped divergences : */
4211  if( (finite(f)==0) || (finite(df)==0) ) {
4212  if( debugfail >= 2 ) {
4213  dualfprintf(fail_file,"general_newton_raphson(): nan encountered in f or df!! \n");
4214  dualfprintf(fail_file,"gnr nan(): f, df, x0, dx0 = %21.15g %21.15g %21.15g %21.15g \n", f,df,x[0],dx[0]);
4215  }
4216  return(1);
4217  }
4218 #endif
4219 
4220 
4221 #if(!OPTIMIZED)
4222  /* Randomly rescale Newton step to break out of iteration cycles: */
4223  if( ((n_iter+1) % CYCLE_BREAK_PERIOD) == 0 ) {
4224  randtmp = ( (1.*rand())/(1.*RAND_MAX) );
4225  for( id = 0; id < n ; id++) dx[id] *= randtmp;
4226  // for( id = 0; id < n ; id++) dx[id] *= ( (1.*rand())/(1.*RAND_MAX) );
4227  }
4228 #endif
4229 
4230  /* Save old values before calculating the new: */
4231  errx_oldest = errx_old;
4232  errx_old = errx;
4233  errx = 0.;
4234  f_old = f;
4235  for( id = 0; id < n ; id++) {
4236  x_olderer[id]=x_older[id]; //x_olderer contains 2 steps ago
4237  x_older[id]=x_old[id]; //x_older contains last step's W used by above funcd() to get resid and dx
4238  x_old[id] = x[id] ; // x_old contains just used W that was used to compute funcd() to get present resid and dx
4239  }
4240 
4241 
4242 
4243 
4244 
4245 
4246 
4247 
4249  //
4251 
4252  /* Make the newton step: */
4253  if( do_line_search == 1 ) {
4254 
4255  /* Compare the residual to its initial value */
4256  if( n_iter == 0 ) {
4257  resid_norm = 0.0e0;
4258  for( id = 0; id < n ; id++) {
4259  resid_norm += fabs(resid[id]);
4260  }
4261  resid_norm /= 1.0*n ;
4262  if( resid_norm == 0.0 ) resid_norm = 1.0;
4263  }
4264 
4265  for( id = 0; id < n ; id++) {
4266  tmp = 0.;
4267  for( jd = 0; jd < n ; jd++) {
4268  tmp += jac[jd][id] * resid[jd];
4269  }
4270  del_f[id] = tmp;
4271  }
4272  for( id = 0; id < n ; id++) {
4273  dn[id] = dx[id];
4274  }
4275 
4276  my_lnsrch(eomtype, n, x_old-1, f_old, del_f-1, dn-1, x-1, &f, TOL_LINE_STEP, SCALEMAX, &retval, res_func, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra);
4277 
4278  /* dx is needed for errx calculation below: */
4279  for( id = 0; id < n ; id++) {
4280  dx[id] = x[id] - x_old[id];
4281  }
4282 
4283 #if(!OPTIMIZED)
4284  if( ltrace ) {
4285  res_func_val = res_func(x, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc );
4286  res_func_old = res_func(x_old, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc );
4287  dualfprintf(fail_file,"gnr(): f_old, f, res_func_old, res_func_val = %21.15g %21.15g %21.15g %21.15g \n",
4288  f_old, f, res_func_old, res_func_val );
4289  dualfprintf(fail_file,"gnr(): x_old = ");
4290  for( id = 0; id < n ; id++) {
4291  dualfprintf(fail_file," %21.15g ",x_old[id]);
4292  }
4293  dualfprintf(fail_file,"\n ");
4294  dualfprintf(fail_file,"gnr(): x = ");
4295  for( id = 0; id < n ; id++) {
4296  dualfprintf(fail_file," %21.15g ",x[id]);
4297  }
4298  dualfprintf(fail_file,"\n ");
4299  dualfprintf(fail_file,"gnr(): dn = ");
4300  for( id = 0; id < n ; id++) {
4301  dualfprintf(fail_file," %21.15g ",dn[id]);
4302  }
4303  dualfprintf(fail_file,"\n ");
4304  dualfprintf(fail_file,"gnr(): del_f = ");
4305  for( id = 0; id < n ; id++) {
4306  dualfprintf(fail_file," %21.15g ",del_f[id]);
4307  }
4308  dualfprintf(fail_file,"\n ");
4309  }
4310 #endif
4311 
4312  /* Check to see if line search problem is because the residual vector is already small enough */
4313  if( retval == 1 ) {
4314  resid_check = 0.0e0;
4315  for( id = 0; id < n ; id++) {
4316  resid_check += fabs(resid[id]);
4317  }
4318  resid_check /= 1.0*n;
4319 
4320  if( resid_check <= resid_norm * NEWT_FUNC_TOL ) {
4321  retval = 0;
4322  }
4323  if( ltrace && retval ) {
4324  dualfprintf(fail_file,"general_newton_raphson(): retval, resid_check = %4i %21.15g \n",retval, resid_check);
4325  }
4326  }
4327  /* If initial Newton step is bad, then try again without line searching: */
4328  if( (retval == 2) && (USE_LINE_SEARCH == do_line_search) ) {
4329  if(debugfail>=2) dualfprintf(fail_file,"bad first step\n");
4330 #if(!OPTIMIZED)
4331  if( ltrace ) {
4332  dualfprintf(fail_file,"gnr(): bad first step: retval, f_old, f = %4i %21.15g %21.15g \n",retval,f_old,f);
4333  dualfprintf(fail_file,"gnr: doing recursive call, retval, errx = %4i %21.15g \n", retval, errx );
4334  }
4335 #endif
4336  retval = general_newton_raphson(showmessages, lpflag, eomtype, x_orig, n, ((do_line_search+1)%2), funcd, res_func, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra,newtonstats);
4337  for( id = 0; id < n ; id++) x[id] = x_orig[id] ;
4338  return( retval );
4339  }
4340 
4341  /* Check to see if it is trapped in a local minimum, i.e. gradient is too small */
4342  if( retval == 1 ) {
4343  grad_check = 0.0e0;
4344  for( id = 0; id < n ; id++) {
4345  resid_check = (x[id] == 0.) ? 1.0 : fabs(x[id]) ;
4346  grad_check += del_f[id] * resid_check ;
4347  }
4348  resid_check = (f == 0.) ? 1.0 : fabs(f) ;
4349  grad_check /= resid_check;
4350 
4351  /* Then we've most likely found a solution: */
4352  if( grad_check > GRADMIN ) {
4353  retval = -1;
4354  }
4355  else if( ltrace ) {
4356  dualfprintf(fail_file,"general_newton_raphson(): retval, grad_check = %4i %21.15g \n",retval, grad_check);
4357  }
4358  }
4359  }
4361  //
4363 
4364 
4365  else {
4366  /* don't use line search : */
4367  diddamp=0;
4368 
4369 #if(JONGENNEWTSTUFF)
4370  if(n_iter==0){
4371  for(id=0;id<n;id++) {
4372  // start history
4373  x_olderer[id]=x[id];
4374  x_older[id]=x[id];
4375  x_old[id]=x[id];
4376  dx_old[id]=dx[id];
4377  // xold[id]=x[id]; // original W
4378  // dxold[id]=dx[id]; // original dW
4379  // residold[id]=resid[id]; // original residual
4380  }
4381  }
4382  else{ // only meaningful to consider goodness of W if compute residual from present W, not previous W.
4383  for(id=0;id<n;id++){
4384  // if(sign(residold[id])!=sign(resid[id])){ // then passed over root -- bad in very relativistic case
4385  if(resid[id] <=-VERYBIG || jac[0][id]==0.0 ){ // then went too far such that v^2<0 or \tilde{u}^2<0 or Ss=NaN
4386  // if(1 || fabs(xold[0])>wglobal[2] || fabs(x_old[0])>wglobal[2]){ // always do now under more careful choice of when returns bad residual indicating problem
4387  if(1 || fabs(x_old[0])>wglobal[2]){ // always do now under more careful choice of when returns bad residual indicating problem
4388  // then ultrarelativistic and need to be more careful
4389  DAMPFACTOR[id]=0.5*DAMPFACTOR[id];
4390  x[id]=x_older[id]; // since current W is what gave bad residual and already saved it into x_old[0]
4391  dx[id]=dx_old[id]; // new dx is probably messed up (didn't yet save dx[0], so use dx_old[0]
4392 #if(PRODUCTION==0)
4393  if(showmessages&&debugfail>=3) dualfprintf(fail_file,"resetW: nstep=%ld steppart=%d :: lntries=%d :: id=%d :: x=%21.15g dx=%21.15g : x/D=%21.15g DAMPFACTOR=%21.15g errx=%21.15g eomtype=%d\n",nstep,steppart,(int)(newtonstats->lntries),id,x[id],dx[id],x[id]/D,DAMPFACTOR[id],newtonstats->lerrx,eomtype);
4394 #endif
4395  // SUPERGODMARK: Noticed that if Mathematica solution can be found but gives utsq<0 then this damping leads to nearly circular loop till max iterations.
4396  diddamp=1;
4397  }
4398  }
4399  } // end loop over id's
4400  if(diddamp){
4401  // Ensure this newly chosen "old" value gives reasonable result
4402  (*funcd) (x, dx, resid, norm, jac, &f, &df, n, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra); /* returns with new dx, f, df */
4403  for(id=0;id<n;id++){
4404  if(resid[id] <=-VERYBIG || jac[0][id]==0.0 ){ // then went too far such that v^2<0 or \tilde{u}^2<0 or Ss=nan or inf
4405 #if(PRODUCTION==0)
4406  if(debugfail>=2) dualfprintf(fail_file,"Old value failed too: resid[%d]=%26.20g jac[0][%d]=%26.20g\n",id,resid[id],id,jac[0][id]);
4407 #endif
4408  return(1); // then old fails to work too
4409  }
4410  // reset old values so will use this good value if immediately hit bad solution
4411  // This restarts history
4412  x_olderer[id]=x[id];
4413  x_older[id]=x[id];
4414  x_old[id]=x[id];
4415  dx_old[id]=dx[id];
4416  }
4417  }
4418 #if(1)
4419  else{
4420  // can only use x_olderer if 3 steps, but add 1 extra.
4421  // Only do if not just doing extra steps after already error is low (so avoid errx=0 case where will cycle).
4422  if(n_iter>4 && doing_extra==0){
4423  // see if cycling Newton steps with no change in errx
4424 #if(0)
4425  // DEBUG:
4426  if(nstep==19 && steppart==1){// && i==0 && j==0 && k==26){
4427  dualfprintf(fail_file,"errx=%26.20g errx_old=%26.20g errx_oldest=%26.20g dx=%26.20g dx_old=%26.20g x=%26.20g x_old=%26.20g x_older=%26.20g x_olderer=%26.20g\n",errx,errx_old,errx_oldest,dx[0],dx_old[0],x[0],x_old[0],x_older[0],x_olderer[0]);
4428  dualfprintf(fail_file,"differrx=%26.20g is=%d diffx=%26.20g is=%d\n",errx_old-errx_oldest,errx_old==errx_oldest,x[0]-x_olderer[0],x[0]==x_olderer[0]);
4429  }
4430 #endif
4431 
4432  if(newt_cyclecheck(n, trueerrortotal, errx, errx_old, errx_oldest, dx, dx_old, x, x_old, x_older, x_olderer, n_iter, xlist,dxlist, errxlist, trueerrorlist)){
4433  // if( newt_cyclecheck(n,errx,errx_old,errx_oldest,dx,dx_old,x,x_old,x_older,x_olderer) ){
4434  if(debugfail>=2) dualfprintf(fail_file,"Caught in cycle, trying to damp: n_iter=%d.\n",n_iter);
4435  // if hits another cycle, damp more
4436  for(id=0;id<n;id++) DAMPFACTOR[id]=0.1*DAMPFACTOR[id]; // severe damp to try to get good step.
4437  diddamp=1; // tells start of loop to not reset DAMPFACTOR
4438  didcycle++; // tells not to reset damp because having trouble and should stay damped
4439  }
4440  }// end if avoiding cycle
4441  }
4442 #endif
4443  }// end else if not first step
4444 
4445 
4446 #endif
4447 
4448 
4450  //
4451  // Step forward with x+dx
4452  //
4454  if(eomtype==EOMENTROPYGRMHD){
4455  // Function to find root of is like log(Wp) , so near root, Newton can jump to negative Wp. But for entropy method, this cannot be allowed without using very small DAMPFACTOR. So instead, somewhat bracket the root between original Wp and 0.0
4456  for( id = 0; id < n ; id++) {
4457  x[id] = MAX((x[id]+DAMPFACTOR[id]*dx[id]),(x[id]-0.0)/100.0);
4458 
4459  // DEBUG:
4460  // dualfprintf(fail_file,"id=%d newx=%21.15g using dx=%21.15g\n",id,x[id],dx[id]);
4461  }
4462  }
4463  else{
4464  for( id = 0; id < n ; id++) {
4465  x[id] +=DAMPFACTOR[id]*dx[id] ;
4466  }
4467  }
4468 
4470  //
4471  // undo damping if no longer damping
4472  // Otherwise too slow eventually, and otherwise error based upon dx thinks error is much smaller than really is.
4473  //
4475  for( id = 0; id < n ; id++) {
4476  if(diddamp==0) DAMPFACTOR[id]=MIN(1.0,2.0*DAMPFACTOR[id]);
4477  }
4478 
4479 
4480 #if(CRAZYDEBUG&&DEBUGINDEX)
4481  if(ifileglobal==0 && jfileglobal==63 && nstep==9 && steppart==2){
4482  dualfprintf(fail_file,"no line search: x[0]=%21.15g dx[0]=%21.15g\n",x[0],dx[0]);
4483  }
4484 #endif
4485 
4486  } /* End of "to do line search or not to do line search..." */
4487 
4488 
4489 
4490 
4491 
4492 
4493 
4495  //
4496  // Calculate the convergence criterion */
4497  //
4499 
4500  /* For the new criterion, always look at error in "W" : */
4501  // METHOD specific:
4502 
4503 #if( NEWCONVERGE == 1 )
4504 
4505 
4506 
4507 #if(JONGENNEWTSTUFF==0)
4508  // SUPERGODMARK: is wglobal always the right normalization for the error?
4509 
4510  // errx = (x[0]==0.) ? fabs(dx[0]) : fabs(dx[0]/x[0]);
4511  // errx = (x[0]==0.) ? fabs(dx[0]) : fabs(dx[0]/wglobal[0]);
4512 
4513  errx = (wglobal[0]==0.) ? fabs(dx[0]) : fabs(dx[0]/MAX(fabs(wglobal[0]),fabs(x[0])));
4514 
4515 #elif(JONGENNEWTSTUFF==1)
4516  // makes no sense to compute error from present dx for present x. The error in x[0] is unknown. errx here will be error in xold
4517  // actually error of previous x[0]
4518  // in this way, final funcd() call is actually an error check only, not an update to x[0]! So final solution returned should be x_old
4519  errx = (wglobal[0]==0.) ? fabs(dx_old[0]) : fabs(dx_old[0]/MAX(fabs(wglobal[0]),fabs(x[0])));
4520  // errx = (wglobal[0]==0.) ? fabs(dx_old[0]) : fabs(dx_old[0]/fabs(wglobal[0]));
4521 #endif
4522 
4523 
4524 
4525 
4526 #if(CRAZYDEBUG&&DEBUGINDEX)
4527  if(ifileglobal==0 && jfileglobal==31 && nstep==9 && steppart==2){
4528  dualfprintf(fail_file,"x[0]=%21.15g dx[0]=%21.15g\n",x[0],dx[0]);
4529  }
4530 #endif
4531 
4532 
4533  /* For the old criterion, look at errors in each indep. variable(s) (except for 5D) : */
4534 #else
4535  for( id = 0; id < n ; id++) {
4536  errx += (x[id]==0.) ? fabs(dx[id]) : fabs(dx[id]/x[id]);
4537  }
4538  errx /= 1.*n;
4539 #endif
4540 
4541 
4542  //dualfprintf(fail_file,"x[0]=%21.15g dx[0]=%21.15g df/f=%21.15g errx=%21.15g out of %21.15g or %21.15g\n",x[0],dx[0],df/f,errx,NEWT_TOL_VAR,newtonstats->tryconv);
4543 
4544 
4545  // save errx for information/debug purposes
4546  newtonstats->lerrx=errx;
4547 
4548 
4549  /****************************************/
4550  /* Make sure that the new x[] is physical : */
4551  /****************************************/
4552  // METHOD specific:
4553 
4554 
4555  validate_x(ptr_validate_x, x, x_old, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra) ;
4556 
4557 
4558  /****************************************/
4559  /* Check to see if we're in a infinite loop with error function: */
4560  /****************************************/
4561 #if( CHECK_FOR_STALL )
4562  if( ( (errx_old == errx) || (errx_oldest == errx) ) && (errx <= MIN_NEWT_TOL_VAR) ) errx = -errx;
4563 #endif
4564 
4565  /****************************************/
4566  /* If there's a problem with line search, then stop iterating: */
4567  /****************************************/
4568  if( (retval == 1) || (retval == -1) ) errx = -errx;
4569 
4570 
4571 #if(!OPTIMIZED)
4572  if( ltrace ) {
4573  dualfprintf(fail_file," general_newton_raphson(): niter,f_old,f,errx_old,errx = %4i %21.15g %21.15g %21.15g %21.15g\n",
4574  n_iter,f_old,f,errx_old,errx );
4575  dualfprintf(fail_file,"gnr(): x_old = ");
4576  for( id = 0; id < n ; id++) {
4577  dualfprintf(fail_file," %21.15g ",x_old[id]);
4578  }
4579  dualfprintf(fail_file,"\n ");
4580  dualfprintf(fail_file,"gnr(): x = ");
4581  for( id = 0; id < n ; id++) {
4582  dualfprintf(fail_file," %21.15g ",x[id]);
4583  }
4584  dualfprintf(fail_file,"\n ");
4585  dualfprintf(fail_file,"gnr(): dx = ");
4586  for( id = 0; id < n ; id++) {
4587  dualfprintf(fail_file," %21.15g ",dx[id]);
4588  }
4589  dualfprintf(fail_file,"\n ");
4590 
4591  }
4592 #endif
4593 
4594  /****************************************/
4595  /* Prepare for the next iteration, set the "old" variables: */
4596  /****************************************/
4597  for( id = 0; id < n ; id++) dx_old[id] = dx[id] ;
4598  f_old = f;
4599  df_old = df;
4600 
4601 
4602  /****************************************/
4603  /* If we've reached the tolerance level, then just do a few extra iterations before stopping */
4604  /****************************************/
4605 
4606 
4607 
4608  // see if reached tolerance
4609  if(newt_errorcheck(n, NEWT_TOL_VAR, NEWT_TOL_ULTRAREL_VAR, errx, trueerrortotal, x[0], dx[0], x, dx, wglobal) && (doing_extra == 0) && (newt_extracheck(n, EXTRA_NEWT_ITER_VAR, EXTRA_NEWT_ITER_ULTRAREL_VAR, errx,x[0],dx[0],x, dx, wglobal) > 0) ) {
4610  doing_extra = 1;
4611  }
4612 
4613  // see if repeating Newton steps with no change in errx
4614  if(doing_extra==0 && newt_repeatcheck(n,errx,errx_old,errx_oldest,dx,dx_old,x,x_old,x_older,x_olderer) ){
4615  if(debugfail>=2) dualfprintf(fail_file,"Repeat found\n");
4616  doing_extra = 1;
4617  }
4618 
4619 
4620  if( doing_extra == 1 ) i_extra++ ;
4621 
4622  if( (newt_errorcheck(n, NEWT_TOL_VAR, NEWT_TOL_ULTRAREL_VAR, errx, trueerrortotal, x[0], dx[0],x,dx, wglobal)&&(doing_extra == 0)) || (i_extra > newt_extracheck(n, EXTRA_NEWT_ITER_VAR, EXTRA_NEWT_ITER_ULTRAREL_VAR, errx,x[0],dx[0],x,dx,wglobal)) || (n_iter >= (MAX_NEWT_ITER_VAR-1)) ) {
4623  keep_iterating = 0;
4624  }
4625 
4626  // if cycling twice, then jump out and assume repeating cycle with no hope of improvement. Do so after CYCLESTOP steps
4627  if(n_iter>CYCLESTOP && didcycle>1){
4628  keep_iterating = 0;
4629  }
4630 
4631 
4632 
4633 #if(0)
4634  // check if error is not dropping sufficiently
4635 #define CHECKDECREASE0 5
4636 #define CHECKDECREASEAVGNUM 3 // should be less than CHECKDECREASE0
4637 #define CHECKDECFACTOR (0.5) // i.e. should drop by this factor compared to older average
4638  if(0&&newtoniter>CHECKDECREASE0){ // too aggressive, causes no solutions more often when could have gotten solution
4639  int ci;
4640  FTYPE avgerror;
4641  avgerror=0.0;
4642  for(ci=0;ci<CHECKDECREASEAVGNUM;ci++) if(errorlist[newtoniter-CHECKDECREASE0+ci]!=1.0) avgerror += errorlist[newtoniter-CHECKDECREASE0+ci]; // GODMARK: uses debug info
4643  avgerror/=(CHECKDECREASEAVGNUM);
4644  FTYPE currenterror;
4645  currenterror=errorlist[newtoniter];
4646  //
4647  // check both errors to ensure they are decreasing
4648  int cond=(currenterror>CHECKDECFACTOR*avgerror);
4649  if(cond){
4650  keep_iterating=0;
4651  if(debugfail>=2) dualfprintf(fail_file,"Error not decreasing properly: currenterror=%g avgerror=%g\n",currenterror,avgerror);
4652  }// end if current error too large compared to older average error
4653  }// end if large enough iterations so can check how error is trending, to avoid many iterations when error is not dropping enough to matter.
4654 #endif
4655 
4656 
4657 
4658 
4659 
4660  n_iter++;
4661 
4662 #if(CRAZYDEBUG&&DEBUGINDEX)
4663  //i=0 j=63 part=2 step=9 :: n_iter=6
4664  // if(ifileglobal==0 && jfileglobal==31 && nstep==9 && steppart==2){
4665  // if(ifileglobal==0 && jfileglobal==63 && nstep==9 && steppart==2){
4666  // if(n_iter>20){
4667  //if(ifileglobal>90 && CRAZYNEWCHECK){
4668  if(CRAZYNEWCHECK){
4669  // if(ifileglobal==0 && jfileglobal==63 && nstep==12 && steppart==0){
4670  // if(ifileglobal==0 && jfileglobal==0 && nstep==0 && steppart==0){
4671  dualfprintf(fail_file,"i=%d j=%d part=%d step=%ld :: n_iter=%d :: errx=%21.15g minerr=%21.15g :: x[0]=%21.15g dx[0]=%21.15g wglobal0=%21.15g\n",ifileglobal,jfileglobal,steppart,nstep,n_iter,errx,MIN_NEWT_TOL_VAR,x[0],dx[0],wglobal[0]);
4672  }
4673  // }
4674 #endif
4675 
4676  } // END of while(keep_iterating)
4677 
4678 
4679 
4680 
4681 
4682 
4683 
4684 
4685 
4686 
4687  // dualfprintf(fail_file,"n_iter=%d\n",n_iter);
4688 
4689 
4690 
4692  //
4693  // Check for bad untrapped divergences
4694  //
4696 
4697  foundnan=0; // assume no nan's
4698  for(id=0;id<n;id++){
4699  if(finite(x[id])==0) foundnan=1;
4700  }
4701 
4702  if( (finite(f)==0) || (finite(df)==0) ){ // deprecated: || (real_dimen_newton>=1)&&(finite(x[0])==0) || (real_dimen_newton>=2)&&(finite(x[1])==0)) {
4703  foundnan=1;
4704  }
4705 
4706 
4707  if(foundnan){
4708 #if(!OPTIMIZED || PRODUCTION==0)
4709  if(debugfail>=1){
4710  dualfprintf(fail_file,"naninf: f=%21.15g df=%21.15g\n",f,df);
4711  for(id=0;id<n;id++){
4712  dualfprintf(fail_file,"naninf: x[%d]=%21.15g\n",id,x[id]);
4713  }
4714  }
4715  if( debugfail >= 2 ) {
4716  dualfprintf(fail_file,"general_newton_raphson(): nan encountered in f or df!! \n");
4717  dualfprintf(fail_file,"gnr nan(): f, df, x0, dx0 = %21.15g %21.15g %21.15g %21.15g \n", f,df,x[0],dx[0]);
4718  }
4719 #endif
4720  *lpflag=UTOPRIMFAILNANRESULT;
4721  return(1);
4722  }
4723 
4724 
4725 
4726 
4727 
4729  //
4730  // Check error
4731  //
4733 
4734  // assign final result
4735 #if(JONGENNEWTSTUFF==1)
4736  for(id=0;id<n;id++) {
4737  x[id]=x_old[id]; // actually evaluated x
4738  }
4739 #endif
4740 
4741 
4742 
4743 
4744 
4745  if(newt_errorcheck(n, NEWT_TOL_VAR, NEWT_TOL_ULTRAREL_VAR, errx, trueerrortotal, x[0], dx[0], x,dx,wglobal) ){
4746 #if(DOHISTOGRAM)
4747  bin_newt_data(MAX_NEWT_ITER_VAR, errx, n_iter, 2, 0 );
4748 #endif
4749 #if(!OPTIMIZED)
4750  if(ltrace2) {
4751  dualfprintf(fail_file," totalcount = %d 2 %d %d %d %21.15g \n",n_iter,retval,do_line_search,i_extra, errx);
4752 
4753  }
4754  // dualfprintf(fail_file,"gnr retval5 = %4i \n", 0);
4755 #endif
4756 
4757 
4758  return(0);
4759  }
4760 
4761 
4762 
4763  if( (fabs(errx) <= MIN_NEWT_TOL_VAR) && newt_errorcheck(n, NEWT_TOL_VAR, NEWT_TOL_ULTRAREL_VAR, errx, trueerrortotal, x[0], dx[0], x,dx,wglobal) ){
4764 #if(DOHISTOGRAM)
4765  bin_newt_data(MAX_NEWT_ITER_VAR, errx, n_iter, 1, 0 );
4766 #endif
4767 #if(!OPTIMIZED)
4768  if(ltrace2) {
4769  dualfprintf(fail_file," totalcount = %d 1 %d %d %d %21.15g \n",n_iter,retval,do_line_search,i_extra,errx);
4770 
4771  }
4772  if(ltrace) {
4773  dualfprintf(fail_file,"general_newton_raphson(): found minimal solution \n");
4774 
4775  }
4776  // dualfprintf(fail_file,"gnr retval4 = %4i \n", 0);
4777 #endif
4778 
4779 
4780 #if(JONGENNEWTSTUFF==1)
4781  for(id=0;id<n;id++) {
4782  x[id]=x_old[id]; // actually evaluated x
4783  }
4784 #endif
4785 
4786  return(0);
4787  }
4788 
4789 
4790 
4791  if( fabs(errx) > MIN_NEWT_TOL_VAR){
4792  if( (do_line_search != USE_LINE_SEARCH) || (USE_LINE_SEARCH < 0) ) {
4793 #if(DOHISTOGRAM)
4794  bin_newt_data(MAX_NEWT_ITER_VAR, errx, n_iter, 0, 0 );
4795 #endif
4796 
4797  if(debugfail>=2){
4798  dualfprintf(fail_file,"fabs(errx)=%21.15g > MIN_NEWT_TOL_VAR=%21.15g n=%d n_iter=%d lntries=%d eomtype=%d\n",fabs(errx),MIN_NEWT_TOL_VAR,n,n_iter,newtonstats->lntries,eomtype);
4799  // if(DEBUGINDEX) dualfprintf(fail_file,"i=%d j=%d part=%d step=%ld :: n_iter=%d :: errx=%21.15g minerr=%21.15g :: x[0]=%21.15g dx[0]=%21.15g wglobal0=%21.15g\n",ifileglobal,jfileglobal,steppart,nstep,n_iter,errx,MIN_NEWT_TOL_VAR,x[0],dx[0],wglobal[0]);
4800  }
4801 
4802 
4803 
4804 #if(!OPTIMIZED)
4805 
4806  if(ltrace2) {
4807  dualfprintf(fail_file," totalcount = %d 0 %d %d %d %21.15g \n",n_iter,retval,do_line_search,i_extra,errx);
4808  }
4809  if(ltrace) {
4810  dualfprintf(fail_file,"general_newton_raphson(): did not find solution \n");
4811  if( retval == -1 ) {
4812  dualfprintf(fail_file,"general_newton_raphson(): lnsrch converged: x = ");
4813  for( id = 0; id < n ; id++) dualfprintf(fail_file," %21.15g ",x[id]);
4814  dualfprintf(fail_file,"\n");
4815  dualfprintf(fail_file,"general_newton_raphson(): lnsrch converged: x_old = ");
4816  for( id = 0; id < n ; id++) dualfprintf(fail_file," %21.15g ",x_old[id]);
4817  dualfprintf(fail_file,"\n");
4818  }
4819 
4820  }
4821  // dualfprintf(fail_file,"gnr retval2 = %4i \n", 1);
4822 #endif
4823  *lpflag=UTOPRIMFAILCONVRET+4;
4824  return(1);
4825  }
4826  else {
4827  /* If bad return and we tried line searching, try it without before giving up: */
4828  // dualfprintf(fail_file,"gnr: doing recursive call, do_line_search, retval, errx = %4i %4i %21.15g \n", do_line_search, retval, errx );
4829  //
4830  retval2 = general_newton_raphson(showmessages, lpflag, eomtype, x_orig, n, ((do_line_search+1)%2), funcd, res_func, wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra,newtonstats);
4831  for( id = 0; id < n ; id++) x[id] = x_orig[id] ;
4832  // dualfprintf(fail_file,"gnr retval3 = %4i \n", retval2);
4833 
4834  if(debugfail>=2) dualfprintf(fail_file,"fabs(errx) > MIN_NEWT_TOL_VAR: n_iter=%d lntries=%d\n",n_iter,newtonstats->lntries);
4835 
4836  return( retval2 );
4837  }
4838  }
4839 
4840 
4841 
4842 
4843 #if(!OPTIMIZED)
4844  dualfprintf(fail_file,"gnr retval6 = %4i \n", 0);
4845 #endif
4846  return(0);
4847 
4848 }
4849 
4850 
4851 
4852 
4853 
4854 
4855 
4856 
4857 
4858 
4859 
4860 
4861 
4863 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)
4864 {
4865 
4866  int ultrarel=0;
4867 
4868  ultrarel=( (fabs(wglobal[1])>wglobal[2]) || (fabs(x0)>wglobal[2]) );
4869 
4870  // see if dx's are too small to matter, despite condition on errx.
4871  int dxsmallcount=0;
4872  int iter;
4873  for(iter=0;iter<n;iter++) dxsmallcount += (fabs(dx[iter])<NUMEPSILON*x[iter]);
4874 
4875  if(!ultrarel){
4876  return(fabs(errx) <= NEWT_TOL_VAR || fabs(trueerrortotal) <= NEWT_TOL_VAR || dxsmallcount==n);
4877  }
4878  else{
4879  return(fabs(errx) <= NEWT_TOL_ULTRAREL_VAR || fabs(trueerrortotal) <= NEWT_TOL_ULTRAREL_VAR || dxsmallcount==n);
4880  }
4881 }
4882 
4883 
4884 
4886 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)
4887 {
4888  int repeat;
4889  int i;
4890 
4891  repeat=0;
4892  for(i=0;i<n;i++){
4893 
4894  if( (x[i]==x_olderer[i] && fabs(dx[i])==fabs(dx_old[i])) && (errx_oldest==errx_old) ){
4895  repeat++;
4896  }
4897  }
4898 
4899  if(repeat==n) return(1);
4900  else return(0);
4901 
4902 
4903 }
4904 
4905 // FTYPE xlist[NEWT_DIM][MAX_NEWT_ITER+EXTRA_NEWT_ITER_ULTRAREL+EXTRA_NEWT_ITER],errxlist[MAX_NEWT_ITER+EXTRA_NEWT_ITER_ULTRAREL+EXTRA_NEWT_ITER];
4906 // FTYPE dxlist[NEWT_DIM][MAX_NEWT_ITER+EXTRA_NEWT_ITER_ULTRAREL+EXTRA_NEWT_ITER],trueerrorlist[MAX_NEWT_ITER+EXTRA_NEWT_ITER_ULTRAREL+EXTRA_NEWT_ITER];
4907 
4908 
4910 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)
4911 {
4912  int i;
4913 
4914  int cycle=0;
4915  for(i=0;i<n;i++){
4916 
4917  // dualfprintf(fail_file,"%d %d %g %g %g %g\n",x[i]==x_olderer[i],errx_oldest==errx_old,x[i],x_olderer[i],errx_oldest,errx_old);
4918  if( (x[i]==x_olderer[i]) && (errx_oldest==errx_old && errx_old!=0.0) ){
4919  cycle++;
4920  }
4921  }
4922 
4923 
4924  int cycle2=0;
4925  int ii;
4926  for(ii=CYCLESTOP;ii<n_iter;ii++){
4927  if(errx_old!=0.0 && errx_old==errxlist[ii]) cycle2++;
4928  if(trueerror!=0.0 && trueerror==trueerrorlist[ii]) cycle2++;
4929  for(i=0;i<n;i++){
4930  if(x[i]==xlist[i][ii]) cycle2++;
4931  }
4932  }
4933 
4934 
4935 
4936 
4937  if(cycle>=n || cycle2>=NUMCYCLES*n) return(1);
4938  else return(0);
4939 
4940 
4941 }
4942 
4943 // see if done with extra iterations
4944 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)
4945 {
4946 
4947  int ultrarel=0;
4948 
4949  ultrarel=( (fabs(wglobal[1])>wglobal[2]) || (fabs(x0)>wglobal[2]) );
4950 
4951  // see if dx's are too small to matter, despite condition on errx.
4952  int dxsmallcount=0;
4953  int iter;
4954  for(iter=0;iter<n;iter++) dxsmallcount += (fabs(dx[iter])<NUMEPSILON*x[iter]);
4955 
4956  if(dxsmallcount) return(0);
4957  else{
4958  if(!ultrarel){
4959  return(EXTRA_NEWT_ITER_VAR);
4960  }
4961  else{
4962  return(EXTRA_NEWT_ITER_ULTRAREL_VAR);
4963  }
4964  }
4965 }
4966 
4967 
4968 #define ALF 1.0e-4
4969 
4970 static void my_lnsrch(int eomtype, int n, FTYPE xold[], FTYPE fold, FTYPE g[], FTYPE p[], FTYPE x[],
4971  FTYPE *f, FTYPE TOLX, FTYPE stpmax, int *check, FTYPE (*func) (FTYPE [], FTYPE *wglobal,FTYPE Bsq,FTYPE QdotB,FTYPE QdotBsq,FTYPE Qtsq,FTYPE Qdotn,FTYPE Qdotnp,FTYPE D,FTYPE Sc, int whicheos, FTYPE *EOSextra), FTYPE *wglobal,FTYPE Bsq,FTYPE QdotB,FTYPE QdotBsq,FTYPE Qtsq,FTYPE Qdotn,FTYPE Qdotnp,FTYPE D,FTYPE Sc, int whicheos, FTYPE *EOSextra)
4972 {
4973  int i;
4974  FTYPE a,alam,alam2,alamin,b,disc,f2,fold2,rhs1,rhs2,slope,sum,temp,test,tmplam;
4975  int icount=0;
4976  int bad_step;
4977  FTYPE bad_step_factor = 2.0;
4978  const int ltrace = 0;
4979  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);
4980 
4981 
4982 
4983  // pick version of validate_x depending upon scheme
4984  pick_validate_x(eomtype, &ptr_validate_x);
4985 
4986 
4987 
4988  *check=0;
4989  for (sum=0.0,i=1;i<=n;i++) sum += p[i]*p[i];
4990  sum=sqrt(sum);
4991  if (sum > stpmax)
4992  for (i=1;i<=n;i++) p[i] *= stpmax/sum;
4993  for (slope=0.0,i=1;i<=n;i++)
4994  slope += g[i]*p[i];
4995  test=0.0;
4996  for (i=1;i<=n;i++) {
4997  // temp=fabs(p[i])/max(fabs(xold[i]),1.0);
4998  temp= (xold[i] == 0.) ? fabs(p[i]) : fabs(p[i]/xold[i]);
4999  if (temp > test) test=temp;
5000  }
5001  alamin=TOLX/test;
5002 
5003 #if(!OPTIMIZED)
5004  if( ltrace ) {
5005  dualfprintf(fail_file,"my_lnsrch(): sum, slope, test, alamin = %21.15g %21.15g %21.15g %21.15g \n",sum,slope,test,alamin);
5006  }
5007 #endif
5008 
5009  alam=1.0;
5010  for (;;) {
5011  for (i=1;i<=n;i++) x[i]=xold[i]+alam*p[i];
5012 
5013  // METHOD specific:
5014  validate_x(ptr_validate_x, (x+1), (xold+1), wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra);
5015 
5016  *f=(*func)(x+1,wglobal,Bsq,QdotB,QdotBsq,Qtsq,Qdotn,Qdotnp,D,Sc,whicheos,EOSextra);
5017 
5018  bad_step = 0;
5019 
5020  if( finite(*f)==0 ) {
5021  bad_step = 1;
5022  }
5023 
5024 
5025  // if( bad_step ) alam /= bad_step_factor;
5026  // if (alam < alamin) bad_step = 0;
5027 
5028  if( bad_step ) {
5029  *check = 2;
5030 #if(!OPTIMIZED)
5031  dualfprintf(fail_file,"my_lnsrch(): bad_step = 1, f = %21.15g , alam = %21.15g \n", *f, alam);
5032  for( i = 1; i <= n; i++ ) {
5033  dualfprintf(fail_file,"my_lnsrch(): (xold[%d], x[%d], p[%d]) = %21.15g , %21.15g , %21.15g \n", i,i,i, xold[i],x[i],p[i]);
5034  }
5035 #endif
5036  return;
5037  }
5038 
5039  if (alam < alamin) {
5040  for (i=1;i<=n;i++) x[i]=xold[i];
5041  *check=1;
5042 #if(!OPTIMIZED)
5043  if( ltrace ) {
5044  dualfprintf(fail_file,"my_lnsrch(): alam < alamin: alam, alamin = %21.15g %21.15g \n", alam,alamin);
5045  }
5046 #endif
5047  return;
5048  }
5049  else if (*f <= fold+ALF*alam*slope) {
5050 #if(!OPTIMIZED)
5051  if( ltrace ) {
5052  dualfprintf(fail_file,"my_lnsrch(): good exit: alam, alamin, f, fold = %21.15g %21.15g %21.15g %21.15g \n", alam,alamin, *f, fold);
5053  }
5054 #endif
5055  return;
5056  }
5057  else {
5058  if (alam == 1.0) {
5059  tmplam = -slope/(2.0*(*f-fold-slope));
5060 #if(!OPTIMIZED)
5061  if( ltrace ) {
5062  dualfprintf(fail_file,"my_lnsrch(): setting tmplam!! tmplam, alam = %21.15g %21.15g !!\n", tmplam, alam);
5063  }
5064 #endif
5065  }
5066  else {
5067  rhs1 = *f-fold-alam*slope;
5068  rhs2=f2-fold2-alam2*slope;
5069  a=(rhs1/(alam*alam)-rhs2/(alam2*alam2))/(alam-alam2);
5070  b=(-alam2*rhs1/(alam*alam)+alam*rhs2/(alam2*alam2))/(alam-alam2);
5071  if (a == 0.0) tmplam = -slope/(2.0*b);
5072  else {
5073  disc=b*b-3.0*a*slope;
5074  if (disc<0.0) {
5075 #if(!OPTIMIZED)
5076  if( disc < -1.e-10 ) {
5077  if( ltrace ) dualfprintf(fail_file,"my_lnsrch(): Big Roundoff problem: disc = %21.15g \n", disc);
5078  }
5079 #endif
5080  disc = 0.;
5081  }
5082  else tmplam=(-b+sqrt(disc))/(3.0*a);
5083  }
5084  if (tmplam>0.5*alam)
5085  tmplam=0.5*alam;
5086 #if(!OPTIMIZED)
5087  if( ltrace ) {
5088  dualfprintf(fail_file,"my_lnsrch(): rhs1, rhs2, a, b, tmplam, alam = %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g !!\n",
5089  rhs1, rhs2, a, b, tmplam, alam );
5090  }
5091 #endif
5092  }
5093  }
5094  alam2=alam;
5095  f2 = *f;
5096  fold2=fold;
5097  alam=max(tmplam,0.1*alam);
5098 #if(!OPTIMIZED)
5099  if( ltrace ) {
5100  dualfprintf(fail_file,"my_lnsrch(): icount, alam, alam2, tmplam = %4i %21.15g %21.15g %21.15g \n",
5101  icount, alam, alam2, tmplam);
5102  }
5103 #endif
5104  icount++;
5105  }
5106 }
5107 #undef ALF
5108 
5109 
5110 #define N_CONV_TYPES 3
5111 #define N_NITER_BINS (MAX_NEWT_ITER + 3) // want to use MAX_NEWT_ITER_VAR, but not fixed, so constrain to be less.
5112 #define NBINS 200
5113 
5114 /****************************************************
5115 
5116  Used to gather statistics for the Newton solver during a disk run:
5117 
5118  -- bins are inclusive on their minimum bounds, and exlusive on their max.
5119 
5120 
5121  lerrx_[min,max] = the min./max. bounds of the errx historgram
5122  d_errx = errx histogram bin width
5123  n_beyond_range = number of solutions with errx max. beyond range speciified by lerrx_[min,max]
5124  n_conv_types = number of types of ways newton procedure can exit;
5125 
5126  print_now = 1 if you want to just print the histograms without counting a solution;
5127 
5128 *****************************************************/
5129 
5130 // OPENMPMARK: Not thread safe due to use of static variables, but not used except for debugging
5131 static void bin_newt_data(FTYPE MAX_NEWT_ITER_VAR, FTYPE errx, int niters, int conv_type, int print_now ) {
5132 
5133 
5134  /* General variables */
5135  int i, j;
5136  static int first_call = 1;
5137  static long int n_bin_calls = 0L;
5138  const long int n_output_freq = 128*128*100;
5139  FTYPE localerrx;
5140 
5141  /* Variables for the errx histogram */
5142  static const FTYPE localerrx_min = -23.;
5143  static const FTYPE localerrx_max = 4.;
5144  static FTYPE d_errx;
5145  static long int n_errx[N_CONV_TYPES][NBINS];
5146  static FTYPE xbin[NBINS];
5147  static long int n_beyond_range = 0L; /* Number of points that lie out of the bounds of our errx histogram */
5148  int ibin;
5149 
5150  /* Variables for the histogram of the number of newton iterations : */
5151  static long int n_niters[N_CONV_TYPES][N_NITER_BINS];
5152 
5153 
5154 #if(USEOPENMP)
5155  if(omp_in_parallel()){
5156  dualfprintf(fail_file,"bin_newt_data() called in parallel region\n");
5157  myexit(93586346);
5158  }
5159 #endif
5160 
5161 
5162  /* Clear arrays, set constants : */
5163  if( first_call ) {
5164  d_errx = ((localerrx_max - localerrx_min)/(1.*NBINS));
5165 
5166  for( i = 0; i < N_CONV_TYPES; i++ ) {
5167  for( j = 0; j < NBINS; j++ ) {
5168  n_errx[i][j] = 0;
5169  }
5170  for( j = 0; j < N_NITER_BINS; j++ ) {
5171  n_niters[i][j] = 0;
5172  }
5173  }
5174 
5175  for( j = 0; j < NBINS; j++ ) {
5176  xbin[j] = localerrx_min + (1.*j + 0.5)*d_errx;
5177  }
5178 
5179  first_call = 0;
5180  }
5181 
5182 
5183  if( print_now != 1 ) {
5184 
5185  /* Check validity of arguments : */
5186  errx = fabs(errx) ;
5187  localerrx = log10(errx + 1.0e-2*pow(10.,localerrx_min));
5188 
5189  if( (niters < 0) || (niters >= N_NITER_BINS) ) {
5190  dualfprintf(fail_file,"bin_newt_data(): bad value for niters = %d \n", niters );
5191  fflush(stdout);
5192  return;
5193  }
5194 
5195  /* Determine ibin */
5196  if( localerrx < localerrx_min ) {
5197  ibin = 0;
5198  n_beyond_range++ ;
5199  }
5200  else if( localerrx >= localerrx_max ) {
5201  ibin = NBINS - 1;
5202  n_beyond_range++ ;
5203  }
5204  else {
5205  ibin = (int) ( ( localerrx - localerrx_min ) / d_errx );
5206  }
5207 
5208  /* Tally this solution */
5209  n_errx[ conv_type][ibin]++;
5210  n_niters[conv_type][niters]++;
5211 
5212  }
5213 
5214 
5215  /* Print out the histograms periodically or when asked to : */
5216  if( print_now || ( (n_bin_calls % n_output_freq) == 0 ) ) {
5217 
5218  dualfprintf(log_file,"t = %21.15g , n_beyond_range = %ld , n_bin_calls = %ld \n", t, n_beyond_range, n_bin_calls);
5219 
5220  /* ERRX */
5221  dualfprintf(log_file,"ERRX-HISTOGRAM--ERRX-HISTOGRAM--ERRX-HISTOGRAM--ERRX-HISTOGRAM--ERRX-HISTOGRAM--ERRX-HISTOGRAM--\n");
5222  dualfprintf(log_file," x");
5223  for( j = 0; j < N_CONV_TYPES; j++ ) {
5224  dualfprintf(log_file," N%d",j);
5225  }
5226  dualfprintf(log_file,"\n");
5227 
5228  for( i = 0; i < NBINS; i++ ) {
5229  dualfprintf(log_file,"%21.15g ",xbin[i]);
5230  for( j = 0; j < N_CONV_TYPES; j++ ) {
5231  dualfprintf(log_file,"%13ld ", n_errx[j][i]);
5232  }
5233  dualfprintf(log_file,"\n");
5234  }
5235  dualfprintf(log_file,"ERRX-HISTOGRAM--ERRX-HISTOGRAM--ERRX-HISTOGRAM--ERRX-HISTOGRAM--ERRX-HISTOGRAM--ERRX-HISTOGRAM--\n");
5236 
5237 
5238  /* NITER */
5239 
5240  dualfprintf(log_file,"NITER-HISTOGRAM--NITER-HISTOGRAM--NITER-HISTOGRAM--NITER-HISTOGRAM--NITER-HISTOGRAM--NITER-HISTOGRAM--\n");
5241  dualfprintf(log_file," niter");
5242  for( j = 0; j < N_CONV_TYPES; j++ ) {
5243  dualfprintf(log_file," N%d",j);
5244  }
5245  dualfprintf(log_file,"\n");
5246 
5247  for( i = 0; i < N_NITER_BINS; i++ ) {
5248  dualfprintf(log_file,"%13d ", i);
5249  for( j = 0; j < N_CONV_TYPES; j++ ) {
5250  dualfprintf(log_file,"%13ld ", n_niters[j][i]);
5251  }
5252  dualfprintf(log_file,"\n");
5253  }
5254  dualfprintf(log_file,"NITER-HISTOGRAM--NITER-HISTOGRAM--NITER-HISTOGRAM--NITER-HISTOGRAM--NITER-HISTOGRAM--NITER-HISTOGRAM--\n");
5255 
5256  }
5257 
5258 
5259  if( print_now != 1 ) n_bin_calls++;
5260 
5261 
5262  return;
5263 
5264 }
5265 
5266 
5267 #undef N_CONV_TYPES
5268 #undef N_NITER_BINS
5269 #undef NBINS
5270 
5271 #undef CRAZYDEBUG
5272 
Generated on Fri May 20 2016 15:52:35 for HARM by   doxygen 1.8.3.1