Cantera  2.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 #include <ctime>
11 
12 using namespace std;
13 
14 namespace Cantera
15 {
16 
17 MultiJac::MultiJac(OneDim& r)
18  : BandMatrix(r.size(),r.bandwidth(),r.bandwidth())
19 {
20  m_size = r.size();
21  m_points = r.points();
22  m_resid = &r;
23  m_r1.resize(m_size);
24  m_ssdiag.resize(m_size);
25  m_mask.resize(m_size);
26  m_elapsed = 0.0;
27  m_nevals = 0;
28  m_age = 100000;
29  doublereal ff = 1.0;
30  while (1.0 + ff != 1.0) {
31  ff *= 0.5;
32  }
33  m_atol = sqrt(ff);
34  m_rtol = 1.0e-5;
35 }
36 
37 void MultiJac::updateTransient(doublereal rdt, integer* mask)
38 {
39  for (size_t n = 0; n < m_size; n++) {
40  value(n,n) = m_ssdiag[n] - mask[n]*rdt;
41  }
42 }
43 
44 void MultiJac::incrementDiagonal(int j, doublereal d)
45 {
46  m_ssdiag[j] += d;
47  value(j,j) = m_ssdiag[j];
48 }
49 
50 void MultiJac::eval(doublereal* x0, doublereal* resid0, doublereal rdt)
51 {
52  m_nevals++;
53  clock_t t0 = clock();
54  bfill(0.0);
55  size_t n, m, ipt=0, j, nv, mv, iloc;
56  doublereal rdx, dx, xsave;
57 
58  for (j = 0; j < m_points; j++) {
59  nv = m_resid->nVars(j);
60  for (n = 0; n < nv; n++) {
61 
62  // perturb x(n)
63  xsave = x0[ipt];
64  dx = m_atol + fabs(xsave)*m_rtol;
65  x0[ipt] = xsave + dx;
66  dx = x0[ipt] - xsave;
67  rdx = 1.0/dx;
68 
69  // calculate perturbed residual
70  m_resid->eval(j, x0, DATA_PTR(m_r1), rdt, 0);
71 
72  // compute nth column of Jacobian
73  for (size_t i = j - 1; i != j+2; i++) {
74  if (i != npos && i < m_points) {
75  mv = m_resid->nVars(i);
76  iloc = m_resid->loc(i);
77  for (m = 0; m < mv; m++) {
78  value(m+iloc,ipt) = (m_r1[m+iloc]
79  - resid0[m+iloc])*rdx;
80  }
81  }
82  }
83  x0[ipt] = xsave;
84  ipt++;
85  }
86 
87  }
88 
89  for (n = 0; n < m_size; n++) {
90  m_ssdiag[n] = value(n,n);
91  }
92 
93  m_elapsed += double(clock() - t0)/CLOCKS_PER_SEC;
94  m_age = 0;
95 }
96 
97 } // 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:232
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:165
size_t loc(size_t jg)
Location in the solution vector of the first component of global point jg.
Definition: OneDim.h:105
OneDim * m_resid
Residual evaluator for this Jacobian.
Definition: MultiJac.h:78
doublereal & value(size_t i, size_t j)
Return a changeable reference to element (i,j).
Definition: BandMatrix.cpp:129
void eval(doublereal *x0, doublereal *resid0, double rdt)
Evaluate the Jacobian at x0.
Definition: MultiJac.cpp:50
void bfill(doublereal v=0.0)
Fill or zero the matrix.
Definition: BandMatrix.cpp:107
#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:97