Skip to content

Commit

Permalink
Fixing bug in fvSource Calculation (#527)
Browse files Browse the repository at this point in the history 
* Adding volume mesh coordinates to fvSource functions

* Updating regression test references
  • Loading branch information
@bernardopacini
bernardopacini authored Nov 23, 2023
1 parent fbbdcd3 commit 1098eef
Show file tree
Hide file tree
Showing 7 changed files with 244 additions and 175 deletions.
40 changes: 30 additions & 10 deletions dafoam/mphys/mphys_dafoam.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1977,6 +1977,8 @@ def setup(self):

self.DASolver = self.options["solver"]

self.discipline = self.DASolver.getOption("discipline")

# loop over all the propNames and check if any of them is active. If yes, add inputs for this prop
for propName in list(self.DASolver.getOption("wingProp").keys()):
if self.DASolver.getOption("wingProp")[propName]["active"]:
Expand All @@ -1999,6 +2001,8 @@ def setup(self):
self.add_input(propName + "_integral_force", distributed=False, shape=2, tags=["mphys_coupling"])
self.add_input(propName + "_prop_center", distributed=False, shape=3, tags=["mphys_coupling"])

self.add_input("%s_vol_coords" % self.discipline, distributed=True, shape_by_conn=True, tags=["mphys_coupling"])

# we have only one output
self.nLocalCells = self.DASolver.solver.getNLocalCells()
self.add_output("fvSource", distributed=True, shape=self.nLocalCells * 3, tags=["mphys_coupling"])
Expand All @@ -2007,6 +2011,9 @@ def compute(self, inputs, outputs):

DASolver = self.DASolver

dafoam_xv = inputs["%s_vol_coords" % self.discipline]
xvVec = DASolver.array2Vec(dafoam_xv)

# initialize output to zeros
fvSourceVec = PETSc.Vec().create(self.comm)
fvSourceVec.setSizes((self.nLocalCells * 3, PETSc.DECIDE), bsize=1)
Expand Down Expand Up @@ -2040,6 +2047,7 @@ def compute(self, inputs, outputs):
radial_location_vec,
integral_force_vec,
prop_center_vec,
xvVec,
fvSourceVec,
)

Expand All @@ -2058,6 +2066,9 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode):
)
return

dafoam_xv = inputs["%s_vol_coords" % self.discipline]
xvVec = DASolver.array2Vec(dafoam_xv)

if "fvSource" in d_outputs:
sBar = d_outputs["fvSource"]
sBarVec = DASolver.array2Vec(sBar)
Expand All @@ -2081,52 +2092,61 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode):
prodVec = PETSc.Vec().createSeq(self.nForceSections, bsize=1, comm=PETSc.COMM_SELF)
prodVec.zeroEntries()
DASolver.solverAD.calcdFvSourcedInputsTPsiAD(
propName.encode(), "aForce".encode(), aVec, tVec, rVec, fVec, cVec, sBarVec, prodVec
propName.encode(), "aForce".encode(), aVec, tVec, rVec, fVec, cVec, xvVec, sBarVec, prodVec
)
aBar = DASolver.vec2ArraySeq(prodVec)
# aBar = self.comm.allreduce(aBar, op=MPI.SUM)
aBar = self.comm.allreduce(aBar, op=MPI.SUM)
d_inputs[propName + "_axial_force"] += aBar

if propName + "_tangential_force" in d_inputs:
prodVec = PETSc.Vec().createSeq(self.nForceSections, bsize=1, comm=PETSc.COMM_SELF)
prodVec.zeroEntries()
DASolver.solverAD.calcdFvSourcedInputsTPsiAD(
propName.encode(), "tForce".encode(), aVec, tVec, rVec, fVec, cVec, sBarVec, prodVec
propName.encode(), "tForce".encode(), aVec, tVec, rVec, fVec, cVec, xvVec, sBarVec, prodVec
)
tBar = DASolver.vec2ArraySeq(prodVec)
# tBar = self.comm.allreduce(tBar, op=MPI.SUM)
tBar = self.comm.allreduce(tBar, op=MPI.SUM)
d_inputs[propName + "_tangential_force"] += tBar

if propName + "_radial_location" in d_inputs:
prodVec = PETSc.Vec().createSeq(self.nForceSections, bsize=1, comm=PETSc.COMM_SELF)
prodVec.zeroEntries()
DASolver.solverAD.calcdFvSourcedInputsTPsiAD(
propName.encode(), "rDist".encode(), aVec, tVec, rVec, fVec, cVec, sBarVec, prodVec
propName.encode(), "rDist".encode(), aVec, tVec, rVec, fVec, cVec, xvVec, sBarVec, prodVec
)
rBar = DASolver.vec2ArraySeq(prodVec)
# rBar = self.comm.allreduce(rBar, op=MPI.SUM)
rBar = self.comm.allreduce(rBar, op=MPI.SUM)
d_inputs[propName + "_radial_location"] += rBar

if propName + "_integral_force" in d_inputs:
prodVec = PETSc.Vec().createSeq(2, bsize=1, comm=PETSc.COMM_SELF)
prodVec.zeroEntries()
DASolver.solverAD.calcdFvSourcedInputsTPsiAD(
propName.encode(), "targetForce".encode(), aVec, tVec, rVec, fVec, cVec, sBarVec, prodVec
propName.encode(), "targetForce".encode(), aVec, tVec, rVec, fVec, cVec, xvVec, sBarVec, prodVec
)
fBar = DASolver.vec2ArraySeq(prodVec)
# fBar = self.comm.allreduce(fBar, op=MPI.SUM)
fBar = self.comm.allreduce(fBar, op=MPI.SUM)
d_inputs[propName + "_integral_force"] += fBar

if propName + "_prop_center" in d_inputs:
prodVec = PETSc.Vec().createSeq(3, bsize=1, comm=PETSc.COMM_SELF)
prodVec.zeroEntries()
DASolver.solverAD.calcdFvSourcedInputsTPsiAD(
propName.encode(), "center".encode(), aVec, tVec, rVec, fVec, cVec, sBarVec, prodVec
propName.encode(), "center".encode(), aVec, tVec, rVec, fVec, cVec, xvVec, sBarVec, prodVec
)
cBar = DASolver.vec2ArraySeq(prodVec)
# cBar = self.comm.allreduce(cBar, op=MPI.SUM)
cBar = self.comm.allreduce(cBar, op=MPI.SUM)
d_inputs[propName + "_prop_center"] += cBar

if "%s_vol_coords" % self.discipline in d_inputs:
prodVec = xvVec.duplicate()
prodVec.zeroEntries()
DASolver.solverAD.calcdFvSourcedInputsTPsiAD(
propName.encode(), "mesh".encode(), aVec, tVec, rVec, fVec, cVec, xvVec, sBarVec, prodVec
)
vBar = DASolver.vec2Array(prodVec)
# vBar = self.comm.allreduce(vBar, op=MPI.SUM)
d_inputs["%s_vol_coords" % self.discipline] += vBar

class DAFoamPropNodes(ExplicitComponent):
"""
Expand Down
45 changes: 45 additions & 0 deletions src/adjoint/DASolver/DASolver.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2462,6 +2462,7 @@ void DASolver::calcFvSource(
Vec rDist,
Vec targetForce,
Vec center,
Vec xvVec,
Vec fvSource)
{
/*
Expand Down Expand Up @@ -2525,6 +2526,8 @@ void DASolver::calcFvSource(
centerTemp[1] = vecArrayCenter[1];
centerTemp[2] = vecArrayCenter[2];

this->updateOFMesh(xvVec);

// Compute fvSource
this->calcFvSourceInternal(propName, aForceTemp, tForceTemp, rDistTemp, targetForceTemp, centerTemp, fvSourceTemp);

Expand Down Expand Up @@ -2565,6 +2568,7 @@ void DASolver::calcdFvSourcedInputsTPsiAD(
Vec rDist,
Vec targetForce,
Vec center,
Vec xvVec,
Vec psi,
Vec dFvSource)
{
Expand All @@ -2580,6 +2584,8 @@ void DASolver::calcdFvSourcedInputsTPsiAD(

VecZeroEntries(dFvSource);

this->updateOFMesh(xvVec);

const dictionary& propSubDict = daOptionPtr_->getAllOptions().subDict("wingProp").subDict(propName);
label nPoints = propSubDict.getLabel("nForceSections");

Expand All @@ -2588,6 +2594,7 @@ void DASolver::calcdFvSourcedInputsTPsiAD(
Field<scalar> rDistField(nPoints);
List<scalar> targetForceList(2);
vector centerVector = vector::zero;
pointField meshPoints = meshPtr_->points();

volVectorField fvSourceVField(
IOobject(
Expand Down Expand Up @@ -2677,6 +2684,18 @@ void DASolver::calcdFvSourcedInputsTPsiAD(
this->globalADTape_.registerInput(centerVector[i]);
}
}
else if (mode == "mesh")
{
forAll(meshPoints, i)
{
for (label j = 0; j < 3; j++)
{
this->globalADTape_.registerInput(meshPoints[i][j]);
}
}
meshPtr_->movePoints(meshPoints);
meshPtr_->moving(false);
}
else
{
FatalErrorIn("calcdFvSourcedInputsTPsiAD") << "mode not valid"
Expand Down Expand Up @@ -2746,6 +2765,20 @@ void DASolver::calcdFvSourcedInputsTPsiAD(
vecArrayProd[i] = centerVector[i].getGradient();
}
}
else if (mode == "mesh")
{
forAll(meshPoints, i)
{
for (label j = 0; j < 3; j++)
{
label rowI = daIndexPtr_->getGlobalXvIndex(i, j);
PetscScalar val = meshPoints[i][j].getGradient();
VecSetValue(dFvSource, rowI, val, INSERT_VALUES);
}
}
VecAssemblyBegin(dFvSource);
VecAssemblyEnd(dFvSource);
}

VecRestoreArray(dFvSource, &vecArrayProd);

Expand Down Expand Up @@ -2791,6 +2824,18 @@ void DASolver::calcdFvSourcedInputsTPsiAD(
this->globalADTape_.deactivateValue(centerVector[i]);
}
}
else if (mode == "mesh")
{
forAll(meshPoints, i)
{
for (label j = 0; j < 3; j++)
{
this->globalADTape_.deactivateValue(meshPoints[i][j]);
}
}
meshPtr_->movePoints(meshPoints);
meshPtr_->moving(false);
}

this->calcFvSourceInternal(propName, aForceField, tForceField, rDistField, targetForceList, centerVector, fvSourceVField);
#endif
Expand Down
2 changes: 2 additions & 0 deletions src/adjoint/DASolver/DASolver.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -961,6 +961,7 @@ public:
Vec rDist,
Vec targetForce,
Vec center,
Vec xvVec,
Vec fvSource);

void calcdFvSourcedInputsTPsiAD(
Expand All @@ -971,6 +972,7 @@ public:
Vec rDist,
Vec targetForce,
Vec center,
Vec xvVec,
Vec psi,
Vec dFvSource);

Expand Down
6 changes: 4 additions & 2 deletions src/pyDASolvers/DASolvers.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -838,9 +838,10 @@ public:
Vec rDist,
Vec targetForce,
Vec center,
Vec xvVec,
Vec fvSource)
{
DASolverPtr_->calcFvSource(propName, aForce, tForce, rDist, targetForce, center, fvSource);
DASolverPtr_->calcFvSource(propName, aForce, tForce, rDist, targetForce, center, xvVec, fvSource);
}

void calcdFvSourcedInputsTPsiAD(
Expand All @@ -851,10 +852,11 @@ public:
Vec rDist,
Vec targetForce,
Vec center,
Vec xvVec,
Vec psi,
Vec dFvSource)
{
DASolverPtr_->calcdFvSourcedInputsTPsiAD(propName, mode, aForce, tForce, rDist, targetForce, center, psi, dFvSource);
DASolverPtr_->calcdFvSourcedInputsTPsiAD(propName, mode, aForce, tForce, rDist, targetForce, center, xvVec, psi, dFvSource);
}

void calcForceProfile(
Expand Down
12 changes: 6 additions & 6 deletions src/pyDASolvers/pyDASolvers.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -126,8 +126,8 @@ cdef extern from "DASolvers.H" namespace "Foam":
double getForwardADDerivVal(char *)
void calcResidualVec(PetscVec)
void setPrimalBoundaryConditions(int)
void calcFvSource(char *, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec)
void calcdFvSourcedInputsTPsiAD(char *, char *, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec)
void calcFvSource(char *, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec)
void calcdFvSourcedInputsTPsiAD(char *, char *, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec)
void calcForceProfile(char *, PetscVec, PetscVec, PetscVec, PetscVec)
void calcdForceProfiledXvWAD(char *, char *, char *, PetscVec, PetscVec, PetscVec, PetscVec)
void calcdForcedStateTPsiAD(char *, PetscVec, PetscVec, PetscVec, PetscVec)
Expand Down Expand Up @@ -547,11 +547,11 @@ cdef class pyDASolvers:
def getUnsteadyObjFuncEndTimeIndex(self):
return self._thisptr.getUnsteadyObjFuncEndTimeIndex()

def calcFvSource(self, propName, Vec aForce, Vec tForce, Vec rDist, Vec targetForce, Vec center, Vec fvSource):
self._thisptr.calcFvSource(propName, aForce.vec, tForce.vec, rDist.vec, targetForce.vec, center.vec, fvSource.vec)
def calcFvSource(self, propName, Vec aForce, Vec tForce, Vec rDist, Vec targetForce, Vec center, Vec xvVec, Vec fvSource):
self._thisptr.calcFvSource(propName, aForce.vec, tForce.vec, rDist.vec, targetForce.vec, center.vec, xvVec.vec, fvSource.vec)

def calcdFvSourcedInputsTPsiAD(self, propName, mode, Vec aForce, Vec tForce, Vec rDist, Vec targetForce, Vec center, Vec psi, Vec dFvSource):
self._thisptr.calcdFvSourcedInputsTPsiAD(propName, mode, aForce.vec, tForce.vec, rDist.vec, targetForce.vec, center.vec, psi.vec, dFvSource.vec)
def calcdFvSourcedInputsTPsiAD(self, propName, mode, Vec aForce, Vec tForce, Vec rDist, Vec targetForce, Vec center, Vec xvVec, Vec psi, Vec dFvSource):
self._thisptr.calcdFvSourcedInputsTPsiAD(propName, mode, aForce.vec, tForce.vec, rDist.vec, targetForce.vec, center.vec, xvVec.vec, psi.vec, dFvSource.vec)

def calcForceProfile(self, propName, Vec aForce, Vec tForce, Vec rDist, Vec integralForce):
self._thisptr.calcForceProfile(propName, aForce.vec, tForce.vec, rDist.vec, integralForce.vec)
Expand Down
Loading

0 comments on commit 1098eef

Please sign in to comment.