CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Bugs (http://www.cfd-online.com/Forums/openfoam-bugs/)
-   -   Possible bug in meshPhiC please take a look (http://www.cfd-online.com/Forums/openfoam-bugs/62521-possible-bug-meshphic-please-take-look.html)

jaswi December 20, 2007 09:30

Dear Developers Wishing yo
 
Dear Developers

Wishing you all a very nice day.

I request you to take a look at this post and recheck whether the meshPhi.C has an error. I am using interFoam with mesh motion and thus need to use meshPhi(rho,U) for making fluxes absolute and relative when required. I followed the trail and below is what I found:


The file /fvc/fvcMeshPhi.H has an overloaded function meshPhi() alongwith other functions.
-----------
meshPhi(U):
-----------

tmp<surfacescalarfield> meshPhi
(
const volVectorField& U
);

---------------
meshPhi(rho,U): when rho - dimensionedScalar
---------------


tmp<surfacescalarfield> meshPhi
(
const dimensionedScalar& rho,
const volVectorField& U
);

---------------
meshPhi(rho,U): when rho - volScalarField
---------------

tmp<surfacescalarfield> meshPhi
(
const volScalarField& rho,
const volVectorField& U
);

Now the corresponding implementation of these functions in /fvc/meshPhi.C are


-----------
meshPhi(U):
-----------
tmp<surfacescalarfield> meshPhi
(
const volVectorField& vf
)
{
Info << "Calling 1st constructor" << endl;

return fv::ddtScheme<vector>::New
(
vf.mesh(),
vf.mesh().ddtScheme("ddt(" + vf.name() + ')')
)().meshPhi(vf);
}


---------------
meshPhi(rho,U): when rho - dimensionedScalar
---------------
tmp<surfacescalarfield> meshPhi
(
const dimensionedScalar& rho,
const volVectorField& vf
)
{
Info << "Calling 2nd constructor" << endl;
return fv::ddtScheme<vector>::New
(
vf.mesh(),
vf.mesh().ddtScheme("ddt(" + rho.name() + ',' + vf.name() + ')')
)().meshPhi(vf);
}




---------------
meshPhi(rho,U): when rho - volScalarField
---------------

tmp<surfacescalarfield> meshPhi
(
const volScalarField& rho,
const volVectorField& vf
)
{
Info << "Calling 3rd constructor" << endl;
return fv::ddtScheme<vector>::New
(
vf.mesh(),
vf.mesh().ddtScheme("ddt(" + rho.name() + ',' + vf.name() + ')')
)().meshPhi(vf);
}


The point is that there is no difference whether I call meshPhi(U) or meshPhi(rho,U).

This issue has been raised before. please take a look at this post:

http://www.cfd-online.com/OpenFOAM_D...es/1/3548.html

In this post by Fredrick, he at least gets an error :
-------------------------------------------------
phi += fvc::meshPhi(U).
For a compressible flow, the fluxes should be:
phi += fvc::meshPhi(rho,U)

This does not work, I got the following error:
--> FOAM FATAL ERROR : Different dimensions for +=
dimensions : [1 0 -1 0 0 0 0] = [0 3 -1 0 0 0 0]
-------------------------------------------------

and Prof. Hrvoje Jasak and Eric Lillberg had a discussion over which is the best method to correct it.


I took the clue from the post and multiplied the return by rho as suggested in that post:
*********************************************
template<class>
tmp<surfacescalarfield> meshPhi
(
const dimensionedScalar& rho,
const GeometricField<type,>& vf
)
{
//HJ, Missing rho. 13/Dec/2006
return rho*fv::ddtScheme<type>::New
(
vf.mesh(),
vf.mesh().ddtScheme("ddt(" + rho.name() + ',' + vf.name() + ')')
)().meshPhi(vf);

}
**********************************************

but attempt to recompile the library gives the following error

openfoam@taifun:~/OpenFOAM/OpenFOAM-1.4.1/src/finiteVolume> wmake libso
SOURCE=finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtSchemes.C ; g++ -m64 -Dlinux64 -DDP -Wall -Wno-strict-aliasing -Wextra -Wno-unused-parameter -Wold-style-cast -march=opteron -O3 -DNoRepository -ftemplate-depth-40 -I/home/openfoam/OpenFOAM/OpenFOAM-1.4.1/src/triSurface/lnInclude -I/home/openfoam/OpenFOAM/OpenFOAM-1.4.1/src/meshTools/lnInclude -IlnInclude -I. -I/home/openfoam/OpenFOAM/OpenFOAM-1.4.1/src/OpenFOAM/lnInclude -fPIC -c $SOURCE -o Make/linux64GccDPOpt/EulerDdtSchemes.o
SOURCE=finiteVolume/fvc/fvcMeshPhi.C ; g++ -m64 -Dlinux64 -DDP -Wall -Wno-strict-aliasing -Wextra -Wno-unused-parameter -Wold-style-cast -march=opteron -O3 -DNoRepository -ftemplate-depth-40 -I/home/openfoam/OpenFOAM/OpenFOAM-1.4.1/src/triSurface/lnInclude -I/home/openfoam/OpenFOAM/OpenFOAM-1.4.1/src/meshTools/lnInclude -IlnInclude -I. -I/home/openfoam/OpenFOAM/OpenFOAM-1.4.1/src/OpenFOAM/lnInclude -fPIC -c $SOURCE -o Make/linux64GccDPOpt/fvcMeshPhi.o
finiteVolume/fvc/fvcMeshPhi.C: In function 'Foam::tmp<foam::geometricfield<double,> > Foam::fvc::meshPhi(const Foam::volScalarField&, const Foam::volVectorField&)':
finiteVolume/fvc/fvcMeshPhi.C:84: error: conversion from 'Foam::tmp<foam::field<double> >' to non-scalar type 'Foam::tmp<foam::geometricfield<double,> >' requested
make: *** [Make/linux64GccDPOpt/fvcMeshPhi.o] Error 1


Please take a look and comment whether the bug is still there or it has been corrected and I am making a wrong assumption.

Thanks alot for your attention.
with best regards

jaswinder

henry December 20, 2007 09:49

Yes it has an error. I have c
 
Yes it has an error. I have corrected it in our development version and here is the fix for 1.4.1: http://www.cfd-online.com/OpenFOAM_D...hment_icon.gif fvcMeshPhi.C which I haven't tested but it should be OK, let me know it is works.

Henry

hjasak December 20, 2007 10:12

In all honesty, this is not qu
 
In all honesty, this is not quite as clear-cut as it seems. Firstly, I need the rho as the name of the argument in order to capture the correct ddt scheme. However, what do you expect meshPhi() to return: as simple volume swept by the face in motion consistent with the ddt term, allowing you to manipulate it afterwards.

The conclusion from the previous round was that meshPhi() should return the mesh flux - if you don't like it, just multiply by phi interpolated in the way you want.

So much from my side,

Hrv

jaswi December 20, 2007 12:03

Thanks a lot "Professors" for
 
Thanks a lot "Professors" for a quick response. :-). I made the changes and compiled the library. it now gives this error.


--> FOAM FATAL ERROR : Different dimensions for -=
dimensions : [0 3 -1 0 0 0 0] = [1 0 -1 0 0 0 0]


From function dimensionSet::operator-=(const dimensionSet& ds) const
in file dimensionSet/dimensionSet.C at line 183.

when i do the following

return fvc::interpolate(rho)*fv::ddtScheme<vector>::New
(
vf.mesh(),
vf.mesh().ddtScheme("ddt(" + rho.name() + ',' + vf.name() + ')')
)().meshPhi(vf);

then flux[m^3/s] is equated to [(kg/m^3)*(m^3/sec)]
and this dimensionally not correct. The interesting point is that the post i referred to mentions this same error before the correction was made.

One more idea is that as I am using mesh motion with interfoam . In interfoam's createField we declare and define a surface scalar field "rhoPhi" which has the same dimensionality what the modified meshPhi(rho,U) returns. But when we correct the fluxes we use statements like
phi+=meshPhi(rho,U) and phi-=meshPhi(rho,U). May be there is something wrong with that approach. I have very less experience with interfoam solver thus cannot say it with confidence.

I have been loosing my sleep over it . the loss of liquid phase volume fraction due to the mesh motion has really boggled me.

Please comment what shall be done to correct this problem.

With Best Regards
Jaswinder

henry December 20, 2007 12:16

Did you recompile interFoam?
 
Did you recompile interFoam?

jaswi December 20, 2007 12:21

Thanks again for the discussio
 
Thanks again for the discussion.


Yes I did that. To make sure that it does I

wclean it and the wmake it. But the error is still there.

henry December 20, 2007 12:25

Scrap that, there is no need t
 
Scrap that, there is no need to recompile interFoam.

I removed mesh motion from the standard interFoam in the 1.4.1 release because it wasn't entirely consistent and had various modes of failure. I have written a completely new interDyMFoam for the 1.5 release which includes all kinds of mesh change both geometric and topological, an early version of which created the animation on the front page of our web-site. We hope to release 1.5 in the 1st quarter 2008.

Henry

jaswi December 20, 2007 12:26

Good Evening Henry I looked
 
Good Evening Henry

I looked through the interFoam again. Maybe we should correct both phi and rhoPhi for the moving mesh.

The correction block can have something like:


//Make fluxes absolute
phi+=meshphi(U)
rhoPhi+=meshPhi(rho,U)


//Make fluxes relative
phi-=meshphi(U)
rhoPhi-=meshPhi(rho,U)

where and when required in the mesh motion procedure. It is completely a arrow in the dark and I might ne completely wrong.

Please Comment

With Best Regards
Jaswinder

jaswi December 20, 2007 12:30

Hi Henry Can I request for
 
Hi Henry

Can I request for a Christmas Gift :-)

Please If i can have that solver now because I have been loosing sleep since a week over this issue and you know how the advisors some time react when there are no results and one can't explain what is wrong.

my email is singhj@in.tum.de

Please let me know if I can have it.

With Best Regards
Jaswinder

henry December 20, 2007 12:43

The new solver relies on a man
 
The new solver relies on a many new developments and will not compile with 1.4.1, sorry. We will get 1.5 out as soon as we can.

Henry

jaswi December 20, 2007 13:03

Thanks alot for considering it
 
Thanks alot for considering it.

It seems I have to explore other options then, as I am running out of time.

By the way, is there some mass transfer modeling coming up as well. I just can't wait to have the new version. The more i understand this stuff the more I am surprised what all it can do.

Thank you for the great work.

Wish you a Merry Christmas and a Happy New Year.
Happy Holidays
Jaswinder


All times are GMT -4. The time now is 07:07.