Commit b4c27879 authored by mattijs's avatar mattijs
Browse files

Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

parents 727b278f 6de3559b
......@@ -164,8 +164,21 @@ void Foam::KRR4::solve
if (x == xTemp)
{
FatalErrorIn("ODES::KRR4")
<< "stepsize not significant"
FatalErrorIn
(
"void Foam::KRR4::solve"
"("
"const ODE&, "
"scalar&, "
"scalarField&, "
"scalarField&, "
"const scalar, "
"const scalarField&, "
"const scalar, "
"scalar&, "
"scalar&"
") const"
) << "solver stalled: step size = 0"
<< exit(FatalError);
}
......@@ -189,8 +202,21 @@ void Foam::KRR4::solve
}
}
FatalErrorIn("ODES::KRR4")
<< "exceeded maxtry"
FatalErrorIn
(
"void Foam::KRR4::solve"
"("
"const ODE&, "
"scalar&, "
"scalarField&, "
"scalarField&, "
"const scalar, "
"const scalarField&, "
"const scalar, "
"scalar&, "
"scalar&"
") const"
) << "Maximum number of solver iterations exceeded"
<< exit(FatalError);
}
......
......@@ -44,6 +44,7 @@ License
makeInjectionModel(CloudType); \
makeInjectionModelType(CellZoneInjection, CloudType); \
makeInjectionModelType(ConeInjection, CloudType); \
makeInjectionModelType(ConeNozzleInjection, CloudType); \
makeInjectionModelType(FieldActivatedInjection, CloudType); \
makeInjectionModelType(ManualInjection, CloudType); \
makeInjectionModelType(NoInjection, CloudType); \
......
......@@ -31,32 +31,6 @@ License
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class CloudType>
void Foam::FacePostProcessing<CloudType>::applyToFace
(
const label faceIn,
label& zoneI,
label& faceI
) const
{
const faceZoneMesh& fzm = this->owner().mesh().faceZones();
forAll(faceZoneIDs_, i)
{
const faceZone& fz = fzm[faceZoneIDs_[i]];
forAll(fz, j)
{
if (fz[j] == faceIn)
{
zoneI = i;
faceI = j;
return;
}
}
}
}
template<class CloudType>
void Foam::FacePostProcessing<CloudType>::makeLogFile
(
......@@ -93,7 +67,7 @@ void Foam::FacePostProcessing<CloudType>::makeLogFile
<< "# Face zone : " << zoneName << nl
<< "# Faces : " << nFaces << nl
<< "# Area : " << totArea << nl
<< "# Time" << tab << "mass" << tab << "massFlux" << endl;
<< "# Time" << tab << "mass" << tab << "massFlowRate" << endl;
}
}
}
......@@ -105,25 +79,27 @@ void Foam::FacePostProcessing<CloudType>::write()
const fvMesh& mesh = this->owner().mesh();
const Time& time = mesh.time();
const faceZoneMesh& fzm = mesh.faceZones();
const scalar dt = time.deltaTValue();
scalar timeNew = time.value();
scalar timeElapsed = timeNew-timeOld_;
totalTime_ += dt;
totalTime_ += timeElapsed;
const scalar alpha = (totalTime_ - dt)/totalTime_;
const scalar beta = dt/totalTime_;
const scalar alpha = (totalTime_ - timeElapsed)/totalTime_;
const scalar beta = timeElapsed/totalTime_; //correct
forAll(faceZoneIDs_, zoneI)
{
massFlowRate_[zoneI] =
alpha*massFlowRate_[zoneI] + beta*mass_[zoneI]/timeElapsed;
massTotal_[zoneI] += mass_[zoneI];
massFlux_[zoneI] = alpha*massFlux_[zoneI] + beta*mass_[zoneI]/dt;
}
const label procI = Pstream::myProcNo();
Info<< "particleFaceFlux output:" << nl;
Info<< "particleFaceFlowRate output:" << nl;
List<scalarField> zoneMassTotal(mass_.size());
List<scalarField> zoneMassFlux(massFlux_.size());
List<scalarField> zoneMassFlowRate(massFlowRate_.size());
forAll(faceZoneIDs_, zoneI)
{
const word& zoneName = fzm[faceZoneIDs_[zoneI]].name();
......@@ -138,26 +114,26 @@ void Foam::FacePostProcessing<CloudType>::write()
);
const scalar sumMassTotal = sum(zoneMassTotal[zoneI]);
scalarListList allProcMassFlux(Pstream::nProcs());
allProcMassFlux[procI] = massFlux_[zoneI];
Pstream::gatherList(allProcMassFlux);
zoneMassFlux[zoneI] =
scalarListList allProcMassFlowRate(Pstream::nProcs());
allProcMassFlowRate[procI] = massFlowRate_[zoneI];
Pstream::gatherList(allProcMassFlowRate);
zoneMassFlowRate[zoneI] =
ListListOps::combine<scalarList>
(
allProcMassFlux, accessOp<scalarList>()
allProcMassFlowRate, accessOp<scalarList>()
);
const scalar sumMassFlux = sum(zoneMassFlux[zoneI]);
const scalar sumMassFlowRate = sum(zoneMassFlowRate[zoneI]);
Info<< " " << zoneName
<< ": total mass = " << sumMassTotal
<< "; average mass flux = " << sumMassFlux
<< "; average mass flow rate = " << sumMassFlowRate
<< nl;
if (outputFilePtr_.set(zoneI))
{
OFstream& os = outputFilePtr_[zoneI];
os << time.timeName() << token::TAB << sumMassTotal << token::TAB
<< sumMassFlux<< endl;
<< sumMassFlowRate<< endl;
}
}
......@@ -235,8 +211,8 @@ void Foam::FacePostProcessing<CloudType>::write()
fZone.name(),
allPoints,
allFaces,
"massFlux",
zoneMassFlux[zoneI],
"massFlowRate",
zoneMassFlowRate[zoneI],
false
);
}
......@@ -248,8 +224,9 @@ void Foam::FacePostProcessing<CloudType>::write()
{
forAll(faceZoneIDs_, zoneI)
{
massFlux_[zoneI] = 0.0;
massFlowRate_[zoneI] = 0.0;
}
timeOld_ = timeNew;
totalTime_ = 0.0;
}
......@@ -278,15 +255,16 @@ Foam::FacePostProcessing<CloudType>::FacePostProcessing
totalTime_(0.0),
mass_(),
massTotal_(),
massFlux_(),
massFlowRate_(),
log_(this->coeffDict().lookup("log")),
outputFilePtr_(),
outputDir_(owner.mesh().time().path())
outputDir_(owner.mesh().time().path()),
timeOld_(owner.mesh().time().value())
{
wordList faceZoneNames(this->coeffDict().lookup("faceZones"));
mass_.setSize(faceZoneNames.size());
massTotal_.setSize(faceZoneNames.size());
massFlux_.setSize(faceZoneNames.size());
massFlowRate_.setSize(faceZoneNames.size());
outputFilePtr_.setSize(faceZoneNames.size());
......@@ -317,7 +295,7 @@ Foam::FacePostProcessing<CloudType>::FacePostProcessing
const faceZone& fz = fzm[zoneI];
mass_[i].setSize(fz.size(), 0.0);
massTotal_[i].setSize(fz.size(), 0.0);
massFlux_[i].setSize(fz.size(), 0.0);
massFlowRate_[i].setSize(fz.size(), 0.0);
label nFaces = returnReduce(fz.size(), sumOp<label>());
Info<< " " << zoneName << " faces: " << nFaces << nl;
......@@ -372,10 +350,11 @@ Foam::FacePostProcessing<CloudType>::FacePostProcessing
totalTime_(pff.totalTime_),
mass_(pff.mass_),
massTotal_(pff.massTotal_),
massFlux_(pff.massFlux_),
massFlowRate_(pff.massFlowRate_),
log_(pff.log_),
outputFilePtr_(),
outputDir_(pff.outputDir_)
outputDir_(pff.outputDir_),
timeOld_(0.0)
{}
......@@ -401,13 +380,26 @@ void Foam::FacePostProcessing<CloudType>::postFace
|| this->owner().solution().transient()
)
{
label zoneId = -1;
label faceId = -1;
applyToFace(faceI, zoneId, faceId);
const faceZoneMesh& fzm = this->owner().mesh().faceZones();
if ((zoneId != -1) && (faceId != -1))
forAll(faceZoneIDs_, i)
{
mass_[zoneId][faceId] += p.mass()*p.nParticle();
const faceZone& fz = fzm[faceZoneIDs_[i]];
label faceId = -1;
forAll(fz, j)
{
if (fz[j] == faceI)
{
faceId = j;
break;
}
}
if (faceId != -1)
{
mass_[i][faceId] += p.mass()*p.nParticle();
}
}
}
}
......
......@@ -82,8 +82,8 @@ class FacePostProcessing
//- Mass total storage
List<scalarField> massTotal_;
//- Mass flux storage
List<scalarField> massFlux_;
//- Mass flow rate storage
List<scalarField> massFlowRate_;
//- Flag to indicate whether data should be written to file
Switch log_;
......@@ -94,6 +94,9 @@ class FacePostProcessing
//- Output directory
fileName outputDir_;
//- Last calculation time
scalar timeOld_;
// Private Member Functions
......@@ -106,14 +109,6 @@ class FacePostProcessing
const scalar totArea
);
//- Return index into storage lists if valid zone and face
void applyToFace
(
const label faceIn,
label& zoneI, label&
faceI
) const;
protected:
......
......@@ -46,6 +46,7 @@ void Foam::faceAreaIntersect::triSliceWithPlane
label nPos = 0;
label posI = -1;
label negI = -1;
label copI = -1;
forAll(tri, i)
{
d[i] = ((tri[i] - p.refPoint()) & p.normal());
......@@ -54,6 +55,7 @@ void Foam::faceAreaIntersect::triSliceWithPlane
{
d[i] = 0.0;
nCoPlanar++;
copI = i;
}
else
{
......@@ -69,12 +71,16 @@ void Foam::faceAreaIntersect::triSliceWithPlane
}
}
if ((nPos == 3) || ((nPos == 2) && (nCoPlanar == 1)))
if
(
(nPos == 3)
|| ((nPos == 2) && (nCoPlanar == 1))
|| ((nPos == 1) && (nCoPlanar == 2)))
{
// all points above cutting plane - add triangle to list
tris[nTris++] = tri;
}
else if ((nPos == 2) || ((nPos == 1) && (nCoPlanar == 1)))
else if ((nPos == 2) && (nCoPlanar == 0))
{
// 2 points above plane, 1 below
// resulting quad above plane split into 2 triangles
......@@ -95,25 +101,48 @@ void Foam::faceAreaIntersect::triSliceWithPlane
setTriPoints(tri[i1], tri[i2], p02, nTris, tris);
setTriPoints(tri[i1], p02, p01, nTris, tris);
}
else if ((nPos == 1) && (nCoPlanar != 1))
else if (nPos == 1)
{
// 1 point above plane, 2 below
// resulting quad below plane split into 2 triangles
// point above the plane
label i0 = posI;
// indices of remaining points
label i1 = d.fcIndex(i0);
label i2 = d.fcIndex(i1);
if (nCoPlanar == 0)
{
// 1 point above plane, 2 below
// determine the two intersection points
point p01 = planeIntersection(d, tri, i0, i1);
point p02 = planeIntersection(d, tri, i0, i2);
// indices of remaining points
label i1 = d.fcIndex(i0);
label i2 = d.fcIndex(i1);
// determine the two intersection points
point p01 = planeIntersection(d, tri, i0, i1);
point p02 = planeIntersection(d, tri, i0, i2);
// forget quad below plane
// - add triangle above plane to list
setTriPoints(tri[i0], p01, p02, nTris, tris);
// forget quad below plane
// - add triangle above plane to list
setTriPoints(tri[i0], p01, p02, nTris, tris);
}
else
{
// 1 point above plane, 1 on plane, 1 below
// point indices
label i1 = negI;
label i2 = copI;
// determine the intersection point
point p01 = planeIntersection(d, tri, i0, i1);
// add triangle above plane to list
if (d.fcIndex(i0) == i1)
{
setTriPoints(tri[i0], p01, tri[i2], nTris, tris);
}
else
{
setTriPoints(tri[i0], tri[i2], p01, nTris, tris);
}
}
}
else
{
......
#!/bin/sh
cd ${0%/*} || exit 1 # run from this directory
makeType=${1:-libso}
set -x
wclean $makeType regionModel
wclean $makeType pyrolysisModels
wclean $makeType surfaceFilmModels
wclean $makeType thermoBaffleModels
wclean $makeType regionCoupling
# ----------------------------------------------------------------- end-of-file
......@@ -35,21 +35,21 @@ void Foam::janafThermo<EquationOfState>::checkInputData() const
{
FatalErrorIn("janafThermo<EquationOfState>::check()")
<< "Tlow(" << Tlow_ << ") >= Thigh(" << Thigh_ << ')'
<< exit(FatalIOError);
<< exit(FatalError);
}
if (Tcommon_ <= Tlow_)
{
FatalErrorIn("janafThermo<EquationOfState>::check()")
<< "Tcommon(" << Tcommon_ << ") <= Tlow(" << Tlow_ << ')'
<< exit(FatalIOError);
<< exit(FatalError);
}
if (Tcommon_ > Thigh_)
{
FatalErrorIn("janafThermo<EquationOfState>::check()")
<< "Tcommon(" << Tcommon_ << ") > Thigh(" << Thigh_ << ')'
<< exit(FatalIOError);
<< exit(FatalError);
}
}
......
......@@ -153,6 +153,25 @@ public:
//- Limit the temperature to be in the range Tlow_ to Thigh_
inline scalar limit(const scalar T) const;
// Access
//- Return const access to the low temperature limit
inline scalar Tlow() const;
//- Return const access to the high temperature limit
inline scalar Thigh() const;
//- Return const access to the common temperature
inline scalar Tcommon() const;
//- Return const access to the high temperature poly coefficients
inline const coeffArray& highCpCoeffs() const;
//- Return const access to the low temperature poly coefficients
inline const coeffArray& lowCpCoeffs() const;
// Fundamental properties
//- Heat capacity at constant pressure [J/(kmol K)]
......
......@@ -119,6 +119,43 @@ inline Foam::scalar Foam::janafThermo<EquationOfState>::limit
}
template<class EquationOfState>
inline Foam::scalar Foam::janafThermo<EquationOfState>::Tlow() const
{
return Tlow_;
}
template<class EquationOfState>
inline Foam::scalar Foam::janafThermo<EquationOfState>::Thigh() const
{
return Thigh_;
}
template<class EquationOfState>
inline Foam::scalar Foam::janafThermo<EquationOfState>::Tcommon() const
{
return Tcommon_;
}
template<class EquationOfState>
inline const typename Foam::janafThermo<EquationOfState>::coeffArray&
Foam::janafThermo<EquationOfState>::highCpCoeffs() const
{
return highCpCoeffs_;
}
template<class EquationOfState>
inline const typename Foam::janafThermo<EquationOfState>::coeffArray&
Foam::janafThermo<EquationOfState>::lowCpCoeffs() const
{
return lowCpCoeffs_;
}
template<class EquationOfState>
inline Foam::scalar Foam::janafThermo<EquationOfState>::cp
(
......
......@@ -17,7 +17,7 @@ FoamFile
dimensions [ 0 -1 0 0 0 0 0 ];
internalField nonuniform List<scalar>
internalField nonuniform List<scalar>
8025
(
18.372496
......
......@@ -17,7 +17,7 @@ FoamFile
dimensions [ 0 0 0 0 0 0 0 ];
internalField nonuniform List<symmTensor>
internalField nonuniform List<symmTensor>
8025
(
(1.785697 0 0 1.785697 0 1.875)
......
......@@ -17,7 +17,7 @@ FoamFile
dimensions [ 0 -1 0 0 0 0 0 ];
internalField nonuniform List<symmTensor>
internalField nonuniform List<symmTensor>
8025
(
(98.997697 0 0 77.711245 0 30.754852)
......
......@@ -17,7 +17,7 @@ FoamFile
dimensions [ 0 -1 0 0 0 0 0 ];
internalField nonuniform List<symmTensor>
internalField nonuniform List<symmTensor>
8025
(
(4.501854 0 0 3.437531 0 1.537743)
......
......@@ -17,7 +17,7 @@ FoamFile
dimensions [ 0 1 0 0 0 0 0 ];
internalField nonuniform List<scalar>
internalField nonuniform List<scalar>
8025
(
0.2
......
......@@ -17,7 +17,7 @@ FoamFile
dimensions [ 0 0 0 0 0 0 0 ];
internalField nonuniform List<scalar>
internalField nonuniform List<scalar>
8025
(
2.019581
......
......@@ -17,7 +17,7 @@ FoamFile
dimensions [ 0 0 0 0 0 0 0 ];
internalField nonuniform List<scalar>
internalField nonuniform List<scalar>
8025
(
0
......
......@@ -17,7 +17,7 @@ FoamFile
dimensions [ 0 0 0 0 0 0 0 ];
internalField nonuniform List<symmTensor>
internalField nonuniform List<symmTensor>
8025
(
(0.634791 0 0 0.634791 0 0.03)
......