How to solve alphaEqn using PythonFlu?
 

Hi foamers,

I am writing twoliquidMixingFoam solver in python using PythonFlu. I have done everything except alphaEqn.H

{
fvScalarMatrix alpha1Eqn
(
fvm::ddt(alpha1)
+ fvm::div(phi, alpha1)
- fvm::laplacian
(
Dab + alphatab*turbulence->nut(), alpha1,
"laplacian(Dab,alpha1)"
)
);

alpha1Eqn.solve();

rhoPhi = alpha1Eqn.flux()*(rho1 - rho2) + phi*rho2;
rho = alpha1*rho1 + (scalar(1) - alpha1)*rho2;

Info<< "Phase 1 volume fraction = "
<< alpha1.weightedAverage(mesh.V()).value()
<< " Min(alpha1) = " << min(alpha1).value()
<< " Max(alpha1) = " << max(alpha1).value()
<< endl;
}

How can I solve this in PythonFlu? In particularly I have doubt in turbulence->nut().


I have tried this. See below code

def alphaEqn( mesh, phi, alpha1, rhoPhi, rho1, rho2, interface):
from Foam import fvm
alpha1Eqn = fvm.ddt(alpha1) + fvm.div(phi, alpha1) - fvm.laplacian( Dab + alphatab*turbulence.nut(), alpha1, "laplacian(Dab,alpha1)")
alpha1Eqn.solve()
rhoPhi.ext_assign(alpha1Eqn.flux()*(rho1 - rho2) + phi*rho2)
rho.ext_assign(alpha1*rho1 + (1.0 - alpha1)*rho2)

from Foam.OpenFOAM import ext_Info, nl
ext_Info() << "Phase1 volume fraction = " << alpha1.weightedAverage( mesh.V() ).value() \
<< " Min(alpha1) = " << alpha1.ext_min().value() \
<< " Max(alpha1) = " << alpha1.ext_max().value() << nl
pass



It shows error on turbulence.nut(). Thanks in advance.

Hi Dhasthagheer,

Could you clarify, what kind of errors do you have? Have you tried to use "turbulence.nu()" instead of "turbulence.nut()". Also, could you print the "turbulence" object itself?

To speed up this investigation I suggest you upload your pythonFlu implementation of "twoLiquidMixingFoam" and your use-case to check on.

Best regards,
Alexey

How to solve alphaEqn using PythonFlu?
 

I have attached my case files and my python file. I am getting following error

Traceback (most recent call last):
File "twoLiquidMixingFoam.py", line 303, in <module>
os._exit( main_standalone( len( argv ), argv ) )
File "twoLiquidMixingFoam.py", line 57, in main_standalone
alphaEqn( mesh, phi, alpha1, rhoPhi, rho1, rho2, interface )
File "twoLiquidMixingFoam.py", line 87, in alphaEqn
alpha1Eqn = fvm.ddt(alpha1) + fvm.div(phi, alpha1) - fvm.laplacian( Dab + alphatab*turbulence.ext_nut(), alpha1, "laplacian(Dab,alpha1")
File "/usr/local/lib/python2.6/dist-packages/Foam/OpenFOAM/OpenFOAM_.py", line 3462, in __mul__
return _OpenFOAM_.dimensionedScalar___mul__(self, *args)
TypeError: in method 'dimensionedScalar___mul__', argument 2 of type 'Foam::dimensioned< Foam::scalar > const &'



The alphaEqn.H header file contains

{
fvScalarMatrix alpha1Eqn
(
fvm::ddt(alpha1)
+ fvm::div(phi, alpha1)
- fvm::laplacian
(
Dab + alphatab*turbulence->nut(), alpha1,
"laplacian(Dab,alpha1)"
)
);

alpha1Eqn.solve();

rhoPhi = alpha1Eqn.flux()*(rho1 - rho2) + phi*rho2;
rho = alpha1*rho1 + (scalar(1) - alpha1)*rho2;

Info<< "Phase 1 volume fraction = "
<< alpha1.weightedAverage(mesh.V()).value()
<< " Min(alpha1) = " << min(alpha1).value()
<< " Max(alpha1) = " << max(alpha1).value()
<< endl;
}

So I tried like following

def alphaEqn( mesh, phi, alpha1, rhoPhi, rho1, rho2, interface):
from Foam.finiteVolume import fvScalarMatrix
from Foam import fvm

alpha1Eqn = fvm.ddt(alpha1) + fvm.div(phi, alpha1) - fvm.laplacian( Dab + alphatab*turbulence.ext_nut(), alpha1, "laplacian(Dab,alpha1")
alpha1Eqn.solve()

rhoPhi.ext_assign(alpha1Eqn.flux()*(rho1 - rho2) + phi*rho2)
rho.ext_assign(alpha1*rho1 + (1.0 - alpha1)*rho2)

from Foam.OpenFOAM import ext_Info, nl
ext_Info() << "Phase1 volume fraction = " << alpha1.weightedAverage( mesh.V() ).value() \
<< " Min(alpha1) = " << alpha1.ext_min().value() \
<< " Max(alpha1) = " << alpha1.ext_max().value() << nl
pass

Thank you Dhasthagheer.

If we swap the two arguments in the target expression (I mean, "alphatab" and "turbulence.ext_nut()") like "turbulence.ext_nut() * alphatab", will it still work for you?

By the way, could you precise what version of OpenFOAM do you use? I have failed to run "blockMesh" for OpenFOAM-1.7.1-deb utility on the attached case, see
Best regards,
Alexey

How to solve alphaEqn using PythonFlu?
 

I am using Openfoam 1.7.1 .

blockMesh working fine for me in this case

$ blockMesh
/*---------------------------------------------------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.7.x |
| \\ / A nd | Web: www.OpenFOAM.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
Build : 1.7.x-4bbf33160caf
Exec : blockMesh
Date : May 15 2011
Time : 22:53:09
Host : Dhastha
PID : 3950
Case : /home/dhastha/New_PythonFlu/tutorial/twoLiquidMixingFoam/twoStream
nProcs : 1
SigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time


Creating block mesh from
"/home/dhastha/New_PythonFlu/tutorial/twoLiquidMixingFoam/twoStream/constant/polyMesh/blockMeshDict"


Creating blockCorners

Creating curved edges

Creating blocks

Creating patches

Creating block mesh topology

Default patch type set to empty
--> FOAM Warning :
From function polyMesh::polyMesh(... construct from shapes...)
in file meshes/polyMesh/polyMeshFromShapeMesh.C at line 576
Found 4 undefined faces in mesh; adding to default patch.

Check block mesh topology

Basic statistics
Number of internal faces : 1
Number of boundary faces : 10
Number of defined boundary faces : 10
Number of undefined boundary faces : 0

Checking patch -> block consistency

Creating block offsets

Creating merge list .

Creating points with scale 0.01

Creating cells

Creating patches

Creating mesh from block mesh

Default patch type set to empty

Writing polyMesh

End

How to solve alphaEqn using PythonFlu?
 

Hi

I swapped the two arguments in the target expression.

see the changes

alpha1Eqn = fvm.ddt(alpha1) + fvm.div(phi, alpha1) - fvm.laplacian( Dab + turbulence.ext_nut()*alphatab, alpha1, "laplacian(Dab,alpha1")

It shows following error

The initial problem with calculation of "alphatab * turbulence.ext_nut()" expression can be solved by the following modifications :

1. swapping of the "alphatab" and "turbulence.ext_nut()" arguments to make expression look "turbulence.ext_nut() * alphatab" (to be able to call already wrapped "volScalarFiled :: operator * ( const dimensionedScalar& )" functionality)
2. changing the call of "turbulence.ext_nut()" (which produces "tmp< volScalarField >") to "turbulence.nu()" (that gives bare "const volScalarField&")

This actually solves the initial problem. After this next appears, the following specialization of "laplacian" function is not wrapped yet :

tmp< scalarMatrix > laplacian( const tmp< volScalarField >& tgamma, volScalarField& vf, const word& name )

So, we got stack on the missed "laplacian" function :(

P.S.

Dhastha, thank you very much for your patience. pythonFlu is still on his way to complete wrapping. And this kind of "lacked" functionality is possible. The next pythonFlu release (which is actually planned in 2 weeks) will solve all the mentioned issues ...

Keep in touch,
Alexey

After some additional investigation the problem was successfully solved without patching the existing pythonFlu version.
In the attachment you can find as the new version of pythonFlu based solver, as well, as modified use-case to test it (there were some important differences with the referenced one).
So, if you are interested, you could find out the exact reason of the failure by comparison ours implementation with your.
If you made your modifications with some deep meaning, you could also use our implementation as a referenced point for your variation.

Good luck,
Alexey