Skip to content

Commit

Permalink
Added heat source parameters as DVs and location snapCenter option.
Browse files Browse the repository at this point in the history
  • Loading branch information
@friedenhe
friedenhe committed Jun 27, 2024
1 parent aa2d351 commit fb5a72a
Show file tree
Hide file tree
Showing 11 changed files with 556 additions and 31 deletions.
50 changes: 50 additions & 0 deletions dafoam/mphys/mphys_dafoam.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -657,6 +657,11 @@ def setup(self):
if "comps" in list(designVariables[dvName].keys()):
nACTDVars = len(designVariables[dvName]["comps"])
self.add_input(dvName, distributed=False, shape=nACTDVars, tags=["mphys_coupling"])
elif dvType == "HSC": # add heat source parameter variables
nHSCVars = 9
if "comps" in list(designVariables[dvName].keys()):
nHSCVars = len(designVariables[dvName]["comps"])
self.add_input(dvName, distributed=False, shape=nHSCVars, tags=["mphys_coupling"])
elif dvType == "Field": # add field variables
# user can prescribe whether the field var is distributed. Default is True
distributed = True
Expand Down Expand Up @@ -888,6 +893,26 @@ def apply_linear(self, inputs, outputs, d_inputs, d_outputs, d_residuals, mode):
d_inputs[inputName] += ACTDBarSub
else:
d_inputs[inputName] += ACTDBar

# compute [dRdHSC]^T*Psi using reverse mode AD
elif self.dvType[inputName] == "HSC":
prodVec = PETSc.Vec().create(self.comm)
prodVec.setSizes((PETSc.DECIDE, 9), bsize=1)
prodVec.setFromOptions()
DASolver.solverAD.calcdRdHSCTPsiAD(
DASolver.xvVec, DASolver.wVec, resBarVec, inputName.encode(), prodVec
)
# we will convert the MPI prodVec to seq array for all procs
HSCBar = DASolver.convertMPIVec2SeqArray(prodVec)
if "comps" in list(designVariables[inputName].keys()):
nHSCVars = len(designVariables[inputName]["comps"])
HSCBarSub = np.zeros(nHSCVars, "d")
for i in range(nHSCVars):
comp = designVariables[inputName]["comps"][i]
HSCBarSub[i] = HSCBar[comp]
d_inputs[inputName] += HSCBarSub
else:
d_inputs[inputName] += HSCBar

# compute dRdFieldT*Psi using reverse mode AD
elif self.dvType[inputName] == "Field":
Expand Down Expand Up @@ -1251,6 +1276,11 @@ def setup(self):
if "comps" in list(designVariables[dvName].keys()):
nACTDVars = len(designVariables[dvName]["comps"])
self.add_input(dvName, distributed=False, shape=nACTDVars, tags=["mphys_coupling"])
elif dvType == "HSC": # add heat source parameter variables
nHSCVars = 9
if "comps" in list(designVariables[dvName].keys()):
nHSCVars = len(designVariables[dvName]["comps"])
self.add_input(dvName, distributed=False, shape=nHSCVars, tags=["mphys_coupling"])
elif dvType == "Field": # add field variables
# user can prescribe whether the field var is distributed. Default is True
distributed = True
Expand Down Expand Up @@ -1447,6 +1477,26 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode):
d_inputs[inputName] += ACTDBarSub * fBar
else:
d_inputs[inputName] += ACTDBar * fBar

# compute dFdHSC
elif self.dvType[inputName] == "HSC":
dFdHSC = PETSc.Vec().create(self.comm)
dFdHSC.setSizes((PETSc.DECIDE, 9), bsize=1)
dFdHSC.setFromOptions()
DASolver.solverAD.calcdFdHSCAD(
DASolver.xvVec, DASolver.wVec, objFuncName.encode(), inputName.encode(), dFdHSC
)
# we will convert the MPI dFdHSC to seq array for all procs
HSCBar = DASolver.convertMPIVec2SeqArray(dFdHSC)
if "comps" in list(designVariables[inputName].keys()):
nHSCVars = len(designVariables[inputName]["comps"])
HSCBarSub = np.zeros(nHSCVars, "d")
for i in range(nHSCVars):
comp = designVariables[inputName]["comps"][i]
HSCBarSub[i] = HSCBar[comp]
d_inputs[inputName] += HSCBarSub * fBar
else:
d_inputs[inputName] += HSCBar * fBar

# compute dFdField
elif self.dvType[inputName] == "Field":
Expand Down
2 changes: 1 addition & 1 deletion dafoam/pyDAFoam.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@
"""

__version__ = "3.1.2"
__version__ = "3.1.3"

import subprocess
import os
Expand Down
37 changes: 37 additions & 0 deletions src/adjoint/DAFvSource/DAFvSource.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -175,6 +175,43 @@ void DAFvSource::syncDAOptionToActuatorDVs()
dvList[11] = diskSubDict.getScalar("expN");
dvList[12] = diskSubDict.getScalar("targetThrust");

// set actuatorDiskDVs_
actuatorDiskDVs_.set(diskName, dvList);
}
}
}
else if (type0 == "heatSource")
{
word source0 = fvSourceSubDict.subDict(diskName0).getWord("source");

if (source0 == "cylinderSmooth")
{
forAll(fvSourceSubDict.toc(), idxI)
{
word diskName = fvSourceSubDict.toc()[idxI];

// sub dictionary with all parameters for this disk
dictionary diskSubDict = fvSourceSubDict.subDict(diskName);

// now read in all parameters for this actuator disk
scalarList centerList;
diskSubDict.readEntry<scalarList>("center", centerList);

scalarList axisList;
diskSubDict.readEntry<scalarList>("axis", axisList);

// we have 13 design variables for each disk
scalarList dvList(9);
dvList[0] = centerList[0];
dvList[1] = centerList[1];
dvList[2] = centerList[2];
dvList[3] = axisList[0];
dvList[4] = axisList[1];
dvList[5] = axisList[2];
dvList[6] = diskSubDict.getScalar("radius");
dvList[7] = diskSubDict.getScalar("length");
dvList[8] = diskSubDict.getScalar("power");

// set actuatorDiskDVs_
actuatorDiskDVs_.set(diskName, dvList);
}
Expand Down
96 changes: 71 additions & 25 deletions src/adjoint/DAFvSource/DAFvSourceHeatSource.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,8 @@ DAFvSourceHeatSource::DAFvSourceHeatSource(
const DAIndex& daIndex)
: DAFvSource(modelType, mesh, daOption, daModel, daIndex)
{
printInterval_ = daOption.getOption<label>("printInterval");

const dictionary& allOptions = daOption_.getAllOptions();

dictionary fvSourceSubDict = allOptions.subDict("fvSource");
Expand All @@ -34,17 +36,18 @@ DAFvSourceHeatSource::DAFvSourceHeatSource(
dictionary sourceSubDict = fvSourceSubDict.subDict(sourceName);
word sourceType = sourceSubDict.getWord("source");

scalarList p1List;
sourceSubDict.readEntry<scalarList>("p1", p1List);
cylinderP1_.set(sourceName, p1List);
scalarList p2List;
sourceSubDict.readEntry<scalarList>("p2", p2List);
cylinderP2_.set(sourceName, p2List);
cylinderRadius_.set(sourceName, sourceSubDict.getScalar("radius"));
power_.set(sourceName, sourceSubDict.getScalar("power"));

if (sourceType == "cylinderToCell")
{
scalarList p1List;
sourceSubDict.readEntry<scalarList>("p1", p1List);
cylinderP1_.set(sourceName, p1List);
scalarList p2List;
sourceSubDict.readEntry<scalarList>("p2", p2List);
cylinderP2_.set(sourceName, p2List);

cylinderRadius_.set(sourceName, sourceSubDict.getScalar("radius"));
power_.set(sourceName, sourceSubDict.getScalar("power"));

fvSourceCellIndices_.set(sourceName, {});
// all available source type are in src/meshTools/sets/cellSources
// Example of IO parameters os in applications/utilities/mesh/manipulation/topoSet
Expand Down Expand Up @@ -92,17 +95,26 @@ DAFvSourceHeatSource::DAFvSourceHeatSource(
Info << "fvSourceCellIndices " << fvSourceCellIndices_ << endl;
}
}
else if (sourceType == "cylinder")
else if (sourceType == "cylinderSmooth")
{
// no need to do special things

// eps is a smoothing parameter, it should be the local mesh cell size in meters
// near the cylinder region
cylinderEps_.set(sourceName, sourceSubDict.getScalar("eps"));
}
else
{
FatalErrorIn("DAFvSourceHeatSource") << "source: " << sourceType << " not supported!"
<< "Options are: cylinderAnnulusToCell and cylinder!"
<< "Options are: cylinderAnnulusToCell and cylinderSmooth!"
<< abort(FatalError);
}
}

// now we need to initialize actuatorDiskDVs_ by synchronizing the values
// defined in fvSource from DAOption to actuatorDiskDVs_
// NOTE: we need to call this function whenever we change the actuator
// design variables during optimization
this->syncDAOptionToActuatorDVs();
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

Expand Down Expand Up @@ -184,20 +196,31 @@ void DAFvSourceHeatSource::calcFvSource(volScalarField& fvSource)

reduce(sourceTotal, sumOp<scalar>());

Info << "Total volume for " << sourceName << ": " << totalV << " m^3" << endl;
Info << "Total heat source for " << sourceName << ": " << sourceTotal << " W." << endl;
if (daOption_.getOption<word>("runStatus") == "solvePrimal")
{
if (mesh_.time().timeIndex() % printInterval_ == 0 || mesh_.time().timeIndex() == 1)
{
Info << "Total volume for " << sourceName << ": " << totalV << " m^3" << endl;
Info << "Total heat source for " << sourceName << ": " << sourceTotal << " W." << endl;
}
}
}
else if (sourceType == "cylinder")
else if (sourceType == "cylinderSmooth")
{
vector p1 = {cylinderP1_[sourceName][0], cylinderP1_[sourceName][1], cylinderP1_[sourceName][2]};
vector p2 = {cylinderP2_[sourceName][0], cylinderP2_[sourceName][1], cylinderP2_[sourceName][2]};
vector cylinderCenter = (p1 + p2) / 2.0;
vector cylinderDir = p2 - p1;
scalar cylinderLen = mag(cylinderDir);
vector cylinderDirNorm = cylinderDir / cylinderLen;
scalar radius = cylinderRadius_[sourceName];
vector cylinderCenter =
{actuatorDiskDVs_[sourceName][0], actuatorDiskDVs_[sourceName][1], actuatorDiskDVs_[sourceName][2]};
vector cylinderDir =
{actuatorDiskDVs_[sourceName][3], actuatorDiskDVs_[sourceName][4], actuatorDiskDVs_[sourceName][5]};
scalar radius = actuatorDiskDVs_[sourceName][6];
scalar cylinderLen = actuatorDiskDVs_[sourceName][7];
scalar power = actuatorDiskDVs_[sourceName][8];
vector cylinderDirNorm = cylinderDir / mag(cylinderDir);
scalar eps = cylinderEps_[sourceName];

scalar totalVolume = constant::mathematical::pi * radius * radius * cylinderLen;

scalar sourceTotal = 0.0;
scalar volumeTotal = 0.0;
forAll(mesh_.cells(), cellI)
{
// the cell center coordinates of this cellI
Expand All @@ -221,16 +244,39 @@ void DAFvSourceHeatSource::calcFvSource(volScalarField& fvSource)
// the magnitude of axial component of cellC2AVecR
scalar cellC2AVecALen = mag(cellC2AVecA);

if (cellC2AVecRLen <= radius && cellC2AVecALen <= cylinderLen / 2.0)
scalar axialSmoothC = 1.0;
scalar radialSmoothC = 1.0;

if (cellC2AVecALen > cylinderLen / 2.0)
{
scalar d2 = (cellC2AVecALen - cylinderLen / 2.0) * (cellC2AVecALen - cylinderLen / 2.0);
axialSmoothC = exp(-d2 / eps / eps);
}

if (cellC2AVecRLen > radius)
{
scalar d2 = (cellC2AVecRLen - radius) * (cellC2AVecRLen - radius);
radialSmoothC = exp(-d2 / eps / eps);
}

fvSource[cellI] += power / totalVolume * axialSmoothC * radialSmoothC;
sourceTotal += power / totalVolume * axialSmoothC * radialSmoothC * mesh_.V()[cellI];
volumeTotal += axialSmoothC * radialSmoothC * mesh_.V()[cellI];
}

if (daOption_.getOption<word>("runStatus") == "solvePrimal")
{
if (mesh_.time().timeIndex() % printInterval_ == 0 || mesh_.time().timeIndex() == 1)
{
fvSource[cellI] += power / totalVolume;
Info << "Total volume for " << sourceName << ": " << volumeTotal << " m^3" << endl;
Info << "Total heat source for " << sourceName << ": " << sourceTotal << " W." << endl;
}
}
}
else
{
FatalErrorIn("DAFvSourceHeatSource") << "source: " << sourceType << " not supported!"
<< "Options are: cylinderAnnulusToCell and cylinder!"
<< "Options are: cylinderAnnulusToCell and cylinderSmooth!"
<< abort(FatalError);
}
}
Expand Down
5 changes: 3 additions & 2 deletions src/adjoint/DAFvSource/DAFvSourceHeatSource.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,9 +35,10 @@ protected:
HashTable<scalarList> cylinderP2_;
HashTable<scalar> cylinderRadius_;
HashTable<scalar> power_;
HashTable<scalar> cylinderEps_;

/// calculate DAFvSourceHeatSource::fvSourceCellIndices_
void calcFvSourceCellIndices(const word sourceName, HashTable<labelList>& fvSourceCellIndices);
/// print interval for primal and adjoint
label printInterval_;

public:
TypeName("heatSource");
Expand Down
36 changes: 33 additions & 3 deletions src/adjoint/DAObjFunc/DAObjFuncLocation.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,6 +65,17 @@ DAObjFuncLocation::DAObjFuncLocation(
center_ = {centerRead[0], centerRead[1], centerRead[2]};
}

snapCenter2Cell_ = objFuncDict_.lookupOrDefault<label>("snapCenter2Cell", 0);
if (snapCenter2Cell_)
{
point centerPoint = {center_[0], center_[1], center_[2]};
snappedCenterCellI_ = mesh_.findCell(centerPoint);
if (snappedCenterCellI_ < 0)
{
FatalErrorIn(" ") << "can not find cell for center " << centerPoint << abort(FatalError);
}
}

if (mode_ == "maxRadius")
{
// we need to identify the patchI and faceI that has maxR
Expand Down Expand Up @@ -158,7 +169,13 @@ void DAObjFuncLocation::calcObjFunc(
const label patchI = daIndex_.bFacePatchI[bFaceI];
const label faceI = daIndex_.bFaceFaceI[bFaceI];

vector faceC = mesh_.Cf().boundaryField()[patchI][faceI] - center_;
vector center = center_;
if (snapCenter2Cell_)
{
center = mesh_.C()[snappedCenterCellI_];
}

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

tensor faceCTensor(tensor::zero);
faceCTensor.xx() = faceC.x();
Expand Down Expand Up @@ -200,7 +217,13 @@ void DAObjFuncLocation::calcObjFunc(
const label patchI = daIndex_.bFacePatchI[bFaceI];
const label faceI = daIndex_.bFaceFaceI[bFaceI];

vector faceC = mesh_.Cf().boundaryField()[patchI][faceI] - center_;
vector center = center_;
if (snapCenter2Cell_)
{
center = mesh_.C()[snappedCenterCellI_];
}

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

tensor faceCTensor(tensor::zero);
faceCTensor.xx() = faceC.x();
Expand Down Expand Up @@ -234,7 +257,14 @@ void DAObjFuncLocation::calcObjFunc(
scalar radius = 0.0;
if (maxRPatchI_ >= 0 && maxRFaceI_ >= 0)
{
vector faceC = mesh_.Cf().boundaryField()[maxRPatchI_][maxRFaceI_] - center_;

vector center = center_;
if (snapCenter2Cell_)
{
center = mesh_.C()[snappedCenterCellI_];
}

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

tensor faceCTensor(tensor::zero);
faceCTensor.xx() = faceC.x();
Expand Down
6 changes: 6 additions & 0 deletions src/adjoint/DAObjFunc/DAObjFuncLocation.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,12 @@ protected:
/// maxRadius patchI
label maxRFaceI_ = -1;

/// whether to snap the center to a cell in the mesh if yes the center will move with the mesh
label snapCenter2Cell_ = 0;

/// the cell index for the center if snapCenter2Cell_ = 1
label snappedCenterCellI_ = -1;

public:
TypeName("location");
// Constructors
Expand Down
Loading

0 comments on commit fb5a72a

Please sign in to comment.