Skip to content

Commit

Permalink
Added more heat sources.
Browse files Browse the repository at this point in the history
  • Loading branch information
@friedenhe
friedenhe committed Jun 22, 2024
1 parent 25bd435 commit b981121
Show file tree
Hide file tree
Showing 4 changed files with 153 additions and 136 deletions.
259 changes: 132 additions & 127 deletions src/adjoint/DAFvSource/DAFvSourceHeatSource.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,85 @@ DAFvSourceHeatSource::DAFvSourceHeatSource(
const DAIndex& daIndex)
: DAFvSource(modelType, mesh, daOption, daModel, daIndex)
{
this->calcFvSourceCellIndices(fvSourceCellIndices_);
const dictionary& allOptions = daOption_.getAllOptions();

dictionary fvSourceSubDict = allOptions.subDict("fvSource");

forAll(fvSourceSubDict.toc(), idxI)
{
word sourceName = fvSourceSubDict.toc()[idxI];
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")
{
fvSourceCellIndices_.set(sourceName, {});
// all available source type are in src/meshTools/sets/cellSources
// Example of IO parameters os in applications/utilities/mesh/manipulation/topoSet

// create a topoSet
autoPtr<topoSet> currentSet(
topoSet::New(
"cellSet",
mesh_,
"set0",
IOobject::NO_READ));
// we need to change the min and max because they need to
// be of type point; however, we can't parse point type
// in pyDict, we need to change them here.

point p1;
point p2;
p1[0] = cylinderP1_[sourceName][0];
p1[1] = cylinderP1_[sourceName][1];
p1[2] = cylinderP1_[sourceName][2];
p2[0] = cylinderP2_[sourceName][0];
p2[1] = cylinderP2_[sourceName][1];
p2[2] = cylinderP2_[sourceName][2];

dictionary tmpDict;
tmpDict.set("p1", p1);
tmpDict.set("p2", p2);
tmpDict.set("radius", cylinderRadius_[sourceName]);

// create the source
autoPtr<topoSetSource> sourceSet(
topoSetSource::New("cylinderToCell", mesh_, tmpDict));

// add the sourceSet to topoSet
sourceSet().applyToSet(topoSetSource::NEW, currentSet());
// get the face index from currentSet, we need to use
// this special for loop
for (const label i : currentSet())
{
fvSourceCellIndices_[sourceName].append(i);
}

if (daOption_.getOption<label>("debug"))
{
Info << "fvSourceCellIndices " << fvSourceCellIndices_ << endl;
}
}
else if (sourceType == "cylinder")
{
// no need to do special things
}
else
{
FatalErrorIn("DAFvSourceHeatSource") << "source: " << sourceType << " not supported!"
<< "Options are: cylinderAnnulusToCell and cylinder!"
<< abort(FatalError);
}
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

Expand All @@ -43,21 +121,19 @@ void DAFvSourceHeatSource::calcFvSource(volScalarField& fvSource)
"source1"
{
"type": "heatSource",
"source": "cylinderAnnulusToCell",
"source": "cylinderToCell",
"p1": [0.5, 0.3, 0.1], # p1 and p2 define the axis and width
"p2": [0.5, 0.3, 0.5], # p2-p1 should be the axis of the cylinder
"innerRadius": 0.1,
"outerRadius": 0.8,
"radius": 0.8,
"power": 101.0, # here we should prescribe the power in W
},
"source2"
{
"type": "heatSource",
"source": "cylinderAnnulusToCell",
"source": "cylinder",
"p1": [0.5, 0.3, 0.1], # p1 and p2 define the axis and width
"p2": [0.5, 0.3, 0.5], # p2-p1 should be the axis of the cylinder
"innerRadius": 0.0,
"outerRadius": 1.5,
"radius": 1.5,
"power": 10.5, # here we should prescribe the power in W
}
}
Expand All @@ -73,25 +149,17 @@ void DAFvSourceHeatSource::calcFvSource(volScalarField& fvSource)

dictionary fvSourceSubDict = allOptions.subDict("fvSource");

word sourceName0 = fvSourceSubDict.toc()[0];
word source0 = fvSourceSubDict.subDict(sourceName0).getWord("source");

if (source0 == "cylinderAnnulusToCell")
forAll(fvSourceSubDict.toc(), idxI)
{

// loop over all the cell indices for all heat sources
forAll(fvSourceCellIndices_.toc(), idxI)
word sourceName = fvSourceSubDict.toc()[idxI];
dictionary sourceSubDict = fvSourceSubDict.subDict(sourceName);
word sourceType = sourceSubDict.getWord("source");
// NOTE: here power should be in W. We will evenly divide the power by the
// total volume of the source
scalar power = power_[sourceName];
if (sourceType == "cylinderToCell")
{

// name of this heat source
word sourceName = fvSourceCellIndices_.toc()[idxI];

// sub dictionary with all parameters for this heat source
dictionary sourceSubDict = fvSourceSubDict.subDict(sourceName);
// NOTE: here power should be in W. We will evenly divide the power by the
// total volume of the source
scalar power = sourceSubDict.getScalar("power");

// loop over all cell indices for this source and assign the source term
// ----- first loop, calculate the total volume
scalar totalV = 0.0;
Expand All @@ -117,120 +185,57 @@ 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;
}
}
else
{
FatalErrorIn("calcFvSourceCells") << "source: " << source0 << " not supported!"
<< "Options are: cylinderAnnulusToCell!"
<< abort(FatalError);
}

fvSource.correctBoundaryConditions();
}

void DAFvSourceHeatSource::calcFvSourceCellIndices(HashTable<labelList>& fvSourceCellIndices)
{
/*
Description:
Calculate the lists of cell indices that are within the source space
NOTE: we support multiple heat sources.
Output:
fvSourceCellIndices: Hash table that contains the lists of cell indices that
are within the source space. The hash key is the name of the source. We support
multiple sources. An example of fvSourceCellIndices can be:
fvSourceCellIndices=
{
"source1": {1,2,3,4,5},
"source2": {6,7,8,9,10,11,12}
Info << "Total heat source for " << sourceName << ": " << sourceTotal << " W." << endl;
}
*/

const dictionary& allOptions = daOption_.getAllOptions();

dictionary fvSourceSubDict = allOptions.subDict("fvSource");

if (fvSourceSubDict.toc().size() == 0)
{
FatalErrorIn("calcSourceCells") << "heatSource is selected as fvSource "
<< " but the options are empty!"
<< abort(FatalError);
}

forAll(fvSourceSubDict.toc(), idxI)
{
word sourceName = fvSourceSubDict.toc()[idxI];

fvSourceCellIndices.set(sourceName, {});

dictionary sourceSubDict = fvSourceSubDict.subDict(sourceName);
word sourceType = sourceSubDict.getWord("source");
// all available source type are in src/meshTools/sets/cellSources
// Example of IO parameters os in applications/utilities/mesh/manipulation/topoSet
if (sourceType == "cylinderAnnulusToCell")
else if (sourceType == "cylinder")
{
// create a topoSet
autoPtr<topoSet> currentSet(
topoSet::New(
"cellSet",
mesh_,
"set0",
IOobject::NO_READ));
// we need to change the min and max because they need to
// be of type point; however, we can't parse point type
// in pyDict, we need to change them here.

scalarList point1;
scalarList point2;
sourceSubDict.readEntry<scalarList>("p1", point1);
sourceSubDict.readEntry<scalarList>("p2", point2);

point p1;
point p2;
p1[0] = point1[0];
p1[1] = point1[1];
p1[2] = point1[2];
p2[0] = point2[0];
p2[1] = point2[1];
p2[2] = point2[2];

scalar outerRadius = sourceSubDict.getScalar("outerRadius");
scalar innerRadius = sourceSubDict.getScalar("innerRadius");

dictionary tmpDict;
tmpDict.set("p1", p1);
tmpDict.set("p2", p2);
tmpDict.set("innerRadius", innerRadius);
tmpDict.set("outerRadius", outerRadius);

// create the source
autoPtr<topoSetSource> sourceSet(
topoSetSource::New(sourceType, mesh_, tmpDict));

// add the sourceSet to topoSet
sourceSet().applyToSet(topoSetSource::NEW, currentSet());
// get the face index from currentSet, we need to use
// this special for loop
for (const label i : currentSet())
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];
scalar totalVolume = constant::mathematical::pi * radius * radius * cylinderLen;

forAll(mesh_.cells(), cellI)
{
fvSourceCellIndices[sourceName].append(i);
// the cell center coordinates of this cellI
vector cellC = mesh_.C()[cellI];
// cell center to disk center vector
vector cellC2AVec = cellC - cylinderCenter;
// tmp tensor for calculating the axial/radial components of cellC2AVec
tensor cellC2AVecE(tensor::zero);
cellC2AVecE.xx() = cellC2AVec.x();
cellC2AVecE.yy() = cellC2AVec.y();
cellC2AVecE.zz() = cellC2AVec.z();

// now we need to decompose cellC2AVec into axial and radial components
// the axial component of cellC2AVec vector
vector cellC2AVecA = cellC2AVecE & cylinderDirNorm;
// the radial component of cellC2AVec vector
vector cellC2AVecR = cellC2AVec - cellC2AVecA;

// the magnitude of radial component of cellC2AVecR
scalar cellC2AVecRLen = mag(cellC2AVecR);
// the magnitude of axial component of cellC2AVecR
scalar cellC2AVecALen = mag(cellC2AVecA);

if (cellC2AVecRLen <= radius && cellC2AVecALen <= cylinderLen / 2.0)
{
fvSource[cellI] += power / totalVolume;
}
}
}
else
{
FatalErrorIn("calcFvSourceCells") << "source: " << sourceType << " not supported!"
<< "Options are: cylinderAnnulusToCell!"
<< abort(FatalError);
FatalErrorIn("DAFvSourceHeatSource") << "source: " << sourceType << " not supported!"
<< "Options are: cylinderAnnulusToCell and cylinder!"
<< abort(FatalError);
}
}

if (daOption_.getOption<label>("debug"))
{
Info << "fvSourceCellIndices " << fvSourceCellIndices << endl;
}
fvSource.correctBoundaryConditions();
}

} // End namespace Foam
Expand Down
7 changes: 6 additions & 1 deletion src/adjoint/DAFvSource/DAFvSourceHeatSource.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,8 +31,13 @@ protected:
/// HashTable that contains lists of cell indices that are within the actuator disk space
HashTable<labelList> fvSourceCellIndices_;

HashTable<scalarList> cylinderP1_;
HashTable<scalarList> cylinderP2_;
HashTable<scalar> cylinderRadius_;
HashTable<scalar> power_;

/// calculate DAFvSourceHeatSource::fvSourceCellIndices_
void calcFvSourceCellIndices(HashTable<labelList>& fvSourceCellIndices);
void calcFvSourceCellIndices(const word sourceName, HashTable<labelList>& fvSourceCellIndices);

public:
TypeName("heatSource");
Expand Down
6 changes: 3 additions & 3 deletions tests/refs/DAFoam_Test_DAHeatTransferFoamRef.txt
View file
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
Dictionary Key: HF_INNER
@value 189999.9999999991 1e-08 1e-10
@value 181261.4809205142 1e-08 1e-10
Dictionary Key: HF_OUTER
@value -209999.9999999985 1e-08 1e-10
@value -218738.5190794858 1e-08 1e-10
Dictionary Key: TNormSum
@value 6400011.348643273 1e-08 1e-10
@value 6400011.350827799 1e-08 1e-10
Dictionary Key: fail
@value 0 1e-08 1e-10
17 changes: 12 additions & 5 deletions tests/runTests_DAHeatTransferFoam.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,11 +53,18 @@
"fvSource": {
"source1": {
"type": "heatSource",
"source": "cylinderAnnulusToCell",
"p1": [0.0, 0.0, 0.0],
"p2": [0.5, 0.0, 0.0],
"innerRadius": 0.0,
"outerRadius": 10.0,
"source": "cylinderToCell",
"p1": [0.6, 0.055, 0.025],
"p2": [1.0, 0.055, 0.025],
"radius": 0.01,
"power": 1000.0,
},
"source2": {
"type": "heatSource",
"source": "cylinder",
"p1": [0.0, 0.055, 0.025],
"p2": [0.4, 0.055, 0.025],
"radius": 0.005,
"power": 1000.0,
},
},
Expand Down

0 comments on commit b981121

Please sign in to comment.