Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
AlievPanfilov.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File AlievPanfilov.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Aliev-Panfilov phenomological cell model.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <iostream>
37 #include <string>
38 
41 
42 namespace Nektar
43 {
44  /**
45  * Registers the class with the Factory.
46  */
49  "AlievPanfilov",
51  "Phenomological model of canine cardiac electrophysiology.");
52 
55  const MultiRegions::ExpListSharedPtr& pField)
56  : CellModel(pSession, pField)
57  {
58  pSession->LoadParameter("k", m_k, 0.0);
59  pSession->LoadParameter("a", m_a, 0.0);
60  pSession->LoadParameter("mu1", m_mu1, 0.0);
61  pSession->LoadParameter("mu2", m_mu2, 0.0);
62  pSession->LoadParameter("eps", m_eps, 0.0);
63 
68 
69  m_nvar = 2;
70  m_concentrations.push_back(1);
71  }
72 
73 
75  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
76  Array<OneD, Array<OneD, NekDouble> >&outarray,
77  const NekDouble time)
78  {
79  // inarray[0] holds initial physical u values throughout
80  // inarray[1] holds initial physical v values throughout
81 
82  // compute u^2: m_u = u*u
83  Vmath::Vmul(m_nq, &inarray[0][0], 1, &inarray[0][0], 1, &m_uu[0], 1);
84 
85  // compute u^3: m_u = u*u*u
86  Vmath::Vmul(m_nq, &inarray[0][0], 1, &m_uu[0], 1, &m_uuu[0], 1);
87 
88  // --------------------------------------
89  // Compute reaction term f(u,v)
90  // --------------------------------------
91 // if (m_spatialParameters->Exists("a"))
92 // {
93 // Vmath::Vmul(m_nq, &m_spatialParameters->GetData("a")->GetPhys()[0], 1,
94 // &inarray[0][0], 1, &m_tmp1[0], 1);
95 //
96 // Vmath::Vvtvm(m_nq, &m_spatialParameters->GetData("a")->GetPhys()[0], 1,
97 // &m_uu[0], 1, &m_tmp1[0], 1, &m_tmp1[0], 1);
98 //
99 // Vmath::Svtvm(m_nq, -1.0, &m_uu[0], 1, &m_tmp1[0], 1, &m_tmp1[0], 1);
100 // }
101 // else
102 // {
103  // Ru = au
104  Vmath::Smul(m_nq, m_a, &inarray[0][0], 1, &m_tmp1[0], 1);
105  // Ru = (-1-a)u*u + au
106  Vmath::Svtvp(m_nq, (-1.0-m_a), &m_uu[0], 1, &m_tmp1[0], 1,
107  &m_tmp1[0], 1);
108 // }
109  // Ru = u*u*u - (1+a)u*u + au
110  Vmath::Vadd(m_nq, &m_uuu[0], 1, &m_tmp1[0], 1, &m_tmp1[0], 1);
111  // Ru = k(u*u*u - (1+a)u*u + au)
112 // if (m_spatialParameters->Exists("k"))
113 // {
114 // Vmath::Vmul(m_nq, &m_spatialParameters->GetData("k")->GetPhys()[0], 1,
115 // &m_tmp1[0], 1, &m_tmp1[0], 1);
116 // }
117 // else
118 // {
119  Vmath::Smul(m_nq, m_k, &m_tmp1[0], 1, &m_tmp1[0], 1);
120 // }
121 
122  // Ru = k(u*u*u - (1+a)u*u + au) + I_stim
123  Vmath::Vadd(m_nq, &outarray[0][0], 1, &m_tmp1[0], 1, &outarray[0][0], 1);
124 
125  // Ru = k(u*u*u - (1+a)u*u + au) + uv + I_stim
126  Vmath::Vvtvp(m_nq, &inarray[0][0], 1, &inarray[1][0], 1, &m_tmp1[0], 1,
127  &outarray[0][0], 1);
128  // Ru = -k(u*u*u - (1+a)u*u + au) - uv - I_stim
129  Vmath::Neg(m_nq, &outarray[0][0], 1);
130 
131 
132  // --------------------------------------
133  // Compute reaction term g(u,v)
134  // --------------------------------------
135  // tmp2 = mu2 + u
136  Vmath::Sadd(m_nq, m_mu2, &inarray[0][0], 1, &m_tmp2[0], 1);
137 
138  // tmp2 = v/(mu2 + u)
139  Vmath::Vdiv(m_nq, &inarray[1][0], 1, &m_tmp2[0], 1, &m_tmp2[0], 1);
140 
141  // tmp2 = mu1*v/(mu2 + u)
142  Vmath::Smul(m_nq, m_mu1, &m_tmp2[0], 1, &m_tmp2[0], 1);
143 
144  // tmp1 = Eps + mu1*v/(mu2+u)
145  Vmath::Sadd(m_nq, m_eps, &m_tmp2[0], 1, &m_tmp2[0], 1);
146 
147  // tmp1 = (-a-1) + u
148 // if (m_spatialParameters->Exists("a"))
149 // {
150 // Vmath::Vsub(m_nq, &inarray[0][0], 1,
151 // &m_spatialParameters->GetData("a")->GetPhys()[0], 1,
152 // &m_tmp1[0], 1);
153 //
154 // Vmath::Sadd(m_nq, -1.0, &inarray[0][0], 1, &m_tmp1[0], 1);
155 // }
156 // else
157 // {
158  Vmath::Sadd(m_nq, (-m_a-1), &inarray[0][0], 1, &m_tmp1[0], 1);
159 // }
160 
161  // tmp1 = k(u-a-1)
162 // if (m_spatialParameters->Exists("k"))
163 // {
164 // Vmath::Vmul(m_nq, &m_spatialParameters->GetData("k")->GetPhys()[0], 1,
165 // &m_tmp1[0], 1, &m_tmp1[0], 1);
166 // }
167 // else
168 // {
169  Vmath::Smul(m_nq, m_k, &m_tmp1[0], 1, &m_tmp1[0], 1);
170 // }
171 
172  // tmp1 = ku(u-a-1) + v
173  Vmath::Vvtvp(m_nq, &inarray[0][0], 1, &m_tmp1[0], 1, &inarray[1][0], 1,
174  &m_tmp1[0], 1);
175 
176  // tmp1 = -ku(u-a-1)-v
177  Vmath::Neg(m_nq, &m_tmp1[0], 1);
178 
179  // outarray = [Eps + mu1*v/(mu2+u)] * [-ku(u-a-1)-v]
180  Vmath::Vmul(m_nq, &m_tmp1[0], 1, &m_tmp2[0], 1, &outarray[1][0], 1);
181  }
182 
183  /**
184  *
185  */
187  {
188  SolverUtils::AddSummaryItem(s, "Cell model","Aliev-Panfilov");
194  }
195 
196 
197  /**
198  *
199  */
201  {
202  Vmath::Fill(m_nq, 0.0, m_cellSol[0], 1);
203  Vmath::Fill(m_nq, 0.0, m_cellSol[1], 1);
204  }
205 }
NekDouble m_a
Trigger parameter a.
Definition: AlievPanfilov.h:75
Array< OneD, NekDouble > m_uu
Temporary space for storing when computing reaction term.
Definition: AlievPanfilov.h:86
Array< OneD, NekDouble > m_tmp2
Workspace for computing reaction term.
Definition: AlievPanfilov.h:92
int m_nq
Number of physical points.
Definition: CellModel.h:117
CellModelAlievPanfilov(const LibUtilities::SessionReaderSharedPtr &pSession, const MultiRegions::ExpListSharedPtr &pField)
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:47
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:46
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:485
Cell model base class.
Definition: CellModel.h:65
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:442
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
Definition: Vmath.cpp:241
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
static CellModelSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const MultiRegions::ExpListSharedPtr &pField)
Creates an instance of this class.
Definition: AlievPanfilov.h:48
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:213
std::vector< int > m_concentrations
Indices of cell model variables which are concentrations.
Definition: CellModel.h:139
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
Definition: Misc.cpp:50
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Array< OneD, Array< OneD, NekDouble > > m_cellSol
Cell model solution variables.
Definition: CellModel.h:126
virtual void v_GenerateSummary(SummaryList &s)
int m_nvar
Number of variables in cell model (inc. transmembrane voltage)
Definition: CellModel.h:119
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:396
double NekDouble
NekDouble m_eps
Restitution parameter .
Definition: AlievPanfilov.h:83
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
Definition: Vmath.cpp:315
CellModelFactory & GetCellModelFactory()
Definition: CellModel.cpp:47
static std::string className
Name of class.
Definition: AlievPanfilov.h:55
NekDouble m_k
Scaling parameter k.
Definition: AlievPanfilov.h:77
NekDouble m_mu1
Restitution parameter .
Definition: AlievPanfilov.h:79
Array< OneD, NekDouble > m_tmp1
Workspace for computing reaction term.
Definition: AlievPanfilov.h:90
Array< OneD, NekDouble > m_uuu
Temporary space for storing when computing reaction term.
Definition: AlievPanfilov.h:88
NekDouble m_mu2
Restitution parameter .
Definition: AlievPanfilov.h:81
virtual void v_Update(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:299
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:183
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215