CellModel.cpp 10.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
///////////////////////////////////////////////////////////////////////////////
//
// File CellModel.cpp
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// License for the specific language governing rights and limitations under
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: Cell model base class.
//
///////////////////////////////////////////////////////////////////////////////

36
37
#include <LibUtilities/BasicUtils/VmathArray.hpp>

38
39
#include <CardiacEPSolver/CellModels/CellModel.h>

40
41
42
#include <StdRegions/StdNodalTriExp.h>
#include <LibUtilities/LinearAlgebra/Blas.hpp>

43
44
namespace Nektar
{
45
46
47
48
49
50
51
52
    CellModelFactory& GetCellModelFactory()
    {
        typedef Loki::SingletonHolder<CellModelFactory,
            Loki::CreateUsingNew,
            Loki::NoDestroy > Type;
        return Type::Instance();
    }

53
54
55
56
57
58
59
60
61
62
63
    /**
     * @class CellModel
     *
     * The CellModel class and derived classes implement a range of cell model
     * ODE systems. A cell model comprises a system of ion concentration
     * variables and zero or more gating variables. Gating variables are
     * time-integrated using the Rush-Larsen method and for each variable y,
     * the corresponding y_inf and tau_y value is computed by Update(). The tau
     * values are stored in separate storage to inarray/outarray, #m_gates_tau.
     */

64
65
66
    /**
     * Cell model base class constructor.
     */
67
68
    CellModel::CellModel(const LibUtilities::SessionReaderSharedPtr& pSession,
                         const MultiRegions::ExpListSharedPtr& pField)
69
    {
70
71
        m_session = pSession;
        m_field = pField;
72
        m_lastTime = 0.0;
73
        m_substeps = pSession->GetParameter("Substeps");
74
        m_nvar = 0;
75
76
77
78
79
80
81
82
83
        m_useNodal = false;

        // Number of points in nodal space is the number of coefficients
        // in modified basis
        std::set<enum StdRegions::ExpansionType> s;
        for (unsigned int i = 0; i < m_field->GetNumElmts(); ++i)
        {
            s.insert(m_field->GetExp(i)->DetExpansionType());
        }
84

85
        // Use nodal projection if only triangles
86
        if (s.size() == 1 && (s.count(StdRegions::eTriangle) == 1 || s.count(StdRegions::eTetrahedron) == 1))
87
88
89
90
91
92
93
94
95
96
97
        {
            m_useNodal = true;
        }

        // ---------------------------
        // Move to nodal points
        if (m_useNodal)
        {
            m_nq = pField->GetNcoeffs();
            int order = m_field->GetExp(0)->GetBasis(0)->GetNumModes();
            int nvar = 1;
98

99
100
101
102
103
104
105
            // Set up a nodal tri
            LibUtilities::BasisKey B0(
                LibUtilities::eModified_A, order,
                LibUtilities::PointsKey(order, LibUtilities::eGaussLobattoLegendre));
            LibUtilities::BasisKey B1(
                LibUtilities::eModified_B, order,
                LibUtilities::PointsKey(order, LibUtilities::eGaussRadauMAlpha1Beta0));
106
107
108
            LibUtilities::BasisKey B2(
                LibUtilities::eModified_C, order,
                LibUtilities::PointsKey(order, LibUtilities::eGaussRadauMAlpha2Beta0));
109

110
111
            m_nodalTri = MemoryManager<StdRegions::StdNodalTriExp>::AllocateSharedPtr(
                        B0, B1, LibUtilities::eNodalTriEvenlySpaced);
112
113
            m_nodalTet = MemoryManager<StdRegions::StdNodalTetExp>::AllocateSharedPtr(
                        B0, B1, B2, LibUtilities::eNodalTetEvenlySpaced);
114
115
116
117
118
        }
        else
        {
            m_nq = pField->GetTotPoints();
        }
119
120
    }

121
122
123
124
125
126
127
128
129
130
131
132

    /**
     * Initialise the cell model. Allocate workspace and variable storage.
     */
    void CellModel::Initialise()
    {
        ASSERTL1(m_nvar > 0, "Cell model must have at least 1 variable.");

        m_cellSol = Array<OneD, Array<OneD, NekDouble> >(m_nvar);
        m_wsp = Array<OneD, Array<OneD, NekDouble> >(m_nvar);
        for (unsigned int i = 0; i < m_nvar; ++i)
        {
133
134
            m_cellSol[i] = Array<OneD, NekDouble>(m_nq);
            m_wsp[i] = Array<OneD, NekDouble>(m_nq);
135
136
137
138
        }
        m_gates_tau = Array<OneD, Array<OneD, NekDouble> >(m_gates.size());
        for (unsigned int i = 0; i < m_gates.size(); ++i)
        {
139
            m_gates_tau[i] = Array<OneD, NekDouble>(m_nq);
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
        }

        v_SetInitialConditions();
    }

    /**
     * Integrates the cell model for one PDE time-step. Cell model is
     * sub-stepped.
     *
     * Ion concentrations and membrane potential are integrated using forward
     * Euler, while gating variables are integrated using the Rush-Larsen
     * scheme.
     */
    void CellModel::TimeIntegrate(
            const Array<OneD, const Array<OneD, NekDouble> > &inarray,
                  Array<OneD,       Array<OneD, NekDouble> > &outarray,
            const NekDouble time)
    {
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
        int phys_offset = 0;
        int coef_offset = 0;
        int nvar = inarray.num_elements();
        Array<OneD, NekDouble> tmp;

        // ---------------------------
        // Check nodal temp array set up
        if (m_useNodal)
        {
            if (!m_nodalTmp.num_elements())
            {
                m_nodalTmp = Array<OneD, Array<OneD, NekDouble> >(nvar);
                for (unsigned int k = 0; k < nvar; ++k)
                {
                    m_nodalTmp[k] = Array<OneD, NekDouble>(m_nq);
                }
            }

            // Move to nodal points
            for (unsigned int k = 0; k < nvar; ++k)
            {
                for (unsigned int i = 0; i < m_field->GetNumElmts(); ++i)
                {
                    phys_offset = m_field->GetPhys_Offset(i);
                    coef_offset = m_field->GetCoeff_Offset(i);
183
184
185
186
187
188
189
190
191
192
                    if (m_field->GetExp(0)->DetExpansionType() == StdRegions::eTriangle)
                    {
                        m_field->GetExp(0)->FwdTrans(inarray[k] + phys_offset, m_nodalTri->UpdateCoeffs());
                        m_nodalTri->ModalToNodal(m_nodalTri->GetCoeffs(), tmp=m_nodalTmp[k]+coef_offset);
                    }
                    else
                    {
                        m_field->GetExp(0)->FwdTrans(inarray[k] + phys_offset, m_nodalTet->UpdateCoeffs());
                        m_nodalTet->ModalToNodal(m_nodalTet->GetCoeffs(), tmp=m_nodalTmp[k]+coef_offset);
                    }
193
194
195
196
197
198
199
200
201
202
203
204
                }
            }
            // Copy new transmembrane potential into cell model
            Vmath::Vcopy(m_nq, m_nodalTmp[0], 1, m_cellSol[0], 1);
        }
        else
        {
            // Copy new transmembrane potential into cell model
            Vmath::Vcopy(m_nq, inarray[0], 1, m_cellSol[0], 1);
        }
        // -------------------------

205
        NekDouble delta_t = (time - m_lastTime)/m_substeps;
206
207
208


        // Perform substepping
209
        for (unsigned int i = 0; i < m_substeps - 1; ++i)
210
211
212
        {
            Update(m_cellSol, m_wsp, time);
            // Voltage
213
            Vmath::Svtvp(m_nq, delta_t, m_wsp[0], 1, m_cellSol[0], 1, m_cellSol[0], 1);
214
215
216
            // Ion concentrations
            for (unsigned int j = 0; j < m_concentrations.size(); ++j)
            {
217
                Vmath::Svtvp(m_nq, delta_t, m_wsp[m_concentrations[j]], 1, m_cellSol[m_concentrations[j]], 1, m_cellSol[m_concentrations[j]], 1);
218
219
220
221
            }
            // Gating variables: Rush-Larsen scheme
            for (unsigned int j = 0; j < m_gates.size(); ++j)
            {
222
223
224
225
                Vmath::Sdiv(m_nq, -delta_t, m_gates_tau[j], 1, m_gates_tau[j], 1);
                Vmath::Vexp(m_nq, m_gates_tau[j], 1, m_gates_tau[j], 1);
                Vmath::Vsub(m_nq, m_cellSol[m_gates[j]], 1, m_wsp[m_gates[j]], 1, m_cellSol[m_gates[j]], 1);
                Vmath::Vvtvp(m_nq, m_cellSol[m_gates[j]], 1, m_gates_tau[j], 1, m_wsp[m_gates[j]], 1, m_cellSol[m_gates[j]], 1);
226
227
228
229
230
231
232
            }
        }

        // Perform final cell model step
        Update(m_cellSol, m_wsp, time);

        // Output dV/dt from last step but integrate remaining cell model vars
233
234
235
236
237
238
239
240
241
        // Transform cell model I_total from nodal to modal space
        if (m_useNodal)
        {
            for (unsigned int k = 0; k < nvar; ++k)
            {
                for (unsigned int i = 0; i < m_field->GetNumElmts(); ++i)
                {
                    int phys_offset = m_field->GetPhys_Offset(i);
                    int coef_offset = m_field->GetCoeff_Offset(i);
242
243
244
245
246
247
248
249
250
251
                    if (m_field->GetExp(0)->DetExpansionType() == StdRegions::eTriangle)
                    {
                        m_nodalTri->NodalToModal(m_wsp[k]+coef_offset, m_nodalTri->UpdateCoeffs());
                        m_field->GetExp(0)->BwdTrans(m_nodalTri->GetCoeffs(), tmp=outarray[k] + phys_offset);
                    }
                    else
                    {
                        m_nodalTet->NodalToModal(m_wsp[k]+coef_offset, m_nodalTet->UpdateCoeffs());
                        m_field->GetExp(0)->BwdTrans(m_nodalTet->GetCoeffs(), tmp=outarray[k] + phys_offset);
                    }
252
253
254
255
256
257
258
                }
            }
        }
        else
        {
            Vmath::Vcopy(m_nq, m_wsp[0], 1, outarray[0], 1);
        }
259
260
261
262

        // Ion concentrations
        for (unsigned int j = 0; j < m_concentrations.size(); ++j)
        {
263
            Vmath::Svtvp(m_nq, delta_t, m_wsp[m_concentrations[j]], 1, m_cellSol[m_concentrations[j]], 1, m_cellSol[m_concentrations[j]], 1);
264
265
266
267
268
        }

        // Gating variables: Rush-Larsen scheme
        for (unsigned int j = 0; j < m_gates.size(); ++j)
        {
269
270
271
272
            Vmath::Sdiv(m_nq, -delta_t, m_gates_tau[j], 1, m_gates_tau[j], 1);
            Vmath::Vexp(m_nq, m_gates_tau[j], 1, m_gates_tau[j], 1);
            Vmath::Vsub(m_nq, m_cellSol[m_gates[j]], 1, m_wsp[m_gates[j]], 1, m_cellSol[m_gates[j]], 1);
            Vmath::Vvtvp(m_nq, m_cellSol[m_gates[j]], 1, m_gates_tau[j], 1, m_wsp[m_gates[j]], 1, m_cellSol[m_gates[j]], 1);
273
        }
274
275

        m_lastTime = time;
276
    }
277
}