[Diego_Bindoni]

Good morning,
i'm writing you to ask for help on the coding environment in OpenFOAMv12.
I need to set up a vertical profile of temperature inside my 3D simulations using setFields, but i struggle into coding it. It doesn't exactly look like when you code for creating new applications or boundary conditions, and the tutorials doesn't give me help into calling the cell points in my code.
Could you please give me help through this?
Thanks in advance

Up to know, i wrote something like this (which clearly doesn't work )


regions
(
boxToCell
{
box (0-1530) (90150);
fieldValues #codeStream
(
codeInclude
#{
#include"pointField.H"
#};

code
#{
const scalar T_min =300; // Minimum temperature (K)
const scalar T_max =500; // Maximum temperature (K)
const scalar z_min =0; // Minimum z-coordinate
const scalar z_max =30; // Maximum z-coordinate

const scalar Tz = T_min + (T_max - T_min) * (pos.z() - z_min) / (z_max - z_min);
os <<"volScalarFieldValue T "<< Tz <<";";
#};
);
}
);

[Tobermory]

That sounds a bit ambitious ... I don't think that setFields will allow that (I could be wrong, but looking at the code I don't think so). You could, however, insert the code in the boundary file (0/T), for the internalField.

[Diego_Bindoni]

Originally Posted by Tobermory

That sounds a bit ambitious ... I don't think that setFields will allow that (I could be wrong, but looking at the code I don't think so). You could, however, insert the code in the boundary file (0/T), for the internalField.

Thanks for the reply. In this case, how should i proceed? I'm really a bit stuck with how to code it practically. Thanks in advance for your time

[Tobermory]

The following video may help:

https://www.youtube.com/watch?v=VDpLpUV1X-k

[Diego_Bindoni]

Originally Posted by Tobermory

The following video may help:

https://www.youtube.com/watch?v=VDpLpUV1X-k

Thanks a lot! I think it could help. I'll get back to this thread if everything went fine (or wrong).

[Diego_Bindoni]

Originally Posted by Diego_Bindoni

Thanks a lot! I think it could help. I'll get back to this thread if everything went fine (or wrong).

Hello, unfortunately it didn't solve the problem. I put here the piece of code related to the internalfield modification that i made, but when i try to run, i have this error

--> FOAM FATAL IO ERROR:
expected keyword 'uniform' or 'nonuniform', found codeInclude

file: /home/dbindoni/OpenFOAM/dbindoni-12/run/AMR/jet_BAMRnut/0/T from line 16 to line 39.

From function Foam::Field<Type>::Field(const Foam::word&, const Foam::unitConversion&, const Foam::dictionary&, Foam::label) [with Type = double; Foam::label = int]
in file /home/ubuntu/OpenFOAM/OpenFOAM-12/src/OpenFOAM/lnInclude/Field.C at line 237.

FOAM exiting

I put here the code

FoamFile
{
formatascii;
classvolScalarField;
objectT;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions [0001000];

internalField #codeStream
{
{
codeInclude
#{
#include "fvCFD.H"
#};

codeOptions
#{
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude
#};

codeLibs
#{
-lfiniteVolume \
-lmeshTools
#};

code
#{
constIOdictionary&d+static_cast<constIOdictionary&>(dict);
constfvMesh&mesh=refCast<constfvMesh>(d.db());

scalarFieldT(mesh.nCells(), 0.);

scalarTup=303;
scalarTbot=273;
scalarhmid=15;

forAll(T, i)
{
constscalarz=mesh.C()[i][2];

if (z<=hmid)
{
T[i] =Tup;
}
else
{
T[i] =Tbot;
}
}

writeEntry(os, "", T);
#};
};
}

[Tobermory]

Is your writeFormat in system/controlDict set to binary? If so, change to ascii and see if that fixes the problem - others on the forum have reported the same error (OF9, Error after #codeStream internalField).


EDIT: the reason for the "expected keyword 'uniform'" error in the example above is that you have too many brackets. You have

internalField #codeStream 
{
  {
     codeInclude etc...
  };
};

but should have instead:

internalField #codeStream 
{
   codeInclude etc...
};

This fixes the original problem ... but the codeStream will still crash since OF12 has tweaked the headers - see my post below for the solution.

[Diego_Bindoni]

Originally Posted by Tobermory

Is your writeFormat in system/controlDict set to binary? If so, change to ascii and see if that fixes the problem - others on the forum have reported the same error (OF9, Error after #codeStream internalField).

Yes i checked that forum and my writeFormat is in ascii. I just think it's not how openfoamv12 accepts code lines now i guess. But there is little to no information for a case like this one, even though it shouldn't be such an exotic condition to have...

[Diego_Bindoni]

Originally Posted by Diego_Bindoni

Yes i checked that forum and my writeFormat is in ascii. I just think it's not how openfoamv12 accepts code lines now i guess. But there is little to no information for a case like this one, even though it shouldn't be such an exotic condition to have...

I partially solved the problem. OpenFOAM allows you to do a coded setFields only through a coded functionObject, that you run with foamPostProcess (ironic) before the simulation starts.


As a reference, in OpenFOAMv12 the case blockedChannel has a generateAlphas function Object that does it. With this you obtain a T field that now i have make it look like a boundary condition


Hope it will be useful for everyone

[Tobermory]

Apologies Diego - OF12 has changed a bunch of the headers ... including the fvCFD.H header that you included. To get Joel's coding from the YT clip to work in OF12 you need something like:

internalField   #codeStream
                {
                    codeInclude
                    #{
                        #include "volFields.H"
                    #};

                    codeOptions
                    #{
                        -I$(LIB_SRC)/finiteVolume/lnInclude \
                        -I$(LIB_SRC)/meshTools/lnInclude
                    #};

                    codeLibs
                    #{
                        -lfiniteVolume \
                        -lmeshTools
                    #};

                    code
                    #{
                        const IOdictionary& d = static_cast<const IOdictionary&>(dict);
                        const fvMesh& mesh = refCast<const fvMesh>(d.db());

                        scalarField T(mesh().nCells(), 0.);
                        scalar Tup = 300;
                        scalar Tbot = 270;
                        scalar hmid = 0;

                        forAll(T, i)
                        {
                            const scalar z = mesh.C()[i][2];

                            if (z <= hmid)
                            {
                                T[i] = Tup;
                            }
                            else
                            {
                                T[i] = Tbot;
                            }
                        }
                        writeEntry(os, "", T);
                    #};
                };

i.e. different include files, and ditch the extra brackets at the start and end of the codestream chunk. This compiles and works fine for me.

[Diego_Bindoni]

Originally Posted by Tobermory

Apologies Diego - OF12 has changed a bunch of the headers ... including the fvCFD.H header that you included. To get Joel's coding from the YT clip to work in OF12 you need something like:

internalField   #codeStream
                {
                    codeInclude
                    #{
                        #include "volFields.H"
                    #};

                    codeOptions
                    #{
                        -I$(LIB_SRC)/finiteVolume/lnInclude \
                        -I$(LIB_SRC)/meshTools/lnInclude
                    #};

                    codeLibs
                    #{
                        -lfiniteVolume \
                        -lmeshTools
                    #};

                    code
                    #{
                        const IOdictionary& d = static_cast<const IOdictionary&>(dict);
                        const fvMesh& mesh = refCast<const fvMesh>(d.db());

                        scalarField T(mesh().nCells(), 0.);
                        scalar Tup = 300;
                        scalar Tbot = 270;
                        scalar hmid = 0;

                        forAll(T, i)
                        {
                            const scalar z = mesh.C()[i][2];

                            if (z <= hmid)
                            {
                                T[i] = Tup;
                            }
                            else
                            {
                                T[i] = Tbot;
                            }
                        }
                        writeEntry(os, "", T);
                    #};
                };

i.e. different include files, and ditch the extra brackets at the start and end of the codestream chunk. This compiles and works fine for me.

Thanks a lot for the help. I tried in my case and it runs only if i have hmid as 0. If i change it in value (like 15, since my domain goes from 0 to 30. It crashes with this error

--> FOAM FATAL IO ERROR:
compound has already been transferred from token
on line 0 the empty compound of type List<scalar>

[Tobermory]

Hmm - that's an odd message. Can you post your 0/T file in full, so that I can see what you've written?

[Diego_Bindoni]

Originally Posted by Tobermory

Hmm - that's an odd message. Can you post your 0/T file in full, so that I can see what you've written?

/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 11
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
formatascii;
classvolScalarField;
objectT;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions [0001000];

internalField #codeStream
{
codeInclude
#{
#include "volFields.H"
#};

codeOptions
#{
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude
#};

codeLibs
#{
-lfiniteVolume \
-lmeshTools
#};

code
#{
constIOdictionary&d=static_cast<constIOdictionary&>(dict);
constfvMesh&mesh=refCast<constfvMesh>(d.db());

scalarFieldT(mesh().nCells(), 0.);
scalarTup=300;
scalarTbot=270;
scalarhmid=5;

forAll(T, i)
{
constscalarz=mesh.C()[i][2];

if (z<=hmid)
{
T[i] =Tup;
}
else
{
T[i] =Tbot;
}
}
writeEntry(os, "", T);
#};
};
boundaryField
{
#includeEtc "caseDicts/setConstraintTypes"

jet_inlet
{
typefixedValue;
valueuniform303;
}

outlet
{
typeinletOutlet;
inletValue $internalField;
value $internalField;
}

top
{
typezeroGradient;
}

bottom
{
typezeroGradient;
}


inlet
{
typezeroGradient;
}


inlet_sides
{
typezeroGradient;
}

sides
{
typeinletOutlet;
inletValue $internalField;
value $internalField;
}

}

[Diego_Bindoni]

Originally Posted by Tobermory

Hmm - that's an odd message. Can you post your 0/T file in full, so that I can see what you've written?

The problem seems to be due to the presence of the inletoutlet condition. it struggle into dynamically linking it with the internalfield coding. Changing it with zerogradient solve th problem, but i would like to keep the inlet outlet condition to still maintain possibility for inflow. I will think about how to keep it, do you have any idea or recommendation?

[Tobermory]

Yes, I thought that might be it. The simplest (but not elegant) way would be to insert the same codeStream coding for the inletValue of the inletOutlet, and stop referring to $internalField.

[Diego_Bindoni]

Originally Posted by Tobermory

Yes, I thought that might be it. The simplest (but not elegant) way would be to insert the same codeStream coding for the inletValue of the inletOutlet, and stop referring to $internalField.

I will add a doubt to the thread. I'm running simulations o a round jet in marine environment, which explains the need for a vertical temperature gradient. The jet, at a depth of 15 meters, injects water at surface temperature, thus creating a positive buoyancy. Since i've added the internalField coded condition however, the temperature is always as the gradient i've provided, thus not allowing any buoyancy effect. I can see the velocity field generated by the jet, but not the associated temperature field generated by it. Do you have any idea of why?

[Tobermory]

If the temperature field is not updating, then check the log file - is an energy equation being solved? is it doing any iterations? are there any warning messages? what solver are you running? Need some more info to understand what's going on.

[Diego_Bindoni]

Originally Posted by Tobermory

If the temperature field is not updating, then check the log file - is an energy equation being solved? is it doing any iterations? are there any warning messages? what solver are you running? Need some more info to understand what's going on.

Ok, i will give a bit more informations about my setup. I attach down a normal print of my log. file. I'm using the solver fluid. No warning messages for what i could read, even at the beginning


Courant Number mean: 0.00117763 max: 0.898341
deltaT = 0.0325333
Time = 827.259251s

diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
smoothSolver: Solving for Ux, Initial residual = 8.89797e-05, Final residual = 1.97563e-07, No Iterations 1
smoothSolver: Solving for Uy, Initial residual = 0.00079219, Final residual = 5.92003e-08, No Iterations 2
smoothSolver: Solving for Uz, Initial residual = 0.00015151, Final residual = 2.89056e-07, No Iterations 1
smoothSolver: Solving for h, Initial residual = 5.07307e-05, Final residual = 1.90883e-07, No Iterations 1
GAMG: Solving for p_rgh, Initial residual = 0.0301946, Final residual = 0.00231544, No Iterations 1
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 7.75695e-09, global = -5.62662e-13, cumulative = 1.37926e-10
GAMG: Solving for p_rgh, Initial residual = 0.00257173, Final residual = 9.14075e-07, No Iterations 10
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 3.06215e-12, global = -3.83321e-14, cumulative = 1.37888e-10
smoothSolver: Solving for k, Initial residual = 0.000132324, Final residual = 7.92914e-07, No Iterations 1
bounding k, min: -1.15002e-05 max: 0.046884 average: 0.000105664
ExecutionTime = 89741.5 s ClockTime = 89759 s

scalarTransport write:
smoothSolver: Solving for s, Initial residual = 4.77571e-05, Final residual = 1.74934e-07, No Iterations 1

[Tobermory]

Okay, that's looking good. So the attention now turns to your model setup and boundary conditions. What equation of state are you using? Can you share your physicalProperties file?

[Diego_Bindoni]

Originally Posted by Tobermory

Okay, that's looking good. So the attention now turns to your model setup and boundary conditions. What equation of state are you using? Can you share your physicalProperties file?

I did some other checks. The equation of state seems to work perfectly fine by activating or deactivating the non uniform internalField (it's a polynomial equation of state btw).

my internalField codeStream is attached below.

                        internalField #codeStream
  {
  codeInclude
  #{
  #include "volFields.H"
  #};
 
 
  codeOptions
  #{
  -I$(LIB_SRC)/finiteVolume/lnInclude \
  -I$(LIB_SRC)/meshTools/lnInclude
  #};
 
 
  code
  #{
  const IOdictionary& d = static_cast<const IOdictionary&>(dict);
  const fvMesh& mesh = refCast<const fvMesh>(d.db());
 

 
  scalarField T(mesh().nCells(), 0.);
  scalar Tup = 300;
  scalar Tbot = 270;
  scalar hmid = 3;
  scalar zup = 0;
  forAll(T, i)
  {
  const scalar z = mesh.C()[i][2];
 
 
  if (z <= hmid)
  {
  T[i] = Tup;
  // T[i] = Tbot + (Tup - Tbot) * (hmid - z) / (hmid - zup);
  // T[i] = Tbot + (Tup - Tbot) * pow((hmid - z) / (hmid - zup), 0.5);
  }
  else
  {
  T[i] = Tbot;
  }
  }
 
 
  writeEntry(os, "", T);
  #};
  };
  p { line-height: 115%; margin-bottom: 0.1in; background: transparent }

The problem is that if i leave the if structure like this, i have two slices at different temperature and a correct diffusion of the temperature from the jet. If instead i adopt the commented options (a linear or a quadratic trend in the temperature) i have the non uniform internalfield, but from the jet no temperature is seen diffusing in the domain (i just have the temperature fixedValue boundary condition on the patch). I can't wrap my head around the code error i put in it. I've also tried to make several if statements to create more slices, but with just another else if i go back to the forced temperature profile, with no diffusion of temperature).

P.S. I've improperly used the term diffusion in this message, but i hope the message still passed.


UPDATE: I've also tried another way around to trick OpenFOAM, buy launching an identical simulation with the jet turned off that could give me a non uniform temperature field, and then i have used mapFields on the original case. Unfortunately it didn't work, which make things even more weird.