#include <FlowModel_impl.h>

Collaboration diagram for FlowModel< T >::Impl:

Public Types
typedef Array< int >	IntArray

typedef Array< T >	TArray

typedef Vector< T, 3 >	VectorT3

typedef VectorTranspose< T, 3 >	VectorT3T

typedef Array< VectorT3 >	VectorT3Array

typedef DiagonalTensor< T, 3 >	DiagTensorT3

typedef CRMatrix< DiagTensorT3, T, VectorT3 >	VVMatrix

typedef VVMatrix::DiagArray	VVDiagArray

typedef CRMatrix< T, T, T >	PPMatrix

typedef PPMatrix::DiagArray	PPDiagArray

typedef PPMatrix::PairWiseAssembler	PPAssembler

typedef Gradient< VectorT3 >	VGradType

typedef Array< Gradient < VectorT3 > >	VGradArray

typedef Gradient< T >	PGradType

typedef Array< PGradType >	PGradArray

typedef FluxJacobianMatrix< T, T >	FMatrix

typedef StressTensor< T >	StressTensorT6

typedef Array< StressTensorT6 >	StressTensorArray

Public Member Functions
	Impl (const GeomFields &geomFields, FlowFields &thermalFields, const MeshList &meshes)

void	init ()

FlowBCMap &	getBCMap ()

FlowVCMap &	getVCMap ()

FlowModelOptions< T > &	getOptions ()

void	updateTime ()

const Field &	getViscosityField () const

void	computeIBFaceVelocity (const StorageSite &particles)

map< string, shared_ptr < ArrayBase > > &	getPersistenceData ()

void	restart ()

void	initMomentumLinearization (LinearSystem &ls)

void	linearizeMomentum (LinearSystem &ls)

MFRPtr	solveMomentum ()

void	interfaceContinuityBC (const Mesh &mesh, const StorageSite &faces, MultiFieldMatrix &matrix, MultiField &rField)

void	correctVelocityBoundary (const Mesh &mesh, const StorageSite &faces, const MultiField &ppField)

void	correctMassFluxBoundary (const StorageSite &faces, const MultiField &ppField)

void	correctPressure (const Mesh &mesh, const MultiField &xField)

void	correctVelocityExplicit (const Mesh &mesh, const MultiField &xField)

void	setDirichlet (MultiFieldMatrix &mfmatrix, MultiField &rField)

void	linearizeContinuity (LinearSystem &ls)

void	setReferencePP (const MultiField &ppField)

void	computeContinuityResidual ()

void	postContinuitySolve (LinearSystem &ls)

void	initContinuityLinearization (LinearSystem &ls)

shared_ptr< LinearSystem >	discretizeContinuity ()

MFRPtr	solveContinuity ()

bool	advance (const int niter)

void	dumpContinuityMatrix (const string fileBase)

void	printBCs ()

VectorT3	getPressureIntegral (const Mesh &mesh, const int faceGroupId)

VectorT3	getPressureIntegralonIBFaces (const Mesh &mesh)

VectorT3	getMomentumFluxIntegral (const Mesh &mesh, const int faceGroupId)

VectorT3	getMomentumDerivativeIntegral (const Mesh &mesh)

boost::shared_ptr< ArrayBase >	getStressTensor (const Mesh &mesh, const ArrayBase &gcellIds)

void	ComputeStressTensorES (const StorageSite &solidFaces)

void	getTraction (const Mesh &mesh)

void	computeSolidSurfaceForce (const StorageSite &solidFaces, bool perUnitArea)

VectorT3	getMomentumFluxIntegralonIBFaces (const Mesh &mesh)

void	printPressureIntegrals ()

void	printMomentumFluxIntegrals ()

void	printMassFluxIntegrals ()

void	fixedPressureMomentumBC (const StorageSite &faces, const Mesh &mesh, MultiFieldMatrix &matrix, const MultiField &xField, MultiField &rField)

T	fixedPressureContinuityBC (const StorageSite &faces, const Mesh &mesh, MultiFieldMatrix &matrix, const MultiField &xField, MultiField &rField)

void	pressureBoundaryPostContinuitySolve (const StorageSite &faces, const Mesh &mesh, const T bp)

void	updateFacePressureBoundary (const Mesh &mesh, const StorageSite &faces)

T	fixedFluxContinuityBC (const StorageSite &faces, const Mesh &mesh, MultiFieldMatrix &matrix, MultiField &xField, MultiField &rField, const FlowBC< T > &bc)

T	discretizeMassFluxInterior (const Mesh &mesh, const StorageSite &faces, MultiFieldMatrix &mfmatrix, const MultiField &xField, MultiField &rField, const bool isSymmetry=false)

void	correctVelocityInterior (const Mesh &mesh, const StorageSite &faces, const MultiField &ppField, const bool isSymmetry=false)

void	updateFacePressureInterior (const Mesh &mesh, const StorageSite &faces, const bool isSymmetry=false)

void	correctMassFluxInterior (const Mesh &mesh, const StorageSite &faces, MultiFieldMatrix &mfmatrix, const MultiField &xField)

void	slipJumpMomentumBC (const StorageSite &faces, const Mesh &mesh, GenericBCS< VectorT3, DiagTensorT3, T > &gbc, const T accomodationCoefficient, const FloatValEvaluator< VectorT3 > &bVelocity)

Private Attributes
const MeshList	_meshes

const GeomFields &	_geomFields

FlowFields &	_flowFields

FlowBCMap	_bcMap

FlowVCMap	_vcMap

FlowModelOptions< T >	_options

GradientModel< VectorT3 >	_velocityGradientModel

GradientModel< T >	_pressureGradientModel

MFRPtr	_initialMomentumNorm

MFRPtr	_initialContinuityNorm

MFRPtr	_initialCoupledNorm

int	_niters

shared_ptr< Field >	_previousVelocity

shared_ptr< Field >	_momApField

bool	_useReferencePressure

int	_globalRefCellID

int	_globalRefProcID

T	_referencePP

map< string, shared_ptr < ArrayBase > >	_persistenceData

Detailed Description

template<class T>
class FlowModel< T >::Impl

Definition at line 41 of file FlowModel_impl.h.

Member Typedef Documentation

template<class T>

typedef DiagonalTensor<T,3> FlowModel< T >::Impl::DiagTensorT3

Definition at line 52 of file FlowModel_impl.h.

template<class T>

typedef FluxJacobianMatrix<T,T> FlowModel< T >::Impl::FMatrix

Definition at line 78 of file FlowModel_impl.h.

template<class T>

typedef Array<int> FlowModel< T >::Impl::IntArray

Definition at line 45 of file FlowModel_impl.h.

template<class T>

typedef Array<PGradType> FlowModel< T >::Impl::PGradArray

Definition at line 76 of file FlowModel_impl.h.

template<class T>

typedef Gradient<T> FlowModel< T >::Impl::PGradType

Definition at line 75 of file FlowModel_impl.h.

template<class T>

typedef PPMatrix::PairWiseAssembler FlowModel< T >::Impl::PPAssembler

Definition at line 59 of file FlowModel_impl.h.

template<class T>

typedef PPMatrix::DiagArray FlowModel< T >::Impl::PPDiagArray

Definition at line 58 of file FlowModel_impl.h.

template<class T>

typedef CRMatrix<T,T,T> FlowModel< T >::Impl::PPMatrix

Definition at line 57 of file FlowModel_impl.h.

template<class T>

typedef Array<StressTensorT6> FlowModel< T >::Impl::StressTensorArray

Definition at line 81 of file FlowModel_impl.h.

template<class T>

typedef StressTensor<T> FlowModel< T >::Impl::StressTensorT6

Definition at line 80 of file FlowModel_impl.h.

template<class T>

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

Definition at line 47 of file FlowModel_impl.h.

template<class T>

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

Definition at line 48 of file FlowModel_impl.h.

template<class T>

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

Definition at line 51 of file FlowModel_impl.h.

template<class T>

typedef VectorTranspose<T,3> FlowModel< T >::Impl::VectorT3T

Definition at line 49 of file FlowModel_impl.h.

template<class T>

typedef Array<Gradient<VectorT3> > FlowModel< T >::Impl::VGradArray

Definition at line 73 of file FlowModel_impl.h.

template<class T>

typedef Gradient<VectorT3> FlowModel< T >::Impl::VGradType

Definition at line 72 of file FlowModel_impl.h.

template<class T>

typedef VVMatrix::DiagArray FlowModel< T >::Impl::VVDiagArray

Definition at line 55 of file FlowModel_impl.h.

template<class T>

typedef CRMatrix<DiagTensorT3,T,VectorT3> FlowModel< T >::Impl::VVMatrix

Definition at line 54 of file FlowModel_impl.h.

Constructor & Destructor Documentation

template<class T>

FlowModel< T >::Impl::Impl	(	const GeomFields &	geomFields,
		FlowFields &	thermalFields,
		const MeshList &	meshes
	)

inline

Definition at line 83 of file FlowModel_impl.h.

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

                               :
     _meshes(meshes),
     _geomFields(geomFields),
     _flowFields(thermalFields),
     _velocityGradientModel(_meshes,_flowFields.velocity,
                            _flowFields.velocityGradient,_geomFields),
     _pressureGradientModel(_meshes,_flowFields.pressure,
                            _flowFields.pressureGradient,_geomFields),
     _initialMomentumNorm(),
     _initialContinuityNorm(),
     _niters(0)
   {
     const int numMeshes = _meshes.size();
     for (int n=0; n<numMeshes; n++)
     {
         const Mesh& mesh = *_meshes[n];
 
         FlowVC<T> *vc(new FlowVC<T>());
         vc->vcType = "flow";
         _vcMap[mesh.getID()] = vc;
         foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
         {
             const FaceGroup& fg = *fgPtr;
             if (_bcMap.find(fg.id) == _bcMap.end())
             {
                 FlowBC<T> *bc(new FlowBC<T>());
                 
                 _bcMap[fg.id] = bc;
                 if ((fg.groupType == "wall"))
                 {
                     bc->bcType = "NoSlipWall";
                 }
                 else if ((fg.groupType == "velocity-inlet"))
                 {
                     bc->bcType = "VelocityBoundary";
                 }
                 else if ((fg.groupType == "pressure-inlet") ||
                          (fg.groupType == "pressure-outlet"))
                 {
                     bc->bcType = "PressureBoundary";
                 }
                 else if ((fg.groupType == "symmetry"))
                 {
                     bc->bcType = "Symmetry";
                 }
                 else
                   throw CException("FlowModel: unknown face group type "
                                    + fg.groupType);
             }
         }
         /*
         foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
           {
             const FaceGroup& fg = *fgPtr;
             if ((fg.groupType == "ESinterface"))
               { bc->bcType = "ESInterfaceBC";}
           }
         */
         
 
     }
   }

Member Function Documentation

template<class T>

bool FlowModel< T >::Impl::advance ( const int niter )

inline

Definition at line 1433 of file FlowModel_impl.h.

   {
 
     for(int n=0; n<niter; n++)
     { 
         MFRPtr mNorm = solveMomentum();
         MFRPtr cNorm = solveContinuity();
 
         if (_niters < 5)
         {
             _initialMomentumNorm->setMax(*mNorm);
             _initialContinuityNorm->setMax(*cNorm);
         }
         
         MFRPtr mNormRatio((*mNorm)/(*_initialMomentumNorm));
         MFRPtr cNormRatio((*cNorm)/(*_initialContinuityNorm));
         
 #ifdef FVM_PARALLEL
         if ( MPI::COMM_WORLD.Get_rank() == 0 ){ // only root process
         if (_options.printNormalizedResiduals)
           {cout << _niters << ": " << *mNormRatio << ";" << *cNormRatio <<  endl;}
         else{
           cout << _niters << ": " << *mNorm << ";" << *cNorm <<  endl;}}
 #endif
 
 #ifndef FVM_PARALLEL
         if (_options.printNormalizedResiduals)
           {cout << _niters << ": " << *mNormRatio << ";" << *cNormRatio <<  endl;}
         else
           {cout << _niters << ": " << *mNorm << ";" << *cNorm <<  endl;}
 #endif
 
         _niters++;
         if ((*mNormRatio < _options.momentumTolerance) &&
             (*cNormRatio < _options.continuityTolerance))
           return true;
     }
     return false;
   }

template<class T>

void FlowModel< T >::Impl::computeContinuityResidual ( )

inline

Definition at line 1235 of file FlowModel_impl.h.

References Model::_meshes, Mesh::getAllFaceCells(), Mesh::getCells(), StorageSite::getCount(), Mesh::getFaces(), and Array< T >::zero().

   {
     const int numMeshes = _meshes.size();
     for (int n=0; n<numMeshes; n++)
     {
         const Mesh& mesh = *_meshes[n];
 
         const StorageSite& cells = mesh.getCells();
         const StorageSite& faces = mesh.getFaces();
 
         TArray& r = dynamic_cast<TArray&>(_flowFields.continuityResidual[cells]);
         const TArray& massFlux = dynamic_cast<const TArray&>(_flowFields.massFlux[faces]);
 
         const CRConnectivity& faceCells = mesh.getAllFaceCells();
         const int nFaces = faces.getCount();
 
         r.zero();
         for(int f=0; f<nFaces; f++)
         {
           const int c0 = faceCells(f,0);
           const int c1 = faceCells(f,1);
 
             r[c0] += massFlux[f];
             r[c1] -= massFlux[f];
         }
     }
   }

template<class T>

void FlowModel< T >::Impl::computeIBFaceVelocity ( const StorageSite & particles )

inline

Definition at line 378 of file FlowModel_impl.h.

References Model::_meshes, Mesh::getCells(), StorageSite::getCount(), Mesh::getFaces(), Mesh::getIBFaceList(), Mesh::getIBFaces(), Mesh::isShell(), and Array< T >::zero().

   {
     typedef CRMatrixTranspose<T,T,T> IMatrix;
     typedef CRMatrixTranspose<T,VectorT3,VectorT3> IMatrixV3;
 
     const VectorT3Array& pV =
       dynamic_cast<const VectorT3Array&>(_flowFields.velocity[particles]);
 
     const int numMeshes = _meshes.size();
     for (int n=0; n<numMeshes; n++)
     {
         const Mesh& mesh = *_meshes[n];
         if (!mesh.isShell() && mesh.getIBFaces().getCount() > 0){
 
            const StorageSite& cells = mesh.getCells();
            const StorageSite& ibFaces = mesh.getIBFaces();
         
            GeomFields::SSPair key1(&ibFaces,&cells);
            const IMatrix& mIC =
            dynamic_cast<const IMatrix&>
             (*_geomFields._interpolationMatrices[key1]);
 
            IMatrixV3 mICV(mIC);
            
 
            GeomFields::SSPair key2(&ibFaces,&particles);
            const IMatrix& mIP =
            dynamic_cast<const IMatrix&>
             (*_geomFields._interpolationMatrices[key2]);
 
            IMatrixV3 mIPV(mIP);
            
 
            shared_ptr<VectorT3Array> ibV(new VectorT3Array(ibFaces.getCount()));
         
            const VectorT3Array& cV =
             dynamic_cast<const VectorT3Array&>(_flowFields.velocity[cells]);
     
 
            ibV->zero();
 
            mICV.multiplyAndAdd(*ibV,cV);
            mIPV.multiplyAndAdd(*ibV,pV);
 
 #if 0
         ofstream debugFile;
         stringstream ss(stringstream::in | stringstream::out);
         ss <<  MPI::COMM_WORLD.Get_rank();
         string  fname1 = "IBVelocity_proc" +  ss.str() + ".dat";
         debugFile.open(fname1.c_str());
         
         //debug use
         const Array<int>& ibFaceList = mesh.getIBFaceList();
         const StorageSite& faces = mesh.getFaces();
         const VectorT3Array& faceCentroid = 
           dynamic_cast<const VectorT3Array&> (_geomFields.coordinate[faces]);
         const double angV = 1.0;
         VectorT3 center;
         center[0]=0.;
         center[1]=0.;
         center[2]=0.;   
 
         for(int f=0; f<ibFaces.getCount();f++){
           int fID = ibFaceList[f];
           debugFile << "f=" <<   f << setw(10) <<  "   fID = " <<  fID << "  faceCentroid = " << faceCentroid[fID] << " ibV = " << (*ibV)[f] << endl;
         }
           
          debugFile.close();
 #endif
 
 
           _flowFields.velocity.addArray(ibFaces,ibV);
 
         }
     }
 
   }

template<class T>

void FlowModel< T >::Impl::computeSolidSurfaceForce	(	const StorageSite &	solidFaces,
		bool	perUnitArea
	)

inline

Definition at line 1939 of file FlowModel_impl.h.

References Model::_meshes, Mesh::getCells(), CRConnectivity::getCol(), Mesh::getConnectivity(), StorageSite::getCount(), CRConnectivity::getRow(), and Array< T >::zero().

   {
     typedef CRMatrixTranspose<T,T,T> IMatrix;
 
     const int nSolidFaces = solidFaces.getCount();
 
     _velocityGradientModel.compute();
 
     boost::shared_ptr<VectorT3Array>
       forcePtr( new VectorT3Array(nSolidFaces));
     VectorT3Array& force = *forcePtr;
 
     force.zero();
     _flowFields.force.addArray(solidFaces,forcePtr);
 
     const VectorT3Array& solidFaceArea =
       dynamic_cast<const VectorT3Array&>(_geomFields.area[solidFaces]);
 
     const TArray& solidFaceAreaMag =
       dynamic_cast<const TArray&>(_geomFields.areaMag[solidFaces]);
     
     const int numMeshes = _meshes.size();
     for (int n=0; n<numMeshes; n++)
     {
         const Mesh& mesh = *_meshes[n];
         const StorageSite& cells = mesh.getCells();
         const VGradArray& vGrad =
           dynamic_cast<const VGradArray&>(_flowFields.velocityGradient[cells]);
 
         const TArray& pCell =
           dynamic_cast<const TArray&>(_flowFields.pressure[cells]);
         
         const TArray& mu =
           dynamic_cast<const TArray&>(getViscosityField()[cells]);
         
         //const FlowVC<T>& vc = *_vcMap[mesh.getID()];
             
         const CRConnectivity& solidFacesToCells
           = mesh.getConnectivity(solidFaces,cells);
         
         const IntArray& sFCRow = solidFacesToCells.getRow();
         const IntArray& sFCCol = solidFacesToCells.getCol();
 
         GeomFields::SSPair key1(&solidFaces,&cells);
         const IMatrix& mIC =
           dynamic_cast<const IMatrix&>
           (*_geomFields._interpolationMatrices[key1]);
 
         const Array<T>& iCoeffs = mIC.getCoeff();
         
         for(int f=0; f<nSolidFaces; f++)
         {
             StressTensor<T> stress = NumTypeTraits<StressTensor<T> >::getZero();
             for(int nc = sFCRow[f]; nc<sFCRow[f+1]; nc++)
             {
                 const int c = sFCCol[nc];
                 const VGradType& vg = vGrad[c];
                 VGradType vgPlusTranspose = vGrad[c];
                 
                 for(int i=0;i<3;i++)
                   for(int j=0;j<3;j++)
                     vgPlusTranspose[i][j] += vg[j][i];
 
                 const T coeff = iCoeffs[nc];
                
                 stress[0] += coeff*(vgPlusTranspose[0][0]*mu[c] - pCell[c]);
                 stress[1] += coeff*(vgPlusTranspose[1][1]*mu[c] - pCell[c]);
                 stress[2] += coeff*(vgPlusTranspose[2][2]*mu[c] - pCell[c]);
                 stress[3] += coeff*(vgPlusTranspose[0][1]*mu[c]);
                 stress[4] += coeff*(vgPlusTranspose[1][2]*mu[c]);
                 stress[5] += coeff*(vgPlusTranspose[2][0]*mu[c]);
             }
 
             const VectorT3& Af = solidFaceArea[f];
             force[f][0] = Af[0]*stress[0] + Af[1]*stress[3] + Af[2]*stress[5];
             force[f][1] = Af[0]*stress[3] + Af[1]*stress[1] + Af[2]*stress[4];
             force[f][2] = Af[0]*stress[5] + Af[1]*stress[4] + Af[2]*stress[2];
             if (perUnitArea)
             {
                 force[f] /= solidFaceAreaMag[f];
             }
         }
     }
   }

template<class T>

void FlowModel< T >::Impl::ComputeStressTensorES ( const StorageSite & solidFaces )

inline

Definition at line 1811 of file FlowModel_impl.h.

References Model::_meshes, Mesh::getCells(), StorageSite::getCount(), and Mesh::getFaceCells().

   {
     typedef Array<StressTensor<T> > StressTensorArray;
     const int nSolidFaces = solidFaces.getCount();
 
 
     //interface
     VectorT3Array& IntVel =
       dynamic_cast<VectorT3Array&>(_flowFields.InterfaceVelocity[solidFaces]);
     TArray& IntPress =
       dynamic_cast<TArray&>(_flowFields.InterfacePressure[solidFaces]); 
     TArray& IntDens =
       dynamic_cast<TArray&>(_flowFields.InterfaceDensity[solidFaces]);
    
     //const TArray& mu = dynamic_cast<const TArray&>(getViscosityField()[cells]);   
     StressTensorArray& stressTensor = dynamic_cast<StressTensorArray &>(_flowFields.InterfaceStress[solidFaces]);
    
 
     const int numMeshes = _meshes.size();
  
     for (int n=0; n<numMeshes; n++)
       {
         const Mesh& mesh = *_meshes[n];
         const StorageSite& cells = mesh.getCells();
         //const int nCells = cells.getCount();
         
         _velocityGradientModel.compute();
 
         const VGradArray& vGrad =
           dynamic_cast<VGradArray&>(_flowFields.velocityGradient[cells]);
         const TArray& pCell =
           dynamic_cast<const TArray&>(_flowFields.pressure[cells]);
         const TArray& dCell =
           dynamic_cast<const TArray&>(_flowFields.density[cells]);
         const VectorT3Array& vCell =
           dynamic_cast<VectorT3Array&>(_flowFields.velocity[cells]);
 
 
         const CRConnectivity& _faceCells=mesh.getFaceCells(solidFaces);
         
         for(int f=0; f<nSolidFaces; f++)
           {
            
             const int c0 = _faceCells(f,0);  
             const int c1 = _faceCells(f,1);     
             
             const VGradType& vg = vGrad[c0];
 
             VGradType vgPlusTranspose = vGrad[c0];
             
             for(int i=0;i<3;i++)
               for(int j=0;j<3;j++)
                 vgPlusTranspose[i][j] += vg[j][i];
             
             stressTensor[f][0] = vgPlusTranspose[0][0];
             stressTensor[f][1] = vgPlusTranspose[1][1];
             stressTensor[f][3] = vgPlusTranspose[0][1];
             stressTensor[f][4] = vgPlusTranspose[1][2];
             stressTensor[f][5] = vgPlusTranspose[2][0];
             
             //copy density pressure and velocity into Interface arrays
             IntDens[f]=dCell[c1];
             IntPress[f]=pCell[c1];
             IntVel[f][0]=vCell[c1][0];
             IntVel[f][1]=vCell[c1][1];
             IntVel[f][2]=vCell[c1][2];
               }
       }
     
     
   }

template<class T>

void FlowModel< T >::Impl::correctMassFluxBoundary	(	const StorageSite &	faces,
		const MultiField &	ppField
	)

inline

Definition at line 835 of file FlowModel_impl.h.

References StorageSite::getCount().

   {
     MultiField::ArrayIndex mfIndex(&_flowFields.massFlux,&faces);
     const TArray& dMassFlux = dynamic_cast<const TArray&>(ppField[mfIndex]);
     TArray& massFlux = dynamic_cast<TArray&>(_flowFields.massFlux[faces]);
     
     const int nFaces = faces.getCount();
     for(int f=0; f<nFaces; f++)
     {
         massFlux[f] -= dMassFlux[f];
     }
   }

template<class T>

void FlowModel< T >::Impl::correctMassFluxInterior	(	const Mesh &	mesh,
		const StorageSite &	faces,
		MultiFieldMatrix &	mfmatrix,
		const MultiField &	xField
	)

inline

Definition at line 372 of file FlowModel_impl.h.

   {
     typedef CRMatrixTranspose<T,T,T> IMatrix;
     typedef CRMatrixTranspose<T,VectorT3,VectorT3> IMatrixV3;
 
     const VectorT3Array& pV =
       dynamic_cast<const VectorT3Array&>(_flowFields.velocity[particles]);
 
     const int numMeshes = _meshes.size();
     for (int n=0; n<numMeshes; n++)
     {
         const Mesh& mesh = *_meshes[n];
         if (!mesh.isShell() && mesh.getIBFaces().getCount() > 0){
 
            const StorageSite& cells = mesh.getCells();
            const StorageSite& ibFaces = mesh.getIBFaces();
         
            GeomFields::SSPair key1(&ibFaces,&cells);
            const IMatrix& mIC =
            dynamic_cast<const IMatrix&>
             (*_geomFields._interpolationMatrices[key1]);
 
            IMatrixV3 mICV(mIC);
            
 
            GeomFields::SSPair key2(&ibFaces,&particles);
            const IMatrix& mIP =
            dynamic_cast<const IMatrix&>
             (*_geomFields._interpolationMatrices[key2]);
 
            IMatrixV3 mIPV(mIP);
            
 
            shared_ptr<VectorT3Array> ibV(new VectorT3Array(ibFaces.getCount()));
         
            const VectorT3Array& cV =
             dynamic_cast<const VectorT3Array&>(_flowFields.velocity[cells]);
     
 
            ibV->zero();
 
            mICV.multiplyAndAdd(*ibV,cV);
            mIPV.multiplyAndAdd(*ibV,pV);
 
 #if 0

template<class T>

void FlowModel< T >::Impl::correctPressure	(	const Mesh &	mesh,
		const MultiField &	xField
	)

inline

Definition at line 849 of file FlowModel_impl.h.

References Mesh::getCells(), and StorageSite::getCountLevel1().

   {
     const StorageSite& cells = mesh.getCells();
     MultiField::ArrayIndex pIndex(&_flowFields.pressure,&cells);
     
     TArray& p = dynamic_cast<TArray&>(_flowFields.pressure[cells]);
     const TArray& pp = dynamic_cast<const TArray&>(xField[pIndex]);
 
     const T pressureURF(_options["pressureURF"]);
       
     const int nCells = cells.getCountLevel1();
     for(int c=0; c<nCells; c++)
     {
         p[c] += pressureURF*(pp[c]-_referencePP);
     }
   }

template<class T>

void FlowModel< T >::Impl::correctVelocityBoundary	(	const Mesh &	mesh,
		const StorageSite &	faces,
		const MultiField &	ppField
	)

inline

Definition at line 804 of file FlowModel_impl.h.

References Mesh::getCells(), StorageSite::getCount(), and Mesh::getFaceCells().

   {
     const StorageSite& cells = mesh.getCells();
 
     const VectorT3Array& faceArea = dynamic_cast<const VectorT3Array&>(_geomFields.area[faces]);
 
     const CRConnectivity& faceCells = mesh.getFaceCells(faces);
 
     MultiField::ArrayIndex pIndex(&_flowFields.pressure,&cells);
     MultiField::ArrayIndex vIndex(&_flowFields.velocity,&cells);
     const VVDiagArray& momAp = dynamic_cast<const VVDiagArray&>((*_momApField)[cells]);
     
     VectorT3Array& V = dynamic_cast<VectorT3Array&>(_flowFields.velocity[cells]);
     const TArray& pp = dynamic_cast<const TArray&>(ppField[pIndex]);
     
     const int nFaces = faces.getCount();
     for(int f=0; f<nFaces; f++)
     {
         const int c0 = faceCells(f,0);
         const int c1 = faceCells(f,1);
 
         const T ppFace = pp[c1];
         const VectorT3 ppA = ppFace*faceArea[f];
 
         V[c0] += ppA/momAp[c0];
     }
   }

template<class T>

void FlowModel< T >::Impl::correctVelocityExplicit	(	const Mesh &	mesh,
		const MultiField &	xField
	)

inline

Definition at line 868 of file FlowModel_impl.h.

References Mesh::getBoundaryFaceGroups(), Mesh::getCells(), StorageSite::getCount(), StorageSite::getCountLevel1(), and FaceGroup::site.

   {
     const StorageSite& cells = mesh.getCells();
     MultiField::ArrayIndex vIndex(&_flowFields.velocity,&cells);
     
     VectorT3Array& V = dynamic_cast<VectorT3Array&>(_flowFields.velocity[cells]);
     const VectorT3Array& Vp = dynamic_cast<const VectorT3Array&>(xField[vIndex]);
 
     const T velocityURF(_options["velocityURF"]);
       
     const int nCells = cells.getCountLevel1();
     for(int c=0; c<nCells; c++)
     {
         V[c] += velocityURF*Vp[c];
     }
 
     // boundary
     foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
     {
         const FaceGroup& fg = *fgPtr;
         const StorageSite& faces = fg.site;
         MultiField::ArrayIndex fluxIndex(&_flowFields.momentumFlux,&faces);
         VectorT3Array& momFlux =
           dynamic_cast<VectorT3Array&>(_flowFields.momentumFlux[faces]);
         const VectorT3Array& dmomFlux =
           dynamic_cast<const VectorT3Array&>(xField[fluxIndex]);
         
         const int nFaces = faces.getCount();
         for(int f=0; f<nFaces; f++)
         {
             momFlux[f] += dmomFlux[f];
         }
     }
   }

template<class T>

void FlowModel< T >::Impl::correctVelocityInterior	(	const Mesh &	mesh,
		const StorageSite &	faces,
		const MultiField &	ppField,
		const bool	isSymmetry = `false`
	)

inline

Definition at line 222 of file FlowModel_impl.h.

         {
             const int c0 = faceCells(f,0);
             const int c1 = faceCells(f,1);
             
             mf[f] = 0.5*(density[c0]*dot(V[c0],faceArea[f]) +
                          density[c1]*dot(V[c1],faceArea[f]));
         }
         _flowFields.massFlux.addArray(faces,mfPtr);
 
         // store momentum flux at interfaces
         foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
         {
             const FaceGroup& fg = *fgPtr;
             const StorageSite& faces = fg.site;
             shared_ptr<VectorT3Array> momFlux(new VectorT3Array(faces.getCount()));
             momFlux->zero();
             _flowFields.momentumFlux.addArray(faces,momFlux);
             
             if(fg.groupType == "ESinterface"){
               cout << "interface init" <<endl;
               const StorageSite& Intfaces = fg.site; 
               
               shared_ptr<VectorT3Array> InterfaceVelFace(new VectorT3Array(Intfaces.getCount()));
               InterfaceVelFace ->zero();
               _flowFields.InterfaceVelocity.addArray(Intfaces,InterfaceVelFace);
               
               shared_ptr<StressTensorArray> InterfaceStressFace(new StressTensorArray(Intfaces.getCount()));
               InterfaceStressFace ->zero();
               _flowFields.InterfaceStress.addArray(Intfaces,InterfaceStressFace);
               
               shared_ptr<TArray> InterfacePressFace(new TArray(Intfaces.getCount()));
               *InterfacePressFace = _options["initialPressure"];
               _flowFields.InterfacePressure.addArray(Intfaces,InterfacePressFace);
               
               shared_ptr<TArray> InterfaceDensityFace(new TArray(Intfaces.getCount()));
               *InterfaceDensityFace =vc["density"];
               _flowFields.InterfaceDensity.addArray(Intfaces,InterfaceDensityFace);
               
               
             }
             
             
         }
         
         foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
         {
             const FaceGroup& fg = *fgPtr;
             const StorageSite& faces = fg.site;
 
             const FlowBC<T>& bc = *_bcMap[fg.id];
 
             FloatValEvaluator<VectorT3>
               bVelocity(bc.getVal("specifiedXVelocity"),
                         bc.getVal("specifiedYVelocity"),
                         bc.getVal("specifiedZVelocity"),
                         faces);
 
             const int nFaces = faces.getCount();
             const CRConnectivity& faceCells = mesh.getFaceCells(faces);
 
             if ((bc.bcType == "NoSlipWall") ||
                 (bc.bcType == "SlipJump") ||
                 (bc.bcType == "VelocityBoundary"))
             {
                 // arrays for this face group
                 TArray& massFlux = dynamic_cast<TArray&>(_flowFields.massFlux[faces]);

template<class T>

shared_ptr<LinearSystem> FlowModel< T >::Impl::discretizeContinuity ( )

inline

Definition at line 1394 of file FlowModel_impl.h.

   {
     shared_ptr<LinearSystem> ls(new LinearSystem());
 
     initContinuityLinearization(*ls);
         
     ls->initAssembly();
     
     linearizeContinuity(*ls);
 
     ls->initSolve();
 
     return ls;
   }

template<class T>

T FlowModel< T >::Impl::discretizeMassFluxInterior	(	const Mesh &	mesh,
		const StorageSite &	faces,
		MultiFieldMatrix &	mfmatrix,
		const MultiField &	xField,
		MultiField &	rField,
		const bool	isSymmetry = `false`
	)

inline

Definition at line 9 of file FlowModel_impl.h.

 {
 public:
 
   typedef Array<int> IntArray;
   
   typedef Array<T> TArray;
   typedef Vector<T,3> VectorT3;
   typedef VectorTranspose<T,3> VectorT3T;
 
   typedef Array<VectorT3> VectorT3Array;
   typedef DiagonalTensor<T,3> DiagTensorT3;
 
   typedef CRMatrix<DiagTensorT3,T,VectorT3> VVMatrix;
   typedef typename VVMatrix::DiagArray VVDiagArray;
   
   typedef CRMatrix<T,T,T> PPMatrix;
   typedef typename PPMatrix::DiagArray PPDiagArray;
   typedef typename PPMatrix::PairWiseAssembler PPAssembler;
 
 #ifdef PV_COUPLED
   typedef CRMatrixRect<VectorT3T,VectorT3,T> PVMatrix;
   typedef typename PVMatrix::DiagArray PVDiagArray;
   typedef typename PVMatrix::PairWiseAssembler PVAssembler;
 
   typedef CRMatrixRect<VectorT3,T,VectorT3> VPMatrix;
   typedef typename VPMatrix::DiagArray VPDiagArray;
   typedef typename VPMatrix::PairWiseAssembler VPAssembler;
 
 #endif
   
   typedef Gradient<VectorT3> VGradType;
   typedef Array<Gradient<VectorT3> > VGradArray;
 
   typedef Gradient<T> PGradType;
   typedef Array<PGradType> PGradArray;
 
   typedef FluxJacobianMatrix<T,T> FMatrix;
   
   typedef StressTensor<T> StressTensorT6;
   typedef Array<StressTensorT6> StressTensorArray;
   
   Impl(const GeomFields& geomFields,
        FlowFields& thermalFields,
        const MeshList& meshes):
     _meshes(meshes),
     _geomFields(geomFields),
     _flowFields(thermalFields),
     _velocityGradientModel(_meshes,_flowFields.velocity,
                            _flowFields.velocityGradient,_geomFields),
     _pressureGradientModel(_meshes,_flowFields.pressure,
                            _flowFields.pressureGradient,_geomFields),
     _initialMomentumNorm(),
     _initialContinuityNorm(),
     _niters(0)
   {
     const int numMeshes = _meshes.size();
     for (int n=0; n<numMeshes; n++)
     {
         const Mesh& mesh = *_meshes[n];
 
         FlowVC<T> *vc(new FlowVC<T>());
         vc->vcType = "flow";
         _vcMap[mesh.getID()] = vc;
         foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
         {
             const FaceGroup& fg = *fgPtr;
             if (_bcMap.find(fg.id) == _bcMap.end())
             {
                 FlowBC<T> *bc(new FlowBC<T>());
                 
                 _bcMap[fg.id] = bc;
                 if ((fg.groupType == "wall"))
                 {
                     bc->bcType = "NoSlipWall";
                 }
                 else if ((fg.groupType == "velocity-inlet"))
                 {
                     bc->bcType = "VelocityBoundary";
                 }
                 else if ((fg.groupType == "pressure-inlet") ||
                          (fg.groupType == "pressure-outlet"))
                 {
                     bc->bcType = "PressureBoundary";
                 }
                 else if ((fg.groupType == "symmetry"))
                 {
                     bc->bcType = "Symmetry";
                 }
                 else
                   throw CException("FlowModel: unknown face group type "
                                    + fg.groupType);
             }
         }
         /*
         foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
           {
             const FaceGroup& fg = *fgPtr;
             if ((fg.groupType == "ESinterface"))
               { bc->bcType = "ESInterfaceBC";}
           }
         */
         
 
     }
   }
 
   void init()
   {
     const int numMeshes = _meshes.size();
     for (int n=0; n<numMeshes; n++)
     {
         const Mesh& mesh = *_meshes[n];
 
         const FlowVC<T>& vc = *_vcMap[mesh.getID()];
             
         const StorageSite& cells = mesh.getCells();
         const StorageSite& faces = mesh.getFaces();
 
         shared_ptr<VectorT3Array> vCell(new VectorT3Array(cells.getCountLevel1()));
 
         VectorT3 initialVelocity;
         initialVelocity[0] = _options["initialXVelocity"];
         initialVelocity[1] = _options["initialYVelocity"];
         initialVelocity[2] = _options["initialZVelocity"];
         *vCell = initialVelocity;
 
         _flowFields.velocity.addArray(cells,vCell);
 
         shared_ptr<VectorT3Array> uparallelCell(new VectorT3Array(cells.getCount()));
         uparallelCell ->zero();
         _flowFields.uparallel.addArray(cells,uparallelCell);
 
       shared_ptr<VectorT3Array> tauCell(new VectorT3Array(faces.getCount()));
         tauCell ->zero();
         _flowFields.tau.addArray(faces,tauCell);
 
      shared_ptr<VectorT3Array> TauWallCell(new VectorT3Array(faces.getCount()));
         TauWallCell ->zero();
         _flowFields.tauwall.addArray(faces,TauWallCell);
 
 
 
 
   
        
 
         if (_options.transient)
         {
             _flowFields.velocityN1.addArray(cells,
                                             dynamic_pointer_cast<ArrayBase>(vCell->newCopy()));
             if (_options.timeDiscretizationOrder > 1)
               _flowFields.velocityN2.addArray(cells,
                                               dynamic_pointer_cast<ArrayBase>(vCell->newCopy()));
 
         }
         
         shared_ptr<TArray> pCell(new TArray(cells.getCountLevel1()));
         shared_ptr<TArray> pFace(new TArray(faces.getCountLevel1()));
         *pCell = _options["initialPressure"];
         *pFace = _options["initialPressure"];
         _flowFields.pressure.addArray(cells,pCell);
         _flowFields.pressure.addArray(faces,pFace);
 
 
         shared_ptr<TArray> rhoCell(new TArray(cells.getCountLevel1()));
         *rhoCell = vc["density"];
         _flowFields.density.addArray(cells,rhoCell);
 
         shared_ptr<TArray> muCell(new TArray(cells.getCountLevel1()));
         *muCell = vc["viscosity"];
         _flowFields.viscosity.addArray(cells,muCell);
 
         shared_ptr<PGradArray> gradp(new PGradArray(cells.getCountLevel1()));
         gradp->zero();
         _flowFields.pressureGradient.addArray(cells,gradp);
 
         shared_ptr<TArray> ci(new TArray(cells.getCountLevel1()));
         ci->zero();

template<class T>

void FlowModel< T >::Impl::dumpContinuityMatrix ( const string fileBase )

inline

Definition at line 1560 of file FlowModel_impl.h.

References Model::_meshes, Mesh::getCells(), CRConnectivity::getCol(), CRMatrix< T_Diag, T_OffDiag, X >::getConnectivity(), CRMatrix< T_Diag, T_OffDiag, X >::getDiag(), MultiFieldMatrix::getMatrix(), CRMatrix< T_Diag, T_OffDiag, X >::getOffDiag(), CRConnectivity::getRow(), and StorageSite::getSelfCount().

   {
     solveMomentum();
     shared_ptr<LinearSystem> ls(discretizeContinuity());
 
     MultiFieldMatrix& matrix = ls->getMatrix();
     MultiField& b = ls->getB();
 
     const Mesh& mesh = *_meshes[0];
     const StorageSite& cells = mesh.getCells();
     
     MultiField::ArrayIndex pIndex(&_flowFields.pressure,&cells);
 
     PPMatrix& ppMatrix =
       dynamic_cast<PPMatrix&>(matrix.getMatrix(pIndex,pIndex));
 
     PPDiagArray& ppDiag = ppMatrix.getDiag();
     PPDiagArray& ppCoeff = ppMatrix.getOffDiag();
 
     TArray& rCell = dynamic_cast<TArray&>(b[pIndex]);
 
     const CRConnectivity& cr = ppMatrix.getConnectivity();
 
     const Array<int>& row = cr.getRow();
     const Array<int>& col = cr.getCol();
     
     const int nCells = cells.getSelfCount();
     int nCoeffs = nCells;
 
     for(int i=0; i<nCells; i++)
       for(int jp=row[i]; jp<row[i+1]; jp++)
       {
           const int j = col[jp];
           if (j<nCells) nCoeffs++;
       }
     
     string matFileName = fileBase + ".mat";
     FILE *matFile = fopen(matFileName.c_str(),"wb");
     
     fprintf(matFile,"%%%%MatrixMarket matrix coordinate real general\n");
     fprintf(matFile,"%d %d %d\n", nCells,nCells,nCoeffs);
 
     for(int i=0; i<nCells; i++)
     {
         fprintf(matFile,"%d %d %lf\n", i+1, i+1, ppDiag[i]);
         for(int jp=row[i]; jp<row[i+1]; jp++)
         {
             const int j = col[jp];
             if (j<nCells)
               fprintf(matFile,"%d %d %lf\n", i+1, j+1, ppCoeff[jp]);
         }
     }
 
     fclose(matFile);
 
     string rhsFileName = fileBase + ".rhs";
     FILE *rhsFile = fopen(rhsFileName.c_str(),"wb");
     
     for(int i=0; i<nCells; i++)
       fprintf(rhsFile,"%lf\n",-rCell[i]);
 
     fclose(rhsFile);
   }

template<class T>

T FlowModel< T >::Impl::fixedFluxContinuityBC	(	const StorageSite &	faces,
		const Mesh &	mesh,
		MultiFieldMatrix &	matrix,
		MultiField &	xField,
		MultiField &	rField,
		const FlowBC< T > &	bc
	)

inline

Definition at line 12 of file FlowModel_impl.h.

 {
 public:
 
   typedef Array<int> IntArray;
   
   typedef Array<T> TArray;
   typedef Vector<T,3> VectorT3;
   typedef VectorTranspose<T,3> VectorT3T;
 
   typedef Array<VectorT3> VectorT3Array;
   typedef DiagonalTensor<T,3> DiagTensorT3;
 
   typedef CRMatrix<DiagTensorT3,T,VectorT3> VVMatrix;
   typedef typename VVMatrix::DiagArray VVDiagArray;
   
   typedef CRMatrix<T,T,T> PPMatrix;
   typedef typename PPMatrix::DiagArray PPDiagArray;
   typedef typename PPMatrix::PairWiseAssembler PPAssembler;
 
 #ifdef PV_COUPLED
   typedef CRMatrixRect<VectorT3T,VectorT3,T> PVMatrix;
   typedef typename PVMatrix::DiagArray PVDiagArray;
   typedef typename PVMatrix::PairWiseAssembler PVAssembler;
 
   typedef CRMatrixRect<VectorT3,T,VectorT3> VPMatrix;
   typedef typename VPMatrix::DiagArray VPDiagArray;
   typedef typename VPMatrix::PairWiseAssembler VPAssembler;
 
 #endif
   
   typedef Gradient<VectorT3> VGradType;
   typedef Array<Gradient<VectorT3> > VGradArray;
 
   typedef Gradient<T> PGradType;
   typedef Array<PGradType> PGradArray;
 
   typedef FluxJacobianMatrix<T,T> FMatrix;
   
   typedef StressTensor<T> StressTensorT6;
   typedef Array<StressTensorT6> StressTensorArray;
   
   Impl(const GeomFields& geomFields,
        FlowFields& thermalFields,
        const MeshList& meshes):
     _meshes(meshes),
     _geomFields(geomFields),
     _flowFields(thermalFields),
     _velocityGradientModel(_meshes,_flowFields.velocity,
                            _flowFields.velocityGradient,_geomFields),
     _pressureGradientModel(_meshes,_flowFields.pressure,
                            _flowFields.pressureGradient,_geomFields),
     _initialMomentumNorm(),
     _initialContinuityNorm(),
     _niters(0)
   {
     const int numMeshes = _meshes.size();
     for (int n=0; n<numMeshes; n++)
     {
         const Mesh& mesh = *_meshes[n];
 
         FlowVC<T> *vc(new FlowVC<T>());
         vc->vcType = "flow";
         _vcMap[mesh.getID()] = vc;

template<class T>

T FlowModel< T >::Impl::fixedPressureContinuityBC	(	const StorageSite &	faces,
		const Mesh &	mesh,
		MultiFieldMatrix &	matrix,
		const MultiField &	xField,
		MultiField &	rField
	)

inline

Definition at line 59 of file FlowModel_impl.h.

                               :
     _meshes(meshes),
     _geomFields(geomFields),
     _flowFields(thermalFields),
     _velocityGradientModel(_meshes,_flowFields.velocity,
                            _flowFields.velocityGradient,_geomFields),
     _pressureGradientModel(_meshes,_flowFields.pressure,
                            _flowFields.pressureGradient,_geomFields),
     _initialMomentumNorm(),
     _initialContinuityNorm(),
     _niters(0)
   {
     const int numMeshes = _meshes.size();
     for (int n=0; n<numMeshes; n++)
     {
         const Mesh& mesh = *_meshes[n];
 
         FlowVC<T> *vc(new FlowVC<T>());
         vc->vcType = "flow";
         _vcMap[mesh.getID()] = vc;
         foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
         {
             const FaceGroup& fg = *fgPtr;
             if (_bcMap.find(fg.id) == _bcMap.end())
             {
                 FlowBC<T> *bc(new FlowBC<T>());
                 
                 _bcMap[fg.id] = bc;
                 if ((fg.groupType == "wall"))
                 {
                     bc->bcType = "NoSlipWall";
                 }
                 else if ((fg.groupType == "velocity-inlet"))
                 {
                     bc->bcType = "VelocityBoundary";
                 }
                 else if ((fg.groupType == "pressure-inlet") ||
                          (fg.groupType == "pressure-outlet"))
                 {
                     bc->bcType = "PressureBoundary";
                 }
                 else if ((fg.groupType == "symmetry"))
                 {
                     bc->bcType = "Symmetry";
                 }
                 else
                   throw CException("FlowModel: unknown face group type "
                                    + fg.groupType);
             }
         }
         /*
         foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
           {
             const FaceGroup& fg = *fgPtr;
             if ((fg.groupType == "ESinterface"))
               { bc->bcType = "ESInterfaceBC";}
           }
         */
         
 
     }
   }
 
   void init()
   {
     const int numMeshes = _meshes.size();
     for (int n=0; n<numMeshes; n++)
     {
         const Mesh& mesh = *_meshes[n];
 
         const FlowVC<T>& vc = *_vcMap[mesh.getID()];
             
         const StorageSite& cells = mesh.getCells();
         const StorageSite& faces = mesh.getFaces();
 
         shared_ptr<VectorT3Array> vCell(new VectorT3Array(cells.getCountLevel1()));
 
         VectorT3 initialVelocity;
         initialVelocity[0] = _options["initialXVelocity"];
         initialVelocity[1] = _options["initialYVelocity"];
         initialVelocity[2] = _options["initialZVelocity"];
         *vCell = initialVelocity;
 
         _flowFields.velocity.addArray(cells,vCell);
 

template<class T>

void FlowModel< T >::Impl::fixedPressureMomentumBC	(	const StorageSite &	faces,
		const Mesh &	mesh,
		MultiFieldMatrix &	matrix,
		const MultiField &	xField,
		MultiField &	rField
	)

inline

Definition at line 12 of file FlowModel_impl.h.

 {
 public:
 
   typedef Array<int> IntArray;
   
   typedef Array<T> TArray;
   typedef Vector<T,3> VectorT3;
   typedef VectorTranspose<T,3> VectorT3T;
 
   typedef Array<VectorT3> VectorT3Array;
   typedef DiagonalTensor<T,3> DiagTensorT3;
 
   typedef CRMatrix<DiagTensorT3,T,VectorT3> VVMatrix;
   typedef typename VVMatrix::DiagArray VVDiagArray;
   
   typedef CRMatrix<T,T,T> PPMatrix;

template<class T>

FlowBCMap& FlowModel< T >::Impl::getBCMap ( )

inline

Definition at line 344 of file FlowModel_impl.h.

344 {return _bcMap;}

FlowModel::Impl::_bcMap

FlowBCMap _bcMap

Definition: FlowModel_impl.h:2125

template<class T>

VectorT3 FlowModel< T >::Impl::getMomentumDerivativeIntegral ( const Mesh & mesh )

inline

Definition at line 1723 of file FlowModel_impl.h.

References Mesh::getCells(), StorageSite::getSelfCount(), and Vector< double, 3 >::getZero().

   {
     VectorT3 r(VectorT3::getZero());
     const StorageSite& cells = mesh.getCells();
     const int nCells = cells.getSelfCount();
 
     const TArray& density =
           dynamic_cast<const TArray&>(_flowFields.density[cells]);
     const VectorT3Array& v =
           dynamic_cast<const VectorT3Array&>(_flowFields.velocity[cells]);
     const VectorT3Array& vN1 =
       dynamic_cast<const VectorT3Array&>(_flowFields.velocityN1[cells]);
 
     const TArray& volume =
       dynamic_cast<const TArray&>(_geomFields.volume[cells]);
     
     const T deltaT = _options["timeStep"];
     
     if (_flowFields.velocityN2.hasArray(cells))
     {
         // second order
         const VectorT3Array& vN2 =
           dynamic_cast<const VectorT3Array&>(_flowFields.velocityN2[cells]);
         T onePointFive(1.5);
         T two(2.0);
         T pointFive(0.5);
         
         for(int c=0; c<nCells; c++)
         {
             const T rhoVbydT = density[c]*volume[c]/deltaT;
             r += rhoVbydT*(onePointFive*v[c]- two*vN1[c]
                            + pointFive*vN2[c]);
         }
     }
     else
     {
         for(int c=0; c<nCells; c++)
         {
             const T rhoVbydT = density[c]*volume[c]/deltaT;
             r += rhoVbydT*(v[c]- vN1[c]);
         }
     }
     return r;
   }

template<class T>

VectorT3 FlowModel< T >::Impl::getMomentumFluxIntegral	(	const Mesh &	mesh,
		const int	faceGroupId
	)

inline

Definition at line 1700 of file FlowModel_impl.h.

References Mesh::getBoundaryFaceGroups(), StorageSite::getCount(), Vector< double, 3 >::getZero(), FaceGroup::id, and FaceGroup::site.

   {
     VectorT3 r(VectorT3::getZero());
     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 VectorT3Array& momFlux =
               dynamic_cast<const VectorT3Array&>(_flowFields.momentumFlux[faces]);
             for(int f=0; f<nFaces; f++)
               r += momFlux[f];
             found=true;
         }
     }
     if (!found)
       throw CException("getMomentumFluxIntegral: invalid faceGroupID");
     return r;
   }

template<class T>

VectorT3 FlowModel< T >::Impl::getMomentumFluxIntegralonIBFaces ( const Mesh & mesh )

inline

Definition at line 2025 of file FlowModel_impl.h.

References StorageSite::getCount(), Mesh::getFaces(), Mesh::getIBFaceList(), Mesh::getIBFaces(), and Vector< double, 3 >::getZero().

   {
     VectorT3 r(VectorT3::getZero());
     const StorageSite& ibFaces = mesh.getIBFaces();
     const StorageSite& faces = mesh.getFaces();
     const Array<int>& ibFaceIndices = mesh.getIBFaceList();
     const int nibf = ibFaces.getCount();
     const VectorT3Array& momFlux =
               dynamic_cast<const VectorT3Array&>(_flowFields.momentumFlux[faces]);
     for ( int f = 0; f < nibf; f ++){
       const int ibFaceIndex = ibFaceIndices[f];
       r += momFlux[ibFaceIndex];
     }
     if (nibf == 0)
       throw CException("getMomentumFluxIntegralonIBFaces:  no ibFaces found!");
     return r;
   }

template<class T>

FlowModelOptions<T>& FlowModel< T >::Impl::getOptions ( )

inline

Definition at line 346 of file FlowModel_impl.h.

346 {return _options;}

FlowModel::Impl::_options

FlowModelOptions< T > _options

Definition: FlowModel_impl.h:2128

template<class T>

map<string,shared_ptr<ArrayBase> >& FlowModel< T >::Impl::getPersistenceData ( )

inline

Definition at line 457 of file FlowModel_impl.h.

References Model::_persistenceData, and Array< T >::zero().

   {
     _persistenceData.clear();
     
     Array<int>* niterArray = new Array<int>(1);
     (*niterArray)[0] = _niters;
     _persistenceData["niters"]=shared_ptr<ArrayBase>(niterArray);
     
     if (_initialMomentumNorm)
     {
         _persistenceData["initialMomentumNorm"] =
           _initialMomentumNorm->getArrayPtr(_flowFields.velocity);
     }
     else
     {
         Array<Vector<T,3> >* xArray = new Array<Vector<T,3> >(1);
         xArray->zero();
         _persistenceData["initialMomentumNorm"]=shared_ptr<ArrayBase>(xArray);
         
     }
     
     if (_initialContinuityNorm)
     {
         _persistenceData["initialContinuityNorm"] =
           _initialContinuityNorm->getArrayPtr(_flowFields.pressure);
     }
     else
     {
         Array<T>* xArray = new Array<T>(1);
         xArray->zero();
         _persistenceData["initialContinuityNorm"]=shared_ptr<ArrayBase>(xArray);
     }
     return _persistenceData;
   }

template<class T>

VectorT3 FlowModel< T >::Impl::getPressureIntegral	(	const Mesh &	mesh,
		const int	faceGroupId
	)

inline

Definition at line 1638 of file FlowModel_impl.h.

References Mesh::getBoundaryFaceGroups(), StorageSite::getCount(), Vector< double, 3 >::getZero(), FaceGroup::id, and FaceGroup::site.

   {
     VectorT3 r(VectorT3::getZero());
     bool found = false;
     foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
     {
         const FaceGroup& fg = *fgPtr;
         if (fg.id == faceGroupId)
         {
             const StorageSite& faces = fg.site;
             const VectorT3Array& faceArea =
               dynamic_cast<const VectorT3Array&>(_geomFields.area[faces]);
             const int nFaces = faces.getCount();
             const TArray& facePressure = dynamic_cast<const TArray&>(_flowFields.pressure[faces]);
             for(int f=0; f<nFaces; f++)
               r += faceArea[f]*facePressure[f];
 
             found=true;
         }
     }
     if (!found)
       throw CException("getPressureIntegral: invalid faceGroupID");
     return r;
   }

template<class T>

VectorT3 FlowModel< T >::Impl::getPressureIntegralonIBFaces ( const Mesh & mesh )

inline

Definition at line 1663 of file FlowModel_impl.h.

References Mesh::getCells(), StorageSite::getCount(), Vector< T, N >::getData(), Mesh::getFaceCells(), Mesh::getFaces(), Mesh::getIBFaceList(), Mesh::getIBFaces(), Vector< double, 3 >::getZero(), and Mesh::IBTYPE_FLUID.

   {
     VectorT3 r(VectorT3::getZero());
     const StorageSite& ibFaces = mesh.getIBFaces();
     const StorageSite& faces = mesh.getFaces();
     const StorageSite& cells = mesh.getCells();
     const Array<int>& ibFaceIndices = mesh.getIBFaceList();
     const int nibf = ibFaces.getCount();
     const CRConnectivity& faceCells = mesh.getFaceCells(faces);
     const IntArray& ibType = dynamic_cast<const IntArray&>(_geomFields.ibType[cells]);
     const VectorT3Array& faceArea =
               dynamic_cast<const VectorT3Array&>(_geomFields.area[faces]);
     const TArray& facePressure = dynamic_cast<const TArray&>(_flowFields.pressure[faces]);
 
     for ( int f = 0; f < nibf; f ++)
     {
         const int ibFaceIndex = ibFaceIndices[f];
         const int c0 = faceCells(ibFaceIndex,0);
 
         // need to check whether the face points in or out of the
         // fluid cell and adjust sign of area accordingly
         
         if (ibType[c0] == Mesh::IBTYPE_FLUID)
           r += faceArea[ibFaceIndex]*facePressure[ibFaceIndex];
         else
           r -= faceArea[ibFaceIndex]*facePressure[ibFaceIndex];
     }
 
 #if FVM_PARALLEL    
      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, r.getData(), 3, MPI::DOUBLE, MPI::SUM);
 #endif     
 
 
     return r;
   }

template<class T>

boost::shared_ptr<ArrayBase> FlowModel< T >::Impl::getStressTensor	(	const Mesh &	mesh,
		const ArrayBase &	gcellIds
	)

inline

Definition at line 1768 of file FlowModel_impl.h.

References Mesh::getCells(), and Array< T >::getLength().

   {
     typedef Array<StressTensor<T> > StressTensorArray;
     
     const StorageSite& cells = mesh.getCells();
     
     const Array<int>& cellIds = dynamic_cast<const Array<int> &>(gcellIds);
     const int nCells = cellIds.getLength();
 
     _velocityGradientModel.compute();
 
     const VGradArray& vGrad =
       dynamic_cast<const VGradArray&>(_flowFields.velocityGradient[cells]);
 
     const TArray& pCell =
       dynamic_cast<const TArray&>(_flowFields.pressure[cells]);
 
     const TArray& mu = dynamic_cast<const TArray&>(getViscosityField()[cells]);
 
     boost::shared_ptr<StressTensorArray> stressTensorPtr( new StressTensorArray(nCells));
     StressTensorArray& stressTensor = *stressTensorPtr;
 
     for(int n=0; n<nCells; n++)
     {
         const int c = cellIds[n];
         const VGradType& vg = vGrad[c];
         VGradType vgPlusTranspose = vGrad[c];
 
         for(int i=0;i<3;i++)
           for(int j=0;j<3;j++)
             vgPlusTranspose[i][j] += vg[j][i];
         
         stressTensor[n][0] = vgPlusTranspose[0][0]*mu[c] - pCell[c];
         stressTensor[n][1] = vgPlusTranspose[1][1]*mu[c] - pCell[c];
         stressTensor[n][2] = vgPlusTranspose[2][2]*mu[c] - pCell[c];
         stressTensor[n][3] = vgPlusTranspose[0][1]*mu[c];
         stressTensor[n][4] = vgPlusTranspose[1][2]*mu[c];
         stressTensor[n][5] = vgPlusTranspose[2][0]*mu[c];
     }
 
     return stressTensorPtr;
   }

template<class T>

void FlowModel< T >::Impl::getTraction ( const Mesh & mesh )

inline

Definition at line 1884 of file FlowModel_impl.h.

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

   {
     const StorageSite& cells = mesh.getCells();
 
     const int nCells = cells.getSelfCount();
 
     shared_ptr<VectorT3Array> tractionXPtr(new VectorT3Array(nCells));
     tractionXPtr->zero();
     _flowFields.tractionX.addArray(cells,tractionXPtr);
     VectorT3Array& tractionX = *tractionXPtr;
 
     shared_ptr<VectorT3Array> tractionYPtr(new VectorT3Array(nCells));
     tractionYPtr->zero();
     _flowFields.tractionY.addArray(cells,tractionYPtr);
     VectorT3Array& tractionY = *tractionYPtr;
 
     shared_ptr<VectorT3Array> tractionZPtr(new VectorT3Array(nCells));
     tractionZPtr->zero();
     _flowFields.tractionZ.addArray(cells,tractionZPtr);
     VectorT3Array& tractionZ = *tractionZPtr;
 
     _velocityGradientModel.compute();
 
     const VGradArray& vGrad =
       dynamic_cast<const VGradArray&>(_flowFields.velocityGradient[cells]);
 
     const TArray& pCell =
       dynamic_cast<const TArray&>(_flowFields.pressure[cells]);
 
     const TArray& mu = dynamic_cast<const TArray&>(getViscosityField()[cells]);
       
     for(int n=0; n<nCells; n++)
     {
         const VGradType& vg = vGrad[n];
         VGradType vgPlusTranspose = vGrad[n];
 
         for(int i=0;i<3;i++)
           for(int j=0;j<3;j++)
             vgPlusTranspose[i][j] += vg[j][i];
           
         tractionX[n][0] = vgPlusTranspose[0][0]*mu[n] - pCell[n];
         tractionX[n][1] = vgPlusTranspose[0][1]*mu[n];
         tractionX[n][2] = vgPlusTranspose[0][2]*mu[n];
 
         tractionY[n][0] = vgPlusTranspose[1][0]*mu[n];
         tractionY[n][1] = vgPlusTranspose[1][1]*mu[n] - pCell[n];
         tractionY[n][2] = vgPlusTranspose[1][2]*mu[n];
 
         tractionZ[n][0] = vgPlusTranspose[2][0]*mu[n];
         tractionZ[n][1] = vgPlusTranspose[2][1]*mu[n];
         tractionZ[n][2] = vgPlusTranspose[2][2]*mu[n] - pCell[n];
     }
   }

template<class T>

FlowVCMap& FlowModel< T >::Impl::getVCMap ( )

inline

Definition at line 345 of file FlowModel_impl.h.

345 {return _vcMap;}

FlowModel::Impl::_vcMap

FlowVCMap _vcMap

Definition: FlowModel_impl.h:2126

template<class T>

const Field& FlowModel< T >::Impl::getViscosityField ( ) const

inline

Definition at line 370 of file FlowModel_impl.h.

   {
     if (_options.turbulent)
       return _flowFields.totalviscosity;
     else
       return _flowFields.viscosity;
   }

template<class T>

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

inline

Definition at line 148 of file FlowModel_impl.h.

References Model::_meshes, Mesh::getCells(), StorageSite::getCount(), StorageSite::getCountLevel1(), Mesh::getFaces(), and Mesh::getID().

   {
     const int numMeshes = _meshes.size();
     for (int n=0; n<numMeshes; n++)
     {
         const Mesh& mesh = *_meshes[n];
 
         const FlowVC<T>& vc = *_vcMap[mesh.getID()];
             
         const StorageSite& cells = mesh.getCells();
         const StorageSite& faces = mesh.getFaces();
 
         shared_ptr<VectorT3Array> vCell(new VectorT3Array(cells.getCountLevel1()));
 
         VectorT3 initialVelocity;
         initialVelocity[0] = _options["initialXVelocity"];
         initialVelocity[1] = _options["initialYVelocity"];
         initialVelocity[2] = _options["initialZVelocity"];
         *vCell = initialVelocity;
 
         _flowFields.velocity.addArray(cells,vCell);
 
         shared_ptr<VectorT3Array> uparallelCell(new VectorT3Array(cells.getCount()));
         uparallelCell ->zero();
         _flowFields.uparallel.addArray(cells,uparallelCell);
 
       shared_ptr<VectorT3Array> tauCell(new VectorT3Array(faces.getCount()));
         tauCell ->zero();
         _flowFields.tau.addArray(faces,tauCell);
 
      shared_ptr<VectorT3Array> TauWallCell(new VectorT3Array(faces.getCount()));
         TauWallCell ->zero();
         _flowFields.tauwall.addArray(faces,TauWallCell);
 
 
 
 
   
        
 
         if (_options.transient)
         {
             _flowFields.velocityN1.addArray(cells,
                                             dynamic_pointer_cast<ArrayBase>(vCell->newCopy()));
             if (_options.timeDiscretizationOrder > 1)
               _flowFields.velocityN2.addArray(cells,
                                               dynamic_pointer_cast<ArrayBase>(vCell->newCopy()));
 
         }
         
         shared_ptr<TArray> pCell(new TArray(cells.getCountLevel1()));
         shared_ptr<TArray> pFace(new TArray(faces.getCountLevel1()));
         *pCell = _options["initialPressure"];
         *pFace = _options["initialPressure"];
         _flowFields.pressure.addArray(cells,pCell);
         _flowFields.pressure.addArray(faces,pFace);
 
 
         shared_ptr<TArray> rhoCell(new TArray(cells.getCountLevel1()));
         *rhoCell = vc["density"];
         _flowFields.density.addArray(cells,rhoCell);
 
         shared_ptr<TArray> muCell(new TArray(cells.getCountLevel1()));
         *muCell = vc["viscosity"];
         _flowFields.viscosity.addArray(cells,muCell);
 
         shared_ptr<PGradArray> gradp(new PGradArray(cells.getCountLevel1()));
         gradp->zero();
         _flowFields.pressureGradient.addArray(cells,gradp);
 
         shared_ptr<TArray> ci(new TArray(cells.getCountLevel1()));
         ci->zero();
         _flowFields.continuityResidual.addArray(cells,ci);
 
         // compute default value of mass flux
         
         const VectorT3Array& faceArea =
           dynamic_cast<const VectorT3Array&>(_geomFields.area[faces]);
         const VectorT3Array& V = *vCell;
         const TArray& density = *rhoCell;
         
         const CRConnectivity& faceCells = mesh.getAllFaceCells();
 
         const int nFaces = faces.getCount();
 
         shared_ptr<TArray> mfPtr(new TArray(faces.getCount()));
         TArray& mf = *mfPtr;
 
         for(int f=0; f<nFaces; f++)
         {
             const int c0 = faceCells(f,0);
             const int c1 = faceCells(f,1);
             
             mf[f] = 0.5*(density[c0]*dot(V[c0],faceArea[f]) +
                          density[c1]*dot(V[c1],faceArea[f]));
         }
         _flowFields.massFlux.addArray(faces,mfPtr);
 
         // store momentum flux at interfaces
         foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
         {
             const FaceGroup& fg = *fgPtr;
             const StorageSite& faces = fg.site;
             shared_ptr<VectorT3Array> momFlux(new VectorT3Array(faces.getCount()));
             momFlux->zero();
             _flowFields.momentumFlux.addArray(faces,momFlux);
             
             if(fg.groupType == "ESinterface"){
               cout << "interface init" <<endl;
               const StorageSite& Intfaces = fg.site; 
               
               shared_ptr<VectorT3Array> InterfaceVelFace(new VectorT3Array(Intfaces.getCount()));
               InterfaceVelFace ->zero();
               _flowFields.InterfaceVelocity.addArray(Intfaces,InterfaceVelFace);
               
               shared_ptr<StressTensorArray> InterfaceStressFace(new StressTensorArray(Intfaces.getCount()));
               InterfaceStressFace ->zero();
               _flowFields.InterfaceStress.addArray(Intfaces,InterfaceStressFace);
               
               shared_ptr<TArray> InterfacePressFace(new TArray(Intfaces.getCount()));
               *InterfacePressFace = _options["initialPressure"];
               _flowFields.InterfacePressure.addArray(Intfaces,InterfacePressFace);
               
               shared_ptr<TArray> InterfaceDensityFace(new TArray(Intfaces.getCount()));
               *InterfaceDensityFace =vc["density"];
               _flowFields.InterfaceDensity.addArray(Intfaces,InterfaceDensityFace);
               
               
             }
             
             
         }
         
         foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
         {
             const FaceGroup& fg = *fgPtr;
             const StorageSite& faces = fg.site;
 
             const FlowBC<T>& bc = *_bcMap[fg.id];
 
             FloatValEvaluator<VectorT3>
               bVelocity(bc.getVal("specifiedXVelocity"),
                         bc.getVal("specifiedYVelocity"),
                         bc.getVal("specifiedZVelocity"),
                         faces);
 
             const int nFaces = faces.getCount();
             const CRConnectivity& faceCells = mesh.getFaceCells(faces);
 
             if ((bc.bcType == "NoSlipWall") ||
                 (bc.bcType == "SlipJump") ||
                 (bc.bcType == "VelocityBoundary"))
             {
                 // arrays for this face group
                 TArray& massFlux = dynamic_cast<TArray&>(_flowFields.massFlux[faces]);
                 
                 const VectorT3Array& faceArea =
                   dynamic_cast<const VectorT3Array&>(_geomFields.area[faces]);
                 
                 for(int f=0; f<nFaces; f++)
                 {
                     const int c0 = faceCells(f,0);
                     massFlux[f] = density[c0]*dot(bVelocity[f],faceArea[f]);
                  }
             }
             else if (bc.bcType == "SpecifiedPressure")
             {
                 T bp = bc["specifiedPressure"];
                 TArray& facePressure = dynamic_cast<TArray&>(_flowFields.pressure[faces]);
                 TArray& cellPressure = dynamic_cast<TArray&>(_flowFields.pressure[cells]);
                 for(int f=0; f<nFaces; f++)
                 {
                     const int c1 = faceCells(f,1);
                     facePressure[f] = cellPressure[c1] = bp;
                 }
             }
             else if ((bc.bcType == "Symmetry"))
             {
                 TArray& massFlux = dynamic_cast<TArray&>(_flowFields.massFlux[faces]);
                 massFlux.zero();
             }
             shared_ptr<VectorT3Array> momFlux(new VectorT3Array(faces.getCount()));
             momFlux->zero();
             _flowFields.momentumFlux.addArray(faces,momFlux);
            
             
             
         }
     }
 
     computeContinuityResidual();
     _niters  =0;
     _initialMomentumNorm = MFRPtr();
     _initialContinuityNorm = MFRPtr();
   }

template<class T>

void FlowModel< T >::Impl::initContinuityLinearization ( LinearSystem & ls )

inline

Definition at line 1341 of file FlowModel_impl.h.

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

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

template<class T>

void FlowModel< T >::Impl::initMomentumLinearization ( LinearSystem & ls )

inline

Definition at line 522 of file FlowModel_impl.h.

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

   {
     const int numMeshes = _meshes.size();
     for (int n=0; n<numMeshes; n++)
     {
         const Mesh& mesh = *_meshes[n];
 
         const StorageSite& cells = mesh.getCells();
         MultiField::ArrayIndex vIndex(&_flowFields.velocity,&cells);
 
         ls.getX().addArray(vIndex,_flowFields.velocity.getArrayPtr(cells));
 
         const CRConnectivity& cellCells = mesh.getCellCells();
         
         shared_ptr<Matrix> m(new CRMatrix<DiagTensorT3,T,VectorT3>(cellCells));
 
         ls.getMatrix().addMatrix(vIndex,vIndex,m);
 
         foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
         {
             const FaceGroup& fg = *fgPtr;
             const StorageSite& faces = fg.site;
 
             MultiField::ArrayIndex fIndex(&_flowFields.momentumFlux,&faces);
             ls.getX().addArray(fIndex,_flowFields.momentumFlux.getArrayPtr(faces));
 
             const CRConnectivity& faceCells = mesh.getFaceCells(faces);
 
             shared_ptr<Matrix> mft(new FluxJacobianMatrix<DiagTensorT3,VectorT3>(faceCells));
             ls.getMatrix().addMatrix(fIndex,vIndex,mft);
 
             shared_ptr<Matrix> mff(new DiagonalMatrix<DiagTensorT3,VectorT3>(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(&_flowFields.momentumFlux,&faces);
             ls.getX().addArray(fIndex,_flowFields.momentumFlux.getArrayPtr(faces));
 
             const CRConnectivity& faceCells = mesh.getFaceCells(faces);
 
             shared_ptr<Matrix> mft(new FluxJacobianMatrix<DiagTensorT3,VectorT3>(faceCells));
             ls.getMatrix().addMatrix(fIndex,vIndex,mft);
 
             shared_ptr<Matrix> mff(new DiagonalMatrix<DiagTensorT3,VectorT3>(faces.getCount()));
             ls.getMatrix().addMatrix(fIndex,fIndex,mff);
         }
     }
   }

template<class T>

void FlowModel< T >::Impl::interfaceContinuityBC	(	const Mesh &	mesh,
		const StorageSite &	faces,
		MultiFieldMatrix &	matrix,
		MultiField &	rField
	)

inline

Definition at line 774 of file FlowModel_impl.h.

References Mesh::getCells(), StorageSite::getCount(), Mesh::getFaceCells(), and StorageSite::getSelfCount().

   {
     const StorageSite& cells = mesh.getCells();
     MultiField::ArrayIndex pIndex(&_flowFields.pressure,&cells);
     const CRConnectivity& faceCells = mesh.getFaceCells(faces);
 
     TArray& rCell = dynamic_cast<TArray&>(rField[pIndex]);
     
     const int nFaces = faces.getCount();
 
     for(int f=0; f<nFaces; f++)
     {
         int c0 = faceCells(f,0);
         int c1 = faceCells(f,1);
         // either c0 or c1 could be the exterior cell
         if (c1 >= cells.getSelfCount())
         {
             rCell[c1] = 0;
         }
         else
         {
             rCell[c0] = 0;
         }
     }
   }

template<class T>

void FlowModel< T >::Impl::linearizeContinuity ( LinearSystem & ls )

inline

Definition at line 998 of file FlowModel_impl.h.

References Model::_meshes, FlowBC< T >::bcType, discretizeMassFluxInterior(), dot(), fixedFluxContinuityBC(), fixedPressureContinuityBC(), LinearSystem::getB(), Mesh::getBoundaryFaceGroups(), Mesh::getCells(), StorageSite::getCount(), Mesh::getFaceCells(), Mesh::getInterfaceGroups(), Mesh::getInteriorFaceGroup(), LinearSystem::getMatrix(), MultiFieldMatrix::getMatrix(), StorageSite::getSelfCount(), LinearSystem::getX(), Mesh::IBTYPE_FLUID, FaceGroup::id, FaceGroup::site, and DiagonalMatrix< Diag, X >::unitize().

   {
     MultiFieldMatrix& matrix = ls.getMatrix();
     MultiField& x = ls.getX();
     MultiField& b = ls.getB();
 
     
     this->_useReferencePressure = true;
     const int numMeshes = _meshes.size();
 
     T netFlux(0.);
     
 
     _flowFields.pressureGradient.syncLocal();
 
     
     for (int n=0; n<numMeshes; n++)
     {
         const Mesh& mesh = *_meshes[n];
         const StorageSite& cells = mesh.getCells();
         const StorageSite& iFaces = mesh.getInteriorFaceGroup().site;
         
         MultiField::ArrayIndex pIndex(&_flowFields.pressure,&cells);
         
         // interior
         
         netFlux += discretizeMassFluxInterior(mesh,iFaces,matrix,x,b);
         
         // interfaces
         foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
         {
             const FaceGroup& fg = *fgPtr;
             const StorageSite& faces = fg.site;
             netFlux += discretizeMassFluxInterior(mesh,faces,matrix,x,b);
 
             // this just sets the residual for the ghost cell to zero
             interfaceContinuityBC(mesh,faces,matrix,b);
         }
         
         // boundaries
         foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
         {
             const FaceGroup& fg = *fgPtr;
             const StorageSite& faces = fg.site;
 
             const FlowBC<T>& bc = *_bcMap[fg.id];
             
             if ((bc.bcType == "NoSlipWall") ||
                 (bc.bcType == "SlipJump") ||
                 (bc.bcType == "Symmetry") ||
                 (bc.bcType == "VelocityBoundary"))
             {
                 
                 netFlux += fixedFluxContinuityBC(faces,mesh,matrix,x,b,bc);
             }
             else if (bc.bcType == "PressureBoundary")
             {
                 netFlux += fixedPressureContinuityBC(faces,mesh,matrix,x,b);
                 this->_useReferencePressure=false;
             }   
             else                
               throw CException(bc.bcType + " not implemented for FlowModel");
             
             MultiField::ArrayIndex mfIndex(&_flowFields.massFlux,&faces);
             typedef DiagonalMatrix<T,T> FFMatrix;
             FFMatrix& dFluxdFlux =
               dynamic_cast<FFMatrix&>(matrix.getMatrix(mfIndex,mfIndex));
             dFluxdFlux.unitize();
         }
 
         // add additional imbalance to netFlux due to transient term
         // when we have deforming mesh
         if (_geomFields.volumeN1.hasArray(cells))
         {
 
             // this block computes the net sum of d(rho)/dt over all the cells
             const TArray& density =
               dynamic_cast<const TArray&>(_flowFields.density[cells]);  
             const TArray& cellVolume =
               dynamic_cast<const TArray&>(_geomFields.volume[cells]);
             const int nCells = cells.getSelfCount();
             T _dT(_options["timeStep"]);
             T onePointFive(1.5);
             T two(2.0);
             T pointFive(0.5);
 
             const TArray& cellVolumeN1 = 
               dynamic_cast<const TArray&>(_geomFields.volumeN1[cells]);
 
             if (_geomFields.volumeN2.hasArray(cells))
             {
                 // second order time discretization
 
                 const TArray& cellVolumeN2 = 
                   dynamic_cast<const TArray&>(_geomFields.volumeN2[cells]);
                 for(int c=0; c<nCells; c++)
                 {
                     const T rhobydT = density[c]/_dT;
                     const T term1 = onePointFive*cellVolume[c];
                     const T term2 = two*cellVolumeN1[c];
                     const T term3 = pointFive*cellVolumeN2[c];
                     netFlux -= rhobydT*(term1 - term2 + term3);
                 }
             }
             else
             {
                 // first order time discretization
                 for(int c=0; c<nCells; c++)
                 {       
                     const T rhobydT = density[c]/_dT;
                     netFlux -= rhobydT*(cellVolume[c] - cellVolumeN1[c]);
                 }
             }
 
             // compute the net integral of the grid flux term 
             foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
             {
                 const FaceGroup& fg = *fgPtr;
                 const StorageSite& faces = fg.site;
                 const int nFaces = faces.getCount();
                 const CRConnectivity& faceCells = mesh.getFaceCells(faces);
                 const VectorT3Array& faceArea =
                   dynamic_cast<const VectorT3Array&>(_geomFields.area[faces]);
                 const VectorT3Array& faceVel =
                   dynamic_cast<const VectorT3Array&>(_geomFields.faceVel[faces]);
                 for(int f=0; f<nFaces; f++)
                 {
                     const int c0 = faceCells(f,0);
                     netFlux += density[c0]*dot(faceVel[f],faceArea[f]);
                 }
             }
         }
     }
   
  
       //sum netflux globalally MPI::
 #ifdef FVM_PARALLEL
       int count = 1;
       MPI::COMM_WORLD.Allreduce( MPI::IN_PLACE, &netFlux, count, MPI::DOUBLE, MPI::SUM);
 
       int useReferencePressureInt = int(this->_useReferencePressure );
       MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &useReferencePressureInt, count, MPI::INT, MPI::PROD);
       this->_useReferencePressure = bool(useReferencePressureInt);
 #endif
         //cout << "net boundary flux = " << netFlux << endl;
     
     if (this->_useReferencePressure)
     {
         // we have all fixed mass flux boundary problem. In case the
         // net mass flux sum over all the boundaries does not add up
         // to zero we distribute that over the cells as a volumetric
         // as a matching source/sink so that the continuity equation
         // can converge to a zero residual. In case of immersed
         // boundary problems we need to use only the fluid cells since
         // all the other cells will always have zero corrections.
         
         T volumeSum(0.);
         
         for (int n=0; n<numMeshes; n++)
         {
             const Mesh& mesh = *_meshes[n];
             const StorageSite& cells = mesh.getCells();
             const IntArray& ibType = dynamic_cast<const IntArray&>(_geomFields.ibType[cells]);
             const TArray& cellVolume = dynamic_cast<const TArray&>(_geomFields.volume[cells]);
             for(int c=0; c<cells.getSelfCount(); c++)
               if (ibType[c] == Mesh::IBTYPE_FLUID)
                 volumeSum += cellVolume[c];
         }
 
 #ifdef FVM_PARALLEL
         MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &volumeSum, 1, MPI::DOUBLE, MPI::SUM);
 #endif  
 
         netFlux /= volumeSum;
 
         for (int n=0; n<numMeshes; n++)
         {
             const Mesh& mesh = *_meshes[n];
             const StorageSite& cells = mesh.getCells();
             const IntArray& ibType = dynamic_cast<const IntArray&>(_geomFields.ibType[cells]);
             const TArray& cellVolume = dynamic_cast<const TArray&>(_geomFields.volume[cells]);
             
             MultiField::ArrayIndex pIndex(&_flowFields.pressure,&cells);
             TArray& rCell = dynamic_cast<TArray&>(b[pIndex]);
             
             for(int c=0; c<cells.getSelfCount(); c++)
             {
               if (ibType[c] == Mesh::IBTYPE_FLUID)
                 rCell[c] += netFlux*cellVolume[c];
             }
         }
 
         // when all the mass fluxes are specified the continuity
         // equation also has an extra degree of freedom. This sets the
         // first cell to have a zero correction to account for this.
         
 
       setDirichlet(matrix,b);
 
     }
   }

template<class T>

void FlowModel< T >::Impl::linearizeMomentum ( LinearSystem & ls )

inline

Definition at line 576 of file FlowModel_impl.h.

References Model::_meshes, BaseGenericBCS< X, Diag, OffDiag >::applyInterfaceBC(), FlowBC< T >::bcType, fixedPressureMomentumBC(), LinearSystem::getB(), Mesh::getBoundaryFaceGroups(), StorageSite::getCount(), Mesh::getInterfaceGroups(), LinearSystem::getMatrix(), FloatVarDict< T >::getVal(), LinearSystem::getX(), FaceGroup::id, Linearizer::linearize(), FaceGroup::site, and slipJumpMomentumBC().

   {
     _velocityGradientModel.compute();
     
     DiscrList discretizations;
     shared_ptr<Discretization>
       dd(new DiffusionDiscretization<VectorT3,DiagTensorT3,T>
          (_meshes,_geomFields,
           _flowFields.velocity,
            getViscosityField(),
           _flowFields.velocityGradient));
 
     shared_ptr<Discretization>
       cd(new ConvectionDiscretization<VectorT3,DiagTensorT3,T>
          (_meshes,_geomFields,
           _flowFields.velocity,
           _flowFields.massFlux,
           _flowFields.continuityResidual,
           _flowFields.velocityGradient));
 
     shared_ptr<Discretization>
       pd(new MomentumPressureGradientDiscretization<T>
          (_meshes,_geomFields,
           _flowFields,
           _pressureGradientModel));
 
     discretizations.push_back(dd);
     discretizations.push_back(cd);
     discretizations.push_back(pd);
 
     if (_options.transient)
     {
         shared_ptr<Discretization>
           td(new TimeDerivativeDiscretization<VectorT3,DiagTensorT3,T>
              (_meshes,_geomFields,
               _flowFields.velocity,
               _flowFields.velocityN1,
               _flowFields.velocityN2,
               _flowFields.density,
               _options["timeStep"]));
         
         discretizations.push_back(td);
     }
 
 
     shared_ptr<Discretization>
       ibm(new GenericIBDiscretization<VectorT3,DiagTensorT3,T>
              (_meshes,_geomFields,_flowFields.velocity));
       
     discretizations.push_back(ibm);
     Linearizer linearizer;
 
     linearizer.linearize(discretizations,_meshes,ls.getMatrix(),
                          ls.getX(), ls.getB());
 
     const int numMeshes = _meshes.size();
     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;
             const int nFaces = faces.getCount();
 
             const FlowBC<T>& bc = *_bcMap[fg.id];
             
 
             GenericBCS<VectorT3,DiagTensorT3,T> gbc(faces,mesh,
                                                     _geomFields,
                                                     _flowFields.velocity,
                                                     _flowFields.momentumFlux,
                                                     ls.getMatrix(),
                                                     ls.getX(),
                                                     ls.getB());
 
             const TArray& massFlux = dynamic_cast<const TArray&>(_flowFields.massFlux[faces]);
             //            const CRConnectivity& faceCells = mesh.getFaceCells(faces);
     
             FloatValEvaluator<VectorT3>
               bVelocity(bc.getVal("specifiedXVelocity"),
                         bc.getVal("specifiedYVelocity"),
                         bc.getVal("specifiedZVelocity"),
                         faces);
 
             if (bc.bcType == "NoSlipWall")
             {
                 gbc.applyDirichletBC(bVelocity);
             }
             else if (bc.bcType == "SlipJump")
             {
                 slipJumpMomentumBC(faces,mesh,
                                    gbc,
                                    bc["accomodationCoefficient"],
                                    bVelocity);
             }
             else if ((bc.bcType == "VelocityBoundary") ||
                      (bc.bcType == "PressureBoundary"))
             {
                 for(int f=0; f<nFaces; f++)
                 {
                     if (massFlux[f] > 0.)
                     {
                         gbc.applyExtrapolationBC(f);
                     }
                     else
                     {
                         gbc.applyDirichletBC(f,bVelocity[f]);
                     }
                 }
                 if (bc.bcType == "PressureBoundary")
                 {
                     fixedPressureMomentumBC(faces,mesh,
                                             ls.getMatrix(), ls.getX(), ls.getB());
                 }
             }
             else if ((bc.bcType == "Symmetry"))
             {
                 gbc.applySymmetryBC();
             }
            
             else
               throw CException(bc.bcType + " not implemented for FlowModel");
         }
 
         foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
         {
             const FaceGroup& fg = *fgPtr;
             const StorageSite& faces = fg.site;
             GenericBCS<VectorT3,DiagTensorT3,T> gbc(faces,mesh,
                                                     _geomFields,
                                                     _flowFields.velocity,
                                                     _flowFields.momentumFlux,
                                                     ls.getMatrix(), ls.getX(), ls.getB());
 
             gbc.applyInterfaceBC();
         }
 
     }
     DiscrList discretizations2;
     shared_ptr<Discretization>
       ud(new Underrelaxer<VectorT3,DiagTensorT3,T>
          (_meshes,_flowFields.velocity,
           _options["momentumURF"]));
 
     discretizations2.push_back(ud);
 
     linearizer.linearize(discretizations2,_meshes,ls.getMatrix(),
                          ls.getX(), ls.getB());
 
   }

template<class T>

void FlowModel< T >::Impl::postContinuitySolve ( LinearSystem & ls )

inline

Definition at line 1263 of file FlowModel_impl.h.

References Model::_meshes, FlowBC< T >::bcType, correctMassFluxInterior(), correctVelocityInterior(), Mesh::getBoundaryFaceGroups(), Mesh::getCells(), LinearSystem::getDelta(), Mesh::getInterfaceGroups(), Mesh::getInteriorFaceGroup(), LinearSystem::getMatrix(), MultiFieldMatrix::hasMatrix(), FaceGroup::id, pressureBoundaryPostContinuitySolve(), FaceGroup::site, updateFacePressureBoundary(), and updateFacePressureInterior().

   {
     MultiFieldMatrix& matrix = ls.getMatrix();
     MultiField& ppField = ls.getDelta();
 
     setReferencePP(ppField);
     
     const int numMeshes = _meshes.size();
     for (int n=0; n<numMeshes; n++)
     {
         const Mesh& mesh = *_meshes[n];
 
         const StorageSite& cells = mesh.getCells();
         const StorageSite& iFaces = mesh.getInteriorFaceGroup().site;
         
         MultiField::ArrayIndex pIndex(&_flowFields.pressure,&cells);
         
         MultiField::ArrayIndex vIndex(&_flowFields.velocity,&cells);
         const bool coupled =  matrix.hasMatrix(pIndex,vIndex);
         
         correctPressure(mesh,ppField);
         
         // interior
         
         correctMassFluxInterior(mesh,iFaces,matrix,ppField);
         
         if (coupled)
           correctVelocityExplicit(mesh,ppField);
         else if (_options.correctVelocity)
           correctVelocityInterior(mesh,iFaces,ppField);
         
         updateFacePressureInterior(mesh,iFaces);
 
         // interfaces
         foreach(const FaceGroupPtr fgPtr, mesh.getInterfaceGroups())
         {
             const FaceGroup& fg = *fgPtr;
             const StorageSite& faces = fg.site;
             correctMassFluxInterior(mesh,faces,matrix,ppField);
             
             if (coupled)
               correctVelocityExplicit(mesh,ppField);
             else if (_options.correctVelocity)
               correctVelocityInterior(mesh,faces,ppField);
             
             updateFacePressureInterior(mesh,faces);
         }
         
         // boundary
         foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
         {
             const FaceGroup& fg = *fgPtr;
             const StorageSite& faces = fg.site;
             const FlowBC<T>& bc = *_bcMap[fg.id];
               
             correctMassFluxBoundary(faces,ppField);
             
             if (!coupled && _options.correctVelocity)
               correctVelocityBoundary(mesh,faces,ppField);
 
             if (bc.bcType == "PressureBoundary")
             {
                 T bp = bc["specifiedPressure"];
                 pressureBoundaryPostContinuitySolve(faces,mesh,bp);
             }
             
             updateFacePressureBoundary(mesh,faces);
         }
 
     }
 
     _flowFields.velocity.syncLocal();
     _flowFields.pressure.syncLocal();
 
     computeContinuityResidual();
     
   }

template<class T>

void FlowModel< T >::Impl::pressureBoundaryPostContinuitySolve	(	const StorageSite &	faces,
		const Mesh &	mesh,
		const T	bp
	)

inline

Definition at line 171 of file FlowModel_impl.h.

         {
             _flowFields.velocityN1.addArray(cells,
                                             dynamic_pointer_cast<ArrayBase>(vCell->newCopy()));
             if (_options.timeDiscretizationOrder > 1)
               _flowFields.velocityN2.addArray(cells,
                                               dynamic_pointer_cast<ArrayBase>(vCell->newCopy()));
 
         }
         
         shared_ptr<TArray> pCell(new TArray(cells.getCountLevel1()));
         shared_ptr<TArray> pFace(new TArray(faces.getCountLevel1()));
         *pCell = _options["initialPressure"];
         *pFace = _options["initialPressure"];
         _flowFields.pressure.addArray(cells,pCell);
         _flowFields.pressure.addArray(faces,pFace);
 
 
         shared_ptr<TArray> rhoCell(new TArray(cells.getCountLevel1()));
         *rhoCell = vc["density"];
         _flowFields.density.addArray(cells,rhoCell);
 
         shared_ptr<TArray> muCell(new TArray(cells.getCountLevel1()));
         *muCell = vc["viscosity"];
         _flowFields.viscosity.addArray(cells,muCell);
 
         shared_ptr<PGradArray> gradp(new PGradArray(cells.getCountLevel1()));

template<class T>

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

inline

Definition at line 1625 of file FlowModel_impl.h.

   {
     foreach(typename FlowBCMap::value_type& pos, _bcMap)
     {
         cout << "Face Group " << pos.first << ":" << endl;
         cout << "    bc type " << pos.second->bcType << endl;
         foreach(typename FlowBC<T>::value_type& vp, *pos.second)
         {
             cout << "   " << vp.first << " "  << vp.second.constant <<  endl;
         }
     }
   }

template<class T>

void FlowModel< T >::Impl::printMassFluxIntegrals ( )

inline

Definition at line 2093 of file FlowModel_impl.h.

References Model::_meshes, Mesh::getBoundaryFaceGroups(), StorageSite::getCount(), Mesh::getID(), FaceGroup::id, and FaceGroup::site.

   {
     const int numMeshes = _meshes.size();
     for (int n=0; n<numMeshes; n++)
     {
         const Mesh& mesh = *_meshes[n];
         foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
         {
             const FaceGroup& fg = *fgPtr;
             
             T r(0.);
             
             const StorageSite& faces = fg.site;
             const int nFaces = faces.getCount();
             const TArray& massFlux = dynamic_cast<const TArray&>(_flowFields.massFlux[faces]);
             for(int f=0; f<nFaces; f++)
               r += massFlux[f];
 
             cout << "Mesh " << mesh.getID() << " faceGroup " << fg.id << " : " << r <<  endl;
         }
     }
   }

template<class T>

void FlowModel< T >::Impl::printMomentumFluxIntegrals ( )

inline

Definition at line 2069 of file FlowModel_impl.h.

References Model::_meshes, Mesh::getBoundaryFaceGroups(), StorageSite::getCount(), Mesh::getID(), Vector< double, 3 >::getZero(), FaceGroup::id, and FaceGroup::site.

   {
     const int numMeshes = _meshes.size();
     for (int n=0; n<numMeshes; n++)
     {
         const Mesh& mesh = *_meshes[n];
         foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
         {
             const FaceGroup& fg = *fgPtr;
             
             VectorT3 r(VectorT3::getZero());
             
             const StorageSite& faces = fg.site;
             const int nFaces = faces.getCount();
             const VectorT3Array& momFlux =
               dynamic_cast<const VectorT3Array&>(_flowFields.momentumFlux[faces]);
             for(int f=0; f<nFaces; f++)
               r += momFlux[f];
 
             cout << "Mesh " << mesh.getID() << " faceGroup " << fg.id << " : " << r <<  endl;
         }
     }
   }

template<class T>

void FlowModel< T >::Impl::printPressureIntegrals ( )

inline

Definition at line 2044 of file FlowModel_impl.h.

References Model::_meshes, Mesh::getBoundaryFaceGroups(), StorageSite::getCount(), Mesh::getID(), Vector< double, 3 >::getZero(), FaceGroup::id, and FaceGroup::site.

   {
     const int numMeshes = _meshes.size();
     for (int n=0; n<numMeshes; n++)
     {
         const Mesh& mesh = *_meshes[n];
         foreach(const FaceGroupPtr fgPtr, mesh.getBoundaryFaceGroups())
         {
             const FaceGroup& fg = *fgPtr;
             
             VectorT3 r(VectorT3::getZero());
             
             const StorageSite& faces = fg.site;
             const VectorT3Array& faceArea =
               dynamic_cast<const VectorT3Array&>(_geomFields.area[faces]);
             const int nFaces = faces.getCount();
             const TArray& facePressure = dynamic_cast<const TArray&>(_flowFields.pressure[faces]);
             for(int f=0; f<nFaces; f++)
               r += faceArea[f]*facePressure[f];
 
             cout << "Mesh " << mesh.getID() << " faceGroup " << fg.id << " : " << r <<  endl;
         }
     }
   }

template<class T>

void FlowModel< T >::Impl::restart ( )

inline

Definition at line 492 of file FlowModel_impl.h.

References Model::_persistenceData.

   {
     if (_persistenceData.find("niters") != _persistenceData.end())
     {
         shared_ptr<ArrayBase> rp = _persistenceData["niters"];
         ArrayBase& r = *rp;
         Array<int>& niterArray = dynamic_cast<Array<int>& >(r);
         _niters = niterArray[0];
     }
 
     if (_persistenceData.find("initialMomentumNorm") != _persistenceData.end())
     {
         shared_ptr<ArrayBase>  r = _persistenceData["initialMomentumNorm"];
         _initialMomentumNorm = MFRPtr(new MultiFieldReduction());
         _initialMomentumNorm->addArray(_flowFields.velocity,r);
     }
 
     if (_persistenceData.find("initialContinuityNorm") != _persistenceData.end())
     {
         shared_ptr<ArrayBase>  r = _persistenceData["initialContinuityNorm"];
         _initialContinuityNorm = MFRPtr(new MultiFieldReduction());
         _initialContinuityNorm->addArray(_flowFields.pressure,r);
     }
     
     computeContinuityResidual();
     
   }

template<class T>

void FlowModel< T >::Impl::setDirichlet	(	MultiFieldMatrix &	mfmatrix,
		MultiField &	rField
	)

inline

Definition at line 906 of file FlowModel_impl.h.

References Model::_meshes, Mesh::getCells(), CRMatrix< T_Diag, T_OffDiag, X >::getConnectivity(), CRMatrix< T_Diag, T_OffDiag, X >::getDiag(), Mesh::getGlobalToLocal(), Mesh::getLocalToGlobal(), MultiFieldMatrix::getMatrix(), CRMatrix< T_Diag, T_OffDiag, X >::getOffDiag(), CRConnectivity::getRow(), StorageSite::getSelfCount(), MultiFieldMatrix::hasMatrix(), Mesh::IBTYPE_FLUID, and max().

   {
     const Mesh& mesh = *_meshes[0];
     
     const StorageSite& cells = mesh.getCells();
     const IntArray& ibType = dynamic_cast<const IntArray&>(_geomFields.ibType[cells]);
     
     MultiField::ArrayIndex pIndex(&_flowFields.pressure,&cells);
     MultiField::ArrayIndex vIndex(&_flowFields.velocity,&cells);
 
     PPMatrix& ppMatrix =
       dynamic_cast<PPMatrix&>(mfmatrix.getMatrix(pIndex,pIndex));
 
     PPDiagArray& ppDiag = ppMatrix.getDiag();
     PPDiagArray& ppCoeff = ppMatrix.getOffDiag();
 
     TArray& rCell = dynamic_cast<TArray&>(rField[pIndex]);
 
     const CRConnectivity& cr = ppMatrix.getConnectivity();
 
     const Array<int>& row = cr.getRow();
 
     const int nCells = cells.getSelfCount();
 
 #ifdef FVM_PARALLEL
     const Array<int>& localToGlobal = mesh.getLocalToGlobal();
 #endif
     
     _globalRefCellID = numeric_limits<int>::max();
     for ( int n = 0; n < nCells; n++ )
     {
         if ( ibType[n] == Mesh::IBTYPE_FLUID )
         {
 #ifdef FVM_PARALLEL
             int glblIndx = localToGlobal[n];
 #else
             int glblIndx = n;
 #endif
               
             if ( glblIndx < _globalRefCellID )
               _globalRefCellID = glblIndx;
         }
     } 
     
 #ifdef FVM_PARALLEL
     MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &_globalRefCellID, 1, MPI::INT, MPI::MIN);
 
     
     _globalRefProcID = -1;
     const map<int,int>& globalToLocal = mesh.getGlobalToLocal();
 
     if ( globalToLocal.count(_globalRefCellID) > 0 ){
        _globalRefProcID = MPI::COMM_WORLD.Get_rank();
     }
     MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &_globalRefProcID, 1, MPI::INT, MPI::MAX);
     
     
     if ( globalToLocal.count(_globalRefCellID) > 0)
     {
        int nc =  globalToLocal.find(_globalRefCellID)->second;
 
 #else
        int nc = _globalRefCellID;
 #endif
        
        ppDiag[nc] = -1.;
        rCell[nc]=0.;
        for(int nb=row[nc]; nb<row[nc+1]; nb++)
          ppCoeff[nb] = 0;
 #ifdef PV_COUPLED
        if (mfmatrix.hasMatrix(pIndex,vIndex))
        {
           PVMatrix& pvMatrix =
             dynamic_cast<PVMatrix&>(mfmatrix.getMatrix(pIndex,vIndex));
         
           PVDiagArray& pvDiag = pvMatrix.getDiag();
           PVDiagArray& pvCoeff = pvMatrix.getOffDiag();
 
           pvDiag[nc] = NumTypeTraits<VectorT3>::getZero();
           for(int nb=row[nc]; nb<row[nc+1]; nb++)
              pvCoeff[nb] = NumTypeTraits<VectorT3>::getZero();
        }
 #endif
        
 #ifdef FVM_PARALLEL
     }
 #endif
   }

template<class T>

void FlowModel< T >::Impl::setReferencePP ( const MultiField & ppField )

inline

Definition at line 1200 of file FlowModel_impl.h.

References Model::_meshes, Mesh::getCells(), and Mesh::getGlobalToLocal().

   {
     if (_useReferencePressure)
     {
         const Mesh& mesh = *_meshes[0];
         
         const StorageSite& cells = mesh.getCells();
         
         MultiField::ArrayIndex pIndex(&_flowFields.pressure,&cells);
         const TArray& pp = dynamic_cast<const TArray&>(ppField[pIndex]);
 
 #ifndef FVM_PARALLEL
         _referencePP = pp[_globalRefCellID];  
 #endif
 
         
 #ifdef FVM_PARALLEL    
         _referencePP = 0.0;
         const map<int,int>& globalToLocal = mesh.getGlobalToLocal();
         if ( MPI::COMM_WORLD.Get_rank() == _globalRefProcID ){
             int localID = globalToLocal.find(_globalRefCellID)->second;
            _referencePP = pp[localID];
         }          
            
         int count = 1;
         MPI::COMM_WORLD.Bcast( &_referencePP, count, MPI::DOUBLE, _globalRefProcID);                    
 
 #endif
        //broadcast this value (referencePP) to all 
 
     }
     else
       _referencePP = 0.0;
   }

template<class T>

void FlowModel< T >::Impl::slipJumpMomentumBC	(	const StorageSite &	faces,
		const Mesh &	mesh,
		GenericBCS< VectorT3, DiagTensorT3, T > &	gbc,
		const T	accomodationCoefficient,
		const FloatValEvaluator< VectorT3 > &	bVelocity
	)

inline

Definition at line 12 of file FlowModel_impl.h.

 {
 public:
 
   typedef Array<int> IntArray;
   
   typedef Array<T> TArray;
   typedef Vector<T,3> VectorT3;
   typedef VectorTranspose<T,3> VectorT3T;
 
   typedef Array<VectorT3> VectorT3Array;
   typedef DiagonalTensor<T,3> DiagTensorT3;
 
   typedef CRMatrix<DiagTensorT3,T,VectorT3> VVMatrix;
   typedef typename VVMatrix::DiagArray VVDiagArray;
   
   typedef CRMatrix<T,T,T> PPMatrix;
   typedef typename PPMatrix::DiagArray PPDiagArray;
   typedef typename PPMatrix::PairWiseAssembler PPAssembler;
 
 #ifdef PV_COUPLED
   typedef CRMatrixRect<VectorT3T,VectorT3,T> PVMatrix;
   typedef typename PVMatrix::DiagArray PVDiagArray;
   typedef typename PVMatrix::PairWiseAssembler PVAssembler;
 
   typedef CRMatrixRect<VectorT3,T,VectorT3> VPMatrix;
   typedef typename VPMatrix::DiagArray VPDiagArray;
   typedef typename VPMatrix::PairWiseAssembler VPAssembler;
 
 #endif
   
   typedef Gradient<VectorT3> VGradType;
   typedef Array<Gradient<VectorT3> > VGradArray;
 
   typedef Gradient<T> PGradType;
   typedef Array<PGradType> PGradArray;
 
   typedef FluxJacobianMatrix<T,T> FMatrix;
   
   typedef StressTensor<T> StressTensorT6;
   typedef Array<StressTensorT6> StressTensorArray;
   
   Impl(const GeomFields& geomFields,
        FlowFields& thermalFields,
        const MeshList& meshes):
     _meshes(meshes),
     _geomFields(geomFields),
     _flowFields(thermalFields),
     _velocityGradientModel(_meshes,_flowFields.velocity,

template<class T>

MFRPtr FlowModel< T >::Impl::solveContinuity ( )

inline

Definition at line 1410 of file FlowModel_impl.h.

   {
     shared_ptr<LinearSystem> ls(discretizeContinuity());
 
     // discard previous velocity
     _previousVelocity = shared_ptr<Field>();
 
     MFRPtr rNorm = _options.getPressureLinearSolver().solve(*ls);
 
     if (!_initialContinuityNorm) _initialContinuityNorm = rNorm;
 
     ls->postSolve();
     _options.pressureLinearSolver->cleanup();
     
     postContinuitySolve(*ls);
 
     // discard the momentum ap coeffficients
     
     _momApField = shared_ptr<Field>();
     return rNorm;
   }

template<class T>

MFRPtr FlowModel< T >::Impl::solveMomentum ( )

inline

Definition at line 730 of file FlowModel_impl.h.

References Model::_meshes, Mesh::getCells(), CRMatrix< T_Diag, T_OffDiag, X >::getDiag(), LinearSystem::getMatrix(), MultiFieldMatrix::getMatrix(), LinearSystem::initAssembly(), LinearSystem::initSolve(), Array< T >::newCopy(), LinearSystem::postSolve(), and LinearSystem::updateSolution().

   {
     LinearSystem ls;
 
     initMomentumLinearization(ls);
     ls.initAssembly();
     linearizeMomentum(ls);
     ls.initSolve();
     
 
     // save current velocity for use in continuity discretization
     _previousVelocity = dynamic_pointer_cast<Field>(_flowFields.velocity.newCopy());
     
     //AMG solver(ls);
     MFRPtr rNorm = _options.getMomentumLinearSolver().solve(ls);
     
     if (!_initialMomentumNorm) _initialMomentumNorm = rNorm;
         
     _options.getMomentumLinearSolver().cleanup();
     
     ls.postSolve();
     
     ls.updateSolution();
    
     // save the momentum ap coeffficients for use in continuity discretization
     _momApField = shared_ptr<Field>(new Field("momAp"));
     const int numMeshes = _meshes.size();
     for (int n=0; n<numMeshes; n++)
     {
         const Mesh& mesh = *_meshes[n];
 
         const StorageSite& cells = mesh.getCells();
         MultiField::ArrayIndex vIndex(&_flowFields.velocity,&cells);
         const VVMatrix& vvMatrix =
           dynamic_cast<const VVMatrix&>(ls.getMatrix().getMatrix(vIndex,vIndex));
         const VVDiagArray& momAp = vvMatrix.getDiag();
         _momApField->addArray(cells,dynamic_pointer_cast<ArrayBase>(momAp.newCopy()));
     }
     _momApField->syncLocal();
     return rNorm;
   }

template<class T>

void FlowModel< T >::Impl::updateFacePressureBoundary	(	const Mesh &	mesh,
		const StorageSite &	faces
	)

inline

Definition at line 219 of file FlowModel_impl.h.

237 {

template<class T>

void FlowModel< T >::Impl::updateFacePressureInterior	(	const Mesh &	mesh,
		const StorageSite &	faces,
		const bool	isSymmetry = `false`
	)

inline

Definition at line 304 of file FlowModel_impl.h.

References dot().

                 {
                     const int c0 = faceCells(f,0);
                     massFlux[f] = density[c0]*dot(bVelocity[f],faceArea[f]);
                  }
             }
             else if (bc.bcType == "SpecifiedPressure")
             {
                 T bp = bc["specifiedPressure"];
                 TArray& facePressure = dynamic_cast<TArray&>(_flowFields.pressure[faces]);
                 TArray& cellPressure = dynamic_cast<TArray&>(_flowFields.pressure[cells]);
                 for(int f=0; f<nFaces; f++)
                 {
                     const int c1 = faceCells(f,1);
                     facePressure[f] = cellPressure[c1] = bp;
                 }
             }
             else if ((bc.bcType == "Symmetry"))
             {
                 TArray& massFlux = dynamic_cast<TArray&>(_flowFields.massFlux[faces]);
                 massFlux.zero();
             }
             shared_ptr<VectorT3Array> momFlux(new VectorT3Array(faces.getCount()));
             momFlux->zero();
             _flowFields.momentumFlux.addArray(faces,momFlux);
            
             
             
         }
     }
 
     computeContinuityResidual();
     _niters  =0;
     _initialMomentumNorm = MFRPtr();
     _initialContinuityNorm = MFRPtr();
   }
   
   FlowBCMap& getBCMap() {return _bcMap;}
   FlowVCMap& getVCMap() {return _vcMap;}
   FlowModelOptions<T>& getOptions() {return _options;}
 
   void updateTime()
   {
     const int numMeshes = _meshes.size();
     for (int n=0; n<numMeshes; n++)
     {
         const Mesh& mesh = *_meshes[n];
 
         const StorageSite& cells = mesh.getCells();
         VectorT3Array& v =
           dynamic_cast<VectorT3Array&>(_flowFields.velocity[cells]);
         VectorT3Array& vN1 =
           dynamic_cast<VectorT3Array&>(_flowFields.velocityN1[cells]);
 
         if (_options.timeDiscretizationOrder > 1)
         {
             VectorT3Array& vN2 =
               dynamic_cast<VectorT3Array&>(_flowFields.velocityN2[cells]);
             vN2 = vN1;
         }
         vN1 = v;
     }
   }
  const Field& getViscosityField() const

template<class T>

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

inline

Definition at line 348 of file FlowModel_impl.h.

References Model::_meshes, and Mesh::getCells().

   {
     const int numMeshes = _meshes.size();
     for (int n=0; n<numMeshes; n++)
     {
         const Mesh& mesh = *_meshes[n];
 
         const StorageSite& cells = mesh.getCells();
         VectorT3Array& v =
           dynamic_cast<VectorT3Array&>(_flowFields.velocity[cells]);
         VectorT3Array& vN1 =
           dynamic_cast<VectorT3Array&>(_flowFields.velocityN1[cells]);
 
         if (_options.timeDiscretizationOrder > 1)
         {
             VectorT3Array& vN2 =
               dynamic_cast<VectorT3Array&>(_flowFields.velocityN2[cells]);
             vN2 = vN1;
         }
         vN1 = v;
     }
   }

Member Data Documentation

template<class T>

FlowBCMap FlowModel< T >::Impl::_bcMap

private

Definition at line 2125 of file FlowModel_impl.h.

template<class T>

FlowFields& FlowModel< T >::Impl::_flowFields

private

Definition at line 2124 of file FlowModel_impl.h.

template<class T>

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

private

Definition at line 2123 of file FlowModel_impl.h.

template<class T>

int FlowModel< T >::Impl::_globalRefCellID

private

Definition at line 2141 of file FlowModel_impl.h.

template<class T>

int FlowModel< T >::Impl::_globalRefProcID

private

Definition at line 2142 of file FlowModel_impl.h.

template<class T>

MFRPtr FlowModel< T >::Impl::_initialContinuityNorm

private

Definition at line 2133 of file FlowModel_impl.h.

template<class T>

MFRPtr FlowModel< T >::Impl::_initialCoupledNorm

private

Definition at line 2134 of file FlowModel_impl.h.

template<class T>

MFRPtr FlowModel< T >::Impl::_initialMomentumNorm

private

Definition at line 2132 of file FlowModel_impl.h.

template<class T>

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

private

Definition at line 2122 of file FlowModel_impl.h.

template<class T>

shared_ptr<Field> FlowModel< T >::Impl::_momApField

private

Definition at line 2138 of file FlowModel_impl.h.

template<class T>

int FlowModel< T >::Impl::_niters

private

Definition at line 2135 of file FlowModel_impl.h.

template<class T>

FlowModelOptions<T> FlowModel< T >::Impl::_options

private

Definition at line 2128 of file FlowModel_impl.h.

template<class T>

map<string,shared_ptr<ArrayBase> > FlowModel< T >::Impl::_persistenceData

private

Definition at line 2146 of file FlowModel_impl.h.

template<class T>

GradientModel<T> FlowModel< T >::Impl::_pressureGradientModel

private

Definition at line 2130 of file FlowModel_impl.h.

template<class T>

shared_ptr<Field> FlowModel< T >::Impl::_previousVelocity

private

Definition at line 2137 of file FlowModel_impl.h.

template<class T>

T FlowModel< T >::Impl::_referencePP

private

Definition at line 2143 of file FlowModel_impl.h.

template<class T>

bool FlowModel< T >::Impl::_useReferencePressure

private

Definition at line 2140 of file FlowModel_impl.h.

template<class T>

FlowVCMap FlowModel< T >::Impl::_vcMap

private

Definition at line 2126 of file FlowModel_impl.h.

template<class T>

GradientModel<VectorT3> FlowModel< T >::Impl::_velocityGradientModel

private

Definition at line 2129 of file FlowModel_impl.h.

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

FlowModel_impl.h

Public Types

Public Member Functions

Private Attributes

Detailed Description

template<class T> class FlowModel< T >::Impl

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation

template<class T>
class FlowModel< T >::Impl