#include "pluto.h"
#include "cooling_defs.h"
Go to the source code of this file.
|
void | Jacobian (double *v, double *rhs, double **dfdy) |
|
void Jacobian |
( |
double * |
v, |
|
|
double * |
rhs, |
|
|
double ** |
dfdy |
|
) |
| |
Compute the jacobian J(k,l) = dfdy
k = row index l = column index
J(0,0) J(0,1) ... J(0, n-1) J(1,0) J(1,1) .... J(1, n-1) . . . . . . . . . J(n-1,0) .... J(n-1, n-1)
or,
+--------------------—+
- | |
- | |
- | |
- dX'/dX | dX'/dp |
- (JXX) | (JXp) |
- | |
- | | +-----------—+-----—+
- dp'/dX | dp'/dp |
- (JpX) | Jpp | +--------------------—+
Definition at line 5 of file jacobian.c.
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.
double MeanMolecularWeight(double *V)
#define ARRAY_2D(nx, ny, type)