#include <SpeciesModel_impl.h>

Collaboration diagram for SpeciesModel< T >::Impl:

Public Types
typedef Array< T >	TArray

typedef Vector< T, 3 >	VectorT3

typedef Array< VectorT3 >	VectorT3Array

typedef Gradient< T >	TGradType

typedef Array< Gradient< T > >	TGradArray

typedef CRMatrix< T, T, T >	T_Matrix

Public Member Functions
	Impl (const GeomFields &geomFields, const MeshList &meshes, const int nSpecies)

void	init ()

SpeciesFields &	getSpeciesFields (const int speciesId)

SpeciesVCMap &	getVCMap (const int speciesId)

SpeciesBCMap &	getBCMap (const int speciesId)

SpeciesBC< T > &	getBC (const int id, const int speciesId)

SpeciesModelOptions< T > &	getOptions ()

void	updateTime ()

void	initLinearization (LinearSystem &ls, const int &m)

void	linearize (LinearSystem &ls, const int &m)

T	getMassFluxIntegral (const Mesh &mesh, const int faceGroupId, const int m)

T	getAverageMassFraction (const Mesh &mesh, const int m)

void	advance (const int niter)

void	printBCs ()

T	getMassFractionResidual (const int speciesId)

Private Attributes
const MeshList	_meshes

const GeomFields &	_geomFields

vector< SpeciesFields * >	_speciesFieldsVector

vector< SpeciesBCMap * >	_bcMapVector

vector< SpeciesVCMap * >	_vcMapVector

SpeciesModelOptions< T >	_options

vector< MFRPtr * >	_initialNormVector

int	_niters

const int	_nSpecies

SpeciesModelFields	_speciesModelFields

vector< MFRPtr * >	_currentResidual

Detailed Description

template<class T>
class SpeciesModel< T >::Impl

Definition at line 33 of file SpeciesModel_impl.h.

Member Typedef Documentation

template<class T>

typedef CRMatrix<T,T,T> SpeciesModel< T >::Impl::T_Matrix

Definition at line 41 of file SpeciesModel_impl.h.

template<class T>

typedef Array<T> SpeciesModel< T >::Impl::TArray

Definition at line 36 of file SpeciesModel_impl.h.

template<class T>

typedef Array<Gradient<T> > SpeciesModel< T >::Impl::TGradArray

Definition at line 40 of file SpeciesModel_impl.h.

template<class T>

typedef Gradient<T> SpeciesModel< T >::Impl::TGradType

Definition at line 39 of file SpeciesModel_impl.h.

template<class T>

typedef Vector<T,3> SpeciesModel< T >::Impl::VectorT3

Definition at line 37 of file SpeciesModel_impl.h.

template<class T>

typedef Array<VectorT3> SpeciesModel< T >::Impl::VectorT3Array

Definition at line 38 of file SpeciesModel_impl.h.

Constructor & Destructor Documentation

template<class T>

SpeciesModel< T >::Impl::Impl	(	const GeomFields &	geomFields,
		const MeshList &	meshes,
		const int	nSpecies
	)

inline

Definition at line 43 of file SpeciesModel_impl.h.

References Model::_meshes, SpeciesBC< T >::bcType, Mesh::getBoundaryFaceGroups(), Mesh::getID(), FaceGroup::groupType, FaceGroup::id, and SpeciesVC< T >::vcType.

                            :
     _meshes(meshes),
     _geomFields(geomFields),
     _niters(0),
     _nSpecies(nSpecies),
     _speciesModelFields("speciesModel")
   {
     
     const int numMeshes = _meshes.size();
 
     for (int m=0; m<_nSpecies; m++)
     {
       SpeciesVCMap *svcmap = new SpeciesVCMap();
       _vcMapVector.push_back(svcmap);
       SpeciesBCMap *sbcmap = new SpeciesBCMap();
       _bcMapVector.push_back(sbcmap);
       SpeciesFields *sFields = new SpeciesFields("species");
       _speciesFieldsVector.push_back(sFields);
       MFRPtr *iNorm = new MFRPtr();
       _initialNormVector.push_back(iNorm);
       MFRPtr *rCurrent = new MFRPtr();
       _currentResidual.push_back(rCurrent); 
         
       for (int n=0; n<numMeshes; n++)
       {
 
         const Mesh& mesh = *_meshes[n];
         SpeciesVC<T> *vc(new SpeciesVC<T>());
         vc->vcType = "flow";
         (*svcmap)[mesh.getID()] = vc;
         
         foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
         {
             const FaceGroup& fg = *fgPtr;
             SpeciesBC<T> *bc(new SpeciesBC<T>());
             
             (*sbcmap)[fg.id] = bc;
 
             if ((fg.groupType == "wall") ||
                 (fg.groupType == "symmetry"))
             {
                 bc->bcType = "SpecifiedMassFlux";
             }
             else if ((fg.groupType == "velocity-inlet") ||
                      (fg.groupType == "pressure-outlet"))
             {
                 bc->bcType = "SpecifiedMassFraction";
             }
             else
               throw CException("SpeciesModel: unknown face group type "
                                + fg.groupType);
         }
       }
     }
   }

Member Function Documentation

template<class T>

void SpeciesModel< T >::Impl::advance ( const int niter )

inline

Definition at line 682 of file SpeciesModel_impl.h.

References LinearSystem::initAssembly(), LinearSystem::initSolve(), LinearSystem::postSolve(), and LinearSystem::updateSolution().

   {
     for(int n=0; n<niter; n++)
     { 
       bool allConverged=true;
       for (int m=0; m<_nSpecies; m++)
       {
         MFRPtr& iNorm = *_initialNormVector[m];
         MFRPtr& rCurrent = *_currentResidual[m];
 
         LinearSystem ls;
         initLinearization(ls, m);
         
         ls.initAssembly();
 
         linearize(ls, m);
 
         ls.initSolve();
 
         MFRPtr rNorm(_options.getLinearSolver().solve(ls));
 
         if (!iNorm) iNorm = rNorm;
         
         MFRPtr normRatio((*rNorm)/(*iNorm));
 
         cout << "Species Number: " << m << endl;
         cout << _niters << ": " << *rNorm << endl;
         
         _options.getLinearSolver().cleanup();
 
         ls.postSolve();
         ls.updateSolution();
 
         _niters++;
 
         rCurrent = rNorm;
  
         if (!(*rNorm < _options.absoluteTolerance ||
             *normRatio < _options.relativeTolerance))
             allConverged=false;
       }
       if (allConverged)
         break;
     }
   }

template<class T>

T SpeciesModel< T >::Impl::getAverageMassFraction	(	const Mesh &	mesh,
		const int	m
	)

inline

Definition at line 665 of file SpeciesModel_impl.h.

References Mesh::getCells(), StorageSite::getSelfCount(), and SpeciesFields::massFraction.

   {
     SpeciesFields& sFields = *_speciesFieldsVector[m];
     const StorageSite& cells = mesh.getCells();
     const TArray& mFCell = dynamic_cast<const TArray&>(sFields.massFraction[cells]);
     const TArray& cellVolume = dynamic_cast<const TArray&>(_geomFields.volume[cells]);
     T totalVolume = 0.0;
     T weightedMassFraction = 0.0;
     for(int c=0; c<cells.getSelfCount(); c++)
       {
         totalVolume += cellVolume[c];
         weightedMassFraction += mFCell[c]*cellVolume[c];
       }
     return weightedMassFraction/totalVolume;
   }

template<class T>

SpeciesBC<T>& SpeciesModel< T >::Impl::getBC	(	const int	id,
		const int	speciesId
	)

inline

Definition at line 337 of file SpeciesModel_impl.h.

337 {return *(*_bcMapVector[speciesId])[id];}

SpeciesModel::Impl::_bcMapVector

vector< SpeciesBCMap * > _bcMapVector

Definition: SpeciesModel_impl.h:767

template<class T>

SpeciesBCMap& SpeciesModel< T >::Impl::getBCMap ( const int speciesId )

inline

Definition at line 335 of file SpeciesModel_impl.h.

335 {return *_bcMapVector[speciesId];}

SpeciesModel::Impl::_bcMapVector

vector< SpeciesBCMap * > _bcMapVector

Definition: SpeciesModel_impl.h:767

template<class T>

T SpeciesModel< T >::Impl::getMassFluxIntegral	(	const Mesh &	mesh,
		const int	faceGroupId,
		const int	m
	)

inline

Definition at line 623 of file SpeciesModel_impl.h.

References Mesh::getBoundaryFaceGroups(), StorageSite::getCount(), Mesh::getInterfaceGroups(), FaceGroup::id, SpeciesFields::massFlux, and FaceGroup::site.

   {
     SpeciesFields& sFields = *_speciesFieldsVector[m];
 
     T r(0.);
     bool found = false;
     foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
     {
         const FaceGroup& fg = *fgPtr;
         if (fg.id == faceGroupId)
         {
             const StorageSite& faces = fg.site;
             const int nFaces = faces.getCount();
             const TArray& massFlux =
               dynamic_cast<const TArray&>(sFields.massFlux[faces]);
             for(int f=0; f<nFaces; f++)
               r += massFlux[f];
             found=true;
         }
     }
 
     foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
     {
         const FaceGroup& fg = *fgPtr;
         if (fg.id == faceGroupId)
           {
             const StorageSite& faces = fg.site;
             const int nFaces = faces.getCount();
             const TArray& massFlux =
               dynamic_cast<const TArray&>(sFields.massFlux[faces]);
             for(int f=0; f<nFaces; f++)
               r += massFlux[f];
             found=true;
           }
     }
 
 
     if (!found)
       throw CException("getMassFluxIntegral: invalid faceGroupID");
     return r;
   }

template<class T>

T SpeciesModel< T >::Impl::getMassFractionResidual ( const int speciesId )

inline

Definition at line 749 of file SpeciesModel_impl.h.

References ArrayBase::getData(), and SpeciesFields::massFraction.

   {
     MFRPtr& rCurrent = *_currentResidual[speciesId];
     SpeciesFields& sFields = *_speciesFieldsVector[speciesId];
     const Field& massFractionField = sFields.massFraction;
     ArrayBase& residualArray = (*rCurrent)[massFractionField];
     T *arrayData = (T*)(residualArray.getData());
     const T residual = arrayData[0];
     return residual;
   }

template<class T>

SpeciesModelOptions<T>& SpeciesModel< T >::Impl::getOptions ( )

inline

Definition at line 339 of file SpeciesModel_impl.h.

339 {return _options;}

SpeciesModel::Impl::_options

SpeciesModelOptions< T > _options

Definition: SpeciesModel_impl.h:770

template<class T>

SpeciesFields& SpeciesModel< T >::Impl::getSpeciesFields ( const int speciesId )

inline

Definition at line 332 of file SpeciesModel_impl.h.

332 {return *_speciesFieldsVector[speciesId];}

SpeciesModel::Impl::_speciesFieldsVector

vector< SpeciesFields * > _speciesFieldsVector

Definition: SpeciesModel_impl.h:765

template<class T>

SpeciesVCMap& SpeciesModel< T >::Impl::getVCMap ( const int speciesId )

inline

Definition at line 334 of file SpeciesModel_impl.h.

334 {return *_vcMapVector[speciesId];}

SpeciesModel::Impl::_vcMapVector

vector< SpeciesVCMap * > _vcMapVector

Definition: SpeciesModel_impl.h:768

template<class T>

void SpeciesModel< T >::Impl::init ( )

inline

Definition at line 101 of file SpeciesModel_impl.h.

References Model::_meshes, Field::addArray(), SpeciesFields::convectionFlux, SpeciesFields::diffusivity, SpeciesFields::elecPotential, Mesh::getBoundaryFaceGroups(), Mesh::getCells(), StorageSite::getCount(), Mesh::getFaces(), Mesh::getID(), Mesh::getInterfaceGroups(), Mesh::getOtherMeshID(), Mesh::getParentMeshID(), Mesh::isDoubleShell(), SpeciesFields::massFlux, SpeciesFields::massFraction, SpeciesFields::massFractionElectricModel, SpeciesFields::massFractionN1, SpeciesFields::massFractionN2, SpeciesFields::one, FaceGroup::site, SpeciesFields::source, Field::syncLocal(), and SpeciesFields::zero.

   {
     const int numMeshes = _meshes.size();
 
     for (int m=0; m<_nSpecies; m++)
     {
       const SpeciesVCMap& svcmap = *_vcMapVector[m];
       const SpeciesBCMap& sbcmap = *_bcMapVector[m];
       SpeciesFields& sFields = *_speciesFieldsVector[m];
       //MFRPtr& iNorm = *_initialNormVector[m];
 
       for (int n=0; n<numMeshes; n++)
       {
         const Mesh& mesh = *_meshes[n];
 
         const StorageSite& cells = mesh.getCells();
         const StorageSite& allFaces = mesh.getFaces();
  
         typename SpeciesVCMap::const_iterator pos = svcmap.find(mesh.getID());
         if (pos==svcmap.end())
         {   
            throw CException("SpeciesModel: Error in VC Map");
         }
 
         const SpeciesVC<T>& vc = *(pos->second);
      
 
         //mass fraction
         
         shared_ptr<TArray> mFCell(new TArray(cells.getCount()));
         
         if (!mesh.isDoubleShell())
           {
             *mFCell = vc["initialMassFraction"];
           }
         else
           {
             // double shell mesh cells take on values of parents
             typename SpeciesVCMap::const_iterator posParent = svcmap.find(mesh.getParentMeshID());
             typename SpeciesVCMap::const_iterator posOther = svcmap.find(mesh.getOtherMeshID());
             const SpeciesVC<T>& vcP = *(posParent->second);
             const SpeciesVC<T>& vcO = *(posOther->second);
 
             const int NumberOfCells = cells.getCount();
             for (int c=0; c<NumberOfCells; c++)
               {
                 if ((c < NumberOfCells/4)||(c >= 3*NumberOfCells/4))
                   {
                     (*mFCell)[c] = vcP["initialMassFraction"];
                   }
                 else
                   {
                     (*mFCell)[c] = vcO["initialMassFraction"];
                   }
               }
           }
         /*
         //Initialize to exp for ln term testing in electric model
         const VectorT3Array& cellCentroid =
           dynamic_cast<const VectorT3Array&>(_geomFields.coordinate[cells]);
 
         for (int c=0; c<cells.getCount(); c++)
         {
           VectorT3 CellCent = cellCentroid[c]; 
           (*mFCell)[c] = exp(pow(CellCent[0],2));
           }*/
 
 
         /*if (n==0)
           {
             *mFCell =_options["initialMassFraction0"];
           }
         else if (n==1)
           {
             *mFCell = _options["initialMassFraction1"];
           }
         else if (n==2)
           {
             *mFCell = _options["initialMassFraction2"];
           }
         else if (n==3)
           {
             *mFCell = _options["initialMassFraction3"];
           }
         else if (n==4)
           {
             *mFCell = _options["initialMassFraction4"];
           }
         else
           {
             *mFCell =_options["initialMassFraction0"];
             }
         
         if (mesh.isDoubleShell())
           {
             char parentInitial[21];
             char otherInitial[21];
             sprintf(parentInitial, "initialMassFraction%d", mesh.getParentMeshID());
             sprintf(otherInitial, "initialMassFraction%d", mesh.getOtherMeshID());
 
             for (int c=0; c<cells.getSelfCount(); c++)
               {
                 if (c < cells.getSelfCount()/2)
                   {
                     (*mFCell)[c] = _options[parentInitial];
                   }
                 else
                   {
                     (*mFCell)[c] = _options[otherInitial];
                   }
               }
           }*/
 
 
         //*mFCell = _options["initialMassFraction"];
 
         /* //Initialize to 1-D linear solution in x direction
         const VectorT3Array& cellCentroid =
           dynamic_cast<const VectorT3Array&>(_geomFields.coordinate[cells]);
 
         for (int c=0; c<cells.getCount(); c++)
         {
           VectorT3 CellCent = cellCentroid[c]; 
           (*mFCell)[c] = -0.05*CellCent[0] + 0.5;
           } */
 
         /*//Initialize to specialized concentrations (for now)
         const VectorT3Array& cellCentroid =
           dynamic_cast<const VectorT3Array&>(_geomFields.coordinate[cells]);
 
         for (int c=0; c<cells.getCount(); c++)
           {
             VectorT3 CellCent = cellCentroid[c]; 
             if (n==0)
               {
                 (*mFCell)[c] = 1000.0;
               }
             elseif (n==1)
               {
                 (*mFCell)[c] = 0.0;
 
                 (*mFCell)[c] = -0.05*CellCent[0] + 0.5;
                 }*/
         
         sFields.massFraction.addArray(cells,mFCell);
 
         if (_options.ButlerVolmer)
           {
             sFields.massFractionElectricModel.addArray(cells,dynamic_pointer_cast<ArrayBase>(mFCell->newCopy()));
           }
 
         if (_options.transient)
         {
             sFields.massFractionN1.addArray(cells,
                                             dynamic_pointer_cast<ArrayBase>(mFCell->newCopy()));
             if (_options.timeDiscretizationOrder > 1)
               sFields.massFractionN2.addArray(cells,
                                               dynamic_pointer_cast<ArrayBase>(mFCell->newCopy()));
 
         }
 
         //diffusivity
         shared_ptr<TArray> diffusivityCell(new TArray(cells.getCount()));
         *diffusivityCell = vc["massDiffusivity"];
         sFields.diffusivity.addArray(cells,diffusivityCell);
         
         //source 
         shared_ptr<TArray> sCell(new TArray(cells.getCount()));
         *sCell = T(0.);
         sFields.source.addArray(cells,sCell);
 
         //create a zero field
         shared_ptr<TArray> zeroCell(new TArray(cells.getCount()));
         *zeroCell = T(0.0);
         sFields.zero.addArray(cells,zeroCell);
 
         //create a one field
         shared_ptr<TArray> oneCell(new TArray(cells.getCount()));
         *oneCell = T(1.0);
         sFields.one.addArray(cells,oneCell);
 
         //initial temparature gradient array
         shared_ptr<TGradArray> gradT(new TGradArray(cells.getCount()));
         gradT->zero();
         _speciesModelFields.speciesGradient.addArray(cells,gradT);
         
         //inital convection flux at faces
         shared_ptr<TArray> convFlux(new TArray(allFaces.getCount()));
         convFlux->zero();
         sFields.convectionFlux.addArray(allFaces,convFlux);
 
         //mass flux at faces
         foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
         {
             const FaceGroup& fg = *fgPtr;
             const StorageSite& faces = fg.site;
           
             shared_ptr<TArray> fluxFace(new TArray(faces.getCount()));
 
             fluxFace->zero();
             sFields.massFlux.addArray(faces,fluxFace);
           
         }
         foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
         {
             const FaceGroup& fg = *fgPtr;
             const StorageSite& faces = fg.site;
           
             shared_ptr<TArray> fluxFace(new TArray(faces.getCount()));
 
             fluxFace->zero();
             sFields.massFlux.addArray(faces,fluxFace);
           
         }
 
         //electric potential (for coupling to ElectricModel)
         // only needed once for each mesh (not depenedent on species)
         if (m==0)
           {
             shared_ptr<TArray> ePotCell(new TArray(cells.getCount()));
             ePotCell->zero();
             sFields.elecPotential.addArray(cells,ePotCell);
           }
       } 
     
     sFields.diffusivity.syncLocal();
     //iNorm = MFRPtr();
     }
     _niters  =0;    
   }

template<class T>

void SpeciesModel< T >::Impl::initLinearization	(	LinearSystem &	ls,
		const int &	m
	)

inline

Definition at line 370 of file SpeciesModel_impl.h.

References Model::_meshes, MultiField::addArray(), MultiFieldMatrix::addMatrix(), Field::getArrayPtr(), Mesh::getBoundaryFaceGroups(), Mesh::getCellCells(), Mesh::getCells(), StorageSite::getCount(), Mesh::getFaceCells(), Mesh::getInterfaceGroups(), LinearSystem::getMatrix(), LinearSystem::getX(), SpeciesFields::massFlux, SpeciesFields::massFraction, and FaceGroup::site.

   {
     SpeciesFields& sFields = *_speciesFieldsVector[m];
     const int numMeshes = _meshes.size();
 
     for (int n=0; n<numMeshes; n++)
     {
         const Mesh& mesh = *_meshes[n];
 
         const StorageSite& cells = mesh.getCells();
         MultiField::ArrayIndex tIndex(&sFields.massFraction,&cells);
 
         ls.getX().addArray(tIndex,sFields.massFraction.getArrayPtr(cells));
 
         const CRConnectivity& cellCells = mesh.getCellCells();
 
         shared_ptr<Matrix> m(new CRMatrix<T,T,T>(cellCells));
 
         ls.getMatrix().addMatrix(tIndex,tIndex,m);
 
         foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
         {
             const FaceGroup& fg = *fgPtr;
             const StorageSite& faces = fg.site;
 
             MultiField::ArrayIndex fIndex(&sFields.massFlux,&faces);
             ls.getX().addArray(fIndex,sFields.massFlux.getArrayPtr(faces));
 
             const CRConnectivity& faceCells = mesh.getFaceCells(faces);
 
             shared_ptr<Matrix> mft(new FluxJacobianMatrix<T,T>(faceCells));
             ls.getMatrix().addMatrix(fIndex,tIndex,mft);
 
             shared_ptr<Matrix> mff(new DiagonalMatrix<T,T>(faces.getCount()));
             ls.getMatrix().addMatrix(fIndex,fIndex,mff);
         }
 
         foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
         {
             const FaceGroup& fg = *fgPtr;
             const StorageSite& faces = fg.site;
 
             MultiField::ArrayIndex fIndex(&sFields.massFlux,&faces);
             ls.getX().addArray(fIndex,sFields.massFlux.getArrayPtr(faces));
 
             const CRConnectivity& faceCells = mesh.getFaceCells(faces);
 
             shared_ptr<Matrix> mft(new FluxJacobianMatrix<T,T>(faceCells));
             ls.getMatrix().addMatrix(fIndex,tIndex,mft);
 
             shared_ptr<Matrix> mff(new DiagonalMatrix<T,T>(faces.getCount()));
             ls.getMatrix().addMatrix(fIndex,fIndex,mff);
         }
      
     }
   }

template<class T>

void SpeciesModel< T >::Impl::linearize	(	LinearSystem &	ls,
		const int &	m
	)

inline

Definition at line 427 of file SpeciesModel_impl.h.

References Model::_meshes, BaseGenericBCS< X, Diag, OffDiag >::applyInterfaceBC(), SpeciesBC< T >::bcType, GradientModel< X >::compute(), SpeciesFields::convectionFlux, SpeciesFields::diffusivity, LinearizeInterfaceJump< X, Diag, OffDiag >::discretize(), LinearizeSpeciesInterface< X, Diag, OffDiag >::discretize(), SpeciesFields::elecPotential, LinearSystem::getB(), Mesh::getBoundaryFaceGroups(), StorageSite::getCount(), Mesh::getInterfaceGroups(), LinearSystem::getMatrix(), Mesh::getOtherMeshID(), Mesh::getParentMeshID(), FloatVarDict< T >::getVal(), LinearSystem::getX(), Field::hasArray(), FaceGroup::id, Mesh::isDoubleShell(), Linearizer::linearize(), SpeciesFields::massFlux, SpeciesFields::massFraction, SpeciesFields::massFractionElectricModel, SpeciesFields::massFractionN1, SpeciesFields::massFractionN2, SpeciesFields::one, FaceGroup::site, SpeciesFields::source, and SpeciesFields::zero.

   {
     const SpeciesVCMap& svcmap = *_vcMapVector[m];
     const SpeciesBCMap& sbcmap = *_bcMapVector[m];
     SpeciesFields& sFields = *_speciesFieldsVector[m];
 
     GradientModel<T> speciesGradientModel(_meshes,sFields.massFraction,
                                           _speciesModelFields.speciesGradient,_geomFields);
     speciesGradientModel.compute();
     
     DiscrList discretizations;
     
     shared_ptr<Discretization>
       dd(new DiffusionDiscretization<T,T,T>
          (_meshes,_geomFields,
           sFields.massFraction,
           sFields.diffusivity,
           _speciesModelFields.speciesGradient));
     discretizations.push_back(dd);
     
     shared_ptr<Discretization>
       cd(new ConvectionDiscretization<T,T,T>
          (_meshes,_geomFields,
           sFields.massFraction,
           sFields.convectionFlux,
           sFields.zero,
           _speciesModelFields.speciesGradient,
           _options.useCentralDifference));
     discretizations.push_back(cd);
     
     
     shared_ptr<Discretization>
       sd(new SourceDiscretization<T>
          (_meshes, 
           _geomFields, 
           sFields.massFraction,
           sFields.source));
     discretizations.push_back(sd);
     
     if (_options.transient)
     {
         shared_ptr<Discretization>
           td(new TimeDerivativeDiscretization<T,T,T>
              (_meshes,_geomFields,
               sFields.massFraction,
               sFields.massFractionN1,
               sFields.massFractionN2,
               sFields.one,
               _options["timeStep"]));
         
         discretizations.push_back(td);
     }
 
     shared_ptr<Discretization>
       ibm(new GenericIBDiscretization<T,T,T>
           (_meshes,_geomFields,sFields.massFraction));
       
     discretizations.push_back(ibm);
 
     Linearizer linearizer;
 
     linearizer.linearize(discretizations,_meshes,ls.getMatrix(),
                          ls.getX(), ls.getB());
 
     const int numMeshes = _meshes.size();
 
     /* linearize shell mesh */
 
     for (int n=0; n<numMeshes; n++)
       {
         const Mesh& mesh = *_meshes[n];
         if (mesh.isDoubleShell())
           {
             const int parentMeshID = mesh.getParentMeshID();
             const int otherMeshID = mesh.getOtherMeshID();
             const Mesh& parentMesh = *_meshes[parentMeshID];
             const Mesh& otherMesh = *_meshes[otherMeshID];
           
             if (_options.ButlerVolmer)
               {
                 bool Cathode = false;
                 bool Anode = false;
                 if (n == _options["ButlerVolmerCathodeShellMeshID"])
                   {
                     Cathode = true;
                   }
                 else if (n == _options["ButlerVolmerAnodeShellMeshID"])
                   {
                     Anode = true;
                   }
                 
                 LinearizeSpeciesInterface<T, T, T> lbv (_geomFields,
                                                         _options["ButlerVolmerRRConstant"],
                                                         _options["A_coeff"],
                                                         _options["B_coeff"],
                                                         _options["interfaceUnderRelax"],
                                                         Anode,
                                                         Cathode,
                                                         sFields.massFraction,
                                                         sFields.massFractionElectricModel,
                                                         sFields.elecPotential);
 
                 lbv.discretize(mesh, parentMesh, otherMesh, ls.getMatrix(), ls.getX(), ls.getB() );
 
               }
             else
               {
                 LinearizeInterfaceJump<T, T, T> lsm (_options["A_coeff"],
                                                      _options["B_coeff"],
                                                      sFields.massFraction);
 
                 lsm.discretize(mesh, parentMesh, otherMesh, ls.getMatrix(), ls.getX(), ls.getB() );
               }
           }
       }
 
     /* boundary and interface condition */
 
     for (int n=0; n<numMeshes; n++)
     {
         const Mesh& mesh = *_meshes[n];
 
         foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
         {
             const FaceGroup& fg = *fgPtr;
             const StorageSite& faces = fg.site;
 
             typename SpeciesBCMap::const_iterator pos = sbcmap.find(fg.id);
             if (pos==sbcmap.end())
             {   
                throw CException("SpeciesModel: Error in BC Map");
             }
 
             const SpeciesBC<T>& bc = *(pos->second);
 
             GenericBCS<T,T,T> gbc(faces,mesh,
                                   _geomFields,
                                   sFields.massFraction,
                                   sFields.massFlux,
                                   ls.getMatrix(), ls.getX(), ls.getB());
 
             if (bc.bcType == "SpecifiedMassFraction")
             {
                 FloatValEvaluator<T>
                   bT(bc.getVal("specifiedMassFraction"),faces);
                 if (sFields.convectionFlux.hasArray(faces))
                 {
                     const TArray& convectingFlux =
                       dynamic_cast<const TArray&>
                       (sFields.convectionFlux[faces]);
                     const int nFaces = faces.getCount();
                                 
                     for(int f=0; f<nFaces; f++)
                     {
                         if (convectingFlux[f] > 0.)
                         {
                             gbc.applyExtrapolationBC(f);
                         }
                         else
                         {
                             gbc.applyDirichletBC(f,bT[f]);
                         }
                     }
                 }
                 else
                   gbc.applyDirichletBC(bT);
             }
             else if (bc.bcType == "SpecifiedMassFlux")
             {
                 const T specifiedFlux(bc["specifiedMassFlux"]);
                 gbc.applyNeumannBC(specifiedFlux);
             }
             else if ((bc.bcType == "Symmetry"))
             {
                  T zeroFlux(NumTypeTraits<T>::getZero());
                  gbc.applyNeumannBC(zeroFlux);
             }
             else
               throw CException(bc.bcType + " not implemented for SpeciesModel");
         }
 
         foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
         {
             const FaceGroup& fg = *fgPtr;
             const StorageSite& faces = fg.site;
             GenericBCS<T,T,T> gbc(faces,mesh,
                                   _geomFields,
                                   sFields.massFraction,
                                   sFields.massFlux,
                                   ls.getMatrix(), ls.getX(), ls.getB());
 
             gbc.applyInterfaceBC();
         }
     }
   }

template<class T>

void SpeciesModel< T >::Impl::printBCs ( )

inline

Definition at line 728 of file SpeciesModel_impl.h.

   {
     for (int m=0; m<_nSpecies; m++)
     {
       const SpeciesBCMap& sbcmap = *_bcMapVector[m];
       cout << "Species Number :" << m << endl;
 
       for(typename SpeciesBCMap::const_iterator pos=sbcmap.begin();
           pos!=sbcmap.end();
           ++pos)
       {
           cout << "Face Group " << pos->first << ":" << endl;
           cout << "    bc type " << pos->second->bcType << endl;
           foreach(typename SpeciesBC<T>::value_type& vp, *(pos->second))
           {
               cout << "   " << vp.first << " "  << vp.second.constant <<  endl;
           }
       }
     }
   }

template<class T>

void SpeciesModel< T >::Impl::updateTime ( )

inline

Definition at line 341 of file SpeciesModel_impl.h.

References Model::_meshes, Mesh::getCells(), SpeciesFields::massFraction, SpeciesFields::massFractionN1, and SpeciesFields::massFractionN2.

   {
     const int numMeshes = _meshes.size();
 
     for (int m=0; m<_nSpecies; m++)
     {
         SpeciesFields& sFields = *_speciesFieldsVector[m];
 
         for (int n=0; n<numMeshes; n++)
         {
             const Mesh& mesh = *_meshes[n];
 
             const StorageSite& cells = mesh.getCells();
             TArray& mF =
               dynamic_cast<TArray&>(sFields.massFraction[cells]);
             TArray& mFN1 =
               dynamic_cast<TArray&>(sFields.massFractionN1[cells]);
 
             if (_options.timeDiscretizationOrder > 1)
             {
                 TArray& mFN2 =
                   dynamic_cast<TArray&>(sFields.massFractionN2[cells]);
                 mFN2 = mFN1;
             }
             mFN1 = mF;
         }
     }
   }

Member Data Documentation

template<class T>

vector<SpeciesBCMap*> SpeciesModel< T >::Impl::_bcMapVector

private

Definition at line 767 of file SpeciesModel_impl.h.

template<class T>

vector<MFRPtr*> SpeciesModel< T >::Impl::_currentResidual

private

Definition at line 780 of file SpeciesModel_impl.h.

template<class T>

const GeomFields& SpeciesModel< T >::Impl::_geomFields

private

Definition at line 763 of file SpeciesModel_impl.h.

template<class T>

vector<MFRPtr*> SpeciesModel< T >::Impl::_initialNormVector

private

Definition at line 773 of file SpeciesModel_impl.h.

template<class T>

const MeshList SpeciesModel< T >::Impl::_meshes

private

Definition at line 762 of file SpeciesModel_impl.h.

template<class T>

int SpeciesModel< T >::Impl::_niters

private

Definition at line 774 of file SpeciesModel_impl.h.

template<class T>

const int SpeciesModel< T >::Impl::_nSpecies

private

Definition at line 775 of file SpeciesModel_impl.h.

template<class T>

SpeciesModelOptions<T> SpeciesModel< T >::Impl::_options

private

Definition at line 770 of file SpeciesModel_impl.h.

template<class T>

vector<SpeciesFields*> SpeciesModel< T >::Impl::_speciesFieldsVector

private

Definition at line 765 of file SpeciesModel_impl.h.

template<class T>

SpeciesModelFields SpeciesModel< T >::Impl::_speciesModelFields

private

Definition at line 777 of file SpeciesModel_impl.h.

template<class T>

vector<SpeciesVCMap*> SpeciesModel< T >::Impl::_vcMapVector

private

Definition at line 768 of file SpeciesModel_impl.h.

The documentation for this class was generated from the following file:

SpeciesModel_impl.h

Public Types

Public Member Functions

Private Attributes

Detailed Description

template<class T> class SpeciesModel< T >::Impl

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation

template<class T>
class SpeciesModel< T >::Impl