Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
BidomainRoth.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File BidomainRoth.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: Bidomain cardiac electrophysiology model - Roth formulation.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <iostream>
37 
40 
41 using namespace std;
42 
43 namespace Nektar
44 {
45 
46 /**
47  * Registers the class with the Factory.
48  */
49 string BidomainRoth::className
51  "BidomainRoth",
52  BidomainRoth::create,
53  "Bidomain Roth model of cardiac electrophysiology.");
54 
55 
56 /**
57  *
58  */
59 BidomainRoth::BidomainRoth(
61  : UnsteadySystem(pSession)
62 {
63 }
64 
65 
66 /**
67  *
68  */
70 {
72 
73  m_session->LoadParameter("Chi", m_chi);
74  m_session->LoadParameter("Cm", m_capMembrane);
75 
76  std::string vCellModel;
77  m_session->LoadSolverInfo("CELLMODEL", vCellModel, "");
78 
79  ASSERTL0(vCellModel != "", "Cell Model not specified.");
80 
82  vCellModel, m_session, m_fields[0]);
83 
84  m_intVariables.push_back(0);
85 
86  // Load variable coefficients
87  StdRegions::VarCoeffType varCoeffEnum[6] = {
94  };
95  std::string varCoeffString[6] = {"xx","xy","yy","xz","yz","zz"};
96  std::string aniso_var[3] = {"fx", "fy", "fz"};
97 
98  const int nq = m_fields[0]->GetNpoints();
99 
100  // Allocate storage for variable coeffs and initialize to 1.
101  for (int i = 0, k = 0; i < m_spacedim; ++i)
102  {
103  for (int j = 0; j < i+1; ++j)
104  {
105  if (i == j)
106  {
107  m_vardiffi[varCoeffEnum[k]] = Array<OneD, NekDouble>(nq, 1.0);
108  m_vardiffe[varCoeffEnum[k]] = Array<OneD, NekDouble>(nq, 1.0);
109  m_vardiffie[varCoeffEnum[k]] = Array<OneD, NekDouble>(nq, 1.0);
110  }
111  else
112  {
113  m_vardiffi[varCoeffEnum[k]] = Array<OneD, NekDouble>(nq, 0.0);
114  m_vardiffe[varCoeffEnum[k]] = Array<OneD, NekDouble>(nq, 0.0);
115  m_vardiffie[varCoeffEnum[k]] = Array<OneD, NekDouble>(nq, 0.0);
116  }
117  ++k;
118  }
119  }
120 
121  // Apply fibre map f \in [0,1], scale to conductivity range
122  // [o_min,o_max], specified by the session parameters o_min and o_max
123  if (m_session->DefinesFunction("ExtracellularAnisotropicConductivity"))
124  {
125  if (m_session->DefinesCmdLineArgument("verbose"))
126  {
127  cout << "Loading Extracellular Anisotropic Fibre map." << endl;
128  }
129 
130  NekDouble o_min = m_session->GetParameter("o_min");
131  NekDouble o_max = m_session->GetParameter("o_max");
132  int k = 0;
133 
134  Array<OneD, NekDouble> vTemp_i;
135  Array<OneD, NekDouble> vTemp_j;
136 
137  /*
138  * Diffusivity matrix D is upper triangular and defined as
139  * d_00 d_01 d_02
140  * d_11 d_12
141  * d_22
142  *
143  * Given a principle fibre direction _f_ the diffusivity is given
144  * by
145  * d_ij = { D_2 + (D_1 - D_2) f_i f_j if i==j
146  * { (D_1 - D_2) f_i f_j if i!=j
147  *
148  * The vector _f_ is given in terms of the variables fx,fy,fz in the
149  * function AnisotropicConductivity. The values of D_1 and D_2 are
150  * the parameters o_max and o_min, respectively.
151  */
152 
153  // Loop through columns of D
154  for (int j = 0; j < m_spacedim; ++j)
155  {
156  ASSERTL0(m_session->DefinesFunction(
157  "ExtracellularAnisotropicConductivity",
158  aniso_var[j]),
159  "Function 'AnisotropicConductivity' not correctly "
160  "defined.");
161  EvaluateFunction(aniso_var[j], vTemp_j,
162  "ExtracellularAnisotropicConductivity");
163 
164  // Loop through rows of D
165  for (int i = 0; i < j + 1; ++i)
166  {
167  ASSERTL0(m_session->DefinesFunction(
168  "ExtracellularAnisotropicConductivity",
169  aniso_var[i]),
170  "Function 'ExtracellularAnisotropicConductivity' not "
171  "correctly defined.");
172 
173  EvaluateFunction(aniso_var[i], vTemp_i,
174  "ExtracellularAnisotropicConductivity");
175 
176  Vmath::Vmul(nq, vTemp_i, 1, vTemp_j, 1,
177  m_vardiffe[varCoeffEnum[k]], 1);
178 
179  Vmath::Smul(nq, o_max-o_min,
180  m_vardiffe[varCoeffEnum[k]], 1,
181  m_vardiffe[varCoeffEnum[k]], 1);
182 
183  if (i == j)
184  {
185  Vmath::Sadd(nq, o_min,
186  m_vardiffe[varCoeffEnum[k]], 1,
187  m_vardiffe[varCoeffEnum[k]], 1);
188  }
189 
190  }
191  }
192  }
193 
194  // Apply fibre map f \in [0,1], scale to conductivity range
195  // [o_min,o_max], specified by the session parameters o_min and o_max
196  if (m_session->DefinesFunction("IntracellularAnisotropicConductivity"))
197  {
198  if (m_session->DefinesCmdLineArgument("verbose"))
199  {
200  cout << "Loading Anisotropic Fibre map." << endl;
201  }
202 
203  NekDouble o_min = m_session->GetParameter("o_min");
204  NekDouble o_max = m_session->GetParameter("o_max");
205  int k = 0;
206 
207  Array<OneD, NekDouble> vTemp_i;
208  Array<OneD, NekDouble> vTemp_j;
209 
210  /*
211  * Diffusivity matrix D is upper triangular and defined as
212  * d_00 d_01 d_02
213  * d_11 d_12
214  * d_22
215  *
216  * Given a principle fibre direction _f_ the diffusivity is given
217  * by
218  * d_ij = { D_2 + (D_1 - D_2) f_i f_j if i==j
219  * { (D_1 - D_2) f_i f_j if i!=j
220  *
221  * The vector _f_ is given in terms of the variables fx,fy,fz in the
222  * function AnisotropicConductivity. The values of D_1 and D_2 are
223  * the parameters o_max and o_min, respectively.
224  */
225 
226  // Loop through columns of D
227  for (int j = 0; j < m_spacedim; ++j)
228  {
229  ASSERTL0(m_session->DefinesFunction(
230  "IntracellularAnisotropicConductivity",
231  aniso_var[j]),
232  "Function 'IntracellularAnisotropicConductivity' not "
233  "correctly defined.");
234 
235  EvaluateFunction(aniso_var[j], vTemp_j,
236  "IntracellularAnisotropicConductivity");
237 
238  // Loop through rows of D
239  for (int i = 0; i < j + 1; ++i)
240  {
241  ASSERTL0(m_session->DefinesFunction(
242  "IntracellularAnisotropicConductivity",
243  aniso_var[i]),
244  "Function 'IntracellularAnisotropicConductivity' not "
245  "correctly defined.");
246  EvaluateFunction(aniso_var[i], vTemp_i,
247  "IntracellularAnisotropicConductivity");
248 
249  Vmath::Vmul(nq, vTemp_i, 1, vTemp_j, 1,
250  m_vardiffi[varCoeffEnum[k]], 1);
251 
252  Vmath::Smul(nq, o_max-o_min,
253  m_vardiffi[varCoeffEnum[k]], 1,
254  m_vardiffi[varCoeffEnum[k]], 1);
255 
256  if (i == j)
257  {
258  Vmath::Sadd(nq, o_min,
259  m_vardiffi[varCoeffEnum[k]], 1,
260  m_vardiffi[varCoeffEnum[k]], 1);
261  }
262 
263  Vmath::Vadd(nq, m_vardiffe[varCoeffEnum[k]], 1,
264  m_vardiffi[varCoeffEnum[k]], 1,
265  m_vardiffie[varCoeffEnum[k]], 1);
266 
267  ++k;
268  }
269  }
270  }
271 
272 
273  // Write out conductivity values
274  for (int j = 0, k = 0; j < m_spacedim; ++j)
275  {
276  // Loop through rows of D
277  for (int i = 0; i < j + 1; ++i)
278  {
279  // Transform variable coefficient and write out to file.
280  m_fields[0]->FwdTrans_IterPerExp(m_vardiffi[varCoeffEnum[k]],
281  m_fields[0]->UpdateCoeffs());
282  std::stringstream filenamei;
283  filenamei << "IConductivity_" << varCoeffString[k] << ".fld";
284  WriteFld(filenamei.str());
285 
286  // Transform variable coefficient and write out to file.
287  m_fields[0]->FwdTrans_IterPerExp(m_vardiffe[varCoeffEnum[k]],
288  m_fields[0]->UpdateCoeffs());
289  std::stringstream filenamee;
290  filenamee << "EConductivity_" << varCoeffString[k] << ".fld";
291  WriteFld(filenamee.str());
292 
293  ++k;
294  }
295  }
296 
297  // Search through the loaded filters and pass the cell model to any
298  // CheckpointCellModel filters loaded.
299  int k = 0;
300  const LibUtilities::FilterMap& f = m_session->GetFilters();
301  LibUtilities::FilterMap::const_iterator x;
302  for (x = f.begin(); x != f.end(); ++x, ++k)
303  {
304  if (x->first == "CheckpointCellModel")
305  {
306  boost::shared_ptr<FilterCheckpointCellModel> c
307  = boost::dynamic_pointer_cast<FilterCheckpointCellModel>(
308  m_filters[k]);
309  c->SetCellModel(m_cell);
310  }
311  }
312 
313  // Load stimuli
315 
316  if (!m_explicitDiffusion)
317  {
319  }
321 }
322 
323 
324 /**
325  *
326  */
328 {
329 
330 }
331 
332 
333 /**
334  * @param inarray Input array.
335  * @param outarray Output array.
336  * @param time Current simulation time.
337  * @param lambda Timestep.
338  */
340  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
341  Array<OneD, Array<OneD, NekDouble> >&outarray,
342  const NekDouble time,
343  const NekDouble lambda)
344 {
345  int nq = m_fields[0]->GetNpoints();
346 
347  StdRegions::ConstFactorMap factorsHelmholtz;
348  // lambda = \Delta t
349  factorsHelmholtz[StdRegions::eFactorLambda]
350  = 1.0/lambda*m_chi*m_capMembrane;
351 
352  // ------------------------------
353  // Solve Helmholtz problem for Vm
354  // ------------------------------
355  // Multiply 1.0/timestep
356  //Vmath::Vadd(nq, inarray[0], 1, ggrad, 1, m_fields[0]->UpdatePhys(), 1);
357  Vmath::Smul(nq, -factorsHelmholtz[StdRegions::eFactorLambda], inarray[0], 1,
358  m_fields[0]->UpdatePhys(), 1);
359 
360  // Solve a system of equations with Helmholtz solver and transform
361  // back into physical space.
362  m_fields[0]->HelmSolve(m_fields[0]->GetPhys(),
363  m_fields[0]->UpdateCoeffs(), NullFlagList,
364  factorsHelmholtz, m_vardiffe);
365 
366  m_fields[0]->BwdTrans( m_fields[0]->GetCoeffs(),
367  m_fields[0]->UpdatePhys());
368  m_fields[0]->SetPhysState(true);
369 
370  // Copy the solution vector (required as m_fields must be set).
371  outarray[0] = m_fields[0]->GetPhys();
372 }
373 
374 
375 /**
376  *
377  */
379  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
380  Array<OneD, Array<OneD, NekDouble> >&outarray,
381  const NekDouble time)
382 {
383  int nq = m_fields[0]->GetNpoints();
384 
385  // Compute I_ion
386  m_cell->TimeIntegrate(inarray, outarray, time);
387 
388  // Compute I_stim
389  for (unsigned int i = 0; i < m_stimulus.size(); ++i)
390  {
391  m_stimulus[i]->Update(outarray, time);
392  }
393 
394  Array<OneD, NekDouble> ggrad0(nq), ggrad1(nq), ggrad2(nq), ggrad(nq);
395  StdRegions::ConstFactorMap factorsPoisson;
396  factorsPoisson[StdRegions::eFactorLambda] = 0.0;
397 
398  // ----------------------------
399  // Compute \nabla g_i \nabla Vm
400  // ----------------------------
401  m_fields[0]->PhysDeriv(inarray[0], ggrad0, ggrad1, ggrad2);
402  m_fields[0]->PhysDeriv(0, ggrad0, ggrad0);
403  m_fields[0]->PhysDeriv(1, ggrad1, ggrad1);
404  m_fields[0]->PhysDeriv(2, ggrad2, ggrad2);
405  if (m_session->DefinesFunction("IntracellularAnisotropicConductivity") &&
406  m_session->DefinesFunction("ExtracellularAnisotropicConductivity"))
407  {
409  1, ggrad0, 1);
411  1, ggrad1, 1);
413  1, ggrad2, 1);
414  }
415  // Add partial derivatives together
416  Vmath::Vadd(nq, ggrad0, 1, ggrad1, 1, ggrad, 1);
417  Vmath::Vadd(nq, ggrad2, 1, ggrad, 1, ggrad, 1);
418 
419  Vmath::Smul(nq, -1.0, ggrad, 1, m_fields[1]->UpdatePhys(), 1);
420 
421  // ----------------------------
422  // Solve Poisson problem for Ve
423  // ----------------------------
424  m_fields[1]->HelmSolve(m_fields[1]->GetPhys(),
425  m_fields[1]->UpdateCoeffs(), NullFlagList, factorsPoisson,
426  m_vardiffie);
427  m_fields[1]->BwdTrans(m_fields[1]->GetCoeffs(),
428  m_fields[1]->UpdatePhys());
429  m_fields[1]->SetPhysState(true);
430 
431  // ------------------------------
432  // Compute Laplacian of Ve (forcing term)
433  // ------------------------------
434  m_fields[1]->PhysDeriv(m_fields[1]->GetPhys(), ggrad0, ggrad1, ggrad2);
435  m_fields[1]->PhysDeriv(0, ggrad0, ggrad0);
436  m_fields[1]->PhysDeriv(1, ggrad1, ggrad1);
437  m_fields[1]->PhysDeriv(2, ggrad2, ggrad2);
438  if (m_session->DefinesFunction("IntracellularAnisotropicConductivity") &&
439  m_session->DefinesFunction("ExtracellularAnisotropicConductivity"))
440  {
442  1, ggrad0, 1);
444  1, ggrad1, 1);
446  1, ggrad2, 1);
447  }
448  // Add partial derivatives together
449  Vmath::Vadd(nq, ggrad0, 1, ggrad1, 1, ggrad, 1);
450  Vmath::Vadd(nq, ggrad2, 1, ggrad, 1, ggrad, 1);
451 
452  Vmath::Vadd(nq, ggrad, 1, outarray[0], 1, outarray[0], 1);
453 }
454 
455 
456 /**
457  *
458  */
460  bool dumpInitialConditions,
461  const int domain)
462 {
464  dumpInitialConditions,
465  domain);
466  m_cell->Initialise();
467 }
468 
469 
470 /**
471  *
472  */
474 {
476  m_cell->GenerateSummary(s);
477 }
478 
479 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Computes the reaction terms and .
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)
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:47
virtual void v_InitObject()
Init object for UnsteadySystem class.
virtual void v_GenerateSummary(SummaryList &s)
Prints a summary of the model parameters.
STL namespace.
void SetCellModel(CellModelSharedPtr &pCellModel)
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:252
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
StdRegions::VarCoeffMap m_vardiffe
Definition: BidomainRoth.h:104
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.
std::vector< StimulusSharedPtr > m_stimulus
Definition: BidomainRoth.h:101
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.
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:213
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
Base class for unsteady solvers.
std::vector< std::pair< std::string, FilterParams > > FilterMap
Definition: SessionReader.h:66
int m_spacedim
Spatial dimension (>= expansion dim).
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:315
StdRegions::VarCoeffMap m_vardiffie
Definition: BidomainRoth.h:105
EquationSystemFactory & GetEquationSystemFactory()
CellModelFactory & GetCellModelFactory()
Definition: CellModel.cpp:47
SOLVER_UTILS_EXPORT void WriteFld(const std::string &outname)
Write field data to the given filename.
CellModelSharedPtr m_cell
Cell model.
Definition: BidomainRoth.h:99
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
StdRegions::VarCoeffMap m_vardiffi
Definition: BidomainRoth.h:103
std::vector< FilterSharedPtr > m_filters
virtual void v_SetInitialConditions(NekDouble initialtime, bool dumpInitialConditions, const int domain)
Sets a custom initial condition.
virtual ~BidomainRoth()
Desctructor.
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:299
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:183
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