Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Monodomain.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File Monodomain.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Monodomain cardiac electrophysiology homogenised model.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <iostream>
37 
41 
42 using namespace std;
43 
44 namespace Nektar
45 {
46  /**
47  * @class Monodomain
48  *
49  * Base model of cardiac electrophysiology of the form
50  * \f{align*}{
51  * \frac{\partial u}{\partial t} = \nabla^2 u + J_{ion},
52  * \f}
53  * where the reaction term, \f$J_{ion}\f$ is defined by a specific cell
54  * model.
55  *
56  * This implementation, at present, treats the reaction terms explicitly
57  * and the diffusive element implicitly.
58  */
59 
60  /**
61  * Registers the class with the Factory.
62  */
63  string Monodomain::className
65  "Monodomain",
66  Monodomain::create,
67  "Monodomain model of cardiac electrophysiology.");
68 
69 
70  /**
71  *
72  */
73  Monodomain::Monodomain(
75  : UnsteadySystem(pSession)
76  {
77  }
78 
79 
80  /**
81  *
82  */
84  {
86 
87  m_session->LoadParameter("Chi", m_chi);
88  m_session->LoadParameter("Cm", m_capMembrane);
89 
90  std::string vCellModel;
91  m_session->LoadSolverInfo("CELLMODEL", vCellModel, "");
92 
93  ASSERTL0(vCellModel != "", "Cell Model not specified.");
94 
96  vCellModel, m_session, m_fields[0]);
97 
98  m_intVariables.push_back(0);
99 
100  // Load variable coefficients
101  StdRegions::VarCoeffType varCoeffEnum[6] = {
108  };
109  std::string varCoeffString[6] = {"xx","xy","yy","xz","yz","zz"};
110  std::string aniso_var[3] = {"fx", "fy", "fz"};
111 
112  const int nq = m_fields[0]->GetNpoints();
113  const int nVarDiffCmpts = m_spacedim * (m_spacedim + 1) / 2;
114 
115  // Allocate storage for variable coeffs and initialize to 1.
116  for (int i = 0, k = 0; i < m_spacedim; ++i)
117  {
118  for (int j = 0; j < i+1; ++j)
119  {
120  if (i == j)
121  {
122  m_vardiff[varCoeffEnum[k]] = Array<OneD, NekDouble>(nq, 1.0);
123  }
124  else
125  {
126  m_vardiff[varCoeffEnum[k]] = Array<OneD, NekDouble>(nq, 0.0);
127  }
128  ++k;
129  }
130  }
131 
132  // Apply fibre map f \in [0,1], scale to conductivity range
133  // [o_min,o_max], specified by the session parameters o_min and o_max
134  if (m_session->DefinesFunction("AnisotropicConductivity"))
135  {
136  if (m_session->DefinesCmdLineArgument("verbose"))
137  {
138  cout << "Loading Anisotropic Fibre map." << endl;
139  }
140 
141  NekDouble o_min = m_session->GetParameter("o_min");
142  NekDouble o_max = m_session->GetParameter("o_max");
143  int k = 0;
144 
145  Array<OneD, NekDouble> vTemp_i;
146  Array<OneD, NekDouble> vTemp_j;
147 
148  /*
149  * Diffusivity matrix D is upper triangular and defined as
150  * d_00 d_01 d_02
151  * d_11 d_12
152  * d_22
153  *
154  * Given a principle fibre direction _f_ the diffusivity is given
155  * by
156  * d_ij = { D_2 + (D_1 - D_2) f_i f_j if i==j
157  * { (D_1 - D_2) f_i f_j if i!=j
158  *
159  * The vector _f_ is given in terms of the variables fx,fy,fz in the
160  * function AnisotropicConductivity. The values of D_1 and D_2 are
161  * the parameters o_max and o_min, respectively.
162  */
163 
164  // Loop through columns of D
165  for (int j = 0; j < m_spacedim; ++j)
166  {
167  ASSERTL0(m_session->DefinesFunction("AnisotropicConductivity",
168  aniso_var[j]),
169  "Function 'AnisotropicConductivity' not correctly "
170  "defined.");
171  EvaluateFunction(aniso_var[j], vTemp_j,
172  "AnisotropicConductivity");
173 
174  // Loop through rows of D
175  for (int i = 0; i < j + 1; ++i)
176  {
177  ASSERTL0(m_session->DefinesFunction(
178  "AnisotropicConductivity",aniso_var[i]),
179  "Function 'AnisotropicConductivity' not correctly "
180  "defined.");
181  EvaluateFunction(aniso_var[i], vTemp_i,
182  "AnisotropicConductivity");
183 
184  Vmath::Vmul(nq, vTemp_i, 1, vTemp_j, 1,
185  m_vardiff[varCoeffEnum[k]], 1);
186 
187  Vmath::Smul(nq, o_max-o_min,
188  m_vardiff[varCoeffEnum[k]], 1,
189  m_vardiff[varCoeffEnum[k]], 1);
190 
191  if (i == j)
192  {
193  Vmath::Sadd(nq, o_min,
194  m_vardiff[varCoeffEnum[k]], 1,
195  m_vardiff[varCoeffEnum[k]], 1);
196  }
197 
198  ++k;
199  }
200  }
201  }
202  else
203  {
204  // Otherwise apply isotropic conductivity value (o_max) to
205  // diagonal components of tensor
206  NekDouble o_max = m_session->GetParameter("o_max");
207  for (int i = 0; i < nVarDiffCmpts; ++i)
208  {
209  Vmath::Smul(nq,o_max,
210  m_vardiff[varCoeffEnum[i]], 1,
211  m_vardiff[varCoeffEnum[i]], 1);
212  }
213  }
214 
215  // Scale by scar map (range 0->1) derived from intensity map
216  // (range d_min -> d_max)
217  if (m_session->DefinesFunction("IsotropicConductivity"))
218  {
219  if (m_session->DefinesCmdLineArgument("verbose"))
220  {
221  cout << "Loading Isotropic Conductivity map." << endl;
222  }
223 
224  const std::string varName = "intensity";
226  EvaluateFunction(varName, vTemp, "IsotropicConductivity");
227 
228  // If the d_min and d_max parameters are defined, then we need to
229  // rescale the isotropic conductivity to convert from the source
230  // domain (e.g. late-gad intensity) to conductivity
231  if ( m_session->DefinesParameter("d_min") ||
232  m_session->DefinesParameter("d_max") ) {
233  const NekDouble f_min = m_session->GetParameter("d_min");
234  const NekDouble f_max = m_session->GetParameter("d_max");
235  const NekDouble scar_min = 0.0;
236  const NekDouble scar_max = 1.0;
237 
238  // Threshold based on d_min, d_max
239  for (int j = 0; j < nq; ++j)
240  {
241  vTemp[j] = (vTemp[j] < f_min ? f_min : vTemp[j]);
242  vTemp[j] = (vTemp[j] > f_max ? f_max : vTemp[j]);
243  }
244 
245  // Rescale to s \in [0,1] (0 maps to d_max, 1 maps to d_min)
246  Vmath::Sadd(nq, -f_min, vTemp, 1, vTemp, 1);
247  Vmath::Smul(nq, -1.0/(f_max-f_min), vTemp, 1, vTemp, 1);
248  Vmath::Sadd(nq, 1.0, vTemp, 1, vTemp, 1);
249  Vmath::Smul(nq, scar_max - scar_min, vTemp, 1, vTemp, 1);
250  Vmath::Sadd(nq, scar_min, vTemp, 1, vTemp, 1);
251  }
252 
253  // Scale anisotropic conductivity values
254  for (int i = 0; i < nVarDiffCmpts; ++i)
255  {
256  Vmath::Vmul(nq, vTemp, 1,
257  m_vardiff[varCoeffEnum[i]], 1,
258  m_vardiff[varCoeffEnum[i]], 1);
259  }
260  }
261 
262 
263  // Write out conductivity values
264  for (int j = 0, k = 0; j < m_spacedim; ++j)
265  {
266  // Loop through rows of D
267  for (int i = 0; i < j + 1; ++i)
268  {
269  // Transform variable coefficient and write out to file.
270  m_fields[0]->FwdTrans_IterPerExp(m_vardiff[varCoeffEnum[k]],
271  m_fields[0]->UpdateCoeffs());
272  std::stringstream filename;
273  filename << "Conductivity_" << varCoeffString[k] << ".fld";
274  WriteFld(filename.str());
275 
276  ++k;
277  }
278  }
279 
280  // Search through the loaded filters and pass the cell model to any
281  // CheckpointCellModel filters loaded.
282  int k = 0;
283  const LibUtilities::FilterMap& f = m_session->GetFilters();
284  LibUtilities::FilterMap::const_iterator x;
285  for (x = f.begin(); x != f.end(); ++x, ++k)
286  {
287  if (x->first == "CheckpointCellModel")
288  {
289  boost::shared_ptr<FilterCheckpointCellModel> c
290  = boost::dynamic_pointer_cast<FilterCheckpointCellModel>(
291  m_filters[k]);
292  c->SetCellModel(m_cell);
293  }
294  if (x->first == "CellHistoryPoints")
295  {
296  boost::shared_ptr<FilterCellHistoryPoints> c
297  = boost::dynamic_pointer_cast<FilterCellHistoryPoints>(
298  m_filters[k]);
299  c->SetCellModel(m_cell);
300  }
301  }
302 
303  // Load stimuli
305 
306  if (!m_explicitDiffusion)
307  {
309  }
311  }
312 
313 
314  /**
315  *
316  */
318  {
319 
320  }
321 
322 
323  /**
324  * @param inarray Input array.
325  * @param outarray Output array.
326  * @param time Current simulation time.
327  * @param lambda Timestep.
328  */
330  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
331  Array<OneD, Array<OneD, NekDouble> >&outarray,
332  const NekDouble time,
333  const NekDouble lambda)
334  {
335  int nvariables = inarray.num_elements();
336  int nq = m_fields[0]->GetNpoints();
338  // lambda = \Delta t
339  factors[StdRegions::eFactorLambda] = 1.0/lambda*m_chi*m_capMembrane;
340 
341  // We solve ( \nabla^2 - HHlambda ) Y[i] = rhs [i]
342  // inarray = input: \hat{rhs} -> output: \hat{Y}
343  // outarray = output: nabla^2 \hat{Y}
344  // where \hat = modal coeffs
345  for (int i = 0; i < nvariables; ++i)
346  {
347  // Multiply 1.0/timestep
348  Vmath::Smul(nq, -factors[StdRegions::eFactorLambda], inarray[i], 1,
349  m_fields[i]->UpdatePhys(), 1);
350 
351  // Solve a system of equations with Helmholtz solver and transform
352  // back into physical space.
353  m_fields[i]->HelmSolve(m_fields[i]->GetPhys(),
354  m_fields[i]->UpdateCoeffs(), NullFlagList,
355  factors, m_vardiff);
356 
357  m_fields[i]->BwdTrans( m_fields[i]->GetCoeffs(),
358  m_fields[i]->UpdatePhys());
359  m_fields[i]->SetPhysState(true);
360 
361  // Copy the solution vector (required as m_fields must be set).
362  outarray[i] = m_fields[i]->GetPhys();
363  }
364  }
365 
366 
367  /**
368  *
369  */
371  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
372  Array<OneD, Array<OneD, NekDouble> >&outarray,
373  const NekDouble time)
374  {
375  // Compute I_ion
376  m_cell->TimeIntegrate(inarray, outarray, time);
377 
378  // Compute I_stim
379  for (unsigned int i = 0; i < m_stimulus.size(); ++i)
380  {
381  m_stimulus[i]->Update(outarray, time);
382  }
383  }
384 
385 
386  /**
387  *
388  */
390  bool dumpInitialConditions,
391  const int domain)
392  {
394  dumpInitialConditions,
395  domain);
396  m_cell->Initialise();
397  }
398 
399 
400  /**
401  *
402  */
404  {
406  if (m_session->DefinesFunction("d00") &&
407  m_session->GetFunctionType("d00", "intensity")
409  {
410  AddSummaryItem(s, "Diffusivity-x",
411  m_session->GetFunction("d00", "intensity")->GetExpression());
412  }
413  if (m_session->DefinesFunction("d11") &&
414  m_session->GetFunctionType("d11", "intensity")
416  {
417  AddSummaryItem(s, "Diffusivity-y",
418  m_session->GetFunction("d11", "intensity")->GetExpression());
419  }
420  if (m_session->DefinesFunction("d22") &&
421  m_session->GetFunctionType("d22", "intensity")
423  {
424  AddSummaryItem(s, "Diffusivity-z",
425  m_session->GetFunction("d22", "intensity")->GetExpression());
426  }
427  m_cell->GenerateSummary(s);
428  }
429 }
void DoImplicitSolve(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, NekDouble time, NekDouble lambda)
Solve for the diffusion term.
Definition: Monodomain.cpp:329
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:162
bool m_explicitDiffusion
Indicates if explicit or implicit treatment of diffusion is used.
void DefineImplicitSolve(FuncPointerT func, ObjectPointerT obj)
virtual void v_InitObject()
Init object for UnsteadySystem class.
Definition: Monodomain.cpp:83
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:47
STL namespace.
StdRegions::VarCoeffMap m_vardiff
Variable diffusivity.
Definition: Monodomain.h:105
void SetCellModel(CellModelSharedPtr &pCellModel)
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:251
NekDouble m_capMembrane
Definition: Monodomain.h:108
void SetCellModel(CellModelSharedPtr &pCellModel)
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
static std::vector< StimulusSharedPtr > LoadStimuli(const LibUtilities::SessionReaderSharedPtr &pSession, const MultiRegions::ExpListSharedPtr &pField)
Definition: Stimulus.cpp:95
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s)
Print a summary of time stepping parameters.
CellModelSharedPtr m_cell
Cell model.
Definition: Monodomain.h:100
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
Base class for unsteady solvers.
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
Definition: Misc.cpp:50
virtual void v_GenerateSummary(SummaryList &s)
Prints a summary of the model parameters.
Definition: Monodomain.cpp:403
std::vector< std::pair< std::string, FilterParams > > FilterMap
Definition: SessionReader.h:66
int m_spacedim
Spatial dimension (>= expansion dim).
virtual ~Monodomain()
Desctructor.
Definition: Monodomain.cpp:317
double NekDouble
virtual SOLVER_UTILS_EXPORT void v_SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
virtual SOLVER_UTILS_EXPORT void v_InitObject()
Init object for UnsteadySystem class.
SOLVER_UTILS_EXPORT void EvaluateFunction(Array< OneD, Array< OneD, NekDouble > > &pArray, std::string pFunctionName, const NekDouble pTime=0.0, const int domain=0)
Evaluates a function as specified in the session file.
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
Definition: Vmath.cpp:301
EquationSystemFactory & GetEquationSystemFactory()
CellModelFactory & GetCellModelFactory()
Definition: CellModel.cpp:47
std::vector< StimulusSharedPtr > m_stimulus
Definition: Monodomain.h:102
SOLVER_UTILS_EXPORT void WriteFld(const std::string &outname)
Write field data to the given filename.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
virtual void v_SetInitialConditions(NekDouble initialtime, bool dumpInitialConditions, const int domain)
Sets a custom initial condition.
Definition: Monodomain.cpp:389
std::vector< FilterSharedPtr > m_filters
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Computes the reaction terms and .
Definition: Monodomain.cpp:370
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169
static FlagList NullFlagList
An empty flag list.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215