Nektar++
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 namespace Nektar
43 {
44  /**
45  * @class Monodomain
46  *
47  * Base model of cardiac electrophysiology of the form
48  * \f{align*}{
49  * \frac{\partial u}{\partial t} = \nabla^2 u + J_{ion},
50  * \f}
51  * where the reaction term, \f$J_{ion}\f$ is defined by a specific cell
52  * model.
53  *
54  * This implementation, at present, treats the reaction terms explicitly
55  * and the diffusive element implicitly.
56  */
57 
58  /**
59  * Registers the class with the Factory.
60  */
63  "Monodomain",
65  "Monodomain model of cardiac electrophysiology.");
66 
67 
68  /**
69  *
70  */
73  : UnsteadySystem(pSession)
74  {
75  }
76 
77 
78  /**
79  *
80  */
82  {
83  UnsteadySystem::v_InitObject();
84 
85  m_session->LoadParameter("Chi", m_chi);
86  m_session->LoadParameter("Cm", m_capMembrane);
87 
88  std::string vCellModel;
89  m_session->LoadSolverInfo("CELLMODEL", vCellModel, "");
90 
91  ASSERTL0(vCellModel != "", "Cell Model not specified.");
92 
94  vCellModel, m_session, m_fields[0]);
95 
96  m_intVariables.push_back(0);
97 
98  // Load variable coefficients
99  StdRegions::VarCoeffType varCoeffEnum[6] = {
106  };
107  std::string varCoeffString[6] = {"xx","xy","yy","xz","yz","zz"};
108  std::string aniso_var[3] = {"fx", "fy", "fz"};
109 
110  const int nq = m_fields[0]->GetNpoints();
111  const int nVarDiffCmpts = m_spacedim * (m_spacedim + 1) / 2;
112 
113  // Allocate storage for variable coeffs and initialize to 1.
114  for (int i = 0, k = 0; i < m_spacedim; ++i)
115  {
116  for (int j = 0; j < i+1; ++j)
117  {
118  if (i == j)
119  {
120  m_vardiff[varCoeffEnum[k]] = Array<OneD, NekDouble>(nq, 1.0);
121  }
122  else
123  {
124  m_vardiff[varCoeffEnum[k]] = Array<OneD, NekDouble>(nq, 0.0);
125  }
126  ++k;
127  }
128  }
129 
130  // Apply fibre map f \in [0,1], scale to conductivity range
131  // [o_min,o_max], specified by the session parameters o_min and o_max
132  if (m_session->DefinesFunction("AnisotropicConductivity"))
133  {
134  if (m_session->DefinesCmdLineArgument("verbose"))
135  {
136  cout << "Loading Anisotropic Fibre map." << endl;
137  }
138 
139  NekDouble o_min = m_session->GetParameter("o_min");
140  NekDouble o_max = m_session->GetParameter("o_max");
141  int k = 0;
142 
143  Array<OneD, NekDouble> vTemp_i;
144  Array<OneD, NekDouble> vTemp_j;
145 
146  /*
147  * Diffusivity matrix D is upper triangular and defined as
148  * d_00 d_01 d_02
149  * d_11 d_12
150  * d_22
151  *
152  * Given a principle fibre direction _f_ the diffusivity is given
153  * by
154  * d_ij = { D_2 + (D_1 - D_2) f_i f_j if i==j
155  * { (D_1 - D_2) f_i f_j if i!=j
156  *
157  * The vector _f_ is given in terms of the variables fx,fy,fz in the
158  * function AnisotropicConductivity. The values of D_1 and D_2 are
159  * the parameters o_max and o_min, respectively.
160  */
161 
162  // Loop through columns of D
163  for (int j = 0; j < m_spacedim; ++j)
164  {
165  ASSERTL0(m_session->DefinesFunction("AnisotropicConductivity",
166  aniso_var[j]),
167  "Function 'AnisotropicConductivity' not correctly "
168  "defined.");
169  EvaluateFunction(aniso_var[j], vTemp_j,
170  "AnisotropicConductivity");
171 
172  // Loop through rows of D
173  for (int i = 0; i < j + 1; ++i)
174  {
175  ASSERTL0(m_session->DefinesFunction(
176  "AnisotropicConductivity",aniso_var[i]),
177  "Function 'AnisotropicConductivity' not correctly "
178  "defined.");
179  EvaluateFunction(aniso_var[i], vTemp_i,
180  "AnisotropicConductivity");
181 
182  Vmath::Vmul(nq, vTemp_i, 1, vTemp_j, 1,
183  m_vardiff[varCoeffEnum[k]], 1);
184 
185  Vmath::Smul(nq, o_max-o_min,
186  m_vardiff[varCoeffEnum[k]], 1,
187  m_vardiff[varCoeffEnum[k]], 1);
188 
189  if (i == j)
190  {
191  Vmath::Sadd(nq, o_min,
192  m_vardiff[varCoeffEnum[k]], 1,
193  m_vardiff[varCoeffEnum[k]], 1);
194  }
195 
196  ++k;
197  }
198  }
199  }
200  else
201  {
202  // Otherwise apply isotropic conductivity value (o_max) to
203  // diagonal components of tensor
204  NekDouble o_max = m_session->GetParameter("o_max");
205  for (int i = 0; i < nVarDiffCmpts; ++i)
206  {
207  Vmath::Smul(nq,o_max,
208  m_vardiff[varCoeffEnum[i]], 1,
209  m_vardiff[varCoeffEnum[i]], 1);
210  }
211  }
212 
213  // Scale by scar map (range 0->1) derived from intensity map
214  // (range d_min -> d_max)
215  if (m_session->DefinesFunction("IsotropicConductivity"))
216  {
217  if (m_session->DefinesCmdLineArgument("verbose"))
218  {
219  cout << "Loading Isotropic Conductivity map." << endl;
220  }
221 
222  const std::string varName = "intensity";
224  EvaluateFunction(varName, vTemp, "IsotropicConductivity");
225 
226  // If the d_min and d_max parameters are defined, then we need to
227  // rescale the isotropic conductivity to convert from the source
228  // domain (e.g. late-gad intensity) to conductivity
229  if ( m_session->DefinesParameter("d_min") ||
230  m_session->DefinesParameter("d_max") ) {
231  const NekDouble f_min = m_session->GetParameter("d_min");
232  const NekDouble f_max = m_session->GetParameter("d_max");
233  const NekDouble scar_min = 0.0;
234  const NekDouble scar_max = 1.0;
235 
236  // Threshold based on d_min, d_max
237  for (int j = 0; j < nq; ++j)
238  {
239  vTemp[j] = (vTemp[j] < f_min ? f_min : vTemp[j]);
240  vTemp[j] = (vTemp[j] > f_max ? f_max : vTemp[j]);
241  }
242 
243  // Rescale to s \in [0,1] (0 maps to d_max, 1 maps to d_min)
244  Vmath::Sadd(nq, -f_min, vTemp, 1, vTemp, 1);
245  Vmath::Smul(nq, -1.0/(f_max-f_min), vTemp, 1, vTemp, 1);
246  Vmath::Sadd(nq, 1.0, vTemp, 1, vTemp, 1);
247  Vmath::Smul(nq, scar_max - scar_min, vTemp, 1, vTemp, 1);
248  Vmath::Sadd(nq, scar_min, vTemp, 1, vTemp, 1);
249  }
250 
251  // Scale anisotropic conductivity values
252  for (int i = 0; i < nVarDiffCmpts; ++i)
253  {
254  Vmath::Vmul(nq, vTemp, 1,
255  m_vardiff[varCoeffEnum[i]], 1,
256  m_vardiff[varCoeffEnum[i]], 1);
257  }
258  }
259 
260 
261  // Write out conductivity values
262  for (int j = 0, k = 0; j < m_spacedim; ++j)
263  {
264  // Loop through rows of D
265  for (int i = 0; i < j + 1; ++i)
266  {
267  // Transform variable coefficient and write out to file.
268  m_fields[0]->FwdTrans_IterPerExp(m_vardiff[varCoeffEnum[k]],
269  m_fields[0]->UpdateCoeffs());
270  std::stringstream filename;
271  filename << "Conductivity_" << varCoeffString[k] << ".fld";
272  WriteFld(filename.str());
273 
274  ++k;
275  }
276  }
277 
278  // Search through the loaded filters and pass the cell model to any
279  // CheckpointCellModel filters loaded.
280  int k = 0;
281  const LibUtilities::FilterMap& f = m_session->GetFilters();
282  LibUtilities::FilterMap::const_iterator x;
283  for (x = f.begin(); x != f.end(); ++x, ++k)
284  {
285  if (x->first == "CheckpointCellModel")
286  {
287  boost::shared_ptr<FilterCheckpointCellModel> c
288  = boost::dynamic_pointer_cast<FilterCheckpointCellModel>(
289  m_filters[k]);
290  c->SetCellModel(m_cell);
291  }
292  if (x->first == "CellHistoryPoints")
293  {
294  boost::shared_ptr<FilterCellHistoryPoints> c
295  = boost::dynamic_pointer_cast<FilterCellHistoryPoints>(
296  m_filters[k]);
297  c->SetCellModel(m_cell);
298  }
299  }
300 
301  // Load stimuli
303 
304  if (!m_explicitDiffusion)
305  {
307  }
309  }
310 
311 
312  /**
313  *
314  */
316  {
317 
318  }
319 
320 
321  /**
322  * @param inarray Input array.
323  * @param outarray Output array.
324  * @param time Current simulation time.
325  * @param lambda Timestep.
326  */
328  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
329  Array<OneD, Array<OneD, NekDouble> >&outarray,
330  const NekDouble time,
331  const NekDouble lambda)
332  {
333  int nvariables = inarray.num_elements();
334  int nq = m_fields[0]->GetNpoints();
336  // lambda = \Delta t
337  factors[StdRegions::eFactorLambda] = 1.0/lambda*m_chi*m_capMembrane;
338 
339  // We solve ( \nabla^2 - HHlambda ) Y[i] = rhs [i]
340  // inarray = input: \hat{rhs} -> output: \hat{Y}
341  // outarray = output: nabla^2 \hat{Y}
342  // where \hat = modal coeffs
343  for (int i = 0; i < nvariables; ++i)
344  {
345  // Multiply 1.0/timestep
346  Vmath::Smul(nq, -factors[StdRegions::eFactorLambda], inarray[i], 1,
347  m_fields[i]->UpdatePhys(), 1);
348 
349  // Solve a system of equations with Helmholtz solver and transform
350  // back into physical space.
351  m_fields[i]->HelmSolve(m_fields[i]->GetPhys(),
352  m_fields[i]->UpdateCoeffs(), NullFlagList,
353  factors, m_vardiff);
354 
355  m_fields[i]->BwdTrans( m_fields[i]->GetCoeffs(),
356  m_fields[i]->UpdatePhys());
357  m_fields[i]->SetPhysState(true);
358 
359  // Copy the solution vector (required as m_fields must be set).
360  outarray[i] = m_fields[i]->GetPhys();
361  }
362  }
363 
364 
365  /**
366  *
367  */
369  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
370  Array<OneD, Array<OneD, NekDouble> >&outarray,
371  const NekDouble time)
372  {
373  // Compute I_ion
374  m_cell->TimeIntegrate(inarray, outarray, time);
375 
376  // Compute I_stim
377  for (unsigned int i = 0; i < m_stimulus.size(); ++i)
378  {
379  m_stimulus[i]->Update(outarray, time);
380  }
381  }
382 
383 
384  /**
385  *
386  */
388  bool dumpInitialConditions,
389  const int domain)
390  {
391  EquationSystem::v_SetInitialConditions(initialtime,
392  dumpInitialConditions,
393  domain);
394  m_cell->Initialise();
395  }
396 
397 
398  /**
399  *
400  */
402  {
403  UnsteadySystem::v_GenerateSummary(s);
404  if (m_session->DefinesFunction("d00") &&
405  m_session->GetFunctionType("d00", "intensity")
407  {
408  AddSummaryItem(s, "Diffusivity-x",
409  m_session->GetFunction("d00", "intensity")->GetExpression());
410  }
411  if (m_session->DefinesFunction("d11") &&
412  m_session->GetFunctionType("d11", "intensity")
414  {
415  AddSummaryItem(s, "Diffusivity-y",
416  m_session->GetFunction("d11", "intensity")->GetExpression());
417  }
418  if (m_session->DefinesFunction("d22") &&
419  m_session->GetFunctionType("d22", "intensity")
421  {
422  AddSummaryItem(s, "Diffusivity-z",
423  m_session->GetFunction("d22", "intensity")->GetExpression());
424  }
425  m_cell->GenerateSummary(s);
426  }
427 }
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:327
#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:81
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
Monodomain(const LibUtilities::SessionReaderSharedPtr &pSession)
Constructor.
Definition: Monodomain.cpp:71
StdRegions::VarCoeffMap m_vardiff
Variable diffusivity.
Definition: Monodomain.h:105
void SetCellModel(CellModelSharedPtr &pCellModel)
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:248
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: CellModel.h:51
NekDouble m_capMembrane
Definition: Monodomain.h:108
void SetCellModel(CellModelSharedPtr &pCellModel)
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
static std::vector< StimulusSharedPtr > LoadStimuli(const LibUtilities::SessionReaderSharedPtr &pSession, const MultiRegions::ExpListSharedPtr &pField)
Definition: Stimulus.cpp:93
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:401
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:315
double NekDouble
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:45
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.
static std::string className
Name of class.
Definition: Monodomain.h:65
virtual void v_SetInitialConditions(NekDouble initialtime, bool dumpInitialConditions, const int domain)
Sets a custom initial condition.
Definition: Monodomain.cpp:387
std::vector< FilterSharedPtr > m_filters
static EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession)
Creates an instance of this class.
Definition: Monodomain.h:56
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:368
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