You are not logged in.
For SpaceNut .... The C language deserves it's own topic. It is in a class by itself.
Post #3 will show a C language program intended to compute thrust produced by an OpenFOAM model.
It was developed by ChatGPT5 (with a bit of effort), and I don't know if it is right or even close to right.
This topic is available for NewMars members to post C language programs or C++ ones.
(th)
Offline
Like button can go here
This post is reserved for an index to posts that may be contributed by NewMars members.
Index:
Post #3: C program for OpenFOAM to compute thrust
(th)
Offline
Like button can go here
This post is for thrustCalc
thrustCalc was written by ChatGPT5 (with a bit of efford) on Monday 2025/11/03
I have NO idea if the program is right or even close to right.
// thrustCalc.C — OpenFOAM v12-compatible, standalone thrust post-processor
// Computes force on the FLUID across a patch:
// T = ∫_patch [ rho*(U·nHat)*U dA + (p - pAmb)*Sf ]
// Vehicle thrust = negative of this vector.
#include "argList.H"
#include "timeSelector.H"
#include "Time.H"
#include "fvMesh.H"
#include "volFields.H"
#include "surfaceFields.H"
#include "fvc.H"
#include "IOdictionary.H"
#include "polyBoundaryMesh.H"
#include "word.H"
#include "vector.H"
#include "OFstream.H"
#include "OSspecific.H" // mkDir
#include <fstream>
#include <iomanip>
using namespace Foam;
int main(int argc, char *argv[])
{
// Enable time-selection flags before args are created by setRootCase.H
timeSelector::addOptions();
// Nice help text with --help
argList::addNote
(
"Compute thrust on a patch:\n"
" T = ∫_patch [ rho*(U·nHat)*U dA + (p - pAmb)*Sf ]\n"
"Force reported is on the FLUID; vehicle thrust is its negative."
);
#include "setRootCase.H"
#include "createTime.H"
instantList timeDirs = timeSelector::select0(runTime, args);
#include "createMesh.H"
// Read settings from system/thrustCalcDict
IOdictionary dict
(
IOobject
(
"thrustCalcDict",
runTime.system(),
runTime,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);
word patchName("outlet");
dict.readIfPresent("patch", patchName);
word inletName(""); // leave absent to ignore
dict.readIfPresent("inletPatch", inletName);
scalar pAmb(0.0);
dict.readIfPresent("pAmb", pAmb);
// Output path
fileName outDir = runTime.path()/"postProcessing"/"thrustCalc"/patchName;
mkDir(outDir);
const fileName outFile = outDir/"thrust.dat";
// Helper: find patch index by name (v12-safe)
auto findPatch = [&](const word& name) -> label
{
if (!name.size()) return -1;
const fvBoundaryMesh& bnd = mesh.boundary();
forAll(bnd, i)
{
if (bnd[i].name() == name) return i;
}
return -1;
};
forAll(timeDirs, ti)
{
runTime.setTime(timeDirs[ti], ti);
mesh.readUpdate();
// Exact time directory name (no rounding)
const word curTName = timeDirs[ti].name();
// Read fields at this time
volVectorField U
(
IOobject("U", curTName, mesh, IOobject::MUST_READ, IOobject::NO_WRITE),
mesh
);
volScalarField p
(
IOobject("p", curTName, mesh, IOobject::MUST_READ, IOobject::NO_WRITE),
mesh
);
volScalarField rho
(
IOobject("rho", curTName, mesh, IOobject::MUST_READ, IOobject::NO_WRITE),
mesh
);
// Compute contribution on a patch (fields are in scope here)
auto computePatch = [&](const word& name, bool subtract) -> vector
{
if (!name.size()) return vector::zero; // silent ignore
vector result(vector::zero);
const label pid = findPatch(name);
if (pid < 0)
{
Info<< "thrustCalc: patch '"<< name <<"' not found (skipping)." << nl;
return result;
}
const vectorField& Sf = mesh.Sf().boundaryField()[pid];
const fvPatchVectorField& Ub = U.boundaryField()[pid];
const fvPatchScalarField& pb = p.boundaryField()[pid];
const fvPatchScalarField& rb = rho.boundaryField()[pid];
const label N = Sf.size();
for (label i=0; i<N; ++i)
{
const scalar magSf = mag(Sf[i]);
const vector nHat = (magSf > VSMALL ? Sf[i]/magSf : vector::zero);
const scalar Un = (Ub[i] & nHat);
const vector mom = rb[i]*Un*Ub[i]*magSf; // rho * Un * U * dA
const vector pTerm = (pb[i] - pAmb)*Sf[i]; // (p - pAmb) * Sf
result += (mom + pTerm);
}
return subtract ? -result : result;
};
const vector thrust = computePatch(patchName, /*subtract*/false)
+ computePatch(inletName, /*subtract*/true);
const scalar t = runTime.value();
const scalar Fmag = mag(thrust);
// Decide if we need a header (file does not exist or empty)
bool needHeader = true;
{
std::ifstream chk(outFile.c_str());
needHeader = !chk.good() || (chk.peek() == std::ifstream::traits_type::eof());
}
// Append one row (with header if needed)
std::ofstream os(outFile.c_str(), std::ios_base::app);
if (needHeader)
{
os << "# timeDir t[s] Fx[N] Fy[N] Fz[N] |F|[N]\n";
}
os.setf(std::ios_base::scientific);
os << curTName << " "
<< std::setprecision(10) << t << " "
<< thrust.x() << " " << thrust.y() << " " << thrust.z() << " "
<< Fmag << "\n";
os.close();
Info<< "thrustCalc: t="<< t
<< " Thrust_on_fluid [N] = " << thrust
<< " |F|=" << Fmag
<< " (vehicle reaction = -value above)" << nl;
}
Info<< "thrustCalc: wrote " << outFile << nl;
return 0;
}See the OpenFOAM topic for output of the program and related discussion.
When I say "with some effort" I am commenting upon the struggle that I observed as ChatGPT5 was attempting to create this program. As I watched the process unfold, I noted that OpenFOAM has multiple versions and the resources ChatGPT5 consults include human generated data which is (apparently) often misleading or wrong, or out of date, or otherwise difficult to work with.
On top of that, the problem to be solved is (from my perspective) extremely complex. The status of many thousands of cells is preserved in each time record. Time records are recorded every .05 seconds of simulation time, so we have five time directories for each 1/4 second of run time.
(th)
Offline
Like button can go here