Skip to content

Commit

Permalink
Fixed a few issues for snapCenter2Cell (#658)
Browse files Browse the repository at this point in the history 
* Fixed bugs for snapCenter2Cell and hasFvSource.

* Fixed a bug.

* More bug fix.

* Fixed bugs in snapCenter.
  • Loading branch information
@friedenhe
friedenhe authored Jun 30, 2024
1 parent c195f8b commit c27805a
Show file tree
Hide file tree
Showing 8 changed files with 86 additions and 33 deletions.
13 changes: 11 additions & 2 deletions src/adjoint/DAFvSource/DAFvSourceHeatSource.C
View file
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
View file
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
View file
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
View file
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
View file
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
View file
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
View file
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
View file
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

0 comments on commit c27805a

Please sign in to comment.