exaFOAM/foam-extend4.1-coherent-io
Archived
1
1
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

FEATURE: Improvements in the coupled solver coupling syntax. First version of the coupled compressible solver. Author: Hrvoje Jasak. Merge: Dominik Christ.

Browse source
This commit is contained in:
Dominik Christ 2014-09-30 16:25:28 +01:00
commit 98114b178a
26 changed files with 345 additions and 484 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

15
applications/solvers/coupled/pUCoupledFoam/couplingTerms.H Normal file
View file

@ -0,0 +1,15 @@
{
// Calculate grad p coupling matrix. Needs to be here if one uses
// gradient schemes with limiters. VV, 9/June/2014
BlockLduSystem<vector, vector> pInU(fvm::grad(p));
// Calculate div U coupling. Could be calculated only once since
// it is only geometry dependent. VV, 9/June/2014
BlockLduSystem<vector, scalar> UInp(fvm::UDiv(U));
// Last argument in insertBlockCoupling says if the column direction
// should be incremented. This is needed for arbitrary positioning
// of U and p in the system. This could be better. VV, 30/April/2014
UpEqn.insertBlockCoupling(0, 3, pInU, true);
UpEqn.insertBlockCoupling(3, 0, UInp, false);
}

16
applications/solvers/coupled/pUCoupledFoam/pUCoupledFoam.C
View file

@ -38,9 +38,9 @@ Authors
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "fvBlockMatrix.H"
#include "singlePhaseTransportModel.H"
#include "RASModel.H"
#include "fvBlockMatrix.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -54,9 +54,6 @@ int main(int argc, char *argv[])
# include "initContinuityErrs.H"
# include "readBlockSolverControls.H"
// Calculate div U coupling only once since it is only geometry dependant.
BlockLduSystem<vector, scalar> UInp(fvm::div(U));
Info<< "\nStarting time loop\n" << endl;
while (runTime.loop())
{
@ -73,15 +70,8 @@ int main(int argc, char *argv[])
// Assemble and insert pressure equation
# include "pEqn.H"
// Calculate grad p coupling matrix. Needs to be here if one uses
// gradient schemes with limiters. VV, 9/June/2014
BlockLduSystem<vector, vector> pInU(fvm::grad(p));
// Last argument in insertBlockCoupling says if the column direction
// should be incremented. This is needed for arbitrary positioning
// of U and p in the system. This could be better. VV, 30/April/2014
UpEqn.insertBlockCoupling(0, 3, pInU, true);
UpEqn.insertBlockCoupling(3, 0, UInp, false);
// Assemble and insert coupling terms
# include "couplingTerms.H"
// Solve the block matrix
UpEqn.solve();

3
applications/solvers/coupled/pUCoupledFullPicardFoam/Make/files
View file

@ -1,3 +0,0 @@
pUCoupledFullPicardFoam.C
EXE = $(FOAM_APPBIN)/pUCoupledFullPicardFoam

12
applications/solvers/coupled/pUCoupledFullPicardFoam/Make/options
View file

@ -1,12 +0,0 @@
EXE_INC = \
-I$(LIB_SRC)/turbulenceModels \
-I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-lincompressibleRASModels \
-lincompressibleTransportModels \
-lfiniteVolume \
-llduSolvers

12
applications/solvers/coupled/pUCoupledFullPicardFoam/UEqn.H
View file

@ -1,12 +0,0 @@
// Momentum equation
fvVectorMatrix UEqn
(
fvm::div(phi, U)
+ fvc::div(-phi, U, "div(phi,U)")
+ turbulence->divDevReff(U)
);
UEqn.relax();
UpEqn.insertEquation(0, UEqn);
UpEqn.insertPicardTensor(0, U, phi);

52
applications/solvers/coupled/pUCoupledFullPicardFoam/createFields.H
View file

@ -1,52 +0,0 @@
Info << "Reading field p\n" << endl;
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info << "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
#include "createPhi.H"
singlePhaseTransportModel laminarTransport(U, phi);
autoPtr<incompressible::RASModel> turbulence
(
incompressible::RASModel::New(U, phi, laminarTransport)
);
// Block vector field for velocity (first entry) and pressure (second
// entry).
Info << "Creating field Up\n" << endl;
volVector4Field Up
(
IOobject
(
"Up",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedVector4("zero", dimless, vector4::zero)
);

23
applications/solvers/coupled/pUCoupledFullPicardFoam/pEqn.H
View file

@ -1,23 +0,0 @@
// Pressure parts of the continuity equation
surfaceScalarField rUAf
(
"rUAf",
fvc::interpolate(1.0/UEqn.A())
);
surfaceScalarField presSource
(
"presSource",
rUAf*(fvc::interpolate(fvc::grad(p, "grad(pSource)")) & mesh.Sf())
);
fvScalarMatrix pEqn
(
- fvm::laplacian(rUAf, p)
==
- fvc::div(presSource)
);
pEqn.setReference(pRefCell, pRefValue);
UpEqn.insertEquation(3, pEqn);

116
applications/solvers/coupled/pUCoupledFullPicardFoam/pUCoupledFullPicardFoam.C
View file

@ -1,116 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Application
pUCoupledFullPicardFoam
Description
Steady-state solver for incompressible, turbulent flow, with implicit
coupling between pressure and velocity achieved by fvBlockMatrix.
Turbulence is in this version solved using the existing turbulence
structure.
Convective term is properly linearised using Picard linearisation and
additional block coefficients. This is full, proper variant!
VV, 23/July/2014.
Authors
Klas Jareteg, Chalmers University of Technology,
Vuko Vukcevic, FMENA Zagreb.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "RASModel.H"
#include "fvBlockMatrix.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
# include "setRootCase.H"
# include "createTime.H"
# include "createMesh.H"
# include "createFields.H"
# include "initContinuityErrs.H"
# include "readBlockSolverControls.H"
// Calculate div U coupling only once since it is only geometry dependant.
BlockLduSystem<vector, scalar> UInp(fvm::div(U));
Info<< "\nStarting time loop\n" << endl;
while (runTime.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;
p.storePrevIter();
// Initialize the Up block system (matrix, source and reference to Up)
fvBlockMatrix<vector4> UpEqn(Up);
// Assemble and insert momentum equation
# include "UEqn.H"
// Assemble and insert pressure equation
# include "pEqn.H"
// Calculate grad p coupling matrix. Needs to be here if one uses
// gradient schemes with limiters. VV, 9/June/2014
BlockLduSystem<vector, vector> pInU(fvm::grad(p));
// Last argument in insertBlockCoupling says if the column direction
// should be incremented. This is needed for arbitrary positioning
// of U and p in the system. This could be better. VV, 30/April/2014
UpEqn.insertBlockCoupling(0, 3, pInU, true);
UpEqn.insertBlockCoupling(3, 0, UInp, false);
// Solve the block matrix
UpEqn.solve();
// Retrieve solution
UpEqn.retrieveSolution(0, U.internalField());
UpEqn.retrieveSolution(3, p.internalField());
U.correctBoundaryConditions();
p.correctBoundaryConditions();
phi = (fvc::interpolate(U) & mesh.Sf()) + pEqn.flux() + presSource;
# include "continuityErrs.H"
p.relax();
turbulence->correct();
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}

3
applications/solvers/coupled/pUCoupledFullPicardFoam/readBlockSolverControls.H
View file

@ -1,3 +0,0 @@
label pRefCell = 0;
scalar pRefValue = 0;
setRefCell(p, mesh.solutionDict().subDict("blockSolver"), pRefCell, pRefValue);

3
applications/solvers/coupled/pUCoupledSemiPicardFoam/Make/files
View file

@ -1,3 +0,0 @@
pUCoupledSemiPicardFoam.C
EXE = $(FOAM_APPBIN)/pUCoupledSemiPicardFoam

12
applications/solvers/coupled/pUCoupledSemiPicardFoam/Make/options
View file

@ -1,12 +0,0 @@
EXE_INC = \
-I$(LIB_SRC)/turbulenceModels \
-I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-lincompressibleRASModels \
-lincompressibleTransportModels \
-lfiniteVolume \
-llduSolvers

11
applications/solvers/coupled/pUCoupledSemiPicardFoam/UEqn.H
View file

@ -1,11 +0,0 @@
// Momentum equation
fvVectorMatrix UEqn
(
fvm::div(2*phi, U,"div(phi,U)")
+ fvc::div(-phi, U,"div(phi,U)")
+ turbulence->divDevReff(U)
);
UEqn.relax();
UpEqn.insertEquation(0, UEqn);

52
applications/solvers/coupled/pUCoupledSemiPicardFoam/createFields.H
View file

@ -1,52 +0,0 @@
Info << "Reading field p\n" << endl;
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info << "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
#include "createPhi.H"
singlePhaseTransportModel laminarTransport(U, phi);
autoPtr<incompressible::RASModel> turbulence
(
incompressible::RASModel::New(U, phi, laminarTransport)
);
// Block vector field for velocity (first entry) and pressure (second
// entry).
Info << "Creating field Up\n" << endl;
volVector4Field Up
(
IOobject
(
"Up",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedVector4("zero", dimless, vector4::zero)
);

23
applications/solvers/coupled/pUCoupledSemiPicardFoam/pEqn.H
View file

@ -1,23 +0,0 @@
// Pressure parts of the continuity equation
surfaceScalarField rUAf
(
"rUAf",
fvc::interpolate(1.0/UEqn.A())
);
surfaceScalarField presSource
(
"presSource",
rUAf*(fvc::interpolate(fvc::grad(p, "grad(pSource)")) & mesh.Sf())
);
fvScalarMatrix pEqn
(
- fvm::laplacian(rUAf, p)
==
- fvc::div(presSource)
);
pEqn.setReference(pRefCell, pRefValue);
UpEqn.insertEquation(3, pEqn);

116
applications/solvers/coupled/pUCoupledSemiPicardFoam/pUCoupledSemiPicardFoam.C
View file

@ -1,116 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Application
pUCoupledSemiPicardFoam
Description
Steady-state solver for incompressible, turbulent flow, with implicit
coupling between pressure and velocity achieved by fvBlockMatrix.
Turbulence is in this version solved using the existing turbulence
structure.
Convective term is properly linearised using Picard linearisation and
additional block coefficients. This is semi variant where the difference
between succesive velocity fields is neglected! VV, 23/July/2014.
Authors
Klas Jareteg, Chalmers University of Technology,
Vuko Vukcevic, FMENA Zagreb.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "RASModel.H"
#include "fvBlockMatrix.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
# include "setRootCase.H"
# include "createTime.H"
# include "createMesh.H"
# include "createFields.H"
# include "initContinuityErrs.H"
# include "readBlockSolverControls.H"
// Calculate div U coupling only once since it is only geometry dependant.
BlockLduSystem<vector, scalar> UInp(fvm::div(U));
Info<< "\nStarting time loop\n" << endl;
while (runTime.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;
p.storePrevIter();
// Initialize the Up block system (matrix, source and reference to Up)
fvBlockMatrix<vector4> UpEqn(Up);
// Assemble and insert momentum equation
# include "UEqn.H"
// Assemble and insert pressure equation
# include "pEqn.H"
// Calculate grad p coupling matrix. Needs to be here if one uses
// gradient schemes with limiters. VV, 9/June/2014
BlockLduSystem<vector, vector> pInU(fvm::grad(p));
// Last argument in insertBlockCoupling says if the column direction
// should be incremented. This is needed for arbitrary positioning
// of U and p in the system. This could be better. VV, 30/April/2014
UpEqn.insertBlockCoupling(0, 3, pInU, true);
UpEqn.insertBlockCoupling(3, 0, UInp, false);
// Solve the block matrix
UpEqn.solve();
// Retrieve solution
UpEqn.retrieveSolution(0, U.internalField());
UpEqn.retrieveSolution(3, p.internalField());
U.correctBoundaryConditions();
p.correctBoundaryConditions();
phi = (fvc::interpolate(U) & mesh.Sf()) + pEqn.flux() + presSource;
# include "continuityErrs.H"
p.relax();
turbulence->correct();
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}

3
applications/solvers/coupled/pUCoupledSemiPicardFoam/readBlockSolverControls.H
View file

@ -1,3 +0,0 @@
label pRefCell = 0;
scalar pRefValue = 0;
setRefCell(p, mesh.solutionDict().subDict("blockSolver"), pRefCell, pRefValue);

34
src/finiteVolume/finiteVolume/divSchemes/divScheme/divScheme.C
View file

@ -103,7 +103,7 @@ tmp
<
BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
>
divScheme<Type>::fvmDiv
divScheme<Type>::fvmUDiv
(
const GeometricField<Type, fvPatchField, volMesh>& vf
) const
@ -128,6 +128,38 @@ divScheme<Type>::fvmDiv
}
template<class Type>
tmp
<
BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
>
divScheme<Type>::fvmUDiv
(
const surfaceScalarField& flux,
const GeometricField<Type, fvPatchField, volMesh>& vf
) const
{
FatalErrorIn
(
"tmp<BlockLduSystem> divScheme<Type>::fvmDiv\n"
"(\n"
" surfaceScalarField&"
" GeometricField<Type, fvPatchField, volMesh>&"
")\n"
) << "Implicit div operator currently defined only for Gauss linear. "
<< abort(FatalError);
typedef typename innerProduct<vector, Type>::type DivType;
tmp<BlockLduSystem<vector, DivType> > tbs
(
new BlockLduSystem<vector, DivType>(vf.mesh())
);
return tbs;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace fv

11
src/finiteVolume/finiteVolume/divSchemes/divScheme/divScheme.H
View file

@ -157,10 +157,19 @@ public:
virtual tmp
<
BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
> fvmDiv
> fvmUDiv
(
const GeometricField<Type, fvPatchField, volMesh>&
) const;
virtual tmp
<
BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
> fvmUDiv
(
const surfaceScalarField&,
const GeometricField<Type, fvPatchField, volMesh>&
) const;
};

35
src/finiteVolume/finiteVolume/divSchemes/gaussDivScheme/gaussDivScheme.C
View file

@ -72,14 +72,14 @@ template<class Type>
tmp
<
BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
> gaussDivScheme<Type>::fvmDiv
> gaussDivScheme<Type>::fvmUDiv
(
const GeometricField<Type, fvPatchField, volMesh>& vf
) const
{
FatalErrorIn
(
"tmp<BlockLduSystem> fvmDiv\n"
"tmp<BlockLduSystem> fvmUDiv\n"
"(\n"
" GeometricField<Type, fvPatchField, volMesh>&"
")\n"
@ -97,6 +97,37 @@ tmp
}
template<class Type>
tmp
<
BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
> gaussDivScheme<Type>::fvmUDiv
(
const surfaceScalarField& flux,
const GeometricField<Type, fvPatchField, volMesh>& vf
) const
{
FatalErrorIn
(
"tmp<BlockLduSystem> fvmUDiv\n"
"(\n"
" const surfaceScalarField& flux"
" const GeometricField<Type, fvPatchField, volMesh>&"
")\n"
) << "Implicit div operator defined only for vector."
<< abort(FatalError);
typedef typename innerProduct<vector, Type>::type DivType;
tmp<BlockLduSystem<vector, DivType> > tbs
(
new BlockLduSystem<vector, DivType>(vf.mesh())
);
return tbs;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace fv

11
src/finiteVolume/finiteVolume/divSchemes/gaussDivScheme/gaussDivScheme.H
View file

@ -102,10 +102,19 @@ public:
tmp
<
BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
> fvmDiv
> fvmUDiv
(
const GeometricField<Type, fvPatchField, volMesh>&
) const;
tmp
<
BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
> fvmUDiv
(
const surfaceScalarField& flux,
const GeometricField<Type, fvPatchField, volMesh>&
) const;
};

101
src/finiteVolume/finiteVolume/divSchemes/vectorGaussDivScheme/vectorGaussDivScheme.C
View file

@ -42,7 +42,7 @@ namespace fv
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<>
tmp<BlockLduSystem<vector, scalar> > gaussDivScheme<vector>::fvmDiv
tmp<BlockLduSystem<vector, scalar> > gaussDivScheme<vector>::fvmUDiv
(
const GeometricField<vector, fvPatchField, volMesh>& vf
) const
@ -60,9 +60,9 @@ tmp<BlockLduSystem<vector, scalar> > gaussDivScheme<vector>::fvmDiv
scalarField& source = bs.source();
// Grab ldu parts of block matrix as linear always
typename CoeffField<vector>::linearTypeField& d = bs.diag().asLinear();
typename CoeffField<vector>::linearTypeField& u = bs.upper().asLinear();
typename CoeffField<vector>::linearTypeField& l = bs.lower().asLinear();
CoeffField<vector>::linearTypeField& d = bs.diag().asLinear();
CoeffField<vector>::linearTypeField& u = bs.upper().asLinear();
CoeffField<vector>::linearTypeField& l = bs.lower().asLinear();
const vectorField& SfIn = mesh.Sf().internalField();
@ -89,9 +89,9 @@ tmp<BlockLduSystem<vector, scalar> > gaussDivScheme<vector>::fvmDiv
if (patch.coupled())
{
typename CoeffField<vector>::linearTypeField& pcoupleUpper =
CoeffField<vector>::linearTypeField& pcoupleUpper =
bs.coupleUpper()[patchI].asLinear();
typename CoeffField<vector>::linearTypeField& pcoupleLower =
CoeffField<vector>::linearTypeField& pcoupleLower =
bs.coupleLower()[patchI].asLinear();
const vectorField pcl = -pw*Sf;
@ -119,6 +119,95 @@ tmp<BlockLduSystem<vector, scalar> > gaussDivScheme<vector>::fvmDiv
}
template<>
tmp<BlockLduSystem<vector, scalar> > gaussDivScheme<vector>::fvmUDiv
(
const surfaceScalarField& flux,
const GeometricField<vector, fvPatchField, volMesh>& vf
) const
{
tmp<surfaceScalarField> tweights = this->tinterpScheme_().weights(vf);
const scalarField& wIn = tweights().internalField();
const fvMesh& mesh = vf.mesh();
tmp<BlockLduSystem<vector, scalar> > tbs
(
new BlockLduSystem<vector, scalar>(mesh)
);
BlockLduSystem<vector, scalar>& bs = tbs();
scalarField& source = bs.source();
// Grab ldu parts of block matrix as linear always
CoeffField<vector>::linearTypeField& d = bs.diag().asLinear();
CoeffField<vector>::linearTypeField& u = bs.upper().asLinear();
CoeffField<vector>::linearTypeField& l = bs.lower().asLinear();
const vectorField& SfIn = mesh.Sf().internalField();
const scalarField& fluxIn = flux.internalField();
l = -wIn*fluxIn*SfIn;
u = l + fluxIn*SfIn;
bs.negSumDiag();
// Boundary contributions
forAll(vf.boundaryField(), patchI)
{
const fvPatchVectorField& pf = vf.boundaryField()[patchI];
const fvPatch& patch = pf.patch();
const vectorField& Sf = patch.Sf();
const fvsPatchScalarField& pw = tweights().boundaryField()[patchI];
const unallocLabelList& fc = patch.faceCells();
const scalarField& pFlux = flux.boundaryField()[patchI];
const vectorField internalCoeffs(pf.valueInternalCoeffs(pw));
// Diag contribution
forAll(pf, faceI)
{
d[fc[faceI]] += cmptMultiply
(
internalCoeffs[faceI],
pFlux[faceI]*Sf[faceI]
);
}
if (patch.coupled())
{
CoeffField<vector>::linearTypeField& pcoupleUpper =
bs.coupleUpper()[patchI].asLinear();
CoeffField<vector>::linearTypeField& pcoupleLower =
bs.coupleLower()[patchI].asLinear();
const vectorField pcl = -pw*pFlux*Sf;
const vectorField pcu = pcl + pFlux*Sf;
// Coupling contributions
pcoupleLower -= pcl;
pcoupleUpper -= pcu;
}
else
{
const vectorField boundaryCoeffs(pf.valueBoundaryCoeffs(pw));
// Boundary contribution
forAll(pf, faceI)
{
source[fc[faceI]] -=
(
boundaryCoeffs[faceI] & (pFlux[faceI]*Sf[faceI])
);
}
}
}
// Interpolation schemes with corrections not accounted for
return tbs;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace fv

9
src/finiteVolume/finiteVolume/divSchemes/vectorGaussDivScheme/vectorGaussDivScheme.H
View file

@ -51,11 +51,18 @@ namespace fv
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<>
tmp<BlockLduSystem<vector, scalar> > gaussDivScheme<vector>::fvmDiv
tmp<BlockLduSystem<vector, scalar> > gaussDivScheme<vector>::fvmUDiv
(
const GeometricField<vector, fvPatchField, volMesh>&
) const;
template<>
tmp<BlockLduSystem<vector, scalar> > gaussDivScheme<vector>::fvmUDiv
(
const surfaceScalarField& flux,
const GeometricField<vector, fvPatchField, volMesh>&
) const;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

87
src/finiteVolume/finiteVolume/fvm/fvmDiv.C
View file

@ -103,7 +103,7 @@ template<class Type>
tmp
<
BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
> div
> UDiv
(
GeometricField<Type, fvPatchField, volMesh>& vf,
const word& name
@ -113,7 +113,7 @@ tmp
(
vf.mesh(),
vf.mesh().schemesDict().divScheme(name)
)().fvmDiv(vf);
)().fvmUDiv(vf);
}
@ -121,12 +121,52 @@ template<class Type>
tmp
<
BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
> div
> UDiv
(
const surfaceScalarField& flux,
GeometricField<Type, fvPatchField, volMesh>& vf,
const word& name
)
{
return fv::divScheme<Type>::New
(
vf.mesh(),
vf.mesh().schemesDict().divScheme(name)
)().fvmUDiv(flux, vf);
}
template<class Type>
tmp
<
BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
> UDiv
(
const tmp<surfaceScalarField>& tflux,
GeometricField<Type, fvPatchField, volMesh>& vf,
const word& name
)
{
tmp
<
BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
>
Div(fvm::UDiv(tflux(), vf, name));
tflux.clear();
return Div;
}
template<class Type>
tmp
<
BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
> UDiv
(
GeometricField<Type, fvPatchField, volMesh>& vf
)
{
return fvm::div
return fvm::UDiv
(
vf,
"div(" + vf.name() + ')'
@ -134,6 +174,45 @@ tmp
}
template<class Type>
tmp
<
BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
> UDiv
(
const surfaceScalarField& flux,
GeometricField<Type, fvPatchField, volMesh>& vf
)
{
return fvm::UDiv
(
flux,
vf,
"div(" + vf.name() + ')'
);
}
template<class Type>
tmp
<
BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
> UDiv
(
const tmp<surfaceScalarField>& tflux,
GeometricField<Type, fvPatchField, volMesh>& vf
)
{
tmp
<
BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
>
Div(fvm::UDiv(tflux(), vf));
tflux.clear();
return Div;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace fvm

49
src/finiteVolume/finiteVolume/fvm/fvmDiv.H
View file

@ -68,7 +68,6 @@ namespace fvm
const word& name
);
template<class Type>
tmp<fvMatrix<Type> > div
(
@ -83,12 +82,12 @@ namespace fvm
GeometricField<Type, fvPatchField, volMesh>&
);
// Implicit div operator for block systems
// Implicit div operators for block systems
template<class Type>
tmp
<
BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
> div
> UDiv
(
GeometricField<Type, fvPatchField, volMesh>&,
const word&
@ -98,8 +97,50 @@ namespace fvm
tmp
<
BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
> div
> UDiv
(
const surfaceScalarField&,
GeometricField<Type, fvPatchField, volMesh>&,
const word&
);
template<class Type>
tmp
<
BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
> UDiv
(
const tmp<surfaceScalarField>&,
GeometricField<Type, fvPatchField, volMesh>&,
const word&
);
template<class Type>
tmp
<
BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
> UDiv
(
GeometricField<Type, fvPatchField, volMesh>&
);
template<class Type>
tmp
<
BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
> UDiv
(
const surfaceScalarField&,
GeometricField<Type, fvPatchField, volMesh>&
);
template<class Type>
tmp
<
BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
> UDiv
(
const tmp<surfaceScalarField>&,
GeometricField<Type, fvPatchField, volMesh>&
);
}

10
src/finiteVolume/finiteVolume/gradSchemes/scalarGaussGrad/scalarGaussGrad.C
View file

@ -60,9 +60,9 @@ tmp<BlockLduSystem<vector, vector> > gaussGrad<scalar>::fvmGrad
vectorField& source = bs.source();
// Grab ldu parts of block matrix as linear always
typename CoeffField<vector>::linearTypeField& d = bs.diag().asLinear();
typename CoeffField<vector>::linearTypeField& u = bs.upper().asLinear();
typename CoeffField<vector>::linearTypeField& l = bs.lower().asLinear();
CoeffField<vector>::linearTypeField& d = bs.diag().asLinear();
CoeffField<vector>::linearTypeField& u = bs.upper().asLinear();
CoeffField<vector>::linearTypeField& l = bs.lower().asLinear();
const vectorField& SfIn = mesh.Sf().internalField();
@ -89,9 +89,9 @@ tmp<BlockLduSystem<vector, vector> > gaussGrad<scalar>::fvmGrad
if (patch.coupled())
{
typename CoeffField<vector>::linearTypeField& pcoupleUpper =
CoeffField<vector>::linearTypeField& pcoupleUpper =
bs.coupleUpper()[patchI].asLinear();
typename CoeffField<vector>::linearTypeField& pcoupleLower =
CoeffField<vector>::linearTypeField& pcoupleLower =
bs.coupleLower()[patchI].asLinear();
const vectorField pcl = -pw*Sf;

10
src/finiteVolume/finiteVolume/gradSchemes/scalarLeastSquaresGrad/scalarLeastSquaresGrad.C
View file

@ -80,9 +80,9 @@ tmp<BlockLduSystem<vector, vector> > leastSquaresGrad<scalar>::fvmGrad
vectorField& source = bs.source();
// Grab ldu parts of block matrix as linear always
typename CoeffField<vector>::linearTypeField& d = bs.diag().asLinear();
typename CoeffField<vector>::linearTypeField& u = bs.upper().asLinear();
typename CoeffField<vector>::linearTypeField& l = bs.lower().asLinear();
CoeffField<vector>::linearTypeField& d = bs.diag().asLinear();
CoeffField<vector>::linearTypeField& u = bs.upper().asLinear();
CoeffField<vector>::linearTypeField& l = bs.lower().asLinear();
// Get references to least square vectors
const leastSquaresVectors& lsv = leastSquaresVectors::New(mesh);
@ -125,9 +125,9 @@ tmp<BlockLduSystem<vector, vector> > leastSquaresGrad<scalar>::fvmGrad
const scalarField cellVInNei =
cellV.boundaryField()[patchI].patchNeighbourField();
typename CoeffField<vector>::linearTypeField& pcoupleUpper =
CoeffField<vector>::linearTypeField& pcoupleUpper =
bs.coupleUpper()[patchI].asLinear();
typename CoeffField<vector>::linearTypeField& pcoupleLower =
CoeffField<vector>::linearTypeField& pcoupleLower =
bs.coupleLower()[patchI].asLinear();
// Coupling and diagonal contributions