Inner product of tensor and vector

jfw_cfd · December 3, 2013, 03:11

Hello!

I try to implement an equation of scalar fluxes (an equation of type fvVectorMatrix) in OpenFOAM-1.6 where terms of an inner product of a tensor and a vector have to be calculated such as

$-\rho C_s u_j z \frac {\partial U_i}{\partial x_j} - \rho u_i u_j \frac {\partial Z}{\partial x_j}$

is a vector (volVectorField), so the first term is the inner product of a vector and a tensor

is the Reynolds stress tensor and Z is a scalar (volScalarField), so the second term is the inner product of a symmetric tensor and a vector
At the moment I wrote this in OpenFOAM as

Code:

 rho*Cs* uz & fvc::grad(U) + rho*turbulence->R() & fvc::grad(Z)

but it results in an error of " no match for 'operator==' ", actually this message

Code:

ZEqns.H:36: error: no match for ‘operator==’ in ‘Foam::operator-(const Foam::tmp<Foam::fvMatrix<Type> >&, const Foam::tmp<Foam::fvMatrix<Type> >&) [with Type = Foam::Vector<double>](((const Foam::tmp<Foam::fvMatrix<Foam::Vector<double> > >&)((const Foam::tmp<Foam::fvMatrix<Foam::Vector<double> > >*)(& Foam::fvm::laplacian(const Foam::tmp<Foam::GeometricField<GType, Foam::fvPatchField, Foam::volMesh> >&, Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>&) [with Type = Foam::Vector<double>, GType = Foam::SymmTensor<double>](((Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>&)(& uz))))))) == Foam::operator*(const Foam::tmp<Foam::GeometricField<double, PatchField, GeoMesh> >&, const Foam::tmp<Foam::GeometricField<Type, PatchField, GeoMesh> >&) [with Type = Foam::SymmTensor<double>, PatchField = Foam::fvPatchField, GeoMesh = Foam::volMesh](((const Foam::tmp<Foam::GeometricField<Foam::SymmTensor<double>, Foam::fvPatchField, Foam::volMesh> >&)((const Foam::tmp<Foam::GeometricField<Foam::SymmTensor<double>, Foam::fvPatchField, Foam::volMesh> >*)(& Foam::compressible::turbulenceModel::R()))))’
/opt/OpenFOAM/OpenFOAM-1.6/src/thermophysicalModels/specie/lnInclude/perfectGasI.H:148: note: candidates are: Foam::perfectGas Foam::operator==(const Foam::perfectGas&, const Foam::perfectGas&)
/opt/OpenFOAM/OpenFOAM-1.6/src/thermophysicalModels/specie/lnInclude/specieI.H:190: note:                 Foam::specie Foam::operator==(const Foam::specie&, const Foam::specie&)
/opt/OpenFOAM/OpenFOAM-1.6/src/OpenFOAM/lnInclude/triFaceI.H:304: note:                 bool Foam::operator==(const Foam::triFace&, const Foam::triFace&)
/opt/OpenFOAM/OpenFOAM-1.6/src/OpenFOAM/lnInclude/objectMapI.H:95: note:                 bool Foam::operator==(const Foam::objectMap&, const Foam::objectMap&)
/opt/OpenFOAM/OpenFOAM-1.6/src/OpenFOAM/lnInclude/objectHit.H:110: note:                 bool Foam::operator==(const Foam::objectHit&, const Foam::objectHit&)
/opt/OpenFOAM/OpenFOAM-1.6/src/OpenFOAM/lnInclude/cellShape.H:158: note:                 bool Foam::operator==(const Foam::cellShape&, const Foam::cellShape&)
/opt/OpenFOAM/OpenFOAM-1.6/src/OpenFOAM/lnInclude/cellModelI.H:124: note:                 bool Foam::operator==(const Foam::cellModel&, const Foam::cellModel&)
/opt/OpenFOAM/OpenFOAM-1.6/src/OpenFOAM/lnInclude/cell.H:134: note:                 bool Foam::operator==(const Foam::cell&, const Foam::cell&)
/opt/OpenFOAM/OpenFOAM-1.6/src/OpenFOAM/lnInclude/faceI.H:138: note:                 bool Foam::operator==(const Foam::face&, const Foam::face&)
/opt/OpenFOAM/OpenFOAM-1.6/src/OpenFOAM/lnInclude/edgeI.H:171: note:                 bool Foam::operator==(const Foam::edge&, const Foam::edge&)
/opt/OpenFOAM/OpenFOAM-1.6/src/OpenFOAM/lnInclude/instant.H:140: note:                 bool Foam::operator==(const Foam::instant&, const Foam::instant&)

Any ideas how to implement this inner product in OpenFOAM-1.6 ? The & operator should be the correct one for an inner product. Do I use the operator the wrong way or on wrong types or isn't it defined for these types?

Thanks in advance for any comment!

Greetings,
Johannes

Bernhard · December 3, 2013, 07:02

The error message you give does not relate to the posted code, because there is no "==" in there.

jfw_cfd · December 3, 2013, 07:11

Hello Bernhard!

Thank you for your reply.
My fault, since I didn't post the whole equation, only two terms of it.

Anyway, I guess I found a solution with which I'll come back a bit later to get an opinion from others.

Greetings and let's "read" again in a while....

jfw_cfd · December 3, 2013, 13:03

As I sad before, I'll come back with a solution which I would like to share with you.

I finally implemented this equation

${\partial \rho u'_iz' \over \partial t} + {\partial \rho u'_iz' U_j \over \partial x_j} - {\partial \over \partial x_j} \left(C_s \rho \frac{k}{\varepsilon} u'_jz' {\partial U_i \over \partial x_j} \right)$
=
$\rho (C_2-1) u'_j z' \frac {\partial U_i}{\partial x_j} - \rho u'_i u'_j \frac {\partial Z}{\partial x_j} - \rho C_1 \frac{\varepsilon}{k} u'_iz'$
where

is the vector (volVectorField) which has to be solved
in OpenFOAM-1.6 as

Code:

    volVectorField R_Z = turbulence->R() & fvc::grad(Z);
    volVectorField uz_gU = uz & fvc::grad(U);

    fvVectorMatrix uzEqn
    (
        fvm::ddt(rho, uz)
      + fvm::div(phi, uz)
      - fvm::laplacian(rho*C_s*turbulence->k()/turbulence->epsilon()*turbulence->R(), uz)
      ==
        rho*(C_2 - 1)*uz_gU
      - rho*R_Z
      - rho*C_1*turbulence->epsilon()/turbulence->k()*uz
    );

The point is, that one has to calculate the two inner products of tensor and vector (which give volVectorFields) outside or before the actual equation. It compiles and runs, but not very stable. I use it together with the LRR turbulence model, but not with the default diffusion term (as explained here http://www.cfd-online.com/Forums/ope...tml#post463684 ).

Any comments? Is there a better way or more "stable form" to implement the equation?

Thank you and kind regards,
Johannes

gxy200992243 · July 7, 2014, 10:04

Hi, Johannes

Have you got the solution to this problem? I am also facing the inner product of a tensor and a vector. Can you post the solution to this problem?

Best regards,

Xiangyu

jfw_cfd · July 7, 2014, 10:33

Hello Xiangyu!

Well, the solution of how to calculate the inner product of a tensor and a vector which I stated in my previous post should actually work (I guess you can see the code snippet). At least it worked for me (compiled and calculated). The instabilities I mentioned above, did not arise from this part of the code as far as I remember (all this happened quite some months ago...).

Unfortunately, I don't have any further news about this topic. So I hope the above code works for you, which should be the case if just the inner product is the problem.

Kind regards,
Johannes

gxy200992243 · July 7, 2014, 10:36

Hi, Johannes!

Thank you very much for your reply! I just tried. Your solution works for my case perfectly! Thank you again for your help!

Best regards,

Xiangyu

jfw_cfd · July 7, 2014, 10:41

Perfect! That's good to hear, Xiangyu!

parthiv1991 · November 23, 2017, 08:39

Hallo Johannes,

I am trying to implement a Rosenbrock-Wanner method in OpenFOAM for incompressible flows by modifying the icoFoam solver.

At one part in my code I have to multiply a volTensorField with volVectorField (i.e. a dot product), for which I use the "&" symbol in OpenFOAM. The code compiles without errors but while running a case it gives the "segmentation fault (core dumped)" error.

Here's the code snippet:

Code:

phi = (fvc::interpolate(U)& mesh.Sf());
rhs.set(m,new volVectorField(runTime.deltaT()*(-fvc::div(phi,U)+ nu * fvc::laplacian(U))));
		
//the next line is where the error is being generated	

    volVectorField slope = kihelp & (st-(g*runTime.deltaT()*symm(J)));
            
			fvVectorMatrix slopeE
			(
			  (fvm::Sp(1,slope))
			    ==
			    rhs[m]
			 );
			 slopeE.solve();

where: st = Identity tensor, kihelp = function of the velocity(volVectorField), J = Jacobian of velocity field (volTensorField).

So the equations I want to solve are in the form of Ax=b,
where: matrix A = Identity - (scalar_const * delta_T * Jacobian)
vector x = kihelp
rhs b = function of Velocity.

And here is how I declared the Identity matrix/tensor :

Code:

GeometricField<symmTensor, fvPatchField, volMesh> st
    (
        IOobject
        (
            "st",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh,
        dimensioned<symmTensor>("st", dimless, symmTensor::one)
        zeroGradientFvPatchSymmTensorField::typeName
    );

Could you or anyone else please suggest me how to go about this?
Any help is appreciated.

Thanks in advance

vitor.geraldes@ist.utl.pt · January 21, 2018, 20:44

I am also facing the same problem. I have tried run this code:

#include "vector.H"
#include "tensor.H"
#include "IOstreams.H"

using namespace Foam;

int main()
{
tensor T(1, 0, 0, 0, 1, 0, 0, 0, 1);
vector d(1, 1, 1);
tensor TT = T*d;

Info<< (TT) << endl;

return 0;
}

but, the compiler displays an error when "tensor TT = T*d;"

According to the programming manual, this should work.

jfw_cfd · January 22, 2018, 04:14

Hi Vitor!

First of all, I have to admit that (unfortunately) I rather lost my close contact with OpenFOAM about two years ago. So I'm probably not the most adequate person to answer your (and Parthiv's) question.

However, when looking at what you want to do and as far as my very limited mathematical background allows, I think a dot product of a tensor and a vector gives you a vector. So shouldn't TT be of type vector?
Side note: T.d is not equal d.T

Maybe this gives you an idea how to go on with your problem.

To Parthiv: Sorry for this very late reply, but I'm not able to give you any solution at the moment. Maybe you fixed that already. If so, please post the solution.

Good luck guys!

December 3, 2013, 13:03	possible solution	#4
jfw_cfd New Member Johannes Join Date: Mar 2011 Location: Austria Posts: 12 Rep Power: 15	As I sad before, I'll come back with a solution which I would like to share with you. I finally implemented this equation ${\partial \rho u'_iz' \over \partial t} + {\partial \rho u'_iz' U_j \over \partial x_j} - {\partial \over \partial x_j} \left(C_s \rho \frac{k}{\varepsilon} u'_jz' {\partial U_i \over \partial x_j} \right)$ = $\rho (C_2-1) u'_j z' \frac {\partial U_i}{\partial x_j} - \rho u'_i u'_j \frac {\partial Z}{\partial x_j} - \rho C_1 \frac{\varepsilon}{k} u'_iz'$ where is the vector (volVectorField) which has to be solved in OpenFOAM-1.6 as Code: volVectorField R_Z = turbulence->R() & fvc::grad(Z); volVectorField uz_gU = uz & fvc::grad(U); fvVectorMatrix uzEqn ( fvm::ddt(rho, uz) + fvm::div(phi, uz) - fvm::laplacian(rhoC_sturbulence->k()/turbulence->epsilon()turbulence->R(), uz) == rho(C_2 - 1)uz_gU - rhoR_Z - rhoC_1turbulence->epsilon()/turbulence->k()uz ); The point is, that one has to calculate the two inner products of tensor and vector (which give volVectorFields) outside or before the actual equation. It compiles and runs, but not very stable. I use it together with the LRR turbulence model, but not with the default diffusion term (as explained here http://www.cfd-online.com/Forums/ope...tml#post463684 ). Any comments? Is there a better way or more "stable form" to implement the equation? Thank you and kind regards, Johannes Last edited by jfw_cfd; December 4, 2013 at 11:59.*

July 7, 2014, 10:04		#5
gxy200992243 New Member Xiangyu Gao Join Date: Sep 2013 Location: West Lafayette, IN, USA Posts: 29 Rep Power: 12	Hi, Johannes Have you got the solution to this problem? I am also facing the inner product of a tensor and a vector. Can you post the solution to this problem? Best regards, Xiangyu vitor.geraldes@ist.utl.pt likes this.

July 7, 2014, 10:33	proposed solution should work	#6
jfw_cfd New Member Johannes Join Date: Mar 2011 Location: Austria Posts: 12 Rep Power: 15	Hello Xiangyu! Well, the solution of how to calculate the inner product of a tensor and a vector which I stated in my previous post should actually work (I guess you can see the code snippet). At least it worked for me (compiled and calculated). The instabilities I mentioned above, did not arise from this part of the code as far as I remember (all this happened quite some months ago...). Unfortunately, I don't have any further news about this topic. So I hope the above code works for you, which should be the case if just the inner product is the problem. Kind regards, Johannes

November 23, 2017, 08:39		#9
parthiv1991 New Member PLD Join Date: Jun 2017 Location: Braunschweig, Germany Posts: 13 Rep Power: 8	Hallo Johannes, I am trying to implement a Rosenbrock-Wanner method in OpenFOAM for incompressible flows by modifying the icoFoam solver. At one part in my code I have to multiply a volTensorField with volVectorField (i.e. a dot product), for which I use the "&" symbol in OpenFOAM. The code compiles without errors but while running a case it gives the "segmentation fault (core dumped)" error. Here's the code snippet: Code: phi = (fvc::interpolate(U)& mesh.Sf()); rhs.set(m,new volVectorField(runTime.deltaT()(-fvc::div(phi,U)+ nu fvc::laplacian(U)))); //the next line is where the error is being generated volVectorField slope = kihelp & (st-(grunTime.deltaT()symm(J))); fvVectorMatrix slopeE ( (fvm::Sp(1,slope)) == rhs[m] ); slopeE.solve(); where: st = Identity tensor, kihelp = function of the velocity(volVectorField), J = Jacobian of velocity field (volTensorField). So the equations I want to solve are in the form of Ax=b, where: matrix A = Identity - (scalar_const * delta_T * Jacobian) vector x = kihelp rhs b = function of Velocity. And here is how I declared the Identity matrix/tensor : Code: GeometricField<symmTensor, fvPatchField, volMesh> st ( IOobject ( "st", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::NO_WRITE ), mesh, dimensioned<symmTensor>("st", dimless, symmTensor::one) zeroGradientFvPatchSymmTensorField::typeName ); Could you or anyone else please suggest me how to go about this? Any help is appreciated. Thanks in advance

January 22, 2018, 04:14		#11
jfw_cfd New Member Johannes Join Date: Mar 2011 Location: Austria Posts: 12 Rep Power: 15	Hi Vitor! First of all, I have to admit that (unfortunately) I rather lost my close contact with OpenFOAM about two years ago. So I'm probably not the most adequate person to answer your (and Parthiv's) question. However, when looking at what you want to do and as far as my very limited mathematical background allows, I think a dot product of a tensor and a vector gives you a vector. So shouldn't TT be of type vector? Side note: T.d is not equal d.T Maybe this gives you an idea how to go on with your problem. To Parthiv: Sorry for this very late reply, but I'm not able to give you any solution at the moment. Maybe you fixed that already. If so, please post the solution. Good luck guys! vitor.geraldes@ist.utl.pt likes this.

December 3, 2013, 07:02		#2
Bernhard Senior Member Bernhard Join Date: Sep 2009 Location: Delft Posts: 790 Rep Power: 21	The error message you give does not relate to the posted code, because there is no "==" in there.

December 3, 2013, 07:11		#3
jfw_cfd New Member Johannes Join Date: Mar 2011 Location: Austria Posts: 12 Rep Power: 15	Hello Bernhard! Thank you for your reply. My fault, since I didn't post the whole equation, only two terms of it. Anyway, I guess I found a solution with which I'll come back a bit later to get an opinion from others. Greetings and let's "read" again in a while....

July 7, 2014, 10:36		#7
gxy200992243 New Member Xiangyu Gao Join Date: Sep 2013 Location: West Lafayette, IN, USA Posts: 29 Rep Power: 12	Hi, Johannes! Thank you very much for your reply! I just tried. Your solution works for my case perfectly! Thank you again for your help! Best regards, Xiangyu

July 7, 2014, 10:41		#8
jfw_cfd New Member Johannes Join Date: Mar 2011 Location: Austria Posts: 12 Rep Power: 15	Perfect! That's good to hear, Xiangyu!

January 21, 2018, 20:44		#10
vitor.geraldes@ist.utl.pt New Member Vitor Geraldes Join Date: Dec 2009 Location: Lisbon, Portugal Posts: 26 Rep Power: 16	I am also facing the same problem. I have tried run this code: #include "vector.H" #include "tensor.H" #include "IOstreams.H" using namespace Foam; int main() { tensor T(1, 0, 0, 0, 1, 0, 0, 0, 1); vector d(1, 1, 1); tensor TT = Td; Info<< (TT) << endl; return 0; } but, the compiler displays an error when "tensor TT = Td;" According to the programming manual, this should work.