Commit adc90617 authored by andy's avatar andy
Browse files

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

parents 09ee2342 dd017458
......@@ -27,7 +27,7 @@
fvScalarMatrix pDDtEqn
(
fvc::ddt(rho) + fvc::div(phiHbyA)
+ correction(fvm::ddt(psi, p) + fvm::div(phid, p))
+ correction(psi*fvm::ddt(p) + fvm::div(phid, p))
);
while (pimple.correctNonOrthogonal())
......
......@@ -68,4 +68,6 @@
}
}
phi.oldTime() = phi;
#include "continuityErrs.H"
{
solve
fvScalarMatrix TEqn
(
fvm::ddt(rho, T)
+ fvm::div(rhoPhi, T)
- fvm::laplacian(twoPhaseProperties.alphaEff(turbulence->mut()), T)
+ (
p*fvc::div(phi)
fvc::div(fvc::absolute(phi, U), p)
+ fvc::ddt(rho, K) + fvc::div(rhoPhi, K)
)
*(
......@@ -14,5 +14,8 @@
)
);
TEqn.relax();
TEqn.solve();
twoPhaseProperties.correct();
}
......@@ -56,6 +56,8 @@
phi -= pcorrEqn.flux();
}
}
#include "continuityErrs.H"
}
phi.oldTime() = phi;
#include "continuityErrs.H"
......@@ -104,8 +104,6 @@ int main(int argc, char *argv[])
}
}
rho = alpha1*rho1 + alpha2*rho2;
runTime.write();
Info<< "ExecutionTime = "
......
......@@ -38,9 +38,9 @@
volScalarField& p = twoPhaseProperties.p();
volScalarField& T = twoPhaseProperties.T();
const volScalarField& rho1 = twoPhaseProperties.thermo1().rho();
volScalarField& rho1 = twoPhaseProperties.thermo1().rho();
const volScalarField& psi1 = twoPhaseProperties.thermo1().psi();
const volScalarField& rho2 = twoPhaseProperties.thermo2().rho();
volScalarField& rho2 = twoPhaseProperties.thermo2().rho();
const volScalarField& psi2 = twoPhaseProperties.thermo2().psi();
volScalarField rho
......@@ -93,18 +93,5 @@
compressible::turbulenceModel::New(rho, U, rhoPhi, twoPhaseProperties)
);
Info<< "Creating field dpdt\n" << endl;
volScalarField dpdt
(
IOobject
(
"dpdt",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
);
Info<< "Creating field kinetic energy K\n" << endl;
volScalarField K("K", 0.5*magSqr(U));
......@@ -26,28 +26,44 @@
tmp<fvScalarMatrix> p_rghEqnComp1;
tmp<fvScalarMatrix> p_rghEqnComp2;
//if (transonic)
//{
//}
//else
if (pimple.transonic())
{
surfaceScalarField phid1("phid1", fvc::interpolate(psi1)*phi);
surfaceScalarField phid2("phid2", fvc::interpolate(psi2)*phi);
p_rghEqnComp1 =
fvc::ddt(rho1) + fvc::div(phi, rho1) - fvc::Sp(fvc::div(phi), rho1)
+ correction
(
psi1*fvm::ddt(p_rgh)
+ fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
);
deleteDemandDrivenData(p_rghEqnComp1().faceFluxCorrectionPtr());
p_rghEqnComp1().relax();
p_rghEqnComp2 =
fvc::ddt(rho2) + fvc::div(phi, rho2) - fvc::Sp(fvc::div(phi), rho2)
+ correction
(
psi2*fvm::ddt(p_rgh)
+ fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
);
deleteDemandDrivenData(p_rghEqnComp2().faceFluxCorrectionPtr());
p_rghEqnComp2().relax();
}
else
{
p_rghEqnComp1 =
fvc::ddt(rho1) + psi1*correction(fvm::ddt(p_rgh))
+ fvc::div(phid1, p_rgh)
- fvc::Sp(fvc::div(phid1), p_rgh);
+ fvc::div(phi, rho1) - fvc::Sp(fvc::div(phi), rho1);
p_rghEqnComp2 =
fvc::ddt(rho2) + psi2*correction(fvm::ddt(p_rgh))
+ fvc::div(phid2, p_rgh)
- fvc::Sp(fvc::div(phid2), p_rgh);
+ fvc::div(phi, rho2) - fvc::Sp(fvc::div(phi), rho2);
}
// Thermodynamic density needs to be updated by psi*d(p) after the
// pressure solution - done in 2 parts. Part 1:
//twoPhaseProperties.rho() -= psi*p_rgh;
// Cache p_rgh prior to solve for density update
volScalarField p_rgh_0(p_rgh);
while (pimple.correctNonOrthogonal())
{
......@@ -69,8 +85,8 @@
if (pimple.finalNonOrthogonalIter())
{
// Second part of thermodynamic density update
//twoPhaseProperties.rho() += psi*p_rgh;
//p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
//p_rgh = p - (alpha1*rho1 + alpha2*rho2)*gh;
dgdt =
(
......@@ -88,15 +104,14 @@
p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
// twoPhaseProperties.correct();
// Update densities from change in p_rgh
rho1 += psi1*(p_rgh - p_rgh_0);
rho2 += psi2*(p_rgh - p_rgh_0);
Info<< "max(U) " << max(mag(U)).value() << endl;
Info<< "min(p_rgh) " << min(p_rgh).value() << endl;
rho = alpha1*rho1 + alpha2*rho2;
K = 0.5*magSqr(U);
if (twoPhaseProperties.dpdt())
{
dpdt = fvc::ddt(p);
}
Info<< "max(U) " << max(mag(U)).value() << endl;
Info<< "min(p_rgh) " << min(p_rgh).value() << endl;
}
......@@ -94,6 +94,16 @@ public:
return thermo2_();
}
rhoThermo& thermo1()
{
return thermo1_();
}
rhoThermo& thermo2()
{
return thermo2_();
}
//- Update properties
virtual void correct();
......
......@@ -62,10 +62,10 @@
volScalarField& p = thermo1.p();
volScalarField rho1("rho" + phase1Name, thermo1.rho());
volScalarField& rho1 = thermo1.rho();
const volScalarField& psi1 = thermo1.psi();
volScalarField rho2("rho" + phase2Name, thermo2.rho());
volScalarField& rho2 = thermo2.rho();
const volScalarField& psi2 = thermo2.psi();
volVectorField U
......
{
rho1 = thermo1.rho();
rho2 = thermo2.rho();
surfaceScalarField alpha1f(fvc::interpolate(alpha1));
surfaceScalarField alpha2f(scalar(1) - alpha1f);
......@@ -91,10 +88,7 @@
tmp<fvScalarMatrix> pEqnComp1;
tmp<fvScalarMatrix> pEqnComp2;
//if (transonic)
//{
//}
//else
if (pimple.transonic())
{
surfaceScalarField phid1
(
......@@ -107,17 +101,42 @@
fvc::interpolate(psi2)*phi2
);
pEqnComp1 =
fvc::ddt(rho1)
+ fvc::div(phi1, rho1) - fvc::Sp(fvc::div(phi1), rho1)
+ correction
(
psi1*fvm::ddt(p)
+ fvm::div(phid1, p) - fvm::Sp(fvc::div(phid1), p)
);
deleteDemandDrivenData(pEqnComp1().faceFluxCorrectionPtr());
pEqnComp1().relax();
pEqnComp2 =
fvc::ddt(rho2)
+ fvc::div(phi2, rho2) - fvc::Sp(fvc::div(phi2), rho2)
+ correction
(
psi2*fvm::ddt(p)
+ fvm::div(phid2, p) - fvm::Sp(fvc::div(phid2), p)
);
deleteDemandDrivenData(pEqnComp2().faceFluxCorrectionPtr());
pEqnComp2().relax();
}
else
{
pEqnComp1 =
fvc::ddt(rho1) + psi1*correction(fvm::ddt(p))
+ fvc::div(phid1, p)
- fvc::Sp(fvc::div(phid1), p);
+ fvc::div(phi1, rho1) - fvc::Sp(fvc::div(phi1), rho1);
pEqnComp2 =
fvc::ddt(rho2) + psi2*correction(fvm::ddt(p))
+ fvc::div(phid2, p)
- fvc::Sp(fvc::div(phid2), p);
+ fvc::div(phi2, rho2) - fvc::Sp(fvc::div(phi2), rho2);
}
// Cache p prior to solve for density update
volScalarField p_0(p);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqnIncomp
......@@ -182,10 +201,9 @@
p = max(p, pMin);
thermo1.correct();
thermo2.correct();
rho1 = thermo1.rho();
rho2 = thermo2.rho();
// Update densities from change in p
rho1 += psi1*(p - p_0);
rho2 += psi2*(p - p_0);
K1 = 0.5*magSqr(U1);
K2 = 0.5*magSqr(U2);
......
......@@ -53,6 +53,8 @@
fvc::makeRelative(phi, U);
}
}
#include "continuityErrs.H"
}
phi.oldTime() = phi;
#include "continuityErrs.H"
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -263,8 +263,8 @@ void Foam::sigFpe::set(const bool verbose)
{
if (supported)
{
Info<< "SetNaN : Initialising allocated memory to NaN"
<< " (FOAM_SETNAN)." << endl;
Info<< "SetNaN : Initialising allocated memory to NaN"
<< " (FOAM_SETNAN)." << endl;
}
else
{
......
......@@ -367,6 +367,8 @@ $(cellShape)/cellShapeIOList.C
meshes/Identifiers/patch/patchIdentifier.C
meshes/MeshObject/meshObject.C
polyMesh = meshes/polyMesh
polyPatches = $(polyMesh)/polyPatches
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -132,7 +132,9 @@ inline const T& Foam::PtrList<T>::operator[](const label i) const
if (!ptrs_[i])
{
FatalErrorIn("PtrList::operator[] const")
<< "hanging pointer, cannot dereference"
<< "hanging pointer at index " << i
<< " (size " << size()
<< "), cannot dereference"
<< abort(FatalError);
}
......@@ -146,7 +148,9 @@ inline T& Foam::PtrList<T>::operator[](const label i)
if (!ptrs_[i])
{
FatalErrorIn("PtrList::operator[]")
<< "hanging pointer, cannot dereference"
<< "hanging pointer at index " << i
<< " (size " << size()
<< "), cannot dereference"
<< abort(FatalError);
}
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -107,7 +107,9 @@ inline const T& Foam::UPtrList<T>::operator[](const label i) const
if (!ptrs_[i])
{
FatalErrorIn("UPtrList::operator[] const")
<< "hanging pointer, cannot dereference"
<< "hanging pointer at index " << i
<< " (size " << size()
<< "), cannot dereference"
<< abort(FatalError);
}
......@@ -121,7 +123,9 @@ inline T& Foam::UPtrList<T>::operator[](const label i)
if (!ptrs_[i])
{
FatalErrorIn("UPtrList::operator[]")
<< "hanging pointer, cannot dereference"
<< "hanging pointer at index " << i
<< " (size " << size()
<< "), cannot dereference"
<< abort(FatalError);
}
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -24,7 +24,6 @@ License
\*---------------------------------------------------------------------------*/
#include "processorPointPatchField.H"
//#include "transformField.H"
#include "processorPolyPatch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -92,73 +91,6 @@ processorPointPatchField<Type>::~processorPointPatchField()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
void processorPointPatchField<Type>::initSwapAddSeparated
(
const Pstream::commsTypes commsType,
Field<Type>& pField
)
const
{
// if (Pstream::parRun())
// {
// // Get internal field into correct order for opposite side
// Field<Type> pf
// (
// this->patchInternalField
// (
// pField,
// procPatch_.reverseMeshPoints()
// )
// );
//
// OPstream::write
// (
// commsType,
// procPatch_.neighbProcNo(),
// reinterpret_cast<const char*>(pf.begin()),
// pf.byteSize(),
// procPatch_.tag()
// );
// }
}
template<class Type>
void processorPointPatchField<Type>::swapAddSeparated
(
const Pstream::commsTypes commsType,
Field<Type>& pField
) const
{
// if (Pstream::parRun())
// {
// Field<Type> pnf(this->size());
//
// IPstream::read
// (
// commsType,
// procPatch_.neighbProcNo(),
// reinterpret_cast<char*>(pnf.begin()),
// pnf.byteSize(),
// procPatch_.tag()
// );
//
// if (doTransform())
// {
// const processorPolyPatch& ppp = procPatch_.procPolyPatch();
// const tensor& forwardT = ppp.forwardT();
//
// transform(pnf, forwardT, pnf);
// }
//
// addToInternalField(pField, pnf, procPatch_.separatedPoints());
// }
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -57,7 +57,6 @@ class processorPointPatchField
//- Local reference to processor patch
const processorPointPatch& procPatch_;
public:
//- Runtime type information
......@@ -176,19 +175,13 @@ public:
)
{}
//- Initialise swap of non-collocated patch point values
virtual void initSwapAddSeparated
(
const Pstream::commsTypes commsType,
Field<Type>&
) const;
//- Complete swap of patch point values and add to local values
//- Assume processor patch always collocated
virtual void swapAddSeparated
(
const Pstream::commsTypes commsType,
Field<Type>&
) const;
) const
{}
};
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -38,7 +38,8 @@ Foam::processorCyclicPointPatchField<Type>::processorCyclicPointPatchField
)
:
coupledPointPatchField<Type>(p, iF),
procPatch_(refCast<const processorCyclicPointPatch>(p))
procPatch_(refCast<const processorCyclicPointPatch>(p)),
receiveBuf_(0)
{}
......@@ -51,7 +52,8 @@ Foam::processorCyclicPointPatchField<Type>::processorCyclicPointPatchField
)
:
coupledPointPatchField<Type>(p, iF, dict),
procPatch_(refCast<const processorCyclicPointPatch>(p))
procPatch_(refCast<const processorCyclicPointPatch>(p)),
receiveBuf_(0)
{}
......@@ -65,7 +67,8 @@ Foam::processorCyclicPointPatchField<Type>::processorCyclicPointPatchField
)
:
coupledPointPatchField<Type>(ptf, p, iF, mapper),
procPatch_(refCast<const processorCyclicPointPatch>(ptf.patch()))
procPatch_(refCast<const processorCyclicPointPatch>(ptf.patch())),
receiveBuf_(0)
{}
......@@ -77,7 +80,8 @@ Foam::processorCyclicPointPatchField<Type>::processorCyclicPointPatchField
)
:
coupledPointPatchField<Type>(ptf, iF),
procPatch_(refCast<const processorCyclicPointPatch>(ptf.patch()))
procPatch_(refCast<const processorCyclicPointPatch>(ptf.patch())),
receiveBuf_(0)
{}
......@@ -109,6 +113,18 @@ void Foam::processorCyclicPointPatchField<Type>::initSwapAddSeparated
)
);
if (commsType == Pstream::nonBlocking)
{
receiveBuf_.setSize(pf.size());
IPstream::read
(
commsType,
procPatch_.neighbProcNo(),
reinterpret_cast<char*>(receiveBuf_.begin()),
receiveBuf_.byteSize(),
procPatch_.tag()
);
}
OPstream::write
(