Skip to content

Commit

Permalink
Added a func to get OFField to the Python layer.
Browse files Browse the repository at this point in the history
  • Loading branch information
friedenhe committed Nov 24, 2023
1 parent e07872a commit 423a8a8
Show file tree
Hide file tree
Showing 6 changed files with 113 additions and 12 deletions.
41 changes: 41 additions & 0 deletions dafoam/pyDAFoam.py
Original file line number Diff line number Diff line change
Expand Up @@ -2962,6 +2962,47 @@ def getAcousticData(self, groupName=None):
# Finally map the vector as required.
return positions, normals, areas, forces

def getOFField(self, fieldName, fieldType, distributed=False):
"""
Get the OpenFOAM field variable in mesh.Db() and return it as a serial or parallel array
Input:
------
fieldName: name of the OF field
fieldType: type of the OF Field, can be either scalar or vector
distributed: whether the array is distributed (parallel)
Output:
------
field: numpy array that saves the field
"""

if fieldType == "scalar":
fieldComp = 1
elif fieldType == "vector":
fieldComp = 3
else:
raise Error("fieldType not valid!")
nLocalCells = self.solver.getNLocalCells()

vecSize = fieldComp * nLocalCells
fieldVec = PETSc.Vec().create(comm=PETSc.COMM_WORLD)
fieldVec.setSizes((vecSize, PETSc.DECIDE), bsize=1)
fieldVec.setFromOptions()
fieldVec.zeroEntries()

self.solver.getOFField(fieldName.encode(), fieldType.encode(), fieldVec)

if distributed:
field = self.vec2Array(fieldVec)
else:
field = self.convertMPIVec2SeqArray(fieldVec)

return field

def calcTotalDeriv(self, dRdX, dFdX, psi, totalDeriv):
"""
Compute total derivative
Expand Down
43 changes: 43 additions & 0 deletions src/adjoint/DASolver/DASolver.C
Original file line number Diff line number Diff line change
Expand Up @@ -1130,6 +1130,49 @@ void DASolver::getThermalAD(
#endif
}

void DASolver::getOFField(
const word fieldName,
const word fieldType,
Vec field) const
{
/*
Description:
assign a OpenFoam layer field variable in mesh.Db() to field
*/

PetscScalar* vecArray;
VecGetArray(field, &vecArray);

if (fieldType == "scalar")
{
const volScalarField& field = meshPtr_->thisDb().lookupObject<volScalarField>(fieldName);
forAll(field, cellI)
{
assignValueCheckAD(vecArray[cellI], field[cellI]);
}
}
else if (fieldType == "vector")
{
const volVectorField& field = meshPtr_->thisDb().lookupObject<volVectorField>(fieldName);
label localIdx = 0;
forAll(field, cellI)
{
for (label comp = 0; comp < 3; comp++)
{
assignValueCheckAD(vecArray[localIdx], field[cellI][comp]);
localIdx++;
}
}
}
else
{
FatalErrorIn("getField") << " fieldType not valid. Options: scalar or vector"
<< abort(FatalError);
}

VecRestoreArray(field, &vecArray);
}

void DASolver::getForces(Vec fX, Vec fY, Vec fZ)
{
/*
Expand Down
6 changes: 6 additions & 0 deletions src/adjoint/DASolver/DASolver.H
Original file line number Diff line number Diff line change
Expand Up @@ -689,6 +689,12 @@ public:
daFieldPtr_->resVec2OFResField(resVec);
}

/// get a variable from OF layer
void getOFField(
const word fieldName,
const word fieldType,
Vec field) const;

/// write the matrix in binary format
void writeMatrixBinary(
const Mat matIn,
Expand Down
9 changes: 9 additions & 0 deletions src/pyDASolvers/DASolvers.H
Original file line number Diff line number Diff line change
Expand Up @@ -576,6 +576,15 @@ public:
DASolverPtr_->calcdRdThermalTPsiAD(volCoords, states, thermal, seeds, product);
}

/// get a variable from OF layer
void getOFField(
const word fieldName,
const word fieldType,
Vec field)
{
DASolverPtr_->getOFField(fieldName, fieldType, field);
}

/// return the acoustic data values
void getAcousticData(Vec x, Vec y, Vec z, Vec nX, Vec nY, Vec nZ, Vec a, Vec fX, Vec fY, Vec fZ, word groupName)
{
Expand Down
4 changes: 4 additions & 0 deletions src/pyDASolvers/pyDASolvers.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -101,6 +101,7 @@ cdef extern from "DASolvers.H" namespace "Foam":
void getThermal(double *, double *, double *)
void getThermalAD(char *, double *, double *, double *, double *)
void setThermal(double *)
void getOFField(char *, char *, PetscVec)
void getAcousticData(PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, char*)
void printAllOptions()
void updateDAOption(object)
Expand Down Expand Up @@ -408,6 +409,9 @@ cdef class pyDASolvers:
def getForces(self, Vec fX, Vec fY, Vec fZ):
self._thisptr.getForces(fX.vec, fY.vec, fZ.vec)

def getOFField(self, fieldName, fieldType, Vec fieldVec):
self._thisptr.getOFField(fieldName, fieldType, fieldVec.vec)

def getThermal(self,
np.ndarray[double, ndim=1, mode="c"] volCoords,
np.ndarray[double, ndim=1, mode="c"] states,
Expand Down
22 changes: 10 additions & 12 deletions tests/runTests_DAPimpleFoamField.py
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,6 @@

twist0 = 30
U0 = 10
nCells = 4032
alpha0 = 0

# test incompressible solvers
Expand Down Expand Up @@ -85,17 +84,6 @@
DVGeo = DVGeometry("./FFD/FFD.xyz")
DVGeo.addRefAxis("bodyAxis", xFraction=0.25, alignIndex="k")


def betaFieldInversion(val, geo):
for idxI, v in enumerate(val):
DASolver.setFieldValue4GlobalCellI(b"betaFI", v, idxI)
DASolver.updateBoundaryConditions(b"betaFI", b"scalar")


beta0 = np.ones(nCells, dtype="d")
beta0[0] = 1.0
DVGeo.addGlobalDV("beta", value=beta0, func=betaFieldInversion, lower=1e-5, upper=10.0, scale=1.0)

# =============================================================================
# DAFoam initialization
# =============================================================================
Expand All @@ -107,6 +95,16 @@ def betaFieldInversion(val, geo):
evalFuncs = []
DASolver.setEvalFuncs(evalFuncs)


def betaFieldInversion(val, geo):
for idxI, v in enumerate(val):
DASolver.setFieldValue4GlobalCellI(b"betaFI", v, idxI)
DASolver.updateBoundaryConditions(b"betaFI", b"scalar")


beta0 = DASolver.getOFField("betaFI", "scalar", False)
DVGeo.addGlobalDV("beta", value=beta0, func=betaFieldInversion, lower=1e-5, upper=10.0, scale=1.0)

# =============================================================================
# Constraint setup
# =============================================================================
Expand Down

0 comments on commit 423a8a8

Please sign in to comment.