Coverage for /opt/hostedtoolcache/Python/3.11.10/x64/lib/python3.11/site-packages/pysagas/sensitivity/calculator.py: 82%
105 statements
« prev ^ index » next coverage.py v7.6.4, created at 2024-10-30 04:27 +0000
1import numpy as np
2import pandas as pd
3from pysagas.flow import FlowState
4from abc import ABC, abstractmethod
5from pysagas.geometry import Vector, Cell
6from pysagas.utilities import add_sens_data
7from pysagas.cfd.solver import SensitivityResults
8from typing import List, Callable, Tuple, Optional, Union, Literal
9from pysagas.sensitivity.models import (
10 piston_sensitivity,
11 van_dyke_sensitivity,
12 isentropic_sensitivity,
13)
16class AbstractSensitivityCalculator(ABC):
17 """Abstract sensitivity calculator class defining the interface."""
19 solver = None
21 @abstractmethod
22 def __init__(self, **kwargs) -> None:
23 pass
25 def __repr__(self) -> str:
26 return f"PySAGAS Sensitivity Calculator for {self.solver}"
28 def __str__(self) -> str:
29 return f"PySAGAS Sensitivity Calculator for {self.solver}"
31 @property
32 @abstractmethod
33 def solver(self):
34 # This is a placeholder for a class variable defining the Sensitivity Calculator's solver
35 pass
37 @abstractmethod
38 def _transcribe_cells(self, parameters: List):
39 """Transcribes the unified Cells from the solver data.
41 Returns
42 --------
43 cells : List[Cell]
44 A list of the Cells from the flow solution.
45 """
46 # This method should be overloaded with the solver-specific method
47 pass
49 @abstractmethod
50 def _extract_parameters(self) -> List[str]:
51 """Extract the geometry parameters.
53 Returns
54 -------
55 parameters : List[str]
56 A list of the parameter names.
57 """
58 # This method should be overloaded with the solver-specific method
59 pass
61 @abstractmethod
62 def to_csv(self):
63 """Dumps the sensitivity data to CSV file."""
64 pass
67class SensitivityCalculator(AbstractSensitivityCalculator):
68 """SensitivityCalculator base class."""
70 def __init__(self, **kwargs) -> None:
71 self.cells: list[Cell] = None
73 def calculate(
74 self,
75 sensitivity_model: Optional[
76 Union[Callable, Literal["piston", "van_dyke", "isentropic"]]
77 ] = "van_dyke",
78 cog: Optional[Vector] = Vector(0, 0, 0),
79 flowstate: Optional[FlowState] = None,
80 **kwargs,
81 ) -> SensitivityResults:
82 """Calculate the force sensitivities of the surface to the
83 parameters.
85 Parameters
86 ----------
87 sensitivity_model : Callable | Literal["piston", "van_dyke", "isentropic"], optional
88 The model used to calculate the pressure/parameter sensitivities.
89 The default is Van Dyke's second-order theory model.
91 cog : Vector, optional
92 The centre of gravity. The default is Vector(0, 0, 0).
94 flowstate : FlowState, optional
95 The flowstate associated with this result. The default is None.
97 Returns
98 -------
99 SensitivityResults : the sensitivity results object.
100 """
101 if self.verbosity > 0:
102 print("\nCalculating aerodynamic sensitivities.")
104 parameters = self._extract_parameters()
106 if self.cells is None:
107 self.cells = self._transcribe_cells(parameters=parameters)
109 params_sens_cols = []
110 for p in parameters:
111 for d in ["x", "y", "z"]:
112 params_sens_cols.append(f"d{d}d_{p}")
114 # Calculate force sensitivity
115 F_sense, M_sense = self.net_sensitivity(
116 cells=self.cells, sensitivity_model=sensitivity_model, cog=cog, **kwargs
117 )
119 # Construct dataframes to return
120 df_f = pd.DataFrame(
121 F_sense, columns=["dFx/dp", "dFy/dp", "dFz/dp"], index=parameters
122 )
123 df_m = pd.DataFrame(
124 M_sense, columns=["dMx/dp", "dMy/dp", "dMz/dp"], index=parameters
125 )
127 if self.verbosity > 0:
128 print("Done.")
130 # Construct results
131 result = SensitivityResults(
132 f_sens=df_f,
133 m_sens=df_m,
134 freestream=flowstate,
135 )
137 return result
139 def to_csv(self):
140 """Dumps the sensitivity data to CSV file."""
141 if self.cells is not None:
142 pass
144 else:
145 raise Exception("No cells have been transcribed.")
147 def _extract_parameters(self):
148 parameters = []
149 for e in self.sensdata.columns:
150 e: str
151 if e.startswith("dxd") or e.startswith("dyd") or e.startswith("dzd"):
152 # Sensitivity coluns
153 if e[3:] not in parameters:
154 parameters.append(e[3:])
155 return parameters
157 @staticmethod
158 def net_sensitivity(
159 cells: List[Cell],
160 sensitivity_model: Optional[
161 Union[Callable, Literal["piston", "van_dyke", "isentropic"]]
162 ] = "van_dyke",
163 cog: Vector = Vector(0, 0, 0),
164 **kwargs,
165 ) -> Tuple[np.array, np.array]:
166 """Calcualtes the net force and moment sensitivities for a list of Cells.
168 Parameters
169 ----------
170 cells : list[Cell]
171 The cells to be analysed.
173 sensitivity_model : Callable | Literal["piston", "van_dyke", "isentropic"], optional
174 The model used to calculate the pressure/parameter sensitivities.
175 The default is Van Dyke's second-order theory model.
177 cog : Vector, optional
178 The reference centre of gravity, used in calculating the moment
179 sensitivities. The defualt is Vector(0,0,0).
181 Returns
182 --------
183 dFdp : np.array
184 The force sensitivity matrix with respect to the parameters.
186 dMdp : np.array
187 The moment sensitivity matrix with respect to the parameters.
188 """
189 # Get sensitivity handle
190 if isinstance(sensitivity_model, Callable):
191 # Use function provided
192 sensitivity_function = sensitivity_model
193 elif isinstance(sensitivity_model, str):
194 # Get function from method specified
195 sensitivity_function = globals().get(f"{sensitivity_model}_sensitivity")
196 if sensitivity_function is None:
197 raise Exception("Invalid sensitivity method specified.")
199 dFdp = 0
200 dMdp = 0
201 for cell in cells:
202 # Calculate force sensitivity
203 dFdp_c, dMdp_c = SensitivityCalculator.cell_sensitivity(
204 cell=cell, sensitivity_function=sensitivity_function, cog=cog, **kwargs
205 )
207 dFdp += dFdp_c
208 dMdp += dMdp_c
210 return dFdp, dMdp
212 @staticmethod
213 def cell_sensitivity(
214 cell: Cell,
215 sensitivity_function: Callable,
216 cog: Vector = Vector(0, 0, 0),
217 **kwargs,
218 ) -> Tuple[np.array, np.array]:
219 """Calculates force and moment sensitivities for a single cell.
221 Parameters
222 ----------
223 cell : Cell
224 The cell.
226 sensitivity_function : Callable
227 The function to use when calculating the pressure sensitivities.
229 cog : Vector, optional
230 The reference centre of gravity, used in calculating the moment
231 sensitivities. The defualt is Vector(0,0,0).
233 Returns
234 --------
235 sensitivities : np.array
236 An array of shape n x 3, for a 3-dimensional cell with
237 n parameters.
239 See Also
240 --------
241 all_dfdp : a wrapper to calculate force sensitivities for many cells
242 """
243 # Initialisation
244 all_directions = [Vector(1, 0, 0), Vector(0, 1, 0), Vector(0, 0, 1)]
245 sensitivities = np.zeros(shape=(cell.dndp.shape[1], 3))
246 moment_sensitivities = np.zeros(shape=(cell.dndp.shape[1], 3))
248 # Calculate moment dependencies
249 r = cell.c - cog
250 F = cell.flowstate.P * cell.A * cell.n.vec
252 # For each parameter
253 for p_i in range(cell.dndp.shape[1]):
254 # Calculate pressure sensitivity
255 dPdp = sensitivity_function(cell=cell, p_i=p_i, **kwargs)
257 # Evaluate for sensitivities for each direction
258 for i, direction in enumerate(all_directions):
259 dF = (
260 dPdp * cell.A * np.dot(cell.n.vec, direction.vec)
261 + cell.flowstate.P
262 * cell.dAdp[p_i]
263 * np.dot(cell.n.vec, direction.vec)
264 + cell.flowstate.P
265 * cell.A
266 * np.dot(-cell.dndp[:, p_i], direction.vec)
267 )
268 sensitivities[p_i, i] = dF
270 # Now evaluate moment sensitivities
271 moment_sensitivities[p_i, :] = np.cross(
272 r.vec, sensitivities[p_i, :]
273 ) + np.cross(cell.dcdp[:, p_i], F)
275 # Append to cell
276 cell.sensitivities = sensitivities
277 cell.moment_sensitivities = moment_sensitivities
279 return sensitivities, moment_sensitivities
282class GenericSensitivityCalculator(SensitivityCalculator):
283 solver = "Generic Flow Solver"
285 def __init__(
286 self,
287 cells: List[Cell],
288 sensitivity_filepath: Optional[str] = None,
289 cells_have_sens_data: Optional[bool] = False,
290 verbosity: Optional[int] = 1,
291 **kwargs,
292 ) -> None:
293 """Instantiate a generic PySAGAS sensitivity calculator.
295 Parameters
296 ----------
297 cells : list[Cell]
298 A list of transcribed Cell objects, containing the nominal flow solution.
300 sensitivity_filepath : str
301 The filepath to the geometry sensitivities. This must be provided, if the
302 cells do not have geometric sensitivity information attached. This must
303 also be provided anyway, to determine the geometric design parameters.
305 cells_have_sens_data : bool, optional
306 The cells already have geometric sensitivity data matched to them. When this
307 is True, the sensitivity_filepath argument (if provided) is ignored. The
308 default is False.
310 verbosity : int, optional
311 The verbosity of the code. The defualt is 1.
312 """
313 super().__init__(**kwargs)
315 if cells_have_sens_data:
316 self.cells = cells
317 else:
318 self._pre_transcribed_cells = cells
320 self.sensdata = pd.read_csv(sensitivity_filepath)
321 self.verbosity = verbosity
323 def _transcribe_cells(self, **kwargs) -> List[Cell]:
324 """This is a dummy method to satisfy the abstract base class. Transcribed cells
325 are provided upon instantiation of the sensitivity calculator.
326 """
327 # Add sensitivity data to _pre_transcribed_cells
328 add_sens_data(
329 cells=self._pre_transcribed_cells,
330 data=self.sensdata,
331 verbosity=self.verbosity,
332 )
334 return self._pre_transcribed_cells