Cantera  2.1.2
MultiJac.cpp
Go to the documentation of this file.
1 /**
2  * @file MultiJac.cpp Implementation file for class MultiJac
3  */
4 
5 /*
6  * Copyright 2002 California Institute of Technology
7  */
8 
10 using namespace std;
11 
12 namespace Cantera
13 {
14 
15 MultiJac::MultiJac(OneDim& r)
16  : BandMatrix(r.size(),r.bandwidth(),r.bandwidth())
17 {
18  m_size = r.size();
19  m_points = r.points();
20  m_resid = &r;
21  m_r1.resize(m_size);
22  m_ssdiag.resize(m_size);
23  m_mask.resize(m_size);
24  m_elapsed = 0.0;
25  m_nevals = 0;
26  m_age = 100000;
27  doublereal ff = 1.0;
28  while (1.0 + ff != 1.0) {
29  ff *= 0.5;
30  }
31  m_atol = sqrt(ff);
32  m_rtol = 1.0e-5;
33 }
34 
35 void MultiJac::updateTransient(doublereal rdt, integer* mask)
36 {
37  for (size_t n = 0; n < m_size; n++) {
38  value(n,n) = m_ssdiag[n] - mask[n]*rdt;
39  }
40 }
41 
42 void MultiJac::incrementDiagonal(int j, doublereal d)
43 {
44  m_ssdiag[j] += d;
45  value(j,j) = m_ssdiag[j];
46 }
47 
48 void MultiJac::eval(doublereal* x0, doublereal* resid0, doublereal rdt)
49 {
50  m_nevals++;
51  clock_t t0 = clock();
52  bfill(0.0);
53  size_t n, m, ipt=0, j, nv, mv, iloc;
54  doublereal rdx, dx, xsave;
55 
56  for (j = 0; j < m_points; j++) {
57  nv = m_resid->nVars(j);
58  for (n = 0; n < nv; n++) {
59 
60  // perturb x(n)
61  xsave = x0[ipt];
62  dx = m_atol + fabs(xsave)*m_rtol;
63  x0[ipt] = xsave + dx;
64  dx = x0[ipt] - xsave;
65  rdx = 1.0/dx;
66 
67  // calculate perturbed residual
68  m_resid->eval(j, x0, DATA_PTR(m_r1), rdt, 0);
69 
70  // compute nth column of Jacobian
71  for (size_t i = j - 1; i != j+2; i++) {
72  if (i != npos && i < m_points) {
73  mv = m_resid->nVars(i);
74  iloc = m_resid->loc(i);
75  for (m = 0; m < mv; m++) {
76  value(m+iloc,ipt) = (m_r1[m+iloc]
77  - resid0[m+iloc])*rdx;
78  }
79  }
80  }
81  x0[ipt] = xsave;
82  ipt++;
83  }
84 
85  }
86 
87  for (n = 0; n < m_size; n++) {
88  m_ssdiag[n] = value(n,n);
89  }
90 
91  m_elapsed += double(clock() - t0)/CLOCKS_PER_SEC;
92  m_age = 0;
93 }
94 
95 } // namespace
void eval(size_t j, double *x, double *r, doublereal rdt=-1.0, int count=1)
Evaluate the multi-domain residual function.
Definition: OneDim.cpp:236
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:173
size_t loc(size_t jg)
Location in the solution vector of the first component of global point jg.
Definition: OneDim.h:106
OneDim * m_resid
Residual evaluator for this jacobian.
Definition: MultiJac.h:79
doublereal & value(size_t i, size_t j)
Return a changeable reference to element (i,j).
Definition: BandMatrix.cpp:138
void eval(doublereal *x0, doublereal *resid0, double rdt)
Evaluate the Jacobian at x0.
Definition: MultiJac.cpp:48
void bfill(doublereal v=0.0)
Fill or zero the matrix.
Definition: BandMatrix.cpp:116
#define DATA_PTR(vec)
Creates a pointer to the start of the raw data for a vector.
Definition: ct_defs.h:36
size_t nVars(size_t jg)
Number of solution components at global point jg.
Definition: OneDim.h:98