Merged volume continuity check files

This commit is contained in:
Hrvoje Jasak 2018-10-12 10:37:05 +01:00
parent 850f11a9fc
commit 6b09536c33
4 changed files with 49 additions and 75 deletions

View file

@ -68,7 +68,7 @@ int main(int argc, char *argv[])
mesh.update();
# include "checkVolContinuity.H"
# include "volContinuity.H"
# include "meshCourantNo.H"
if

View file

@ -61,7 +61,7 @@ int main(int argc, char *argv[])
mesh.update();
# include "checkVolContinuity.H"
# include "volContinuity.H"
# include "meshCourantNo.H"
if (runTime.timeIndex() % checkFrequency == 0)

View file

@ -1,55 +0,0 @@
if (mesh.moving())
{
// The ddt term constructed by hand because it must be Euler
const scalar dt = runTime.deltaT().value();
const scalarField& V = mesh.V().field();
const scalarField& V0 = mesh.V0().field();
scalar sumV = gSum(V);
scalar sumV0 = gSum(V0);
// Rewrite of mesh motion check based on round-off error problems
// Note precision around V and avoiding division with dt
// HJ, 4/Oct/2011
scalarField volChange = (V - V0)/V;
scalarField divFlux =
- fvc::div(mesh.phi())().internalField()*dt;
scalarField conserve = volChange + divFlux;
scalar sumLocalContErr = dt*gSum(mag(conserve*mesh.V()))/sumV;
scalar globalContErr = dt*gSum(conserve*mesh.V())/sumV;
label nMotionErrors = 0;
scalar maxMotionError = 0;
forAll (conserve, cellI)
{
maxMotionError = Foam::max(maxMotionError, mag(conserve[cellI]));
if (mag(conserve[cellI]) > 1e-8)
{
Info<< "Motion conservation error in cell " << cellI << ": "
<< conserve[cellI]
<< " V: " << V[cellI] << " V0: " << V0[cellI]
<< " volume change: " << volChange[cellI]
<< " div flux: " << divFlux[cellI] << endl;
nMotionErrors++;
}
}
if (nMotionErrors > 0)
{
Info<< "Number of motion errors: " << nMotionErrors
<< " out of nCells = " << mesh.nCells() << endl;
}
Info<< "volume continuity errors : "
<< "volume = " << sumV
<< ", old volume = " << sumV0
<< ", max error = " << maxMotionError
<< ", sum local = " << sumLocalContErr
<< ", global = " << globalContErr << endl;
}

View file

@ -1,26 +1,55 @@
if (mesh.moving())
{
dimensionedScalar one("one", dimless, 1.0);
volScalarField volConservation
(
"volConservation",
-fvc::div(mesh.phi())
);
// The ddt term constructed by hand because it must be Euler
volConservation.internalField() +=
(1.0 - mesh.V0()/mesh.V())/runTime.deltaT().value();
const scalar dt = runTime.deltaT().value();
scalar sumLocalVolContErr = runTime.deltaT().value()*
mag(volConservation)().weightedAverage(mesh.V()).value();
const scalarField& V = mesh.V().field();
const scalarField& V0 = mesh.V0().field();
scalar sumV = gSum(V);
scalar sumV0 = gSum(V0);
scalar globalVolContErr = runTime.deltaT().value()*
volConservation.weightedAverage(mesh.V()).value();
// Rewrite of mesh motion check based on round-off error problems
// Note precision around V and avoiding division with dt
// HJ, 4/Oct/2011
scalarField volChange = (V - V0)/V;
scalarField divFlux =
- fvc::div(mesh.phi())().internalField()*dt;
scalarField conserve = volChange + divFlux;
scalar sumLocalContErr = dt*gSum(mag(conserve*mesh.V()))/sumV;
scalar globalContErr = dt*gSum(conserve*mesh.V())/sumV;
label nMotionErrors = 0;
scalar maxMotionError = 0;
forAll (conserve, cellI)
{
maxMotionError = Foam::max(maxMotionError, mag(conserve[cellI]));
if (mag(conserve[cellI]) > 1e-8)
{
Info<< "Motion conservation error in cell " << cellI << ": "
<< conserve[cellI]
<< " V: " << V[cellI] << " V0: " << V0[cellI]
<< " volume change: " << volChange[cellI]
<< " div flux: " << divFlux[cellI] << endl;
nMotionErrors++;
}
}
if (nMotionErrors > 0)
{
Info<< "Number of motion errors: " << nMotionErrors
<< " out of nCells = " << mesh.nCells() << endl;
}
Info<< "volume continuity errors : "
<< "volume = " << sum(mesh.V()).value()
<< ", max error = " << max(volConservation.internalField())
<< ", sum local = " << sumLocalVolContErr
<< ", global = " << globalVolContErr << endl;
<< "volume = " << sumV
<< ", old volume = " << sumV0
<< ", max error = " << maxMotionError
<< ", sum local = " << sumLocalContErr
<< ", global = " << globalContErr << endl;
}