Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
CellModel.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File CellModel.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: Cell model base class.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
37 
39 
41 //#include <LibUtilities/LinearAlgebra/Blas.hpp>
42 
43 using namespace std;
44 
45 namespace Nektar
46 {
48  {
49  typedef Loki::SingletonHolder<CellModelFactory,
50  Loki::CreateUsingNew,
51  Loki::NoDestroy > Type;
52  return Type::Instance();
53  }
54 
55  /**
56  * @class CellModel
57  *
58  * The CellModel class and derived classes implement a range of cell model
59  * ODE systems. A cell model comprises a system of ion concentration
60  * variables and zero or more gating variables. Gating variables are
61  * time-integrated using the Rush-Larsen method and for each variable y,
62  * the corresponding y_inf and tau_y value is computed by Update(). The tau
63  * values are stored in separate storage to inarray/outarray, #m_gates_tau.
64  */
65 
66  /**
67  * Cell model base class constructor.
68  */
69  CellModel::CellModel(const LibUtilities::SessionReaderSharedPtr& pSession,
70  const MultiRegions::ExpListSharedPtr& pField)
71  {
72  m_session = pSession;
73  m_field = pField;
74  m_lastTime = 0.0;
75  m_substeps = pSession->GetParameter("Substeps");
76  m_nvar = 0;
77  m_useNodal = false;
78 
79  // Number of points in nodal space is the number of coefficients
80  // in modified basis
81  std::set<enum LibUtilities::ShapeType> s;
82  for (unsigned int i = 0; i < m_field->GetNumElmts(); ++i)
83  {
84  s.insert(m_field->GetExp(i)->DetShapeType());
85  }
86 
87  // Use nodal projection if only triangles
88  if (s.size() == 1 && (s.count(LibUtilities::eTriangle) == 1 ||
89  s.count(LibUtilities::eTetrahedron) == 1))
90  {
91  // This is disabled for now as it causes problems at high order.
92  // m_useNodal = true;
93  }
94 
95  // ---------------------------
96  // Move to nodal points
97  if (m_useNodal)
98  {
99  m_nq = pField->GetNcoeffs();
100  int order = m_field->GetExp(0)->GetBasis(0)->GetNumModes();
101 
102  // Set up a nodal tri
112 
117  }
118  else
119  {
120  m_nq = pField->GetTotPoints();
121  }
122  }
123 
124 
125  /**
126  * Initialise the cell model. Allocate workspace and variable storage.
127  */
128  void CellModel::Initialise()
129  {
130  ASSERTL1(m_nvar > 0, "Cell model must have at least 1 variable.");
131 
132  m_cellSol = Array<OneD, Array<OneD, NekDouble> >(m_nvar);
133  m_wsp = Array<OneD, Array<OneD, NekDouble> >(m_nvar);
134  for (unsigned int i = 0; i < m_nvar; ++i)
135  {
136  m_cellSol[i] = Array<OneD, NekDouble>(m_nq);
137  m_wsp[i] = Array<OneD, NekDouble>(m_nq);
138  }
139  m_gates_tau = Array<OneD, Array<OneD, NekDouble> >(m_gates.size());
140  for (unsigned int i = 0; i < m_gates.size(); ++i)
141  {
142  m_gates_tau[i] = Array<OneD, NekDouble>(m_nq);
143  }
144 
145  if (m_session->DefinesFunction("CellModelInitialConditions"))
146  {
147  LoadCellModel();
148  }
149  else
150  {
151  v_SetInitialConditions();
152  }
153  }
154 
155  /**
156  * Integrates the cell model for one PDE time-step. Cell model is
157  * sub-stepped.
158  *
159  * Ion concentrations and membrane potential are integrated using forward
160  * Euler, while gating variables are integrated using the Rush-Larsen
161  * scheme.
162  */
163  void CellModel::TimeIntegrate(
164  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
165  Array<OneD, Array<OneD, NekDouble> > &outarray,
166  const NekDouble time)
167  {
168  int phys_offset = 0;
169  int coef_offset = 0;
170  int nvar = inarray.num_elements();
172 
173  // ---------------------------
174  // Check nodal temp array set up
175  if (m_useNodal)
176  {
177  if (!m_nodalTmp.num_elements())
178  {
179  m_nodalTmp = Array<OneD, Array<OneD, NekDouble> >(nvar);
180  for (unsigned int k = 0; k < nvar; ++k)
181  {
182  m_nodalTmp[k] = Array<OneD, NekDouble>(m_nq);
183  }
184  }
185 
186  // Move to nodal points
187  Array<OneD, NekDouble> tmpCoeffs(max(m_nodalTri->GetNcoeffs(), m_nodalTet->GetNcoeffs()));
188 
189  for (unsigned int k = 0; k < nvar; ++k)
190  {
191  for (unsigned int i = 0; i < m_field->GetNumElmts(); ++i)
192  {
193  phys_offset = m_field->GetPhys_Offset(i);
194  coef_offset = m_field->GetCoeff_Offset(i);
195  if (m_field->GetExp(0)->DetShapeType() == LibUtilities::eTriangle)
196  {
197  m_field->GetExp(0)->FwdTrans(inarray[k] + phys_offset, tmpCoeffs);
198  m_nodalTri->ModalToNodal(tmpCoeffs, tmp=m_nodalTmp[k]+coef_offset);
199  }
200  else
201  {
202  m_field->GetExp(0)->FwdTrans(inarray[k] + phys_offset, tmpCoeffs);
203  m_nodalTet->ModalToNodal(tmpCoeffs, tmp=m_nodalTmp[k]+coef_offset);
204  }
205  }
206  }
207  // Copy new transmembrane potential into cell model
208  Vmath::Vcopy(m_nq, m_nodalTmp[0], 1, m_cellSol[0], 1);
209  }
210  else
211  {
212  // Copy new transmembrane potential into cell model
213  Vmath::Vcopy(m_nq, inarray[0], 1, m_cellSol[0], 1);
214  }
215  // -------------------------
216 
217  NekDouble delta_t = (time - m_lastTime)/m_substeps;
218 
219 
220  // Perform substepping
221  for (unsigned int i = 0; i < m_substeps - 1; ++i)
222  {
223  Update(m_cellSol, m_wsp, time);
224  // Voltage
225  Vmath::Svtvp(m_nq, delta_t, m_wsp[0], 1, m_cellSol[0], 1, m_cellSol[0], 1);
226  // Ion concentrations
227  for (unsigned int j = 0; j < m_concentrations.size(); ++j)
228  {
229  Vmath::Svtvp(m_nq, delta_t, m_wsp[m_concentrations[j]], 1, m_cellSol[m_concentrations[j]], 1, m_cellSol[m_concentrations[j]], 1);
230  }
231  // Gating variables: Rush-Larsen scheme
232  for (unsigned int j = 0; j < m_gates.size(); ++j)
233  {
234  Vmath::Sdiv(m_nq, -delta_t, m_gates_tau[j], 1, m_gates_tau[j], 1);
235  Vmath::Vexp(m_nq, m_gates_tau[j], 1, m_gates_tau[j], 1);
236  Vmath::Vsub(m_nq, m_cellSol[m_gates[j]], 1, m_wsp[m_gates[j]], 1, m_cellSol[m_gates[j]], 1);
237  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);
238  }
239  }
240 
241  // Perform final cell model step
242  Update(m_cellSol, m_wsp, time);
243 
244  // Output dV/dt from last step but integrate remaining cell model vars
245  // Transform cell model I_total from nodal to modal space
246  if (m_useNodal)
247  {
248  Array<OneD, NekDouble> tmpCoeffs(max(m_nodalTri->GetNcoeffs(), m_nodalTet->GetNcoeffs()));
249 
250  for (unsigned int k = 0; k < nvar; ++k)
251  {
252  for (unsigned int i = 0; i < m_field->GetNumElmts(); ++i)
253  {
254  int phys_offset = m_field->GetPhys_Offset(i);
255  int coef_offset = m_field->GetCoeff_Offset(i);
256  if (m_field->GetExp(0)->DetShapeType() == LibUtilities::eTriangle)
257  {
258  m_nodalTri->NodalToModal(m_wsp[k]+coef_offset, tmpCoeffs);
259  m_field->GetExp(0)->BwdTrans(tmpCoeffs, tmp=outarray[k] + phys_offset);
260  }
261  else
262  {
263  m_nodalTet->NodalToModal(m_wsp[k]+coef_offset, tmpCoeffs);
264  m_field->GetExp(0)->BwdTrans(tmpCoeffs, tmp=outarray[k] + phys_offset);
265  }
266  }
267  }
268  }
269  else
270  {
271  Vmath::Vcopy(m_nq, m_wsp[0], 1, outarray[0], 1);
272  }
273 
274  // Ion concentrations
275  for (unsigned int j = 0; j < m_concentrations.size(); ++j)
276  {
277  Vmath::Svtvp(m_nq, delta_t, m_wsp[m_concentrations[j]], 1, m_cellSol[m_concentrations[j]], 1, m_cellSol[m_concentrations[j]], 1);
278  }
279 
280  // Gating variables: Rush-Larsen scheme
281  for (unsigned int j = 0; j < m_gates.size(); ++j)
282  {
283  Vmath::Sdiv(m_nq, -delta_t, m_gates_tau[j], 1, m_gates_tau[j], 1);
284  Vmath::Vexp(m_nq, m_gates_tau[j], 1, m_gates_tau[j], 1);
285  Vmath::Vsub(m_nq, m_cellSol[m_gates[j]], 1, m_wsp[m_gates[j]], 1, m_cellSol[m_gates[j]], 1);
286  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);
287  }
288 
289  m_lastTime = time;
290  }
291 
292  Array<OneD, NekDouble> CellModel::GetCellSolutionCoeffs(unsigned int idx)
293  {
294  ASSERTL0(idx < m_nvar, "Index out of range for cell model.");
295 
296  Array<OneD, NekDouble> outarray(m_field->GetNcoeffs());
298 
299  if (m_useNodal)
300  {
301  for (unsigned int i = 0; i < m_field->GetNumElmts(); ++i)
302  {
303  int coef_offset = m_field->GetCoeff_Offset(i);
304  if (m_field->GetExp(0)->DetShapeType() == LibUtilities::eTriangle)
305  {
306  m_nodalTri->NodalToModal(m_cellSol[idx]+coef_offset, tmp=outarray+coef_offset);
307  }
308  else
309  {
310  m_nodalTet->NodalToModal(m_cellSol[idx]+coef_offset, tmp=outarray+coef_offset);
311  }
312  }
313  }
314  else
315  {
316  m_field->FwdTrans_IterPerExp(m_cellSol[idx], outarray);
317  }
318 
319  return outarray;
320  }
321 
322  Array<OneD, NekDouble> CellModel::GetCellSolution(unsigned int idx)
323  {
324  return m_cellSol[idx];
325  }
326 
327  void CellModel::LoadCellModel()
328  {
329  const bool root = (m_session->GetComm()->GetRank() == 0);
330  const std::string fncName = "CellModelInitialConditions";
331  const int nvar = m_cellSol[0].num_elements();
332  std::string varName;
333  Array<OneD, NekDouble> coeffs(m_field->GetNcoeffs());
335  int j = 0;
336 
337  SpatialDomains::MeshGraphSharedPtr vGraph = m_field->GetGraph();
338 
339  if (root)
340  {
341  cout << "Cell model initial conditions: " << endl;
342  }
343 
344  // First determine all the files we need to load
345  std::set<std::string> filelist;
346  for (j = 1; j < nvar; ++j)
347  {
348  // Get the name of the jth variable
349  varName = GetCellVarName(j);
350 
351  if (m_session->GetFunctionType(fncName, varName) ==
353  {
354  filelist.insert(m_session->GetFunctionFilename(fncName,
355  varName));
356  }
357  }
358 
359  // Read files
360  typedef std::vector<LibUtilities::FieldDefinitionsSharedPtr> FDef;
361  typedef std::vector<std::vector<NekDouble> > FData;
362  std::map<std::string, FDef> FieldDef;
363  std::map<std::string, FData> FieldData;
364  LibUtilities::FieldMetaDataMap fieldMetaDataMap;
366  std::set<std::string>::const_iterator setIt;
369  AllocateSharedPtr(m_session->GetComm());
370  for (setIt = filelist.begin(); setIt != filelist.end(); ++setIt)
371  {
372  if (root)
373  {
374  cout << " - Reading file: " << *setIt << endl;
375  }
376  FieldDef[*setIt] = FDef(0);
377  FieldData[*setIt] = FData(0);
378  fld->Import(*setIt, FieldDef[*setIt], FieldData[*setIt],
379  fieldMetaDataMap);
380  }
381 
382  // Get time of checkpoint from file if available
383  iter = fieldMetaDataMap.find("Time");
384  if(iter != fieldMetaDataMap.end())
385  {
386  m_lastTime = boost::lexical_cast<NekDouble>(iter->second);
387  }
388 
389  // Load each cell model variable
390  // j=0 and j=1 are for transmembrane or intra/extra-cellular volt.
391  Vmath::Zero(m_nq, m_cellSol[0], 1);
392  for(j = 1; j < m_cellSol.num_elements(); ++j)
393  {
394  // Get the name of the jth variable
395  varName = GetCellVarName(j);
396 
397  // Check if this variable is defined in a file or analytically
398  if (m_session->GetFunctionType(fncName, varName) ==
400  {
401  const std::string file =
402  m_session->GetFunctionFilename(fncName, varName);
403 
404  if (root)
405  {
406  cout << " - Field " << varName << ": from file "
407  << file << endl;
408  }
409 
410  // Extract the data into the modal coefficients
411  for(int i = 0; i < FieldDef[file].size(); ++i)
412  {
413  m_field->ExtractDataToCoeffs(FieldDef[file][i],
414  FieldData[file][i],
415  varName,
416  coeffs);
417  }
418 
419  // If using nodal cell model then we do a modal->nodal transform
420  // otherwise we do a backward transform onto physical points.
421  if (m_useNodal)
422  {
423  for (unsigned int i = 0; i < m_field->GetNumElmts(); ++i)
424  {
425  int coef_offset = m_field->GetCoeff_Offset(i);
426  if (m_field->GetExp(0)->DetShapeType() ==
428  {
429  m_nodalTri->ModalToNodal(coeffs+coef_offset,
430  tmp=m_cellSol[j]+coef_offset);
431  }
432  else
433  {
434  m_nodalTet->ModalToNodal(coeffs+coef_offset,
435  tmp=m_cellSol[j]+coef_offset);
436  }
437  }
438  }
439  else
440  {
441  m_field->BwdTrans(coeffs, m_cellSol[j]);
442  }
443  }
444  else if (m_session->GetFunctionType(fncName, varName) ==
446  {
448  m_session->GetFunction(fncName, varName);
449 
450  if (root)
451  {
452  cout << " - Field " << varName << ": "
453  << equ->GetExpression() << endl;
454  }
455 
456  const unsigned int nphys = m_field->GetNpoints();
457  Array<OneD, NekDouble> x0(nphys);
458  Array<OneD, NekDouble> x1(nphys);
459  Array<OneD, NekDouble> x2(nphys);
460  m_field->GetCoords(x0,x1,x2);
461 
462  if (m_useNodal)
463  {
464  Array<OneD, NekDouble> phys(nphys);
465  Array<OneD, NekDouble> tmpCoeffs(max(m_nodalTri->GetNcoeffs(), m_nodalTet->GetNcoeffs()));
466 
467  equ->Evaluate(x0, x1, x2, phys);
468  for (unsigned int i = 0; i < m_field->GetNumElmts(); ++i)
469  {
470  int phys_offset = m_field->GetPhys_Offset(i);
471  int coef_offset = m_field->GetCoeff_Offset(i);
472  if (m_field->GetExp(0)->DetShapeType() ==
474  {
475  m_field->GetExp(0)->FwdTrans(
476  phys + phys_offset, tmpCoeffs);
477  m_nodalTri->ModalToNodal(
478  tmpCoeffs,
479  tmp = m_cellSol[j] + coef_offset);
480  }
481  else
482  {
483  m_field->GetExp(0)->FwdTrans(
484  phys + phys_offset, tmpCoeffs);
485  m_nodalTet->ModalToNodal(
486  tmpCoeffs,
487  tmp = m_cellSol[j] + coef_offset);
488  }
489  }
490  }
491  else
492  {
493  equ->Evaluate(x0, x1, x2, m_cellSol[j]);
494  }
495  }
496  }
497  }
498 }
LibUtilities::NekFactory< std::string, CellModel, const LibUtilities::SessionReaderSharedPtr &, const MultiRegions::ExpListSharedPtr & > CellModelFactory
Datatype of the NekFactory used to instantiate classes derived from the EquationSystem class...
Definition: CellModel.h:61
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Principle Modified Functions .
Definition: BasisType.h:51
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
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:471
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:428
Principle Modified Functions .
Definition: BasisType.h:49
STL namespace.
void Sdiv(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha/y.
Definition: Vmath.cpp:257
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
std::map< std::string, std::string > FieldMetaDataMap
Definition: FieldIO.h:53
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:57
Principle Modified Functions .
Definition: BasisType.h:50
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
void Vexp(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:107
Defines a specification for a set of points.
Definition: Points.h:58
boost::shared_ptr< FieldIO > FieldIOSharedPtr
Definition: FieldIO.h:225
double NekDouble
3D Evenly-spaced points on a Tetrahedron
Definition: PointsType.h:71
boost::shared_ptr< Equation > EquationSharedPtr
CellModelFactory & GetCellModelFactory()
Definition: CellModel.cpp:47
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:329
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:58
2D Evenly-spaced points on a Triangle
Definition: PointsType.h:70
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:442
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
Describes the specification for a Basis.
Definition: Basis.h:50
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:50
Provides a generic Factory class.
Definition: NekFactory.hpp:116