Merge remote-tracking branch 'origin/master' into develop

This commit is contained in:
Mark Olesen 2024-12-16 10:07:54 +01:00
commit 4c21ad3d3f
13 changed files with 249 additions and 159 deletions

View File

@ -119,10 +119,10 @@ public:
#if __cplusplus >= 201703L
//- Construct (shallow copy) from std::string_view content
explicit ispanstream(std::string_view s)
{
:
buffer_type(const_cast<char*>(s.data()), s.size()),
stream_type(static_cast<buffer_type*>(this));
}
stream_type(static_cast<buffer_type*>(this))
{}
#endif
//- Construct (shallow copy) from span character content

View File

@ -157,7 +157,10 @@ Foam::error::error(const error& err)
throwing_(err.throwing_),
messageStreamPtr_(nullptr)
{
if (err.messageStreamPtr_ && (err.messageStreamPtr_->count() > 0))
// FIXME: OStringStream copy construct does not adjust tellp and
// thus the count is wrong! (#3281)
// // if (err.messageStreamPtr_ && (err.messageStreamPtr_->count() > 0))
if (err.messageStreamPtr_)
{
messageStreamPtr_.reset(new OStringStream(*err.messageStreamPtr_));
}

View File

@ -888,7 +888,7 @@ boundaryInternalField() const
*this
);
auto& result = tresult;
auto& result = tresult.ref();
forAll(result, patchi)
{

View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2017-2023 OpenCFD Ltd.
Copyright (C) 2017-2024 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -136,7 +136,8 @@ std::string Foam::fileFormats::NASCore::nextNasField
(
const std::string& str,
std::string::size_type& pos,
std::string::size_type len
const std::string::size_type width,
const bool free_format
)
{
const auto beg = pos;
@ -144,15 +145,23 @@ std::string Foam::fileFormats::NASCore::nextNasField
if (end == std::string::npos)
{
pos = beg + len; // Continue after field width
if (free_format)
{
// Nothing left
pos = str.size();
return str.substr(beg);
}
// Fixed format - continue after field width
pos = beg + width;
return str.substr(beg, width);
}
else
{
len = (end - beg); // Efffective width
pos = end + 1; // Continue after comma
// Free format - continue after comma
pos = end + 1;
return str.substr(beg, (end - beg));
}
return str.substr(beg, len);
}
@ -249,8 +258,8 @@ void Foam::fileFormats::NASCore::writeCoord
// 2 ID : point ID - requires starting index of 1
// 3 CP : coordinate system ID (blank)
// 4 X1 : point x coordinate
// 5 X2 : point x coordinate
// 6 X3 : point x coordinate
// 5 X2 : point y coordinate
// 6 X3 : point z coordinate
// 7 CD : coordinate system for displacements (blank)
// 8 PS : single point constraints (blank)
// 9 SEID : super-element ID

View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2017-2023 OpenCFD Ltd.
Copyright (C) 2017-2024 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -47,7 +47,6 @@ SourceFiles
namespace Foam
{
namespace fileFormats
{
@ -79,18 +78,18 @@ public:
//- Output load format
enum loadFormat
{
PLOAD2,
PLOAD4
PLOAD2, //!< Face load (eg, pressure)
PLOAD4 //!< Vertex load
};
//- Selection names for the NASTRAN file field formats
//- Selection names for the NASTRAN load formats
static const Enum<loadFormat> loadFormatNames;
// Constructors
//- Default construct
NASCore() = default;
NASCore() noexcept = default;
// Public Static Member Functions
@ -98,18 +97,20 @@ public:
//- Extract numbers from things like "-2.358-8" (same as "-2.358e-8")
static scalar readNasScalar(const std::string& str);
//- A string::substr() to handle fixed-format and free-format NASTRAN.
// Returns the substr to the next comma (if found) or the given length
//
// \param str The string to extract from
// \param pos On input, the position of the first character of the
// substring. On output, advances to the next position to use.
// \param len The fixed-format length to use if a comma is not found.
//- A std::string::substr() variant to handle fixed-format and
//- free-format NASTRAN.
// Returns the substr until the next comma (if found)
// or the given fixed width
static std::string nextNasField
(
//! The string to extract from
const std::string& str,
//! [in,out] The parse position within \p str
std::string::size_type& pos,
std::string::size_type len
//! The fixed-format width to use (if comma is not found)
const std::string::size_type width,
//! The input is known to be free-format
const bool free_format = false
);

View File

@ -62,13 +62,17 @@ bool Foam::fileFormats::NASedgeFormat::read
while (is.good())
{
string::size_type linei = 0; // parsing position within current line
string line;
is.getLine(line);
if (line.empty() || line[0] == '$')
if (line.empty())
{
continue; // Skip empty or comment
continue; // Ignore empty
}
else if (line[0] == '$')
{
// Ignore comment
continue;
}
// Check if character 72 is continuation
@ -94,41 +98,69 @@ bool Foam::fileFormats::NASedgeFormat::read
}
// Parsing position within current line
std::string::size_type linei = 0;
// Is free format if line contains a comma
const bool freeFormat = line.contains(',');
// First word (column 0-8)
const word cmd(word::validate(nextNasField(line, linei, 8)));
if (cmd == "CBEAM" || cmd == "CROD")
{
// discard elementId (8-16)
(void) nextNasField(line, linei, 8); // 8-16
// discard groupId (16-24)
(void) nextNasField(line, linei, 8); // 16-24
// Fixed format:
// 8-16 : element id
// 16-24 : group id
// 24-32 : vertex
// 32-40 : vertex
label a = readLabel(nextNasField(line, linei, 8)); // 24-32
label b = readLabel(nextNasField(line, linei, 8)); // 32-40
// discard elementId
(void) nextNasField(line, linei, 8, freeFormat);
// discard groupId
(void) nextNasField(line, linei, 8, freeFormat);
dynEdges.append(edge(a,b));
label a = readLabel(nextNasField(line, linei, 8, freeFormat));
label b = readLabel(nextNasField(line, linei, 8, freeFormat));
dynEdges.emplace_back(a,b);
}
else if (cmd == "PLOTEL")
{
// Fixed format:
// 8-16 : element id
// 16-24 : vertex
// 24-32 : vertex
// 32-40 : vertex
// discard elementId (8-16)
(void) nextNasField(line, linei, 8); // 8-16
(void) nextNasField(line, linei, 8, freeFormat);
label a = readLabel(nextNasField(line, linei, 8)); // 16-24
label b = readLabel(nextNasField(line, linei, 8)); // 24-32
label a = readLabel(nextNasField(line, linei, 8, freeFormat));
label b = readLabel(nextNasField(line, linei, 8, freeFormat));
dynEdges.append(edge(a,b));
dynEdges.emplace_back(a,b);
}
else if (cmd == "GRID")
{
label index = readLabel(nextNasField(line, linei, 8)); // 8-16
(void) nextNasField(line, linei, 8); // 16-24
scalar x = readNasScalar(nextNasField(line, linei, 8)); // 24-32
scalar y = readNasScalar(nextNasField(line, linei, 8)); // 32-40
scalar z = readNasScalar(nextNasField(line, linei, 8)); // 40-48
// Fixed (short) format:
// 8-16 : point id
// 16-24 : coordinate system (unsupported)
// 24-32 : point x coordinate
// 32-40 : point y coordinate
// 40-48 : point z coordinate
// 48-56 : displacement coordinate system (optional, unsupported)
// 56-64 : single point constraints (optional, unsupported)
// 64-70 : super-element id (optional, unsupported)
pointId.append(index);
dynPoints.append(point(x, y, z));
label index = readLabel(nextNasField(line, linei, 8, freeFormat));
(void) nextNasField(line, linei, 8, freeFormat);
scalar x = readNasScalar(nextNasField(line, linei, 8, freeFormat));
scalar y = readNasScalar(nextNasField(line, linei, 8, freeFormat));
scalar z = readNasScalar(nextNasField(line, linei, 8, freeFormat));
pointId.push_back(index);
dynPoints.emplace_back(x, y, z);
}
else if (cmd == "GRID*")
{
@ -138,6 +170,8 @@ bool Foam::fileFormats::NASedgeFormat::read
// GRID* 126 0 -5.55999875E+02 -5.68730474E+02
// * 2.14897901E+02
// Cannot be long format and free format at the same time!
label index = readLabel(nextNasField(line, linei, 16)); // 8-24
(void) nextNasField(line, linei, 16); // 24-40
scalar x = readNasScalar(nextNasField(line, linei, 16)); // 40-56
@ -157,8 +191,8 @@ bool Foam::fileFormats::NASedgeFormat::read
(void) nextNasField(line, linei, 8); // 0-8
scalar z = readNasScalar(nextNasField(line, linei, 16)); // 8-16
pointId.append(index);
dynPoints.append(point(x, y, z));
pointId.push_back(index);
dynPoints.emplace_back(x, y, z);
}
}

View File

@ -136,82 +136,6 @@ void Foam::designVariablesUpdate::writeToCostFile(bool zeroAdjointSolns)
}
void Foam::designVariablesUpdate::checkConvergence
(
const scalarField& oldCorrection
)
{
bool converged(false);
// Design variables convergence check
if (designVarsThreshold_)
{
const labelList& activeVarIDs =
designVars_->activeDesignVariables();
const scalarField correction(oldCorrection, activeVarIDs);
const scalarField activeVars(designVars_->getVars(), activeVarIDs);
const scalar scaledCorrection =
gMax(mag(correction)/(mag(activeVars) + SMALL));
DebugInfo
<< "Current correction " << correction << nl
<< "Active vars " << activeVars << endl;
Info<< "Max. scaled correction of the design variables = "
<< scaledCorrection
<< endl;
if (scaledCorrection < designVarsThreshold_())
{
Info<< tab << "Design variables have converged " << endl;
converged = true;
}
}
// Objective convergence check
if (objectiveThreshold_)
{
const scalar newObjective = computeObjectiveValue();
const scalar oldObjective = updateMethod_->getObjectiveValueOld();
const scalar relativeUpdate =
mag(newObjective - oldObjective)/(mag(oldObjective) + SMALL);
Info<< "Relative change of the objective value = "
<< relativeUpdate
<< endl;
if (relativeUpdate < objectiveThreshold_())
{
Info<< tab << "Objective function has converged " << endl;
converged = true;
}
}
// Feasibility check
const scalarField& constraints = updateMethod_->getConstraintValues();
const scalar feasibility = sum(pos(constraints)*constraints);
Info<< "Feasibility = " << feasibility << endl;
if (converged && feasibility < feasibilityThreshold_)
{
Info<< "Stopping criteria met and all constraints satisfied." << nl
<< "Optimisation has converged, stopping ..." << nl << nl
<< "End" << nl
<< endl;
// Force writing of all objective and constraint functions, to get
// the converged results to files
for (adjointSolverManager& am : adjointSolvManagers_)
{
for (adjointSolver& as : am.adjointSolvers())
{
// Use dummy weighted objective
as.getObjectiveManager().writeObjectives();
}
}
writeToCostFile(true);
if (UPstream::parRun())
{
UPstream::exit(0);
}
else
{
std::exit(0);
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::designVariablesUpdate::designVariablesUpdate
@ -570,10 +494,86 @@ void Foam::designVariablesUpdate::postUpdate(const scalarField& oldCorrection)
}
}
}
checkConvergence();
}
if (convergenceCriteriaDefined_)
}
void Foam::designVariablesUpdate::checkConvergence()
{
if (!convergenceCriteriaDefined_)
{
checkConvergence(oldCorrection);
return;
}
bool converged(false);
const scalarField& oldCorrection = updateMethod_->returnCorrection();
// Design variables convergence check
if (designVarsThreshold_)
{
const labelList& activeVarIDs =
designVars_->activeDesignVariables();
const scalarField correction(oldCorrection, activeVarIDs);
const scalarField activeVars(designVars_->getVars(), activeVarIDs);
const scalar scaledCorrection =
gMax(mag(correction)/(mag(activeVars) + SMALL));
DebugInfo
<< "Current correction " << correction << nl
<< "Active vars " << activeVars << endl;
Info<< "Max. scaled correction of the design variables = "
<< scaledCorrection
<< endl;
if (scaledCorrection < designVarsThreshold_())
{
Info<< tab << "Design variables have converged " << endl;
converged = true;
}
}
// Objective convergence check
if (objectiveThreshold_ && updateMethod_->getObjectiveValueOld())
{
const scalar newObjective = computeObjectiveValue();
const scalar oldObjective = updateMethod_->getObjectiveValueOld()();
const scalar relativeUpdate =
mag(newObjective - oldObjective)/(mag(oldObjective) + SMALL);
Info<< "Relative change of the objective value = "
<< relativeUpdate
<< endl;
if (relativeUpdate < objectiveThreshold_())
{
Info<< tab << "Objective function has converged " << endl;
converged = true;
}
}
// Feasibility check
const scalarField& constraints = updateMethod_->getConstraintValues();
const scalar feasibility = sum(pos(constraints)*constraints);
Info<< "Feasibility = " << feasibility << endl;
if (converged && feasibility < feasibilityThreshold_)
{
Info<< "Stopping criteria met and all constraints satisfied." << nl
<< "Optimisation has converged, stopping ..." << nl << nl
<< "End" << nl
<< endl;
// Force writing of all objective and constraint functions, to get
// the converged results to files
for (adjointSolverManager& am : adjointSolvManagers_)
{
for (adjointSolver& as : am.adjointSolvers())
{
// Use dummy weighted objective
as.getObjectiveManager().writeObjectives();
}
}
writeToCostFile(true);
if (UPstream::parRun())
{
UPstream::exit(0);
}
else
{
std::exit(0);
}
}
}

View File

@ -138,11 +138,6 @@ protected:
//- Write to cost file
void writeToCostFile(bool zeroAdjointSolns = false);
//- Check if the optimisation loop has converged based on the provided
//- criteria
// May terminate the program
void checkConvergence(const scalarField& oldCorrection);
private:
@ -216,6 +211,11 @@ public:
//- in the fixedStep approach
void postUpdate(const scalarField& oldCorrection);
//- Check if the optimisation loop has converged based on the provided
//- criteria
// May terminate the program
void checkConvergence();
//- Get access to design variables
inline autoPtr<designVariables>& getDesignVariables()
{

View File

@ -172,6 +172,9 @@ void Foam::optimisationManager::fixedStepUpdate()
// Solve primal equations
solvePrimalEquations();
// Check the convergence criteria of the optimisation loop
dvUpdate_->checkConvergence();
// Reset adjoint sensitivities in all adjoint solver managers
clearSensitivities();

View File

@ -237,7 +237,7 @@ Foam::updateMethod::updateMethod
objectiveDerivatives_(designVars().getVars().size(), Zero),
constraintDerivatives_(0),
objectiveValue_(0),
objectiveValueOld_(0),
objectiveValueOld_(nullptr),
cValues_(0),
correction_(readOrZeroField("correction", designVars().getVars().size())),
cumulativeCorrection_(0),
@ -334,7 +334,11 @@ void Foam::updateMethod::setObjectiveValue(const scalar value)
void Foam::updateMethod::setObjectiveValueOld(const scalar value)
{
objectiveValueOld_ = value;
if (!objectiveValueOld_)
{
objectiveValueOld_.reset(new scalar(Zero));
}
objectiveValueOld_.ref() = value;
}
@ -350,7 +354,8 @@ Foam::scalar Foam::updateMethod::getObjectiveValue() const
}
Foam::scalar Foam::updateMethod::getObjectiveValueOld() const
const Foam::autoPtr<Foam::scalar>&
Foam::updateMethod::getObjectiveValueOld() const
{
return objectiveValueOld_;
}

View File

@ -84,7 +84,7 @@ protected:
//- Old objective value
// Used for convergence checking
scalar objectiveValueOld_;
autoPtr<scalar> objectiveValueOld_;
//- Constraint values
scalarField cValues_;
@ -239,7 +239,7 @@ public:
scalar getObjectiveValue() const;
//- Get old objective value
scalar getObjectiveValueOld() const;
const autoPtr<scalar>& getObjectiveValueOld() const;
//- Get values of constraints
const scalarField& getConstraintValues() const;

View File

@ -183,8 +183,6 @@ bool Foam::fileFormats::NASsurfaceFormat<Face>::read
while (is.good())
{
// Parsing position within current line
std::string::size_type linei = 0;
is.getLine(line);
if (NASCore::debug > 1) Info<< "Process: " << line << nl;
@ -320,16 +318,30 @@ bool Foam::fileFormats::NASsurfaceFormat<Face>::read
}
}
// Parsing position within current line
std::string::size_type linei = 0;
// Is free format if line contains a comma
const bool freeFormat = line.contains(',');
// First word (column 0-8)
const word cmd(word::validate(nextNasField(line, linei, 8)));
if (cmd == "CTRIA3")
{
label elemId = readLabel(nextNasField(line, linei, 8)); // 8-16
label groupId = readLabel(nextNasField(line, linei, 8)); // 16-24
const auto a = readLabel(nextNasField(line, linei, 8)); // 24-32
const auto b = readLabel(nextNasField(line, linei, 8)); // 32-40
const auto c = readLabel(nextNasField(line, linei, 8)); // 40-48
// Fixed format:
// 8-16 : element id
// 16-24 : group id
// 24-32 : vertex
// 32-40 : vertex
// 40-48 : vertex
label elemId = readLabel(nextNasField(line, linei, 8, freeFormat));
label groupId = readLabel(nextNasField(line, linei, 8, freeFormat));
const auto a = readLabel(nextNasField(line, linei, 8, freeFormat));
const auto b = readLabel(nextNasField(line, linei, 8, freeFormat));
const auto c = readLabel(nextNasField(line, linei, 8, freeFormat));
// Convert groupId into zoneId
const auto iterZone = zoneLookup.cfind(groupId);
@ -358,12 +370,20 @@ bool Foam::fileFormats::NASsurfaceFormat<Face>::read
}
else if (cmd == "CQUAD4")
{
label elemId = readLabel(nextNasField(line, linei, 8)); // 8-16
label groupId = readLabel(nextNasField(line, linei, 8)); // 16-24
const auto a = readLabel(nextNasField(line, linei, 8)); // 24-32
const auto b = readLabel(nextNasField(line, linei, 8)); // 32-40
const auto c = readLabel(nextNasField(line, linei, 8)); // 40-48
const auto d = readLabel(nextNasField(line, linei, 8)); // 48-56
// Fixed format:
// 8-16 : element id
// 16-24 : group id
// 24-32 : vertex
// 32-40 : vertex
// 40-48 : vertex
// 48-56 : vertex
label elemId = readLabel(nextNasField(line, linei, 8, freeFormat));
label groupId = readLabel(nextNasField(line, linei, 8, freeFormat));
const auto a = readLabel(nextNasField(line, linei, 8, freeFormat));
const auto b = readLabel(nextNasField(line, linei, 8, freeFormat));
const auto c = readLabel(nextNasField(line, linei, 8, freeFormat));
const auto d = readLabel(nextNasField(line, linei, 8, freeFormat));
// Convert groupId into zoneId
const auto iterZone = zoneLookup.cfind(groupId);
@ -407,11 +427,21 @@ bool Foam::fileFormats::NASsurfaceFormat<Face>::read
}
else if (cmd == "GRID")
{
label index = readLabel(nextNasField(line, linei, 8)); // 8-16
(void) nextNasField(line, linei, 8); // 16-24
scalar x = readNasScalar(nextNasField(line, linei, 8)); // 24-32
scalar y = readNasScalar(nextNasField(line, linei, 8)); // 32-40
scalar z = readNasScalar(nextNasField(line, linei, 8)); // 40-48
// Fixed (short) format:
// 8-16 : point id
// 16-24 : coordinate system (not supported)
// 24-32 : point x coordinate
// 32-40 : point y coordinate
// 40-48 : point z coordinate
// 48-56 : displacement coordinate system (optional, unsupported)
// 56-64 : single point constraints (optional, unsupported)
// 64-70 : super-element id (optional, unsupported)
label index = readLabel(nextNasField(line, linei, 8, freeFormat));
(void) nextNasField(line, linei, 8, freeFormat);
scalar x = readNasScalar(nextNasField(line, linei, 8, freeFormat));
scalar y = readNasScalar(nextNasField(line, linei, 8, freeFormat));
scalar z = readNasScalar(nextNasField(line, linei, 8, freeFormat));
pointId.push_back(index);
dynPoints.emplace_back(x, y, z);
@ -424,6 +454,8 @@ bool Foam::fileFormats::NASsurfaceFormat<Face>::read
// GRID* 126 0 -5.55999875E+02 -5.68730474E+02
// * 2.14897901E+02
// Cannot be long format and free format at the same time!
label index = readLabel(nextNasField(line, linei, 16)); // 8-24
(void) nextNasField(line, linei, 16); // 24-40
scalar x = readNasScalar(nextNasField(line, linei, 16)); // 40-56
@ -453,7 +485,10 @@ bool Foam::fileFormats::NASsurfaceFormat<Face>::read
// have the 'weird' format where the immediately preceeding
// comment contains the information.
label groupId = readLabel(nextNasField(line, linei, 8)); // 8-16
// Fixed format:
// 8-16 : pshell id
label groupId = readLabel(nextNasField(line, linei, 8, freeFormat));
if (lastComment.size() > 1 && !nameLookup.contains(groupId))
{

View File

@ -1,7 +1,7 @@
#------------------------------------------------------------------------------
include $(GENERAL_RULES)/Icx/c++
c++ARCH := -fp-model=precise
c++ARCH := -pthread -fp-model=precise
ifneq (,$(strip $(WM_COMPILE_OPTION)))
sinclude $(DEFAULT_RULES)/c++$(WM_COMPILE_OPTION)