Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | Protected Member Functions | Protected Attributes | List of all members
Nektar::CellModel Class Referenceabstract

Cell model base class. More...

#include <CellModel.h>

Inheritance diagram for Nektar::CellModel:
Inheritance graph
[legend]
Collaboration diagram for Nektar::CellModel:
Collaboration graph
[legend]

Public Member Functions

 CellModel (const LibUtilities::SessionReaderSharedPtr &pSession, const MultiRegions::ExpListSharedPtr &pField)
 
virtual ~CellModel ()
 
void Initialise ()
 Initialise the cell model storage and set initial conditions. More...
 
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. More...
 
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. More...
 
void GenerateSummary (SummaryList &s)
 Print a summary of the cell model. More...
 
unsigned int GetNumCellVariables ()
 
std::string GetCellVarName (unsigned int idx)
 
Array< OneD, NekDoubleGetCellSolutionCoeffs (unsigned int idx)
 
Array< OneD, NekDoubleGetCellSolution (unsigned int idx)
 

Protected Member Functions

virtual void v_Update (const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)=0
 
virtual void v_GenerateSummary (SummaryList &s)=0
 
virtual std::string v_GetCellVarName (unsigned int idx)
 
virtual void v_SetInitialConditions ()=0
 
void LoadCellModel ()
 

Protected Attributes

LibUtilities::SessionReaderSharedPtr m_session
 Session. More...
 
MultiRegions::ExpListSharedPtr m_field
 Transmembrane potential field from PDE system. More...
 
int m_nq
 Number of physical points. More...
 
int m_nvar
 Number of variables in cell model (inc. transmembrane voltage) More...
 
NekDouble m_lastTime
 Timestep for pde model. More...
 
int m_substeps
 Number of substeps to take. More...
 
Array< OneD, Array< OneD,
NekDouble > > 
m_cellSol
 Cell model solution variables. More...
 
Array< OneD, Array< OneD,
NekDouble > > 
m_wsp
 Cell model integration workspace. More...
 
bool m_useNodal
 Flag indicating whether nodal projection in use. More...
 
StdRegions::StdNodalTriExpSharedPtr m_nodalTri
 StdNodalTri for cell model calculations. More...
 
StdRegions::StdNodalTetExpSharedPtr m_nodalTet
 
Array< OneD, Array< OneD,
NekDouble > > 
m_nodalTmp
 Temporary array for nodal projection. More...
 
std::vector< int > m_concentrations
 Indices of cell model variables which are concentrations. More...
 
std::vector< int > m_gates
 Indices of cell model variables which are gates. More...
 
Array< OneD, Array< OneD,
NekDouble > > 
m_gates_tau
 Storage for gate tau values. More...
 

Detailed Description

Cell model base class.

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.

Definition at line 65 of file CellModel.h.

Constructor & Destructor Documentation

Nektar::CellModel::CellModel ( const LibUtilities::SessionReaderSharedPtr pSession,
const MultiRegions::ExpListSharedPtr pField 
)

Cell model base class constructor.

Definition at line 69 of file CellModel.cpp.

References Nektar::LibUtilities::eGaussLobattoLegendre, Nektar::LibUtilities::eGaussRadauMAlpha1Beta0, Nektar::LibUtilities::eGaussRadauMAlpha2Beta0, Nektar::LibUtilities::eModified_A, Nektar::LibUtilities::eModified_B, Nektar::LibUtilities::eModified_C, Nektar::LibUtilities::eNodalTetEvenlySpaced, Nektar::LibUtilities::eNodalTriEvenlySpaced, Nektar::LibUtilities::eTetrahedron, and Nektar::LibUtilities::eTriangle.

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
103  LibUtilities::BasisKey B0(
105  LibUtilities::PointsKey(order, LibUtilities::eGaussLobattoLegendre));
106  LibUtilities::BasisKey B1(
108  LibUtilities::PointsKey(order, LibUtilities::eGaussRadauMAlpha1Beta0));
109  LibUtilities::BasisKey B2(
111  LibUtilities::PointsKey(order, LibUtilities::eGaussRadauMAlpha2Beta0));
112 
117  }
118  else
119  {
120  m_nq = pField->GetTotPoints();
121  }
122  }
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.
Principle Modified Functions .
Definition: BasisType.h:49
StdRegions::StdNodalTriExpSharedPtr m_nodalTri
StdNodalTri for cell model calculations.
Definition: CellModel.h:133
MultiRegions::ExpListSharedPtr m_field
Transmembrane potential field from PDE system.
Definition: CellModel.h:115
StdRegions::StdNodalTetExpSharedPtr m_nodalTet
Definition: CellModel.h:134
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:57
NekDouble m_lastTime
Timestep for pde model.
Definition: CellModel.h:121
Principle Modified Functions .
Definition: BasisType.h:50
int m_nvar
Number of variables in cell model (inc. transmembrane voltage)
Definition: CellModel.h:119
3D Evenly-spaced points on a Tetrahedron
Definition: PointsType.h:71
bool m_useNodal
Flag indicating whether nodal projection in use.
Definition: CellModel.h:131
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:58
2D Evenly-spaced points on a Triangle
Definition: PointsType.h:70
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
virtual Nektar::CellModel::~CellModel ( )
inlinevirtual

Definition at line 71 of file CellModel.h.

71 {}

Member Function Documentation

void Nektar::CellModel::GenerateSummary ( SummaryList s)
inline

Print a summary of the cell model.

Definition at line 92 of file CellModel.h.

References v_GenerateSummary().

93  {
95  }
virtual void v_GenerateSummary(SummaryList &s)=0
Array< OneD, NekDouble > Nektar::CellModel::GetCellSolution ( unsigned int  idx)

Definition at line 322 of file CellModel.cpp.

323  {
324  return m_cellSol[idx];
325  }
Array< OneD, Array< OneD, NekDouble > > m_cellSol
Cell model solution variables.
Definition: CellModel.h:126
Array< OneD, NekDouble > Nektar::CellModel::GetCellSolutionCoeffs ( unsigned int  idx)

Definition at line 292 of file CellModel.cpp.

References ASSERTL0, and Nektar::LibUtilities::eTriangle.

293  {
294  ASSERTL0(idx < m_nvar, "Index out of range for cell model.");
295 
296  Array<OneD, NekDouble> outarray(m_field->GetNcoeffs());
297  Array<OneD, NekDouble> tmp;
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  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
StdRegions::StdNodalTriExpSharedPtr m_nodalTri
StdNodalTri for cell model calculations.
Definition: CellModel.h:133
MultiRegions::ExpListSharedPtr m_field
Transmembrane potential field from PDE system.
Definition: CellModel.h:115
StdRegions::StdNodalTetExpSharedPtr m_nodalTet
Definition: CellModel.h:134
Array< OneD, Array< OneD, NekDouble > > m_cellSol
Cell model solution variables.
Definition: CellModel.h:126
int m_nvar
Number of variables in cell model (inc. transmembrane voltage)
Definition: CellModel.h:119
bool m_useNodal
Flag indicating whether nodal projection in use.
Definition: CellModel.h:131
std::string Nektar::CellModel::GetCellVarName ( unsigned int  idx)
inline

Definition at line 102 of file CellModel.h.

References v_GetCellVarName().

103  {
104  return v_GetCellVarName(idx);
105  }
virtual std::string v_GetCellVarName(unsigned int idx)
Definition: CellModel.h:152
unsigned int Nektar::CellModel::GetNumCellVariables ( )
inline

Definition at line 97 of file CellModel.h.

References m_nvar.

98  {
99  return m_nvar;
100  }
int m_nvar
Number of variables in cell model (inc. transmembrane voltage)
Definition: CellModel.h:119
void Nektar::CellModel::Initialise ( )

Initialise the cell model storage and set initial conditions.

Initialise the cell model. Allocate workspace and variable storage.

Definition at line 128 of file CellModel.cpp.

References ASSERTL1.

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  {
152  }
153  }
int m_nq
Number of physical points.
Definition: CellModel.h:117
Array< OneD, Array< OneD, NekDouble > > m_wsp
Cell model integration workspace.
Definition: CellModel.h:128
virtual void v_SetInitialConditions()=0
Array< OneD, Array< OneD, NekDouble > > m_gates_tau
Storage for gate tau values.
Definition: CellModel.h:143
Array< OneD, Array< OneD, NekDouble > > m_cellSol
Cell model solution variables.
Definition: CellModel.h:126
int m_nvar
Number of variables in cell model (inc. transmembrane voltage)
Definition: CellModel.h:119
std::vector< int > m_gates
Indices of cell model variables which are gates.
Definition: CellModel.h:141
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:218
LibUtilities::SessionReaderSharedPtr m_session
Session.
Definition: CellModel.h:113
void Nektar::CellModel::LoadCellModel ( )
protected

Definition at line 327 of file CellModel.cpp.

References Nektar::LibUtilities::eFunctionTypeExpression, Nektar::LibUtilities::eFunctionTypeFile, Nektar::LibUtilities::eTriangle, Nektar::iterator, and Vmath::Zero().

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());
334  Array<OneD, NekDouble> tmp;
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  }
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
StdRegions::StdNodalTriExpSharedPtr m_nodalTri
StdNodalTri for cell model calculations.
Definition: CellModel.h:133
MultiRegions::ExpListSharedPtr m_field
Transmembrane potential field from PDE system.
Definition: CellModel.h:115
std::map< std::string, std::string > FieldMetaDataMap
Definition: FieldIO.h:53
StdRegions::StdNodalTetExpSharedPtr m_nodalTet
Definition: CellModel.h:134
NekDouble m_lastTime
Timestep for pde model.
Definition: CellModel.h:121
Array< OneD, Array< OneD, NekDouble > > m_cellSol
Cell model solution variables.
Definition: CellModel.h:126
boost::shared_ptr< FieldIO > FieldIOSharedPtr
Definition: FieldIO.h:225
double NekDouble
boost::shared_ptr< Equation > EquationSharedPtr
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
bool m_useNodal
Flag indicating whether nodal projection in use.
Definition: CellModel.h:131
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:442
LibUtilities::SessionReaderSharedPtr m_session
Session.
Definition: CellModel.h:113
void Nektar::CellModel::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.

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.

Definition at line 163 of file CellModel.cpp.

References Nektar::LibUtilities::eTriangle, Vmath::Sdiv(), Vmath::Svtvp(), Vmath::Vcopy(), Vmath::Vexp(), Vmath::Vsub(), and Vmath::Vvtvp().

167  {
168  int phys_offset = 0;
169  int coef_offset = 0;
170  int nvar = inarray.num_elements();
171  Array<OneD, NekDouble> tmp;
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  {
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  {
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  }
int m_nq
Number of physical points.
Definition: CellModel.h:117
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
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
MultiRegions::ExpListSharedPtr m_field
Transmembrane potential field from PDE system.
Definition: CellModel.h:115
StdRegions::StdNodalTetExpSharedPtr m_nodalTet
Definition: CellModel.h:134
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
std::vector< int > m_concentrations
Indices of cell model variables which are concentrations.
Definition: CellModel.h:139
Array< OneD, Array< OneD, NekDouble > > m_gates_tau
Storage for gate tau values.
Definition: CellModel.h:143
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
std::vector< int > m_gates
Indices of cell model variables which are gates.
Definition: CellModel.h:141
double NekDouble
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
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
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
int m_substeps
Number of substeps to take.
Definition: CellModel.h:123
void Nektar::CellModel::Update ( const Array< OneD, const Array< OneD, NekDouble > > &  inarray,
Array< OneD, Array< OneD, NekDouble > > &  outarray,
const NekDouble  time 
)
inline

Compute the derivatives of cell model variables.

Definition at line 83 of file CellModel.h.

References v_Update().

87  {
88  v_Update(inarray, outarray, time);
89  }
virtual void v_Update(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)=0
virtual void Nektar::CellModel::v_GenerateSummary ( SummaryList s)
protectedpure virtual
virtual std::string Nektar::CellModel::v_GetCellVarName ( unsigned int  idx)
inlineprotectedvirtual

Reimplemented in Nektar::CourtemancheRamirezNattel98, and Nektar::FentonKarma.

Definition at line 152 of file CellModel.h.

Referenced by GetCellVarName().

153  {
154  return "Var" + boost::lexical_cast<std::string>(idx);
155  }
virtual void Nektar::CellModel::v_SetInitialConditions ( )
protectedpure virtual
virtual void Nektar::CellModel::v_Update ( const Array< OneD, const Array< OneD, NekDouble > > &  inarray,
Array< OneD, Array< OneD, NekDouble > > &  outarray,
const NekDouble  time 
)
protectedpure virtual

Member Data Documentation

Array<OneD, Array<OneD, NekDouble> > Nektar::CellModel::m_cellSol
protected
std::vector<int> Nektar::CellModel::m_concentrations
protected
MultiRegions::ExpListSharedPtr Nektar::CellModel::m_field
protected

Transmembrane potential field from PDE system.

Definition at line 115 of file CellModel.h.

std::vector<int> Nektar::CellModel::m_gates
protected
Array<OneD, Array<OneD, NekDouble> > Nektar::CellModel::m_gates_tau
protected
NekDouble Nektar::CellModel::m_lastTime
protected

Timestep for pde model.

Definition at line 121 of file CellModel.h.

StdRegions::StdNodalTetExpSharedPtr Nektar::CellModel::m_nodalTet
protected

Definition at line 134 of file CellModel.h.

Array<OneD, Array<OneD, NekDouble> > Nektar::CellModel::m_nodalTmp
protected

Temporary array for nodal projection.

Definition at line 136 of file CellModel.h.

StdRegions::StdNodalTriExpSharedPtr Nektar::CellModel::m_nodalTri
protected

StdNodalTri for cell model calculations.

Definition at line 133 of file CellModel.h.

int Nektar::CellModel::m_nq
protected
int Nektar::CellModel::m_nvar
protected
LibUtilities::SessionReaderSharedPtr Nektar::CellModel::m_session
protected

Session.

Definition at line 113 of file CellModel.h.

int Nektar::CellModel::m_substeps
protected

Number of substeps to take.

Definition at line 123 of file CellModel.h.

bool Nektar::CellModel::m_useNodal
protected

Flag indicating whether nodal projection in use.

Definition at line 131 of file CellModel.h.

Array<OneD, Array<OneD, NekDouble> > Nektar::CellModel::m_wsp
protected

Cell model integration workspace.

Definition at line 128 of file CellModel.h.