Cantera  2.5.1
MultiJac.cpp
Go to the documentation of this file.
1 //! @file MultiJac.cpp Implementation file for class MultiJac
2 
3 // This file is part of Cantera. See License.txt in the top-level directory or
4 // at https://cantera.org/license.txt for license and copyright information.
5 
7 #include <ctime>
8 
9 using namespace std;
10 
11 namespace Cantera
12 {
13 
14 MultiJac::MultiJac(OneDim& r)
15  : BandMatrix(r.size(),r.bandwidth(),r.bandwidth())
16 {
17  m_size = r.size();
18  m_points = r.points();
19  m_resid = &r;
20  m_r1.resize(m_size);
21  m_ssdiag.resize(m_size);
22  m_mask.resize(m_size);
23  m_elapsed = 0.0;
24  m_nevals = 0;
25  m_age = 100000;
26  m_atol = sqrt(std::numeric_limits<double>::epsilon());
27  m_rtol = 1.0e-5;
28 }
29 
30 void MultiJac::updateTransient(doublereal rdt, integer* mask)
31 {
32  for (size_t n = 0; n < m_size; n++) {
33  value(n,n) = m_ssdiag[n] - mask[n]*rdt;
34  }
35 }
36 
37 void MultiJac::incrementDiagonal(int j, doublereal d)
38 {
39  m_ssdiag[j] += d;
40  value(j,j) = m_ssdiag[j];
41 }
42 
43 void MultiJac::eval(doublereal* x0, doublereal* resid0, doublereal rdt)
44 {
45  m_nevals++;
46  clock_t t0 = clock();
47  bfill(0.0);
48  size_t ipt=0;
49 
50  for (size_t j = 0; j < m_points; j++) {
51  size_t nv = m_resid->nVars(j);
52  for (size_t n = 0; n < nv; n++) {
53  // perturb x(n); preserve sign(x(n))
54  double xsave = x0[ipt];
55  double dx;
56  if (xsave >= 0) {
57  dx = xsave*m_rtol + m_atol;
58  } else {
59  dx = xsave*m_rtol - m_atol;
60  }
61  x0[ipt] = xsave + dx;
62  dx = x0[ipt] - xsave;
63  double rdx = 1.0/dx;
64 
65  // calculate perturbed residual
66  m_resid->eval(j, x0, m_r1.data(), rdt, 0);
67 
68  // compute nth column of Jacobian
69  for (size_t i = j - 1; i != j+2; i++) {
70  if (i != npos && i < m_points) {
71  size_t mv = m_resid->nVars(i);
72  size_t iloc = m_resid->loc(i);
73  for (size_t m = 0; m < mv; m++) {
74  value(m+iloc,ipt) = (m_r1[m+iloc] - resid0[m+iloc])*rdx;
75  }
76  }
77  }
78  x0[ipt] = xsave;
79  ipt++;
80  }
81  }
82 
83  for (size_t n = 0; n < m_size; n++) {
84  m_ssdiag[n] = value(n,n);
85  }
86 
87  m_elapsed += double(clock() - t0)/CLOCKS_PER_SEC;
88  m_age = 0;
89 }
90 
91 } // namespace
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:188
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:264