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