Skip to content

Commit

Permalink
Allowed heat source to move with derivs (mdolab#656)
Browse files Browse the repository at this point in the history 
* Added heat source parameters as DVs and location snapCenter option.

* Fixed some bugs and added the derivs for the heatSource pars.

* Fixed a bug.
  • Loading branch information
@friedenhe
friedenhe authored Jun 27, 2024
1 parent aa2d351 commit a941eaf
Show file tree
Hide file tree
Showing 20 changed files with 774 additions and 112 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
96 changes: 57 additions & 39 deletions src/adjoint/DAFvSource/DAFvSource.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -135,49 +135,67 @@ void DAFvSource::syncDAOptionToActuatorDVs()

// now we need to initialize actuatorDiskDVs_
dictionary fvSourceSubDict = daOption_.getAllOptions().subDict("fvSource");
word diskName0 = fvSourceSubDict.toc()[0];

word type0 = fvSourceSubDict.subDict(diskName0).getWord("type");

if (type0 == "actuatorDisk")
forAll(fvSourceSubDict.toc(), idxI)
{
word source0 = fvSourceSubDict.subDict(diskName0).getWord("source");
word diskName = fvSourceSubDict.toc()[idxI];
// sub dictionary with all parameters for this disk
dictionary diskSubDict = fvSourceSubDict.subDict(diskName);
word type = diskSubDict.getWord("type");
word source = diskSubDict.getWord("source");

if (type == "actuatorDisk" && source == "cylinderAnnulusSmooth")
{

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

scalarList dirList;
diskSubDict.readEntry<scalarList>("direction", dirList);

// we have 13 design variables for each disk
scalarList dvList(13);
dvList[0] = centerList[0];
dvList[1] = centerList[1];
dvList[2] = centerList[2];
dvList[3] = dirList[0];
dvList[4] = dirList[1];
dvList[5] = dirList[2];
dvList[6] = diskSubDict.getScalar("innerRadius");
dvList[7] = diskSubDict.getScalar("outerRadius");
dvList[8] = diskSubDict.getScalar("scale");
dvList[9] = diskSubDict.getScalar("POD");
dvList[10] = diskSubDict.getScalar("expM");
dvList[11] = diskSubDict.getScalar("expN");
dvList[12] = diskSubDict.getScalar("targetThrust");

// set actuatorDiskDVs_
actuatorDiskDVs_.set(diskName, dvList);
}
else if (type == "heatSource" && source == "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 dirList;
diskSubDict.readEntry<scalarList>("direction", dirList);

// we have 13 design variables for each disk
scalarList dvList(13);
dvList[0] = centerList[0];
dvList[1] = centerList[1];
dvList[2] = centerList[2];
dvList[3] = dirList[0];
dvList[4] = dirList[1];
dvList[5] = dirList[2];
dvList[6] = diskSubDict.getScalar("innerRadius");
dvList[7] = diskSubDict.getScalar("outerRadius");
dvList[8] = diskSubDict.getScalar("scale");
dvList[9] = diskSubDict.getScalar("POD");
dvList[10] = diskSubDict.getScalar("expM");
dvList[11] = diskSubDict.getScalar("expN");
dvList[12] = diskSubDict.getScalar("targetThrust");

// set actuatorDiskDVs_
actuatorDiskDVs_.set(diskName, dvList);
}
// 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
22 changes: 19 additions & 3 deletions src/adjoint/DAFvSource/DAFvSource.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -128,9 +128,25 @@ public:
/// calculate fvSource based on the latest actuatorDVs
void updateFvSource()
{
volVectorField& fvSource = const_cast<volVectorField&>(
mesh_.thisDb().lookupObject<volVectorField>("fvSource"));
this->calcFvSource(fvSource);
if (mesh_.thisDb().foundObject<volVectorField>("fvSource"))
{
volVectorField& fvSource = const_cast<volVectorField&>(
mesh_.thisDb().lookupObject<volVectorField>("fvSource"));

this->calcFvSource(fvSource);
}
else if (mesh_.thisDb().foundObject<volScalarField>("fvSource"))
{
volScalarField& fvSource = const_cast<volScalarField&>(
mesh_.thisDb().lookupObject<volScalarField>("fvSource"));

this->calcFvSource(fvSource);
}
else
{
FatalErrorIn("") << "fvSource not found! "
<< abort(FatalError);
}
}

/// synchronize the values in DAOption and actuatorDiskDVs_
Expand Down
38 changes: 13 additions & 25 deletions src/adjoint/DAFvSource/DAFvSourceActuatorDisk.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -103,21 +103,15 @@ void DAFvSourceActuatorDisk::calcFvSource(volVectorField& fvSource)

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

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

if (source0 == "cylinderAnnulusToCell")
forAll(fvSourceSubDict.toc(), idxI)
{
word diskName = fvSourceSubDict.toc()[idxI];
dictionary diskSubDict = fvSourceSubDict.subDict(diskName);
word source = diskSubDict.getWord("source");

// loop over all the cell indices for all actuator disks
forAll(fvSourceCellIndices_.toc(), idxI)
if (source == "cylinderAnnulusToCell")
{

// name of this disk
word diskName = fvSourceCellIndices_.toc()[idxI];

// sub dictionary with all parameters for this disk
dictionary diskSubDict = fvSourceSubDict.subDict(diskName);
// loop over all the cell indices for all actuator disks

// now read in all parameters for this actuator disk
scalarList point1;
Expand Down Expand Up @@ -213,14 +207,8 @@ void DAFvSourceActuatorDisk::calcFvSource(volVectorField& fvSource)
}
}
}
}
else if (source0 == "cylinderAnnulusSmooth")
{

forAll(fvSourceSubDict.toc(), idxI)
else if (source == "cylinderAnnulusSmooth")
{
word diskName = fvSourceSubDict.toc()[idxI];
dictionary diskSubDict = fvSourceSubDict.subDict(diskName);

vector center = {actuatorDiskDVs_[diskName][0], actuatorDiskDVs_[diskName][1], actuatorDiskDVs_[diskName][2]};
vector dirNorm = {actuatorDiskDVs_[diskName][3], actuatorDiskDVs_[diskName][4], actuatorDiskDVs_[diskName][5]};
Expand Down Expand Up @@ -432,12 +420,12 @@ void DAFvSourceActuatorDisk::calcFvSource(volVectorField& fvSource)
}
}
}
}
else
{
FatalErrorIn("calcFvSourceCells") << "source: " << source0 << " not supported!"
<< "Options are: cylinderAnnulusToCell and cylinderAnnulusSmooth!"
<< abort(FatalError);
else
{
FatalErrorIn("calcFvSourceCells") << "source: " << source << " not supported!"
<< "Options are: cylinderAnnulusToCell and cylinderAnnulusSmooth!"
<< abort(FatalError);
}
}

fvSource.correctBoundaryConditions();
Expand Down
Loading

0 comments on commit a941eaf

Please sign in to comment.