Coverage for /opt/hostedtoolcache/Python/3.11.10/x64/lib/python3.11/site-packages/pysagas/cfd/solver.py: 31%
143 statements
« prev ^ index » next coverage.py v7.6.4, created at 2024-10-30 04:27 +0000
1from __future__ import annotations
2import os
3import copy
4import meshio
5import numpy as np
6import pandas as pd
7from pysagas import banner
8from abc import ABC, abstractmethod
9from typing import List, Optional, Dict
10from pysagas import Cell, FlowState, Vector
13class AbstractFlowSolver(ABC):
14 """Abstract interface for flow solvers."""
16 method = None
18 @abstractmethod
19 def __init__(self, **kwargs) -> None:
20 pass
22 def __repr__(self) -> str:
23 return f"PySAGAS {self.method} flow solver"
25 def __str__(self) -> str:
26 return f"PySAGAS {self.method} flow solver"
28 @property
29 @abstractmethod
30 def method(self):
31 # This is a placeholder for a class variable defining the solver method
32 pass
34 @abstractmethod
35 def solve(self, *args, **kwargs):
36 pass
39class FlowSolver(AbstractFlowSolver):
40 """Base class for CFD flow solver."""
42 def __init__(
43 self,
44 cells: List[Cell],
45 freestream: Optional[FlowState] = None,
46 verbosity: Optional[int] = 1,
47 ) -> None:
48 """Instantiates the flow solver.
50 Parameters
51 ----------
52 cells : list[Cell]
53 A list of the cells describing the geometry.
55 freestream : Flowstate, optional
56 The free-stream flow state. The default is the freestream provided
57 upon instantiation of the solver.
59 verbosity : int, optional
60 The verbosity of the solver. The default is 1.
61 """
63 # TODO - allow providing a geometry parser, eg. tri file parser for Cart3D,
64 # or a (yet to be implemented) STL parser.
66 self.cells = cells
67 self.freestream = freestream
68 self.verbosity = verbosity
69 self._last_solve_freestream: FlowState = None
70 self._last_sens_freestream: FlowState = None
72 # Results
73 self.flow_result: FlowResults = None
74 self.flow_sensitivity: SensitivityResults = None
76 # Print banner
77 if self.verbosity > 0:
78 banner()
79 print(f"\033[4m{self.__repr__()}\033[0m".center(50, " "))
81 def solve(
82 self,
83 freestream: Optional[FlowState] = None,
84 mach: Optional[float] = None,
85 aoa: Optional[float] = None,
86 ) -> bool:
87 """Run the flow solver.
89 Parameters
90 ----------
91 freestream : Flowstate, optional
92 The free-stream flow state. The default is the freestream provided
93 upon instantiation of the solver.
95 mach : float, optional
96 The free-stream Mach number. The default is that specified in
97 the freestream flow state.
99 aoa : float, optional
100 The free-stream angle of attack. The default is that specified in
101 the freestream flow state.
103 Returns
104 -------
105 result : FlowResults
106 The flow results.
108 Raises
109 ------
110 Exception : when no freestream can be found.
111 """
113 if not freestream:
114 # No freestream conditions specified
115 if not self.freestream:
116 # No conditions provided on instantiation either
117 raise Exception("Please specify freestream conditions.")
118 else:
119 # Use nominal freestream as base
120 freestream = copy.copy(self.freestream)
121 if mach:
122 # Update mach number
123 freestream._M = mach
125 if aoa:
126 # Update flow direction
127 flow_direction = Vector(1, 1 * np.tan(np.deg2rad(aoa)), 0).unit
128 freestream.direction = flow_direction
130 else:
131 # Solve-specific freestream conditions provided
132 if mach or aoa:
133 print("Using freestream provided; ignoring Mach/aoa provided.")
135 # Check if already solved
136 if self._last_solve_freestream and freestream == self._last_solve_freestream:
137 return True
138 else:
139 # Update last solve freestream and continue
140 self._last_solve_freestream = freestream
141 return False
143 def solve_sens(
144 self,
145 freestream: Optional[FlowState] = None,
146 mach: Optional[float] = None,
147 aoa: Optional[float] = None,
148 ) -> SensitivityResults:
149 """Run the flow solver to obtain sensitivity information.
151 Parameters
152 ----------
153 freestream : Flowstate, optional
154 The free-stream flow state. The default is the freestream provided
155 upon instantiation of the solver.
157 Mach : float, optional
158 The free-stream Mach number. The default is that specified in
159 the freestream flow state.
161 aoa : float, optional
162 The free-stream angle of attack. The default is that specified in
163 the freestream flow state.
165 Returns
166 -------
167 SensitivityResults
169 Raises
170 ------
171 Exception : when no freestream can be found.
172 """
173 if not freestream:
174 # No freestream conditions specified
175 if not self.freestream:
176 # No conditions provided on instantiation either
177 raise Exception("Please specify freestream conditions.")
178 else:
179 # Use nominal freestream as base
180 freestream = copy.copy(self.freestream)
181 if mach:
182 # Update mach number
183 freestream._M = mach
185 if aoa:
186 # Update flow direction
187 flow_direction = Vector(1, 1 * np.tan(np.deg2rad(aoa)), 0).unit
188 freestream.direction = flow_direction
190 else:
191 # Solve-specific freestream conditions provided
192 if mach or aoa:
193 print("Using freestream provided; ignoring Mach/aoa provided.")
195 # Check if already solved
196 if self._last_sens_freestream and freestream == self._last_sens_freestream:
197 return True
198 else:
199 # Update last solve freestream and continue
200 self._last_sens_freestream = freestream
201 return False
203 def save(self, name: str, attributes: Dict[str, list]):
204 """Save the solution to VTK file format. Note that the
205 geometry mesh must have been loaded from a parser which
206 includes connectivity information (eg. PyMesh or TRI).
208 Parameters
209 ----------
210 name : str
211 The filename prefix of the desired output file.
213 attributes : Dict[str, list]
214 The attributes dictionary used to initialise the file writer.
215 Note that this is not required when calling `solve` from a
216 child FlowSolver.
217 """
218 # Construct vertices and faces
219 faces = []
220 vertex_faces = {}
221 for cell in self.cells:
222 # Get faces
223 faces.append(cell._face_ids)
225 # Get attributes
226 for a in attributes:
227 attributes[a].append(cell.attributes.get(a))
229 # Get vertices
230 for i, vid in enumerate(cell._face_ids):
231 vertex_faces[vid] = cell.vertices[i]
233 # Sort vertices
234 sorted_vertices = dict(sorted(vertex_faces.items()))
236 # Note: vertex attributes can be treated the same way as
237 # the vertices - they must be sorted.
239 # Finally
240 vertices = np.array([v for v in sorted_vertices.values()])
241 faces = np.array(faces)
243 # Convert mesh object with attributes
244 mesh_obj = meshio.Mesh(
245 points=vertices,
246 cells=[("triangle", faces)],
247 cell_data={k: [v] for k, v in attributes.items()},
248 )
249 mesh_obj.write(f"{name}.vtk")
251 @staticmethod
252 def body_to_wind(v: Vector, aoa: float):
253 """Converts a vector from body axes to wind axes.
255 Parameters
256 ----------
257 v : Vector
258 The result to transform.
260 aoa : float
261 The angle of attack, specified in degrees.
263 Returns
264 --------
265 r : Vector
266 The transformed result.
267 """
268 # TODO - use euler matrices
269 A = v.x # axial
270 N = v.y # normal
272 L = N * np.cos(np.deg2rad(aoa)) - A * np.sin(np.deg2rad(aoa))
273 D = N * np.sin(np.deg2rad(aoa)) + A * np.cos(np.deg2rad(aoa))
275 # Construct new result
276 r = Vector(x=D, y=L, z=v.z)
278 return r
280 def sweep(self, aoa_range: list[float], mach_range: list[float]) -> pd.DataFrame:
281 """Evaluate the aerodynamic coefficients over a sweep of angle
282 of attack and mach numbers.
284 Parameters
285 -----------
286 aoa_range : list[float]
287 The angles of attack to evaluate at, specified in degrees.
289 mach_range : list[float]
290 The Mach numbers to evaluate at.
292 Returns
293 --------
294 results : pd.DataFrame
295 A Pandas DataFrame of the aerodynamic coefficients at each angle of
296 attack and Mach number combination.
297 """
298 # TODO - add bool to evaluate gradients on the sweep
299 # Run sweep
300 i = 0
301 results = pd.DataFrame()
302 for aoa in aoa_range:
303 for mach in mach_range:
304 flowresult: FlowResults = self.solve(aoa=aoa, Mach=mach)
305 CL, CD, Cm = flowresult.coefficients()
307 coefficients_df = pd.DataFrame(
308 data={"aoa": aoa, "mach": mach, "CL": CL, "CD": CD, "Cm": Cm},
309 index=[i],
310 )
312 # Append results
313 results = pd.concat(
314 [
315 results,
316 coefficients_df,
317 ]
318 )
320 # Iterate index
321 i += 1
323 return results
326class FlowResults:
327 """A class containing the aerodynamic force and moment
328 information with respect to design parameters.
330 Attributes
331 ----------
332 freestream : FlowState
333 The freestream flow state.
335 net_force : pd.DataFrame
336 The net force in cartesian coordinate frame (x,y,z).
338 m_sense : pd.DataFrame
339 The net moment in cartesian coordinate frame (x,y,z).
340 """
342 def __init__(
343 self, freestream: FlowState, net_force: Vector, net_moment: Vector
344 ) -> None:
345 self.freestream = freestream
346 self.net_force = net_force
347 self.net_moment = net_moment
349 # Calculate angle of attack
350 self.aoa = np.rad2deg(
351 np.arctan(freestream.direction.y / freestream.direction.x)
352 )
354 def __str__(self) -> str:
355 return f"Net force = {self.net_force} N"
357 def __repr__(self) -> str:
358 return self.__str__()
360 def coefficients(self, A_ref: Optional[float] = 1.0, c_ref: Optional[float] = 1.0):
361 """Calculate the aerodynamic coefficients CL, CD and Cm.
363 Parameters
364 -----------
365 A_ref : float, optional
366 The reference area (m^2). The default is 1 m^2.
368 c_ref : float, optional
369 The reference length (m). The default is 1 m.
371 Returns
372 -------
373 C_L, C_D, C_m
374 """
375 w = FlowSolver.body_to_wind(v=self.net_force, aoa=self.aoa)
376 C_L = w.y / (self.freestream.q * A_ref)
377 C_D = w.x / (self.freestream.q * A_ref)
379 mw = FlowSolver.body_to_wind(v=self.net_moment, aoa=self.aoa)
380 C_m = -mw.z / (self.freestream.q * A_ref * c_ref)
381 return C_L, C_D, C_m
383 @property
384 def coef(self):
385 CL, CD, Cm = self.coefficients()
386 print(f"CL = {CL:.3f}\nCD = {CD:.3f}\nCm = {Cm:.3f}")
389class SensitivityResults:
390 """A class containing the aerodynamic force and moment sensitivity
391 information with respect to design parameters.
393 Attributes
394 ----------
395 freestream : FlowState
396 The freestream flow state.
398 f_sense : pd.DataFrame
399 The force sensitivities in cartesian coordinate frame (x,y,z).
401 m_sense : pd.DataFrame
402 The moment sensitivities in cartesian coordinate frame (x,y,z).
403 """
405 def __init__(
406 self,
407 f_sens: pd.DataFrame,
408 m_sens: pd.DataFrame,
409 freestream: Optional[FlowState] = None,
410 ) -> None:
411 self.f_sens = f_sens
412 self.m_sens = m_sens
413 self.freestream = freestream
415 def __str__(self) -> str:
416 s1 = f"Force sensitivties:\n{self.f_sens}"
417 s2 = f"Moment sensitivties:\n{self.m_sens}"
419 return f"{s1}\n\n{s2}"
421 def __repr__(self) -> str:
422 return self.__str__()
424 def coefficients(
425 self, A_ref: Optional[float] = 1.0, c_ref: Optional[float] = 1.0
426 ) -> tuple[pd.DataFrame, pd.DataFrame]:
427 """Calculate the aerodynamic coefficients CL, CD and Cm.
429 Parameters
430 -----------
431 A_ref : float, optional
432 The reference area (m^2). The default is 1 m^2.
434 c_ref : float, optional
435 The reference length (m). The default is 1 m.
437 Returns
438 -------
439 A tuple (cf_sens, cm_sens) containing the force and moment coefficient sensitivities.
440 """
441 # Non-dimensionalise
442 cf_sens: pd.DataFrame = self.f_sens / (self.freestream.q * A_ref)
443 cm_sens: pd.DataFrame = self.m_sens / (self.freestream.q * A_ref * c_ref)
445 # Translate to aero frame
446 # TODO - review if this might swap value order on overwrite of columns
447 aoa = self.freestream.aoa
448 cf_sens["dFy/dp"] = cf_sens["dFy/dp"] * np.cos(np.deg2rad(aoa)) - cf_sens[
449 "dFx/dp"
450 ] * np.sin(np.deg2rad(aoa))
451 cf_sens["dFx/dp"] = cf_sens["dFy/dp"] * np.sin(np.deg2rad(aoa)) + cf_sens[
452 "dFx/dp"
453 ] * np.cos(np.deg2rad(aoa))
455 # Rename columns
456 cf_sens.rename(
457 columns={"dFx/dp": "dCD", "dFy/dp": "dCL", "dFz/dp": "dCZ"}, inplace=True
458 )
459 # cm_sens.rename(columns={"dMx/dp": "dCl", "dMy/dp": "dCn", "dMz/dp": "dCm"}, inplace=True)
461 return cf_sens, cm_sens