Menu
▾
▴
#87 Interpolation Error on FAM Mesh with Cyclic BC
Hi
I remembered the bug tracker after submitting the bug to the cfd-online forum, so please follow the link:
http://www.cfd-online.com/Forums/openfoam-bugs/91245-interpolation-error-fam-mesh-cyclic-bcs.html
I have made a small solver and test case, which I have attached to this bug report. Just run the AllwmakeRun-script.
Kind regards,
Niels
I have been looking for an explanation, and I have found the possible bug in
OpenFOAM/OpenFOAM-1.6-ext/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolationScheme/edgeInterpolationScheme.C
Two types of interpolations are defined, namely "interpolate" and "euclidianInterpolate". The latter uses standard interpolation over the coupled patch, namely:
tsf().boundaryField()[pi] =
pLambda*vf.boundaryField()[pi].patchInternalField()
+ (1 - pLambda)*vf.boundaryField()[pi].patchNeighbourField();
where pLambda are the interpolation weights. I have confirmed that they are correct (see earlier comment in the bug report).
The "interpolate" routine, however, applies transformation tensors on the areaField<Type>, which is no problem when it is performs on a scalar, however, it gives wrong results on vectors and tensors. Interpolation routine:
for (label i=0; i<size; i++)
{
const tensorField& curT =
mesh.edgeTransformTensors()[start + i];
const tensor& Te = curT[0];
const tensor& TP = curT[1];
const tensor& TN = curT[2];
pSf[i] =
transform
(
Te.T(),
pLambda[i]*transform(TP, pOwnVf[i])
+ (1 - pLambda[i])*transform(TN, pNgbVf[i])
);
}
Furthermore the first interpolation routine has been used previously, but has been commented and replaced by the transformation routine.
Hope this has helped shedding a bit of a light on the problem.
/ Niels
I have further dug into the code, and I have narrowed it down to the computation of the transformation tensors on the boundary, thus interpolation by the use of transformation vectors are actually OK!
The problems arise in the method "void faMesh::calcEdgeTransformTensors() const" in faMeshDemandDrivenData.C on line 857. The current implementation is:
// ----------------------
il = ngbCf[edgeI] - E;
// ----------------------
which is correct for coupled patches that are NOT cyclic, e.g. parallel patches, as the face centre of the neighbouring face is on the correct side of the side. For cyclic patches, on the other hand, both the neighbour and owner faces are on the same side of the boundary edge, hence the code in that case should read
// ----------------------
il = E - ngbCf[edgeI];
// ----------------------
This means that the line 857 should be replaced by something like
// ----------------------
if ( patch == cyclicTypePatch )
il = E - ngbCf[edgeI];
else
il = ngbCf[edgeI] - E;
// ----------------------
Is there an easy way to make the check " patch == cyclicTypePatch " ?
/ Niels
A suggestion to fix the problem, however, it would require someone with a broad knowledge of all boundary conditions, to see, if it could potentially be problematic. As stated above, the sign of "il" on line 857 depends on the type of coupled boundary condition. My suggestion is to change line 857 from
// ----------------------
il = ngbCf[edgeI] - E;
// ----------------------
to
// ----------------------
il = Foam::sign( ( ngbCf[edgeI] - E ) & ( E - Cf.internalField()[edgeFaces[edgeI]] ) ) * (ngbCf[edgeI] - E );
// ----------------------
which should automatically take care of the sign change; at least for processor and cyclic boundaries. Problems with other types of coupled BCs, which could arise is currently beyond my insight.
Kind regards,
Niels