Convert the conservative variables {D, m, sigma_c} (where sigma_c = D*sigma is the conserved entropy) to primitive variable using a Newton-Raphson/Bisection scheme.
26 double v2, v2max, v2min, acc = 1.e-13;
27 double h, W, W2, W3, S2_W2;
29 double lor, lor2, rho, sigma, th;
30 double dv2_dW, dW_dlor;
31 double fv2, dfv2, dv2, dv2old;
32 double m2,
B2, S, S2,
D;
46 v2 = 0.5*(v2max+v2min);
50 lor2 = 1.0/(1.0 - v2);
61 th = sigma*x/sqrt(1.0 + 3.0*sigma*x);
62 h = 2.5*th + sqrt(2.25*th*th + 1.0);
65 scrh = -2.0*x/(3.0*lor)*sigma*(2.0*sigma*x - 3.0*th*th);
66 scrh /= 2.0*th*(1.0 + 3.0*sigma*
x);
67 dW_dlor = W/lor + D*lor*(5.0*h - 8.0*th)/(2.0*h - 5.0*th)*
scrh;
75 fv2 = v2 - (m2 + S2_W2*(scrh + W))/(scrh*
scrh);
77 dv2_dW = S2*(3.0*W*scrh + B2*
B2) + m2*W3;
78 dv2_dW *= -2.0/(W3*scrh*scrh*
scrh);
80 dfv2 = 1.0 - 0.5*dv2_dW*dW_dlor*lor2*lor;
83 if ((((v2-v2max)*dfv2-fv2)*((v2-v2min)*dfv2-fv2) > 0.0)
84 || (fabs(2.*fv2) > fabs(dv2old*dfv2))) {
86 dv2 = 0.5*(v2max-v2min);
88 if (v2min == v2) done = 1;
94 if (temp == v2) done = 1;
97 if (fabs(dv2) < acc*v2 || fabs(fv2) < acc) done = 1;
99 lor2 = 1.0/(1.0 - v2);
110 x = pow(rho,2.0/3.0);
111 th = sigma*x/sqrt(1.0 + 3.0*sigma*x);
112 h = 2.5*th + sqrt(2.25*th*th + 1.0);
115 scrh = -2.0*x/(3.0*lor)*sigma*(2.0*sigma*x - 3.0*th*th);
116 scrh /= 2.0*th*(1.0 + 3.0*sigma*
x);
117 dW_dlor = W/lor + D*lor*(5.0*h-8.0*th)/(2.0*h-5.0*th)*
scrh;
128 fv2 = v2 - (m2 + S2_W2*(scrh + W))/(scrh*
scrh);
130 if (fv2 < 0.0) v2min = v2;
133 dv2_dW = S2*(3.0*W*scrh + B2*
B2) + m2*W3;
134 dv2_dW *= -2.0/(W3*scrh*scrh*
scrh);
136 dfv2 = 1.-0.5*dv2_dW*dW_dlor*lor2*lor;
142 par->
prs = sigma*x*rho;
149 print (
"! EntropySolve: too many iterations,%d, ",k);
154 if (par->
prs < 0.0 || sigma < 0.0 || W < 0.0) {
156 print (
"! EntropySolve: negative pressure, p = %12.6e\n",par->
prs);
167 #if RMHD_REDUCED_ENERGY == YES
170 + 0.5*(1.0 + v2)*B2 - 0.5*S2_W2;
172 scrh = 2.25*th*th + v2;
173 par->
E = par->
prs*(2.5*lor2 - 1.0)
174 + D*lor2/(lor*sqrt(2.25*th*th + 1.0) + 1.0)*scrh
175 + 0.5*(1.0 + v2)*B2 - 0.5*S2_W2;
178 par->
E = W - par->
prs + 0.5*(1.0 + v2)*B2 - 0.5*S2_W2;
double m2
Square of total momentum (input).
double D
Lab density (input).
void print(const char *fmt,...)
double prs
Thermal pressure (output).
double sigma_c
Conserved entropy (input).
double B2
Square of magnetic field (input).
double E
Total energy (input).