Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fixed a few issues for snapCenter2Cell #658

Merged
merged 4 commits into from
Jun 30, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
13 changes: 11 additions & 2 deletions src/adjoint/DAFvSource/DAFvSourceHeatSource.C
Original file line number Diff line number Diff line change
Expand Up @@ -112,11 +112,16 @@ DAFvSourceHeatSource::DAFvSourceHeatSource(
if (snapCenter2Cell_[sourceName])
{
point centerPoint = {actuatorDiskDVs_[sourceName][0], actuatorDiskDVs_[sourceName][1], actuatorDiskDVs_[sourceName][2]};
snappedCenterCellI_.set(sourceName, mesh_.findCell(centerPoint));

// NOTE: we need to call a self-defined findCell func to make it work correctly in ADR
label myCellI = DAUtility::myFindCell(mesh_, centerPoint);

snappedCenterCellI_.set(sourceName, myCellI);
label foundCellI = 0;
if (snappedCenterCellI_[sourceName] >= 0)
{
foundCellI = 1;
//Pout << "snap source " << sourceName << " to center " << mesh_.C()[snappedCenterCellI_[sourceName]] << endl;
}
reduce(foundCellI, sumOp<label>());
if (foundCellI != 1)
Expand All @@ -128,6 +133,10 @@ DAFvSourceHeatSource::DAFvSourceHeatSource(
<< " be outside of the mesh domain or on a mesh face "
<< abort(FatalError);
}

vector snappedCenter = vector::zero;
this->findGlobalSnappedCenter(snappedCenterCellI_[sourceName], snappedCenter);
Info << "heat source " << sourceName << " snap to center " << snappedCenter << endl;
}
}
else
Expand Down Expand Up @@ -232,7 +241,7 @@ void DAFvSourceHeatSource::calcFvSource(volScalarField& fvSource)
{
vector cylinderCenter =
{actuatorDiskDVs_[sourceName][0], actuatorDiskDVs_[sourceName][1], actuatorDiskDVs_[sourceName][2]};

if (snapCenter2Cell_[sourceName])
{
this->findGlobalSnappedCenter(snappedCenterCellI_[sourceName], cylinderCenter);
Expand Down
45 changes: 26 additions & 19 deletions src/adjoint/DAObjFunc/DAObjFuncLocation.C
Original file line number Diff line number Diff line change
Expand Up @@ -69,7 +69,10 @@ DAObjFuncLocation::DAObjFuncLocation(
if (snapCenter2Cell_)
{
point centerPoint = {center_[0], center_[1], center_[2]};
snappedCenterCellI_ = mesh_.findCell(centerPoint);

// NOTE: we need to call a self-defined findCell func to make it work correctly in ADR
snappedCenterCellI_ = DAUtility::myFindCell(mesh_, centerPoint);

label foundCellI = 0;
if (snappedCenterCellI_ >= 0)
{
Expand All @@ -85,6 +88,9 @@ DAObjFuncLocation::DAObjFuncLocation(
<< " be outside of the mesh domain or on a mesh face "
<< abort(FatalError);
}
vector snappedCenter = vector::zero;
this->findGlobalSnappedCenter(snappedCenterCellI_, snappedCenter);
Info << "snap to center " << snappedCenter << endl;
}

if (mode_ == "maxRadius")
Expand Down Expand Up @@ -196,19 +202,19 @@ void DAObjFuncLocation::calcObjFunc(
// calculate Location
scalar objValTmp = 0.0;

vector center = center_;
if (snapCenter2Cell_)
{
this->findGlobalSnappedCenter(snappedCenterCellI_, center);
}

forAll(objFuncFaceSources, idxI)
{
const label& objFuncFaceI = objFuncFaceSources[idxI];
label bFaceI = objFuncFaceI - daIndex_.nLocalInternalFaces;
const label patchI = daIndex_.bFacePatchI[bFaceI];
const label faceI = daIndex_.bFaceFaceI[bFaceI];

vector center = center_;
if (snapCenter2Cell_)
{
this->findGlobalSnappedCenter(snappedCenterCellI_, center);
}

vector faceC = mesh_.Cf().boundaryField()[patchI][faceI] - center;

tensor faceCTensor(tensor::zero);
Expand Down Expand Up @@ -244,19 +250,19 @@ void DAObjFuncLocation::calcObjFunc(
// calculate Location
scalar objValTmp = 0.0;

vector center = center_;
if (snapCenter2Cell_)
{
this->findGlobalSnappedCenter(snappedCenterCellI_, center);
}

forAll(objFuncFaceSources, idxI)
{
const label& objFuncFaceI = objFuncFaceSources[idxI];
label bFaceI = objFuncFaceI - daIndex_.nLocalInternalFaces;
const label patchI = daIndex_.bFacePatchI[bFaceI];
const label faceI = daIndex_.bFaceFaceI[bFaceI];

vector center = center_;
if (snapCenter2Cell_)
{
this->findGlobalSnappedCenter(snappedCenterCellI_, center);
}

vector faceC = mesh_.Cf().boundaryField()[patchI][faceI] - center;

tensor faceCTensor(tensor::zero);
Expand Down Expand Up @@ -289,14 +295,15 @@ void DAObjFuncLocation::calcObjFunc(
else if (mode_ == "maxRadius")
{
scalar radius = 0.0;
if (maxRPatchI_ >= 0 && maxRFaceI_ >= 0)

vector center = center_;
if (snapCenter2Cell_)
{
this->findGlobalSnappedCenter(snappedCenterCellI_, center);
}

vector center = center_;
if (snapCenter2Cell_)
{
this->findGlobalSnappedCenter(snappedCenterCellI_, center);
}
if (maxRPatchI_ >= 0 && maxRFaceI_ >= 0)
{

vector faceC = mesh_.Cf().boundaryField()[maxRPatchI_][maxRFaceI_] - center;

Expand Down
30 changes: 23 additions & 7 deletions src/adjoint/DAObjFunc/DAObjFuncMass.C
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,7 @@ DAObjFuncMass::DAObjFuncMass(

objFuncDict_.readEntry<scalar>("scale", scale_);

rho_ = objFuncDict_.lookupOrDefault<scalar>("rho", -1.0);
}

/// calculate the value of objective function
Expand Down Expand Up @@ -82,15 +83,30 @@ void DAObjFuncMass::calcObjFunc(
objFuncValue = 0.0;

const objectRegistry& db = mesh_.thisDb();
const volScalarField& rho = db.lookupObject<volScalarField>("solid:rho");

// calculate mass
forAll(objFuncCellSources, idxI)
if (rho_ < 0.0)
{
const label& cellI = objFuncCellSources[idxI];
scalar volume = mesh_.V()[cellI];
objFuncCellValues[idxI] = scale_ * volume * rho[cellI];
objFuncValue += objFuncCellValues[idxI];
const volScalarField& rho = db.lookupObject<volScalarField>("solid:rho");

// calculate mass
forAll(objFuncCellSources, idxI)
{
const label& cellI = objFuncCellSources[idxI];
scalar volume = mesh_.V()[cellI];
objFuncCellValues[idxI] = scale_ * volume * rho[cellI];
objFuncValue += objFuncCellValues[idxI];
}
}
else
{
// calculate mass
forAll(objFuncCellSources, idxI)
{
const label& cellI = objFuncCellSources[idxI];
scalar volume = mesh_.V()[cellI];
objFuncCellValues[idxI] = scale_ * volume * rho_;
objFuncValue += objFuncCellValues[idxI];
}
}

// need to reduce the sum of force across all processors
Expand Down
3 changes: 3 additions & 0 deletions src/adjoint/DAObjFunc/DAObjFuncMass.H
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,9 @@ public:
TypeName("mass");
// Constructors

/// user-defined density
scalar rho_ = -1.0;

//- Construct from components
DAObjFuncMass(
const fvMesh& mesh,
Expand Down
4 changes: 2 additions & 2 deletions src/adjoint/DAObjFunc/DAObjFuncVariance.C
Original file line number Diff line number Diff line change
Expand Up @@ -207,7 +207,7 @@ DAObjFuncVariance::DAObjFuncVariance(
forAll(probePointCoords_, idxI)
{
point pointCoord = {probePointCoords_[idxI][0], probePointCoords_[idxI][1], probePointCoords_[idxI][2]};
label cellI = mesh_.findCell(pointCoord);
label cellI = DAUtility::myFindCell(mesh_, pointCoord);
if (cellI >= 0)
{
probeCellIndex_.append(cellI);
Expand Down Expand Up @@ -286,7 +286,7 @@ DAObjFuncVariance::DAObjFuncVariance(
forAll(probePointCoords_, idxI)
{
point pointCoord = {probePointCoords_[idxI][0], probePointCoords_[idxI][1], probePointCoords_[idxI][2]};
label cellI = mesh_.findCell(pointCoord);
label cellI = DAUtility::myFindCell(mesh_, pointCoord);
if (cellI >= 0)
{
probeCellIndex_.append(cellI);
Expand Down
2 changes: 1 addition & 1 deletion src/adjoint/DAResidual/DAResidualHeatTransferFoam.H
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,7 @@ protected:

autoPtr<dimensionedScalar> kPtr_;

label hasFvSource_;
label hasFvSource_ = 0;

public:
TypeName("DAHeatTransferFoam");
Expand Down
16 changes: 15 additions & 1 deletion src/adjoint/DAUtility/DAUtility.C
Original file line number Diff line number Diff line change
Expand Up @@ -770,7 +770,7 @@ void DAUtility::primalResidualControl(
// calculate the initial residual mag and set it to primalResidualNorms_

// for vectors, we need to use the median value for the residual
// this is because we often need to run 2D simulations with symmetry
// this is because we often need to run 2D simulations with symmetry
// BC, so one component of the residual vector, which is related to the symmetry BC,
// may be high while the other two components' residuals are low.
// In this case, we can use the median value for the residual vector, which better
Expand All @@ -792,6 +792,20 @@ void DAUtility::primalResidualControl(
}
}

label DAUtility::myFindCell(
const primitiveMesh& mesh,
const point& point)
{
/*
A self-defined findCell function. We will need to cast fvMesh to primitiveMesh
and then call primitiveMesh's findCell. For some reasons, the fvMesh's findCell
did not work correctly in ADR mode...
*/

label cellI = mesh.findCell(point);
return cellI;
}

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

} // End namespace Foam
Expand Down
6 changes: 5 additions & 1 deletion src/adjoint/DAUtility/DAUtility.H
Original file line number Diff line number Diff line change
Expand Up @@ -144,12 +144,16 @@ public:
const SolverPerformance<scalar>& solverP,
const label printToScreen,
const word varName);

/// control when to print the residual and also compute the maxInitRes
static void primalResidualControl(
const SolverPerformance<vector>& solverP,
const label printToScreen,
const word varName);

static label myFindCell(
const primitiveMesh& mesh,
const point& point);
};

template<class classType>
Expand Down
Loading