CFD Online Logo CFD Online URL
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Running, Solving & CFD

Read-in variable source term from files

Register Blogs Members List Search Today's Posts Mark Forums Read

LinkBack Thread Tools Search this Thread Display Modes
Old   May 25, 2016, 10:48
Default Read-in variable source term from files
New Member
Join Date: Jul 2015
Posts: 5
Rep Power: 7
FlorisvdBeek is on a distinguished road
Hi folks,

For a project, Im trying to simulate the effects of a body force created by a plasma actuator on the flow. The solver is based on rhoPimpleFoam (OpenFOAM 3.0.1), with an additional term in the momentum equation, added to the UEqn.H file.

I have the body force fields as external files (from measurements), so I would like to read in the force fields from these files. As exactly this is happening at the initialization (createFields.H) in folder 0/, I created a readForceField.H:

// Read in the body force

Info<< "Reading field bodyForce\n" << endl; // bodyForce term
volVectorField bodyForce
Which is executed right before entering the PIMPLE loop:

  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
     \\/     M anipulation  |
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <>.


    Transient solver for laminar or turbulent flow of compressible fluids
    for HVAC and similar applications.

    Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
    pseudo-transient simulations.


#include "fvCFD.H"
#include "psiThermo.H"
#include "turbulentFluidThermoModel.H"
#include "bound.H"
#include "pimpleControl.H"
#include "fvIOoptionList.H"
#include "localEulerDdtScheme.H"
#include "fvcSmooth.H"

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

int main(int argc, char *argv[])
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"

    pimpleControl pimple(mesh);

    // initialization	
    scalar currentTime=0;
//     scalar currentForcePhase=0;
//     scalar currentEnergyPhase=0;
    scalar period=0;
    scalar globalBodyForceVal=0;
    scalar globalEnergyTermVal=0;
    #include "createTimeControls.H"
    #include "createRDeltaT.H"
    #include "initContinuityErrs.H"
    #include "createFields.H"
    #include "createMRF.H"
    #include "createFvOptions.H"
    // two extra header files that are needed  mk

    #include "createBodyForce4.H"
    scalar auxTime=auxTime0;

    #include "mathematicalConstants.H"

    if (!LTS)
        #include "compressibleCourantNo.H"
        #include "setInitialDeltaT.H"

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


    Info<< "\nStarting time loop\n" << endl;

    while (
        #include "readTimeControls.H"

        if (LTS)
            #include "setRDeltaT.H"
            #include "compressibleCourantNo.H"
            #include "setDeltaT.H"


        Info<< "Time = " << runTime.timeName() << nl << endl;
	if (  currentTime == auxTime ) 
	#include "readForceField.H"    
	  auxTime +=auxTime0; 

        if (pimple.nCorrPIMPLE() <= 1)

	  #include "rhoEqn.H"

        // --- Pressure-velocity PIMPLE corrector loop
        while (pimple.loop())
            #include "UEqn.H"
            #include "EEqn.H"

            // --- Pressure corrector loop
            while (pimple.correct())
                if (pimple.consistent())
                    #include "pcEqn.H"
                    #include "pEqn.H"

            if (pimple.turbCorr())


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

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

    return 0;

// ************************************************************************* //
The files are located in folders which would be otherwise created (i.e. 0/, 0.001/, 0.002/ etc..). The read-in is controlled by the if-statement, which checks whether the currentTime is equal to a multitude of auxTime (which initially has the value of writeInterval, the time step after which the solver should write its output, and is increased once currentTime has exceded it). This is necessary as the timestep between the body force measurements is much larger than the timestep of the solver, so the solver should read in the files only once every "writeInterval" timestep. Now my problems:
  • The first is the if-statement: "if ( currentTime == auxTime ) ", which does not seem to work.. Both values are scalars, but the solver is skipping the loop even if I monitor that currentTime is equal to auxTime.
  • My second problem (if I remove the if-statement), the solver does not change the values of the body force after 'reading in' the file at the new timestep: nothing is happening.

Any ideas/clues/alternative ways to do what Id like to do?

Thanks in advance !
FlorisvdBeek is offline   Reply With Quote

Old   October 25, 2016, 16:10
Default Floating point equality comparison is harmful and often wrong
Join Date: Oct 2013
Posts: 92
Rep Power: 9
fedvasu is on a distinguished road
Hi Florisvdbeek,

I am pretty sure, you have solved this issue by now...

but just for documentation.

if ( currentTime == auxTime )
This won't and shouldn't work, because comparing floating points is both pointless and perilious.

either convert the times into strings and make sure that each time that directory has the correct time(string value) or if you really want to compare times, make times as label (int).

I am not sure about the second part.
fedvasu is offline   Reply With Quote


Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On

Similar Threads
Thread Thread Starter Forum Replies Last Post
[swak4Foam] funkyDoCalc with OF2.3 massflow NiFl OpenFOAM Community Contributions 14 November 25, 2020 04:30
[] Patches to compile OpenFOAM 2.2 on Mac OS X gschaider OpenFOAM Installation 136 October 10, 2017 18:25
[Other] OpenFOAM Installation for navalFoam sachinlb OpenFOAM Community Contributions 22 July 28, 2017 06:26
Read txt file and import values to source term. Constantine Fluent UDF and Scheme Programming 40 July 21, 2016 18:20

All times are GMT -4. The time now is 20:36.