5 void Jacobian (
double *v,
double *rhs,
double **dfdy)
40 double dmu_dX[
NIONS], N, mu;
41 double T, dLdX, dCdX, dRdX,
eps,
scrh;
44 double *L, *La, *Lb, *Lc;
45 double *C, *Ca, *Cb, *Cc;
46 double *R, *Ra, *Rb, *Rc;
47 double *
X, *dnel_dX, *de;
48 double Unit_Time, E_cost;
49 static double **delta;
60 for (k = 0; k <
NIONS; k++){
61 for (l = 0; l <
NIONS; l++){
62 delta[
k][l] = (k == l ? 1.0:0.0);
84 for (k = 0; k < n - 1; k++){
98 for (l = 0; l < n - 1; l++) {
99 dfdy[0][l] = - (R[0] + C[0])*delta[0][l] +
104 for (k = 1; k < n - 1; k++) {
105 for (l = 0; l < n - 1; l++) {
106 dLdX = La[
k]*dnel_dX[l] + Lb[
k]*delta[0][l] + Lc[
k]*delta[1][l];
107 dCdX = Ca[
k]*dnel_dX[l] + Cb[
k]*delta[0][l] + Cc[
k]*delta[1][l];
108 dfdy[
k][l] = L[
k]*delta[k - 1][l] + dLdX*X[k - 1] - C[
k]*delta[
k][l] - dCdX*X[
k];
111 for (k = 1; k < n - 2; k++) {
112 for (l = 0; l < n - 1; l++) {
113 dRdX = Ra[
k]*dnel_dX[l] + Rb[
k]*delta[0][l] + Rc[
k]*delta[1][l];
114 dfdy[
k][l] += R[
k]*delta[k + 1][l] + dRdX*X[k + 1];
117 for (k = 0; k < n - 1; k++) {
118 for (l = 0; l < n - 1; l++) {
119 dfdy[
k][l] *= Unit_Time;
126 for (l = 0; l < n - 1; l++){
131 for (nv = 0; nv <
NIONS; nv++){
150 for (nv = 0; nv <
NVAR; nv++){
151 vm[nv] = vp[nv] = v[nv];
153 vp[PRS] = v[PRS]*(1.0 +
eps);
154 vm[PRS] = v[PRS]*(1.0 -
eps);
161 for (k = 0; k < n - 1; k++){
162 dfdp[
k] = (rhs_p[k +
NFLX] - rhs_m[k +
NFLX])/(2.0*eps*v[PRS]);
163 dfdy[
k][n - 1] = dfdp[
k];
165 dfdp[n - 1] = (rhs_p[PRS] - rhs_m[PRS])/(2.0*eps*v[PRS]);
166 dfdy[n - 1][n - 1] = dfdp[n - 1];
171 for (k = 0; k <
n ; k++) {
172 for (l = 0; l < n - 1; l++) {
173 dfdy[
k][l] += dfdp[
k]*scrh*dmu_dX[l];
#define UNIT_DENSITY
Unit density in gr/cm^3.
void Radiat(double *, double *)
#define UNIT_VELOCITY
Unit velocity in cm/sec.
#define UNIT_LENGTH
Unit Length in cm.
void Jacobian(real *v, real *rhs, real **dfdy)
double MeanMolecularWeight(double *V)
#define ARRAY_2D(nx, ny, type)