CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Community Contributions (https://www.cfd-online.com/Forums/openfoam-community-contributions/)
-   -   [pythonFlu] How to solve alphaEqn using PythonFlu? (https://www.cfd-online.com/Forums/openfoam-community-contributions/88347-how-solve-alphaeqn-using-pythonflu.html)

 dhasthagheer May 14, 2011 16:00

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.

 alexey2petrov May 15, 2011 08:25

Hi Dhasthagheer,

Quote:
 Originally Posted by dhasthagheer (Post 307645) I am writing twoliquidMixingFoam solver in python using PythonFlu ... It shows error on turbulence.nut()
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

 dhasthagheer May 15, 2011 09:33

How to solve alphaEqn using PythonFlu?

1 Attachment(s)
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

 alexey2petrov May 15, 2011 10:38

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
Code:

```/*---------------------------------------------------------------------------*\ | =========                |                                                | | \\      /  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-131caa989cd3 Exec  : blockMesh Date  : May 15 2011 Time  : 18:21:09 Host  : maverik-amd64 PID    : 25957 Case  : /home/alexey/tmp nProcs : 1 SigFpe : Enabling floating point exception trapping (FOAM_SIGFPE). // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Create time Creating block mesh from     "/home/alexey/tmp/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 --> FOAM FATAL ERROR: Failed writing polyMesh.     From function blockMesh     in file blockMeshApp.C at line 345. FOAM exiting```
Best regards,
Alexey

 dhasthagheer May 15, 2011 13:26

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

 dhasthagheer May 15, 2011 13:40

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

Code:

```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 + turbulence.ext_nut()*alphatab, alpha1, "laplacian(Dab,alpha1") TypeError: unsupported operand type(s) for *: 'ext_tmp_volScalarField' and 'dimensionedScalar```

 alexey2petrov May 15, 2011 16:24

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

 alexey2petrov May 19, 2011 05:08

1 Attachment(s)
Quote:
 Originally Posted by alexey2petrov (Post 307706) So, we got stack ...
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

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