Skip to content

Commit

Permalink
Conservative GGI bridging updates and fixes
Browse files Browse the repository at this point in the history 
Scaling the fluxes and killing non-orthogonal correction vectors in case
bridging is used, both in order to ensure global conservation across bridged GGI
with potentially partially overlapping faces.
  • Loading branch information
Vuko Vukcevic committed Mar 15, 2018
1 parent 412eaba commit 88cdcdc
Show file tree
Hide file tree
Showing 6 changed files with 101 additions and 74 deletions.
1 change: 1 addition & 0 deletions src/finiteVolume/Make/files
View file
Original file line number Diff line number Diff line change
Expand Up @@ -119,6 +119,7 @@ $(constraintFvPatchFields)/symmetry/symmetryFvPatchFields.C
$(constraintFvPatchFields)/wedge/wedgeFvPatchFields.C
$(constraintFvPatchFields)/wedge/wedgeFvPatchScalarField.C
$(constraintFvPatchFields)/wedge/wedgeFvPatchVectorNFields.C
$(constraintFvPatchFields)/ggi/ggiFvPatchScalarField.C
$(constraintFvPatchFields)/ggi/ggiFvPatchFields.C
$(constraintFvPatchFields)/ggi/ggiFvPatchVectorNFields.C
$(constraintFvPatchFields)/jumpGgi/jumpGgiFvPatchFields.C
Expand Down
25 changes: 1 addition & 24 deletions src/finiteVolume/fields/fvPatchFields/constraint/ggi/ggiFvPatchField.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,7 @@ Note on parallelisation
#include "ggiFvPatchField.H"
#include "symmTransformField.H"
#include "coeffFields.H"
#include "fvMatrices.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

Expand Down Expand Up @@ -520,12 +521,6 @@ void Foam::ggiFvPatchField<Type>::manipulateValueCoeffs

// Scale partially overlapping boundary coeffs to ensure conservation
ggiPatch_.shadow().scalePartialFaces(slaveBC);

// Info<< "MANIPULATING VALUE COEFFS" << endl;
// Info<< "Slave internal coeffs: " << slaveIC << nl
// << "Slave boundary coeffs: " << slaveBC << nl
// << "Master internal coeffs: " << masterIC << nl
// << "Master boundary coeffs: " << masterBC << endl;
}
}

Expand All @@ -552,12 +547,6 @@ void Foam::ggiFvPatchField<Type>::manipulateGradientCoeffs
Field<Type>& slaveIC = matrix.internalCoeffs()[sPatchI];
Field<Type>& slaveBC = matrix.boundaryCoeffs()[sPatchI];

// Info<< "ORIGINAL COEFFS" << endl;
// Info<< "Slave internal coeffs: " << slaveIC << nl
// << "Slave boundary coeffs: " << slaveBC << nl
// << "Master internal coeffs: " << masterIC << nl
// << "Master boundary coeffs: " << masterBC << endl;

// Get surface area magnitudes
const scalarField& magSfMaster = ggiPatch_.magSf();
const scalarField& magSfSlave = ggiPatch_.shadow().magSf();
Expand All @@ -579,12 +568,6 @@ void Foam::ggiFvPatchField<Type>::manipulateGradientCoeffs
const Field<Type> slaveInterpolatedBC =
magSfMaster*ggiPatch_.interpolate(slaveBC/magSfSlave);

// Info<< "INTERPOLATED COEFFS" << endl;
// Info<< "Slave internal coeffs: " << slaveInterpolatedIC << nl
// << "Slave boundary coeffs: " << slaveInterpolatedBC << nl
// << "Master internal coeffs: " << masterInterpolatedIC << nl
// << "Master boundary coeffs: " << masterInterpolatedBC << endl;


// Set partially covered master coeffs using slave data
ggiPatch_.setPartialFaces(slaveInterpolatedBC, masterIC);
Expand All @@ -610,12 +593,6 @@ void Foam::ggiFvPatchField<Type>::manipulateGradientCoeffs
// Add to partially overlapping slave internal coeffs to ensure
// conservation. Note: boundary coeffs must be negated
ggiPatch_.shadow().addToPartialFaces((-slaveBC)(), slaveIC);

// Info<< "MANIPULATING GRADIENT COEFFS" << endl;
// Info<< "Slave internal coeffs: " << slaveIC << nl
// << "Slave boundary coeffs: " << slaveBC << nl
// << "Master internal coeffs: " << masterIC << nl
// << "Master boundary coeffs: " << masterBC << endl;
}
}

Expand Down
2 changes: 2 additions & 0 deletions src/finiteVolume/fields/fvPatchFields/constraint/ggi/ggiFvPatchField.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -290,6 +290,8 @@ public:

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#include "ggiFvPatchScalarField.H"

#ifdef NoRepository
# include "ggiFvPatchField.C"
#endif
Expand Down
115 changes: 76 additions & 39 deletions src/finiteVolume/fields/fvPatchFields/constraint/ggi/ggiFvPatchScalarField.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,13 +25,14 @@ Author
Vuko Vukcevic, Wikki Ltd. All rights reserved
Description
Specilalisation of manipulateMatrix member function for scalars needed to
Specilalisation of patchFlux member function for scalars needed to
ensure conservation across GGI patches that are partially covered if the
bridge overlap switch is on.
\*---------------------------------------------------------------------------*/

#include "ggiFvPatchField.H"
#include "surfaceFields.H"
#include "fvScalarMatrix.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Expand All @@ -42,51 +43,87 @@ namespace Foam
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //

template<>
void Foam::ggiFvPatchField<scalar>::manipulateMatrix(fvMatrix<scalar>& matrix)
void Foam::ggiFvPatchField<scalar>::patchFlux
(
GeometricField<scalar, fvsPatchField, surfaceMesh>& flux,
const fvMatrix<scalar>& matrix
) const
{
Info<< "IN MANIPULATE MATRIX!" << endl;
// Since we have adjusted the internal/boundary coefficients in the
// manipulateMatrix member function below, we must not use
// patchNeighbourField for reconstructing the flux. We only need to use
// interpolated shadow field. VV, 6/Mar/2018.

// Get patch ID
const label patchI = this->patch().index();

// Get internal field
const scalarField& iField = this->internalField();

// Get shadow face-cells and assemble shadow field
const unallocLabelList& sfc = ggiPatch_.shadow().faceCells();

scalarField sField(sfc.size());
forAll (sField, i)
{
sField[i] = iField[sfc[i]];
}

// Interpolate shadow to this side. Note: must not bridge since internal
// coeffs and boundary coeffs take it into account
scalarField neighbourField(ggiPatch_.interpolate(sField));

// Conservative treatment for bridged overlap
if (ggiPatch_.bridgeOverlap())
{
// Get this patch and shadow patch index
const label patchI = this->patch().index();
const label sPatchI = ggiPatch_.shadowIndex();
// Only set fully uncovered faces. Partially covered faces taken into
// account by manipulating value and gradient matrix coefficients. Note:
// mirror field is the same as patch internal field
ggiPatch_.setUncoveredFaces
(
this->patchInternalField()(),
neighbourField
);
}

// Calculate the flux with correct neighbour field (fully uncovered faces
// bridged, while partially uncovered faces taken into account by
// manipulating value and gradient matrix coefficients in order to ensure
// conservation for both convection and diffusion part across partially
// overlapping faces). VV, 14/Mar/2018.
flux.boundaryField()[patchI] =
matrix.internalCoeffs()[patchI]*this->patchInternalField()
- matrix.boundaryCoeffs()[patchI]*neighbourField;

// Scale the flux on slave patch to ensure global conservation across this
// partially overlapping GGI patch. The slight disbalance can happen since
// we interpolate the matrix coeffs and field separately, and we should do
// it together to ensure conservation. Current code design does not easily
// allow this, so this is a meaningful temporary solution when we have
// partially overlapping faces. VV, 14/Mar/2018.
if (ggiPatch_.bridgeOverlap() && !ggiPatch_.master())
{
// This is slave patch, master already updated the fluxes and we can use
// that info to scale the fluxes on this side

// Get all matrix coefficients
scalarField& thisIC = matrix.internalCoeffs()[patchI];
scalarField& thisBC = matrix.boundaryCoeffs()[patchI];
scalarField& shadowIC = matrix.internalCoeffs()[sPatchI];
scalarField& shadowBC = matrix.boundaryCoeffs()[sPatchI];
// Get the total flux through master patch
const scalar masterFlux =
gSum(flux.boundaryField()[ggiPatch_.shadowIndex()]);

// if (!ggiPatch_.master())
{
// This is master, interpolate shadow on this side
const scalarField shadowInterpolatedIC = ggiPatch_.interpolate(shadowIC);
const scalarField shadowInterpolatedBC = ggiPatch_.interpolate(shadowBC);

// Set new internal coeffs using the data from the other side and
// taking into account partially and fully uncovered faces
const scalarField origIC = thisIC;
thisIC = shadowInterpolatedBC;
ggiPatch_.scaleForPartialCoverage(origIC, thisIC);
ggiPatch_.scaleForPartialCoverage(origIC, thisIC);

// Set new boundary coeffs using the data from the other side and
// taking into account partially and fully uncovered faces
const scalarField origBC = thisBC;
thisBC = shadowInterpolatedIC;
ggiPatch_.scaleForPartialCoverage(origBC, thisBC);
ggiPatch_.scaleForPartialCoverage(origBC, thisBC);

Info<< "This new internal coeffs: " << thisIC << nl
<< "This old internal coeffs: " << origIC << nl
<< "This shadow boundary coeffs" << shadowBC << nl << endl;

Info<< "This new boundary coeffs: " << thisBC << nl
<< "This old boundary coeffs: " << origBC << nl
<< "This shadow internale coeffs" << shadowIC << nl << endl;
// Get the total flux through slave patch
const scalar slaveFlux = gSum(flux.boundaryField()[patchI]);

// Calculate the scaling factor. Note: negative sign since master and
// slave fluxes have opposite sign
const scalar scalingFactor = -masterFlux/(slaveFlux + SMALL);

// Scale the slave flux to ensure global patch conservation
flux.boundaryField()[patchI] *= scalingFactor;

if (debug)
{
Info<< "Scaling flux on patch: " << ggiPatch_.name()
<< " with " << scalingFactor
<< " to ensure conservation." << endl;
}
}
}
Expand Down
9 changes: 5 additions & 4 deletions src/finiteVolume/fields/fvPatchFields/constraint/ggi/ggiFvPatchScalarField.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@ Author
Vuko Vukcevic, Wikki Ltd. All rights reserved
Description
Specilalisation of manipulateMatrix member function for scalars needed to
Specilalisation of patchFlux member function for scalars needed to
ensure conservation across GGI patches that are partially covered if the
bridge overlap switch is on.
Expand All @@ -47,10 +47,11 @@ namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

template<>
void Foam::ggiFvPatchField<scalar>::manipulateMatrix
void Foam::ggiFvPatchField<scalar>::patchFlux
(
fvMatrix<scalar>& matrix
);
GeometricField<scalar, fvsPatchField, surfaceMesh>& flux,
const fvMatrix<scalar>& matrix
) const;


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Expand Down
23 changes: 16 additions & 7 deletions src/finiteVolume/fvMesh/fvPatches/constraint/ggi/ggiFvPatch.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -156,15 +156,24 @@ void Foam::ggiFvPatch::makeCorrVecs(vectorField& cv) const
// Non-orthogonality correction on a ggi interface
// MB, 7/April/2009

// Calculate correction vectors on coupled patches
const scalarField& patchDeltaCoeffs = deltaCoeffs();
// No non-orthognal correction if the bridge overlap is switched on to
// ensure conservative interpolation for partially overlapping faces
if (bridgeOverlap())
{
cv = vector::zero;
}
else
{
// Calculate correction vectors on coupled patches
const scalarField& patchDeltaCoeffs = deltaCoeffs();

const vectorField patchDeltas = delta();
const vectorField n = nf();
const vectorField patchDeltas = delta();
const vectorField n = nf();

// If non-orthogonality is over 90 deg, kill correction vector
// HJ, 6/Jan/2011
cv = pos(patchDeltas & n)*(n - patchDeltas*patchDeltaCoeffs);
// If non-orthogonality is over 90 deg, kill correction vector
// HJ, 6/Jan/2011
cv = pos(patchDeltas & n)*(n - patchDeltas*patchDeltaCoeffs);
}
}


Expand Down

0 comments on commit 88cdcdc

Please sign in to comment.