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