Cantera  3.0.0
Loading...
Searching...
No Matches
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
9namespace Cantera
10{
11
12MultiJac::MultiJac(OneDim& r)
13 : BandMatrix(r.size(),r.bandwidth(),r.bandwidth())
14{
15 m_size = r.size();
16 m_points = r.points();
17 m_resid = &r;
18 m_r1.resize(m_size);
19 m_ssdiag.resize(m_size);
20 m_mask.resize(m_size);
21}
22
23void MultiJac::updateTransient(double rdt, integer* mask)
24{
25 for (size_t n = 0; n < m_size; n++) {
26 value(n,n) = m_ssdiag[n] - mask[n]*rdt;
27 }
28}
29
30void MultiJac::incrementDiagonal(int j, double d)
31{
32 m_ssdiag[j] += d;
33 value(j,j) = m_ssdiag[j];
34}
35
36void MultiJac::eval(double* x0, double* resid0, double rdt)
37{
38 m_nevals++;
39 clock_t t0 = clock();
40 bfill(0.0);
41 size_t ipt=0;
42
43 for (size_t j = 0; j < m_points; j++) {
44 size_t nv = m_resid->nVars(j);
45 for (size_t n = 0; n < nv; n++) {
46 // perturb x(n); preserve sign(x(n))
47 double xsave = x0[ipt];
48 double dx;
49 if (xsave >= 0) {
50 dx = xsave*m_rtol + m_atol;
51 } else {
52 dx = xsave*m_rtol - m_atol;
53 }
54 x0[ipt] = xsave + dx;
55 dx = x0[ipt] - xsave;
56 double rdx = 1.0/dx;
57
58 // calculate perturbed residual
59 m_resid->eval(j, x0, m_r1.data(), rdt, 0);
60
61 // compute nth column of Jacobian
62 for (size_t i = j - 1; i != j+2; i++) {
63 if (i != npos && i < m_points) {
64 size_t mv = m_resid->nVars(i);
65 size_t iloc = m_resid->loc(i);
66 for (size_t m = 0; m < mv; m++) {
67 value(m+iloc,ipt) = (m_r1[m+iloc] - resid0[m+iloc])*rdx;
68 }
69 }
70 }
71 x0[ipt] = xsave;
72 ipt++;
73 }
74 }
75
76 for (size_t n = 0; n < m_size; n++) {
77 m_ssdiag[n] = value(n,n);
78 }
79
80 m_elapsed += double(clock() - t0)/CLOCKS_PER_SEC;
81 m_age = 0;
82}
83
84} // namespace
Namespace for the Cantera kernel.
Definition AnyMap.cpp:564
const size_t npos
index returned by functions to indicate "no position"
Definition ct_defs.h:195