#include <FentonKarma.h>

Inheritance diagram for Nektar::FentonKarma:

Collaboration diagram for Nektar::FentonKarma:

Public Member Functions
	FentonKarma (const LibUtilities::SessionReaderSharedPtr &pSession, const MultiRegions::ExpListSharedPtr &pField)
	Constructor. More...

virtual	~FentonKarma ()
	Destructor. More...

Public Member Functions inherited from Nektar::CellModel
	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)

Static Public Member Functions
static CellModelSharedPtr	create (const LibUtilities::SessionReaderSharedPtr &pSession, const MultiRegions::ExpListSharedPtr &pField)
	Creates an instance of this class. More...

Static Public Attributes
static std::string	className
	Name of class. More...

Protected Member Functions
virtual void	v_Update (const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
	Computes the reaction terms $f(u,v)$ and $g(u,v)$. More...

virtual void	v_GenerateSummary (SummaryList &s)
	Prints a summary of the model parameters. More...

virtual void	v_SetInitialConditions ()

virtual std::string	v_GetCellVarName (unsigned int idx)

Protected Member Functions inherited from Nektar::CellModel
void	LoadCellModel ()

Private Types
enum	Variants { eBR, eMBR, eMLR1, eGP, eCF1, eCF2a, eCF2b, eCF2c, eCF3a, eCF3b }

Private Attributes
NekDouble	C_m

NekDouble	V_0

NekDouble	u_fi

NekDouble	u_c

NekDouble	u_v

NekDouble	u_r

NekDouble	g_fi_max

NekDouble	tau_d

NekDouble	tau_v1_minus

NekDouble	tau_v2_minus

NekDouble	tau_v_plus

NekDouble	tau_0

NekDouble	tau_r

NekDouble	tau_si

NekDouble	tau_y_plus

NekDouble	tau_y_minus

NekDouble	u_csi

NekDouble	k1

NekDouble	k2

NekDouble	tau_w_minus

NekDouble	tau_w_plus

bool	isCF3

enum Variants	model_variant

Static Private Attributes
static std::string	lookupIds []

static std::string	def

Additional Inherited Members
Protected Attributes inherited from Nektar::CellModel
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

Definition at line 43 of file FentonKarma.h.

Member Enumeration Documentation

enum Nektar::FentonKarma::Variants

private

Enumerator
eBR
eMBR
eMLR1
eGP
eCF1
eCF2a
eCF2b
eCF2c
eCF3a
eCF3b

Definition at line 105 of file FentonKarma.h.

                       {
             eBR,
             eMBR,
             eMLR1,
             eGP,
             eCF1,
             eCF2a,
             eCF2b,
             eCF2c,
             eCF3a,
             eCF3b
         };

Constructor & Destructor Documentation

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

Constructor.

Initialise Fenton-Karma model parameters. k1 is k in the original model. k2 is an additional parameter introduced in Cherry&Fenton 2004. u_r and u_fi are introduced by Cherry&Fenton 2004 and is the same as u_c in the original model.

Definition at line 88 of file FentonKarma.cpp.

References ASSERTL0, C_m, eBR, eCF1, eCF2a, eCF2b, eCF2c, eCF3a, eCF3b, eGP, eMBR, eMLR1, g_fi_max, isCF3, k1, k2, Nektar::CellModel::m_gates, Nektar::CellModel::m_nvar, model_variant, tau_0, tau_d, tau_r, tau_si, tau_v1_minus, tau_v2_minus, tau_v_plus, tau_w_minus, tau_w_plus, tau_y_minus, tau_y_plus, u_c, u_csi, u_fi, u_r, u_v, and V_0.

         : CellModel(pSession, pField)
     {
         model_variant = pSession->GetSolverInfoAsEnum<FentonKarma::Variants>(
                                                             "CellModelVariant");
 
         C_m  =  1; // picoF
         V_0  = -85;
 
         switch (model_variant) {
             case eBR:
                 g_fi_max     = 4;
                 tau_r        = 33.33;
                 tau_si       = 29;
                 tau_0        = 12.5;
                 tau_v_plus   = 3.33;
                 tau_v1_minus = 1250;
                 tau_v2_minus = 19.6;
                 tau_w_plus   = 870;
                 tau_w_minus  = 41;
                 u_c          = 0.13;
                 u_v          = 0.04;
                 u_r          = 0.13;
                 u_fi         = 0.13;
                 u_csi        = 0.85;
                 k1           = 10;
                 k2           = 0.0;
                 break;
             case eMBR:
                 g_fi_max     = 4;
                 tau_r        = 50;
                 tau_si       = 44.84;
                 tau_0        = 8.3;
                 tau_v_plus   = 3.33;
                 tau_v1_minus = 1000;
                 tau_v2_minus = 19.2;
                 tau_w_plus   = 667;
                 tau_w_minus  = 11;
                 u_c          = 0.13;
                 u_v          = 0.04;
                 u_r          = 0.13;
                 u_fi         = 0.13;
                 u_csi        = 0.85;
                 k1           = 10;
                 k2           = 0.0;
                 break;
             case eMLR1:
                 g_fi_max     = 5.8;
                 tau_r        = 130;
                 tau_si       = 127;
                 tau_0        = 12.5;
                 tau_v_plus   = 10;
                 tau_v1_minus = 18.2;
                 tau_v2_minus = 18.2;
                 tau_w_plus   = 1020;
                 tau_w_minus  = 80;
                 u_c          = 0.13;
                 u_v          = 0.0;
                 u_r          = 0.13;
                 u_fi         = 0.13;
                 u_csi        = 0.85;
                 k1           = 10;
                 k2           = 0.0;
                 break;
             case eGP:
                 g_fi_max     = 8.7;
                 tau_r        = 25;
                 tau_si       = 22.22;
                 tau_0        = 12.5;
                 tau_v_plus   = 10;
                 tau_v1_minus = 333;
                 tau_v2_minus = 40;
                 tau_w_plus   = 1000;
                 tau_w_minus  = 65;
                 u_c          = 0.13;
                 u_v          = 0.025;
                 u_r          = 0.13;
                 u_fi         = 0.13;
                 u_csi        = 0.85;
                 k1           = 10;
                 k2           = 0.0;
                 break;
             case eCF1:
                 g_fi_max     = 6.6666;
                 tau_r        = 12.5;
                 tau_si       = 10;
                 tau_0        = 1.5;
                 tau_v_plus   = 10;
                 tau_v1_minus = 350;
                 tau_v2_minus = 80;
                 tau_w_plus   = 562;
                 tau_w_minus  = 48.5;
                 u_c          = 0.25;
                 u_v          = 0.001;
                 u_r          = 0.25;
                 u_fi         = 0.15;
                 u_csi        = 0.2;
                 k1           = 15;
                 k2           = 0;
                 break;
             case eCF2a:
                 g_fi_max     = 6.6666;
                 tau_r        = 31;
                 tau_si       = 26.5;
                 tau_0        = 1.5;
                 tau_v_plus   = 10;
                 tau_v1_minus = 20;
                 tau_v2_minus = 20;
                 tau_w_plus   = 800;
                 tau_w_minus  = 45;
                 u_c          = 0.25;
                 u_v          = 0.05;
                 u_r          = 0.6;
                 u_fi         = 0.11;
                 u_csi        = 0.7;
                 k1           = 10;
                 k2           = 1;
                 break;
             case eCF2b:
                 g_fi_max     = 6.6666;
                 tau_r        = 31;
                 tau_si       = 26.5;
                 tau_0        = 1.5;
                 tau_v_plus   = 10;
                 tau_v1_minus = 100;
                 tau_v2_minus = 20;
                 tau_w_plus   = 800;
                 tau_w_minus  = 45;
                 u_c          = 0.25;
                 u_v          = 0.05;
                 u_r          = 0.6;
                 u_fi         = 0.11;
                 u_csi        = 0.7;
                 k1           = 10;
                 k2           = 1;
                 break;
             case eCF2c:
                 g_fi_max     = 6.6666;
                 tau_r        = 31;
                 tau_si       = 26.5;
                 tau_0        = 1.5;
                 tau_v_plus   = 10;
                 tau_v1_minus = 150;
                 tau_v2_minus = 20;
                 tau_w_plus   = 800;
                 tau_w_minus  = 45;
                 u_c          = 0.25;
                 u_v          = 0.05;
                 u_r          = 0.6;
                 u_fi         = 0.11;
                 u_csi        = 0.7;
                 k1           = 10;
                 k2           = 1;
                 break;
             case eCF3a:
                 g_fi_max     = 13.3333;
                 tau_r        = 38;
                 tau_si       = 127;
                 tau_0        = 8.3;
                 tau_v_plus   = 3.33;
                 tau_v1_minus = 45;
                 tau_v2_minus = 300;
                 tau_w_plus   = 600;
                 tau_w_minus  = 40;
                 tau_y_plus   = 1000;
                 tau_y_minus  = 230;
                 u_c          = 0.25;
                 u_v          = 0.5;
                 u_r          = 0.25;
                 u_fi         = 0.25;
                 u_csi        = 0.7;
                 k1           = 60;
                 k2           = 0;
                 break;
             case eCF3b:
                 g_fi_max     = 13.3333;
                 tau_r        = 38;
                 tau_si       = 127;
                 tau_0        = 8.3;
                 tau_v_plus   = 3.33;
                 tau_v1_minus = 20;
                 tau_v2_minus = 300;
                 tau_w_plus   = 600;
                 tau_w_minus  = 40;
                 tau_y_plus   = 1000;
                 tau_y_minus  = 230;
                 u_c          = 0.25;
                 u_v          = 0.5;
                 u_r          = 0.25;
                 u_fi         = 0.25;
                 u_csi        = 0.7;
                 k1           = 60;
                 k2           = 0;
                 break;
         }
 
         tau_d  = C_m/g_fi_max;
         
         isCF3 = (model_variant == eCF3a || model_variant == eCF3b);
 
         // List gates and concentrations
         m_gates.push_back(1);
         m_gates.push_back(2);
         m_nvar = 3;
 
         // Cherry-Fenton 2004 Model 3 has extra gating variable
         if (isCF3)
         {
             m_gates.push_back(3);
             m_nvar = 4;
         }
 
         ASSERTL0(!isCF3, "Cherry-Fenton model 3 not implemented yet.");
     }

Nektar::FentonKarma::~FentonKarma ( )

virtual

Destructor.

Definition at line 310 of file FentonKarma.cpp.

     {
         
     }

Member Function Documentation

static CellModelSharedPtr Nektar::FentonKarma::create	(	const LibUtilities::SessionReaderSharedPtr &	pSession,
		const MultiRegions::ExpListSharedPtr &	pField
	)

inlinestatic

Creates an instance of this class.

Definition at line 47 of file FentonKarma.h.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr().

         {
             return MemoryManager<FentonKarma>
             ::AllocateSharedPtr(pSession, pField);
         }

void Nektar::FentonKarma::v_GenerateSummary ( SummaryList & s )

protectedvirtual

Prints a summary of the model parameters.

Implements Nektar::CellModel.

Definition at line 392 of file FentonKarma.cpp.

References Nektar::SolverUtils::AddSummaryItem(), lookupIds, and model_variant.

     {
         SolverUtils::AddSummaryItem(s, "Cell model", "Fenton Karma");
         SolverUtils::AddSummaryItem(s, "Cell model var.", lookupIds[model_variant]);
     }

std::string Nektar::FentonKarma::v_GetCellVarName ( unsigned int idx )

protectedvirtual

Reimplemented from Nektar::CellModel.

Definition at line 411 of file FentonKarma.cpp.

     {
         switch (idx)
         {
             case 0:  return "u";
             case 1:  return "v";
             case 2:  return "w";
             case 3:  return "y";
             default: return "unknown";
         }
     }

void Nektar::FentonKarma::v_SetInitialConditions ( )

protectedvirtual

Implements Nektar::CellModel.

Definition at line 399 of file FentonKarma.cpp.

References Vmath::Fill(), isCF3, Nektar::CellModel::m_cellSol, and Nektar::CellModel::m_nq.

     {
         Vmath::Fill(m_nq, 0.0,  m_cellSol[0],  1);
         Vmath::Fill(m_nq, 1.0,  m_cellSol[1],  1);
         Vmath::Fill(m_nq, 1.0,  m_cellSol[2],  1);
         if (isCF3)
         {
             Vmath::Fill(m_nq, 0.1,  m_cellSol[2],  1);
         }
         
     }

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

protectedvirtual

Computes the reaction terms $f(u,v)$ and $g(u,v)$.

Implements Nektar::CellModel.

Definition at line 317 of file FentonKarma.cpp.

References ASSERTL0, isCF3, k1, k2, Nektar::CellModel::m_gates_tau, Nektar::CellModel::m_nq, tau_0, tau_d, tau_r, tau_si, tau_v1_minus, tau_v2_minus, tau_v_plus, tau_w_minus, tau_w_plus, u_c, u_csi, u_fi, u_r, and u_v.

     {
         ASSERTL0(inarray.get() != outarray.get(),
                  "Must have different arrays for input and output.");
         
         // Variables
         //  0   u    membrane potential
         //  1   v    v gate
         //  2   w    w gate
         int n = m_nq;
         int i = 0;
 
         // Declare pointers
         const NekDouble *u = &inarray[0][0];
         const NekDouble *v = &inarray[1][0];
         const NekDouble *w = &inarray[2][0];
         const NekDouble *y = isCF3 ? &inarray[3][0] : 0;
         NekDouble *u_new   = &outarray[0][0];
         NekDouble *v_new   = &outarray[1][0];
         NekDouble *w_new   = &outarray[2][0];
         //NekDouble *y_new   = isCF3 ? &outarray[3][0] : 0;
         NekDouble *v_tau   = &m_gates_tau[0][0];
         NekDouble *w_tau   = &m_gates_tau[1][0];
         //NekDouble *y_tau   = isCF3 ? &m_gates_tau[2][0] : 0;
 
         // Temporary variables
         NekDouble J_fi, J_so, J_si, h1, h2, h3, alpha, beta;
 
         // Compute rates for each point in domain
         for (i = 0; i < n; ++i)
         {
             // Heavyside functions
             h1 = (*u < u_c) ? 0.0 : 1.0;
             h2 = (*u < u_v) ? 0.0 : 1.0;
             h3 = (*u < u_r) ? 0.0 : 1.0;
 
             // w-gate
             alpha = (1-h1)/tau_w_minus;
             beta = h1/tau_w_plus;
             *w_tau = 1.0 / (alpha + beta);
             *w_new = alpha * (*w_tau);
 
             // v-gate
             alpha = (1-h1)/(h2*tau_v1_minus + (1-h2)*tau_v2_minus);
             beta = h1/tau_v_plus;
             *v_tau = 1.0 / (alpha + beta);
             *v_new = alpha * (*v_tau);
 
             // y-gate
             if (isCF3)
             {
                 // TODO: implementation for y_tau and y_new
             }
 
             // J_fi
             J_fi = -(*v)*h1*(1 - *u)*(*u - u_fi)/tau_d;
 
             // J_so
             // added extra (1-k2*v) term from Cherry&Fenton 2004
             J_so = (*u)*(1-h3)*(1-k2*(*v))/tau_0 +
                     (isCF3 ? h3*(*u)*(*y)/tau_r : h3/tau_r);
 
             // J_si
             J_si = -(*w)*(1 + tanh(k1*(*u - u_csi)))/(2.0*tau_si);
 
             // u
             *u_new = -J_fi - J_so - J_si;
 
             ++u, ++v, ++w, ++u_new, ++v_new, ++w_new, ++v_tau, ++w_tau;
         }
     }

Member Data Documentation

NekDouble Nektar::FentonKarma::C_m

private

Definition at line 81 of file FentonKarma.h.

Referenced by FentonKarma().

std::string Nektar::FentonKarma::className

static

Initial value:

=
            GetCellModelFactory().RegisterCreatorFunction(
                                                    "FentonKarma",
                                                    FentonKarma::create,
                                                    "Phenomenological Model.")

Name of class.

Definition at line 56 of file FentonKarma.h.

std::string Nektar::FentonKarma::def

staticprivate

Initial value:

=
            LibUtilities::SessionReader::RegisterDefaultSolverInfo(
                    "CellModelVariant","BR")

Definition at line 120 of file FentonKarma.h.

NekDouble Nektar::FentonKarma::g_fi_max

private

Definition at line 87 of file FentonKarma.h.

Referenced by FentonKarma().

bool Nektar::FentonKarma::isCF3

private

Definition at line 103 of file FentonKarma.h.

Referenced by FentonKarma(), v_SetInitialConditions(), and v_Update().

NekDouble Nektar::FentonKarma::k1

private

Definition at line 98 of file FentonKarma.h.

Referenced by FentonKarma(), and v_Update().

NekDouble Nektar::FentonKarma::k2

private

Definition at line 99 of file FentonKarma.h.

Referenced by FentonKarma(), and v_Update().

std::string Nektar::FentonKarma::lookupIds

staticprivate

Initial value:

= {
            LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
                    "BR",   FentonKarma::eBR),
            LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
                    "MBR",  FentonKarma::eMBR),
            LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
                    "MLR1", FentonKarma::eMLR1),
            LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
                    "GP",   FentonKarma::eGP),
            LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
                    "CF1",  FentonKarma::eCF1),
            LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
                    "CF2a",  FentonKarma::eCF2a),
            LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
                    "CF2b",  FentonKarma::eCF2b),
            LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
                    "CF2c",  FentonKarma::eCF2c),
            LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
                    "CF3a",  FentonKarma::eCF3a),
            LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
                    "CF3b",  FentonKarma::eCF3b)
    }