Modified laplacian solver

stockzahn · February 6, 2019, 05:15

Dear all,

I tried to modify laplacianFoam to

- be able to use temperature dependent material properties and

- obtain the generated heat flux at the boundaries.

I'm sure there already exist many solutions, function objects etc. to get this done, but for exercise I wanted to try that by my own. My code compiles and the solver in principle works. But during coding some questiosn came up I wanted to ask. But to start with the codes of the three modified files:

laplacianFoam.C - modified

Code:

#include "fvCFD.H"
#include "fvOptions.H"
#include "simpleControl.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    #include "setRootCaseLists.H"

    #include "createTime.H"
    #include "createMesh.H"

    simpleControl simple(mesh);
    bool varProps = false;

    #include "createFields.H"

    //First try: code snippet example to read information from controlDict - not working
    if (runTime.controlDict().lookup("variable") == "yes" || runTime.controlDict().lookup("variable") == "true")
    {
        varProps = true;
    } 

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    Info<< "\nCalculating temperature distribution\n" << endl;

    while (simple.loop(runTime))
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        while (simple.correctNonOrthogonal())
        {
            if(varProps == true)    //variable material properties
            {
                volScalarField k = c0_k * pow(T,0) + c1_k * pow(T,1) + c2_k * pow(T,2) + c3_k * pow(T,3) + c4_k * pow(T,4) + c5_k * pow(T,5) + c6_k * pow(T,6) + c7_k * pow(T,7);

                volScalarField c = c0_c * pow(T,0) + c1_c * pow(T,1) + c2_c * pow(T,2) + c3_c * pow(T,3) + c4_c * pow(T,4) + c5_c * pow(T,5) + c6_c * pow(T,6) + c7_c * pow(T,7);

                volScalarField rho = c0_rho * pow(T,0) + c1_rho * pow(T,1) + c2_rho * pow(T,2) + c3_rho * pow(T,3) + c4_rho * pow(T,4) + c5_rho * pow(T,5) + c6_rho * pow(T,6) + c7_rho * pow(T,7);

                volScalarField thermDiff = k / (c * rho);

                fvScalarMatrix TEqn
                (
                    fvm::ddt(T) - fvm::laplacian(thermDiff, T)
                    ==
                    fvOptions(T)
                );
                fvOptions.constrain(TEqn);
                TEqn.solve();
                    fvOptions.correct(T);
            }
            else    //constant material properties
            {
                fvScalarMatrix TEqn
                (
                    fvm::ddt(T) - fvm::laplacian(DT, T)
                    ==
                    fvOptions(T)
                );
                fvOptions.constrain(TEqn);
                TEqn.solve();
                fvOptions.correct(T);
            }
        }

        #include "write.H"

        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }

    Info<< "End\n" << endl;

    return 0;
}

createFields.H - modified

Code:

Info<< "Reading field T\n" << endl;

volScalarField T
(
    IOobject
    (
        "T",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

Info<< "Reading transportProperties\n" << endl;

IOdictionary transportProperties
(
    IOobject
    (
        "transportProperties",
        runTime.constant(),
        mesh,
        IOobject::MUST_READ_IF_MODIFIED,
        IOobject::NO_WRITE
    )
);

//Second try: code snippet example to read information from transportDict - not working 
if(transportProperties.lookup("variable")=="yes" || transportProperties.lookup("variable")=="true") 
{
    varProps = true;
} 

Info<< "Determining material properties" << nl << endl;
Info<< "Material properties variable (phi = f(T))" 
<< "    Reading polynomial coefficients for thermal conductivity" << endl;

dimensionedScalar c0_k (transportProperties.lookup("c0_k"));    
dimensionedScalar c1_k (transportProperties.lookup("c1_k"));    
dimensionedScalar c2_k (transportProperties.lookup("c2_k"));    
dimensionedScalar c3_k (transportProperties.lookup("c3_k"));    
dimensionedScalar c4_k (transportProperties.lookup("c4_k"));    
dimensionedScalar c5_k (transportProperties.lookup("c5_k"));    
dimensionedScalar c6_k (transportProperties.lookup("c6_k"));    
dimensionedScalar c7_k (transportProperties.lookup("c7_k"));

Info<< "    Reading polynomial coefficients for specific heat capacity" << endl;

dimensionedScalar c0_c (transportProperties.lookup("c0_c"));    
dimensionedScalar c1_c (transportProperties.lookup("c1_c"));    
dimensionedScalar c2_c (transportProperties.lookup("c2_c"));    
dimensionedScalar c3_c (transportProperties.lookup("c3_c"));    
dimensionedScalar c4_c (transportProperties.lookup("c4_c"));    
dimensionedScalar c5_c (transportProperties.lookup("c5_c"));    
dimensionedScalar c6_c (transportProperties.lookup("c6_c"));    
dimensionedScalar c7_c (transportProperties.lookup("c7_c"));        

Info<< "    Reading polynomial coefficients for density" << nl << endl;

dimensionedScalar c0_rho (transportProperties.lookup("c0_rho"));    
dimensionedScalar c1_rho (transportProperties.lookup("c1_rho"));    
dimensionedScalar c2_rho (transportProperties.lookup("c2_rho"));    
dimensionedScalar c3_rho (transportProperties.lookup("c3_rho"));    
dimensionedScalar c4_rho (transportProperties.lookup("c4_rho"));    
dimensionedScalar c5_rho (transportProperties.lookup("c5_rho"));    
dimensionedScalar c6_rho (transportProperties.lookup("c6_rho"));    
dimensionedScalar c7_rho (transportProperties.lookup("c7_rho"));

Info<< "Thermal properties constant" << nl
<< "    Reading thermal diffusivity DT\n" << endl;

dimensionedScalar DT (transportProperties.lookup("DT"));

#include "createFvOptions.H"

write.H - modified

Code:

if (runTime.writeTime())
    {
        volVectorField gradT(fvc::grad(T));

        volScalarField gradTx
        (
            IOobject
            (
                "gradTx",
                runTime.timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::AUTO_WRITE
            ),
            gradT.component(vector::X)
        );

        volScalarField gradTy
        (
            IOobject
            (
                "gradTy",
                runTime.timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::AUTO_WRITE
            ),
            gradT.component(vector::Y)
        );

        volScalarField gradTz
        (
            IOobject
            (
                "gradTz",
                runTime.timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::AUTO_WRITE
            ),
            gradT.component(vector::Z)
        );
    
        Info << nl << endl;
        
        //CALCULATE AND SAVE BOUNDARY HEAT FLUXES:
        forAll(mesh.boundary(), iPatch)
        {
    
        //CELL CENTERS:
    
            //list containing the indices of the boundary neighbour cells (patch cells)
            labelList neighbourCells = mesh.boundary()[iPatch].faceCells();
            
            //scalar field containing the inverse values of the cell center-face distances of the patch cells
            scalarField dCentreCell = mesh.boundary()[iPatch].deltaCoeffs();
    
            //scalar field containing the temperatures of the cell centers
            scalarField cellTemperatures;
            cellTemperatures.resize(T.boundaryField()[iPatch].size());
    
            //scalar field containing the thermal conductivities of the patch cells
            scalarField thermConductivity;
            thermConductivity.resize(T.boundaryField()[iPatch].size());    
    
            //loop over patch cells to obtain cell centre temperatures
            forAll(mesh.boundary()[iPatch], iFace)
            {        
                label iCell = neighbourCells[iFace];                    //index of patch cell
                
                cellTemperatures[iFace] = T.internalField()[iCell];     //temperature of patch cell
                scalar TC = cellTemperatures[iFace];                    //temperature of patch cell
    
                //scalar TC = T.internalField()[iCell];                   //should be identical with the two code lines above - does not work
                
                thermConductivity[iFace] = c0_k.value() * pow(TC ,0) + c1_k.value() * pow(TC ,1) + c2_k.value() * pow(TC ,2) + c3_k.value() * pow(TC ,3) + c4_k.value() * pow(TC ,4) + c5_k.value() * pow(TC ,5) + c6_k.value() * pow(TC ,6) + c7_k.value() * pow(TC ,7);  //th. conductivity of patch cell
            }
    
        //CELL FACES:
            //scalar field containing the areas of the faces on the patch
            scalarField faceAreas = mesh.boundary()[iPatch].magSf();
    
            //scalar field containing the temperatures of the cell faces
            scalarField faceTemperatures = T.boundaryField()[iPatch];
    
        //BOUNDARY HEAT FLUX:
            //scalar field containing the face normals of the face cells
            scalarField patchHeatFlux = thermConductivity * (faceTemperatures - cellTemperatures) * dCentreCell * faceAreas;
    
            //overall heat flux over patch
            scalar Q = gSum(patchHeatFlux);
             Info<< "Total heat flux passing patch " << mesh.boundary()[iPatch].name() << " = " << Q << endl;
    
        //GENERATE IO-OBJECT WITH BOUNDARY HEAT FLUX:
            scalarField boundaryHeatFlux
            (
                IOobject
                (
                    mesh.boundary()[iPatch].name(),
                    runTime.timeName(),
                    T.boundaryField()[iPatch].mesh,        //problem argument - what to put here to make that work?
                    IOobject::NO_READ,
                    IOobject::AUTO_WRITE
                ),
                patchHeatFlux
            );
    
        }
        Info << nl << endl;
        runTime.write();
    }

Generally, I would like to ask for improvements of the code (to decrease calculation time, ...). Every hint is welcome.

Now to ask some specific questions:

1) How do I implement and read a new boolean statement in the controlDict and/or the transportProperties? I tried both and my trials are coloured in red in the code snippets (laplacian.C & createFields.H). Unfortunately none of them work.

2) I would like to create an output file with a list of the heat fluxes of all boundary cells for each patch. Unfortunately this partt of the code doesn't compile, the problem is the registryObject and I don't get what I have to put there to make it work (coloured in blue in write.H). Could somebody briefly explain what I have to put in there?

3) For me the two upper lines to obtain the cell center temperatures coloured in green (write.H) are identical with the single (commented) line below. However, if I use the single line command, the results differ (and are obviously wrong). Does anyone know why that is?

Thank you very much in advance,
stockzahn

snak · February 9, 2019, 10:37

Hi stockzahn,

Quote:

Originally Posted by stockzahn

1) How do I implement and read a new boolean statement in the controlDict and/or the transportProperties? I tried both and my trials are coloured in red in the code snippets (laplacian.C & createFields.H). Unfortunately none of them work.

There will be some confusion of classes and types

If you use a bool type,
- you can use only “true” or “false”, not “yes|no|on|off”,
- there is no need to use an equal-to operator (==) because the bool variable itself will be the answer (true or false).

The function of a dictionary class returns ITstream, not bool.
https://github.com/OpenFOAM/OpenFOAM...tionary.H#L320

It is impossible
- to compare different class and type, and
- to assign a type to another class.

If you prefer to use a bool type, lookupOrDefault function with a bool argument returns bool. (what you give will be returned as a template class)
The following line will work as what you want (I did NOT tested it , sorry…)

Code:

bool varProps = runTime.controlDict().lookupOrDefault(”variable",  false)

When you need the flexible words like yes|no|on|off, a Switch class is needed.
https://github.com/OpenFOAM/OpenFOAM...witch/Switch.H

“runTimeModifiable” entry in controlDict is declared as Swith.
“adjustTimeStep” entry is dealt as Swith by explicitly specifying it.

Code:

lookupOrDefault<Switch>("adjustTimeStep", false)

https://github.com/OpenFOAM/OpenFOAM...me/Time.H#L191

The code of Time class will be a good example of dictionary handling.

stockzahn · February 13, 2019, 10:59

Thanks for your response snak,

unfortunately I don't quite understand your explanation. I tried to get the code in the links you send me, but I wasn't able to correct my mistake.

Quote:

Originally Posted by snak

If you use a bool type,
- you can use only “true” or “false”, not “yes|no|on|off”,
- there is no need to use an equal-to operator (==) because the bool variable itself will be the answer (true or false).

The function of a dictionary class returns ITstream, not bool.

Thanks for these clarifications, but to which equal-to operator (==) you are referring to? In my code

Code:

if (runTime.controlDict().lookup("variable") == "yes" || runTime.controlDict().lookup("variable") == "true")
    {
        varProps = true;
    }

I want to change my bool variable "varProps" to true, if in the dictionary written. As I understood it, the "lookup" function for a dictionary returns the value next to the key word ("variable" in my case). In the condition for the if-statement I just want to check this value, but I do not assign it to a bool variable. It is similar to get the values for thermal diffusivity (or the polynomial coefficients) from the transport properties.

Quote:

Originally Posted by snak

https://github.com/OpenFOAM/OpenFOAM...tionary.H#L320

It is impossible
- to compare different class and type, and
- to assign a type to another class.

So the comparison of the return value of "lookup"-function with a word/string ("yes", "true") seems to be a problem? What kind of types are possible to return with the "lookup"-function?

Quote:

Originally Posted by snak

If you prefer to use a bool type, lookupOrDefault function with a bool argument returns bool. (what you give will be returned as a template class)
The following line will work as what you want (I did NOT tested it , sorry…)

Code:

bool varProps = runTime.controlDict().lookupOrDefault(”variable",   false)

When you need the flexible words like yes|no|on|off, a Switch class is needed.
https://github.com/OpenFOAM/OpenFOAM...witch/Switch.H

“runTimeModifiable” entry in controlDict is declared as Swith.
“adjustTimeStep” entry is dealt as Swith by explicitly specifying it.

Code:

lookupOrDefault<Switch>("adjustTimeStep", false)

https://github.com/OpenFOAM/OpenFOAM...me/Time.H#L191

The code of Time class will be a good example of dictionary handling.

Actually I don't need to return a bool, I just want to make a new entry into the dictionary, like I did for the polynomial coefficients, and dependent on the value/word next to it (preferably a word like "yes" or "true" or "on" instead of numbers 1 or 0), the solver should decide, if variable thermal properties are used or not. If a certain keyword is found, then the bool "varProp" is switched to "true".

Since you pointed out that in the controlDict the option "runTimeModifiable" is defined as Switch: Do all the possible option in the controlDict have to be defined (like in Time.H)? In the transportDict that obviously is not necessary, therefore there should be a possibilty to define keywords. I suppose.

Thanks,
stockzahn

stockzahn · February 13, 2019, 11:03

I've just got a hint from Mowgli, for all of you who calculate polynoms like I do it in the first code snippet of my post #1: Have a look to the Horner scheme! Decreased the calculation time distinctly.

snak · February 13, 2019, 23:42

Hi, stockzahn

Quote:

bool varProps = false;

if (runTime.controlDict().lookup("variable") == "yes" || runTime.controlDict().lookup("variable") == "true")
{
varProps = true;
}

this part can be

Quote:

bool varProps = runTime.controlDict().lookupOrDefault(”variable", false)

If a "variable" entry exists in your controlDict, the value of varProps will be the value written there (true or false).
If there is no entry, the value is false.

lookupOrDefault function can read the entry or can set the default value even if there is no entry.
lookup function will cause an error with message when the entry does not exist.

https://github.com/OpenFOAM/OpenFOAM...emplates.C#L32

---
The other example (not tested)

Code:

Switch varProps = runTime.controlDict().lookupOrDefault<Switch>(”variable", false)

if (varProps)
{ 
    //do something
}

---

Quote:

What kind of types are possible to return with the "lookup"-function?

lookup function returns ITsream& type. it is declared in the following
https://github.com/OpenFOAM/OpenFOAM...tionary.H#L320

---
lookup-func can be used (not tested)

Code:

bool varProps = false;
varProps << runTime.controlDict().lookup("variable")
if (varProps)
{
     // do something
}

lookup-func will cause FatalIOError if the entry varible does not exist.

stockzahn · February 15, 2019, 08:03

Thanks again snak,

that one:

Code:

bool varProps = runTime.controlDict().lookupOrDefault(”variable",   false)

worked perefctly - as I wanted it. I put it in the transportProperties though:

Code:

bool varProps = transportProperties.lookupOrDefault(”variable",   false);

Thanks a lot again!

stockzahn

snak · February 16, 2019, 07:08

Congrats on your success, stockzahn.

bool can only handle yes or no.

Using Switch class will improve the usability with on/off and yes/no.

If you have any chance, could you test the Switch version?

Code:

Switch varProps = transportProperties.lookupOrDefault<Switch>("variable", false);

This will work just same as your code, I guess.

stockzahn · February 18, 2019, 03:41

Dear snak,

Both versions (bool and Switch) work with all keyword pairs (on/off, true/false, yes/no).

Switch:

Code:

Switch varProps =  transportProperties.lookupOrDefault<Switch>("variable",  false);

bool:

Code:

bool varProps =  transportProperties.lookupOrDefault("variable",  false);

stockzahn

February 15, 2019, 08:03		#6
stockzahn Member Join Date: Oct 2018 Location: France Posts: 34 Rep Power: 7	Thanks again snak, that one: Code: bool varProps = runTime.controlDict().lookupOrDefault(”variable", false) worked perefctly - as I wanted it. I put it in the transportProperties though: Code: bool varProps = transportProperties.lookupOrDefault(”variable", false); Thanks a lot again! stockzahn

February 16, 2019, 07:08		#7
snak Senior Member shinji nakagawa Join Date: Mar 2009 Location: Japan Posts: 113 Blog Entries: 1 Rep Power: 18	Congrats on your success, stockzahn. bool can only handle yes or no. Using Switch class will improve the usability with on/off and yes/no. If you have any chance, could you test the Switch version? Code: Switch varProps = transportProperties.lookupOrDefault<Switch>("variable", false); This will work just same as your code, I guess.

February 18, 2019, 03:41		#8
stockzahn Member Join Date: Oct 2018 Location: France Posts: 34 Rep Power: 7	Dear snak, Both versions (bool and Switch) work with all keyword pairs (on/off, true/false, yes/no). Switch: Code: Switch varProps = transportProperties.lookupOrDefault<Switch>("variable", false); bool: Code: bool varProps = transportProperties.lookupOrDefault("variable", false); stockzahn Last edited by stockzahn; February 18, 2019 at 08:23.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
New Solver for a modified pUCoupledFoam	noiseMaker	OpenFOAM	0	August 2, 2016 06:14
Error with modified interFoam solver	robertmaier9	OpenFOAM	3	September 30, 2011 06:26
Working directory via command line	Luiz	CFX	4	March 6, 2011 20:02
why the solver reject it? Anyone with experience?	bearcat	CFX	6	April 28, 2008 14:08
Laplacian viscous stress term in compressible solver	jelmer	OpenFOAM Running, Solving & CFD	3	June 23, 2006 07:31

February 13, 2019, 11:03		#4
stockzahn Member Join Date: Oct 2018 Location: France Posts: 34 Rep Power: 7	I've just got a hint from Mowgli, for all of you who calculate polynoms like I do it in the first code snippet of my post #1: Have a look to the Horner scheme! Decreased the calculation time distinctly.