2 # distutils: language = c++
3 # distutils: sources = DASolvers.C
6 DAFoam : Discrete Adjoint with OpenFOAM
10 Cython wrapper functions that call OpenFOAM libraries defined
11 in the *.C and *.H files. The python naming convention is to
12 add "py" before the C++ class name
16 from petsc4py.PETSc cimport Vec, PetscVec, Mat, PetscMat, KSP, PetscKSP
18 np.import_array() # initialize C API to call PyArray_SimpleNewFromData
20 cdef public api CPointerToPyArray(const double* data, int size) with gil:
21 if not (data and size >= 0): raise ValueError
22 cdef np.npy_intp dims = size
23 return np.PyArray_SimpleNewFromData(1, &dims, np.NPY_DOUBLE, <void*>data)
25 ctypedef void (*pyComputeInterface)(const double *, int, double *, int, void *)
26 ctypedef void (*pyJacVecProdInterface)(const double *, double *, int, const double *, const double *, int, void *)
27 ctypedef void (*pySetCharInterface)(const char *, void *)
29 cdef void pyCalcBetaCallBack(const double* inputs, int n, double* outputs, int m, void *func):
30 inputs_data = CPointerToPyArray(inputs, n)
31 outputs_data = CPointerToPyArray(outputs, m)
32 (<object>func)(inputs_data, n, outputs_data, m)
34 cdef void pyCalcBetaJacVecProdCallBack(const double* inputs, double* inputs_b, int n, const double* outputs, const double* outputs_b, int m, void *func):
35 inputs_data = CPointerToPyArray(inputs, n)
36 inputs_b_data = CPointerToPyArray(inputs_b, n)
37 outputs_data = CPointerToPyArray(outputs, m)
38 outputs_b_data = CPointerToPyArray(outputs_b, m)
39 (<object>func)(inputs_data, inputs_b_data, n, outputs_data, outputs_b_data, m)
41 cdef void pySetModelNameCallBack(const char* modelName, void *func):
42 (<object>func)(modelName)
44 # declare cpp functions
45 cdef extern from "DASolvers.H" namespace "Foam":
47 DASolvers(char *, object) except +
51 void calcJacTVecProduct(char *, char *, double *, char *, char *, double *, double *)
52 int getInputSize(char *, char *)
53 int getOutputSize(char *, char *)
54 void calcOutput(char *, char *, double *)
55 int getInputDistributed(char *, char *)
56 int getOutputDistributed(char *, char *)
57 void setSolverInput(char *, char *, int, double *, double *)
58 void calcdRdWT(int, PetscMat)
59 void initializedRdWTMatrixFree()
60 void destroydRdWTMatrixFree()
61 void createMLRKSPMatrixFree(PetscMat, PetscKSP)
62 void updateKSPPCMat(PetscMat, PetscKSP)
63 int solveLinearEqn(PetscKSP, PetscVec, PetscVec)
64 void calcdRdWOldTPsiAD(int, double *, double *)
65 void updateOFFields(double *)
66 void getOFFields(double *)
67 void getOFField(char *, char *, double *)
68 void getOFMeshPoints(double *)
69 void updateOFMesh(double *)
70 int getGlobalXvIndex(int, int)
71 int getNLocalAdjointStates()
72 int getNLocalAdjointBoundaryStates()
76 double getdFScaling(char *, int)
77 double getTimeOpFuncVal(char *)
78 double getElapsedClockTime()
79 double getElapsedCpuTime()
80 void calcCouplingFaceCoords(double *, double *)
81 int getNRegressionParameters(char *)
82 void printAllOptions()
83 void updateDAOption(object)
84 double getPrevPrimalSolTime()
85 void writeFailedMesh()
86 void updateBoundaryConditions(char *, char *)
87 void updateStateBoundaryConditions()
88 void calcPrimalResidualStatistics(char *)
89 void setPrimalBoundaryConditions(int)
90 int runFPAdj(PetscVec, PetscVec)
91 int solveAdjointFP(PetscVec, PetscVec)
92 void initTensorFlowFuncs(pyComputeInterface, void *, pyJacVecProdInterface, void *, pySetCharInterface, void *)
93 void readStateVars(double, int)
94 void readMeshPoints(double)
95 void writeMeshPoints(double *, double)
96 void calcPCMatWithFvMatrix(PetscMat, int)
99 void setTime(double, int)
100 int getDdtSchemeOrder()
101 void writeSensMapSurface(char *, double *, double *, int, double)
102 void writeSensMapField(char *, double *, char *, double)
103 double getLatestTime()
104 void writeAdjointFields(char *, double, double *)
105 int hasVolCoordInput()
106 void meanStatesToStates()
107 void updateInputFieldUnsteady()
109 # create python wrappers that call cpp functions
110 cdef class pyDASolvers:
112 # define a class pointer for cpp functions
116 # initialize this class pointer with NULL
120 # deallocate the class pointer, and
121 # make sure we don't have memory leak
122 def __dealloc__(self):
123 if self._thisptr != NULL:
126 # point the class pointer to the cpp class constructor
127 def __init__(self, argsAll, pyOptions):
133 Chars that contains all the arguments
134 for running OpenFOAM solvers, including
135 the name of the solver.
138 Dictionary that defines all the options
143 solver = pyDASolvers("DASolvers -parallel -python", aeroOptions)
145 self._thisptr = new DASolvers(argsAll, pyOptions)
147 # wrap all the other member functions in the cpp class
148 def initSolver(self):
149 self._thisptr.initSolver()
151 def solvePrimal(self):
152 return self._thisptr.solvePrimal()
154 def runColoring(self):
155 self._thisptr.runColoring()
157 def setSolverInput(self,
161 np.ndarray[double, ndim=1, mode="c"] inputs,
162 np.ndarray[double, ndim=1, mode="c"] seeds):
164 assert len(inputs) == inputSize, "invalid input array size!"
165 assert len(seeds) == inputSize, "invalid seed array size!"
167 cdef double *inputs_data = <double*>inputs.data
168 cdef double *seeds_data = <double*>seeds.data
170 self._thisptr.setSolverInput(
177 def getInputSize(self, inputName, inputType):
178 return self._thisptr.getInputSize(inputName.encode(), inputType.encode())
180 def getOutputSize(self, outputName, outputType):
181 return self._thisptr.getOutputSize(outputName.encode(), outputType.encode())
183 def getInputDistributed(self, inputName, inputType):
184 return self._thisptr.getInputDistributed(inputName.encode(), inputType.encode())
186 def getOutputDistributed(self, outputName, outputType):
187 return self._thisptr.getOutputDistributed(outputName.encode(), outputType.encode())
189 def calcOutput(self, outputName, outputType, np.ndarray[double, ndim=1, mode="c"] output):
190 cdef double *output_data = <double*>output.data
191 self._thisptr.calcOutput(outputName.encode(), outputType.encode(), output_data)
193 def calcJacTVecProduct(self,
196 np.ndarray[double, ndim=1, mode="c"] inputs,
199 np.ndarray[double, ndim=1, mode="c"] seeds,
200 np.ndarray[double, ndim=1, mode="c"] product):
202 inputSize = self.getInputSize(inputName, inputType)
203 outputSize = self.getOutputSize(outputName, outputType)
205 assert len(inputs) == inputSize, "invalid input array size!"
206 assert len(seeds) == outputSize, "invalid seed array size!"
207 assert len(product) == inputSize, "invalid product array size!"
209 cdef double *inputs_data = <double*>inputs.data
210 cdef double *seeds_data = <double*>seeds.data
211 cdef double *product_data = <double*>product.data
213 self._thisptr.calcJacTVecProduct(
222 def calcdRdWT(self, isPC, Mat dRdWT):
223 self._thisptr.calcdRdWT(isPC, dRdWT.mat)
225 def calcdRdWOldTPsiAD(self,
227 np.ndarray[double, ndim=1, mode="c"] psi,
228 np.ndarray[double, ndim=1, mode="c"] dRdWOldTPsi):
230 assert len(psi) == self.getNLocalAdjointStates(), "invalid input array size!"
231 assert len(dRdWOldTPsi) == self.getNLocalAdjointStates(), "invalid seed array size!"
233 cdef double *psi_data = <double*>psi.data
234 cdef double *dRdWOldTPsi_data = <double*>dRdWOldTPsi.data
236 self._thisptr.calcdRdWOldTPsiAD(oldTimeLevel, psi_data, dRdWOldTPsi_data)
238 def initializedRdWTMatrixFree(self):
239 self._thisptr.initializedRdWTMatrixFree()
241 def destroydRdWTMatrixFree(self):
242 self._thisptr.destroydRdWTMatrixFree()
244 def createMLRKSPMatrixFree(self, Mat jacPCMat, KSP myKSP):
245 self._thisptr.createMLRKSPMatrixFree(jacPCMat.mat, myKSP.ksp)
247 def updateKSPPCMat(self, Mat PCMat, KSP myKSP):
248 self._thisptr.updateKSPPCMat(PCMat.mat, myKSP.ksp)
250 def solveLinearEqn(self, KSP myKSP, Vec rhsVec, Vec solVec):
251 return self._thisptr.solveLinearEqn(myKSP.ksp, rhsVec.vec, solVec.vec)
253 def updateOFFields(self, np.ndarray[double, ndim=1, mode="c"] states):
254 assert len(states) == self.getNLocalAdjointStates(), "invalid array size!"
255 cdef double *states_data = <double*>states.data
256 self._thisptr.updateOFFields(states_data)
258 def getOFFields(self, np.ndarray[double, ndim=1, mode="c"] states):
259 assert len(states) == self.getNLocalAdjointStates(), "invalid array size!"
260 cdef double *states_data = <double*>states.data
261 self._thisptr.getOFFields(states_data)
263 def getOFMeshPoints(self, np.ndarray[double, ndim=1, mode="c"] points):
264 assert len(points) == self.getNLocalPoints() * 3, "invalid array size!"
265 cdef double *points_data = <double*>points.data
266 self._thisptr.getOFMeshPoints(points_data)
268 def getOFField(self, fieldName, fieldType, np.ndarray[double, ndim=1, mode="c"] field):
269 if fieldType == "scalar":
270 assert len(field) == self.getNLocalCells(), "invalid array size!"
271 elif fieldType == "vector":
272 assert len(field) == self.getNLocalCells() * 3, "invalid array size!"
273 cdef double *field_data = <double*>field.data
274 self._thisptr.getOFField(fieldName.encode(), fieldType.encode(), field_data)
276 def updateOFMesh(self, np.ndarray[double, ndim=1, mode="c"] vol_coords):
277 assert len(vol_coords) == self.getNLocalPoints() * 3, "invalid array size!"
278 cdef double *vol_coords_data = <double*>vol_coords.data
279 self._thisptr.updateOFMesh(vol_coords_data)
281 def getGlobalXvIndex(self, pointI, coordI):
282 return self._thisptr.getGlobalXvIndex(pointI, coordI)
284 def getNLocalAdjointStates(self):
285 return self._thisptr.getNLocalAdjointStates()
287 def getNLocalAdjointBoundaryStates(self):
288 return self._thisptr.getNLocalAdjointBoundaryStates()
290 def getNLocalCells(self):
291 return self._thisptr.getNLocalCells()
293 def getNLocalPoints(self):
294 return self._thisptr.getNLocalPoints()
297 return self._thisptr.checkMesh()
299 def getTimeOpFuncVal(self, functionName):
300 return self._thisptr.getTimeOpFuncVal(functionName.encode())
302 def getdFScaling(self, functionName, timeIdx=-1):
303 return self._thisptr.getdFScaling(functionName.encode(), timeIdx)
305 def getElapsedClockTime(self):
306 return self._thisptr.getElapsedClockTime()
308 def getElapsedCpuTime(self):
309 return self._thisptr.getElapsedCpuTime()
311 def calcCouplingFaceCoords(self,
312 np.ndarray[double, ndim=1, mode="c"] volCoords,
313 np.ndarray[double, ndim=1, mode="c"] surfCoords):
315 assert len(volCoords) == self.getNLocalPoints() * 3, "invalid array size!"
317 cdef double *volCoords_data = <double*>volCoords.data
318 cdef double *surfCoords_data = <double*>surfCoords.data
320 self._thisptr.calcCouplingFaceCoords(volCoords_data, surfCoords_data)
322 def getNRegressionParameters(self, modelName):
323 return self._thisptr.getNRegressionParameters(modelName)
325 def printAllOptions(self):
326 self._thisptr.printAllOptions()
328 def updateDAOption(self, pyOptions):
329 self._thisptr.updateDAOption(pyOptions)
331 def getPrevPrimalSolTime(self):
332 return self._thisptr.getPrevPrimalSolTime()
334 def writeFailedMesh(self):
335 self._thisptr.writeFailedMesh()
337 def updateBoundaryConditions(self, fieldName, fieldType):
338 self._thisptr.updateBoundaryConditions(fieldName.encode(), fieldType.encode())
340 def updateStateBoundaryConditions(self):
341 self._thisptr.updateStateBoundaryConditions()
343 def calcPrimalResidualStatistics(self, mode):
344 self._thisptr.calcPrimalResidualStatistics(mode.encode())
346 def setPrimalBoundaryConditions(self, printInfo):
347 self._thisptr.setPrimalBoundaryConditions(printInfo)
349 def readStateVars(self, timeVal, timeLevel):
350 self._thisptr.readStateVars(timeVal, timeLevel)
352 def readMeshPoints(self, timeVal):
353 self._thisptr.readMeshPoints(timeVal)
355 def writeMeshPoints(self, np.ndarray[double, ndim=1, mode="c"] points, timeVal):
356 assert len(points) == self.getNLocalPoints() * 3, "invalid array size!"
357 cdef double *points_data = <double*>points.data
359 self._thisptr.writeMeshPoints(points_data, timeVal)
361 def calcPCMatWithFvMatrix(self, Mat PCMat, turbOnly=0):
362 self._thisptr.calcPCMatWithFvMatrix(PCMat.mat, turbOnly)
364 def setTime(self, time, timeIndex):
365 self._thisptr.setTime(time, timeIndex)
367 def getDdtSchemeOrder(self):
368 return self._thisptr.getDdtSchemeOrder()
370 def getEndTime(self):
371 return self._thisptr.getEndTime()
374 return self._thisptr.getDeltaT()
376 def runFPAdj(self, Vec dFdW, Vec psi):
377 return self._thisptr.runFPAdj(dFdW.vec, psi.vec)
379 def solveAdjointFP(self, Vec dFdW, Vec psi):
380 return self._thisptr.solveAdjointFP(dFdW.vec, psi.vec)
382 def initTensorFlowFuncs(self, compute, jacVecProd, setModelName):
383 self._thisptr.initTensorFlowFuncs(pyCalcBetaCallBack, <void*>compute, pyCalcBetaJacVecProdCallBack, <void*>jacVecProd, pySetModelNameCallBack, <void*>setModelName)
385 def writeSensMapSurface(self,
387 np.ndarray[double, ndim=1, mode="c"] dFdXs,
388 np.ndarray[double, ndim=1, mode="c"] Xs,
392 assert len(dFdXs) == size, "invalid array size!"
393 assert len(Xs) == size, "invalid array size!"
395 cdef double *dFdXs_data = <double*>dFdXs.data
396 cdef double *Xs_data = <double*>Xs.data
398 self._thisptr.writeSensMapSurface(
405 def writeSensMapField(self,
407 np.ndarray[double, ndim=1, mode="c"] dFdField,
411 nCells = self.getNLocalCells()
412 if fieldType == "scalar":
413 assert len(dFdField) == nCells, "invalid array size!"
414 elif fieldType == "vector":
415 assert len(dFdField) == 3 * nCells, "invalid array size!"
417 print("fieldType can be either scalar or vector")
420 cdef double *dFdField_data = <double*>dFdField.data
422 self._thisptr.writeSensMapField(
428 def getLatestTime(self):
429 return self._thisptr.getLatestTime()
431 def hasVolCoordInput(self):
432 return self._thisptr.hasVolCoordInput()
434 def writeAdjointFields(self, function, writeTime, np.ndarray[double, ndim=1, mode="c"] psi):
435 nAdjStates = self.getNLocalAdjointStates()
437 assert len(psi) == nAdjStates, "invalid array size!"
439 cdef double *psi_data = <double*>psi.data
441 return self._thisptr.writeAdjointFields(function.encode(), writeTime, psi_data)
443 def meanStatesToStates(self):
444 self._thisptr.meanStatesToStates()
446 def updateInputFieldUnsteady(self):
447 self._thisptr.updateInputFieldUnsteady()