Skip to content

Commit

Permalink
Fixed a bug for the variance and added a util for extracting data. (#521
Browse files Browse the repository at this point in the history
)
  • Loading branch information
friedenhe authored Nov 20, 2023
1 parent f1495b9 commit 11caf24
Show file tree
Hide file tree
Showing 4 changed files with 159 additions and 5 deletions.
13 changes: 8 additions & 5 deletions src/adjoint/DAObjFunc/DAObjFuncVariance.C
Original file line number Diff line number Diff line change
Expand Up @@ -83,6 +83,7 @@ void DAObjFuncVariance::calcObjFunc(

// initialize objFunValue
objFuncValue = 0.0;
label nRefPoints = 0;

const objectRegistry& db = mesh_.thisDb();

Expand All @@ -99,16 +100,15 @@ void DAObjFuncVariance::calcObjFunc(
IOobject::NO_WRITE),
mesh_);

label nRefPoints = 0;
forAll(var, cellI)
{
if (varData[cellI] < 1e16)
{
objFuncValue += scale_ * (var[cellI] - varData[cellI]) * (var[cellI] - varData[cellI]);
scalar varDif = (var[cellI] - varData[cellI]);
objFuncValue += scale_ * varDif * varDif;
nRefPoints++;
}
}
objFuncValue /= nRefPoints;
}
else if (varType_ == "vector")
{
Expand All @@ -123,7 +123,6 @@ void DAObjFuncVariance::calcObjFunc(
IOobject::NO_WRITE),
mesh_);

label nRefPoints = 0;
forAll(var, cellI)
{
for (label comp = 0; comp < 3; comp++)
Expand All @@ -136,7 +135,6 @@ void DAObjFuncVariance::calcObjFunc(
}
}
}
objFuncValue /= nRefPoints;
}
else
{
Expand All @@ -145,9 +143,14 @@ void DAObjFuncVariance::calcObjFunc(
<< abort(FatalError);
}

// reduce the sum of all the ref points for averaging
reduce(nRefPoints, sumOp<label>());

// need to reduce the sum of force across all processors
reduce(objFuncValue, sumOp<scalar>());

objFuncValue /= nRefPoints;

return;
}

Expand Down
3 changes: 3 additions & 0 deletions src/utilities/postProcessing/getFIData/Make/files
Original file line number Diff line number Diff line change
@@ -0,0 +1,3 @@
getFIData.C

EXE = $(DAFOAM_ROOT_PATH)/OpenFOAM/sharedBins/getFIData
12 changes: 12 additions & 0 deletions src/utilities/postProcessing/getFIData/Make/options
Original file line number Diff line number Diff line change
@@ -0,0 +1,12 @@
EXE_INC = \
-std=c++11 \
-Wno-old-style-cast \
-Wno-conversion-null \
-Wno-deprecated-copy \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude

EXE_LIBS = \
-lfiniteVolume \
-lmeshTools

136 changes: 136 additions & 0 deletions src/utilities/postProcessing/getFIData/getFIData.C
Original file line number Diff line number Diff line change
@@ -0,0 +1,136 @@
/*---------------------------------------------------------------------------*\
DAFoam : Discrete Adjoint with OpenFOAM
Version : v3
Description
Extract the reference data for field inversion
\*---------------------------------------------------------------------------*/

#include "fvCFD.H"
#include "argList.H"
#include "Time.H"
#include "fvMesh.H"
#include "OFstream.H"

using namespace Foam;

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

int main(int argc, char* argv[])
{

Info << "Get ref data for field inversion...." << endl;

#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"

// Read FIDataDict
IOdictionary FIDataDict(
IOobject(
"FIDataDict",
mesh.time().system(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE));

scalar endTime = runTime.endTime().value();
scalar deltaT = runTime.deltaT().value();
label nInstances = round(endTime / deltaT);

forAll(FIDataDict.toc(), idxI)
{
word varName = FIDataDict.toc()[idxI];
dictionary dataSubDict = FIDataDict.subDict(varName);
word dataType = dataSubDict.getWord("dataType");
word fieldType = dataSubDict.getWord("fieldType");

if (dataType == "probePoints")
{
List<List<scalar>> probePointCoords = {{0, 0, 0}};
dataSubDict.readEntry<List<List<scalar>>>("pointCoords", probePointCoords);
// Info << "Reading probePointCoords from FIDataDict. probePointCoords = " << probePointCoords << endl;

// hard coded probe point info.
labelList probePointCellIList(probePointCoords.size(), -999);
forAll(probePointCoords, idxI)
{
point point = {probePointCoords[idxI][0], probePointCoords[idxI][1], probePointCoords[idxI][2]};
probePointCellIList[idxI] = mesh.findNearestCell(point);
}

for (label n = 1; n < nInstances + 1; n++)
{
scalar t = n * deltaT;
runTime.setTime(t, n);

if (fieldType == "scalar")
{
volScalarField varRead(
IOobject(
varName,
mesh.time().timeName(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE),
mesh);

forAll(varRead, cellI)
{
if (!probePointCellIList.found(cellI))
{
varRead[cellI] = 1e17;
}
}

varRead.rename(varName + "Data");
varRead.write();
}
else if (fieldType == "vector")
{
volVectorField varRead(
IOobject(
varName,
mesh.time().timeName(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE),
mesh);

forAll(varRead, cellI)
{
if (!probePointCellIList.found(cellI))
{
for (label i = 0; i < 3; i++)
{
varRead[cellI][i] = 1e17;
}
}
}

varRead.rename(varName + "Data");
varRead.write();
}
else
{
FatalErrorIn("")
<< "fieldType" << fieldType << " not supported! Options are: scalar or vector"
<< abort(FatalError);
}
}
}
else
{
FatalErrorIn("")
<< "dataType" << dataType << " not supported! Options are: probePoints"
<< abort(FatalError);
}
}
Info << "Finished!" << endl;

return 0;
}

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

0 comments on commit 11caf24

Please sign in to comment.