pyDASolvers.pyx
Go to the documentation of this file.
1 
2 # distutils: language = c++
3 # distutils: sources = DASolvers.C
4 
5 """
6  DAFoam : Discrete Adjoint with OpenFOAM
7  Version : v4
8 
9  Description:
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
13 """
14 
15 # for using Petsc
16 from petsc4py.PETSc cimport Vec, PetscVec, Mat, PetscMat, KSP, PetscKSP
17 cimport numpy as np
18 np.import_array() # initialize C API to call PyArray_SimpleNewFromData
19 
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)
24 
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 *)
28 
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)
33 
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)
40 
41 cdef void pySetModelNameCallBack(const char* modelName, void *func):
42  (<object>func)(modelName)
43 
44 # declare cpp functions
45 cdef extern from "DASolvers.H" namespace "Foam":
46  cppclass DASolvers:
47  DASolvers(char *, object) except +
48  void initSolver()
49  int solvePrimal()
50  void runColoring()
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()
73  int getNLocalCells()
74  int getNLocalPoints()
75  int checkMesh()
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)
97  double getEndTime()
98  double getDeltaT()
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()
108 
109 # create python wrappers that call cpp functions
110 cdef class pyDASolvers:
111 
112  # define a class pointer for cpp functions
113  cdef:
114  DASolvers * _thisptr
115 
116  # initialize this class pointer with NULL
117  def __cinit__(self):
118  self._thisptr = NULL
119 
120  # deallocate the class pointer, and
121  # make sure we don't have memory leak
122  def __dealloc__(self):
123  if self._thisptr != NULL:
124  del self._thisptr
125 
126  # point the class pointer to the cpp class constructor
127  def __init__(self, argsAll, pyOptions):
128  """
129  Parameters
130  ----------
131 
132  argsAll: char
133  Chars that contains all the arguments
134  for running OpenFOAM solvers, including
135  the name of the solver.
136 
137  pyOptions: dict
138  Dictionary that defines all the options
139  in DAFoam
140 
141  Examples
142  --------
143  solver = pyDASolvers("DASolvers -parallel -python", aeroOptions)
144  """
145  self._thisptr = new DASolvers(argsAll, pyOptions)
146 
147  # wrap all the other member functions in the cpp class
148  def initSolver(self):
149  self._thisptr.initSolver()
150 
151  def solvePrimal(self):
152  return self._thisptr.solvePrimal()
153 
154  def runColoring(self):
155  self._thisptr.runColoring()
156 
157  def setSolverInput(self,
158  inputName,
159  inputType,
160  inputSize,
161  np.ndarray[double, ndim=1, mode="c"] inputs,
162  np.ndarray[double, ndim=1, mode="c"] seeds):
163 
164  assert len(inputs) == inputSize, "invalid input array size!"
165  assert len(seeds) == inputSize, "invalid seed array size!"
166 
167  cdef double *inputs_data = <double*>inputs.data
168  cdef double *seeds_data = <double*>seeds.data
169 
170  self._thisptr.setSolverInput(
171  inputName.encode(),
172  inputType.encode(),
173  inputSize,
174  inputs_data,
175  seeds_data)
176 
177  def getInputSize(self, inputName, inputType):
178  return self._thisptr.getInputSize(inputName.encode(), inputType.encode())
179 
180  def getOutputSize(self, outputName, outputType):
181  return self._thisptr.getOutputSize(outputName.encode(), outputType.encode())
182 
183  def getInputDistributed(self, inputName, inputType):
184  return self._thisptr.getInputDistributed(inputName.encode(), inputType.encode())
185 
186  def getOutputDistributed(self, outputName, outputType):
187  return self._thisptr.getOutputDistributed(outputName.encode(), outputType.encode())
188 
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)
192 
193  def calcJacTVecProduct(self,
194  inputName,
195  inputType,
196  np.ndarray[double, ndim=1, mode="c"] inputs,
197  outputName,
198  outputType,
199  np.ndarray[double, ndim=1, mode="c"] seeds,
200  np.ndarray[double, ndim=1, mode="c"] product):
201 
202  inputSize = self.getInputSize(inputName, inputType)
203  outputSize = self.getOutputSize(outputName, outputType)
204 
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!"
208 
209  cdef double *inputs_data = <double*>inputs.data
210  cdef double *seeds_data = <double*>seeds.data
211  cdef double *product_data = <double*>product.data
212 
213  self._thisptr.calcJacTVecProduct(
214  inputName.encode(),
215  inputType.encode(),
216  inputs_data,
217  outputName.encode(),
218  outputType.encode(),
219  seeds_data,
220  product_data)
221 
222  def calcdRdWT(self, isPC, Mat dRdWT):
223  self._thisptr.calcdRdWT(isPC, dRdWT.mat)
224 
225  def calcdRdWOldTPsiAD(self,
226  oldTimeLevel,
227  np.ndarray[double, ndim=1, mode="c"] psi,
228  np.ndarray[double, ndim=1, mode="c"] dRdWOldTPsi):
229 
230  assert len(psi) == self.getNLocalAdjointStates(), "invalid input array size!"
231  assert len(dRdWOldTPsi) == self.getNLocalAdjointStates(), "invalid seed array size!"
232 
233  cdef double *psi_data = <double*>psi.data
234  cdef double *dRdWOldTPsi_data = <double*>dRdWOldTPsi.data
235 
236  self._thisptr.calcdRdWOldTPsiAD(oldTimeLevel, psi_data, dRdWOldTPsi_data)
237 
238  def initializedRdWTMatrixFree(self):
239  self._thisptr.initializedRdWTMatrixFree()
240 
241  def destroydRdWTMatrixFree(self):
242  self._thisptr.destroydRdWTMatrixFree()
243 
244  def createMLRKSPMatrixFree(self, Mat jacPCMat, KSP myKSP):
245  self._thisptr.createMLRKSPMatrixFree(jacPCMat.mat, myKSP.ksp)
246 
247  def updateKSPPCMat(self, Mat PCMat, KSP myKSP):
248  self._thisptr.updateKSPPCMat(PCMat.mat, myKSP.ksp)
249 
250  def solveLinearEqn(self, KSP myKSP, Vec rhsVec, Vec solVec):
251  return self._thisptr.solveLinearEqn(myKSP.ksp, rhsVec.vec, solVec.vec)
252 
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)
257 
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)
262 
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)
267 
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)
275 
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)
280 
281  def getGlobalXvIndex(self, pointI, coordI):
282  return self._thisptr.getGlobalXvIndex(pointI, coordI)
283 
284  def getNLocalAdjointStates(self):
285  return self._thisptr.getNLocalAdjointStates()
286 
287  def getNLocalAdjointBoundaryStates(self):
288  return self._thisptr.getNLocalAdjointBoundaryStates()
289 
290  def getNLocalCells(self):
291  return self._thisptr.getNLocalCells()
292 
293  def getNLocalPoints(self):
294  return self._thisptr.getNLocalPoints()
295 
296  def checkMesh(self):
297  return self._thisptr.checkMesh()
298 
299  def getTimeOpFuncVal(self, functionName):
300  return self._thisptr.getTimeOpFuncVal(functionName.encode())
301 
302  def getdFScaling(self, functionName, timeIdx=-1):
303  return self._thisptr.getdFScaling(functionName.encode(), timeIdx)
304 
305  def getElapsedClockTime(self):
306  return self._thisptr.getElapsedClockTime()
307 
308  def getElapsedCpuTime(self):
309  return self._thisptr.getElapsedCpuTime()
310 
311  def calcCouplingFaceCoords(self,
312  np.ndarray[double, ndim=1, mode="c"] volCoords,
313  np.ndarray[double, ndim=1, mode="c"] surfCoords):
314 
315  assert len(volCoords) == self.getNLocalPoints() * 3, "invalid array size!"
316 
317  cdef double *volCoords_data = <double*>volCoords.data
318  cdef double *surfCoords_data = <double*>surfCoords.data
319 
320  self._thisptr.calcCouplingFaceCoords(volCoords_data, surfCoords_data)
321 
322  def getNRegressionParameters(self, modelName):
323  return self._thisptr.getNRegressionParameters(modelName)
324 
325  def printAllOptions(self):
326  self._thisptr.printAllOptions()
327 
328  def updateDAOption(self, pyOptions):
329  self._thisptr.updateDAOption(pyOptions)
330 
331  def getPrevPrimalSolTime(self):
332  return self._thisptr.getPrevPrimalSolTime()
333 
334  def writeFailedMesh(self):
335  self._thisptr.writeFailedMesh()
336 
337  def updateBoundaryConditions(self, fieldName, fieldType):
338  self._thisptr.updateBoundaryConditions(fieldName.encode(), fieldType.encode())
339 
340  def updateStateBoundaryConditions(self):
341  self._thisptr.updateStateBoundaryConditions()
342 
343  def calcPrimalResidualStatistics(self, mode):
344  self._thisptr.calcPrimalResidualStatistics(mode.encode())
345 
346  def setPrimalBoundaryConditions(self, printInfo):
347  self._thisptr.setPrimalBoundaryConditions(printInfo)
348 
349  def readStateVars(self, timeVal, timeLevel):
350  self._thisptr.readStateVars(timeVal, timeLevel)
351 
352  def readMeshPoints(self, timeVal):
353  self._thisptr.readMeshPoints(timeVal)
354 
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
358 
359  self._thisptr.writeMeshPoints(points_data, timeVal)
360 
361  def calcPCMatWithFvMatrix(self, Mat PCMat, turbOnly=0):
362  self._thisptr.calcPCMatWithFvMatrix(PCMat.mat, turbOnly)
363 
364  def setTime(self, time, timeIndex):
365  self._thisptr.setTime(time, timeIndex)
366 
367  def getDdtSchemeOrder(self):
368  return self._thisptr.getDdtSchemeOrder()
369 
370  def getEndTime(self):
371  return self._thisptr.getEndTime()
372 
373  def getDeltaT(self):
374  return self._thisptr.getDeltaT()
375 
376  def runFPAdj(self, Vec dFdW, Vec psi):
377  return self._thisptr.runFPAdj(dFdW.vec, psi.vec)
378 
379  def solveAdjointFP(self, Vec dFdW, Vec psi):
380  return self._thisptr.solveAdjointFP(dFdW.vec, psi.vec)
381 
382  def initTensorFlowFuncs(self, compute, jacVecProd, setModelName):
383  self._thisptr.initTensorFlowFuncs(pyCalcBetaCallBack, <void*>compute, pyCalcBetaJacVecProdCallBack, <void*>jacVecProd, pySetModelNameCallBack, <void*>setModelName)
384 
385  def writeSensMapSurface(self,
386  name,
387  np.ndarray[double, ndim=1, mode="c"] dFdXs,
388  np.ndarray[double, ndim=1, mode="c"] Xs,
389  size,
390  timeName):
391 
392  assert len(dFdXs) == size, "invalid array size!"
393  assert len(Xs) == size, "invalid array size!"
394 
395  cdef double *dFdXs_data = <double*>dFdXs.data
396  cdef double *Xs_data = <double*>Xs.data
397 
398  self._thisptr.writeSensMapSurface(
399  name.encode(),
400  dFdXs_data,
401  Xs_data,
402  size,
403  timeName)
404 
405  def writeSensMapField(self,
406  name,
407  np.ndarray[double, ndim=1, mode="c"] dFdField,
408  fieldType,
409  timeName):
410 
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!"
416  else:
417  print("fieldType can be either scalar or vector")
418  exit(1)
419 
420  cdef double *dFdField_data = <double*>dFdField.data
421 
422  self._thisptr.writeSensMapField(
423  name.encode(),
424  dFdField_data,
425  fieldType.encode(),
426  timeName)
427 
428  def getLatestTime(self):
429  return self._thisptr.getLatestTime()
430 
431  def hasVolCoordInput(self):
432  return self._thisptr.hasVolCoordInput()
433 
434  def writeAdjointFields(self, function, writeTime, np.ndarray[double, ndim=1, mode="c"] psi):
435  nAdjStates = self.getNLocalAdjointStates()
436 
437  assert len(psi) == nAdjStates, "invalid array size!"
438 
439  cdef double *psi_data = <double*>psi.data
440 
441  return self._thisptr.writeAdjointFields(function.encode(), writeTime, psi_data)
442 
443  def meanStatesToStates(self):
444  self._thisptr.meanStatesToStates()
445 
446  def updateInputFieldUnsteady(self):
447  self._thisptr.updateInputFieldUnsteady()