[cosimobianchini]

Hi all,
I want to report a bug found in version 1.4.1-patched regarding cyclic patches.
Is there anyone else who found the same misbehave of cyclic patches in the new version or it is just me doing something wrong?

http://brun.de.unifi.it/docpub/cyclic.tgz

This is the link to download a case in which version 1.3 and 1.4.1 are effectively very different.
It is a quarter of a circular pipe with a swirled inlet velocity. It is an axialsymmetric test.
The old version is performing good while the 1.4.1 patched got something wrong.
In the cyclic.tgz file you found both the results obtained with simpleFoam and the two simulations and the log files.
I could not understand which step of the process is messed up.
Please give it a look,

Cosimo

[henry]

Thank you for the test case and example results. I have fixed the bug and tested it on your case which now behaves very similarly to the 1.3 version (actually converges slightly faster due to the improved under-relaxation practice in 1.4.1). To apply the fix replace OpenFOAM-1.4.1/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C with fvMatrixSolve.C and recompile all the applications which use the solvers. To test the fix you only need to recompile simpleFoam which will include the new solver directly as it is a template.

[msrinath80]

Hi Henry,

Does the current OpenFOAM (1.4.1) in sourceforge have all the bug fixes done so far?

[henry]

The 1.4.1 source pack we release via sourceforge is still the original 1.4.1 release. We will consider releasing another patched version as 1.4.2 before the next full release version 1.5.

[luca]

Hi Henry,

thank you very much for your quick reply. We will do further tests in the next days based on the new fvMatrixSolve patch you posted.

keep in touch

Luca

[gtg627e]

Dear Henry,

If I replace fvMatrixSolve.C with the new one you provide, do I need to recompile finiteVolume, as well as the applications that call vMatrixSolve.C?
Specifically,
OpenFOAM/OpenFOAM-1.4.1/src/finiteVolume/wclean && wmake ?

OpenFOAM/OpenFOAM-1.4.1/applications/solvers/incompressible/simpleFoam/wclean && wmake?

Thank you,

Alessandro

[henry]

Surprisingly you don't have to recompile finiteVolume because fvMatrixSolve.C contains the templated solvers that are not actually used in the finiteVolume library but the finiteVolume applications. Also by replacing fvMatrixSolve.C with a newer file means that the dependency system will notice so you should not need to do wclean before wmake.

[rmagnan]

Hello All,

Currently, cyclic patches are restricted to planar faces. This limitation seems related to the way the transform matrix from one side of the cyclic to its sibling is computed. Currently, it is computed by taking the first pair of matching faces and computing a rotation matrix that will map their normals into one another. Implicitly, this computation assumes the rotation axis is perpendicular to both normals. When this is not the case (most curved patches fall into this category), we get solutions where scalar quantities are correctly periodic, magnitude of vector quantities are also OK, but individual components of the vectors are not periodic.

I introduced a few changes to compute the rotation tensor by also taking into consideration the face centers of the first pair of cyclic faces. I now have a working cyclic interface on non-planar patches. This fix was tested in 1.4.1-dev. In the OpenCFD release of 1.4.1 (with cyclic patch - see cyclic bug), the same changes still exhibit some (yet) unexplained non-periodicity

The following figure shows before/after comparison of the z component of velocity. Attached is a test case I have been using.



test_cyclic2.tgz

in src/OpenFOAM/primitives/transform/transform.H, the following overloaded function was added:


inline tensor rotationTensor
(
const vector& p1, // center of face 1
const vector& p2, // center of face 2
const vector& n1, // normal to face 1
const vector& n2 // normal to face 2
)
{
// compute two orthonormal basis
// for each coord. system:
// p is an arbitrary vector defined from the origin
// to the base point where the normal n is attached
// z is aligned with the normal n at point p
// x is computed by cross product of z with p
// y is orthogonal to both x and z to make a right hand system
vector z1 = n1/mag(n1);
vector x1 = z1^p1;
scalar magx1 = mag(x1);
if ( magx1 <= VSMALL )
{
// points and vectors are in the same plane
// the "default" computation is valid
return rotationTensor( n1, n2 );
}
x1 /= magx1;
vector y1 = z1^x1;
tensor E1( x1.x(), y1.x(), z1.x(),
x1.y(), y1.y(), z1.y(),
x1.z(), y1.z(), z1.z() );

vector z2 = n2/mag(n2);
vector x2 = z2^p2; x2 /= mag(x2);
vector y2 = z2^x2;
tensor E2( x2.x(), y2.x(), z2.x(),
x2.y(), y2.y(), z2.y(),
x2.z(), y2.z(), z2.z() );

// the rotation tensor defining a rotation from E1 to E2
// is computed by writing base2 with respect to base1
return E2 & E1.T();
}


The code in coupledPolyPatch::calcTransformTensors and cyclicPolyPatch::order has been modified to use the new interface - replacing the calls to the original rotationTensor by calls to this new one adding appropriate face centers as the first two arguments.

[sek]

I have a case for which I have non-planar cyclic boundaries. And I am getting is a very unreasonable solution. It appears that the fix your suggested has NOT made into 1.5 or 1.5-dev. What was the verdict?

[sek]

Originally Posted by sek

I have a case for which I have non-planar cyclic boundaries. And I am getting is a very unreasonable solution. It appears that the fix your suggested has NOT made into 1.5 or 1.5-dev. What was the verdict?

I have a case for which I have non-planar cyclic boundaries. And I am getting is a very unreasonable solution. It appears that the fix your suggested has NOT made into 1.5 or 1.5-dev. What was the verdict?

[stevin10]

Originally Posted by rmagnan

Hello All,

Currently, cyclic patches are restricted to planar faces. This limitation seems related to the way the transform matrix from one side of the cyclic to its sibling is computed. Currently, it is computed by taking the first pair of matching faces and computing a rotation matrix that will map their normals into one another. Implicitly, this computation assumes the rotation axis is perpendicular to both normals. When this is not the case (most curved patches fall into this category), we get solutions where scalar quantities are correctly periodic, magnitude of vector quantities are also OK, but individual components of the vectors are not periodic.

I introduced a few changes to compute the rotation tensor by also taking into consideration the face centers of the first pair of cyclic faces. I now have a working cyclic interface on non-planar patches. This fix was tested in 1.4.1-dev. In the OpenCFD release of 1.4.1 (with cyclic patch - see cyclic bug), the same changes still exhibit some (yet) unexplained non-periodicity

The following figure shows before/after comparison of the z component of velocity. Attached is a test case I have been using.



test_cyclic2.tgz

in src/OpenFOAM/primitives/transform/transform.H, the following overloaded function was added:


inline tensor rotationTensor
(
const vector& p1, // center of face 1
const vector& p2, // center of face 2
const vector& n1, // normal to face 1
const vector& n2 // normal to face 2
)
{
// compute two orthonormal basis
// for each coord. system:
// p is an arbitrary vector defined from the origin
// to the base point where the normal n is attached
// z is aligned with the normal n at point p
// x is computed by cross product of z with p
// y is orthogonal to both x and z to make a right hand system
vector z1 = n1/mag(n1);
vector x1 = z1^p1;
scalar magx1 = mag(x1);
if ( magx1 <= VSMALL )
{
// points and vectors are in the same plane
// the "default" computation is valid
return rotationTensor( n1, n2 );
}
x1 /= magx1;
vector y1 = z1^x1;
tensor E1( x1.x(), y1.x(), z1.x(),
x1.y(), y1.y(), z1.y(),
x1.z(), y1.z(), z1.z() );

vector z2 = n2/mag(n2);
vector x2 = z2^p2; x2 /= mag(x2);
vector y2 = z2^x2;
tensor E2( x2.x(), y2.x(), z2.x(),
x2.y(), y2.y(), z2.y(),
x2.z(), y2.z(), z2.z() );

// the rotation tensor defining a rotation from E1 to E2
// is computed by writing base2 with respect to base1
return E2 & E1.T();
}


The code in coupledPolyPatch::calcTransformTensors and cyclicPolyPatch:rder has been modified to use the new interface - replacing the calls to the original rotationTensor by calls to this new one adding appropriate face centers as the first two arguments.

Hi,

I happened to come across this post of yours where you describe a solution to the vector component wise cyclic boundary issue. It seems that this still has not been addressed in the newest versions of OpenFOAM-1.6 & 1.7.1. I am currently trying to model laminar, multiphase periodic pipe flows and constantly have problems with the radial velocity components at the boundaries.

I know this post was a long time ago, but I was wondering if you still had a copy of these codes that I could simply implement. I have no idea how to modify the coupledPolyPatch & cyclicPolyPatch codes, that you describe at the end of your post.

I am hoping you can still be of help. Will be very much appreciated!

Regards Stevin