Commit 691e7186 authored by Henry Weller's avatar Henry Weller
Browse files

ACMI: Corrected conservation issue

Patch contributed by Mattijs Janssens
Resolves bug-report http://bugs.openfoam.org/view.php?id=2057
parent e34763b3
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013-2015 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012-2015 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2012-2016 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012-2015 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2012-2016 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......
......@@ -997,7 +997,7 @@ bool Foam::polyMesh::upToDatePoints(const regIOobject& io) const
void Foam::polyMesh::setUpToDatePoints(regIOobject& io) const
{
io.eventNo() = points_.eventNo();
io.eventNo() = points_.eventNo()+1;
}
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......
......@@ -137,22 +137,15 @@ Foam::tmp<Foam::Field<Type>>
Foam::cyclicACMIFvPatchField<Type>::patchNeighbourField() const
{
const Field<Type>& iField = this->primitiveField();
const labelUList& nbrFaceCellsCoupled =
cyclicACMIPatch_.cyclicACMIPatch().neighbPatch().faceCells();
const labelUList& faceCellsNonOverlap =
cyclicACMIPatch_.cyclicACMIPatch().nonOverlapPatch().faceCells();
Field<Type> pnfCoupled(iField, nbrFaceCellsCoupled);
Field<Type> pfNonOverlap(iField, faceCellsNonOverlap);
const cyclicACMIPolyPatch& cpp = cyclicACMIPatch_.cyclicACMIPatch();
tmp<Field<Type>> tpnf
(
new Field<Type>
cyclicACMIPatch_.interpolate
(
cyclicACMIPatch_.interpolate
Field<Type>
(
pnfCoupled,
pfNonOverlap
iField,
cpp.neighbPatch().faceCells()
)
)
);
......@@ -207,10 +200,12 @@ void Foam::cyclicACMIFvPatchField<Type>::updateInterfaceMatrix
const Pstream::commsTypes
) const
{
const cyclicACMIPolyPatch& cpp = cyclicACMIPatch_.cyclicACMIPatch();
// note: only applying coupled contribution
const labelUList& nbrFaceCellsCoupled =
cyclicACMIPatch_.cyclicACMIPatch().neighbPatch().faceCells();
cpp.neighbPatch().faceCells();
scalarField pnf(psiInternal, nbrFaceCellsCoupled);
......@@ -237,10 +232,11 @@ void Foam::cyclicACMIFvPatchField<Type>::updateInterfaceMatrix
const Pstream::commsTypes
) const
{
const cyclicACMIPolyPatch& cpp = cyclicACMIPatch_.cyclicACMIPatch();
// note: only applying coupled contribution
const labelUList& nbrFaceCellsCoupled =
cyclicACMIPatch_.cyclicACMIPatch().neighbPatch().faceCells();
const labelUList& nbrFaceCellsCoupled = cpp.neighbPatch().faceCells();
Field<Type> pnf(psiInternal, nbrFaceCellsCoupled);
......@@ -258,149 +254,6 @@ void Foam::cyclicACMIFvPatchField<Type>::updateInterfaceMatrix
}
template<class Type>
Foam::tmp<Foam::Field<Type>> Foam::cyclicACMIFvPatchField<Type>::snGrad
(
const scalarField& deltaCoeffs
) const
{
// note: only applying coupled contribution
return coupledFvPatchField<Type>::snGrad(deltaCoeffs);
}
template<class Type>
void Foam::cyclicACMIFvPatchField<Type>::updateCoeffs()
{
// update non-overlap patch - some will implement updateCoeffs, and
// others will implement evaluate
// scale neighbour field by (1 - mask)
const scalarField& mask = cyclicACMIPatch_.cyclicACMIPatch().mask();
const fvPatchField<Type>& npf = nonOverlapPatchField();
const_cast<fvPatchField<Type>&>(npf).updateCoeffs(1.0 - mask);
}
template<class Type>
void Foam::cyclicACMIFvPatchField<Type>::initEvaluate
(
const Pstream::commsTypes comms
)
{
// update non-overlap patch (if not already updated by updateCoeffs)
// scale neighbour field by (1 - mask)
fvPatchField<Type>& npf =
const_cast<fvPatchField<Type>&>(nonOverlapPatchField());
if (!npf.updated())
{
const scalarField& mask = cyclicACMIPatch_.cyclicACMIPatch().mask();
npf.evaluate(comms);
npf *= 1.0 - mask;
}
}
template<class Type>
void Foam::cyclicACMIFvPatchField<Type>::evaluate
(
const Pstream::commsTypes comms
)
{
// blend contributions from the coupled and non-overlap patches
// neighbour patch field is updated via updateCoeffs or initEvaluate
// and is already scaled by (1 - mask)
const fvPatchField<Type>& npf = nonOverlapPatchField();
coupledFvPatchField<Type>::evaluate(comms);
const Field<Type>& cpf = *this;
const scalarField& mask = cyclicACMIPatch_.cyclicACMIPatch().mask();
Field<Type>::operator=(mask*cpf + npf);
fvPatchField<Type>::evaluate();
}
template<class Type>
Foam::tmp<Foam::Field<Type>>
Foam::cyclicACMIFvPatchField<Type>::valueInternalCoeffs
(
const tmp<scalarField>& w
) const
{
// note: do not blend based on mask field
// - when applied this is scaled by the areas which are already scaled
return coupledFvPatchField<Type>::valueInternalCoeffs(w);
}
template<class Type>
Foam::tmp<Foam::Field<Type>>
Foam::cyclicACMIFvPatchField<Type>::valueBoundaryCoeffs
(
const tmp<scalarField>& w
) const
{
// note: do not blend based on mask field
// - when applied this is scaled by the areas which are already scaled
return coupledFvPatchField<Type>::valueBoundaryCoeffs(w);
}
template<class Type>
Foam::tmp<Foam::Field<Type>>
Foam::cyclicACMIFvPatchField<Type>::gradientInternalCoeffs
(
const scalarField& deltaCoeffs
) const
{
// note: do not blend based on mask field
// - when applied this is scaled by the areas which are already scaled
return coupledFvPatchField<Type>::gradientInternalCoeffs(deltaCoeffs);
}
template<class Type>
Foam::tmp<Foam::Field<Type>>
Foam::cyclicACMIFvPatchField<Type>::gradientInternalCoeffs() const
{
// note: do not blend based on mask field
// - when applied this is scaled by the areas which are already scaled
return coupledFvPatchField<Type>::gradientInternalCoeffs();
}
template<class Type>
Foam::tmp<Foam::Field<Type>>
Foam::cyclicACMIFvPatchField<Type>::gradientBoundaryCoeffs
(
const scalarField& deltaCoeffs
) const
{
// note: do not blend based on mask field
// - when applied this is scaled by the areas which are already scaled
return coupledFvPatchField<Type>::gradientBoundaryCoeffs(deltaCoeffs);
}
template<class Type>
Foam::tmp<Foam::Field<Type>>
Foam::cyclicACMIFvPatchField<Type>::gradientBoundaryCoeffs() const
{
// note: do not blend based on mask field
// - when applied this is scaled by the areas which are already scaled
return coupledFvPatchField<Type>::gradientBoundaryCoeffs();
}
template<class Type>
void Foam::cyclicACMIFvPatchField<Type>::manipulateMatrix
(
......
......@@ -218,63 +218,6 @@ public:
const Pstream::commsTypes commsType
) const;
//- Return patch-normal gradient
virtual tmp<Field<Type>> snGrad
(
const scalarField& deltaCoeffs
) const;
//- Update the coefficients associated with the patch field
void updateCoeffs();
//- Initialise the evaluation of the patch field
virtual void initEvaluate
(
const Pstream::commsTypes commsType
);
//- Evaluate the patch field
virtual void evaluate
(
const Pstream::commsTypes commsType
);
//- Return the matrix diagonal coefficients corresponding to the
// evaluation of the value of this patchField with given weights
virtual tmp<Field<Type>> valueInternalCoeffs
(
const tmp<scalarField>&
) const;
//- Return the matrix source coefficients corresponding to the
// evaluation of the value of this patchField with given weights
virtual tmp<Field<Type>> valueBoundaryCoeffs
(
const tmp<scalarField>&
) const;
//- Return the matrix diagonal coefficients corresponding to the
// evaluation of the gradient of this patchField
virtual tmp<Field<Type>> gradientInternalCoeffs
(
const scalarField& deltaCoeffs
) const;
//- Return the matrix diagonal coefficients corresponding to the
// evaluation of the gradient of this patchField
virtual tmp<Field<Type>> gradientInternalCoeffs() const;
//- Return the matrix source coefficients corresponding to the
// evaluation of the gradient of this patchField
virtual tmp<Field<Type>> gradientBoundaryCoeffs
(
const scalarField& deltaCoeffs
) const;
//- Return the matrix source coefficients corresponding to the
// evaluation of the gradient of this patchField
virtual tmp<Field<Type>> gradientBoundaryCoeffs() const;
//- Manipulate matrix
virtual void manipulateMatrix(fvMatrix<Type>& matrix);
......
......@@ -43,7 +43,12 @@ void Foam::cyclicACMIFvPatch::updateAreas() const
{
if (cyclicACMIPolyPatch_.updated())
{
// Set Sf and magSf for both sides' coupled and non-overlapping patches
if (debug)
{
Pout<< "cyclicACMIFvPatch::updateAreas() : updating fv areas for "
<< name() << " and " << this->nonOverlapPatch().name()
<< endl;
}
// owner couple
const_cast<vectorField&>(Sf()) = patch().faceAreas();
......@@ -81,25 +86,37 @@ void Foam::cyclicACMIFvPatch::makeWeights(scalarField& w) const
if (coupled())
{
const cyclicACMIFvPatch& nbrPatch = neighbFvPatch();
const fvPatch& nbrPatchNonOverlap = nonOverlapPatch();
const scalarField deltas(nf() & coupledFvPatch::delta());
const scalarField nbrDeltas
// These deltas are of the cyclic part alone - they are
// not affected by the amount of overlap with the nonOverlapPatch
scalarField nbrDeltas
(
interpolate
(
nbrPatch.nf() & nbrPatch.coupledFvPatch::delta(),
nbrPatchNonOverlap.nf() & nbrPatchNonOverlap.delta()
nbrPatch.nf() & nbrPatch.coupledFvPatch::delta()
)
);
scalar tol = cyclicACMIPolyPatch::tolerance();
forAll(deltas, facei)
{
scalar di = deltas[facei];
scalar dni = nbrDeltas[facei];
w[facei] = dni/(di + dni);
if (dni < tol)
{
// Avoid zero weights on disconnected faces. This value
// will be weighted with the (zero) face area so will not
// influence calculations.
w[facei] = 1.0;
}
else
{
w[facei] = dni/(di + dni);
}
}
}
else
......@@ -122,30 +139,12 @@ Foam::tmp<Foam::vectorField> Foam::cyclicACMIFvPatch::delta() const
{
if (coupled())
{
const cyclicACMIFvPatch& nbrPatchCoupled = neighbFvPatch();
const fvPatch& nbrPatchNonOverlap = nonOverlapPatch();
const cyclicACMIFvPatch& nbrPatch = neighbFvPatch();
const vectorField patchD(coupledFvPatch::delta());
vectorField nbrPatchD
(
interpolate
(
nbrPatchCoupled.coupledFvPatch::delta(),
nbrPatchNonOverlap.delta()
)
);
const vectorField nbrPatchD0
(
interpolate
(
vectorField(nbrPatchCoupled.size(), Zero),
nbrPatchNonOverlap.delta()()
)
);
vectorField nbrPatchD(interpolate(nbrPatch.coupledFvPatch::delta()));
nbrPatchD -= nbrPatchD0;
tmp<vectorField> tpdv(new vectorField(patchD.size()));
vectorField& pdv = tpdv.ref();
......
......@@ -174,10 +174,11 @@ public:
//- Return delta (P to N) vectors across coupled patch
virtual tmp<vectorField> delta() const;
//- Interpolate (make sure to have uptodate areas)
template<class Type>
tmp<Field<Type>> interpolate
(
const Field<Type>& fldCoupled
const Field<Type>& fld
) const
{
updateAreas();
......@@ -185,57 +186,18 @@ public:
return
cyclicACMIPolyPatch_.cyclicAMIPolyPatch::interpolate
(
fldCoupled
fld
);
}
//- Interpolate (make sure to have uptodate areas)
template<class Type>
tmp<Field<Type>> interpolate
(
const tmp<Field<Type>>& tfldCoupled
const tmp<Field<Type>>& tfld
) const
{
updateAreas();
return
cyclicACMIPolyPatch_.cyclicAMIPolyPatch::interpolate
(
tfldCoupled
);
}
template<class Type>
tmp<Field<Type>> interpolate
(
const Field<Type>& fldCoupled,
const Field<Type>& fldNonOverlap
) const
{
updateAreas();
return
cyclicACMIPolyPatch_.interpolate
(
fldCoupled,
fldNonOverlap
);
}
template<class Type>
tmp<Field<Type>> interpolate
(
const tmp<Field<Type>>& tFldCoupled,
const tmp<Field<Type>>& tFldNonOverlap
) const
{
updateAreas();
return
cyclicACMIPolyPatch_.interpolate
(
tFldCoupled,
tFldNonOverlap
);
return interpolate(tfld());
}
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......
......@@ -378,10 +378,17 @@ public:
//- Return const access to source patch weights
inline const scalarListList& srcWeights() const;
//- Return access to source patch weights
inline scalarListList& srcWeights();
//- Return const access to normalisation factor of source
// patch weights (i.e. the sum before normalisation)
inline const scalarField& srcWeightsSum() const;
//- Return access to normalisation factor of source
// patch weights (i.e. the sum before normalisation)
inline scalarField& srcWeightsSum();
//- Source map pointer - valid only if singlePatchProc = -1
// This gets source data into a form to be consumed by
// tgtAddress, tgtWeights
......@@ -399,10 +406,17 @@ public:
//- Return const access to target patch weights
inline const scalarListList& tgtWeights() const;
//- Return access to target patch weights
inline scalarListList& tgtWeights();
//- Return const access to normalisation factor of target
// patch weights (i.e. the sum before normalisation)
inline const scalarField& tgtWeightsSum() const;
//- Return access to normalisation factor of target
// patch weights (i.e. the sum before normalisation)
inline scalarField& tgtWeightsSum();
//- Target map pointer - valid only if singlePatchProc=-1.
// This gets target data into a form to be consumed by
// srcAddress, srcWeights
......
......@@ -72,6 +72,14 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::srcWeights() const
}
template<class SourcePatch, class TargetPatch>
inline Foam::scalarListList&
Foam::AMIInterpolation<SourcePatch, TargetPatch>::srcWeights()
{
return srcWeights_;
}
template<class SourcePatch, class TargetPatch>
inline const Foam::scalarField&
Foam::AMIInterpolation<SourcePatch, TargetPatch>::srcWeightsSum() const
......@@ -80,6 +88,14 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::srcWeightsSum() const
}
template<class SourcePatch, class TargetPatch>
inline Foam::scalarField&
Foam::AMIInterpolation<SourcePatch, TargetPatch>::srcWeightsSum()
{
return srcWeightsSum_;
}
template<class SourcePatch, class TargetPatch>
inline const Foam::mapDistribute&
Foam::AMIInterpolation<SourcePatch, TargetPatch>::srcMap() const
......@@ -112,6 +128,14 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::tgtWeights() const
}
template<class SourcePatch, class TargetPatch>
inline Foam::scalarListList&
Foam::AMIInterpolation<SourcePatch, TargetPatch>::tgtWeights()
{
return tgtWeights_;
}