Cell model base class. More...

#include <CellModel.h>

Inheritance diagram for Nektar::CellModel:

Collaboration diagram for Nektar::CellModel:

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, NekDouble >	GetCellSolutionCoeffs (unsigned int idx)

Array< OneD, NekDouble >	GetCellSolution (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.

     {
         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 || 
                               s.count(LibUtilities::eTetrahedron) == 1))
         {
             // This is disabled for now as it causes problems at high order.
             // m_useNodal = true;
         }
 
         // ---------------------------
         // Move to nodal points
         if (m_useNodal)
         {
             m_nq = pField->GetNcoeffs();
             int order = m_field->GetExp(0)->GetBasis(0)->GetNumModes();
 
             // Set up a nodal tri
             LibUtilities::BasisKey B0(
                 LibUtilities::eModified_A, order,
                 LibUtilities::PointsKey(order, LibUtilities::eGaussLobattoLegendre));
             LibUtilities::BasisKey B1(
                 LibUtilities::eModified_B, order,
                 LibUtilities::PointsKey(order, LibUtilities::eGaussRadauMAlpha1Beta0));
             LibUtilities::BasisKey B2(
                 LibUtilities::eModified_C, order,
                 LibUtilities::PointsKey(order, LibUtilities::eGaussRadauMAlpha2Beta0));
 
             m_nodalTri = MemoryManager<StdRegions::StdNodalTriExp>::AllocateSharedPtr(
                         B0, B1, LibUtilities::eNodalTriEvenlySpaced);
             m_nodalTet = MemoryManager<StdRegions::StdNodalTetExp>::AllocateSharedPtr(
                         B0, B1, B2, LibUtilities::eNodalTetEvenlySpaced);
         }
         else
         {
             m_nq = pField->GetTotPoints();
         }
     }

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

         {
             v_GenerateSummary(s);
         }

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;
 
         if (m_useNodal)
         {
             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"))
         {
             LoadCellModel();
         }
         else
         {
             v_SetInitialConditions();
         }
     }

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

     {
         const bool root = (m_session->GetComm()->GetRank() == 0);
         const std::string fncName = "CellModelInitialConditions";
         const int nvar = m_cellSol[0].num_elements();
         std::string varName;
         Array<OneD, NekDouble> coeffs(m_field->GetNcoeffs());
         Array<OneD, NekDouble> tmp;
         int j = 0;
 
         SpatialDomains::MeshGraphSharedPtr vGraph = m_field->GetGraph();
 
         if (root)
         {
             cout << "Cell model initial conditions: " << endl;
         }
 
         // First determine all the files we need to load
         std::set<std::string> filelist;
         for (j = 1; j < nvar; ++j)
         {
             // Get the name of the jth variable
             varName = GetCellVarName(j);
 
             if (m_session->GetFunctionType(fncName, varName) ==
                     LibUtilities::eFunctionTypeFile)
             {
                 filelist.insert(m_session->GetFunctionFilename(fncName,
                                                                varName));
             }
         }
 
         // Read files
         typedef std::vector<LibUtilities::FieldDefinitionsSharedPtr> FDef;
         typedef std::vector<std::vector<NekDouble> > FData;
         std::map<std::string, FDef>  FieldDef;
         std::map<std::string, FData> FieldData;
         LibUtilities::FieldMetaDataMap fieldMetaDataMap;
         LibUtilities::FieldMetaDataMap::iterator iter;
         std::set<std::string>::const_iterator setIt;
         LibUtilities::FieldIOSharedPtr fld =
                 MemoryManager<LibUtilities::FieldIO>::
                                 AllocateSharedPtr(m_session->GetComm());
         for (setIt = filelist.begin(); setIt != filelist.end(); ++setIt)
         {
             if (root)
             {
                 cout << "  - Reading file: " << *setIt << endl;
             }
             FieldDef[*setIt] = FDef(0);
             FieldData[*setIt] = FData(0);
             fld->Import(*setIt, FieldDef[*setIt], FieldData[*setIt],
                         fieldMetaDataMap);
         }
 
         // Get time of checkpoint from file if available
         iter = fieldMetaDataMap.find("Time");
         if(iter != fieldMetaDataMap.end())
         {
              m_lastTime = boost::lexical_cast<NekDouble>(iter->second);
         }
 
         // Load each cell model variable
         // j=0 and j=1 are for transmembrane or intra/extra-cellular volt.
         Vmath::Zero(m_nq, m_cellSol[0], 1);
         for(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) ==
                     LibUtilities::eFunctionTypeFile)
             {
                 const std::string file =
                         m_session->GetFunctionFilename(fncName, varName);
 
                 if (root)
                 {
                     cout << "  - Field " << varName << ": from file "
                          << file << endl;
                 }
 
                 // Extract the data into the modal coefficients
                 for(int i = 0; i < FieldDef[file].size(); ++i)
                 {
                     m_field->ExtractDataToCoeffs(FieldDef[file][i],
                                                  FieldData[file][i],
                                                  varName,
                                                  coeffs);
                 }
 
                 // If using nodal cell model then we do a modal->nodal transform
                 // otherwise we do a backward transform onto physical points.
                 if (m_useNodal)
                 {
                     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->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) ==
                     LibUtilities::eFunctionTypeExpression)
             {
                 LibUtilities::EquationSharedPtr equ =
                         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);
 
                 if (m_useNodal)
                 {
                     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() ==
                                 LibUtilities::eTriangle)
                         {
                             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_useNodal)
         {
             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
             Vmath::Vcopy(m_nq, m_nodalTmp[0], 1, m_cellSol[0], 1);
         }
         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)
         {
             Update(m_cellSol, m_wsp, time);
             // 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)
             {
                 Vmath::Svtvp(m_nq, delta_t, m_wsp[m_concentrations[j]], 1, m_cellSol[m_concentrations[j]], 1, m_cellSol[m_concentrations[j]], 1);
             }
             // 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);
                 Vmath::Vexp(m_nq, m_gates_tau[j], 1, m_gates_tau[j], 1);
                 Vmath::Vsub(m_nq, m_cellSol[m_gates[j]], 1, m_wsp[m_gates[j]], 1, m_cellSol[m_gates[j]], 1);
                 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);
             }
         }
 
         // Perform final cell model step
         Update(m_cellSol, m_wsp, time);
 
         // Output dV/dt from last step but integrate remaining cell model vars
         // Transform cell model I_total from nodal to modal space
         if (m_useNodal)
         {
             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)
         {
             Vmath::Svtvp(m_nq, delta_t, m_wsp[m_concentrations[j]], 1, m_cellSol[m_concentrations[j]], 1, m_cellSol[m_concentrations[j]], 1);
         }
 
         // 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);
             Vmath::Vexp(m_nq, m_gates_tau[j], 1, m_gates_tau[j], 1);
             Vmath::Vsub(m_nq, m_cellSol[m_gates[j]], 1, m_wsp[m_gates[j]], 1, m_cellSol[m_gates[j]], 1);
             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);
         }
 
         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

Implemented in Nektar::CourtemancheRamirezNattel98, Nektar::FentonKarma, Nektar::Winslow99, Nektar::LuoRudy91, Nektar::Fox02, Nektar::TenTusscher06, Nektar::CellModelFitzHughNagumo, Nektar::CellModelAlievPanfilov, and Nektar::PanditGilesDemir03.

Referenced by GenerateSummary().

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

Implemented in Nektar::CourtemancheRamirezNattel98, Nektar::Winslow99, Nektar::FentonKarma, Nektar::LuoRudy91, Nektar::Fox02, Nektar::TenTusscher06, Nektar::CellModelFitzHughNagumo, and Nektar::CellModelAlievPanfilov.

Referenced by Initialise().

virtual void Nektar::CellModel::v_Update	(	const Array< OneD, const Array< OneD, NekDouble > > &	inarray,
		Array< OneD, Array< OneD, NekDouble > > &	outarray,
		const NekDouble	time
	)