songtianlun/openfoam
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

ENH: additional in-place clamp_min(), clamp_max() field methods

Browse Source 
- these already existed for a single value, but now handle the full
  field. This is more memory-friendly.

      fld.clamp_min(lower);  OLD: fld = max(fld, lower);
      fld.clamp_max(upper);  OLD: fld = min(fld, upper);
This commit is contained in:
Mark Olesen 2025-04-02 11:19:04 +02:00
parent b8a0706e72
commit 8efef93a8c
18 changed files with 258 additions and 103 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

46
applications/solvers/combustion/reactingFoam/setRDeltaT.H
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2013-2016 OpenFOAM Foundation
Copyright (C) 2020 OpenCFD Ltd.
Copyright (C) 2020,2025 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -57,10 +57,22 @@ License
// (relative to reference value)
scalar alphaY(pimpleDict.getOrDefault<scalar>("alphaY", 1.0));
Info<< "Time scales min/max:" << endl;
// Cache old reciprocal time scale field
volScalarField rDeltaT0("rDeltaT0", rDeltaT);
// The old reciprocal time scale field, with any damping factor
tmp<volScalarField> rDeltaT0_damped;
// Calculate damped value before applying any other changes
if
(
rDeltaTDampingCoeff < 1
&& runTime.timeIndex() > runTime.startTimeIndex() + 1
)
{
rDeltaT0_damped = (scalar(1) - rDeltaTDampingCoeff)*(rDeltaT);
}
Info<< "Time scales min/max:" << endl;
// Flow time scale
{
@ -70,8 +82,8 @@ License
/((2*maxCo)*mesh.V()*rho())
);
// Limit the largest time scale
rDeltaT.max(1/maxDeltaT);
// Limit the largest time scale (=> smallest reciprocal time)
rDeltaT.clamp_min(1/maxDeltaT);
Info<< " Flow = "
<< 1/gMax(rDeltaT.primitiveField()) << ", "
@ -86,11 +98,11 @@ License
mag(Qdot)/(alphaTemp*rho*thermo.Cp()*T)
);
rDeltaT.primitiveFieldRef().clamp_min(rDeltaTT);
Info<< " Temperature = "
<< 1/(gMax(rDeltaTT.field()) + VSMALL) << ", "
<< 1/(gMin(rDeltaTT.field()) + VSMALL) << endl;
rDeltaT.ref() = max(rDeltaT(), rDeltaTT);
}
// Reaction rate time scale
@ -138,11 +150,11 @@ License
if (foundY)
{
rDeltaT.primitiveFieldRef().clamp_min(rDeltaTY);
Info<< " Composition = "
<< 1/(gMax(rDeltaTY.field()) + VSMALL) << ", "
<< 1/(gMin(rDeltaTY.field()) + VSMALL) << endl;
rDeltaT.ref() = max(rDeltaT(), rDeltaTY);
}
else
{
@ -161,20 +173,12 @@ License
fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
}
// Limit rate of change of time scale
// Limit rate of change of time scale (=> smallest reciprocal time)
// - reduce as much as required
// - only increase at a fraction of old time scale
if
(
rDeltaTDampingCoeff < 1
&& runTime.timeIndex() > runTime.startTimeIndex() + 1
)
if (rDeltaT0_damped)
{
rDeltaT = max
(
rDeltaT,
(scalar(1) - rDeltaTDampingCoeff)*rDeltaT0
);
rDeltaT.clamp_min(rDeltaT0_damped());
}
// Update tho boundary values of the reciprocal time-step

4
applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
View File

@ -78,8 +78,8 @@
}
rho = thermo.rho();
rho = max(rho, rhoMin[i]);
rho = min(rho, rhoMax[i]);
rho.clamp_range(rhoMin[i], rhoMax[i]);
rho.relax();
Info<< "Min/max rho:" << min(rho).value() << ' '

45
applications/solvers/lagrangian/coalChemistryFoam/setRDeltaT.H
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2020 OpenCFD Ltd.
Copyright (C) 2020,2025 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -54,10 +54,21 @@ License
scalar alphaTemp(pimpleDict.getOrDefault("alphaTemp", 0.05));
Info<< "Time scales min/max:" << endl;
// The old reciprocal time scale field, with any damping factor
tmp<volScalarField> rDeltaT0_damped;
// Cache old reciprocal time scale field
volScalarField rDeltaT0("rDeltaT0", rDeltaT);
// Calculate damped value before applying any other changes
if
(
rDeltaTDampingCoeff < 1
&& runTime.timeIndex() > runTime.startTimeIndex() + 1
)
{
rDeltaT0_damped = (scalar(1) - rDeltaTDampingCoeff)*(rDeltaT);
}
Info<< "Time scales min/max:" << endl;
// Flow time scale
{
@ -67,8 +78,8 @@ License
/((2*maxCo)*mesh.V()*rho())
);
// Limit the largest time scale
rDeltaT.max(1/maxDeltaT);
// Limit the largest time scale (=> smallest reciprocal time)
rDeltaT.clamp_min(1/maxDeltaT);
Info<< " Flow = "
<< gMin(1/rDeltaT.primitiveField()) << ", "
@ -93,15 +104,11 @@ License
)
);
rDeltaT.primitiveFieldRef().clamp_min(rDeltaTT);
Info<< " Temperature = "
<< gMin(1/(rDeltaTT.field() + VSMALL)) << ", "
<< gMax(1/(rDeltaTT.field() + VSMALL)) << endl;
rDeltaT.ref() = max
(
rDeltaT(),
rDeltaTT
);
}
// Update tho boundary values of the reciprocal time-step
@ -113,20 +120,12 @@ License
fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
}
// Limit rate of change of time scale
// Limit rate of change of time scale (=> smallest reciprocal time)
// - reduce as much as required
// - only increase at a fraction of old time scale
if
(
rDeltaTDampingCoeff < 1.0
&& runTime.timeIndex() > runTime.startTimeIndex() + 1
)
if (rDeltaT0_damped)
{
rDeltaT = max
(
rDeltaT,
(scalar(1) - rDeltaTDampingCoeff)*rDeltaT0
);
rDeltaT.clamp_min(rDeltaT0_damped());
}
Info<< " Overall = "

48
applications/solvers/lagrangian/reactingParcelFoam/setRDeltaT.H
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2020 OpenCFD Ltd.
Copyright (C) 2020,2025 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -54,10 +54,21 @@ License
scalar alphaTemp(pimpleDict.getOrDefault("alphaTemp", 0.05));
Info<< "Time scales min/max:" << endl;
// The old reciprocal time scale field, with any damping factor
tmp<volScalarField> rDeltaT0_damped;
// Cache old reciprocal time scale field
volScalarField rDeltaT0("rDeltaT0", rDeltaT);
// Calculate damped value before applying any other changes
if
(
rDeltaTDampingCoeff < 1
&& runTime.timeIndex() > runTime.startTimeIndex() + 1
)
{
rDeltaT0_damped = (scalar(1) - rDeltaTDampingCoeff)*(rDeltaT);
}
Info<< "Time scales min/max:" << endl;
// Flow time scale
{
@ -67,8 +78,8 @@ License
/((2*maxCo)*mesh.V()*rho())
);
// Limit the largest time scale
rDeltaT.max(1/maxDeltaT);
// Limit the largest time scale (=> smallest reciprocal time)
rDeltaT.clamp_min(1/maxDeltaT);
Info<< " Flow = "
<< gMin(1/rDeltaT.primitiveField()) << ", "
@ -92,15 +103,11 @@ License
)
);
rDeltaT.primitiveFieldRef().clamp_min(rDeltaTT);
Info<< " Temperature = "
<< gMin(1/(rDeltaTT.field() + VSMALL)) << ", "
<< gMax(1/(rDeltaTT.field() + VSMALL)) << endl;
rDeltaT.ref() = max
(
rDeltaT(),
rDeltaTT
);
}
// Update the boundary values of the reciprocal time-step
@ -112,22 +119,17 @@ License
fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
}
// Limit rate of change of time scale
// Limit rate of change of time scale (=> smallest reciprocal time)
// - reduce as much as required
// - only increase at a fraction of old time scale
if
(
rDeltaTDampingCoeff < 1.0
&& runTime.timeIndex() > runTime.startTimeIndex() + 1
)
if (rDeltaT0_damped)
{
rDeltaT = max
(
rDeltaT,
(scalar(1) - rDeltaTDampingCoeff)*rDeltaT0
);
rDeltaT.clamp_min(rDeltaT0_damped());
}
// Update the boundary values of the reciprocal time-step
rDeltaT.correctBoundaryConditions();
Info<< " Overall = "
<< gMin(1/rDeltaT.primitiveField())
<< ", " << gMax(1/rDeltaT.primitiveField()) << endl;

3
applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H
View File

@ -48,8 +48,7 @@ U.correctBoundaryConditions();
fvOptions.correct(U);
rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);
rho.clamp_range(rhoMin, rhoMax);
rho.relax();
Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;

3
applications/solvers/lagrangian/simpleCoalParcelFoam/pEqn.H
View File

@ -49,8 +49,7 @@
fvOptions.correct(U);
rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);
rho.clamp_range(rhoMin, rhoMax);
rho.relax();
Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;

6
applications/solvers/lagrangian/sprayFoam/pEqn.H
View File

@ -1,6 +1,5 @@
rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);
rho.clamp_range(rhoMin, rhoMax);
rho.relax();
volScalarField rAU(1.0/UEqn.A());
@ -94,8 +93,7 @@ p.relax();
// Recalculate density from the relaxed pressure
rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);
rho.clamp_range(rhoMin, rhoMax);
rho.relax();
Info<< "rho min/max : " << min(rho).value() << " " << max(rho).value() << endl;

6
applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/pEqn.H
View File

@ -1,6 +1,5 @@
rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);
rho.clamp_range(rhoMin, rhoMax);
rho.relax();
volScalarField rAU(1.0/UEqn.A());
@ -94,8 +93,7 @@ p.relax();
// Recalculate density from the relaxed pressure
rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);
rho.clamp_range(rhoMin, rhoMax);
rho.relax();
Info<< "rho min/max : " << min(rho).value() << " " << max(rho).value() << endl;

29
applications/solvers/multiphase/VoF/setRDeltaT.H
View File

@ -53,6 +53,21 @@
pimpleDict.getOrDefault<scalar>("maxDeltaT", GREAT)
);
// The old reciprocal time scale field, with any damping factor
tmp<volScalarField> rDeltaT0_damped;
// Calculate damped value before applying any other changes
if
(
rDeltaTDampingCoeff < 1
&& runTime.timeIndex() > runTime.startTimeIndex() + 1
)
{
rDeltaT0_damped = (scalar(1) - rDeltaTDampingCoeff)*(rDeltaT);
}
volScalarField rDeltaT0("rDeltaT0", rDeltaT);
// Set the reciprocal time-step from the local Courant number
@ -114,20 +129,12 @@
<< gMin(1/rDeltaT.primitiveField())
<< ", " << gMax(1/rDeltaT.primitiveField()) << endl;
// Limit rate of change of time scale
// Limit rate of change of time scale (=> smallest reciprocal time)
// - reduce as much as required
// - only increase at a fraction of old time scale
if
(
rDeltaTDampingCoeff < 1.0
&& runTime.timeIndex() > runTime.startTimeIndex() + 1
)
if (rDeltaT0_damped)
{
rDeltaT = max
(
rDeltaT,
(scalar(1) - rDeltaTDampingCoeff)*rDeltaT0
);
rDeltaT.clamp_min(rDeltaT0_damped());
Info<< "Damped flow time scale min/max = "
<< gMin(1/rDeltaT.primitiveField())

32
applications/test/minMax2/Test-minMax2.cxx
View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2019-2023 OpenCFD Ltd.
Copyright (C) 2019-2025 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -175,6 +175,36 @@ int main(int argc, char *argv[])
<< nl;
}
{
scalarField field1(10);
scalarField field2(10);
scalarField work;
Random rnd(4567);
for (scalar& val : field1)
{
val = rnd.position(scalar(-0.2), scalar(1.2));
}
for (scalar& val : field2)
{
val = rnd.position(scalar(-0.1), scalar(1.1));
}
Info<< nl
<< "field1: " << flatOutput(field1) << nl
<< "field2: " << flatOutput(field2) << nl;
work = field1;
work.clamp_min(field2);
Info<< "clamp_min: " << flatOutput(work) << nl;
work = field1;
work.clamp_max(field2);
Info<< "clamp_max: " << flatOutput(work) << nl;
}
Info<< nl << "\nDone\n" << endl;
return 0;
}

32
src/OpenFOAM/fields/FieldFields/FieldField/FieldField.C
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2019-2023 OpenCFD Ltd.
Copyright (C) 2019-2025 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -308,6 +308,36 @@ void FieldField<Field, Type>::clamp_max
}
template<template<class> class Field, class Type>
void FieldField<Field, Type>::clamp_min
(
const FieldField<Field, Type>& lower
)
{
const label loopLen = this->size();
for (label i = 0; i < loopLen; ++i)
{
(*this)[i].clamp_min(lower[i]);
}
}
template<template<class> class Field, class Type>
void FieldField<Field, Type>::clamp_max
(
const FieldField<Field, Type>& upper
)
{
const label loopLen = this->size();
for (label i = 0; i < loopLen; ++i)
{
(*this)[i].clamp_max(upper[i]);
}
}
template<template<class> class Field, class Type>
void FieldField<Field, Type>::clamp_range
(

8
src/OpenFOAM/fields/FieldFields/FieldField/FieldField.H
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2022-2023 OpenCFD Ltd.
Copyright (C) 2022-2025 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -152,9 +152,15 @@ public:
//- Impose lower (floor) clamp on the field values (in-place)
void clamp_min(const Type& lower);
//- Impose lower (floor) clamp on the field values (in-place)
void clamp_min(const FieldField<Field, Type>& lower);
//- Impose upper (ceiling) clamp on the field values (in-place)
void clamp_max(const Type& upper);
//- Impose upper (ceiling) clamp on the field values (in-place)
void clamp_max(const FieldField<Field, Type>& upper);
//- Clamp field values (in-place) to the specified range.
// Does not check if range is valid or not.
void clamp_range(const Type& lower, const Type& upper);

35
src/OpenFOAM/fields/Fields/Field/Field.C
View File

@ -669,6 +669,7 @@ void Foam::Field<Type>::clamp_min(const Type& lower)
}
}
template<class Type>
void Foam::Field<Type>::clamp_max(const Type& upper)
{
@ -681,6 +682,40 @@ void Foam::Field<Type>::clamp_max(const Type& upper)
}
template<class Type>
void Foam::Field<Type>::clamp_min(const UList<Type>& lower)
{
// Use free function max() [sic] to impose component-wise clamp_min
std::transform
(
// Can use (std::execution::par_unseq | std::execution::unseq)
this->begin(),
// this->end() but with some extra range safety
this->begin(lower.size()),
lower.begin(),
this->begin(),
maxOp<Type>()
);
}
template<class Type>
void Foam::Field<Type>::clamp_max(const UList<Type>& upper)
{
// Use free function min() [sic] to impose component-wise clamp_max
std::transform
(
// Can use (std::execution::par_unseq | std::execution::unseq)
this->begin(),
// this->end() but with some extra range safety
this->begin(upper.size()),
upper.begin(),
this->begin(),
minOp<Type>()
);
}
template<class Type>
void Foam::Field<Type>::clamp_range(const Type& lower, const Type& upper)
{

9
src/OpenFOAM/fields/Fields/Field/Field.H
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2015-2023 OpenCFD Ltd.
Copyright (C) 2015-2025 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -50,6 +50,7 @@ SourceFiles
#include "direction.H"
#include "labelList.H"
#include "keyType.H"
#include "ops.H"
#include "scalarList.H"
#include "VectorSpace.H"
#include "IOobjectOption.H"
@ -441,9 +442,15 @@ public:
//- Impose lower (floor) clamp on the field values (in-place)
void clamp_min(const Type& lower);
//- Impose lower (floor) clamp on the field values (in-place)
void clamp_min(const UList<Type>& lower);
//- Impose upper (ceiling) clamp on the field values (in-place)
void clamp_max(const Type& upper);
//- Impose upper (ceiling) clamp on the field values (in-place)
void clamp_max(const UList<Type>& upper);
//- Clamp field values (in-place) to the specified range.
// Does not check if range is valid or not.
void clamp_range(const Type& lower, const Type& upper);

32
src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
View File

@ -1276,6 +1276,38 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::clamp_max
}
template<class Type, template<class> class PatchField, class GeoMesh>
void Foam::GeometricField<Type, PatchField, GeoMesh>::clamp_min
(
const GeometricField<Type, PatchField, GeoMesh>& lower
)
{
primitiveFieldRef().clamp_min(lower.primitiveField());
boundaryFieldRef().clamp_min(lower.boundaryField());
correctLocalBoundaryConditions();
if (GeometricBoundaryField<Type, PatchField, GeoMesh>::debug)
{
boundaryField().check();
}
}
template<class Type, template<class> class PatchField, class GeoMesh>
void Foam::GeometricField<Type, PatchField, GeoMesh>::clamp_max
(
const GeometricField<Type, PatchField, GeoMesh>& upper
)
{
primitiveFieldRef().clamp_max(upper.primitiveField());
boundaryFieldRef().clamp_max(upper.boundaryField());
correctLocalBoundaryConditions();
if (GeometricBoundaryField<Type, PatchField, GeoMesh>::debug)
{
boundaryField().check();
}
}
template<class Type, template<class> class PatchField, class GeoMesh>
void Foam::GeometricField<Type, PatchField, GeoMesh>::clamp_range
(

16
src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H
View File

@ -936,17 +936,25 @@ public:
//- Impose lower (floor) clamp on the field values (in-place)
void clamp_min(const Type& lower);
//- Impose upper (ceiling) clamp on the field values (in-place)
void clamp_max(const Type& upper);
//- Impose lower (floor) clamp on the field values (in-place)
// No dimension checking
void clamp_min(const dimensioned<Type>& lower);
//- Impose lower (floor) clamp on the field values (in-place)
// No dimension checking
void clamp_min(const GeometricField<Type, PatchField, GeoMesh>& lower);
//- Impose upper (ceiling) clamp on the field values (in-place)
void clamp_max(const Type& upper);
//- Impose upper (ceiling) clamp on the field values (in-place)
// No dimension checking
void clamp_max(const dimensioned<Type>& upper);
//- Impose upper (ceiling) clamp on the field values (in-place)
// No dimension checking
void clamp_max(const GeometricField<Type, PatchField, GeoMesh>& upper);
//- Clamp field values (in-place) to the specified range.
// Does not check if range is valid or not. No dimension checking.
void clamp_range(const dimensioned<MinMax<Type>>& range);
@ -1055,11 +1063,13 @@ public:
//- Use minimum of the field and specified value. ie, clamp_max().
// This sets the \em ceiling on the field values
// \deprecated(2023-01) prefer clamp_max()
FOAM_DEPRECATED_STRICT(2023-01, "clamp_max() method")
void min(const dimensioned<Type>& upper) { this->clamp_max(upper); }
//- Use maximum of the field and specified value. ie, clamp_min().
// This sets the \em floor on the field values
// \deprecated(2023-01) prefer clamp_min()
FOAM_DEPRECATED_STRICT(2023-01, "clamp_min() method")
void max(const dimensioned<Type>& lower) { this->clamp_min(lower); }
//- Deprecated(2019-01) identical to clamp_range()

4
src/finiteArea/faMatrices/faScalarMatrix/faScalarMatrix.C
View File

@ -85,7 +85,7 @@ Foam::solverPerformance Foam::faMatrix<Foam::scalar>::solve
internalCoeffs_,
psi_.boundaryField().scalarInterfaces(),
solverControls
)->solve(psi.ref(), totalSource);
)->solve(psi.primitiveFieldRef(), totalSource);
if (logLevel)
{
@ -146,7 +146,7 @@ Foam::tmp<Foam::areaScalarField> Foam::faMatrix<Foam::scalar>::H() const
Hphi.primitiveFieldRef() = (lduMatrix::H(psi_.primitiveField()) + source_);
addBoundarySource(Hphi.primitiveFieldRef());
Hphi.ref() /= psi_.mesh().S();
Hphi.primitiveFieldRef() /= psi_.mesh().S();
Hphi.correctBoundaryConditions();
return tHphi;

3
src/thermophysicalModels/basic/rhoThermo/rhoThermo.C
View File

@ -256,8 +256,7 @@ void Foam::rhoThermo::correctRho
)
{
rho_ += deltaRho;
rho_ = max(rho_, rhoMin);
rho_ = min(rho_, rhoMax);
rho_.clamp_range(rhoMin, rhoMax);
}
void Foam::rhoThermo::correctRho(const Foam::volScalarField& deltaRho)