CorrectPhi: Change the divU argument to autoPtr<volScalarField>

If divU is valid the velocity divergence is included in the pcorr equation.
This simplifies the logic in multiphase moveMesh functions supporting
incompressible (with or without mass sources) and compressible fluids.

applications/solvers/lagrangian/denseParticleFoam/correctPhic.H

@@ -9,7 +9,7 @@ CorrectPhi
     Uc,
     p,
     dimensionedScalar("rAUf", dimTime, 1),
-    geometricZeroField(),
+    autoPtr<volScalarField>(),
     pressureReference,
     pimple
 );

applications/solvers/modules/VoFSolver/moveMesh.C

@@ -67,32 +67,16 @@ void Foam::solvers::VoFSolver::moveMesh()
 
                 if (incompressible())
                 {
-                    if (divU.valid())
+                    CorrectPhi
-                    {
+                    (
-                        CorrectPhi
+                        phi,
-                        (
+                        U,
-                            phi,
+                        p_rgh,
-                            U,
+                        surfaceScalarField("rAUf", fvc::interpolate(rAU())),
-                            p_rgh,
+                        divU,
-                            surfaceScalarField("rAUf", fvc::interpolate(rAU())),
+                        pressureReference(),
-                            divU(),
+                        pimple
-                            pressureReference(),
+                    );
-                            pimple
-                        );
-                    }
-                    else
-                    {
-                        CorrectPhi
-                        (
-                            phi,
-                            U,
-                            p_rgh,
-                            surfaceScalarField("rAUf", fvc::interpolate(rAU())),
-                            geometricZeroField(),
-                            pressureReference(),
-                            pimple
-                        );
-                    }
                 }
                 else
                 {

applications/solvers/modules/incompressibleFluid/moveMesh.C

@@ -54,7 +54,7 @@ void Foam::solvers::incompressibleFluid::moveMesh()
                     U,
                     p,
                     dimensionedScalar("rAUf", dimTime, 1),
-                    geometricZeroField(),
+                    autoPtr<volScalarField>(),
                     pressureReference,
                     pimple
                 );

applications/solvers/modules/incompressibleMultiphaseVoF/incompressibleMultiphaseVoF.C

@@ -128,7 +128,7 @@ Foam::solvers::incompressibleMultiphaseVoF::incompressibleMultiphaseVoF
                 U,
                 p_rgh,
                 surfaceScalarField("rAUf", fvc::interpolate(rAU())),
-                geometricZeroField(),
+                autoPtr<volScalarField>(),
                 pressureReference(),
                 pimple
             );
@@ -143,7 +143,7 @@ Foam::solvers::incompressibleMultiphaseVoF::incompressibleMultiphaseVoF
                 U,
                 p_rgh,
                 dimensionedScalar(dimTime/rho.dimensions(), 1),
-                geometricZeroField(),
+                autoPtr<volScalarField>(),
                 pressureReference(),
                 pimple
             );

applications/solvers/modules/incompressibleVoF/incompressibleVoF.C

@@ -113,7 +113,7 @@ Foam::solvers::incompressibleVoF::incompressibleVoF(fvMesh& mesh)
                 U,
                 p_rgh,
                 surfaceScalarField("rAUf", fvc::interpolate(rAU())),
-                geometricZeroField(),
+                autoPtr<volScalarField>(),
                 pressureReference(),
                 pimple
             );
@@ -128,7 +128,7 @@ Foam::solvers::incompressibleVoF::incompressibleVoF(fvMesh& mesh)
                 U,
                 p_rgh,
                 dimensionedScalar(dimTime/rho.dimensions(), 1),
-                geometricZeroField(),
+                autoPtr<volScalarField>(),
                 pressureReference(),
                 pimple
             );

applications/solvers/modules/multiphaseEuler/multiphaseEuler/multiphaseEuler.H

@@ -153,7 +153,7 @@ protected:
         //- Stored divU from the previous mesh so that it can be
         //  mapped and used in correctPhi to ensure the corrected phi
         //  has the same divergence
-        tmp<volScalarField> divU;
+        autoPtr<volScalarField> divU;
 
         //- Read controls
         void readControls();

applications/solvers/modules/multiphaseEuler/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C

@@ -519,10 +519,10 @@ Foam::MovingPhaseModel<BasePhaseModel>::K() const
 
 
 template<class BasePhaseModel>
-Foam::tmp<Foam::volScalarField>
+const Foam::autoPtr<Foam::volScalarField>&
 Foam::MovingPhaseModel<BasePhaseModel>::divU() const
 {
-    return divU_.valid() ? tmp<volScalarField>(divU_()) : tmp<volScalarField>();
+    return divU_;
 }
 
 
@@ -531,13 +531,13 @@ void Foam::MovingPhaseModel<BasePhaseModel>::divU(tmp<volScalarField> divU)
 {
     if (!divU_.valid())
     {
-        divU_ = divU;
+        divU_ = divU.ptr();
-        divU_.ref().rename(IOobject::groupName("divU", this->name()));
+        divU_().rename(IOobject::groupName("divU", this->name()));
-        divU_.ref().checkIn();
+        divU_().checkIn();
     }
     else
     {
-        divU_.ref() = divU;
+        divU_() = divU;
     }
 }
 

applications/solvers/modules/multiphaseEuler/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H

@@ -119,7 +119,7 @@ protected:
         mutable tmp<surfaceScalarField> DUDtf_;
 
         //- Dilatation rate
-        tmp<volScalarField> divU_;
+        autoPtr<volScalarField> divU_;
 
         //- Turbulence model
         autoPtr<phaseCompressible::momentumTransportModel> momentumTransport_;
@@ -253,7 +253,7 @@ public:
         // Compressibility (variable density)
 
             //- Return the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
-            virtual tmp<volScalarField> divU() const;
+            virtual const autoPtr<volScalarField>& divU() const;
 
             //- Set the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
             virtual void divU(tmp<volScalarField> divU);

applications/solvers/modules/multiphaseEuler/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C

@@ -259,10 +259,11 @@ Foam::StationaryPhaseModel<BasePhaseModel>::K() const
 
 
 template<class BasePhaseModel>
-Foam::tmp<Foam::volScalarField>
+const Foam::autoPtr<Foam::volScalarField>&
 Foam::StationaryPhaseModel<BasePhaseModel>::divU() const
 {
-    return tmp<volScalarField>();
+    static autoPtr<volScalarField> divU_;
+    return divU_;
 }
 
 

applications/solvers/modules/multiphaseEuler/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.H

@@ -136,7 +136,7 @@ public:
         // Compressibility (variable density)
 
             //- Return the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
-            virtual tmp<volScalarField> divU() const;
+            virtual const autoPtr<volScalarField>& divU() const;
 
             //- Set the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
             virtual void divU(tmp<volScalarField> divU);

applications/solvers/modules/multiphaseEuler/phaseSystems/phaseModel/phaseModel/phaseModel.H

@@ -246,7 +246,7 @@ public:
             virtual bool isochoric() const = 0;
 
             //- Return the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
-            virtual tmp<volScalarField> divU() const = 0;
+            virtual const autoPtr<volScalarField>& divU() const = 0;
 
             //- Set the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
             virtual void divU(tmp<volScalarField> divU) = 0;

applications/solvers/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.C

@@ -728,7 +728,7 @@ void Foam::phaseSystem::correctBoundaryFlux()
 void Foam::phaseSystem::correctPhi
 (
     const volScalarField& p_rgh,
-    const tmp<volScalarField>& divU,
+    const autoPtr<volScalarField>& divU,
     const pressureReference& pressureReference,
     nonOrthogonalSolutionControl& pimple
 )
@@ -778,32 +778,16 @@ void Foam::phaseSystem::correctPhi
 
     if (incompressible())
     {
-        if (divU.valid())
+        CorrectPhi
-        {
+        (
-            CorrectPhi
+            phi_,
-            (
+            movingPhases()[0].U(),
-                phi_,
+            p_rgh,
-                movingPhases()[0].U(),
+            dimensionedScalar(dimTime/dimDensity, 1),
-                p_rgh,
+            divU,
-                dimensionedScalar(dimTime/dimDensity, 1),
+            pressureReference,
-                divU(),
+            pimple
-                pressureReference,
+        );
-                pimple
-            );
-        }
-        else
-        {
-            CorrectPhi
-            (
-                phi_,
-                movingPhases()[0].U(),
-                p_rgh,
-                dimensionedScalar(dimTime/dimDensity, 1),
-                geometricZeroField(),
-                pressureReference,
-                pimple
-            );
-        }
     }
     else
     {

applications/solvers/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.H

@@ -656,7 +656,7 @@ public:
             void correctPhi
             (
                 const volScalarField& p_rgh,
-                const tmp<volScalarField>& divU,
+                const autoPtr<volScalarField>& divU,
                 const pressureReference& pressureReference,
                 nonOrthogonalSolutionControl& pimple
             );

src/finiteVolume/cfdTools/general/CorrectPhi/CorrectPhi.C

@@ -39,14 +39,14 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-template<class RAUfType, class DivUType>
+template<class RAUfType>
 void Foam::CorrectPhi
 (
     surfaceScalarField& phi,
     const volVectorField& U,
     const volScalarField& p,
     const RAUfType& rAUf,
-    const DivUType& divU,
+    const autoPtr<volScalarField>& divU,
     const pressureReference& pressureReference,
     nonOrthogonalSolutionControl& pcorrControl
 )
@@ -98,7 +98,13 @@ void Foam::CorrectPhi
         // matches the divU provided (from previous iteration, time-step...)
         fvScalarMatrix pcorrEqn
         (
-            fvm::laplacian(rAUf, pcorr) == fvc::div(phi) - divU
+            fvm::laplacian(rAUf, pcorr)
+         ==
+            (
+                divU.valid()
+              ? fvc::div(phi) - divU()
+              : fvc::div(phi)
+            )
         );
 
         pcorrEqn.setReference(pressureReference.refCell(), 0);
@@ -113,14 +119,14 @@ void Foam::CorrectPhi
 }
 
 
-template<class RAUfType, class DivRhoUType>
+template<class RAUfType>
 void Foam::CorrectPhi
 (
     surfaceScalarField& phi,
     const volScalarField& p,
     const volScalarField& psi,
     const RAUfType& rAUf,
-    const DivRhoUType& divRhoU,
+    const volScalarField& divRhoU,
     nonOrthogonalSolutionControl& pcorrControl
 )
 {

src/finiteVolume/cfdTools/general/CorrectPhi/CorrectPhi.H

@@ -70,26 +70,26 @@ namespace Foam
         const bool evaluateUBCs
     );
 
-    template<class RAUfType, class DivUType>
+    template<class RAUfType>
     void CorrectPhi
     (
         surfaceScalarField& phi,
         const volVectorField& U,
         const volScalarField& p,
         const RAUfType& rAUf,
-        const DivUType& divU,
+        const autoPtr<volScalarField>& divU,
         const pressureReference& pressureReference,
         nonOrthogonalSolutionControl& pcorrControl
     );
 
-    template<class RAUfType, class DivRhoUType>
+    template<class RAUfType>
     void CorrectPhi
     (
         surfaceScalarField& phi,
         const volScalarField& p,
         const volScalarField& psi,
         const RAUfType& rAUf,
-        const DivRhoUType& divRhoU,
+        const volScalarField& divRhoU,
         nonOrthogonalSolutionControl& pcorrControl
     );
 }