Memosa-FVM  0.2
BatteryModel_impl.h
Go to the documentation of this file.
1 // This file os part of FVM
2 // Copyright (c) 2012 FVM Authors
3 // See LICENSE file for terms.
4 
5 #include "Mesh.h"
6 #include <sstream>
7 
8 #include "NumType.h"
9 #include "Array.h"
10 #include "Field.h"
11 #include "CRConnectivity.h"
12 #include "LinearSystem.h"
13 //#include "FieldSet.h"
14 #include "StorageSite.h"
15 #include "MultiFieldMatrix.h"
16 #include "CRMatrix.h"
17 #include "FluxJacobianMatrix.h"
18 #include "DiagonalMatrix.h"
19 #include "GenericBCS.h"
20 #include "Vector.h"
21 #include "DiffusionDiscretization.h"
22 #include "ConvectionDiscretization.h"
23 #include "TimeDerivativeDiscretization.h"
24 #include "AMG.h"
25 #include "Linearizer.h"
26 #include "GradientModel.h"
27 #include "GenericIBDiscretization.h"
28 #include "SourceDiscretization.h"
29 #include "LinearizeInterfaceJump.h"
30 #include "BatteryLinearizeSpeciesInterface.h"
31 #include "BatteryLinearizePotentialInterface.h"
32 #include "BatteryLinearizeThermalInterface.h"
33 #include "BatteryPCLinearizeInterface_BV.h"
34 #include "BatteryBinaryElectrolyteDiscretization.h"
35 #include "BatteryPCDiffusionDiscretization.h"
36 #include "BatteryPC_BCS.h"
37 #include "BatteryPCTimeDerivativeDiscretization.h"
38 #include "BatteryPCBinaryElectrolyteDiscretization.h"
39 #include "BatteryPCHeatSourceDiscretization.h"
40 #include "BatteryFixInterfaceGhost.h"
41 #include "LinearizeInterfaceJumpUnconnected.h"
42 
43 template<class T>
44 class BatteryModel<T>::Impl
45 {
46 public:
47  typedef Array<T> TArray;
48  typedef Gradient<T> TGradType;
49  typedef Array<Gradient<T> > TGradArray;
50  typedef CRMatrix<T,T,T> T_Matrix;
51 
52  // PC Thermal included
53  typedef Vector<T,3> VectorT3;
54  typedef Array<VectorT3> VectorT3Array;
55  typedef SquareTensor<T,3> SquareTensorT3;
56  typedef Array<Gradient<VectorT3> > VectorT3GradArray;
57  typedef Gradient<VectorT3> VectorT3Grad;
58  typedef DiagonalTensor<T,3> DiagTensorT3;
59 
60  // PC No Thermal included
61  typedef Vector<T,2> VectorT2;
62  typedef Array<VectorT2> VectorT2Array;
63  typedef SquareTensor<T,2> SquareTensorT2;
64  typedef Array<Gradient<VectorT2> > VectorT2GradArray;
65  typedef Gradient<VectorT2> VectorT2Grad;
66  typedef DiagonalTensor<T,2> DiagTensorT2;
67 
68  Impl(GeomFields& geomFields,
69  const MeshList& realMeshes,
70  const MeshList& meshes,
71  const int nSpecies) :
72  _realMeshes(realMeshes),
73  _meshes(meshes),
74  _geomFields(geomFields),
75  _niters(0),
76  _nSpecies(nSpecies),
77  _batteryModelFields("batteryModel")
78  {
79 
80  const int numMeshes = _meshes.size();
81 
82  for (int n=0; n<numMeshes; n++)
83  {
84  const Mesh& mesh = *_meshes[n];
85 
86  BatteryPotentialVC<T> *pvc(new BatteryPotentialVC<T>());
87  pvc->vcType = "dielectric";
88  _pvcMap[mesh.getID()] = pvc;
89 
90  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
91  {
92  const FaceGroup& fg = *fgPtr;
93  BatteryPotentialBC<T> *pbc(new BatteryPotentialBC<T>());
94 
95  _pbcMap[fg.id] = pbc;
96  if (fg.groupType == "wall")
97  {
98  //pbc->bcType = "SpecifiedPotential";
99  pbc->bcType = "Symmetry";
100  }
101  else if (fg.groupType == "symmetry")
102  {
103  pbc->bcType = "Symmetry";
104  }
105  else
106  throw CException("ElectricModel: unknown face group type "
107  + fg.groupType);
108  }
109  }
110 
111  for (int n=0; n<numMeshes; n++)
112  {
113  const Mesh& mesh = *_meshes[n];
114  BatteryThermalVC<T> *tvc(new BatteryThermalVC<T>());
115  tvc->vcType = "flow";
116  _tvcMap[mesh.getID()] = tvc;
117 
118  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
119  {
120  const FaceGroup& fg = *fgPtr;
121  BatteryThermalBC<T> *tbc(new BatteryThermalBC<T>());
122 
123  _tbcMap[fg.id] = tbc;
124 
125  if ((fg.groupType == "wall") ||
126  (fg.groupType == "symmetry"))
127  {
128  //tbc->bcType = "SpecifiedHeatFlux";
129  tbc->bcType = "Symmetry";
130  }
131  else if ((fg.groupType == "velocity-inlet") ||
132  (fg.groupType == "pressure-outlet"))
133  {
134  tbc->bcType = "SpecifiedTemperature";
135  }
136  else
137  throw CException("ThermalModel: unknown face group type "
138  + fg.groupType);
139  }
140  }
141 
142  for (int m=0; m<_nSpecies; m++)
143  {
144  BatterySpeciesVCMap *svcmap = new BatterySpeciesVCMap();
145  _svcMapVector.push_back(svcmap);
146  BatterySpeciesBCMap *sbcmap = new BatterySpeciesBCMap();
147  _sbcMapVector.push_back(sbcmap);
148  BatterySpeciesFields *sFields = new BatterySpeciesFields("species");
149  _speciesFieldsVector.push_back(sFields);
150  MFRPtr *iNorm = new MFRPtr();
151  _initialSpeciesNormVector.push_back(iNorm);
152  MFRPtr *rCurrent = new MFRPtr();
153  _currentSpeciesResidual.push_back(rCurrent);
154 
155  for (int n=0; n<numMeshes; n++)
156  {
157 
158  const Mesh& mesh = *_meshes[n];
159  BatterySpeciesVC<T> *svc(new BatterySpeciesVC<T>());
160  svc->vcType = "flow";
161  (*svcmap)[mesh.getID()] = svc;
162 
163  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
164  {
165  const FaceGroup& fg = *fgPtr;
166  BatterySpeciesBC<T> *sbc(new BatterySpeciesBC<T>());
167 
168  (*sbcmap)[fg.id] = sbc;
169 
170  if ((fg.groupType == "wall") ||
171  (fg.groupType == "symmetry"))
172  {
173  //sbc->bcType = "SpecifiedMassFlux";
174  sbc->bcType = "Symmetry";
175  }
176  else
177  throw CException("SpeciesModel: unknown face group type "
178  + fg.groupType);
179  }
180  }
181  }
182  }
183 
184  void init()
185  {
186  const int numMeshes = _meshes.size();
187 
188  // initialize fields for battery model
189  for (int n=0; n<numMeshes; n++)
190  {
191  const Mesh& mesh = *_meshes[n];
192  const BatteryPotentialVC<T>& pvc = *_pvcMap[mesh.getID()];
193  const BatteryThermalVC<T>& tvc = *_tvcMap[mesh.getID()];
194  const StorageSite& cells = mesh.getCells();
195 
196  const int nCells = cells.getCount();
197 
198  // print cell centroids for test case
199  //const VectorT3Array& cellCentroid = dynamic_cast<const VectorT3Array&>(_geomFields.coordinate[cells]);
200  VectorT3Array& cellCentroid = dynamic_cast<VectorT3Array&>(_geomFields.coordinate[cells]);
201 
202  cout << "Mesh: " << n << endl;
203  for (int c=0; c<nCells; c++)
204  {
205  if (((cellCentroid[c])[2] < -3.33)&&((cellCentroid[c])[2] > -3.34))
206  cout << c << ": " << (cellCentroid[c])[0] << " " << (cellCentroid[c])[1] << endl;
207  }
208 
209  //initial potential setup
210  shared_ptr<TArray> pCell(new TArray(nCells));
211 
212  if (!mesh.isDoubleShell())
213  {
214  *pCell = pvc["initialPotential"];
215  }
216  else
217  {
218  // double shell mesh cells take on values of parents
219  typename BatteryPotentialVCMap::const_iterator posParent = _pvcMap.find(mesh.getParentMeshID());
220  typename BatteryPotentialVCMap::const_iterator posOther = _pvcMap.find(mesh.getOtherMeshID());
221  const BatteryPotentialVC<T>& pvcP = *(posParent->second);
222  const BatteryPotentialVC<T>& pvcO = *(posOther->second);
223 
224  const int NumberOfCells = cells.getCount();
225  for (int c=0; c<NumberOfCells; c++)
226  {
227  if ((c < NumberOfCells/4)||(c >= 3*NumberOfCells/4))
228  {
229  (*pCell)[c] = pvcP["initialPotential"];
230  }
231  else
232  {
233  (*pCell)[c] = pvcO["initialPotential"];
234  }
235  }
236  }
237 
238  _batteryModelFields.potential.addArray(cells,pCell);
239 
240  if (_options.transient)
241  {
242  // not actually discretized as an unsteady term
243  // only used to recover values from beginning of given time step
244  // in case a new technique needs to be tried after a solution failure
245  _batteryModelFields.potentialN1.addArray(cells, dynamic_pointer_cast<ArrayBase>(pCell->newCopy()));
246  }
247 
248  //initial temperature setup
249  shared_ptr<TArray> tCell(new TArray(nCells));
250 
251  if (!mesh.isDoubleShell())
252  {
253  *tCell = tvc["initialTemperature"];
254  }
255  else
256  {
257  // double shell mesh cells take on values of parents
258  typename BatteryThermalVCMap::const_iterator posParent = _tvcMap.find(mesh.getParentMeshID());
259  typename BatteryThermalVCMap::const_iterator posOther = _tvcMap.find(mesh.getOtherMeshID());
260  const BatteryThermalVC<T>& tvcP = *(posParent->second);
261  const BatteryThermalVC<T>& tvcO = *(posOther->second);
262 
263  const int NumberOfCells = cells.getCount();
264  for (int c=0; c<NumberOfCells; c++)
265  {
266  if ((c < NumberOfCells/4)||(c >= 3*NumberOfCells/4))
267  {
268  (*tCell)[c] = tvcP["initialTemperature"];
269  }
270  else
271  {
272  (*tCell)[c] = tvcO["initialTemperature"];
273  }
274  }
275  }
276 
277  _batteryModelFields.temperature.addArray(cells,tCell);
278 
279  if (_options.transient)
280  {
281  _batteryModelFields.temperatureN1.addArray(cells, dynamic_pointer_cast<ArrayBase>(tCell->newCopy()));
282  if (_options.timeDiscretizationOrder > 1)
283  _batteryModelFields.temperatureN2.addArray(cells, dynamic_pointer_cast<ArrayBase>(tCell->newCopy()));
284  }
285 
286  //electric/ionic conductivity setup
287  shared_ptr<TArray> condCell(new TArray(nCells));
288  *condCell = pvc["conductivity"];
289  _batteryModelFields.conductivity.addArray(cells,condCell);
290 
291  //thermal conductivity setup
292  shared_ptr<TArray> thermCondCell(new TArray(nCells));
293  *thermCondCell = tvc["thermalConductivity"];
294  _batteryModelFields.thermalConductivity.addArray(cells,thermCondCell);
295 
296  // species gradient setup
297  shared_ptr<TGradArray> gradS(new TGradArray(cells.getCount()));
298  gradS->zero();
299  _batteryModelFields.speciesGradient.addArray(cells,gradS);
300 
301  //potential gradient setup
302  shared_ptr<TGradArray> gradp(new TGradArray(nCells));
303  gradp->zero();
304  _batteryModelFields.potential_gradient.addArray(cells,gradp);
305 
306  //temperature gradient setup
307  shared_ptr<TGradArray> gradt(new TGradArray(nCells));
308  gradt->zero();
309  _batteryModelFields.temperatureGradient.addArray(cells,gradt);
310 
311  //create rho *specific heat field rho*Cp
312  shared_ptr<TArray> rhoCp(new TArray(nCells));
313  *rhoCp = tvc["density"] * tvc["specificHeat"];
314  _batteryModelFields.rhoCp.addArray(cells, rhoCp);
315 
316  //heat source
317  shared_ptr<TArray> sCell(new TArray(nCells));
318  *sCell = T(0.);
319  _batteryModelFields.heatSource.addArray(cells,sCell);
320 
321  //initial potential flux; Note: potential_flux only stored on boundary faces
322  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
323  {
324  const FaceGroup& fg = *fgPtr;
325  const StorageSite& faces = fg.site;
326 
327  shared_ptr<TArray> pFlux(new TArray(faces.getCount()));
328  pFlux->zero();
329  _batteryModelFields.potential_flux.addArray(faces,pFlux);
330 
331  }
332  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
333  {
334  const FaceGroup& fg = *fgPtr;
335  const StorageSite& faces = fg.site;
336 
337  shared_ptr<TArray> pFlux(new TArray(faces.getCount()));
338  pFlux->zero();
339  _batteryModelFields.potential_flux.addArray(faces,pFlux);
340 
341  }
342 
343  //initial heat flux; Note: heatflux only stored on boundary faces
344  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
345  {
346  const FaceGroup& fg = *fgPtr;
347  const StorageSite& faces = fg.site;
348 
349  shared_ptr<TArray> tFlux(new TArray(faces.getCount()));
350  tFlux->zero();
351  _batteryModelFields.heatFlux.addArray(faces,tFlux);
352 
353  }
354  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
355  {
356  const FaceGroup& fg = *fgPtr;
357  const StorageSite& faces = fg.site;
358 
359  shared_ptr<TArray> tFlux(new TArray(faces.getCount()));
360  tFlux->zero();
361  _batteryModelFields.heatFlux.addArray(faces,tFlux);
362 
363  }
364 
365  shared_ptr<TArray> lnLCCell(new TArray(nCells));
366  lnLCCell->zero();
367  _batteryModelFields.lnLithiumConcentration.addArray(cells,lnLCCell);
368 
369  shared_ptr<TGradArray> gradLnLC(new TGradArray(nCells));
370  gradLnLC->zero();
371  _batteryModelFields.lnLithiumConcentrationGradient.addArray(cells,gradLnLC);
372 
373  //set up combined fields for point-coupled solve
374  if (_options.thermalModelPC)
375  {
376  shared_ptr<VectorT3Array> pstCell(new VectorT3Array(nCells));
377  pstCell->zero();
378  _batteryModelFields.potentialSpeciesTemp.addArray(cells,pstCell);
379 
380  if (_options.transient)
381  {
382  _batteryModelFields.potentialSpeciesTempN1.addArray(cells,
383  dynamic_pointer_cast<ArrayBase>(pstCell->newCopy()));
384  if (_options.timeDiscretizationOrder > 1)
385  {
386  //throw CException("BatteryModel: Time Discretization Order > 1 not implimented.");
387  _batteryModelFields.potentialSpeciesTempN2.addArray(cells,
388  dynamic_pointer_cast<ArrayBase>(pstCell->newCopy()));
389  }
390 
391  }
392 
393  //combined diffusivity field
394  shared_ptr<VectorT3Array> pstDiffCell(new VectorT3Array(nCells));
395  pstDiffCell->zero();
396  _batteryModelFields.potentialSpeciesTempDiffusivity.addArray(cells,pstDiffCell);
397 
398  //combined gradient setup
399  shared_ptr<VectorT3GradArray> gradpst(new VectorT3GradArray(nCells));
400  gradpst->zero();
401  _batteryModelFields.potentialSpeciesTempGradient.addArray(cells,gradpst);
402 
403  //initial combined flux flied; Note: flux only stored on boundary faces
404  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
405  {
406  const FaceGroup& fg = *fgPtr;
407  const StorageSite& faces = fg.site;
408 
409  shared_ptr<VectorT3Array> pstFlux(new VectorT3Array(faces.getCount()));
410  pstFlux->zero();
411  _batteryModelFields.potentialSpeciesTempFlux.addArray(faces,pstFlux);
412 
413  }
414  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
415  {
416  const FaceGroup& fg = *fgPtr;
417  const StorageSite& faces = fg.site;
418 
419  shared_ptr<VectorT3Array> pstFlux(new VectorT3Array(faces.getCount()));
420  pstFlux->zero();
421  _batteryModelFields.potentialSpeciesTempFlux.addArray(faces,pstFlux);
422  }
423  }
424  else
425  {
426  shared_ptr<VectorT2Array> pstCell(new VectorT2Array(nCells));
427  pstCell->zero();
428  _batteryModelFields.potentialSpeciesTemp.addArray(cells,pstCell);
429 
430  if (_options.transient)
431  {
432  _batteryModelFields.potentialSpeciesTempN1.addArray(cells,
433  dynamic_pointer_cast<ArrayBase>(pstCell->newCopy()));
434  if (_options.timeDiscretizationOrder > 1)
435  {
436  //throw CException("BatteryModel: Time Discretization Order > 1 not implimented.");
437  _batteryModelFields.potentialSpeciesTempN2.addArray(cells,
438  dynamic_pointer_cast<ArrayBase>(pstCell->newCopy()));
439  }
440 
441  }
442 
443  //combined diffusivity field
444  shared_ptr<VectorT2Array> pstDiffCell(new VectorT2Array(nCells));
445  pstDiffCell->zero();
446  _batteryModelFields.potentialSpeciesTempDiffusivity.addArray(cells,pstDiffCell);
447 
448  //combined gradient setup
449  shared_ptr<VectorT2GradArray> gradpst(new VectorT2GradArray(nCells));
450  gradpst->zero();
451  _batteryModelFields.potentialSpeciesTempGradient.addArray(cells,gradpst);
452 
453  //initial combined flux flied; Note: flux only stored on boundary faces
454  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
455  {
456  const FaceGroup& fg = *fgPtr;
457  const StorageSite& faces = fg.site;
458 
459  shared_ptr<VectorT2Array> pstFlux(new VectorT2Array(faces.getCount()));
460  pstFlux->zero();
461  _batteryModelFields.potentialSpeciesTempFlux.addArray(faces,pstFlux);
462 
463  }
464  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
465  {
466  const FaceGroup& fg = *fgPtr;
467  const StorageSite& faces = fg.site;
468 
469  shared_ptr<VectorT2Array> pstFlux(new VectorT2Array(faces.getCount()));
470  pstFlux->zero();
471  _batteryModelFields.potentialSpeciesTempFlux.addArray(faces,pstFlux);
472  }
473  }
474  }
475 
476  for (int m=0; m<_nSpecies; m++)
477  {
478  const BatterySpeciesVCMap& svcmap = *_svcMapVector[m];
479  BatterySpeciesFields& sFields = *_speciesFieldsVector[m];
480  //MFRPtr& iNorm = *_initialNormVector[m];
481 
482  for (int n=0; n<numMeshes; n++)
483  {
484  const Mesh& mesh = *_meshes[n];
485 
486  const StorageSite& cells = mesh.getCells();
487  const StorageSite& allFaces = mesh.getFaces();
488 
489  typename BatterySpeciesVCMap::const_iterator pos = svcmap.find(mesh.getID());
490  if (pos==svcmap.end())
491  {
492  throw CException("BatteryModel: Error in Species VC Map");
493  }
494 
495  const BatterySpeciesVC<T>& svc = *(pos->second);
496 
497 
498  //concentration
499 
500  shared_ptr<TArray> concCell(new TArray(cells.getCount()));
501 
502  if (!mesh.isDoubleShell())
503  {
504  *concCell = svc["initialConcentration"];
505  }
506  else
507  {
508  // double shell mesh cells take on values of parents
509  typename BatterySpeciesVCMap::const_iterator posParent = svcmap.find(mesh.getParentMeshID());
510  typename BatterySpeciesVCMap::const_iterator posOther = svcmap.find(mesh.getOtherMeshID());
511  const BatterySpeciesVC<T>& svcP = *(posParent->second);
512  const BatterySpeciesVC<T>& svcO = *(posOther->second);
513 
514  const int NumberOfCells = cells.getCount();
515  for (int c=0; c<NumberOfCells; c++)
516  {
517  if ((c < NumberOfCells/4)||(c >= 3*NumberOfCells/4))
518  {
519  (*concCell)[c] = svcP["initialConcentration"];
520  }
521  else
522  {
523  (*concCell)[c] = svcO["initialConcentration"];
524  }
525  }
526  }
527 
528  sFields.concentration.addArray(cells,concCell);
529 
530  if (_options.transient)
531  {
532  sFields.concentrationN1.addArray(cells,
533  dynamic_pointer_cast<ArrayBase>(concCell->newCopy()));
534  if (_options.timeDiscretizationOrder > 1)
535  sFields.concentrationN2.addArray(cells,
536  dynamic_pointer_cast<ArrayBase>(concCell->newCopy()));
537 
538  }
539 
540  //diffusivity
541  shared_ptr<TArray> diffusivityCell(new TArray(cells.getCount()));
542  *diffusivityCell = svc["massDiffusivity"];
543  sFields.diffusivity.addArray(cells,diffusivityCell);
544 
545  //create a one field
546  shared_ptr<TArray> oneCell(new TArray(cells.getCount()));
547  *oneCell = T(1.0);
548  sFields.one.addArray(cells,oneCell);
549 
550  //inital convection flux at faces
551  shared_ptr<TArray> convFlux(new TArray(allFaces.getCount()));
552  convFlux->zero();
553  sFields.convectionFlux.addArray(allFaces,convFlux);
554 
555  //mass flux at faces
556  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
557  {
558  const FaceGroup& fg = *fgPtr;
559  const StorageSite& faces = fg.site;
560 
561  shared_ptr<TArray> fluxFace(new TArray(faces.getCount()));
562 
563  fluxFace->zero();
564  sFields.massFlux.addArray(faces,fluxFace);
565 
566  }
567  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
568  {
569  const FaceGroup& fg = *fgPtr;
570  const StorageSite& faces = fg.site;
571 
572  shared_ptr<TArray> fluxFace(new TArray(faces.getCount()));
573 
574  fluxFace->zero();
575  sFields.massFlux.addArray(faces,fluxFace);
576 
577  }
578 
579 
580  }
581 
582  sFields.diffusivity.syncLocal();
583  }
584  _niters =0;
585  _initialPotentialNorm = MFRPtr();
586  _currentPotentialResidual = MFRPtr();
587  _initialThermalNorm = MFRPtr();
588  _currentThermalResidual = MFRPtr();
589  _initialPCNorm = MFRPtr();
590  _currentPCResidual = MFRPtr();
591  _batteryModelFields.conductivity.syncLocal();
592  _batteryModelFields.thermalConductivity.syncLocal();
593  copyPCDiffusivity();
594  _batteryModelFields.potentialSpeciesTempDiffusivity.syncLocal();
595  }
596 
597  BatterySpeciesFields& getBatterySpeciesFields(const int speciesId) {return *_speciesFieldsVector[speciesId];}
598  BatteryModelFields& getBatteryModelFields() {return _batteryModelFields;}
599  BatterySpeciesVCMap& getSpeciesVCMap(const int speciesId) {return *_svcMapVector[speciesId];}
600  BatterySpeciesBCMap& getSpeciesBCMap(const int speciesId) {return *_sbcMapVector[speciesId];}
601  BatteryPotentialVCMap& getPotentialVCMap() {return _pvcMap;}
602  BatteryPotentialBCMap& getPotentialBCMap() {return _pbcMap;}
603  BatteryThermalVCMap& getThermalVCMap() {return _tvcMap;}
604  BatteryThermalBCMap& getThermalBCMap() {return _tbcMap;}
605 
606  BatteryModelOptions<T>& getOptions() {return _options;}
607 
608  void updateTime()
609  {
610  const int numMeshes = _meshes.size();
611 
612  for (int m=0; m<_nSpecies; m++)
613  {
614  BatterySpeciesFields& sFields = *_speciesFieldsVector[m];
615 
616  for (int n=0; n<numMeshes; n++)
617  {
618  const Mesh& mesh = *_meshes[n];
619 
620  const StorageSite& cells = mesh.getCells();
621  TArray& mF =
622  dynamic_cast<TArray&>(sFields.concentration[cells]);
623  TArray& mFN1 =
624  dynamic_cast<TArray&>(sFields.concentrationN1[cells]);
625 
626  if (_options.timeDiscretizationOrder > 1)
627  {
628  TArray& mFN2 =
629  dynamic_cast<TArray&>(sFields.concentrationN2[cells]);
630  mFN2 = mFN1;
631  }
632  mFN1 = mF;
633  }
634  }
635 
636  for (int n=0; n<numMeshes; n++)
637  {
638  const Mesh& mesh = *_meshes[n];
639  const StorageSite& cells = mesh.getCells();
640 
641  if (_options.thermalModelPC)
642  {
643  VectorT3Array& pST =
644  dynamic_cast<VectorT3Array&>(_batteryModelFields.potentialSpeciesTemp[cells]);
645  VectorT3Array& pSTN1 =
646  dynamic_cast<VectorT3Array&>(_batteryModelFields.potentialSpeciesTempN1[cells]);
647 
648  if (_options.timeDiscretizationOrder > 1)
649  {
650  //throw CException("BatteryModel: Time Discretization Order > 1 not implimented.");
651  VectorT3Array& pSTN2 = dynamic_cast<VectorT3Array&>(_batteryModelFields.potentialSpeciesTempN2[cells]);
652  pSTN2 = pSTN1;
653  }
654  pSTN1 = pST;
655  }
656  else
657  {
658  VectorT2Array& pST =
659  dynamic_cast<VectorT2Array&>(_batteryModelFields.potentialSpeciesTemp[cells]);
660  VectorT2Array& pSTN1 =
661  dynamic_cast<VectorT2Array&>(_batteryModelFields.potentialSpeciesTempN1[cells]);
662 
663  if (_options.timeDiscretizationOrder > 1)
664  {
665  //throw CException("BatteryModel: Time Discretization Order > 1 not implimented.");
666  VectorT2Array& pSTN2 = dynamic_cast<VectorT2Array&>(_batteryModelFields.potentialSpeciesTempN2[cells]);
667  pSTN2 = pSTN1;
668  }
669  pSTN1 = pST;
670  }
671 
672  TArray& temp =
673  dynamic_cast<TArray&>(_batteryModelFields.temperature[cells]);
674  TArray& tempN1 =
675  dynamic_cast<TArray&>(_batteryModelFields.temperatureN1[cells]);
676 
677  if (_options.timeDiscretizationOrder > 1)
678  {
679  TArray& tempN2 =
680  dynamic_cast<TArray&>(_batteryModelFields.temperatureN2[cells]);
681  tempN2 = tempN1;
682  }
683  tempN1 = temp;
684 
685  // only used for recoverLastTimeStep, no unsteady term in potential equation
686  TArray& potential =
687  dynamic_cast<TArray&>(_batteryModelFields.potential[cells]);
688  TArray& potentialN1 =
689  dynamic_cast<TArray&>(_batteryModelFields.potentialN1[cells]);
690  potentialN1 = potential;
691 
692  }
693  }
694 
695 void recoverLastTimestep()
696  {
697  const int numMeshes = _meshes.size();
698 
699  for (int m=0; m<_nSpecies; m++)
700  {
701  BatterySpeciesFields& sFields = *_speciesFieldsVector[m];
702 
703  for (int n=0; n<numMeshes; n++)
704  {
705  const Mesh& mesh = *_meshes[n];
706 
707  const StorageSite& cells = mesh.getCells();
708  TArray& mF =
709  dynamic_cast<TArray&>(sFields.concentration[cells]);
710  TArray& mFN1 =
711  dynamic_cast<TArray&>(sFields.concentrationN1[cells]);
712 
713  mF = mFN1;
714 
715  /*if (_options.timeDiscretizationOrder > 1)
716  {
717  TArray& mFN2 =
718  dynamic_cast<TArray&>(sFields.concentrationN2[cells]);
719  mFN1 = mFN2;
720  }*/
721  }
722  }
723 
724  for (int n=0; n<numMeshes; n++)
725  {
726  const Mesh& mesh = *_meshes[n];
727  const StorageSite& cells = mesh.getCells();
728 
729  if (_options.thermalModelPC)
730  {
731  VectorT3Array& pST =
732  dynamic_cast<VectorT3Array&>(_batteryModelFields.potentialSpeciesTemp[cells]);
733  VectorT3Array& pSTN1 =
734  dynamic_cast<VectorT3Array&>(_batteryModelFields.potentialSpeciesTempN1[cells]);
735  pST = pSTN1;
736 
737  /*if (_options.timeDiscretizationOrder > 1)
738  {
739  VectorT3Array& pSTN2 = dynamic_cast<VectorT3Array&>(_batteryModelFields.potentialSpeciesTempN2[cells]);
740  pSTN1 = pSTN2;
741  }*/
742  }
743  else
744  {
745  VectorT2Array& pST =
746  dynamic_cast<VectorT2Array&>(_batteryModelFields.potentialSpeciesTemp[cells]);
747  VectorT2Array& pSTN1 =
748  dynamic_cast<VectorT2Array&>(_batteryModelFields.potentialSpeciesTempN1[cells]);
749  pST = pSTN1;
750  /*if (_options.timeDiscretizationOrder > 1)
751  {
752  VectorT2Array& pSTN2 = dynamic_cast<VectorT2Array&>(_batteryModelFields.potentialSpeciesTempN2[cells]);
753  pSTN1 = pSTN2;
754  }*/
755  }
756 
757  TArray& temp =
758  dynamic_cast<TArray&>(_batteryModelFields.temperature[cells]);
759  TArray& tempN1 =
760  dynamic_cast<TArray&>(_batteryModelFields.temperatureN1[cells]);
761  temp = tempN1;
762 
763  /*if (_options.timeDiscretizationOrder > 1)
764  {
765  TArray& tempN2 =
766  dynamic_cast<TArray&>(_batteryModelFields.temperatureN2[cells]);
767  tempN1 = tempN2;
768  }*/
769 
770  TArray& potential =
771  dynamic_cast<TArray&>(_batteryModelFields.potential[cells]);
772  TArray& potentialN1 =
773  dynamic_cast<TArray&>(_batteryModelFields.potentialN1[cells]);
774  potential = potentialN1;
775  }
776  }
777 
778  void initSpeciesLinearization(LinearSystem& ls, LinearSystem& lsShell, const int& SpeciesNumber)
779  {
780  BatterySpeciesFields& sFields = *_speciesFieldsVector[SpeciesNumber];
781  const int numMeshes = _meshes.size();
782 
783  for (int n=0; n<numMeshes; n++)
784  {
785  const Mesh& mesh = *_meshes[n];
786 
787  if ((!(mesh.isDoubleShell()))||((mesh.isDoubleShell())&&(mesh.isConnectedShell())))
788  {
789  const StorageSite& cells = mesh.getCells();
790  MultiField::ArrayIndex tIndex(&sFields.concentration,&cells);
791 
792  ls.getX().addArray(tIndex,sFields.concentration.getArrayPtr(cells));
793 
794  const CRConnectivity& cellCells = mesh.getCellCells();
795 
796  shared_ptr<Matrix> m(new CRMatrix<T,T,T>(cellCells));
797 
798  ls.getMatrix().addMatrix(tIndex,tIndex,m);
799  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
800  {
801  const FaceGroup& fg = *fgPtr;
802  const StorageSite& faces = fg.site;
803 
804  MultiField::ArrayIndex fIndex(&sFields.massFlux,&faces);
805  ls.getX().addArray(fIndex,sFields.massFlux.getArrayPtr(faces));
806 
807  const CRConnectivity& faceCells = mesh.getFaceCells(faces);
808 
809  shared_ptr<Matrix> mft(new FluxJacobianMatrix<T,T>(faceCells));
810  ls.getMatrix().addMatrix(fIndex,tIndex,mft);
811 
812  shared_ptr<Matrix> mff(new DiagonalMatrix<T,T>(faces.getCount()));
813  ls.getMatrix().addMatrix(fIndex,fIndex,mff);
814  }
815 
816  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
817  {
818  const FaceGroup& fg = *fgPtr;
819  const StorageSite& faces = fg.site;
820 
821  MultiField::ArrayIndex fIndex(&sFields.massFlux,&faces);
822  ls.getX().addArray(fIndex,sFields.massFlux.getArrayPtr(faces));
823 
824  const CRConnectivity& faceCells = mesh.getFaceCells(faces);
825 
826  shared_ptr<Matrix> mft(new FluxJacobianMatrix<T,T>(faceCells));
827  ls.getMatrix().addMatrix(fIndex,tIndex,mft);
828 
829  shared_ptr<Matrix> mff(new DiagonalMatrix<T,T>(faces.getCount()));
830  ls.getMatrix().addMatrix(fIndex,fIndex,mff);
831  }
832  }
833  else
834  {
835  // must be unconnected shell, add to shell LS
836  const StorageSite& cells = mesh.getCells();
837  MultiField::ArrayIndex tIndex(&sFields.concentration,&cells);
838 
839  lsShell.getX().addArray(tIndex,sFields.concentration.getArrayPtr(cells));
840 
841  const CRConnectivity& cellCells = mesh.getCellCells();
842 
843  shared_ptr<Matrix> m(new CRMatrix<T,T,T>(cellCells));
844 
845  lsShell.getMatrix().addMatrix(tIndex,tIndex,m);
846  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
847  {
848  const FaceGroup& fg = *fgPtr;
849  const StorageSite& faces = fg.site;
850 
851  MultiField::ArrayIndex fIndex(&sFields.massFlux,&faces);
852  lsShell.getX().addArray(fIndex,sFields.massFlux.getArrayPtr(faces));
853 
854  const CRConnectivity& faceCells = mesh.getFaceCells(faces);
855 
856  shared_ptr<Matrix> mft(new FluxJacobianMatrix<T,T>(faceCells));
857  lsShell.getMatrix().addMatrix(fIndex,tIndex,mft);
858 
859  shared_ptr<Matrix> mff(new DiagonalMatrix<T,T>(faces.getCount()));
860  lsShell.getMatrix().addMatrix(fIndex,fIndex,mff);
861  }
862 
863  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
864  {
865  const FaceGroup& fg = *fgPtr;
866  const StorageSite& faces = fg.site;
867 
868  MultiField::ArrayIndex fIndex(&sFields.massFlux,&faces);
869  lsShell.getX().addArray(fIndex,sFields.massFlux.getArrayPtr(faces));
870 
871  const CRConnectivity& faceCells = mesh.getFaceCells(faces);
872 
873  shared_ptr<Matrix> mft(new FluxJacobianMatrix<T,T>(faceCells));
874  lsShell.getMatrix().addMatrix(fIndex,tIndex,mft);
875 
876  shared_ptr<Matrix> mff(new DiagonalMatrix<T,T>(faces.getCount()));
877  lsShell.getMatrix().addMatrix(fIndex,fIndex,mff);
878  }
879  }
880  }
881  }
882 
883  void initPCLinearization(LinearSystem& ls)
884  {
885  const int numMeshes = _meshes.size();
886 
887  for (int n=0; n<numMeshes; n++)
888  {
889  const Mesh& mesh = *_meshes[n];
890 
891  if (mesh.isDoubleShell())
892  {
893 
894  const StorageSite& cells = mesh.getCells();
895  MultiField::ArrayIndex tIndex(&_batteryModelFields.potentialSpeciesTemp,&cells);
896 
897  ls.getX().addArray(tIndex,_batteryModelFields.potentialSpeciesTemp.getArrayPtr(cells));
898 
899  const CRConnectivity& cellCells = mesh.getCellCells();
900 
901  if (_options.thermalModelPC)
902  {
903  shared_ptr<Matrix> m(new CRMatrix<SquareTensorT3,SquareTensorT3,VectorT3>(cellCells));
904 
905  ls.getMatrix().addMatrix(tIndex,tIndex,m);
906 
907  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
908  {
909  const FaceGroup& fg = *fgPtr;
910  const StorageSite& faces = fg.site;
911 
912  MultiField::ArrayIndex fIndex(&_batteryModelFields.potentialSpeciesTempFlux,&faces);
913  ls.getX().addArray(fIndex,_batteryModelFields.potentialSpeciesTempFlux.getArrayPtr(faces));
914 
915  const CRConnectivity& faceCells = mesh.getFaceCells(faces);
916 
917  shared_ptr<Matrix> mft(new FluxJacobianMatrix<SquareTensorT3,VectorT3>(faceCells));
918  ls.getMatrix().addMatrix(fIndex,tIndex,mft);
919 
920  shared_ptr<Matrix> mff(new DiagonalMatrix<SquareTensorT3,VectorT3>(faces.getCount()));
921  ls.getMatrix().addMatrix(fIndex,fIndex,mff);
922  }
923 
924  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
925  {
926  const FaceGroup& fg = *fgPtr;
927  const StorageSite& faces = fg.site;
928 
929  MultiField::ArrayIndex fIndex(&_batteryModelFields.potentialSpeciesTempFlux,&faces);
930  ls.getX().addArray(fIndex,_batteryModelFields.potentialSpeciesTempFlux.getArrayPtr(faces));
931 
932  const CRConnectivity& faceCells = mesh.getFaceCells(faces);
933 
934  shared_ptr<Matrix> mft(new FluxJacobianMatrix<SquareTensorT3,VectorT3>(faceCells));
935  ls.getMatrix().addMatrix(fIndex,tIndex,mft);
936 
937  shared_ptr<Matrix> mff(new DiagonalMatrix<SquareTensorT3,VectorT3>(faces.getCount()));
938  ls.getMatrix().addMatrix(fIndex,fIndex,mff);
939  }
940  }
941  else
942  {
943  shared_ptr<Matrix> m(new CRMatrix<SquareTensorT2,SquareTensorT2,VectorT2>(cellCells));
944 
945  ls.getMatrix().addMatrix(tIndex,tIndex,m);
946 
947  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
948  {
949  const FaceGroup& fg = *fgPtr;
950  const StorageSite& faces = fg.site;
951 
952  MultiField::ArrayIndex fIndex(&_batteryModelFields.potentialSpeciesTempFlux,&faces);
953  ls.getX().addArray(fIndex,_batteryModelFields.potentialSpeciesTempFlux.getArrayPtr(faces));
954 
955  const CRConnectivity& faceCells = mesh.getFaceCells(faces);
956 
957  shared_ptr<Matrix> mft(new FluxJacobianMatrix<SquareTensorT2,VectorT2>(faceCells));
958  ls.getMatrix().addMatrix(fIndex,tIndex,mft);
959 
960  shared_ptr<Matrix> mff(new DiagonalMatrix<SquareTensorT2,VectorT2>(faces.getCount()));
961  ls.getMatrix().addMatrix(fIndex,fIndex,mff);
962  }
963 
964  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
965  {
966  const FaceGroup& fg = *fgPtr;
967  const StorageSite& faces = fg.site;
968 
969  MultiField::ArrayIndex fIndex(&_batteryModelFields.potentialSpeciesTempFlux,&faces);
970  ls.getX().addArray(fIndex,_batteryModelFields.potentialSpeciesTempFlux.getArrayPtr(faces));
971 
972  const CRConnectivity& faceCells = mesh.getFaceCells(faces);
973 
974  shared_ptr<Matrix> mft(new FluxJacobianMatrix<SquareTensorT2,VectorT2>(faceCells));
975  ls.getMatrix().addMatrix(fIndex,tIndex,mft);
976 
977  shared_ptr<Matrix> mff(new DiagonalMatrix<SquareTensorT2,VectorT2>(faces.getCount()));
978  ls.getMatrix().addMatrix(fIndex,fIndex,mff);
979  }
980  }
981  }
982  else
983  {
984  const StorageSite& cells = mesh.getCells();
985  MultiField::ArrayIndex tIndex(&_batteryModelFields.potentialSpeciesTemp,&cells);
986 
987  ls.getX().addArray(tIndex,_batteryModelFields.potentialSpeciesTemp.getArrayPtr(cells));
988 
989  const CRConnectivity& cellCells = mesh.getCellCells();
990 
991  if (_options.thermalModelPC)
992  {
993  shared_ptr<Matrix> m(new CRMatrix<DiagTensorT3,DiagTensorT3,VectorT3>(cellCells));
994 
995  ls.getMatrix().addMatrix(tIndex,tIndex,m);
996 
997  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
998  {
999  const FaceGroup& fg = *fgPtr;
1000  const StorageSite& faces = fg.site;
1001 
1002  MultiField::ArrayIndex fIndex(&_batteryModelFields.potentialSpeciesTempFlux,&faces);
1003  ls.getX().addArray(fIndex,_batteryModelFields.potentialSpeciesTempFlux.getArrayPtr(faces));
1004 
1005  const CRConnectivity& faceCells = mesh.getFaceCells(faces);
1006 
1007  shared_ptr<Matrix> mft(new FluxJacobianMatrix<DiagTensorT3,VectorT3>(faceCells));
1008  ls.getMatrix().addMatrix(fIndex,tIndex,mft);
1009 
1010  shared_ptr<Matrix> mff(new DiagonalMatrix<DiagTensorT3,VectorT3>(faces.getCount()));
1011  ls.getMatrix().addMatrix(fIndex,fIndex,mff);
1012  }
1013 
1014  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
1015  {
1016  const FaceGroup& fg = *fgPtr;
1017  const StorageSite& faces = fg.site;
1018 
1019  MultiField::ArrayIndex fIndex(&_batteryModelFields.potentialSpeciesTempFlux,&faces);
1020  ls.getX().addArray(fIndex,_batteryModelFields.potentialSpeciesTempFlux.getArrayPtr(faces));
1021 
1022  const CRConnectivity& faceCells = mesh.getFaceCells(faces);
1023 
1024  shared_ptr<Matrix> mft(new FluxJacobianMatrix<DiagTensorT3,VectorT3>(faceCells));
1025  ls.getMatrix().addMatrix(fIndex,tIndex,mft);
1026 
1027  shared_ptr<Matrix> mff(new DiagonalMatrix<DiagTensorT3,VectorT3>(faces.getCount()));
1028  ls.getMatrix().addMatrix(fIndex,fIndex,mff);
1029  }
1030  }
1031  else
1032  {
1033  shared_ptr<Matrix> m(new CRMatrix<DiagTensorT2,DiagTensorT2,VectorT2>(cellCells));
1034 
1035  ls.getMatrix().addMatrix(tIndex,tIndex,m);
1036 
1037  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
1038  {
1039  const FaceGroup& fg = *fgPtr;
1040  const StorageSite& faces = fg.site;
1041 
1042  MultiField::ArrayIndex fIndex(&_batteryModelFields.potentialSpeciesTempFlux,&faces);
1043  ls.getX().addArray(fIndex,_batteryModelFields.potentialSpeciesTempFlux.getArrayPtr(faces));
1044 
1045  const CRConnectivity& faceCells = mesh.getFaceCells(faces);
1046 
1047  shared_ptr<Matrix> mft(new FluxJacobianMatrix<DiagTensorT2,VectorT2>(faceCells));
1048  ls.getMatrix().addMatrix(fIndex,tIndex,mft);
1049 
1050  shared_ptr<Matrix> mff(new DiagonalMatrix<DiagTensorT2,VectorT2>(faces.getCount()));
1051  ls.getMatrix().addMatrix(fIndex,fIndex,mff);
1052  }
1053 
1054  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
1055  {
1056  const FaceGroup& fg = *fgPtr;
1057  const StorageSite& faces = fg.site;
1058 
1059  MultiField::ArrayIndex fIndex(&_batteryModelFields.potentialSpeciesTempFlux,&faces);
1060  ls.getX().addArray(fIndex,_batteryModelFields.potentialSpeciesTempFlux.getArrayPtr(faces));
1061 
1062  const CRConnectivity& faceCells = mesh.getFaceCells(faces);
1063 
1064  shared_ptr<Matrix> mft(new FluxJacobianMatrix<DiagTensorT2,VectorT2>(faceCells));
1065  ls.getMatrix().addMatrix(fIndex,tIndex,mft);
1066 
1067  shared_ptr<Matrix> mff(new DiagonalMatrix<DiagTensorT2,VectorT2>(faces.getCount()));
1068  ls.getMatrix().addMatrix(fIndex,fIndex,mff);
1069  }
1070  }
1071  }
1072 
1073  }
1074  }
1075 
1076  void initPotentialLinearization(LinearSystem& ls)
1077  {
1078  const int numMeshes = _meshes.size();
1079  for (int n=0; n<numMeshes; n++)
1080  {
1081  const Mesh& mesh = *_meshes[n];
1082 
1083  const StorageSite& cells = mesh.getCells();
1084  MultiField::ArrayIndex tIndex(&_batteryModelFields.potential,&cells);
1085 
1086  ls.getX().addArray(tIndex,_batteryModelFields.potential.getArrayPtr(cells));
1087 
1088  const CRConnectivity& cellCells = mesh.getCellCells();
1089 
1090  shared_ptr<Matrix> m(new CRMatrix<T,T,T>(cellCells));
1091 
1092  ls.getMatrix().addMatrix(tIndex,tIndex,m);
1093 
1094  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
1095  {
1096  const FaceGroup& fg = *fgPtr;
1097  const StorageSite& faces = fg.site;
1098 
1099  MultiField::ArrayIndex fIndex(&_batteryModelFields.potential_flux,&faces);
1100  ls.getX().addArray(fIndex,_batteryModelFields.potential_flux.getArrayPtr(faces));
1101 
1102  const CRConnectivity& faceCells = mesh.getFaceCells(faces);
1103 
1104  shared_ptr<Matrix> mft(new FluxJacobianMatrix<T,T>(faceCells));
1105  ls.getMatrix().addMatrix(fIndex,tIndex,mft);
1106 
1107  shared_ptr<Matrix> mff(new DiagonalMatrix<T,T>(faces.getCount()));
1108  ls.getMatrix().addMatrix(fIndex,fIndex,mff);
1109  }
1110 
1111  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
1112  {
1113  const FaceGroup& fg = *fgPtr;
1114  const StorageSite& faces = fg.site;
1115 
1116  MultiField::ArrayIndex fIndex(&_batteryModelFields.potential_flux,&faces);
1117  ls.getX().addArray(fIndex,_batteryModelFields.potential_flux.getArrayPtr(faces));
1118 
1119  const CRConnectivity& faceCells = mesh.getFaceCells(faces);
1120 
1121  shared_ptr<Matrix> mft(new FluxJacobianMatrix<T,T>(faceCells));
1122  ls.getMatrix().addMatrix(fIndex,tIndex,mft);
1123 
1124  shared_ptr<Matrix> mff(new DiagonalMatrix<T,T>(faces.getCount()));
1125  ls.getMatrix().addMatrix(fIndex,fIndex,mff);
1126  }
1127 
1128  }
1129  }
1130 
1131 void initThermalLinearization(LinearSystem& ls)
1132  {
1133  const int numMeshes = _meshes.size();
1134  for (int n=0; n<numMeshes; n++)
1135  {
1136  const Mesh& mesh = *_meshes[n];
1137 
1138  const StorageSite& cells = mesh.getCells();
1139  MultiField::ArrayIndex tIndex(&_batteryModelFields.temperature,&cells);
1140 
1141  ls.getX().addArray(tIndex,_batteryModelFields.temperature.getArrayPtr(cells));
1142 
1143  const CRConnectivity& cellCells = mesh.getCellCells();
1144 
1145  shared_ptr<Matrix> m(new CRMatrix<T,T,T>(cellCells));
1146 
1147  ls.getMatrix().addMatrix(tIndex,tIndex,m);
1148 
1149  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
1150  {
1151  const FaceGroup& fg = *fgPtr;
1152  const StorageSite& faces = fg.site;
1153 
1154  MultiField::ArrayIndex fIndex(&_batteryModelFields.heatFlux,&faces);
1155  ls.getX().addArray(fIndex,_batteryModelFields.heatFlux.getArrayPtr(faces));
1156 
1157  const CRConnectivity& faceCells = mesh.getFaceCells(faces);
1158 
1159  shared_ptr<Matrix> mft(new FluxJacobianMatrix<T,T>(faceCells));
1160  ls.getMatrix().addMatrix(fIndex,tIndex,mft);
1161 
1162  shared_ptr<Matrix> mff(new DiagonalMatrix<T,T>(faces.getCount()));
1163  ls.getMatrix().addMatrix(fIndex,fIndex,mff);
1164  }
1165 
1166  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
1167  {
1168  const FaceGroup& fg = *fgPtr;
1169  const StorageSite& faces = fg.site;
1170 
1171  MultiField::ArrayIndex fIndex(&_batteryModelFields.heatFlux,&faces);
1172  ls.getX().addArray(fIndex,_batteryModelFields.heatFlux.getArrayPtr(faces));
1173 
1174  const CRConnectivity& faceCells = mesh.getFaceCells(faces);
1175 
1176  shared_ptr<Matrix> mft(new FluxJacobianMatrix<T,T>(faceCells));
1177  ls.getMatrix().addMatrix(fIndex,tIndex,mft);
1178 
1179  shared_ptr<Matrix> mff(new DiagonalMatrix<T,T>(faces.getCount()));
1180  ls.getMatrix().addMatrix(fIndex,fIndex,mff);
1181  }
1182 
1183  }
1184  }
1185 
1186  void linearizeSpecies(LinearSystem& ls, LinearSystem& lsShell, const int& m)
1187  {
1188  const BatterySpeciesBCMap& sbcmap = *_sbcMapVector[m];
1189  BatterySpeciesFields& sFields = *_speciesFieldsVector[m];
1190 
1191  sFields.concentration.syncLocal();
1192 
1193  //fix interface if it is an unconnected shell mesh
1194  // (tested inside so that multiple interfaces can be handled)
1195  BatteryFixInterfaceGhost<T, T, T> fig (_meshes,
1196  _geomFields,
1197  sFields.concentration,
1198  sFields.diffusivity);
1199 
1200  fig.fixInterfaces();
1201 
1202 
1203  GradientModel<T> speciesGradientModel(_meshes,sFields.concentration,
1204  _batteryModelFields.speciesGradient,_geomFields);
1205  speciesGradientModel.compute();
1206 
1207  DiscrList discretizations;
1208 
1209  /*shared_ptr<Discretization>
1210  dd(new DiffusionDiscretization<T,T,T>
1211  (_meshes,_geomFields,
1212  sFields.concentration,
1213  sFields.diffusivity,
1214  _batteryModelFields.speciesGradient));
1215  discretizations.push_back(dd);
1216 
1217  if (_options.transient)
1218  {
1219  shared_ptr<Discretization>
1220  td(new TimeDerivativeDiscretization<T,T,T>
1221  (_meshes,_geomFields,
1222  sFields.concentration,
1223  sFields.concentrationN1,
1224  sFields.concentrationN2,
1225  sFields.one,
1226  _options["timeStep"]));
1227 
1228  discretizations.push_back(td);
1229  }
1230 
1231 
1232  Linearizer linearizer;
1233 
1234  linearizer.linearize(discretizations,_meshes,ls.getMatrix(),
1235  ls.getX(), ls.getB());
1236  */
1237 
1238  shared_ptr<Discretization>
1239  dd(new DiffusionDiscretization<T,T,T>
1240  (_realMeshes,_geomFields,
1241  sFields.concentration,
1242  sFields.diffusivity,
1243  _batteryModelFields.speciesGradient));
1244  discretizations.push_back(dd);
1245 
1246  if (_options.transient)
1247  {
1248  shared_ptr<Discretization>
1249  td(new TimeDerivativeDiscretization<T,T,T>
1250  (_realMeshes,_geomFields,
1251  sFields.concentration,
1252  sFields.concentrationN1,
1253  sFields.concentrationN2,
1254  sFields.one,
1255  _options["timeStep"]));
1256 
1257  discretizations.push_back(td);
1258  }
1259 
1260 
1261  Linearizer linearizer;
1262 
1263  linearizer.linearize(discretizations,_realMeshes,ls.getMatrix(),
1264  ls.getX(), ls.getB());
1265 
1266  const int numMeshes = _meshes.size();
1267 
1268  // linearize shell mesh
1269 
1270  for (int n=0; n<numMeshes; n++)
1271  {
1272  const Mesh& mesh = *_meshes[n];
1273  if ((mesh.isDoubleShell())&&(mesh.isConnectedShell()))
1274  {
1275  const int parentMeshID = mesh.getParentMeshID();
1276  const int otherMeshID = mesh.getOtherMeshID();
1277  const Mesh& parentMesh = *_meshes[parentMeshID];
1278  const Mesh& otherMesh = *_meshes[otherMeshID];
1279 
1280  if (_options.ButlerVolmer)
1281  {
1282  bool Cathode = false;
1283  bool Anode = false;
1284  if (n == _options["ButlerVolmerCathodeShellMeshID"])
1285  {
1286  Cathode = true;
1287  }
1288  else if (n == _options["ButlerVolmerAnodeShellMeshID"])
1289  {
1290  Anode = true;
1291  }
1292 
1293  BatteryLinearizeSpeciesInterface<T, T, T> lbv (_geomFields,
1294  _options["ButlerVolmerRRConstant"],
1295  _options["interfaceSpeciesUnderRelax"],
1296  Anode,
1297  Cathode,
1298  sFields.concentration,
1299  sFields.concentration,
1300  _batteryModelFields.potential,
1301  _batteryModelFields.temperature);
1302 
1303  lbv.discretize(mesh, parentMesh, otherMesh, ls.getMatrix(), ls.getX(), ls.getB() );
1304 
1305  }
1306  else
1307  {
1308  LinearizeInterfaceJump<T, T, T> lsm (T(1.0),
1309  T(0.0),
1310  sFields.concentration);
1311 
1312  lsm.discretize(mesh, parentMesh, otherMesh, ls.getMatrix(), ls.getX(), ls.getB() );
1313  }
1314  }
1315  /*
1316  if (n==0)
1317  {
1318  cout << "BEFORE" << endl;
1319  printMatrixElementsOnFace(mesh, 8, ls, sFields.concentration);
1320  }
1321  else if (n==1)
1322  printMatrixElementsOnFace(mesh, 11, ls, sFields.concentration);
1323  */
1324  if ((mesh.isDoubleShell())&&(!(mesh.isConnectedShell())))
1325  {
1326  const int parentMeshID = mesh.getParentMeshID();
1327  const int otherMeshID = mesh.getOtherMeshID();
1328  const Mesh& parentMesh = *_meshes[parentMeshID];
1329  const Mesh& otherMesh = *_meshes[otherMeshID];
1330  LinearizeInterfaceJumpUnconnected<T, T, T> lsm (T(1.0),
1331  T(0.0),
1332  sFields.concentration);
1333 
1334  lsm.discretize(mesh, parentMesh, otherMesh, ls.getMatrix(), ls.getX(), ls.getB() , lsShell);
1335 
1336  //sync for debugging
1337  //sync doesn't effect result, but matches what the matrix/residual now contains
1338  sFields.concentration.syncLocal();
1339  }
1340  }
1341  /*
1342  cout << "AFTER" << endl;
1343  printMatrixElementsOnFace(*_meshes[0], 8, ls, sFields.concentration);
1344  printMatrixElementsOnFace(*_meshes[1], 11, ls, sFields.concentration);
1345  */
1346 
1348  const int numRealMeshes = _realMeshes.size();
1349  for (int n=0; n<numRealMeshes; n++)
1350  {
1351  const Mesh& mesh = *_meshes[n];
1352 
1353  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
1354  {
1355  const FaceGroup& fg = *fgPtr;
1356  const StorageSite& faces = fg.site;
1357 
1358  typename BatterySpeciesBCMap::const_iterator pos = sbcmap.find(fg.id);
1359  if (pos==sbcmap.end())
1360  {
1361  throw CException("BatteryModel: Error in Species BC Map");
1362  }
1363 
1364  const BatterySpeciesBC<T>& bc = *(pos->second);
1365 
1366  GenericBCS<T,T,T> gbc(faces,mesh,
1367  _geomFields,
1368  sFields.concentration,
1369  sFields.massFlux,
1370  ls.getMatrix(), ls.getX(), ls.getB());
1371 
1372  if (bc.bcType == "SpecifiedConcentration")
1373  {
1374  FloatValEvaluator<T>
1375  bT(bc.getVal("specifiedConcentration"),faces);
1376  if (sFields.convectionFlux.hasArray(faces))
1377  {
1378  const TArray& convectingFlux =
1379  dynamic_cast<const TArray&>
1380  (sFields.convectionFlux[faces]);
1381  const int nFaces = faces.getCount();
1382 
1383  for(int f=0; f<nFaces; f++)
1384  {
1385  if (convectingFlux[f] > 0.)
1386  {
1387  gbc.applyExtrapolationBC(f);
1388  }
1389  else
1390  {
1391  gbc.applyDirichletBC(f,bT[f]);
1392  }
1393  }
1394  }
1395  else
1396  gbc.applyDirichletBC(bT);
1397  }
1398  else if (bc.bcType == "SpecifiedMassFlux")
1399  {
1400  const T specifiedFlux(bc["specifiedMassFlux"]);
1401  gbc.applyNeumannBC(specifiedFlux);
1402  }
1403  else if ((bc.bcType == "Symmetry"))
1404  {
1405  T zeroFlux(NumTypeTraits<T>::getZero());
1406  gbc.applyNeumannBC(zeroFlux);
1407  }
1408  else
1409  throw CException(bc.bcType + " not implemented for Species in BatteryModel");
1410  }
1411 
1412 
1413  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
1414  {
1415  const FaceGroup& fg = *fgPtr;
1416  const StorageSite& faces = fg.site;
1417  GenericBCS<T,T,T> gbc(faces,mesh,
1418  _geomFields,
1419  sFields.concentration,
1420  sFields.massFlux,
1421  ls.getMatrix(), ls.getX(), ls.getB());
1422 
1423  gbc.applyInterfaceBC();
1424  }
1425 
1426  }
1427  }
1428 
1429 void linearizePotential(LinearSystem& ls)
1430  {
1431  const int numMeshes = _meshes.size();
1432 
1433  GradientModel<T> potentialGradientModel(_meshes,_batteryModelFields.potential,
1434  _batteryModelFields.potential_gradient,_geomFields);
1435  potentialGradientModel.compute();
1436 
1437  DiscrList discretizations;
1438 
1439  shared_ptr<Discretization>
1440  dd(new DiffusionDiscretization<T,T,T>(_meshes,_geomFields,
1441  _batteryModelFields.potential,
1442  _batteryModelFields.conductivity,
1443  _batteryModelFields.potential_gradient));
1444  discretizations.push_back(dd);
1445 
1446  if(_options.ButlerVolmer)
1447  {
1448  bool foundElectrolyte = false;
1449  //populate lnSpeciesConc Field
1450  for (int n=0; n<numMeshes; n++)
1451  {
1452  if (n==_options["BatteryElectrolyteMeshID"])
1453  {
1454  foundElectrolyte = true;
1455  const Mesh& mesh = *_meshes[n];
1456  const StorageSite& cells = mesh.getCells();
1457  BatterySpeciesFields& sFields = *_speciesFieldsVector[0];
1458  const TArray& speciesConc = dynamic_cast<const TArray&>(sFields.concentration[cells]);
1459  TArray& lnSpeciesConc = dynamic_cast<TArray&>(_batteryModelFields.lnLithiumConcentration[cells]);
1460  for (int c=0; c<cells.getCount(); c++)
1461  {
1462  T CellConc = speciesConc[c];
1463  if (CellConc <= 0.0)
1464  {
1465  cout << "Error: Cell Concentration <= 0 MeshID: " << n << " CellNum: " << c << endl;
1466  CellConc = 0.01;
1467  }
1468  lnSpeciesConc[c] = log(CellConc);
1469  }
1470  }
1471  }
1472 
1473  if ((!(foundElectrolyte))&&(_options.ButlerVolmer))
1474  cout << "Warning: Electrolyte Mesh ID not set." << endl;
1475 
1476  //compute gradient for ln term discretization
1477  GradientModel<T> lnSpeciesGradientModel(_meshes,_batteryModelFields.lnLithiumConcentration,
1478  _batteryModelFields.lnLithiumConcentrationGradient,_geomFields);
1479  lnSpeciesGradientModel.compute();
1480 
1481  //discretize ln term (only affects electrolyte mesh)
1482  shared_ptr<Discretization>
1483  bedd(new BatteryBinaryElectrolyteDiscretization<T,T,T>(_meshes,_geomFields,
1484  _batteryModelFields.potential,
1485  _batteryModelFields.conductivity,
1486  _batteryModelFields.lnLithiumConcentration,
1487  _batteryModelFields.lnLithiumConcentrationGradient,
1488  _batteryModelFields.temperature,
1489  _options["BatteryElectrolyteMeshID"]));
1490  discretizations.push_back(bedd);
1491 
1492  }
1493 
1494  Linearizer linearizer;
1495 
1496  linearizer.linearize(discretizations,_meshes,ls.getMatrix(),
1497  ls.getX(), ls.getB());
1498 
1499 
1500 
1501 
1502 
1503 
1504  /* linearize double shell mesh for interface potential jump*/
1505 
1506  for (int n=0; n<numMeshes; n++)
1507  {
1508  const Mesh& mesh = *_meshes[n];
1509  if (mesh.isDoubleShell())
1510  {
1511  const int parentMeshID = mesh.getParentMeshID();
1512  const int otherMeshID = mesh.getOtherMeshID();
1513  const Mesh& parentMesh = *_meshes[parentMeshID];
1514  const Mesh& otherMesh = *_meshes[otherMeshID];
1515 
1516  if (_options.ButlerVolmer)
1517  {
1518  bool Cathode = false;
1519  bool Anode = false;
1520  if (n == _options["ButlerVolmerCathodeShellMeshID"])
1521  {
1522  Cathode = true;
1523  }
1524  else if (n == _options["ButlerVolmerAnodeShellMeshID"])
1525  {
1526  Anode = true;
1527  }
1528 
1529  BatterySpeciesFields& sFields = *_speciesFieldsVector[0];
1530  BatteryLinearizePotentialInterface<T, T, T> lbv (_geomFields,
1531  _batteryModelFields.potential,
1532  sFields.concentration,
1533  _batteryModelFields.temperature,
1534  _options["ButlerVolmerRRConstant"],
1535  Anode,
1536  Cathode);
1537 
1538  lbv.discretize(mesh, parentMesh, otherMesh, ls.getMatrix(), ls.getX(), ls.getB() );
1539  }
1540  else
1541  {
1542 
1543  LinearizeInterfaceJump<T, T, T> lsm (T(1.0),
1544  T(0.0),
1545  _batteryModelFields.potential);
1546 
1547  lsm.discretize(mesh, parentMesh, otherMesh, ls.getMatrix(), ls.getX(), ls.getB() );
1548  }
1549  }
1550  }
1551 
1552  /* boundary and interface condition */
1553 
1554  for (int n=0; n<numMeshes; n++)
1555  {
1556  const Mesh& mesh = *_meshes[n];
1557 
1558  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
1559  {
1560  const FaceGroup& fg = *fgPtr;
1561  const StorageSite& faces = fg.site;
1562  const int nFaces = faces.getCount();
1563  const BatteryPotentialBC<T>& bc = *_pbcMap[fg.id];
1564 
1565 
1566  GenericBCS<T,T,T> gbc(faces,mesh,
1567  _geomFields,
1568  _batteryModelFields.potential,
1569  _batteryModelFields.potential_flux,
1570  ls.getMatrix(), ls.getX(), ls.getB());
1571 
1572  if (bc.bcType == "SpecifiedPotential")
1573  {
1574  FloatValEvaluator<T> bT(bc.getVal("specifiedPotential"), faces);
1575  for(int f=0; f<nFaces; f++)
1576  {
1577  gbc.applyDirichletBC(f, bT[f]);
1578  }
1579  }
1580  else if (bc.bcType == "SpecifiedPotentialFlux")
1581  {
1582  const T specifiedFlux(bc["specifiedPotentialFlux"]);
1583  gbc.applyNeumannBC(specifiedFlux);
1584  }
1585  else if (bc.bcType == "Symmetry")
1586  {
1587  T zeroFlux(NumTypeTraits<T>::getZero());
1588  gbc.applyNeumannBC(zeroFlux);
1589  }
1590  else
1591  throw CException(bc.bcType + " not implemented for Potential in BatteryModel");
1592  }
1593 
1594  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
1595  {
1596 
1597  const FaceGroup& fg = *fgPtr;
1598  const StorageSite& faces = fg.site;
1599 
1600 
1601  GenericBCS<T,T,T> gbc(faces,mesh,
1602  _geomFields,
1603  _batteryModelFields.potential,
1604  _batteryModelFields.potential_flux,
1605  ls.getMatrix(), ls.getX(), ls.getB());
1606 
1607  gbc.applyInterfaceBC();
1608  }
1609  }
1610  }
1611 
1612 void linearizeThermal(LinearSystem& ls)
1613  {
1614 
1615  const int numMeshes = _meshes.size();
1616 
1617  // calculate new thermal gradient
1618  GradientModel<T> thermalGradientModel(_meshes,_batteryModelFields.temperature,
1619  _batteryModelFields.temperatureGradient,_geomFields);
1620  thermalGradientModel.compute();
1621 
1622  // recalculate other gradients for use in heat source term
1623  GradientModel<T> potentialGradientModel(_meshes,_batteryModelFields.potential,
1624  _batteryModelFields.potential_gradient,_geomFields);
1625  potentialGradientModel.compute();
1626 
1627  BatterySpeciesFields& sFields = *_speciesFieldsVector[0];
1628  GradientModel<T> speciesGradientModel(_meshes,sFields.concentration,
1629  _batteryModelFields.speciesGradient,_geomFields);
1630  speciesGradientModel.compute();
1631 
1632 
1633  // fill source field with dot product of gradients from species and potential
1634  for (int n=0; n<numMeshes; n++)
1635  {
1636  const Mesh& mesh = *_meshes[n];
1637  const StorageSite& cells = mesh.getCells();
1638 
1639  const BatterySpeciesVCMap& svcmap = *_svcMapVector[0];
1640  typename BatterySpeciesVCMap::const_iterator pos = svcmap.find(n);
1641  if (pos==svcmap.end())
1642  {
1643  throw CException("BatteryModel: Error in Species VC Map");
1644  }
1645  const BatterySpeciesVC<T>& svc = *(pos->second);
1646  const T massDiffusivity = svc["massDiffusivity"];
1647 
1648  const TGradArray& potentialGradCell = dynamic_cast<const TGradArray&>(_batteryModelFields.potential_gradient[cells]);
1649  const TGradArray& speciesGradCell = dynamic_cast<const TGradArray&>(_batteryModelFields.speciesGradient[cells]);
1650  TArray& heatSource = dynamic_cast<TArray&>(_batteryModelFields.heatSource[cells]);
1651  for (int c=0; c<cells.getCount(); c++)
1652  {
1653  const TGradType pGrad = potentialGradCell[c];
1654  const TGradType sGrad = speciesGradCell[c];
1655  T CellSource = pGrad[0]*sGrad[0] + pGrad[1]*sGrad[1];
1656  if ((*_meshes[0]).getDimension() == 3)
1657  CellSource += pGrad[2]*sGrad[2];
1658  //heatSource[c] = CellSource;
1659  heatSource[c] = CellSource*massDiffusivity*96485.0;
1660  }
1661  }
1662 
1663  DiscrList discretizations;
1664 
1665  shared_ptr<Discretization>
1666  dd(new DiffusionDiscretization<T,T,T>(_meshes,_geomFields,
1667  _batteryModelFields.temperature,
1668  _batteryModelFields.thermalConductivity,
1669  _batteryModelFields.temperatureGradient));
1670  discretizations.push_back(dd);
1671 
1672  if (_options.transient)
1673  {
1674  shared_ptr<Discretization>
1675  td(new TimeDerivativeDiscretization<T,T,T>
1676  (_meshes,_geomFields,
1677  _batteryModelFields.temperature,
1678  _batteryModelFields.temperatureN1,
1679  _batteryModelFields.temperatureN2,
1680  _batteryModelFields.rhoCp,
1681  _options["timeStep"]));
1682 
1683  discretizations.push_back(td);
1684  }
1685 
1686  shared_ptr<Discretization>
1687  sd(new SourceDiscretization<T>
1688  (_meshes,
1689  _geomFields,
1690  _batteryModelFields.temperature,
1691  _batteryModelFields.heatSource));
1692  discretizations.push_back(sd);
1693 
1694  Linearizer linearizer;
1695 
1696  linearizer.linearize(discretizations,_meshes,ls.getMatrix(),
1697  ls.getX(), ls.getB());
1698 
1699 
1700 
1701 
1702 
1703  /* linearize double shell mesh for thermal special interface*/
1704 
1705  for (int n=0; n<numMeshes; n++)
1706  {
1707  const Mesh& mesh = *_meshes[n];
1708  if (mesh.isDoubleShell())
1709  {
1710  const int parentMeshID = mesh.getParentMeshID();
1711  const int otherMeshID = mesh.getOtherMeshID();
1712  const Mesh& parentMesh = *_meshes[parentMeshID];
1713  const Mesh& otherMesh = *_meshes[otherMeshID];
1714 
1715  if (_options.ButlerVolmer)
1716  {
1717 
1718  bool Cathode = false;
1719  bool Anode = false;
1720  if (n == _options["ButlerVolmerCathodeShellMeshID"])
1721  {
1722  Cathode = true;
1723  }
1724  else if (n == _options["ButlerVolmerAnodeShellMeshID"])
1725  {
1726  Anode = true;
1727  }
1728 
1729  BatterySpeciesFields& sFields = *_speciesFieldsVector[0];
1730  BatteryLinearizeThermalInterface<T, T, T> lbv (_geomFields,
1731  _batteryModelFields.temperature,
1732  sFields.concentration,
1733  _batteryModelFields.potential,
1734  _options["ButlerVolmerRRConstant"],
1735  Anode,
1736  Cathode);
1737 
1738  lbv.discretize(mesh, parentMesh, otherMesh, ls.getMatrix(), ls.getX(), ls.getB() );
1739 
1740  }
1741  else
1742  {
1743 
1744  LinearizeInterfaceJump<T, T, T> lsm (T(1.0),
1745  T(0.0),
1746  _batteryModelFields.temperature);
1747 
1748  lsm.discretize(mesh, parentMesh, otherMesh, ls.getMatrix(), ls.getX(), ls.getB() );
1749  }
1750  }
1751  }
1752 
1753  /* boundary and interface condition */
1754 
1755  for (int n=0; n<numMeshes; n++)
1756  {
1757  const Mesh& mesh = *_meshes[n];
1758 
1759  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
1760  {
1761  const FaceGroup& fg = *fgPtr;
1762  const StorageSite& faces = fg.site;
1763  const int nFaces = faces.getCount();
1764  const BatteryThermalBC<T>& bc = *_tbcMap[fg.id];
1765 
1766 
1767  GenericBCS<T,T,T> gbc(faces,mesh,
1768  _geomFields,
1769  _batteryModelFields.temperature,
1770  _batteryModelFields.heatFlux,
1771  ls.getMatrix(), ls.getX(), ls.getB());
1772 
1773  if (bc.bcType == "SpecifiedTemperature")
1774  {
1775  FloatValEvaluator<T> bT(bc.getVal("specifiedTemperature"), faces);
1776  for(int f=0; f<nFaces; f++)
1777  {
1778  gbc.applyDirichletBC(f, bT[f]);
1779  }
1780  }
1781  else if (bc.bcType == "SpecifiedHeatFlux")
1782  {
1783  const T specifiedFlux(bc["specifiedHeatFlux"]);
1784  gbc.applyNeumannBC(specifiedFlux);
1785  }
1786  else if (bc.bcType == "Symmetry")
1787  {
1788  T zeroFlux(NumTypeTraits<T>::getZero());
1789  gbc.applyNeumannBC(zeroFlux);
1790  }
1791  else
1792  throw CException(bc.bcType + " not implemented for Thermal in BatteryModel");
1793  }
1794 
1795  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
1796  {
1797 
1798  const FaceGroup& fg = *fgPtr;
1799  const StorageSite& faces = fg.site;
1800 
1801 
1802  GenericBCS<T,T,T> gbc(faces,mesh,
1803  _geomFields,
1804  _batteryModelFields.temperature,
1805  _batteryModelFields.heatFlux,
1806  ls.getMatrix(), ls.getX(), ls.getB());
1807 
1808  gbc.applyInterfaceBC();
1809  }
1810  }
1811  }
1812 
1813 
1814 void linearizePC_Thermal(LinearSystem& ls)
1815  {
1816 
1817  // only lithium (species 0) for now
1818  const BatterySpeciesVCMap& svcmap = *_svcMapVector[0];
1819  const BatterySpeciesBCMap& sbcmap = *_sbcMapVector[0];
1820 
1821  // get combined gradients of potential and species
1822  GradientModel<VectorT3> pstGradientModel(_meshes,_batteryModelFields.potentialSpeciesTemp,
1823  _batteryModelFields.potentialSpeciesTempGradient,_geomFields);
1824  pstGradientModel.compute();
1825 
1826  //populate lnSpeciesConc fields for inclusion of ln term in potential equation
1827  bool foundElectrolyte = false;
1828 
1829  const int numMeshes = _meshes.size();
1830  for (int n=0; n<numMeshes; n++)
1831  {
1832  if (n==_options["BatteryElectrolyteMeshID"])
1833  {
1834  foundElectrolyte = true;
1835  const Mesh& mesh = *_meshes[n];
1836  const StorageSite& cells = mesh.getCells();
1837  const VectorT3Array& speciesPotentialTemp = dynamic_cast<const VectorT3Array&>(_batteryModelFields.potentialSpeciesTemp[cells]);
1838  TArray& lnSpeciesConc = dynamic_cast<TArray&>(_batteryModelFields.lnLithiumConcentration[cells]);
1839  for (int c=0; c<cells.getCount(); c++)
1840  {
1841  T CellConc = (speciesPotentialTemp[c])[1];
1842  if (CellConc <= 0.0)
1843  {
1844  cout << "Error: Cell Concentration <= 0 MeshID: " << n << " CellNum: " << c << endl;
1845  CellConc = 0.01;
1846  }
1847  lnSpeciesConc[c] = log(CellConc);
1848  }
1849  }
1850  }
1851 
1852  if ((!(foundElectrolyte))&&(_options.ButlerVolmer))
1853  cout << "Warning: Electrolyte Mesh ID not set." << endl;
1854 
1855  //compute gradient for ln term discretization
1856  GradientModel<T> lnSpeciesGradientModel(_meshes,_batteryModelFields.lnLithiumConcentration,
1857  _batteryModelFields.lnLithiumConcentrationGradient,_geomFields);
1858  lnSpeciesGradientModel.compute();
1859 
1860  // fill source field with dot product of gradients from species and potential
1861  for (int n=0; n<numMeshes; n++)
1862  {
1863  const Mesh& mesh = *_meshes[n];
1864  const StorageSite& cells = mesh.getCells();
1865 
1866  typename BatterySpeciesVCMap::const_iterator pos = svcmap.find(n);
1867  if (pos==svcmap.end())
1868  {
1869  throw CException("BatteryModel: Error in Species VC Map");
1870  }
1871  const BatterySpeciesVC<T>& svc = *(pos->second);
1872  const T massDiffusivity = svc["massDiffusivity"];
1873 
1874  const VectorT3GradArray& pstGradCell = dynamic_cast<const VectorT3GradArray&>(_batteryModelFields.potentialSpeciesTempGradient[cells]);
1875  TArray& heatSource = dynamic_cast<TArray&>(_batteryModelFields.heatSource[cells]);
1876  for (int c=0; c<cells.getCount(); c++)
1877  {
1878  const VectorT3Grad combinedCellGradient = pstGradCell[c];
1879  TGradType pGrad(NumTypeTraits<TGradType>::getZero());
1880  TGradType sGrad(NumTypeTraits<TGradType>::getZero());
1881  pGrad[0] = combinedCellGradient[0][0];
1882  pGrad[1] = combinedCellGradient[1][0];
1883  sGrad[0] = combinedCellGradient[0][1];
1884  sGrad[1] = combinedCellGradient[1][1];
1885  T CellSource = pGrad[0]*sGrad[0] + pGrad[1]*sGrad[1];
1886  if ((*_meshes[0]).getDimension() == 3)
1887  {
1888  pGrad[2] = combinedCellGradient[2][0];
1889  sGrad[2] = combinedCellGradient[2][1];
1890  CellSource += pGrad[2]*sGrad[2];
1891  }
1892  heatSource[c] = CellSource*massDiffusivity*96485.0;
1893  }
1894  }
1895 
1896 
1897  DiscrList discretizations;
1898 
1899  shared_ptr<Discretization>
1900  dd(new BatteryPCDiffusionDiscretization<VectorT3,DiagTensorT3,DiagTensorT3>
1901  (_realMeshes,_geomFields,
1902  _batteryModelFields.potentialSpeciesTemp,
1903  _batteryModelFields.potentialSpeciesTempDiffusivity,
1904  _batteryModelFields.potentialSpeciesTempGradient));
1905  discretizations.push_back(dd);
1906 
1907  //discretize ln term (only affects electrolyte mesh of potential model)
1908  shared_ptr<Discretization>
1909  bedd(new BatteryPCBinaryElectrolyteDiscretization<VectorT3,DiagTensorT3,DiagTensorT3>(_realMeshes,_geomFields,
1910  _batteryModelFields.potentialSpeciesTemp,
1911  _batteryModelFields.potentialSpeciesTempDiffusivity,
1912  _batteryModelFields.lnLithiumConcentration,
1913  _batteryModelFields.lnLithiumConcentrationGradient,
1914  _options.thermalModelPC,
1915  _options["BatteryElectrolyteMeshID"]));
1916  discretizations.push_back(bedd);
1917 
1918  if (_options.transient)
1919  {
1920  shared_ptr<Discretization>
1921  td(new BatteryPCTimeDerivativeDiscretization<VectorT3,DiagTensorT3,DiagTensorT3>
1922  (_realMeshes,_geomFields,
1923  _batteryModelFields.potentialSpeciesTemp,
1924  _batteryModelFields.potentialSpeciesTempN1,
1925  _batteryModelFields.potentialSpeciesTempN2,
1926  _batteryModelFields.rhoCp,
1927  _options.thermalModelPC,
1928  _options["timeStep"]));
1929 
1930  discretizations.push_back(td);
1931  }
1932 
1933  if (1)
1934  {
1935  shared_ptr<Discretization>
1936  sd(new BatteryPCHeatSourceDiscretization<VectorT3>
1937  (_realMeshes,_geomFields,
1938  _batteryModelFields.potentialSpeciesTemp,
1939  _batteryModelFields.heatSource));
1940 
1941  discretizations.push_back(sd);
1942  }
1943 
1944  Linearizer linearizer;
1945 
1946  linearizer.linearize(discretizations,_realMeshes,ls.getMatrix(),
1947  ls.getX(), ls.getB());
1948 
1949  // linearize shell mesh
1950 
1951  for (int n=0; n<numMeshes; n++)
1952  {
1953  const Mesh& mesh = *_meshes[n];
1954  if (mesh.isDoubleShell())
1955  {
1956  const int parentMeshID = mesh.getParentMeshID();
1957  const int otherMeshID = mesh.getOtherMeshID();
1958  const Mesh& parentMesh = *_meshes[parentMeshID];
1959  const Mesh& otherMesh = *_meshes[otherMeshID];
1960 
1961  if (_options.ButlerVolmer)
1962  {
1963  bool Cathode = false;
1964  bool Anode = false;
1965  if (n == _options["ButlerVolmerCathodeShellMeshID"])
1966  {
1967  Cathode = true;
1968  }
1969  else if (n == _options["ButlerVolmerAnodeShellMeshID"])
1970  {
1971  Anode = true;
1972  }
1973 
1974  BatteryPCLinearizeInterface_BV<VectorT3, SquareTensorT3, SquareTensorT3, DiagTensorT3> lbv (_geomFields,
1975  _options["ButlerVolmerRRConstant"],
1976  _options["interfaceSpeciesUnderRelax"],
1977  Anode,
1978  Cathode,
1979  _options.thermalModelPC,
1980  _batteryModelFields.potentialSpeciesTemp);
1981 
1982  lbv.discretize(mesh, parentMesh, otherMesh, ls.getMatrix(), ls.getX(), ls.getB() );
1983 
1984  }
1985  else
1986  {
1987  LinearizeInterfaceJump<VectorT3, SquareTensorT3, SquareTensorT3> lsm (T(1.0),
1988  NumTypeTraits<VectorT3>::getZero(),
1989  _batteryModelFields.potentialSpeciesTemp);
1990 
1991  lsm.discretize(mesh, parentMesh, otherMesh, ls.getMatrix(), ls.getX(), ls.getB() );
1992  }
1993  }
1994  }
1995 
1997 
1998  for (int n=0; n<numMeshes; n++)
1999  {
2000  const Mesh& mesh = *_meshes[n];
2001 
2002  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
2003  {
2004  const FaceGroup& fg = *fgPtr;
2005  const StorageSite& faces = fg.site;
2006 
2007  typename BatterySpeciesBCMap::const_iterator pos = sbcmap.find(fg.id);
2008  if (pos==sbcmap.end())
2009  {
2010  throw CException("BatteryModel: Error in Species BC Map");
2011  }
2012 
2013  const BatterySpeciesBC<T>& sbc = *(pos->second);
2014 
2015  const BatteryPotentialBC<T>& pbc = *_pbcMap[fg.id];
2016 
2017  const BatteryThermalBC<T>& tbc = *_tbcMap[fg.id];
2018 
2019  BatteryPC_BCS<VectorT3,DiagTensorT3,DiagTensorT3> bbc(faces,mesh,
2020  _geomFields,
2021  _batteryModelFields.potentialSpeciesTemp,
2022  _batteryModelFields.potentialSpeciesTempFlux,
2023  ls.getMatrix(), ls.getX(), ls.getB());
2024 
2025  if (sbc.bcType == "SpecifiedConcentration")
2026  {
2027  FloatValEvaluator<T>
2028  sbT(sbc.getVal("specifiedConcentration"),faces);
2029 
2030  bbc.applySingleEquationDirichletBC(sbT,1);
2031  }
2032  else if (sbc.bcType == "SpecifiedMassFlux")
2033  {
2034  const T specifiedFlux(sbc["specifiedMassFlux"]);
2035  bbc.applySingleEquationNeumannBC(specifiedFlux,1);
2036  }
2037  else if ((sbc.bcType == "Symmetry"))
2038  {
2039  T zeroFlux(NumTypeTraits<T>::getZero());
2040  bbc.applySingleEquationNeumannBC(zeroFlux,1);
2041  }
2042  else
2043  throw CException(sbc.bcType + " not implemented for Species in BatteryModel");
2044 
2045  if (pbc.bcType == "SpecifiedPotential")
2046  {
2047  FloatValEvaluator<T> pbT(pbc.getVal("specifiedPotential"),faces);
2048 
2049  bbc.applySingleEquationDirichletBC(pbT,0);
2050  }
2051  else if (pbc.bcType == "SpecifiedPotentialFlux")
2052  {
2053  const T specifiedFlux(pbc["specifiedPotentialFlux"]);
2054  bbc.applySingleEquationNeumannBC(specifiedFlux,0);
2055  }
2056  else if ((pbc.bcType == "Symmetry"))
2057  {
2058  T zeroFlux(NumTypeTraits<T>::getZero());
2059  bbc.applySingleEquationNeumannBC(zeroFlux,0);
2060  }
2061  else
2062  throw CException(pbc.bcType + " not implemented for potential in BatteryModel");
2063 
2064  if (tbc.bcType == "SpecifiedTemperature")
2065  {
2066  FloatValEvaluator<T> tbT(tbc.getVal("specifiedTemperature"),faces);
2067 
2068  bbc.applySingleEquationDirichletBC(tbT,2);
2069  }
2070  else if (tbc.bcType == "SpecifiedHeatFlux")
2071  {
2072  const T specifiedFlux(tbc["specifiedHeatFlux"]);
2073  bbc.applySingleEquationNeumannBC(specifiedFlux,2);
2074  }
2075  else if ((tbc.bcType == "Symmetry"))
2076  {
2077  T zeroFlux(NumTypeTraits<T>::getZero());
2078  bbc.applySingleEquationNeumannBC(zeroFlux,2);
2079  }
2080  else
2081  throw CException(tbc.bcType + " not implemented for temperature in BatteryModel");
2082 
2083  }
2084 
2085  if (mesh.isDoubleShell())
2086  {
2087  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
2088  {
2089  const FaceGroup& fg = *fgPtr;
2090  const StorageSite& faces = fg.site;
2091  GenericBCS<VectorT3,SquareTensorT3,SquareTensorT3> gbc(faces,mesh,
2092  _geomFields,
2093  _batteryModelFields.potentialSpeciesTemp,
2094  _batteryModelFields.potentialSpeciesTempFlux,
2095  ls.getMatrix(), ls.getX(), ls.getB());
2096 
2097  gbc.applyInterfaceBC();
2098  }
2099  }
2100  else
2101  {
2102  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
2103  {
2104  const FaceGroup& fg = *fgPtr;
2105  const StorageSite& faces = fg.site;
2106  GenericBCS<VectorT3,DiagTensorT3,DiagTensorT3> gbc(faces,mesh,
2107  _geomFields,
2108  _batteryModelFields.potentialSpeciesTemp,
2109  _batteryModelFields.potentialSpeciesTempFlux,
2110  ls.getMatrix(), ls.getX(), ls.getB());
2111 
2112  gbc.applyInterfaceBC();
2113  }
2114  }
2115  //set all temperature residuals to zero if no thermal model PC (comes in to play if trying to still do advanceThermal separatly)
2116 /*
2117  if (!_options.thermalModelPC)
2118  {
2119  for (int n=0; n<numMeshes; n++)
2120  {
2121  const Mesh& mesh = *_meshes[n];
2122  const StorageSite& cells = mesh.getCells();
2123  const MultiField::ArrayIndex cVarIndex(&_batteryModelFields.potentialSpeciesTemp,&cells);
2124  VectorT3Array& rCell = dynamic_cast<VectorT3Array&>((ls.getB())[cVarIndex]);
2125  for (int c=0; c<cells.getCount(); c++)
2126  {
2127  (rCell[c])[2] = 0.0;
2128  }
2129  }
2130  }*/
2131  }
2132  }
2133 
2134 void linearizePC(LinearSystem& ls)
2135  {
2136  // only lithium (species 0) for now
2137  const BatterySpeciesBCMap& sbcmap = *_sbcMapVector[0];
2138 
2139  // get combined gradients of potential and species
2140  GradientModel<VectorT2> pstGradientModel(_meshes,_batteryModelFields.potentialSpeciesTemp,
2141  _batteryModelFields.potentialSpeciesTempGradient,_geomFields);
2142  pstGradientModel.compute();
2143 
2144  //populate lnSpeciesConc fields for inclusion of ln term in potential equation
2145  bool foundElectrolyte = false;
2146 
2147  const int numMeshes = _meshes.size();
2148  for (int n=0; n<numMeshes; n++)
2149  {
2150  if (n==_options["BatteryElectrolyteMeshID"])
2151  {
2152  foundElectrolyte = true;
2153  const Mesh& mesh = *_meshes[n];
2154  const StorageSite& cells = mesh.getCells();
2155  const VectorT2Array& speciesPotentialTemp = dynamic_cast<const VectorT2Array&>(_batteryModelFields.potentialSpeciesTemp[cells]);
2156  TArray& lnSpeciesConc = dynamic_cast<TArray&>(_batteryModelFields.lnLithiumConcentration[cells]);
2157  for (int c=0; c<cells.getCount(); c++)
2158  {
2159  T CellConc = (speciesPotentialTemp[c])[1];
2160  if (CellConc <= 0.0)
2161  {
2162  cout << "Error: Cell Concentration <= 0 MeshID: " << n << " CellNum: " << c << endl;
2163  CellConc = 0.01;
2164  }
2165  lnSpeciesConc[c] = log(CellConc);
2166  }
2167  }
2168  }
2169 
2170  if ((!(foundElectrolyte))&&(_options.ButlerVolmer))
2171  cout << "Warning: Electrolyte Mesh ID not set." << endl;
2172 
2173  //compute gradient for ln term discretization
2174  GradientModel<T> lnSpeciesGradientModel(_meshes,_batteryModelFields.lnLithiumConcentration,
2175  _batteryModelFields.lnLithiumConcentrationGradient,_geomFields);
2176  lnSpeciesGradientModel.compute();
2177 
2178 
2179  DiscrList discretizations;
2180 
2181  shared_ptr<Discretization>
2182  dd(new BatteryPCDiffusionDiscretization<VectorT2,DiagTensorT2,DiagTensorT2>
2183  (_realMeshes,_geomFields,
2184  _batteryModelFields.potentialSpeciesTemp,
2185  _batteryModelFields.potentialSpeciesTempDiffusivity,
2186  _batteryModelFields.potentialSpeciesTempGradient));
2187  discretizations.push_back(dd);
2188 
2189  //discretize ln term (only affects electrolyte mesh of potential model)
2190  shared_ptr<Discretization>
2191  bedd(new BatteryPCBinaryElectrolyteDiscretization<VectorT2,DiagTensorT2,DiagTensorT2>(_realMeshes,_geomFields,
2192  _batteryModelFields.potentialSpeciesTemp,
2193  _batteryModelFields.potentialSpeciesTempDiffusivity,
2194  _batteryModelFields.lnLithiumConcentration,
2195  _batteryModelFields.lnLithiumConcentrationGradient,
2196  _options.thermalModelPC,
2197  _options["BatteryElectrolyteMeshID"]));
2198  discretizations.push_back(bedd);
2199 
2200  if (_options.transient)
2201  {
2202  shared_ptr<Discretization>
2203  td(new BatteryPCTimeDerivativeDiscretization<VectorT2,DiagTensorT2,DiagTensorT2>
2204  (_realMeshes,_geomFields,
2205  _batteryModelFields.potentialSpeciesTemp,
2206  _batteryModelFields.potentialSpeciesTempN1,
2207  _batteryModelFields.potentialSpeciesTempN2,
2208  _batteryModelFields.rhoCp,
2209  _options.thermalModelPC,
2210  _options["timeStep"]));
2211 
2212  discretizations.push_back(td);
2213  }
2214 
2215  Linearizer linearizer;
2216 
2217  linearizer.linearize(discretizations,_realMeshes,ls.getMatrix(),
2218  ls.getX(), ls.getB());
2219 
2220  // linearize shell mesh
2221 
2222  for (int n=0; n<numMeshes; n++)
2223  {
2224  const Mesh& mesh = *_meshes[n];
2225  if (mesh.isDoubleShell())
2226  {
2227  const int parentMeshID = mesh.getParentMeshID();
2228  const int otherMeshID = mesh.getOtherMeshID();
2229  const Mesh& parentMesh = *_meshes[parentMeshID];
2230  const Mesh& otherMesh = *_meshes[otherMeshID];
2231 
2232  if (_options.ButlerVolmer)
2233  {
2234  bool Cathode = false;
2235  bool Anode = false;
2236  if (n == _options["ButlerVolmerCathodeShellMeshID"])
2237  {
2238  Cathode = true;
2239  }
2240  else if (n == _options["ButlerVolmerAnodeShellMeshID"])
2241  {
2242  Anode = true;
2243  }
2244 
2245  BatteryPCLinearizeInterface_BV<VectorT2, SquareTensorT2, SquareTensorT2, DiagTensorT2> lbv (_geomFields,
2246  _options["ButlerVolmerRRConstant"],
2247  _options["interfaceSpeciesUnderRelax"],
2248  Anode,
2249  Cathode,
2250  _options.thermalModelPC,
2251  _batteryModelFields.potentialSpeciesTemp);
2252 
2253  lbv.discretize(mesh, parentMesh, otherMesh, ls.getMatrix(), ls.getX(), ls.getB() );
2254 
2255  }
2256  else
2257  {
2258  LinearizeInterfaceJump<VectorT2, SquareTensorT2, SquareTensorT2> lsm (T(1.0),
2259  NumTypeTraits<VectorT2>::getZero(),
2260  _batteryModelFields.potentialSpeciesTemp);
2261 
2262  lsm.discretize(mesh, parentMesh, otherMesh, ls.getMatrix(), ls.getX(), ls.getB() );
2263  }
2264  }
2265  }
2266 
2268 
2269  for (int n=0; n<numMeshes; n++)
2270  {
2271  const Mesh& mesh = *_meshes[n];
2272 
2273  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
2274  {
2275  const FaceGroup& fg = *fgPtr;
2276  const StorageSite& faces = fg.site;
2277 
2278  typename BatterySpeciesBCMap::const_iterator pos = sbcmap.find(fg.id);
2279  if (pos==sbcmap.end())
2280  {
2281  throw CException("BatteryModel: Error in Species BC Map");
2282  }
2283 
2284  const BatterySpeciesBC<T>& sbc = *(pos->second);
2285 
2286  const BatteryPotentialBC<T>& pbc = *_pbcMap[fg.id];
2287 
2288  BatteryPC_BCS<VectorT2,DiagTensorT2,DiagTensorT2> bbc(faces,mesh,
2289  _geomFields,
2290  _batteryModelFields.potentialSpeciesTemp,
2291  _batteryModelFields.potentialSpeciesTempFlux,
2292  ls.getMatrix(), ls.getX(), ls.getB());
2293 
2294  if (sbc.bcType == "SpecifiedConcentration")
2295  {
2296  FloatValEvaluator<T>
2297  sbT(sbc.getVal("specifiedConcentration"),faces);
2298 
2299  bbc.applySingleEquationDirichletBC(sbT,1);
2300  }
2301  else if (sbc.bcType == "SpecifiedMassFlux")
2302  {
2303  const T specifiedFlux(sbc["specifiedMassFlux"]);
2304  bbc.applySingleEquationNeumannBC(specifiedFlux,1);
2305  }
2306  else if ((sbc.bcType == "Symmetry"))
2307  {
2308  T zeroFlux(NumTypeTraits<T>::getZero());
2309  bbc.applySingleEquationNeumannBC(zeroFlux,1);
2310  }
2311  else
2312  throw CException(sbc.bcType + " not implemented for Species in BatteryModel");
2313 
2314  if (pbc.bcType == "SpecifiedPotential")
2315  {
2316  FloatValEvaluator<T> pbT(pbc.getVal("specifiedPotential"),faces);
2317 
2318  bbc.applySingleEquationDirichletBC(pbT,0);
2319  }
2320  else if (pbc.bcType == "SpecifiedPotentialFlux")
2321  {
2322  const T specifiedFlux(pbc["specifiedPotentialFlux"]);
2323  bbc.applySingleEquationNeumannBC(specifiedFlux,0);
2324  }
2325  else if ((pbc.bcType == "Symmetry"))
2326  {
2327  T zeroFlux(NumTypeTraits<T>::getZero());
2328  bbc.applySingleEquationNeumannBC(zeroFlux,0);
2329  }
2330  else
2331  throw CException(pbc.bcType + " not implemented for potential in BatteryModel");
2332 
2333  }
2334 
2335  if (mesh.isDoubleShell())
2336  {
2337  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
2338  {
2339  const FaceGroup& fg = *fgPtr;
2340  const StorageSite& faces = fg.site;
2341  GenericBCS<VectorT2,SquareTensorT2,SquareTensorT2> gbc(faces,mesh,
2342  _geomFields,
2343  _batteryModelFields.potentialSpeciesTemp,
2344  _batteryModelFields.potentialSpeciesTempFlux,
2345  ls.getMatrix(), ls.getX(), ls.getB());
2346 
2347  gbc.applyInterfaceBC();
2348  }
2349  }
2350  else
2351  {
2352  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
2353  {
2354  const FaceGroup& fg = *fgPtr;
2355  const StorageSite& faces = fg.site;
2356  GenericBCS<VectorT2,DiagTensorT2,DiagTensorT2> gbc(faces,mesh,
2357  _geomFields,
2358  _batteryModelFields.potentialSpeciesTemp,
2359  _batteryModelFields.potentialSpeciesTempFlux,
2360  ls.getMatrix(), ls.getX(), ls.getB());
2361 
2362  gbc.applyInterfaceBC();
2363  }
2364  }
2365  }
2366  }
2367 
2368 
2369  T getMassFluxIntegral(const Mesh& mesh, const int faceGroupId, const int m)
2370  {
2371  BatterySpeciesFields& sFields = *_speciesFieldsVector[m];
2372 
2373  T r(0.);
2374  bool found = false;
2375  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
2376  {
2377  const FaceGroup& fg = *fgPtr;
2378  if (fg.id == faceGroupId)
2379  {
2380  const StorageSite& faces = fg.site;
2381  const int nFaces = faces.getCount();
2382  const TArray& massFlux =
2383  dynamic_cast<const TArray&>(sFields.massFlux[faces]);
2384  for(int f=0; f<nFaces; f++)
2385  r += massFlux[f];
2386  found=true;
2387  }
2388  }
2389 
2390  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
2391  {
2392  const FaceGroup& fg = *fgPtr;
2393  if (fg.id == faceGroupId)
2394  {
2395  const StorageSite& faces = fg.site;
2396  const int nFaces = faces.getCount();
2397  const TArray& massFlux =
2398  dynamic_cast<const TArray&>(sFields.massFlux[faces]);
2399  for(int f=0; f<nFaces; f++)
2400  r += massFlux[f];
2401  found=true;
2402  }
2403  }
2404 
2405 
2406  if (!found)
2407  throw CException("getMassFluxIntegral: invalid faceGroupID");
2408  return r;
2409  }
2410 
2411 T getPotentialFluxIntegral(const Mesh& mesh, const int faceGroupId)
2412  {
2413 
2414  T r(0.);
2415  bool found = false;
2416  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
2417  {
2418  const FaceGroup& fg = *fgPtr;
2419  if (fg.id == faceGroupId)
2420  {
2421  const StorageSite& faces = fg.site;
2422  const int nFaces = faces.getCount();
2423  const TArray& potentialFlux =
2424  dynamic_cast<const TArray&>(_batteryModelFields.potential_flux[faces]);
2425  for(int f=0; f<nFaces; f++)
2426  r += potentialFlux[f];
2427  found=true;
2428  }
2429  }
2430 
2431  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
2432  {
2433  const FaceGroup& fg = *fgPtr;
2434  if (fg.id == faceGroupId)
2435  {
2436  const StorageSite& faces = fg.site;
2437  const int nFaces = faces.getCount();
2438  const TArray& potentialFlux =
2439  dynamic_cast<const TArray&>(_batteryModelFields.potential_flux[faces]);
2440  for(int f=0; f<nFaces; f++)
2441  r += potentialFlux[f];
2442  found=true;
2443  }
2444  }
2445 
2446 
2447  if (!found)
2448  throw CException("getPotentialFluxIntegral: invalid faceGroupID");
2449  return r;
2450  }
2451 
2452 T getHeatFluxIntegral(const Mesh& mesh, const int faceGroupId)
2453  {
2454 
2455  T r(0.);
2456  bool found = false;
2457  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
2458  {
2459  const FaceGroup& fg = *fgPtr;
2460  if (fg.id == faceGroupId)
2461  {
2462  const StorageSite& faces = fg.site;
2463  const int nFaces = faces.getCount();
2464  const TArray& heatFlux =
2465  dynamic_cast<const TArray&>(_batteryModelFields.heatFlux[faces]);
2466  for(int f=0; f<nFaces; f++)
2467  r += heatFlux[f];
2468  found=true;
2469  }
2470  }
2471 
2472  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
2473  {
2474  const FaceGroup& fg = *fgPtr;
2475  if (fg.id == faceGroupId)
2476  {
2477  const StorageSite& faces = fg.site;
2478  const int nFaces = faces.getCount();
2479  const TArray& heatFlux =
2480  dynamic_cast<const TArray&>(_batteryModelFields.heatFlux[faces]);
2481  for(int f=0; f<nFaces; f++)
2482  r += heatFlux[f];
2483  found=true;
2484  }
2485  }
2486 
2487 
2488  if (!found)
2489  throw CException("getPotentialFluxIntegral: invalid faceGroupID");
2490  return r;
2491  }
2492 
2493 T getAverageConcentration(const Mesh& mesh, const int m)
2494  {
2495  BatterySpeciesFields& sFields = *_speciesFieldsVector[m];
2496  const StorageSite& cells = mesh.getCells();
2497  const TArray& concCell = dynamic_cast<const TArray&>(sFields.concentration[cells]);
2498  const TArray& cellVolume = dynamic_cast<const TArray&>(_geomFields.volume[cells]);
2499  T totalVolume = 0.0;
2500  T weightedConcentration = 0.0;
2501  for(int c=0; c<cells.getSelfCount(); c++)
2502  {
2503  totalVolume += cellVolume[c];
2504  weightedConcentration += concCell[c]*cellVolume[c];
2505  }
2506  return weightedConcentration/totalVolume;
2507  }
2508 
2509 
2510  void printMatrixElementsOnFace(const Mesh& mesh, const int fgId, LinearSystem& ls, Field& varField)
2511  {
2512 
2513  const FaceGroup& fg = mesh.getFaceGroup(fgId);
2514  const StorageSite& faces = fg.site;
2515  const StorageSite& cells = mesh.getCells();
2516  const int nFaces = faces.getCount();
2517  const CRConnectivity& faceCells = mesh.getFaceCells(faces);
2518 
2519  MultiFieldMatrix& mfmatrix = ls.getMatrix();
2520  MultiField& rField = ls.getB();
2521  MultiField& xField = ls.getX();
2522 
2523  const MultiField::ArrayIndex cVarIndex(&varField,&cells);
2524  CRMatrix<T,T,T>& matrix = dynamic_cast<CRMatrix<T,T,T>&>(mfmatrix.getMatrix(cVarIndex,cVarIndex));
2525  TArray& rCell = dynamic_cast<TArray&>(rField[cVarIndex]);
2526  TArray& xCell = dynamic_cast<TArray&>(xField[cVarIndex]);
2527  typename CRMatrix<T,T,T>::DiagArray& diag = matrix.getDiag();
2528 
2529  for(int f=0; f<nFaces; f++)
2530  {
2531  const int c0 = faceCells(f,0);
2532  const int c1 = faceCells(f,1);
2533  T& offdiagC0_C1 = matrix.getCoeff(c0, c1);
2534  T& offdiagC1_C0 = matrix.getCoeff(c1, c0);
2535  const T rC0= rCell[c0];
2536  const T rC1= rCell[c1];
2537  const T diagC0= diag[c0];
2538  const T diagC1= diag[c1];
2539  const T xC0 = xCell[c0];
2540  const T xC1 = xCell[c1];
2541 
2542  if (f==0)
2543  {
2544  cout << "--------------" << endl;
2545  cout << mesh.getID() << ":" << f << endl;
2546  cout << offdiagC0_C1 << endl;
2547  cout << offdiagC1_C0 << endl;
2548  cout << rC0 << endl;
2549  cout << rC1 << endl;
2550  cout << diagC0 << endl;
2551  cout << diagC1 << endl;
2552  cout << xC0 << endl;
2553  cout << xC1 << endl;
2554 
2555  cout << "--------------" << endl;
2556 
2557  }
2558 
2559  }
2560  return;
2561  }
2562 
2563  LinearSystem& matrixMerge(LinearSystem& ls, const int m)
2564  {
2565  /*
2566  typedef CRMatrix<double,double,double> CombinedMatrix;
2567  const MultiFieldMatrix& mfMatrix = ls.getMatrix();
2568  const MultiField& bField = ls.getB();
2569  const MultiField& deltaField = ls.getDelta();
2570  if (bField.getLength() == 1)
2571  return ls;
2572 
2573  const int numMeshes = _meshes.size();
2574  int totalBoundaryCells = 0;
2575  int totalInteriorCells = 0;
2576  for (int n1=0; n1<numMeshes; n1++)
2577  {
2578  Mesh& mesh = *_meshes[n1];
2579  StorageSite& cells = mesh.getCells();
2580  //const CRConnectivity& cellCells = mesh.getCellCells();
2581 
2582 
2583  int n1InterfaceCells = 0;
2584  for (int n2=0; n2<numMeshes; n2++)
2585  {
2586  StorageSite& otherCells = (*_meshes[n2]).getCells();
2587  StorageSite::ScatterMap& n1ScatterMap = cells.getScatterMap();
2588 
2589  typename StorageSite::ScatterMap::const_iterator pos = n1ScatterMap.find(&otherCells);
2590  if (pos==n1ScatterMap.end())
2591  {
2592  //no interfaces with mesh n2
2593  }
2594  else
2595  {
2596  shared_ptr< Array<int> > n1n2ScatterOrigPtr = n1ScatterMap[&otherCells];
2597  n1InterfaceCells += n1n2ScatterOrigPtr->getLength();
2598  }
2599 
2600  //shared_ptr<CRConnectivity> flatConnPtr = origConn.getMultiTranspose(blockSize);
2601 
2602  //CombinedMatrix combinedMatrix(*combinedConnPtr);
2603 
2604  // origMatrix.setFlatMatrix(flatMatrix);
2605  //ls.getMatrix().addMatrix(tIndex,tIndex,m);
2606  }
2607  const int n1InteriorCells = cells.getSelfCount();
2608  const int n1BoundaryCells = cells.getCount() - n1InterfaceCells - n1InteriorCells;
2609  totalBoundaryCells += n1BoundaryCells;
2610  totalInteriorCells += n1InteriorCells;
2611  }
2612 
2613  cout << totalBoundaryCells << endl;
2614  cout << totalInteriorCells << endl;
2615 
2616  const MultiField::ArrayIndex rowIndex = bField.getArrayIndex(0);
2617  const Matrix& origMatrix = mfMatrix.getMatrix(rowIndex,rowIndex);
2618  const CRConnectivity& origConn = origMatrix.getConnectivity();
2619  const ArrayBase& origB = bField[rowIndex];
2620 
2621  LinearSystem lsCombined;
2622 
2623  shared_ptr<CRConnectivity> trPtr(new CRConnectivity(origConn.getColSite(),origConn.getRowSite()));
2624  CRConnectivity& tr = *trPtr;
2625 
2626  tr.getrowDim *= varSize;
2627  tr._colDim *= varSize;
2628 
2629  tr.initCount();
2630 
2631  const Array<int>& myRow = *_row;
2632  const Array<int>& myCol = *_col;*/
2633 
2634 
2635  return ls;
2636  }
2637 
2638  void advanceSpecies(const int niter)
2639  {
2640  for(int n=0; n<niter; n++)
2641  {
2642  bool allConverged=true;
2643  for (int m=0; m<_nSpecies; m++)
2644  {
2645  MFRPtr& iNorm = *_initialSpeciesNormVector[m];
2646  MFRPtr& rCurrent = *_currentSpeciesResidual[m];
2647 
2648  LinearSystem ls,lsShell;
2649  //LinearSystem ls;
2650 
2651  initSpeciesLinearization(ls, lsShell, m);
2652  //initSpeciesLinearization(ls, ls, m);
2653 
2654  ls.initAssembly();
2655  lsShell.initAssembly();
2656 
2657  linearizeSpecies(ls, lsShell, m);
2658  //linearizeSpecies(ls, ls, m);
2659 
2660  ls.initSolve();
2661 
2662  //cout << "BEFORE LS" << endl;
2664  //printMatrixElementsOnFace(*_meshes[0], 8, ls, sFields.concentration);
2665  //printMatrixElementsOnFace(*_meshes[1], 11, ls, sFields.concentration);
2666 
2667 
2668  MFRPtr rNorm(_options.getLinearSolverSpecies().solve(ls));
2669 
2670  if (!iNorm) iNorm = rNorm;
2671  MFRPtr normRatio((*rNorm)/(*iNorm));
2672 
2673  //cout << "Species Number: " << m << endl;
2674  if (((_niters+1) % _options.advanceVerbosity)==0)
2675  cout << _niters << ": " << *rNorm << endl;
2676 
2677  _options.getLinearSolverSpecies().cleanup();
2678 
2679  ls.postSolve();
2680  ls.updateSolution();
2681 
2682  //cout << "AFTER LS" << endl;
2683  //printMatrixElementsOnFace(*_meshes[0], 8, ls, sFields.concentration);
2684  //printMatrixElementsOnFace(*_meshes[1], 11, ls, sFields.concentration);
2685 
2686 
2687  //updateShellGhosts(); not useful right now
2688  lsShell.initSolve();
2689  MFRPtr rNorm2(_options.getLinearSolverSpecies().solve(lsShell));
2690  _options.getLinearSolverSpecies().cleanup();
2691  lsShell.postSolve();
2692  lsShell.updateSolution();
2693  cout << _niters << ": " << *rNorm2 << endl;
2694 
2695 
2696 
2697  _niters++;
2698 
2699  rCurrent = rNorm;
2700 
2701  if (!(*rNorm < _options.absoluteSpeciesTolerance ||
2702  *normRatio < _options.relativeSpeciesTolerance))
2703  allConverged=false;
2704  }
2705  if (allConverged)
2706  break;
2707  }
2708  }
2709 
2710 void advancePotential(const int niter)
2711  {
2712  for(int n=0; n<niter; n++)
2713  {
2714  LinearSystem ls;
2715 
2716  initPotentialLinearization(ls);
2717 
2718  ls.initAssembly();
2719 
2720  linearizePotential(ls);
2721 
2722  ls.initSolve();
2723 
2724  MFRPtr rNorm(_options.getLinearSolverPotential().solve(ls));
2725 
2726  if (!_initialPotentialNorm) _initialPotentialNorm = rNorm;
2727 
2728  _currentPotentialResidual = rNorm;
2729 
2730  MFRPtr normRatio((*rNorm)/(*_initialPotentialNorm));
2731 
2732  if ((_niters % _options.advanceVerbosity)==0)
2733  cout << "Potential: " << _niters << ": " << *rNorm << endl;
2734 
2735  _options.getLinearSolverPotential().cleanup();
2736 
2737  ls.postSolve();
2738  ls.updateSolution();
2739 
2740  _niters++;
2741 
2742  if (*rNorm < _options.absolutePotentialTolerance ||
2743  *normRatio < _options.relativePotentialTolerance)
2744  break;
2745  }
2746  }
2747 
2748 void advanceThermal(const int niter)
2749  {
2750  for(int n=0; n<niter; n++)
2751  {
2752  LinearSystem ls;
2753 
2754  initThermalLinearization(ls);
2755 
2756  ls.initAssembly();
2757 
2758  linearizeThermal(ls);
2759 
2760  ls.initSolve();
2761 
2762  MFRPtr rNorm(_options.getLinearSolverThermal().solve(ls));
2763 
2764  if (!_initialThermalNorm) _initialThermalNorm = rNorm;
2765 
2766  _currentThermalResidual = rNorm;
2767 
2768  MFRPtr normRatio((*rNorm)/(*_initialThermalNorm));
2769 
2770  if ((_niters % _options.advanceVerbosity)==0)
2771  cout << "Thermal: " << _niters << ": " << *rNorm << endl;
2772 
2773  _options.getLinearSolverThermal().cleanup();
2774 
2775  ls.postSolve();
2776  ls.updateSolution();
2777 
2778  _niters++;
2779 
2780  if (*rNorm < _options.absoluteThermalTolerance ||
2781  *normRatio < _options.relativeThermalTolerance)
2782  break;
2783  }
2784  }
2785 
2786  void advanceCoupled(const int niter)
2787  {
2788  for(int n=0; n<niter; n++)
2789  {
2790 
2791  LinearSystem ls;
2792 
2793  initPCLinearization(ls);
2794 
2795  ls.initAssembly();
2796 
2797  if (_options.thermalModelPC)
2798  linearizePC_Thermal(ls);
2799  else
2800  linearizePC(ls);
2801 
2802  ls.initSolve();
2803 
2804  MFRPtr rNorm = _options.getLinearSolverPC().solve(ls);
2805 
2806  if (!_initialPCNorm) _initialPCNorm = rNorm;
2807 
2808  _currentPCResidual = rNorm;
2809 
2810  MFRPtr normRatio((*rNorm)/(*_initialPCNorm));
2811 
2812  _options.getLinearSolverPC().cleanup();
2813 
2814  ls.postSolve();
2815 
2816  ls.updateSolution();
2817 
2818  if ((_niters % _options.advanceVerbosity)==0)
2819  cout << "Point-Coupled: " << _niters << ": " << *rNorm << endl;
2820 
2821  _niters++;
2822 
2823  if (*rNorm < _options.absolutePCTolerance ||
2824  *normRatio < _options.relativePCTolerance)
2825  break;
2826  }
2827  }
2828 
2829  T getSpeciesResidual(const int speciesId)
2830  {
2831  MFRPtr& rCurrent = *_currentSpeciesResidual[speciesId];
2832  BatterySpeciesFields& sFields = *_speciesFieldsVector[speciesId];
2833  const Field& concentrationField = sFields.concentration;
2834  ArrayBase& residualArray = (*rCurrent)[concentrationField];
2835  T *arrayData = (T*)(residualArray.getData());
2836  const T residual = arrayData[0];
2837  return residual;
2838  }
2839 
2840  T getPotentialResidual()
2841  {
2842  const Field& potentialField = _batteryModelFields.potential;
2843  ArrayBase& residualArray = (*_currentPotentialResidual)[potentialField];
2844  T *arrayData = (T*)(residualArray.getData());
2845  const T residual = arrayData[0];
2846  return residual;
2847  }
2848 
2849  T getThermalResidual()
2850  {
2851  const Field& thermalField = _batteryModelFields.temperature;
2852  ArrayBase& residualArray = (*_currentThermalResidual)[thermalField];
2853  T *arrayData = (T*)(residualArray.getData());
2854  const T residual = arrayData[0];
2855  return residual;
2856  }
2857 
2858  T getPCResidual(const int v)
2859  {
2860  const Field& potentialSpeciesTempField = _batteryModelFields.potentialSpeciesTemp;
2861  ArrayBase& residualArray = (*_currentPCResidual)[potentialSpeciesTempField];
2862  T *arrayData = (T*)(residualArray.getData());
2863  if (v<2)
2864  {
2865  const T residual = arrayData[v];
2866  return residual;
2867  }
2868  else if ((v==2)&&(_options.thermalModelPC))
2869  {
2870  const T residual = arrayData[v];
2871  return residual;
2872  }
2873  else
2874  {
2875  const T residual = 0.0;
2876  return residual;
2877  }
2878  }
2879 
2880  void updateShellGhosts()
2881  {
2882  const int numMeshes = _meshes.size();
2883  for (int n=0; n<numMeshes; n++)
2884  {
2885  const Mesh& mesh = *_meshes[n];
2886  if ((mesh.isDoubleShell())&&(!(mesh.isConnectedShell())))
2887  {
2888  BatterySpeciesFields& sFields = *_speciesFieldsVector[0]; //first species is lithium
2889 
2890  const StorageSite& cells = mesh.getCells();
2891  const StorageSite& faces = mesh.getParentFaceGroupSite();
2892  const CRConnectivity& cellCells = mesh.getCellCells();
2893  TArray& speciesCell = dynamic_cast<TArray&>(sFields.concentration[cells]);
2894 
2895  const int parentMeshID = mesh.getParentMeshID();
2896  const int otherMeshID = mesh.getOtherMeshID();
2897  const Mesh& parentMesh = *_meshes[parentMeshID];
2898  const Mesh& otherMesh = *_meshes[otherMeshID];
2899 
2900  const CRConnectivity& parentFaceCells = parentMesh.getFaceCells(faces);
2901  const StorageSite& parentCells = parentMesh.getCells();
2902  const StorageSite& otherFaces = mesh.getOtherFaceGroupSite();
2903  const CRConnectivity& otherFaceCells = otherMesh.getFaceCells(otherFaces);
2904  const StorageSite& otherCells = otherMesh.getCells();
2905  const TArray& speciesCellParent = dynamic_cast< const TArray&>(sFields.concentration[parentCells]);
2906  const TArray& speciesCellOther = dynamic_cast< const TArray&>(sFields.concentration[otherCells]);
2907 
2908 
2909  for (int f=0; f<faces.getCount(); f++)
2910  {
2911  int c0p = parentFaceCells(f,0);
2912  int c1p = parentFaceCells(f,1);
2913  if (c1p < parentCells.getSelfCount())
2914  {
2915  int temp = c0p;
2916  c0p = c1p;
2917  c1p = temp;
2918  }
2919 
2920  int c0o = otherFaceCells(f,0);
2921  int c1o = otherFaceCells(f,1);
2922  if (c1o < otherCells.getSelfCount())
2923  {
2924  int temp = c0o;
2925  c0o = c1o;
2926  c1o = temp;
2927  }
2928 
2929  const int c0 = f;
2930  const int c1 = cellCells(f,0);
2931  const int c2 = cellCells(f,1);
2932  const int c3 = cellCells(f,2);
2933 
2934  if (c0==0)
2935  {
2936  cout << speciesCell[c0] << " " << speciesCell[c1] << " " << speciesCell[c2] << " " << speciesCell[c3] << endl;
2937  }
2938 
2939  // copy sln varialbe values from interior cells of meshes to ghost cells of shell mesh
2940  speciesCell[c3] = speciesCellParent[c0p];
2941  speciesCell[c2] = speciesCellOther[c0o];
2942 
2943  if (c0==0)
2944  {
2945  cout << speciesCell[c0] << " " << speciesCell[c1] << " " << speciesCell[c2] << " " << speciesCell[c3] << endl;
2946  }
2947  }
2948  }
2949  }
2950  }
2951 
2952  void copySeparateToCoupled()
2953  {
2954  const int numMeshes = _meshes.size();
2955  for (int n=0; n<numMeshes; n++)
2956  {
2957  //copy concentration and potential and thermal
2958  const Mesh& mesh = *_meshes[n];
2959  const StorageSite& cells = mesh.getCells();
2960  BatterySpeciesFields& sFields = *_speciesFieldsVector[0]; //first species is lithium
2961  const TArray& mFSeparate = dynamic_cast<const TArray&>(sFields.concentration[cells]);
2962  const TArray& pSeparate = dynamic_cast<const TArray&>(_batteryModelFields.potential[cells]);
2963  const TArray& tSeparate = dynamic_cast<const TArray&>(_batteryModelFields.temperature[cells]);
2964 
2965  if (_options.thermalModelPC)
2966  {
2967  VectorT3Array& coupledValues = dynamic_cast<VectorT3Array&>(_batteryModelFields.potentialSpeciesTemp[cells]);
2968  for (int c=0; c<cells.getCount(); c++)
2969  {
2970  T concentration = mFSeparate[c];
2971  T potential = pSeparate[c];
2972  T temperature = tSeparate[c];
2973  VectorT3& coupledCellVector = coupledValues[c];
2974 
2975  //populate coupled field with separate field data
2976  coupledCellVector[0] = potential;
2977  coupledCellVector[1] = concentration;
2978  coupledCellVector[2] = temperature;
2979  }
2980  }
2981  else
2982  {
2983  VectorT2Array& coupledValues = dynamic_cast<VectorT2Array&>(_batteryModelFields.potentialSpeciesTemp[cells]);
2984  for (int c=0; c<cells.getCount(); c++)
2985  {
2986  T concentration = mFSeparate[c];
2987  T potential = pSeparate[c];
2988  VectorT2& coupledCellVector = coupledValues[c];
2989 
2990  //populate coupled field with separate field data
2991  coupledCellVector[0] = potential;
2992  coupledCellVector[1] = concentration;
2993  }
2994  }
2995 
2996  //copy fluxes
2997  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
2998  {
2999  const FaceGroup& fg = *fgPtr;
3000  const StorageSite& faces = fg.site;
3001  const TArray& mFFluxSeparate = dynamic_cast<const TArray&>(sFields.massFlux[faces]);
3002  const TArray& pFluxSeparate = dynamic_cast<const TArray&>(_batteryModelFields.potential_flux[faces]);
3003  const TArray& tFluxSeparate = dynamic_cast<const TArray&>(_batteryModelFields.heatFlux[faces]);
3004  if (_options.thermalModelPC)
3005  {
3006  VectorT3Array& coupledFluxValues = dynamic_cast<VectorT3Array&>(_batteryModelFields.potentialSpeciesTempFlux[faces]);
3007  for (int f=0; f<faces.getCount(); f++)
3008  {
3009  T concentrationFlux = mFFluxSeparate[f];
3010  T potentialFlux = pFluxSeparate[f];
3011  T heatFlux = tFluxSeparate[f];
3012  VectorT3& coupledCellFluxVector = coupledFluxValues[f];
3013 
3014  //populate coupled field with separate field data
3015  coupledCellFluxVector[0] = potentialFlux;
3016  coupledCellFluxVector[1] = concentrationFlux;
3017  coupledCellFluxVector[2] = heatFlux;
3018  }
3019  }
3020  else
3021  {
3022  VectorT2Array& coupledFluxValues = dynamic_cast<VectorT2Array&>(_batteryModelFields.potentialSpeciesTempFlux[faces]);
3023  for (int f=0; f<faces.getCount(); f++)
3024  {
3025  T concentrationFlux = mFFluxSeparate[f];
3026  T potentialFlux = pFluxSeparate[f];
3027  VectorT2& coupledCellFluxVector = coupledFluxValues[f];
3028 
3029  //populate coupled field with separate field data
3030  coupledCellFluxVector[0] = potentialFlux;
3031  coupledCellFluxVector[1] = concentrationFlux;
3032  }
3033  }
3034  }
3035  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
3036  {
3037  const FaceGroup& fg = *fgPtr;
3038  const StorageSite& faces = fg.site;
3039  const TArray& mFFluxSeparate = dynamic_cast<const TArray&>(sFields.massFlux[faces]);
3040  const TArray& pFluxSeparate = dynamic_cast<const TArray&>(_batteryModelFields.potential_flux[faces]);
3041  const TArray& tFluxSeparate = dynamic_cast<const TArray&>(_batteryModelFields.heatFlux[faces]);
3042  if (_options.thermalModelPC)
3043  {
3044  VectorT3Array& coupledFluxValues = dynamic_cast<VectorT3Array&>(_batteryModelFields.potentialSpeciesTempFlux[faces]);
3045  for (int f=0; f<faces.getCount(); f++)
3046  {
3047  T concentrationFlux = mFFluxSeparate[f];
3048  T potentialFlux = pFluxSeparate[f];
3049  T heatFlux = tFluxSeparate[f];
3050  VectorT3& coupledCellFluxVector = coupledFluxValues[f];
3051 
3052  //populate coupled field with separate field data
3053  coupledCellFluxVector[0] = potentialFlux;
3054  coupledCellFluxVector[1] = concentrationFlux;
3055  coupledCellFluxVector[2] = heatFlux;
3056  }
3057  }
3058  else
3059  {
3060  VectorT2Array& coupledFluxValues = dynamic_cast<VectorT2Array&>(_batteryModelFields.potentialSpeciesTempFlux[faces]);
3061  for (int f=0; f<faces.getCount(); f++)
3062  {
3063  T concentrationFlux = mFFluxSeparate[f];
3064  T potentialFlux = pFluxSeparate[f];
3065  VectorT2& coupledCellFluxVector = coupledFluxValues[f];
3066 
3067  //populate coupled field with separate field data
3068  coupledCellFluxVector[0] = potentialFlux;
3069  coupledCellFluxVector[1] = concentrationFlux;
3070  }
3071  }
3072  }
3073 
3074  // copy to N1 and N2 if transient
3075  if (_options.transient)
3076  {
3077  const TArray& potSeparateN1 = dynamic_cast<const TArray&>(_batteryModelFields.potentialN1[cells]);
3078  const TArray& mFSeparateN1 = dynamic_cast<const TArray&>(sFields.concentrationN1[cells]);
3079  const TArray& tSeparateN1 = dynamic_cast<const TArray&>(_batteryModelFields.temperatureN1[cells]);
3080  if (_options.thermalModelPC)
3081  {
3082  VectorT3Array& coupledValuesN1 = dynamic_cast<VectorT3Array&>(_batteryModelFields.potentialSpeciesTempN1[cells]);
3083 
3084  for (int c=0; c<cells.getCount(); c++)
3085  {
3086  const T potentialN1 = potSeparateN1[c];
3087  const T concentrationN1 = mFSeparateN1[c];
3088  const T temperatureN1 = tSeparateN1[c];
3089  VectorT3& coupledCellVectorN1 = coupledValuesN1[c];
3090 
3091  //populate coupled field with separate field data
3092  coupledCellVectorN1[1] = concentrationN1;
3093  coupledCellVectorN1[0] = potentialN1;
3094  coupledCellVectorN1[2] = temperatureN1;
3095  }
3096 
3097  if (_options.timeDiscretizationOrder > 1)
3098  {
3099  VectorT3Array& coupledValuesN2 = dynamic_cast<VectorT3Array&>(_batteryModelFields.potentialSpeciesTempN2[cells]);
3100  const TArray& mFSeparateN2 = dynamic_cast<const TArray&>(sFields.concentrationN2[cells]);
3101  const TArray& tSeparateN2 = dynamic_cast<const TArray&>(_batteryModelFields.temperatureN2[cells]);
3102  for (int c=0; c<cells.getCount(); c++)
3103  {
3104  const T concentrationN2 = mFSeparateN2[c];
3105  const T temperatureN2 = tSeparateN2[c];
3106  VectorT3& coupledCellVectorN2 = coupledValuesN2[c];
3107 
3108  //populate coupled field with separate field data
3109  coupledCellVectorN2[1] = concentrationN2;
3110  coupledCellVectorN2[2] = temperatureN2;
3111  }
3112  }
3113  }
3114  else
3115  {
3116  VectorT2Array& coupledValuesN1 = dynamic_cast<VectorT2Array&>(_batteryModelFields.potentialSpeciesTempN1[cells]);
3117 
3118  for (int c=0; c<cells.getCount(); c++)
3119  {
3120  const T potentialN1 = potSeparateN1[c];
3121  const T concentrationN1 = mFSeparateN1[c];
3122  VectorT2& coupledCellVectorN1 = coupledValuesN1[c];
3123 
3124  //populate coupled field with separate field data
3125  coupledCellVectorN1[1] = concentrationN1;
3126  coupledCellVectorN1[0] = potentialN1;
3127  }
3128 
3129  if (_options.timeDiscretizationOrder > 1)
3130  {
3131  VectorT2Array& coupledValuesN2 = dynamic_cast<VectorT2Array&>(_batteryModelFields.potentialSpeciesTempN2[cells]);
3132  const TArray& mFSeparateN2 = dynamic_cast<const TArray&>(sFields.concentrationN2[cells]);
3133  for (int c=0; c<cells.getCount(); c++)
3134  {
3135  const T concentrationN2 = mFSeparateN2[c];
3136  VectorT2& coupledCellVectorN2 = coupledValuesN2[c];
3137 
3138  //populate coupled field with separate field data
3139  coupledCellVectorN2[1] = concentrationN2;
3140  }
3141  }
3142  }
3143  }
3144  }
3145  }
3146 
3147  void copyCoupledToSeparate()
3148  {
3149  const int numMeshes = _meshes.size();
3150  for (int n=0; n<numMeshes; n++)
3151  {
3152  const Mesh& mesh = *_meshes[n];
3153  const StorageSite& cells = mesh.getCells();
3154  BatterySpeciesFields& sFields = *_speciesFieldsVector[0]; //first species is lithium
3155 
3156  // copy conc and potential
3157  TArray& mFSeparate = dynamic_cast<TArray&>(sFields.concentration[cells]);
3158  TArray& pSeparate = dynamic_cast<TArray&>(_batteryModelFields.potential[cells]);
3159  TArray& tSeparate = dynamic_cast<TArray&>(_batteryModelFields.temperature[cells]);
3160  if (_options.thermalModelPC)
3161  {
3162  const VectorT3Array& coupledValues = dynamic_cast<const VectorT3Array&>(_batteryModelFields.potentialSpeciesTemp[cells]);
3163  for (int c=0; c<cells.getCount(); c++)
3164  {
3165  VectorT3 coupledCellVector = coupledValues[c];
3166  T potential = coupledCellVector[0];
3167  T concentration = coupledCellVector[1];
3168  T temperature = coupledCellVector[2];
3169 
3170  //populate separate fields with coupled field data
3171  pSeparate[c] = potential;
3172  mFSeparate[c] = concentration;
3173  tSeparate[c] = temperature;
3174  }
3175  }
3176  else
3177  {
3178  const VectorT2Array& coupledValues = dynamic_cast<const VectorT2Array&>(_batteryModelFields.potentialSpeciesTemp[cells]);
3179  for (int c=0; c<cells.getCount(); c++)
3180  {
3181  VectorT2 coupledCellVector = coupledValues[c];
3182  T potential = coupledCellVector[0];
3183  T concentration = coupledCellVector[1];
3184 
3185  //populate separate fields with coupled field data
3186  pSeparate[c] = potential;
3187  mFSeparate[c] = concentration;
3188  }
3189  }
3190  // copy fluxes
3191  foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
3192  {
3193  const FaceGroup& fg = *fgPtr;
3194  const StorageSite& faces = fg.site;
3195  TArray& mFFluxSeparate = dynamic_cast<TArray&>(sFields.massFlux[faces]);
3196  TArray& pFluxSeparate = dynamic_cast<TArray&>(_batteryModelFields.potential_flux[faces]);
3197  TArray& tFluxSeparate = dynamic_cast<TArray&>(_batteryModelFields.heatFlux[faces]);
3198  if (_options.thermalModelPC)
3199  {
3200  const VectorT3Array& coupledFluxValues = dynamic_cast<const VectorT3Array&>(_batteryModelFields.potentialSpeciesTempFlux[faces]);
3201  for (int f=0; f<faces.getCount(); f++)
3202  {
3203  VectorT3 coupledCellFluxVector = coupledFluxValues[f];
3204  T potentialFlux = coupledCellFluxVector[0];
3205  T concentrationFlux = coupledCellFluxVector[1];
3206  T heatFlux = coupledCellFluxVector[2];
3207 
3208  //populate separate fields with coupled field data
3209  pFluxSeparate[f] = potentialFlux;
3210  mFFluxSeparate[f] = concentrationFlux;
3211  tFluxSeparate[f] = heatFlux;
3212  }
3213  }
3214  else
3215  {
3216  const VectorT2Array& coupledFluxValues = dynamic_cast<const VectorT2Array&>(_batteryModelFields.potentialSpeciesTempFlux[faces]);
3217  for (int f=0; f<faces.getCount(); f++)
3218  {
3219  VectorT2 coupledCellFluxVector = coupledFluxValues[f];
3220  T potentialFlux = coupledCellFluxVector[0];
3221  T concentrationFlux = coupledCellFluxVector[1];
3222 
3223  //populate separate fields with coupled field data
3224  pFluxSeparate[f] = potentialFlux;
3225  mFFluxSeparate[f] = concentrationFlux;
3226  }
3227  }
3228  }
3229  foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
3230  {
3231  const FaceGroup& fg = *fgPtr;
3232  const StorageSite& faces = fg.site;
3233  TArray& mFFluxSeparate = dynamic_cast<TArray&>(sFields.massFlux[faces]);
3234  TArray& pFluxSeparate = dynamic_cast<TArray&>(_batteryModelFields.potential_flux[faces]);
3235  TArray& tFluxSeparate = dynamic_cast<TArray&>(_batteryModelFields.heatFlux[faces]);
3236  if (_options.thermalModelPC)
3237  {
3238  const VectorT3Array& coupledFluxValues = dynamic_cast<const VectorT3Array&>(_batteryModelFields.potentialSpeciesTempFlux[faces]);
3239  for (int f=0; f<faces.getCount(); f++)
3240  {
3241  VectorT3 coupledCellFluxVector = coupledFluxValues[f];
3242  T potentialFlux = coupledCellFluxVector[0];
3243  T concentrationFlux = coupledCellFluxVector[1];
3244  T heatFlux = coupledCellFluxVector[2];
3245 
3246  //populate separate fields with coupled field data
3247  pFluxSeparate[f] = potentialFlux;
3248  mFFluxSeparate[f] = concentrationFlux;
3249  tFluxSeparate[f] = heatFlux;
3250 
3251  }
3252  }
3253  else
3254  {
3255  const VectorT2Array& coupledFluxValues = dynamic_cast<const VectorT2Array&>(_batteryModelFields.potentialSpeciesTempFlux[faces]);
3256  for (int f=0; f<faces.getCount(); f++)
3257  {
3258  VectorT2 coupledCellFluxVector = coupledFluxValues[f];
3259  T potentialFlux = coupledCellFluxVector[0];
3260  T concentrationFlux = coupledCellFluxVector[1];
3261 
3262  //populate separate fields with coupled field data
3263  pFluxSeparate[f] = potentialFlux;
3264  mFFluxSeparate[f] = concentrationFlux;
3265  }
3266  }
3267  }
3268 
3269  // copy to N1 and N2 if transient
3270  if (_options.transient)
3271  {
3272  TArray& potSeparateN1 = dynamic_cast<TArray&>(_batteryModelFields.potentialN1[cells]);
3273  TArray& mFSeparateN1 = dynamic_cast<TArray&>(sFields.concentrationN1[cells]);
3274  TArray& tSeparateN1 = dynamic_cast<TArray&>(_batteryModelFields.temperatureN1[cells]);
3275 
3276  if (_options.thermalModelPC)
3277  {
3278  const VectorT3Array& coupledValuesN1 = dynamic_cast< const VectorT3Array&>(_batteryModelFields.potentialSpeciesTempN1[cells]);
3279 
3280  for (int c=0; c<cells.getCount(); c++)
3281  {
3282  const VectorT3& coupledCellVectorN1 = coupledValuesN1[c];
3283 
3284  //populate separate field with coupled field data
3285  mFSeparateN1[c] = coupledCellVectorN1[1];
3286  potSeparateN1[c] = coupledCellVectorN1[0];
3287  tSeparateN1[c] = coupledCellVectorN1[2];
3288  }
3289 
3290  if (_options.timeDiscretizationOrder > 1)
3291  {
3292  const VectorT3Array& coupledValuesN2 = dynamic_cast<const VectorT3Array&>(_batteryModelFields.potentialSpeciesTempN2[cells]);
3293  TArray& mFSeparateN2 = dynamic_cast<TArray&>(sFields.concentrationN2[cells]);
3294  TArray& tSeparateN2 = dynamic_cast<TArray&>(_batteryModelFields.temperatureN2[cells]);
3295  for (int c=0; c<cells.getCount(); c++)
3296  {
3297  const VectorT3& coupledCellVectorN2 = coupledValuesN2[c];
3298 
3299  //populate separate field with coupled field data
3300  mFSeparateN2[c] = coupledCellVectorN2[1];
3301  tSeparateN2[c] = coupledCellVectorN2[2];
3302  }
3303  }
3304  }
3305  else
3306  {
3307  const VectorT2Array& coupledValuesN1 = dynamic_cast< const VectorT2Array&>(_batteryModelFields.potentialSpeciesTempN1[cells]);
3308 
3309  for (int c=0; c<cells.getCount(); c++)
3310  {
3311  const VectorT2& coupledCellVectorN1 = coupledValuesN1[c];
3312 
3313  //populate separate field with coupled field data
3314  mFSeparateN1[c] = coupledCellVectorN1[1];
3315  potSeparateN1[c] = coupledCellVectorN1[0];
3316  }
3317 
3318  if (_options.timeDiscretizationOrder > 1)
3319  {
3320  const VectorT2Array& coupledValuesN2 = dynamic_cast<const VectorT2Array&>(_batteryModelFields.potentialSpeciesTempN2[cells]);
3321  TArray& mFSeparateN2 = dynamic_cast<TArray&>(sFields.concentrationN2[cells]);
3322  for (int c=0; c<cells.getCount(); c++)
3323  {
3324  const VectorT2& coupledCellVectorN2 = coupledValuesN2[c];
3325 
3326  //populate separate field with coupled field data
3327  mFSeparateN2[c] = coupledCellVectorN2[1];
3328  }
3329  }
3330  }
3331  }
3332  }
3333  }
3334 
3335 void copyPCDiffusivity()
3336  {
3337  const int numMeshes = _meshes.size();
3338  for (int n=0; n<numMeshes; n++)
3339  {
3340  //copy concentration diff and potential cond and heat thermal cond
3341  const Mesh& mesh = *_meshes[n];
3342  const StorageSite& cells = mesh.getCells();
3343  BatterySpeciesFields& sFields = *_speciesFieldsVector[0]; //first species is lithium
3344  const TArray& mFDiffSeparate = dynamic_cast<const TArray&>(sFields.diffusivity[cells]);
3345  const TArray& pDiffSeparate = dynamic_cast<const TArray&>(_batteryModelFields.conductivity[cells]);
3346  const TArray& tDiffSeparate = dynamic_cast<const TArray&>(_batteryModelFields.thermalConductivity[cells]);
3347 
3348  if (_options.thermalModelPC)
3349  {
3350  VectorT3Array& coupledDiffValues = dynamic_cast<VectorT3Array&>(_batteryModelFields.potentialSpeciesTempDiffusivity[cells]);
3351 
3352  for (int c=0; c<cells.getCount(); c++)
3353  {
3354  T massDiff = mFDiffSeparate[c];
3355  T potentialDiff = pDiffSeparate[c];
3356  T thermalDiff = tDiffSeparate[c];
3357  VectorT3& coupledCellDiffVector = coupledDiffValues[c];
3358 
3359  //populate coupled field with separate field data
3360  coupledCellDiffVector[0] = potentialDiff;
3361  coupledCellDiffVector[1] = massDiff;
3362  coupledCellDiffVector[2] = thermalDiff;
3363  }
3364  }
3365  else
3366  {
3367  VectorT2Array& coupledDiffValues = dynamic_cast<VectorT2Array&>(_batteryModelFields.potentialSpeciesTempDiffusivity[cells]);
3368 
3369  for (int c=0; c<cells.getCount(); c++)
3370  {
3371  T massDiff = mFDiffSeparate[c];
3372  T potentialDiff = pDiffSeparate[c];
3373  VectorT2& coupledCellDiffVector = coupledDiffValues[c];
3374 
3375  //populate coupled field with separate field data
3376  coupledCellDiffVector[0] = potentialDiff;
3377  coupledCellDiffVector[1] = massDiff;
3378  }
3379  }
3380  }
3381  }
3382 
3383 T getFaceGroupArea(const Mesh& mesh, const int fgID)
3384  {
3385  foreach(const FaceGroupPtr fgPtr, mesh.getAllFaceGroups())
3386  {
3387  const FaceGroup& fg = *fgPtr;
3388  if (fg.id == fgID)
3389  {
3390  const StorageSite& faces = fg.site;
3391  const int nFaces = faces.getCount();
3392  const TArray& faceAreaMag = dynamic_cast<const TArray&>(_geomFields.areaMag[faces]);
3393  T totalArea = 0.0;
3394  for(int f=0; f<nFaces; f++)
3395  {
3396  totalArea += faceAreaMag[f];
3397  if (f%1000==0)
3398  cout << "Percent Done: " << (f*100.0/nFaces) << "%" << endl;
3399  }
3400  return totalArea;
3401  }
3402  }
3403  throw CException("getFaceGroupArea: No face group with given id");
3404  }
3405 
3406 T getFaceGroupVoltage(const Mesh& mesh, const int fgID)
3407  {
3408  foreach(const FaceGroupPtr fgPtr, mesh.getAllFaceGroups())
3409  {
3410  const FaceGroup& fg = *fgPtr;
3411  if (fg.id == fgID)
3412  {
3413  const StorageSite& faces = fg.site;
3414  const int nFaces = faces.getCount();
3415  const TArray& faceAreaMag = dynamic_cast<const TArray&>(_geomFields.areaMag[faces]);
3416  const StorageSite& cells = mesh.getCells();
3417  const CRConnectivity& faceCells = mesh.getFaceCells(faces);
3418  const TArray& cellPotential = dynamic_cast<const TArray&>(_batteryModelFields.potential[cells]);
3419  T totalArea = 0.0;
3420  T voltageSum = 0.0;
3421  for(int f=0; f<nFaces; f++)
3422  {
3423  //const T faceVoltage = harmonicAverage(cellPotential[faceCells(f,0)],cellPotential[faceCells(f,1)]);
3424  const T faceVoltage = cellPotential[faceCells(f,1)];
3425  totalArea += faceAreaMag[f];
3426  voltageSum += faceVoltage*faceAreaMag[f];
3427  }
3428  return (voltageSum/totalArea);
3429  }
3430  }
3431  throw CException("getFaceGroupVoltage: No face group with given id");
3432  }
3433 
3434 T getMeshVolume(const Mesh& mesh)
3435  {
3436  const StorageSite& cells = mesh.getCells();
3437  const TArray& cellVolume = dynamic_cast<const TArray&>(_geomFields.volume[cells]);
3438  T totalVolume = 0.0;
3439  for(int c=0; c<cells.getSelfCount(); c++)
3440  {
3441  totalVolume += cellVolume[c];
3442  }
3443  return totalVolume;
3444  }
3445 
3446 private:
3447  const MeshList _realMeshes;
3448  const MeshList _meshes;
3449  GeomFields& _geomFields;
3450 
3451  vector<BatterySpeciesFields*> _speciesFieldsVector;
3452 
3453  vector<BatterySpeciesBCMap*> _sbcMapVector;
3454  vector<BatterySpeciesVCMap*> _svcMapVector;
3455  BatteryPotentialBCMap _pbcMap;
3456  BatteryPotentialVCMap _pvcMap;
3457  BatteryThermalBCMap _tbcMap;
3458  BatteryThermalVCMap _tvcMap;
3459 
3460  BatteryModelOptions<T> _options;
3461 
3462  MFRPtr _initialPotentialNorm;
3463  MFRPtr _initialThermalNorm;
3464  MFRPtr _initialPCNorm;
3465  vector<MFRPtr*> _initialSpeciesNormVector;
3466  int _niters;
3467  const int _nSpecies;
3468 
3469  BatteryModelFields _batteryModelFields;
3470 
3471  MFRPtr _currentPotentialResidual;
3472  MFRPtr _currentThermalResidual;
3473  MFRPtr _currentPCResidual;
3474  vector<MFRPtr*> _currentSpeciesResidual;
3475 };
3476 
3477 template<class T>
3478 BatteryModel<T>::BatteryModel(GeomFields& geomFields,
3479  const MeshList& realMeshes,
3480  const MeshList& meshes,
3481  const int nSpecies) :
3482  Model(meshes),
3483  _impl(new Impl(geomFields,realMeshes,meshes,nSpecies))
3484 {
3485  logCtor();
3486 }
3487 
3488 template<class T>
3489 BatteryModel<T>::~BatteryModel()
3490 {
3491  logDtor();
3492 }
3493 
3494 template<class T>
3495 void
3496 BatteryModel<T>::init()
3497 {
3498  _impl->init();
3499 }
3500 template<class T>
3501 BatterySpeciesFields&
3502 BatteryModel<T>::getBatterySpeciesFields(const int speciesId) {return _impl->getBatterySpeciesFields(speciesId);}
3503 
3504 template<class T>
3505 typename BatteryModel<T>::BatterySpeciesBCMap&
3506 BatteryModel<T>::getSpeciesBCMap(const int speciesId) {return _impl->getSpeciesBCMap(speciesId);}
3507 
3508 template<class T>
3509 typename BatteryModel<T>::BatterySpeciesVCMap&
3510 BatteryModel<T>::getSpeciesVCMap(const int speciesId) {return _impl->getSpeciesVCMap(speciesId);}
3511 
3512 template<class T>
3513 BatteryModelFields&
3514 BatteryModel<T>::getBatteryModelFields() {return _impl->getBatteryModelFields();}
3515 
3516 template<class T>
3517 typename BatteryModel<T>::BatteryPotentialBCMap&
3518 BatteryModel<T>::getPotentialBCMap() {return _impl->getPotentialBCMap();}
3519 
3520 template<class T>
3521 typename BatteryModel<T>::BatteryPotentialVCMap&
3522 BatteryModel<T>::getPotentialVCMap() {return _impl->getPotentialVCMap();}
3523 
3524 template<class T>
3525 typename BatteryModel<T>::BatteryThermalBCMap&
3526 BatteryModel<T>::getThermalBCMap() {return _impl->getThermalBCMap();}
3527 
3528 template<class T>
3529 typename BatteryModel<T>::BatteryThermalVCMap&
3530 BatteryModel<T>::getThermalVCMap() {return _impl->getThermalVCMap();}
3531 
3532 template<class T>
3533 BatteryModelOptions<T>&
3534 BatteryModel<T>::getOptions() {return _impl->getOptions();}
3535 
3536 template<class T>
3537 void
3538 BatteryModel<T>::advanceSpecies(const int niter)
3539 {
3540  _impl->advanceSpecies(niter);
3541 }
3542 
3543 template<class T>
3544 void
3545 BatteryModel<T>::advancePotential(const int niter)
3546 {
3547  _impl->advancePotential(niter);
3548 }
3549 
3550 template<class T>
3551 void
3552 BatteryModel<T>::advanceThermal(const int niter)
3553 {
3554  _impl->advanceThermal(niter);
3555 }
3556 
3557 template<class T>
3558 void
3559 BatteryModel<T>::advanceCoupled(const int niter)
3560 {
3561  _impl->advanceCoupled(niter);
3562 }
3563 
3564 template<class T>
3565 void
3566 BatteryModel<T>::updateTime()
3567 {
3568  _impl->updateTime();
3569 }
3570 
3571 template<class T>
3572 void
3573 BatteryModel<T>::recoverLastTimestep()
3574 {
3575  _impl->recoverLastTimestep();
3576 }
3577 
3578 
3579 template<class T>
3580 T
3581 BatteryModel<T>::getMassFluxIntegral(const Mesh& mesh, const int faceGroupId, const int m)
3582 {
3583  return _impl->getMassFluxIntegral(mesh, faceGroupId, m);
3584 }
3585 
3586 template<class T>
3587 T
3588 BatteryModel<T>::getPotentialFluxIntegral(const Mesh& mesh, const int faceGroupId)
3589 {
3590  return _impl->getPotentialFluxIntegral(mesh, faceGroupId);
3591 }
3592 
3593 template<class T>
3594 T
3595 BatteryModel<T>::getHeatFluxIntegral(const Mesh& mesh, const int faceGroupId)
3596 {
3597  return _impl->getHeatFluxIntegral(mesh, faceGroupId);
3598 }
3599 
3600 template<class T>
3601 T
3602 BatteryModel<T>::getAverageConcentration(const Mesh& mesh, const int m)
3603 {
3604  return _impl->getAverageConcentration(mesh, m);
3605 }
3606 
3607 template<class T>
3608 T
3609 BatteryModel<T>::getSpeciesResidual(const int speciesId)
3610 {
3611  return _impl->getSpeciesResidual(speciesId);
3612 }
3613 
3614 template<class T>
3615 T
3616 BatteryModel<T>::getPotentialResidual()
3617 {
3618  return _impl->getPotentialResidual();
3619 }
3620 
3621 template<class T>
3622 T
3623 BatteryModel<T>::getThermalResidual()
3624 {
3625  return _impl->getThermalResidual();
3626 }
3627 
3628 template<class T>
3629 T
3630 BatteryModel<T>::getPCResidual(const int v)
3631 {
3632  return _impl->getPCResidual(v);
3633 }
3634 
3635 template<class T>
3636 void
3637 BatteryModel<T>::copyCoupledToSeparate()
3638 {
3639  _impl->copyCoupledToSeparate();
3640 }
3641 
3642 template<class T>
3643 void
3644 BatteryModel<T>::copySeparateToCoupled()
3645 {
3646  _impl->copySeparateToCoupled();
3647 }
3648 
3649 template<class T>
3650 T
3651 BatteryModel<T>::getFaceGroupArea(const Mesh& mesh, const int fgID)
3652 {
3653  return _impl->getFaceGroupArea(mesh, fgID);
3654 }
3655 
3656 template<class T>
3657 T
3658 BatteryModel<T>::getFaceGroupVoltage(const Mesh& mesh, const int fgID)
3659 {
3660  return _impl->getFaceGroupVoltage(mesh, fgID);
3661 }
3662 
3663 template<class T>
3664 T
3665 BatteryModel<T>::getMeshVolume(const Mesh& mesh)
3666 {
3667  return _impl->getMeshVolume(mesh);
3668 }
3669 
BatteryModel::Impl::_batteryModelFields
BatteryModelFields _batteryModelFields
Definition: BatteryModel_impl.h:3469
BatteryModel::Impl::_speciesFieldsVector
vector< BatterySpeciesFields * > _speciesFieldsVector
Definition: BatteryModel_impl.h:3451
BatteryModel::Impl::_currentSpeciesResidual
vector< MFRPtr * > _currentSpeciesResidual
Definition: BatteryModel_impl.h:3474
BatteryModel::Impl::getFaceGroupVoltage
T getFaceGroupVoltage(const Mesh &mesh, const int fgID)
Definition: BatteryModel_impl.h:3406
Mesh::getBoundaryFaceGroups
const FaceGroupList & getBoundaryFaceGroups() const
Definition: Mesh.h:187
BatteryModel::Impl::VectorT3GradArray
Array< Gradient< VectorT3 > > VectorT3GradArray
Definition: BatteryModel_impl.h:56
BatteryModel::Impl::TGradType
Gradient< T > TGradType
Definition: BatteryModel_impl.h:48
Mesh::isDoubleShell
bool isDoubleShell() const
Definition: Mesh.h:324
LinearSystem::initAssembly
void initAssembly()
Definition: LinearSystem.cpp:32
BatteryModel::advanceThermal
void advanceThermal(const int niter)
Definition: BatteryModel_impl.h:3552
BatteryModel::getSpeciesBCMap
BatterySpeciesBCMap & getSpeciesBCMap(const int speciesId)
Definition: BatteryModel_impl.h:3506
MultiFieldMatrix::getMatrix
Matrix & getMatrix(const Index &rowIndex, const Index &colIndex)
Definition: MultiFieldMatrix.cpp:42
FaceGroupPtr
shared_ptr< FaceGroup > FaceGroupPtr
Definition: Mesh.h:46
StorageSite::getSelfCount
int getSelfCount() const
Definition: StorageSite.h:40
BatteryModel::Impl::getBatteryModelFields
BatteryModelFields & getBatteryModelFields()
Definition: BatteryModel_impl.h:598
BatteryModel::Impl::getSpeciesBCMap
BatterySpeciesBCMap & getSpeciesBCMap(const int speciesId)
Definition: BatteryModel_impl.h:600
BatteryModel::getPotentialVCMap
BatteryPotentialVCMap & getPotentialVCMap()
Definition: BatteryModel_impl.h:3522
BatteryModel::Impl::_tbcMap
BatteryThermalBCMap _tbcMap
Definition: BatteryModel_impl.h:3457
BatteryModel::Impl::linearizeSpecies
void linearizeSpecies(LinearSystem &ls, LinearSystem &lsShell, const int &m)
Definition: BatteryModel_impl.h:1186
BatteryModel::getOptions
BatteryModelOptions< T > & getOptions()
Definition: BatteryModel_impl.h:3534
Field::hasArray
bool hasArray(const StorageSite &s) const
Definition: Field.cpp:37
Mesh::getParentFaceGroupSite
const StorageSite & getParentFaceGroupSite() const
Definition: Mesh.h:329
BatteryModel::getPotentialBCMap
BatteryPotentialBCMap & getPotentialBCMap()
Definition: BatteryModel_impl.h:3518
BatteryModel::getMassFluxIntegral
T getMassFluxIntegral(const Mesh &mesh, const int faceGroupId, const int m)
Definition: BatteryModel_impl.h:3581
FaceGroup
Definition: Mesh.h:28
BatteryPCLinearizeInterface_BV::discretize
void discretize(const Mesh &mesh, const Mesh &parentMesh, const Mesh &otherMesh, MultiFieldMatrix &mfmatrix, MultiField &xField, MultiField &rField)
Definition: BatteryPCLinearizeInterface_BV.h:56
BatteryModel::Impl::_initialSpeciesNormVector
vector< MFRPtr * > _initialSpeciesNormVector
Definition: BatteryModel_impl.h:3465
BatteryModel::Impl::printMatrixElementsOnFace
void printMatrixElementsOnFace(const Mesh &mesh, const int fgId, LinearSystem &ls, Field &varField)
Definition: BatteryModel_impl.h:2510
Mesh::getAllFaceGroups
const FaceGroupList & getAllFaceGroups() const
Definition: Mesh.h:193
BatteryModel::copySeparateToCoupled
void copySeparateToCoupled()
Definition: BatteryModel_impl.h:3644
BatteryModel::Impl::initPotentialLinearization
void initPotentialLinearization(LinearSystem &ls)
Definition: BatteryModel_impl.h:1076
BatteryModel::Impl::getSpeciesResidual
T getSpeciesResidual(const int speciesId)
Definition: BatteryModel_impl.h:2829
Field
Definition: Field.h:14
BatteryModel::Impl::Impl
Impl(GeomFields &geomFields, const MeshList &realMeshes, const MeshList &meshes, const int nSpecies)
Definition: BatteryModel_impl.h:68
Mesh::isConnectedShell
bool isConnectedShell() const
Definition: Mesh.h:325
CRMatrix::getCoeff
OffDiag & getCoeff(const int i, const int j)
Definition: CRMatrix.h:836
BatteryModel::getAverageConcentration
T getAverageConcentration(const Mesh &mesh, const int m)
Definition: BatteryModel_impl.h:3602
BatteryModel::Impl::getBatterySpeciesFields
BatterySpeciesFields & getBatterySpeciesFields(const int speciesId)
Definition: BatteryModel_impl.h:597
BatteryModel::getSpeciesVCMap
BatterySpeciesVCMap & getSpeciesVCMap(const int speciesId)
Definition: BatteryModel_impl.h:3510
BatteryModel::advanceCoupled
void advanceCoupled(const int niter)
Definition: BatteryModel_impl.h:3559
BatteryModel::Impl::matrixMerge
LinearSystem & matrixMerge(LinearSystem &ls, const int m)
Definition: BatteryModel_impl.h:2563
Mesh::getFaceGroup
const FaceGroup & getFaceGroup(const int fgId) const
Definition: Mesh.cpp:1570
Mesh
Definition: Mesh.h:49
BatteryModel::getPotentialFluxIntegral
T getPotentialFluxIntegral(const Mesh &mesh, const int faceGroupId)
Definition: BatteryModel_impl.h:3588
BatteryModel::~BatteryModel
virtual ~BatteryModel()
Definition: BatteryModel_impl.h:3489
BatteryModel::getBatteryModelFields
BatteryModelFields & getBatteryModelFields()
Definition: BatteryModel_impl.h:3514
BatteryModel::Impl::getThermalVCMap
BatteryThermalVCMap & getThermalVCMap()
Definition: BatteryModel_impl.h:603
CRMatrix::getDiag
Array< Diag > & getDiag()
Definition: CRMatrix.h:856
BatteryModel::Impl::_meshes
const MeshList _meshes
Definition: BatteryModel_impl.h:3448
BatteryModel::Impl::DiagTensorT3
DiagonalTensor< T, 3 > DiagTensorT3
Definition: BatteryModel_impl.h:58
logCtor
#define logCtor()
Definition: RLogInterface.h:26
BatteryLinearizeSpeciesInterface::discretize
void discretize(const Mesh &mesh, const Mesh &parentMesh, const Mesh &otherMesh, MultiFieldMatrix &mfmatrix, MultiField &xField, MultiField &rField)
Definition: BatteryLinearizeSpeciesInterface.h:57
BatteryModel::Impl::TGradArray
Array< Gradient< T > > TGradArray
Definition: BatteryModel_impl.h:49
BatteryModel::advanceSpecies
void advanceSpecies(const int niter)
Definition: BatteryModel_impl.h:3538
BatteryModel::Impl::getPCResidual
T getPCResidual(const int v)
Definition: BatteryModel_impl.h:2858
FaceGroup::groupType
string groupType
Definition: Mesh.h:42
FloatVarDict::getVal
FloatVal< T > getVal(const string varName) const
Definition: FloatVarDict.h:85
VectorT3Array
Array< Vector< double, 3 > > VectorT3Array
Definition: CellMark.cpp:7
Mesh::getOtherMeshID
int getOtherMeshID() const
Definition: Mesh.h:327
BatteryModel::Impl::linearizeThermal
void linearizeThermal(LinearSystem &ls)
Definition: BatteryModel_impl.h:1612
BatteryModel::Impl::VectorT3Array
Array< VectorT3 > VectorT3Array
Definition: BatteryModel_impl.h:54
BatteryModel::getFaceGroupVoltage
T getFaceGroupVoltage(const Mesh &mesh, const int fgID)
Definition: BatteryModel_impl.h:3658
Mesh::getOtherFaceGroupSite
const StorageSite & getOtherFaceGroupSite() const
Definition: Mesh.h:332
BatteryPC_BCS.h
BatteryModel::advancePotential
void advancePotential(const int niter)
Definition: BatteryModel_impl.h:3545
BatteryModel::Impl::SquareTensorT2
SquareTensor< T, 2 > SquareTensorT2
Definition: BatteryModel_impl.h:63
BatteryModel::Impl::getPotentialFluxIntegral
T getPotentialFluxIntegral(const Mesh &mesh, const int faceGroupId)
Definition: BatteryModel_impl.h:2411
BatteryModel::BatterySpeciesVCMap
std::map< int, BatterySpeciesVC< T > * > BatterySpeciesVCMap
Definition: BatteryModel.h:24
FaceGroup::id
const int id
Definition: Mesh.h:41
LinearSystem::updateSolution
void updateSolution()
Definition: LinearSystem.cpp:265
BatteryModel::BatteryModel
BatteryModel(GeomFields &geomFields, const MeshList &realMeshes, const MeshList &meshes, const int nSpecies)
Definition: BatteryModel_impl.h:3478
BatteryModel::Impl::initThermalLinearization
void initThermalLinearization(LinearSystem &ls)
Definition: BatteryModel_impl.h:1131
BatteryModel::Impl::getMeshVolume
T getMeshVolume(const Mesh &mesh)
Definition: BatteryModel_impl.h:3434
BatteryModel::Impl::getAverageConcentration
T getAverageConcentration(const Mesh &mesh, const int m)
Definition: BatteryModel_impl.h:2493
LinearizeInterfaceJumpUnconnected::discretize
void discretize(const Mesh &mesh, const Mesh &parentMesh, const Mesh &otherMesh, MultiFieldMatrix &mfmatrix, MultiField &xField, MultiField &rField, LinearSystem &lsShell)
Definition: LinearizeInterfaceJumpUnconnected.h:46
MultiField::ArrayIndex
pair< const Field *, const StorageSite * > ArrayIndex
Definition: MultiField.h:21
BatteryModel::Impl::_sbcMapVector
vector< BatterySpeciesBCMap * > _sbcMapVector
Definition: BatteryModel_impl.h:3453
BatteryModel::Impl::advancePotential
void advancePotential(const int niter)
Definition: BatteryModel_impl.h:2710
BatteryPCLinearizeInterface_BV.h
BatteryModel::Impl::VectorT2Grad
Gradient< VectorT2 > VectorT2Grad
Definition: BatteryModel_impl.h:65
BatteryModel::Impl::linearizePC_Thermal
void linearizePC_Thermal(LinearSystem &ls)
Definition: BatteryModel_impl.h:1814
LinearSystem::initSolve
void initSolve()
Definition: LinearSystem.cpp:40
BatteryModel::getSpeciesResidual
T getSpeciesResidual(const int speciesId)
Definition: BatteryModel_impl.h:3609
BatteryModel::Impl::getSpeciesVCMap
BatterySpeciesVCMap & getSpeciesVCMap(const int speciesId)
Definition: BatteryModel_impl.h:599
BatteryModel::recoverLastTimestep
void recoverLastTimestep()
Definition: BatteryModel_impl.h:3573
BatteryPCLinearizeInterface_BV
Definition: BatteryPCLinearizeInterface_BV.h:24
DiscrList
vector< shared_ptr< Discretization > > DiscrList
Definition: Discretization.h:31
Mesh::getFaces
const StorageSite & getFaces() const
Definition: Mesh.h:108
BatteryModel::BatterySpeciesBCMap
std::map< int, BatterySpeciesBC< T > * > BatterySpeciesBCMap
Definition: BatteryModel.h:23
BatteryModel::Impl::_realMeshes
const MeshList _realMeshes
Definition: BatteryModel_impl.h:3447
Model::_meshes
const MeshList _meshes
Definition: Model.h:29
Mesh::getCellCells
const CRConnectivity & getCellCells() const
Definition: Mesh.cpp:480
Model
Definition: Model.h:13
Mesh::getCells
const StorageSite & getCells() const
Definition: Mesh.h:109
BaseGenericBCS::applyInterfaceBC
void applyInterfaceBC(const int f) const
Definition: GenericBCS.h:325
BatteryModel::Impl::advanceSpecies
void advanceSpecies(const int niter)
Definition: BatteryModel_impl.h:2638
BatteryModel::Impl::VectorT3Grad
Gradient< VectorT3 > VectorT3Grad
Definition: BatteryModel_impl.h:57
BatteryModel::getThermalVCMap
BatteryThermalVCMap & getThermalVCMap()
Definition: BatteryModel_impl.h:3530
logDtor
#define logDtor()
Definition: RLogInterface.h:33
BatteryModel::Impl::VectorT2
Vector< T, 2 > VectorT2
Definition: BatteryModel_impl.h:61
BatteryModel::Impl::getHeatFluxIntegral
T getHeatFluxIntegral(const Mesh &mesh, const int faceGroupId)
Definition: BatteryModel_impl.h:2452
Field::addArray
void addArray(const StorageSite &, shared_ptr< ArrayBase > a)
Definition: Field.cpp:72
Mesh::getFaceCells
const CRConnectivity & getFaceCells(const StorageSite &site) const
Definition: Mesh.cpp:388
BatteryModel::Impl::VectorT2GradArray
Array< Gradient< VectorT2 > > VectorT2GradArray
Definition: BatteryModel_impl.h:64
BatteryPC_BCS
Definition: BatteryPC_BCS.h:22
Array
Definition: Array.h:14
Mesh::getParentMeshID
int getParentMeshID() const
Definition: Mesh.h:326
BatteryModel::Impl::SquareTensorT3
SquareTensor< T, 3 > SquareTensorT3
Definition: BatteryModel_impl.h:55
BatteryModel::BatteryPotentialBCMap
std::map< int, BatteryPotentialBC< T > * > BatteryPotentialBCMap
Definition: BatteryModel.h:25
BatteryModel::Impl::getThermalBCMap
BatteryThermalBCMap & getThermalBCMap()
Definition: BatteryModel_impl.h:604
MultiFieldMatrix::addMatrix
void addMatrix(const Index &rowI, const Index &colI, shared_ptr< Matrix > m)
Definition: MultiFieldMatrix.h:137
BatteryModel::Impl::getPotentialVCMap
BatteryPotentialVCMap & getPotentialVCMap()
Definition: BatteryModel_impl.h:601
BatteryModel::Impl::DiagTensorT2
DiagonalTensor< T, 2 > DiagTensorT2
Definition: BatteryModel_impl.h:66
BatteryModel::getPCResidual
T getPCResidual(const int v)
Definition: BatteryModel_impl.h:3630
BatteryModel::Impl::VectorT3
Vector< T, 3 > VectorT3
Definition: BatteryModel_impl.h:53
ArrayBase::getData
virtual void * getData() const =0
BatteryModel::Impl::initPCLinearization
void initPCLinearization(LinearSystem &ls)
Definition: BatteryModel_impl.h:883
BatteryModel::Impl::advanceThermal
void advanceThermal(const int niter)
Definition: BatteryModel_impl.h:2748
StorageSite::getCount
int getCount() const
Definition: StorageSite.h:39
LinearSystem::getB
MultiField & getB()
Definition: LinearSystem.h:33
BatteryModel::Impl::T_Matrix
CRMatrix< T, T, T > T_Matrix
Definition: BatteryModel_impl.h:50
BatteryModel::Impl::linearizePC
void linearizePC(LinearSystem &ls)
Definition: BatteryModel_impl.h:2134
BatteryModel::Impl::_svcMapVector
vector< BatterySpeciesVCMap * > _svcMapVector
Definition: BatteryModel_impl.h:3454
BatteryModel::Impl::_pvcMap
BatteryPotentialVCMap _pvcMap
Definition: BatteryModel_impl.h:3456
BatteryModel::getFaceGroupArea
T getFaceGroupArea(const Mesh &mesh, const int fgID)
Definition: BatteryModel_impl.h:3651
BatteryModel::copyCoupledToSeparate
void copyCoupledToSeparate()
Definition: BatteryModel_impl.h:3637
MFRPtr
shared_ptr< MultiFieldReduction > MFRPtr
Definition: MultiFieldReduction.h:56
MultiField::addArray
void addArray(const ArrayIndex &aIndex, shared_ptr< ArrayBase > a)
Definition: MultiField.cpp:270
BatteryModel::getBatterySpeciesFields
BatterySpeciesFields & getBatterySpeciesFields(const int speciesId)
Definition: BatteryModel_impl.h:3502
BatteryModel::getMeshVolume
T getMeshVolume(const Mesh &mesh)
Definition: BatteryModel_impl.h:3665
Mesh::getInterfaceGroups
const FaceGroupList & getInterfaceGroups() const
Definition: Mesh.h:190
LinearSystem::getX
MultiField & getX()
Definition: LinearSystem.h:32
BatteryModel::Impl::getFaceGroupArea
T getFaceGroupArea(const Mesh &mesh, const int fgID)
Definition: BatteryModel_impl.h:3383
BatteryLinearizePotentialInterface::discretize
void discretize(const Mesh &mesh, const Mesh &parentMesh, const Mesh &otherMesh, MultiFieldMatrix &mfmatrix, MultiField &xField, MultiField &rField)
Definition: BatteryLinearizePotentialInterface.h:53
BatteryModel::init
virtual void init()
Definition: BatteryModel_impl.h:3496
BatteryModel::Impl::VectorT2Array
Array< VectorT2 > VectorT2Array
Definition: BatteryModel_impl.h:62
BatteryLinearizeThermalInterface::discretize
void discretize(const Mesh &mesh, const Mesh &parentMesh, const Mesh &otherMesh, MultiFieldMatrix &mfmatrix, MultiField &xField, MultiField &rField)
Definition: BatteryLinearizeThermalInterface.h:53
BatteryModel::Impl::getMassFluxIntegral
T getMassFluxIntegral(const Mesh &mesh, const int faceGroupId, const int m)
Definition: BatteryModel_impl.h:2369
BatteryModel::Impl::linearizePotential
void linearizePotential(LinearSystem &ls)
Definition: BatteryModel_impl.h:1429
BatteryModel::BatteryPotentialVCMap
std::map< int, BatteryPotentialVC< T > * > BatteryPotentialVCMap
Definition: BatteryModel.h:26
BatteryModel::BatteryThermalBCMap
std::map< int, BatteryThermalBC< T > * > BatteryThermalBCMap
Definition: BatteryModel.h:27
Field::getArrayPtr
shared_ptr< ArrayBase > getArrayPtr(const StorageSite &)
Definition: Field.cpp:63
BatteryModel::Impl::getPotentialBCMap
BatteryPotentialBCMap & getPotentialBCMap()
Definition: BatteryModel_impl.h:602
BatteryModel::Impl::initSpeciesLinearization
void initSpeciesLinearization(LinearSystem &ls, LinearSystem &lsShell, const int &SpeciesNumber)
Definition: BatteryModel_impl.h:778
BatteryModel::Impl::advanceCoupled
void advanceCoupled(const int niter)
Definition: BatteryModel_impl.h:2786
BatteryModel::getHeatFluxIntegral
T getHeatFluxIntegral(const Mesh &mesh, const int faceGroupId)
Definition: BatteryModel_impl.h:3595
Field::syncLocal
void syncLocal()
Definition: Field.cpp:334
Mesh::getID
int getID() const
Definition: Mesh.h:106
BatteryModel::Impl::_tvcMap
BatteryThermalVCMap _tvcMap
Definition: BatteryModel_impl.h:3458
Linearizer::linearize
virtual void linearize(DiscrList &discretizations, const MeshList &meshes, MultiFieldMatrix &matrix, MultiField &x, MultiField &b)
Definition: Linearizer.cpp:17
LinearizeInterfaceJump::discretize
void discretize(const Mesh &mesh, const Mesh &parentMesh, const Mesh &otherMesh, MultiFieldMatrix &mfmatrix, MultiField &xField, MultiField &rField)
Definition: LinearizeInterfaceJump.h:46
BatteryModel::BatteryThermalVCMap
std::map< int, BatteryThermalVC< T > * > BatteryThermalVCMap
Definition: BatteryModel.h:28
BatteryModel::getThermalBCMap
BatteryThermalBCMap & getThermalBCMap()
Definition: BatteryModel_impl.h:3526
LinearSystem::getMatrix
MultiFieldMatrix & getMatrix()
Definition: LinearSystem.h:37
MeshList
vector< Mesh * > MeshList
Definition: Mesh.h:439
BatteryModel::Impl::getOptions
BatteryModelOptions< T > & getOptions()
Definition: BatteryModel_impl.h:606
BatteryModel::Impl::_pbcMap
BatteryPotentialBCMap _pbcMap
Definition: BatteryModel_impl.h:3455
FaceGroup::site
StorageSite site
Definition: Mesh.h:40
BatteryModel::Impl::_options
BatteryModelOptions< T > _options
Definition: BatteryModel_impl.h:3460