Nektar++
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 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: Bidomain cardiac electrophysiology model - Roth formulation.
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <iostream>
36 
39 
40 using namespace std;
41 
42 namespace Nektar
43 {
44 
45 /**
46  * Registers the class with the Factory.
47  */
48 string BidomainRoth::className =
50  "BidomainRoth", BidomainRoth::create,
51  "Bidomain Roth model of cardiac electrophysiology.");
52 
53 /**
54  *
55  */
56 BidomainRoth::BidomainRoth(const LibUtilities::SessionReaderSharedPtr &pSession,
58  : UnsteadySystem(pSession, pGraph)
59 {
60 }
61 
62 /**
63  *
64  */
66 {
68 
69  m_session->LoadParameter("Chi", m_chi);
70  m_session->LoadParameter("Cm", m_capMembrane);
71 
72  std::string vCellModel;
73  m_session->LoadSolverInfo("CELLMODEL", vCellModel, "");
74 
75  ASSERTL0(vCellModel != "", "Cell Model not specified.");
76 
78  m_fields[0]);
79 
80  m_intVariables.push_back(0);
81 
82  // Load variable coefficients
83  StdRegions::VarCoeffType varCoeffEnum[6] = {
87  std::string varCoeffString[6] = {"xx", "xy", "yy", "xz", "yz", "zz"};
88  std::string aniso_var[3] = {"fx", "fy", "fz"};
89 
90  const int nq = m_fields[0]->GetNpoints();
91 
92  // Allocate storage for variable coeffs and initialize to 1.
93  for (int i = 0, k = 0; i < m_spacedim; ++i)
94  {
95  for (int j = 0; j < i + 1; ++j)
96  {
97  if (i == j)
98  {
99  m_vardiffi[varCoeffEnum[k]] = Array<OneD, NekDouble>(nq, 1.0);
100  m_vardiffe[varCoeffEnum[k]] = Array<OneD, NekDouble>(nq, 1.0);
101  m_vardiffie[varCoeffEnum[k]] = Array<OneD, NekDouble>(nq, 1.0);
102  }
103  else
104  {
105  m_vardiffi[varCoeffEnum[k]] = Array<OneD, NekDouble>(nq, 0.0);
106  m_vardiffe[varCoeffEnum[k]] = Array<OneD, NekDouble>(nq, 0.0);
107  m_vardiffie[varCoeffEnum[k]] = Array<OneD, NekDouble>(nq, 0.0);
108  }
109  ++k;
110  }
111  }
112 
113  // Apply fibre map f \in [0,1], scale to conductivity range
114  // [o_min,o_max], specified by the session parameters o_min and o_max
115  if (m_session->DefinesFunction("ExtracellularAnisotropicConductivity"))
116  {
117  if (m_session->DefinesCmdLineArgument("verbose"))
118  {
119  cout << "Loading Extracellular Anisotropic Fibre map." << endl;
120  }
121 
122  NekDouble o_min = m_session->GetParameter("o_min");
123  NekDouble o_max = m_session->GetParameter("o_max");
124  int k = 0;
125 
126  Array<OneD, NekDouble> vTemp_i;
127  Array<OneD, NekDouble> vTemp_j;
128 
129  /*
130  * Diffusivity matrix D is upper triangular and defined as
131  * d_00 d_01 d_02
132  * d_11 d_12
133  * d_22
134  *
135  * Given a principle fibre direction _f_ the diffusivity is given
136  * by
137  * d_ij = { D_2 + (D_1 - D_2) f_i f_j if i==j
138  * { (D_1 - D_2) f_i f_j if i!=j
139  *
140  * The vector _f_ is given in terms of the variables fx,fy,fz in the
141  * function AnisotropicConductivity. The values of D_1 and D_2 are
142  * the parameters o_max and o_min, respectively.
143  */
144 
145  // Loop through columns of D
146  for (int j = 0; j < m_spacedim; ++j)
147  {
148  ASSERTL0(m_session->DefinesFunction(
149  "ExtracellularAnisotropicConductivity", aniso_var[j]),
150  "Function 'AnisotropicConductivity' not correctly "
151  "defined.");
152 
153  GetFunction("ExtracellularAnisotropicConductivity")
154  ->Evaluate(aniso_var[j], vTemp_j);
155 
156  // Loop through rows of D
157  for (int i = 0; i < j + 1; ++i)
158  {
159  ASSERTL0(
160  m_session->DefinesFunction(
161  "ExtracellularAnisotropicConductivity", aniso_var[i]),
162  "Function 'ExtracellularAnisotropicConductivity' not "
163  "correctly defined.");
164 
165  GetFunction("ExtracellularAnisotropicConductivity")
166  ->Evaluate(aniso_var[i], vTemp_i);
167 
168  Vmath::Vmul(nq, vTemp_i, 1, vTemp_j, 1,
169  m_vardiffe[varCoeffEnum[k]], 1);
170 
171  Vmath::Smul(nq, o_max - o_min, m_vardiffe[varCoeffEnum[k]], 1,
172  m_vardiffe[varCoeffEnum[k]], 1);
173 
174  if (i == j)
175  {
176  Vmath::Sadd(nq, o_min, m_vardiffe[varCoeffEnum[k]], 1,
177  m_vardiffe[varCoeffEnum[k]], 1);
178  }
179  }
180  }
181  }
182 
183  // Apply fibre map f \in [0,1], scale to conductivity range
184  // [o_min,o_max], specified by the session parameters o_min and o_max
185  if (m_session->DefinesFunction("IntracellularAnisotropicConductivity"))
186  {
187  if (m_session->DefinesCmdLineArgument("verbose"))
188  {
189  cout << "Loading Anisotropic Fibre map." << endl;
190  }
191 
192  NekDouble o_min = m_session->GetParameter("o_min");
193  NekDouble o_max = m_session->GetParameter("o_max");
194  int k = 0;
195 
196  Array<OneD, NekDouble> vTemp_i;
197  Array<OneD, NekDouble> vTemp_j;
198 
199  /*
200  * Diffusivity matrix D is upper triangular and defined as
201  * d_00 d_01 d_02
202  * d_11 d_12
203  * d_22
204  *
205  * Given a principle fibre direction _f_ the diffusivity is given
206  * by
207  * d_ij = { D_2 + (D_1 - D_2) f_i f_j if i==j
208  * { (D_1 - D_2) f_i f_j if i!=j
209  *
210  * The vector _f_ is given in terms of the variables fx,fy,fz in the
211  * function AnisotropicConductivity. The values of D_1 and D_2 are
212  * the parameters o_max and o_min, respectively.
213  */
214 
215  // Loop through columns of D
216  for (int j = 0; j < m_spacedim; ++j)
217  {
218  ASSERTL0(m_session->DefinesFunction(
219  "IntracellularAnisotropicConductivity", aniso_var[j]),
220  "Function 'IntracellularAnisotropicConductivity' not "
221  "correctly defined.");
222 
223  GetFunction("IntracellularAnisotropicConductivity")
224  ->Evaluate(aniso_var[j], vTemp_j);
225 
226  // Loop through rows of D
227  for (int i = 0; i < j + 1; ++i)
228  {
229  ASSERTL0(
230  m_session->DefinesFunction(
231  "IntracellularAnisotropicConductivity", aniso_var[i]),
232  "Function 'IntracellularAnisotropicConductivity' not "
233  "correctly defined.");
234  GetFunction("IntracellularAnisotropicConductivity")
235  ->Evaluate(aniso_var[i], vTemp_i);
236 
237  Vmath::Vmul(nq, vTemp_i, 1, vTemp_j, 1,
238  m_vardiffi[varCoeffEnum[k]], 1);
239 
240  Vmath::Smul(nq, o_max - o_min, m_vardiffi[varCoeffEnum[k]], 1,
241  m_vardiffi[varCoeffEnum[k]], 1);
242 
243  if (i == j)
244  {
245  Vmath::Sadd(nq, o_min, m_vardiffi[varCoeffEnum[k]], 1,
246  m_vardiffi[varCoeffEnum[k]], 1);
247  }
248 
249  Vmath::Vadd(nq, m_vardiffe[varCoeffEnum[k]], 1,
250  m_vardiffi[varCoeffEnum[k]], 1,
251  m_vardiffie[varCoeffEnum[k]], 1);
252 
253  ++k;
254  }
255  }
256  }
257 
258  // Write out conductivity values
259  for (int j = 0, k = 0; j < m_spacedim; ++j)
260  {
261  // Loop through rows of D
262  for (int i = 0; i < j + 1; ++i)
263  {
264  // Transform variable coefficient and write out to file.
265  m_fields[0]->FwdTransLocalElmt(m_vardiffi[varCoeffEnum[k]],
266  m_fields[0]->UpdateCoeffs());
267  std::stringstream filenamei;
268  filenamei << "IConductivity_" << varCoeffString[k] << ".fld";
269  WriteFld(filenamei.str());
270 
271  // Transform variable coefficient and write out to file.
272  m_fields[0]->FwdTransLocalElmt(m_vardiffe[varCoeffEnum[k]],
273  m_fields[0]->UpdateCoeffs());
274  std::stringstream filenamee;
275  filenamee << "EConductivity_" << varCoeffString[k] << ".fld";
276  WriteFld(filenamee.str());
277 
278  ++k;
279  }
280  }
281 
282  // Search through the loaded filters and pass the cell model to any
283  // CheckpointCellModel filters loaded.
284  for (auto &x : m_filters)
285  {
286  if (x.first == "CheckpointCellModel")
287  {
288  std::shared_ptr<FilterCheckpointCellModel> c =
289  std::dynamic_pointer_cast<FilterCheckpointCellModel>(x.second);
290  c->SetCellModel(m_cell);
291  }
292  }
293  // Load stimuli
295 
296  if (!m_explicitDiffusion)
297  {
299  }
301 }
302 
303 /**
304  *
305  */
307 {
308 }
309 
310 /**
311  * @param inarray Input array.
312  * @param outarray Output array.
313  * @param time Current simulation time.
314  * @param lambda Timestep.
315  */
317  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
318  Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time,
319  const NekDouble lambda)
320 {
321  int nq = m_fields[0]->GetNpoints();
322 
323  StdRegions::ConstFactorMap factorsHelmholtz;
324  // lambda = \Delta t
325  factorsHelmholtz[StdRegions::eFactorLambda] =
326  1.0 / lambda * m_chi * m_capMembrane;
327 
328  // ------------------------------
329  // Solve Helmholtz problem for Vm
330  // ------------------------------
331  // Multiply 1.0/timestep
332  // Vmath::Vadd(nq, inarray[0], 1, ggrad, 1, m_fields[0]->UpdatePhys(), 1);
333  Vmath::Smul(nq, -factorsHelmholtz[StdRegions::eFactorLambda], inarray[0], 1,
334  m_fields[0]->UpdatePhys(), 1);
335 
336  // Solve a system of equations with Helmholtz solver and transform
337  // back into physical space.
338  m_fields[0]->HelmSolve(m_fields[0]->GetPhys(), m_fields[0]->UpdateCoeffs(),
339  factorsHelmholtz, m_vardiffe);
340 
341  m_fields[0]->BwdTrans(m_fields[0]->GetCoeffs(), m_fields[0]->UpdatePhys());
342  m_fields[0]->SetPhysState(true);
343 
344  // Copy the solution vector (required as m_fields must be set).
345  outarray[0] = m_fields[0]->GetPhys();
346 }
347 
348 /**
349  *
350  */
352  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
353  Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time)
354 {
355  int nq = m_fields[0]->GetNpoints();
356 
357  // Compute I_ion
358  m_cell->TimeIntegrate(inarray, outarray, time);
359 
360  // Compute I_stim
361  for (unsigned int i = 0; i < m_stimulus.size(); ++i)
362  {
363  m_stimulus[i]->Update(outarray, time);
364  }
365 
366  Array<OneD, NekDouble> ggrad0(nq), ggrad1(nq), ggrad2(nq), ggrad(nq);
367  StdRegions::ConstFactorMap factorsPoisson;
368  factorsPoisson[StdRegions::eFactorLambda] = 0.0;
369 
370  // ----------------------------
371  // Compute \nabla g_i \nabla Vm
372  // ----------------------------
373  m_fields[0]->PhysDeriv(inarray[0], ggrad0, ggrad1, ggrad2);
374  m_fields[0]->PhysDeriv(0, ggrad0, ggrad0);
375  m_fields[0]->PhysDeriv(1, ggrad1, ggrad1);
376  m_fields[0]->PhysDeriv(2, ggrad2, ggrad2);
377  if (m_session->DefinesFunction("IntracellularAnisotropicConductivity") &&
378  m_session->DefinesFunction("ExtracellularAnisotropicConductivity"))
379  {
381  ggrad0, 1);
383  ggrad1, 1);
385  ggrad2, 1);
386  }
387  // Add partial derivatives together
388  Vmath::Vadd(nq, ggrad0, 1, ggrad1, 1, ggrad, 1);
389  Vmath::Vadd(nq, ggrad2, 1, ggrad, 1, ggrad, 1);
390 
391  Vmath::Smul(nq, -1.0, ggrad, 1, m_fields[1]->UpdatePhys(), 1);
392 
393  // ----------------------------
394  // Solve Poisson problem for Ve
395  // ----------------------------
396  m_fields[1]->HelmSolve(m_fields[1]->GetPhys(), m_fields[1]->UpdateCoeffs(),
397  factorsPoisson, m_vardiffie);
398  m_fields[1]->BwdTrans(m_fields[1]->GetCoeffs(), m_fields[1]->UpdatePhys());
399  m_fields[1]->SetPhysState(true);
400 
401  // ------------------------------
402  // Compute Laplacian of Ve (forcing term)
403  // ------------------------------
404  m_fields[1]->PhysDeriv(m_fields[1]->GetPhys(), ggrad0, ggrad1, ggrad2);
405  m_fields[1]->PhysDeriv(0, ggrad0, ggrad0);
406  m_fields[1]->PhysDeriv(1, ggrad1, ggrad1);
407  m_fields[1]->PhysDeriv(2, ggrad2, ggrad2);
408  if (m_session->DefinesFunction("IntracellularAnisotropicConductivity") &&
409  m_session->DefinesFunction("ExtracellularAnisotropicConductivity"))
410  {
412  ggrad0, 1);
414  ggrad1, 1);
416  ggrad2, 1);
417  }
418  // Add partial derivatives together
419  Vmath::Vadd(nq, ggrad0, 1, ggrad1, 1, ggrad, 1);
420  Vmath::Vadd(nq, ggrad2, 1, ggrad, 1, ggrad, 1);
421 
422  Vmath::Vadd(nq, ggrad, 1, outarray[0], 1, outarray[0], 1);
423 }
424 
425 /**
426  *
427  */
429  bool dumpInitialConditions,
430  const int domain)
431 {
432  EquationSystem::v_SetInitialConditions(initialtime, dumpInitialConditions,
433  domain);
434  m_cell->Initialise();
435 }
436 
437 /**
438  *
439  */
441 {
443  m_cell->GenerateSummary(s);
444 }
445 
446 } // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
virtual ~BidomainRoth()
Desctructor.
virtual void v_GenerateSummary(SummaryList &s)
Prints a summary of the model parameters.
StdRegions::VarCoeffMap m_vardiffi
Definition: BidomainRoth.h:100
std::vector< StimulusSharedPtr > m_stimulus
Definition: BidomainRoth.h:98
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.
CellModelSharedPtr m_cell
Cell model.
Definition: BidomainRoth.h:96
virtual void v_SetInitialConditions(NekDouble initialtime, bool dumpInitialConditions, const int domain)
Sets a custom initial condition.
StdRegions::VarCoeffMap m_vardiffie
Definition: BidomainRoth.h:102
virtual void v_InitObject()
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
Computes the reaction terms and .
StdRegions::VarCoeffMap m_vardiffe
Definition: BidomainRoth.h:101
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:198
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:144
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
void DefineImplicitSolve(FuncPointerT func, ObjectPointerT obj)
int m_spacedim
Spatial dimension (>= expansion dim).
virtual SOLVER_UTILS_EXPORT void v_SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
SOLVER_UTILS_EXPORT void WriteFld(const std::string &outname)
Write field data to the given filename.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction(std::string name, const MultiRegions::ExpListSharedPtr &field=MultiRegions::NullExpListSharedPtr, bool cache=false)
Get a SessionFunction by name.
Base class for unsteady solvers.
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
std::vector< std::pair< std::string, FilterSharedPtr > > m_filters
bool m_explicitDiffusion
Indicates if explicit or implicit treatment of diffusion is used.
virtual SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareField=true)
Init object for UnsteadySystem class.
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s)
Print a summary of time stepping parameters.
static std::vector< StimulusSharedPtr > LoadStimuli(const LibUtilities::SessionReaderSharedPtr &pSession, const MultiRegions::ExpListSharedPtr &pField)
Definition: Stimulus.cpp:89
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:48
EquationSystemFactory & GetEquationSystemFactory()
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:172
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:282
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
CellModelFactory & GetCellModelFactory()
Definition: CellModel.cpp:46
double NekDouble
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:209
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:359
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:248
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:384