From 6b09536c330a622da4deb49826efb44f4485ba43 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Fri, 12 Oct 2018 10:37:05 +0100 Subject: [PATCH] Merged volume continuity check files --- .../moveDyMEngineMesh/moveDyMEngineMesh.C | 2 +- .../moveDynamicMesh/moveDynamicMesh.C | 2 +- .../general/include/checkVolContinuity.H | 55 ---------------- .../cfdTools/general/include/volContinuity.H | 65 ++++++++++++++----- 4 files changed, 49 insertions(+), 75 deletions(-) delete mode 100644 src/finiteVolume/cfdTools/general/include/checkVolContinuity.H diff --git a/applications/utilities/mesh/manipulation/moveDyMEngineMesh/moveDyMEngineMesh.C b/applications/utilities/mesh/manipulation/moveDyMEngineMesh/moveDyMEngineMesh.C index 5e45cdca8..bcbd026a8 100644 --- a/applications/utilities/mesh/manipulation/moveDyMEngineMesh/moveDyMEngineMesh.C +++ b/applications/utilities/mesh/manipulation/moveDyMEngineMesh/moveDyMEngineMesh.C @@ -68,7 +68,7 @@ int main(int argc, char *argv[]) mesh.update(); -# include "checkVolContinuity.H" +# include "volContinuity.H" # include "meshCourantNo.H" if diff --git a/applications/utilities/mesh/manipulation/moveDynamicMesh/moveDynamicMesh.C b/applications/utilities/mesh/manipulation/moveDynamicMesh/moveDynamicMesh.C index ee8f78f0d..2b4aae1f1 100644 --- a/applications/utilities/mesh/manipulation/moveDynamicMesh/moveDynamicMesh.C +++ b/applications/utilities/mesh/manipulation/moveDynamicMesh/moveDynamicMesh.C @@ -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) diff --git a/src/finiteVolume/cfdTools/general/include/checkVolContinuity.H b/src/finiteVolume/cfdTools/general/include/checkVolContinuity.H deleted file mode 100644 index c79cbb9a6..000000000 --- a/src/finiteVolume/cfdTools/general/include/checkVolContinuity.H +++ /dev/null @@ -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; -} diff --git a/src/finiteVolume/cfdTools/general/include/volContinuity.H b/src/finiteVolume/cfdTools/general/include/volContinuity.H index f6ff4f0f3..c79cbb9a6 100644 --- a/src/finiteVolume/cfdTools/general/include/volContinuity.H +++ b/src/finiteVolume/cfdTools/general/include/volContinuity.H @@ -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; }