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 Reference

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.
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.
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.
void GenerateSummary (SummaryList &s)
 Print a summary of the cell model.
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.
MultiRegions::ExpListSharedPtr m_field
 Transmembrane potential field from PDE system.
int m_nq
 Number of physical points.
int m_nvar
 Number of variables in cell model (inc. transmembrane voltage)
NekDouble m_lastTime
 Timestep for pde model.
int m_substeps
 Number of substeps to take.
Array< OneD, Array< OneD,
NekDouble > > 
m_cellSol
 Cell model solution variables.
Array< OneD, Array< OneD,
NekDouble > > 
m_wsp
 Cell model integration workspace.
bool m_useNodal
 Flag indicating whether nodal projection in use.
StdRegions::StdNodalTriExpSharedPtr m_nodalTri
 StdNodalTri for cell model calculations.
StdRegions::StdNodalTetExpSharedPtr m_nodalTet
Array< OneD, Array< OneD,
NekDouble > > 
m_nodalTmp
 Temporary array for nodal projection.
std::vector< int > m_concentrations
 Indices of cell model variables which are concentrations.
std::vector< int > m_gates
 Indices of cell model variables which are gates.
Array< OneD, Array< OneD,
NekDouble > > 
m_gates_tau
 Storage for gate tau values.

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::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.

{
m_session = pSession;
m_field = pField;
m_lastTime = 0.0;
m_substeps = pSession->GetParameter("Substeps");
m_nvar = 0;
m_useNodal = false;
// Number of points in nodal space is the number of coefficients
// in modified basis
std::set<enum LibUtilities::ShapeType> s;
for (unsigned int i = 0; i < m_field->GetNumElmts(); ++i)
{
s.insert(m_field->GetExp(i)->DetShapeType());
}
// Use nodal projection if only triangles
if (s.size() == 1 && (s.count(LibUtilities::eTriangle) == 1 ||
{
// This is disabled for now as it causes problems at high order.
// m_useNodal = true;
}
// ---------------------------
// Move to nodal points
{
m_nq = pField->GetNcoeffs();
int order = m_field->GetExp(0)->GetBasis(0)->GetNumModes();
// Set up a nodal tri
}
else
{
m_nq = pField->GetTotPoints();
}
}
virtual Nektar::CellModel::~CellModel ( )
inlinevirtual

Definition at line 71 of file CellModel.h.

{}

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().

Array< OneD, NekDouble > Nektar::CellModel::GetCellSolution ( unsigned int  idx)

Definition at line 320 of file CellModel.cpp.

References m_cellSol.

{
return m_cellSol[idx];
}
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.

{
ASSERTL0(idx < m_nvar, "Index out of range for cell model.");
Array<OneD, NekDouble> outarray(m_field->GetNcoeffs());
Array<OneD, NekDouble> tmp;
{
for (unsigned int i = 0; i < m_field->GetNumElmts(); ++i)
{
int coef_offset = m_field->GetCoeff_Offset(i);
if (m_field->GetExp(0)->DetShapeType() == LibUtilities::eTriangle)
{
m_nodalTri->NodalToModal(m_cellSol[idx]+coef_offset, tmp=outarray+coef_offset);
}
else
{
m_nodalTet->NodalToModal(m_cellSol[idx]+coef_offset, tmp=outarray+coef_offset);
}
}
}
else
{
m_field->FwdTrans_IterPerExp(m_cellSol[idx], outarray);
}
return outarray;
}
std::string Nektar::CellModel::GetCellVarName ( unsigned int  idx)
inline

Definition at line 102 of file CellModel.h.

References v_GetCellVarName().

Referenced by LoadCellModel().

{
return v_GetCellVarName(idx);
}
unsigned int Nektar::CellModel::GetNumCellVariables ( )
inline

Definition at line 97 of file CellModel.h.

References m_nvar.

{
return m_nvar;
}
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().

{
ASSERTL1(m_nvar > 0, "Cell model must have at least 1 variable.");
m_cellSol = Array<OneD, Array<OneD, NekDouble> >(m_nvar);
m_wsp = Array<OneD, Array<OneD, NekDouble> >(m_nvar);
for (unsigned int i = 0; i < m_nvar; ++i)
{
m_cellSol[i] = Array<OneD, NekDouble>(m_nq);
m_wsp[i] = Array<OneD, NekDouble>(m_nq);
}
m_gates_tau = Array<OneD, Array<OneD, NekDouble> >(m_gates.size());
for (unsigned int i = 0; i < m_gates.size(); ++i)
{
m_gates_tau[i] = Array<OneD, NekDouble>(m_nq);
}
if (m_session->DefinesFunction("CellModelInitialConditions"))
{
}
else
{
}
}
void Nektar::CellModel::LoadCellModel ( )
protected

Definition at line 325 of file CellModel.cpp.

References 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().

{
const bool root = (m_session->GetComm()->GetRank() == 0);
const std::string fncName = "CellModelInitialConditions";
std::string varName;
Array<OneD, NekDouble> coeffs(m_field->GetNcoeffs());
Array<OneD, NekDouble> tmp;
if (root)
{
cout << "Cell model initial conditions: " << endl;
}
// Load each cell model variable
// j=0 and j=1 are for transmembrane or intra/extra-cellular volt.
for(int j = 1; j < m_cellSol.num_elements(); ++j)
{
// Get the name of the jth variable
varName = GetCellVarName(j);
// Check if this variable is defined in a file or analytically
if (m_session->GetFunctionType(fncName, varName) ==
{
const std::string file =
m_session->GetFunctionFilename(fncName, varName);
if (root)
{
cout << " - Field " << varName << ": from file "
<< file << endl;
}
std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
std::vector<std::vector<NekDouble> > FieldData;
// Read the restart file containing this variable
fld->Import(file, FieldDef, FieldData);
fld->ImportFieldMetaData(file,fieldMetaDataMap);
iter = fieldMetaDataMap.find("Time");
if(iter != fieldMetaDataMap.end())
{
m_lastTime = boost::lexical_cast<NekDouble>(iter->second);
}
// Extract the data into the modal coefficients
for(int i = 0; i < FieldDef.size(); ++i)
{
m_field->ExtractDataToCoeffs(FieldDef[i],
FieldData[i],
varName,
coeffs);
}
// If using nodal cell model then we do a modal->nodal transform
// otherwise we do a backward transform onto physical points.
{
for (unsigned int i = 0; i < m_field->GetNumElmts(); ++i)
{
int coef_offset = m_field->GetCoeff_Offset(i);
if (m_field->GetExp(0)->DetShapeType() ==
{
m_nodalTri->ModalToNodal(coeffs+coef_offset,
tmp=m_cellSol[j]+coef_offset);
}
else
{
m_nodalTet->ModalToNodal(coeffs+coef_offset,
tmp=m_cellSol[j]+coef_offset);
}
}
}
else
{
m_field->BwdTrans(coeffs, m_cellSol[j]);
}
}
else if (m_session->GetFunctionType(fncName, varName) ==
{
m_session->GetFunction(fncName, varName);
if (root)
{
cout << " - Field " << varName << ": "
<< equ->GetExpression() << endl;
}
const unsigned int nphys = m_field->GetNpoints();
Array<OneD, NekDouble> x0(nphys);
Array<OneD, NekDouble> x1(nphys);
Array<OneD, NekDouble> x2(nphys);
m_field->GetCoords(x0,x1,x2);
{
Array<OneD, NekDouble> phys(nphys);
Array<OneD, NekDouble> tmpCoeffs(max(m_nodalTri->GetNcoeffs(), m_nodalTet->GetNcoeffs()));
equ->Evaluate(x0, x1, x2, phys);
for (unsigned int i = 0; i < m_field->GetNumElmts(); ++i)
{
int phys_offset = m_field->GetPhys_Offset(i);
int coef_offset = m_field->GetCoeff_Offset(i);
if (m_field->GetExp(0)->DetShapeType() ==
{
m_field->GetExp(0)->FwdTrans(
phys + phys_offset, tmpCoeffs);
m_nodalTri->ModalToNodal(
tmpCoeffs,
tmp = m_cellSol[j] + coef_offset);
}
else
{
m_field->GetExp(0)->FwdTrans(
phys + phys_offset, tmpCoeffs);
m_nodalTet->ModalToNodal(
tmpCoeffs,
tmp = m_cellSol[j] + coef_offset);
}
}
}
else
{
equ->Evaluate(x0, x1, x2, m_cellSol[j]);
}
}
}
}
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().

{
int phys_offset = 0;
int coef_offset = 0;
int nvar = inarray.num_elements();
Array<OneD, NekDouble> tmp;
// ---------------------------
// Check nodal temp array set up
{
if (!m_nodalTmp.num_elements())
{
m_nodalTmp = Array<OneD, Array<OneD, NekDouble> >(nvar);
for (unsigned int k = 0; k < nvar; ++k)
{
m_nodalTmp[k] = Array<OneD, NekDouble>(m_nq);
}
}
// Move to nodal points
Array<OneD, NekDouble> tmpCoeffs(max(m_nodalTri->GetNcoeffs(), m_nodalTet->GetNcoeffs()));
for (unsigned int k = 0; k < nvar; ++k)
{
for (unsigned int i = 0; i < m_field->GetNumElmts(); ++i)
{
phys_offset = m_field->GetPhys_Offset(i);
coef_offset = m_field->GetCoeff_Offset(i);
if (m_field->GetExp(0)->DetShapeType() == LibUtilities::eTriangle)
{
m_field->GetExp(0)->FwdTrans(inarray[k] + phys_offset, tmpCoeffs);
m_nodalTri->ModalToNodal(tmpCoeffs, tmp=m_nodalTmp[k]+coef_offset);
}
else
{
m_field->GetExp(0)->FwdTrans(inarray[k] + phys_offset, tmpCoeffs);
m_nodalTet->ModalToNodal(tmpCoeffs, tmp=m_nodalTmp[k]+coef_offset);
}
}
}
// Copy new transmembrane potential into cell model
}
else
{
// Copy new transmembrane potential into cell model
Vmath::Vcopy(m_nq, inarray[0], 1, m_cellSol[0], 1);
}
// -------------------------
NekDouble delta_t = (time - m_lastTime)/m_substeps;
// Perform substepping
for (unsigned int i = 0; i < m_substeps - 1; ++i)
{
// Voltage
Vmath::Svtvp(m_nq, delta_t, m_wsp[0], 1, m_cellSol[0], 1, m_cellSol[0], 1);
// Ion concentrations
for (unsigned int j = 0; j < m_concentrations.size(); ++j)
{
}
// Gating variables: Rush-Larsen scheme
for (unsigned int j = 0; j < m_gates.size(); ++j)
{
Vmath::Sdiv(m_nq, -delta_t, m_gates_tau[j], 1, m_gates_tau[j], 1);
}
}
// Perform final cell model step
// Output dV/dt from last step but integrate remaining cell model vars
// Transform cell model I_total from nodal to modal space
{
Array<OneD, NekDouble> tmpCoeffs(max(m_nodalTri->GetNcoeffs(), m_nodalTet->GetNcoeffs()));
for (unsigned int k = 0; k < nvar; ++k)
{
for (unsigned int i = 0; i < m_field->GetNumElmts(); ++i)
{
int phys_offset = m_field->GetPhys_Offset(i);
int coef_offset = m_field->GetCoeff_Offset(i);
if (m_field->GetExp(0)->DetShapeType() == LibUtilities::eTriangle)
{
m_nodalTri->NodalToModal(m_wsp[k]+coef_offset, tmpCoeffs);
m_field->GetExp(0)->BwdTrans(tmpCoeffs, tmp=outarray[k] + phys_offset);
}
else
{
m_nodalTet->NodalToModal(m_wsp[k]+coef_offset, tmpCoeffs);
m_field->GetExp(0)->BwdTrans(tmpCoeffs, tmp=outarray[k] + phys_offset);
}
}
}
}
else
{
Vmath::Vcopy(m_nq, m_wsp[0], 1, outarray[0], 1);
}
// Ion concentrations
for (unsigned int j = 0; j < m_concentrations.size(); ++j)
{
}
// Gating variables: Rush-Larsen scheme
for (unsigned int j = 0; j < m_gates.size(); ++j)
{
Vmath::Sdiv(m_nq, -delta_t, m_gates_tau[j], 1, m_gates_tau[j], 1);
}
m_lastTime = time;
}
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().

{
v_Update(inarray, outarray, time);
}
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().

{
return "Var" + boost::lexical_cast<std::string>(idx);
}
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().