#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.
virtual	~FentonKarma ()
	Destructor.
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.
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, 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.

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

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)$.
virtual void	v_GenerateSummary (SummaryList &s)
	Prints a summary of the model parameters.
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.
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

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.

        {
            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)
    }