[FlyBob91]

Hello to all,
i'm creating a new adjoint solver and i have to implement new boundary conditions.
i made all but i still have to add a dimensioned scalar constant to make the boundary condition correct, dimensionally talking.


i added these lines to the adjointOutletVelocityFvPatchVectorField.C

scalar Ke = readScalar(db().lookupObject<IOdictionary>("transportProperties").lookup("Ke"));

and i compiled correctly the solver.

Here my transportProperties file

transportModel  Newtonian;

nu              nu [0 2 -1 0 0 0 0] 1.02e-06;

Ke             Ke [0 7 -5 0 0 0 0] 1;

lambda          lambda   [0 -2 1 0 0 0 0] 1e5;
alphaMax        alphaMax [0 0 -1 0 0 0 0] 200.0;

Unfortutely i have the following error when i run my case

wrong token type - expected Scalar, found on line 22 the word 'Ke'

file: /home/roby/Scrivania/HAGAGA/constant/transportProperties.Ke at line 22.

    From function operator>>(Istream&, Scalar&)
    in file lnInclude/Scalar.C at line 93.

Any idea?

Thanks for help

[shereez234]

You are defining a scalar and giving it dimensions. If you want a dimensioned scalar then the correct way to declare it is:

const dimensioned<double> ke;

[FlyBob91]

Originally Posted by shereez234

You are defining a scalar and giving it dimensions. If you want a dimensioned scalar then the correct way to declare it is:

const dimensioned<double> ke;

Hi shereez,
thank you for your reply. I corrected my code and that line doesn't give me error anymore.
However another very strange error has appeared.
Here part of my code

   // GETTING NU
      const incompressible::RASModel& ras
    = db().lookupObject<incompressible::RASModel>("RASProperties");
    
    scalarField nuw
    = ras.nuEff()().boundaryField()[patch().index()];

    // DISTANCE^(-1)
    const scalarField& deltainv = patch().deltaCoeffs();

    // nu / h
    const scalarField& nd =  - nuw * deltainv;

    // NORMAL VECTOR
    vectorField n = patch().nf();

    // TANGENT VECTOR
    vectorField r = patch().Cf() - vector(10,11,12);

    // SCALAR NORMAL VELOCITY COMPONENT
    scalarField Un = Up & patch().nf();

    // DIMENSIONED CORRECTION COEFFICIENT

    const dimensioned<double> Ke = readScalar(db().lookupObject<IOdictionary>("transportProperties").lookup("Ke"));  

    // PATCH-ADIACENT TANGENT ADJOINT VELOCITY (VECTOR)
    vectorField Uac   = Uap.patchInternalField();
    vectorField Uac_n = ( Uac & patch().nf() ) * patch().nf();
    vectorField Uac_t = Uac - Uac_n;

    // TANGENT VELOCITY COMPONENT (VECTOR)
    vectorField U_n = ( Up & patch().nf() ) * patch().nf();
    vectorField U_t = Up - U_n;


    vectorField :: operator = (
                              phiap*patch().Sf()/sqr(patch().magSf())
                              + (1/(Un + nd))
                              * (nd * Uac_t + Ke *

 (n ^ r) / (pow(mag(r^U_t),3) + SMALL))
                              );

and here part of my error

/home/roby/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/lnInclude/GeometricFieldFunctionsM.C:119:55: note: template<template<class> class PatchField, class GeoMesh> Foam::tmp<Foam::GeometricField<Foam::Tensor<double>, PatchField, GeoMesh> > Foam::operator*(const Foam::GeometricField<Foam::Vector<double>, PatchField, GeoMesh>&)
 tmp<GeometricField<ReturnType, PatchField, GeoMesh> > operator Op             \
                                                       ^
/home/roby/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/lnInclude/GeometricTensorField.C:60:1: note: in expansion of macro ‘UNARY_OPERATOR’
 UNARY_OPERATOR(tensor, vector, *, hdual, transform)
 ^
/home/roby/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/lnInclude/GeometricFieldFunctionsM.C:119:55: note:   template argument deduction/substitution failed:
 tmp<GeometricField<ReturnType, PatchField, GeoMesh> > operator Op             \
                                                       ^
/home/roby/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/lnInclude/GeometricTensorField.C:60:1: note: in expansion of macro ‘UNARY_OPERATOR’
 UNARY_OPERATOR(tensor, vector, *, hdual, transform)
 ^
adjointOutletVelocity/adjointOutletVelocityFvPatchVectorField.C:143:58: note:   ‘const Foam::dimensioned<double>’ is not derived from ‘const Foam::GeometricField<Foam::Vector<double>, PatchField, GeoMesh>’
                               * (nd * Uac_t + Ke * (n ^ r) / (pow(mag(r^U_t),3) + SMALL))
                                                          ^
In file included from /home/roby/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/lnInclude/GeometricTensorField.C:30:0,
                 from /home/roby/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/lnInclude/GeometricTensorField.H:88,
                 from /home/roby/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/lnInclude/GeometricFields.H:36,
                 from /home/roby/OpenFOAM/OpenFOAM-2.3.0/src/finiteVolume/lnInclude/volFields.H:37,
                 from adjointOutletVelocity/adjointOutletVelocityFvPatchVectorField.C:27:
/home/roby/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/lnInclude/GeometricFieldFunctionsM.C:147:55: note: template<template<class> class PatchField, class GeoMesh> Foam::tmp<Foam::GeometricField<Foam::Tensor<double>, PatchField, GeoMesh> > Foam::operator*(const Foam::tmp<Foam::GeometricField<Foam::Vector<double>, PatchField, GeoMesh> >&)
 tmp<GeometricField<ReturnType, PatchField, GeoMesh> > operator Op             \
                                                       ^
/home/roby/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/lnInclude/GeometricTensorField.C:60:1: note: in expansion of macro ‘UNARY_OPERATOR’
 UNARY_OPERATOR(tensor, vector, *, hdual, transform)
 ^
/home/roby/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/lnInclude/GeometricFieldFunctionsM.C:147:55: note:   template argument deduction/substitution failed:
 tmp<GeometricField<ReturnType, PatchField, GeoMesh> > operator Op             \
                                                       ^
/home/roby/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/lnInclude/GeometricTensorField.C:60:1: note: in expansion of macro ‘UNARY_OPERATOR’
 UNARY_OPERATOR(tensor, vector, *, hdual, transform)
 ^
adjointOutletVelocity/adjointOutletVelocityFvPatchVectorField.C:143:58: note:   ‘const Foam::dimensioned<double>’ is not derived from ‘const Foam::tmp<Foam::GeometricField<Foam::Vector<double>, PatchField, GeoMesh> >’
                               * (nd * Uac_t + Ke * (n ^ r) / (pow(mag(r^U_t),3) + SMALL))

Any idea?

[rhythm]

Hi,

I'm not sure whether this is correct advice. But you are using operator ^ in multiple areas. As I understand you are trying to raise to the power, eg. n ^ r . But my guess is that this operator is not defined, and you should use function:

double pow(double x, double y)
x - base 
y - power

Feel free to correct me, if my interpretation of your code is incorrect.
Cheers,
Ben

[FlyBob91]

Originally Posted by rhythm

Hi,

I'm not sure whether this is correct advice. But you are using operator ^ in multiple areas. As I understand you are trying to raise to the power, eg. n ^ r . But my guess is that this operator is not defined, and you should use function:

double pow(double x, double y)
x - base 
y - power

Feel free to correct me, if my interpretation of your code is incorrect.
Cheers,
Ben

Hi rhythm, thank you for your reply.
It's strange because when i don't include the Ke term, the solvers compile correctly but when i run my case the solvers stops.
I firstly thought it was up to the fact the formula wasn't dimensionally correct, so i try to insert the Ke term to make it balanced (but i'm not sure this is the problem).

In the image the boudary condition i have to implement. I took it from a report so it should be correct.

n represents the normal vector of the outlet, while r the distance vector between the outlet center and the point considered.


Best Regards,
Roberto

[shereez234]

   // GETTING NU
      const incompressible::RASModel& ras
    = db().lookupObject<incompressible::RASModel>("RASProperties");
    
    scalarField nuw
    = ras.nuEff()().boundaryField()[patch().index()];

    // DISTANCE^(-1)
    const scalarField& deltainv = patch().deltaCoeffs();

    // nu / h
    const scalarField& nd =  - nuw * deltainv;

    // NORMAL VECTOR
    vectorField n = patch().nf();

    // TANGENT VECTOR
    vectorField r = patch().Cf() - vector(10,11,12);

    // SCALAR NORMAL VELOCITY COMPONENT
    scalarField Un = Up & patch().nf();

    // DIMENSIONED CORRECTION COEFFICIENT

    const dimensioned<double> Ke = readScalar(db().lookupObject<IOdictionary>("transportProperties").lookup("Ke"));  

    // PATCH-ADIACENT TANGENT ADJOINT VELOCITY (VECTOR)
    vectorField Uac   = Uap.patchInternalField();
    vectorField Uac_n = ( Uac & patch().nf() ) * patch().nf();
    vectorField Uac_t = Uac - Uac_n;

    // TANGENT VELOCITY COMPONENT (VECTOR)
    vectorField U_n = ( Up & patch().nf() ) * patch().nf();
    vectorField U_t = Up - U_n;


    vectorField :: operator = (
                              phiap*patch().Sf()/sqr(patch().magSf())
                              + (1/(Un + nd))
                              * (nd * Uac_t + Ke * pow(n, r) / (pow(mag(pow(r,U_t),3) + SMALL))
                              );

[rhythm]

Hi FlyBob91,

Ah.. My apologies I wasn't sure that the hat represented a cross product between vectorFields.
(P.S.: shereez234 sorry about me causing confusion in the thread. As FlyBob91 wants to perform cross product operation between n and r then his use of 'hat' operator is correct, and it is not a power operator)

I firstly thought it was up to the fact the formula wasn't dimensionally correct, so i try to insert the Ke term to make it balanced (but i'm not sure this is the problem).

If it was a dimensionality issue, then the code should explicitly say that it is the case.

As I am not entirely sure what is the issue, would you mind sending the code, for me to play around with and see if I can fix the issue?

Cheers,
Ben

[FlyBob91]

Originally Posted by shereez234

    vectorField :: operator = (
                              phiap*patch().Sf()/sqr(patch().magSf())
                              + (1/(Un + nd))
                              * (nd * Uac_t + Ke * pow(n, r) / (pow(mag(pow(r,U_t),3) + SMALL))
                              );

Hi shereez, thanks again for your reply.
I don't understand why you wrote

pow(n, r)

and

pow(r,U_t),3

.
I thought the cross product should be written with the operator ^

Best Regards,
Roberto

[FlyBob91]

Originally Posted by rhythm

As I am not entirely sure what is the issue, would you mind sending the code, for me to play around with and see if I can fix the issue?

Many thanks ben.


P.S. i'm working on OpenFoam 2.3.0

[shereez234]

Originally Posted by FlyBob91

Hi shereez, thanks again for your reply.
I don't understand why you wrote

pow(n, r)

and

pow(r,U_t),3

.
I thought the cross product should be written with the operator ^

Best Regards,
Roberto

That is correct. I am sorry I thought you were trying to use the power. Many Apologies

[rhythm]

Hi,

With very minor change I got it to compile. Try running it and see if it works.

There is a comment in the code explaining the change and why I believe it should be the correct way of doing it.

Cheers,
Ben

[FlyBob91]

Originally Posted by rhythm

Hi,

With very minor change I got it to compile. Try running it and see if it works.

There is a comment in the code explaining the change and why I believe it should be the correct way of doing it.

Cheers,
Ben

Hi Ben,
I have no words for thanking you and sorry for my late answer but i have been really busy in these days.

Your comment in the boundary condition file makes me think

do not carry dimension values with them

In fact i wanted to introduce the Ke term because i thought the boudary condition should be dimensionally balanced. In fact i wanted to equalize Ke to one and give the dimension necessary. I should have verificated it firstly.
So, trying with the old but good method of trial and error, i found the problematic therm was this one

const scalarField& nd =  - nuw * deltainv;

but i don't really understand why. Anyway i found the solution by decoupling it and writing nuw and deltainv directly in the boundary condition

 phiap*patch().Sf()/sqr(patch().magSf())
                              + (1/(Un - nuw * deltainv))
                              * (-nuw * deltainv * Uac_t + 2 * (n ^ r) / (pow(mag(r^U_t),3) + SMALL))
                              );

In this manner the solver compiled and the case run correctly.

However i have a last question in the following therm

// TANGENT VECTOR
    vectorField r = patch().Cf() - vector(10,11,12);

the vector

vector(10,11,12)

represents the outlet center. There is manner to extract it with some command without writing it directly?


Thank you another time for you help, you played a key role for the solution.


Best Regards,
Roberto

[rhythm]

Hi,

Regarding the dimensions of the boundary conditions.
I have had a discussion with my colleagues last year and it seems that in OpenFOAM the implementation is currently such that the boundary condition files (such as your adjointOutletPressureFvPatchScalarField) do not carry dimensions. So far I haven't seen any of the boundary conditions that are written differently. We presumed that later when the solver is running the Field just picks numerical values from the boundary condition and uses it's own set dimension values.

My guess is that this is done to avoid the need to have different boundary conditions for each type of field with different dimensions, or it might just be a less involved way to write it, but that is just a guess.

Anyway, when you are writing a boundary condition file, you only need to work with primitives like scalarFields, vectorFields, doubles, floats, etc. which you've been doing for the most part (apart from that Ke parameter).


Regarding the question about vector value, I'm afraid I don't have an answer, though I'm pretty sure there should be a function in the code that does that.

Cheers,
Ben