Nektar++
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 67 of file CellModel.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), 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, Nektar::LibUtilities::eTriangle, m_field, m_lastTime, m_nodalTet, m_nodalTri, m_nq, m_nvar, m_session, m_substeps, and m_useNodal.

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  }
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
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
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
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
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 320 of file CellModel.cpp.

References m_cellSol.

321  {
322  return m_cellSol[idx];
323  }
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 290 of file CellModel.cpp.

References ASSERTL0, Nektar::LibUtilities::eTriangle, m_cellSol, m_field, m_nodalTet, m_nodalTri, m_nvar, and m_useNodal.

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  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
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().

Referenced by LoadCellModel().

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 126 of file CellModel.cpp.

References ASSERTL1, LoadCellModel(), m_cellSol, m_gates, m_gates_tau, m_nq, m_nvar, m_session, m_wsp, and v_SetInitialConditions().

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  }
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:191
LibUtilities::SessionReaderSharedPtr m_session
Session.
Definition: CellModel.h:113
void Nektar::CellModel::LoadCellModel ( )
protected

Definition at line 325 of file CellModel.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::LibUtilities::eFunctionTypeExpression, Nektar::LibUtilities::eFunctionTypeFile, Nektar::LibUtilities::eTriangle, GetCellVarName(), Nektar::iterator, m_cellSol, m_field, m_lastTime, m_nodalTet, m_nodalTri, m_nq, m_session, m_useNodal, and Vmath::Zero().

Referenced by Initialise().

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  }
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:65
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:236
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:432
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 161 of file CellModel.cpp.

References Nektar::LibUtilities::eTriangle, m_cellSol, m_concentrations, m_field, m_gates, m_gates_tau, m_lastTime, m_nodalTet, m_nodalTmp, m_nodalTri, m_nq, m_substeps, m_useNodal, m_wsp, Vmath::Sdiv(), Vmath::Svtvp(), Update(), Vmath::Vcopy(), Vmath::Vexp(), Vmath::Vsub(), and Vmath::Vvtvp().

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  }
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:1038
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().

Referenced by TimeIntegrate().

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.

Referenced by CellModel(), GetCellSolutionCoeffs(), LoadCellModel(), and TimeIntegrate().

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.

Referenced by CellModel(), LoadCellModel(), and TimeIntegrate().

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

Definition at line 134 of file CellModel.h.

Referenced by CellModel(), GetCellSolutionCoeffs(), LoadCellModel(), and TimeIntegrate().

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

Temporary array for nodal projection.

Definition at line 136 of file CellModel.h.

Referenced by TimeIntegrate().

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

StdNodalTri for cell model calculations.

Definition at line 133 of file CellModel.h.

Referenced by CellModel(), GetCellSolutionCoeffs(), LoadCellModel(), and TimeIntegrate().

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.

Referenced by CellModel(), Initialise(), and LoadCellModel().

int Nektar::CellModel::m_substeps
protected

Number of substeps to take.

Definition at line 123 of file CellModel.h.

Referenced by CellModel(), and TimeIntegrate().

bool Nektar::CellModel::m_useNodal
protected

Flag indicating whether nodal projection in use.

Definition at line 131 of file CellModel.h.

Referenced by CellModel(), GetCellSolutionCoeffs(), LoadCellModel(), and TimeIntegrate().

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

Cell model integration workspace.

Definition at line 128 of file CellModel.h.

Referenced by Initialise(), and TimeIntegrate().