Coverage for /opt/hostedtoolcache/Python/3.11.10/x64/lib/python3.11/site-packages/pysagas/cfd/oblique_prandtl_meyer.py: 36%
258 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 tqdm import tqdm
4from scipy.optimize import bisect
5from typing import Optional, Tuple
6from pysagas.flow import FlowState
7from pysagas.geometry import Vector
8from scipy.optimize import root_scalar
9from pysagas.utilities import add_sens_data
10from pysagas.cfd.solver import FlowSolver, FlowResults, SensitivityResults
13class OPM(FlowSolver):
14 """Oblique shock thoery and Prandtl-Meyer expansion theory flow solver.
16 This implementation uses oblique shock theory for flow-facing cell
17 elements, and Prandtl-Meyer expansion theory for rearward-facing elements.
19 Extended Summary
20 ----------------
21 Data attribute 'method' refers to which method was used for a particular
22 cell, according to:
23 - -1 : invalid / skipped (eg. 90 degree face)
24 - 0 : parallel face, do nothing
25 - 1 : Prandlt-Meyer
26 - 2 : normal shock
27 - 3 : oblique shock
28 """
30 method = "Oblique/Prandtl-Meyer combination"
31 PM_ANGLE_THRESHOLD = -20 # degrees
33 def solve(
34 self,
35 freestream: Optional[FlowState] = None,
36 mach: Optional[float] = None,
37 aoa: Optional[float] = None,
38 cog: Vector = Vector(0, 0, 0),
39 ) -> FlowResults:
40 already_run = super().solve(freestream=freestream, mach=mach, aoa=aoa)
41 if already_run:
42 # Already have a result
43 if self.verbosity > 1:
44 print("Attention: this flowstate has already been solved.")
45 result = self.flow_result
47 else:
48 # Get flow
49 flow = self._last_solve_freestream
51 # Construct progress bar
52 if self.verbosity > 0:
53 print()
54 desc = f"Running OPM solver at AoA = {flow.aoa:.2f} and Mach = {flow.M:.2f}"
55 pbar = tqdm(
56 total=len(self.cells),
57 desc=desc,
58 position=0,
59 leave=True,
60 )
62 # Iterate over all cells
63 net_force = Vector(0, 0, 0)
64 net_moment = Vector(0, 0, 0)
65 bad = 0
66 total = 0
67 for cell in self.cells:
68 # Calculate orientation of cell to flow
69 theta = np.pi / 2 - np.arccos(
70 np.dot(flow.direction.vec, cell.n.vec)
71 / (cell.n.norm * flow.direction.norm)
72 )
74 # Also calculate vector to COG from cell centroid
75 r = cell.c - cog
77 # Solve flow for this cell
78 if theta < 0:
79 # Rearward facing cell
80 if theta < np.deg2rad(self.PM_ANGLE_THRESHOLD):
81 M2, p2, T2 = (flow.M, 0.0, flow.T)
82 method = -1
83 bad += 1
85 else:
86 # Use Prandtl-Meyer expansion theory
87 method = 1
88 M2, p2, T2 = self._solve_pm(
89 abs(theta), flow.M, flow.P, flow.T, flow.gamma
90 )
92 elif theta > 0:
93 # Use shock theory
94 beta_max = OPM.beta_max(M=flow.M, gamma=flow.gamma)
95 theta_max = OPM.theta_from_beta(
96 M1=flow.M, beta=beta_max, gamma=flow.gamma
97 )
98 if theta > theta_max:
99 # Detached shock
100 method = 2
101 M2, p2, T2 = self._solve_normal(
102 flow.M, flow.P, flow.T, flow.gamma
103 )
105 else:
106 # Use oblique shock theory
107 method = 3
108 M2, p2, T2 = self._solve_oblique(
109 theta, flow.M, flow.P, flow.T, flow.gamma
110 )
112 else:
113 # Cell is parallel to flow, assume no change
114 method = 0
115 M2, p2, T2 = (flow.M, flow.P, flow.T)
117 total += 1
119 # Save results for this cell
120 cell.attributes["pressure"] = p2
121 cell.attributes["Mach"] = M2
122 cell.attributes["temperature"] = T2
123 cell.attributes["method"] = method
125 # Add flowstate info to Cell
126 # TODO - calculate flowstate direction
127 cell.flowstate = FlowState(
128 mach=M2,
129 pressure=p2,
130 temperature=T2,
131 direction=None,
132 )
134 # Calculate force vector
135 F = cell.n * p2 * cell.A
136 net_force += F
137 net_moment += Vector.from_coordinates(np.cross(r.vec, F.vec))
139 # Update progress bar
140 if self.verbosity > 0:
141 pbar.update(1)
143 if self.verbosity > 0:
144 pbar.close()
146 if bad / total > 0.25:
147 print(
148 f"WARNING: {100*bad/total:.2f}% of cells were not "
149 "solved due to PM threshold."
150 )
152 # Construct results
153 result = FlowResults(
154 freestream=flow, net_force=net_force, net_moment=net_moment
155 )
157 # Save
158 self.flow_result = result
160 # Print result
161 if self.verbosity > 1:
162 print(result)
164 return result
166 @staticmethod
167 def pm(M: float, gamma: float = 1.4):
168 """Solves the Prandtl-Meyer function and returns the Prandtl-Meyer angle
169 in radians."""
170 v = ((gamma + 1) / (gamma - 1)) ** 0.5 * np.arctan(
171 ((M**2 - 1) * (gamma - 1) / (gamma + 1)) ** 0.5
172 ) - np.arctan((M**2 - 1) ** 0.5)
173 return v
175 @staticmethod
176 def inv_pm(angle: float, gamma: float = 1.4):
177 """Solves the inverse Prandtl-Meyer function using a bisection algorithm."""
178 func = lambda M: OPM.pm(M, gamma) - angle
179 pm = bisect(func, 1.0, 42.0)
180 return pm
182 @staticmethod
183 def _solve_pm(
184 theta: float, M1: float, p1: float = 1.0, T1: float = 1.0, gamma: float = 1.4
185 ):
186 """Solves for the Mach number, pressure and temperature after a Prandtl-Meyer
187 expansion fan.
189 Parameters
190 ----------
191 theta : float
192 The deflection angle, specified in radians.
194 M1 : float
195 The pre-expansion Mach number.
197 p1 : float, optional
198 The pre-expansion pressure (Pa). The default is 1.0.
200 T1 : float, optional
201 The pre-expansion temperature (K). The default is 1.0.
203 gamma : float, optional
204 The ratio of specific heats. The default is 1.4.
206 Returns
207 --------
208 M2 : float
209 The post-expansion Mach number.
211 p2 : float
212 The post-expansion pressure (Pa).
214 T2 : float
215 The post-expansion temperature (K).
216 """
217 # Solve for M2
218 v_M1 = OPM.pm(M=M1, gamma=gamma)
219 v_M2 = theta + v_M1
220 M2 = OPM.inv_pm(v_M2)
222 a = (gamma - 1) / 2
223 n = 1 + a * M1**2
224 d = 1 + a * M2**2
226 p2 = p1 * (n / d) ** (gamma / (gamma - 1))
228 T2 = T1 * (n / d)
230 return M2, p2, T2
232 @staticmethod
233 def _solve_oblique(
234 theta: float, M1: float, p1: float = 1.0, T1: float = 1.0, gamma: float = 1.4
235 ):
236 """Solves the flow using oblique shock theory.
238 Parameters
239 ----------
240 theta : float
241 The flow deflection angle, specified in radians.
243 M1 : float
244 The pre-shock Mach number.
246 p1 : float, optional
247 The pre-shock pressure (Pa). The default is 1.0.
249 T1 : float, optional
250 The pre-expansion temperature (K). The default is 1.0.
252 gamma : float, optional
253 The ratio of specific heats. The default is 1.4.
255 Returns
256 --------
257 M2 : float
258 The post-shock Mach number.
260 p2 : float
261 The post-shock pressure (Pa).
263 T2 : float
264 The post-shock temperature (K).
265 """
266 # Calculate angle and ratios
267 beta = OPM.oblique_beta(M1, theta, gamma)
268 p2_p1 = OPM.oblique_p2_p1(M1, beta, gamma)
269 rho2_rho1 = OPM.oblique_rho2_rho1(M1, beta, gamma)
270 T2_T1 = p2_p1 / rho2_rho1
272 # Calculate properties
273 M2 = OPM.oblique_M2(M1, beta, theta, gamma)
274 p2 = p2_p1 * p1
275 T2 = T1 * T2_T1
277 return M2, p2, T2
279 @staticmethod
280 def _solve_normal(M1: float, p1: float = 1.0, T1: float = 1.0, gamma: float = 1.4):
281 """Solves the flow using normal shock theory.
283 Parameters
284 ----------
285 M1 : float
286 The pre-shock Mach number.
288 p1 : float, optional
289 The pre-shock pressure (Pa). The default is 1.0.
291 T1 : float, optional
292 The pre-expansion temperature (K). The default is 1.0.
294 gamma : float, optional
295 The ratio of specific heats. The default is 1.4.
297 Returns
298 --------
299 M2 : float
300 The post-shock Mach number.
302 p2 : float
303 The post-shock pressure (Pa).
305 T2 : float
306 The post-shock temperature (K).
307 """
308 rho2_rho1 = (gamma + 1) * M1**2 / (2 + (gamma - 1) * M1**2)
309 p2_p1 = 1 + 2 * gamma * (M1**2 - 1) / (gamma + 1)
311 M2 = ((1 + M1 * (gamma - 1) / 2) / (gamma * M1**2 - (gamma - 1) / 2)) ** 0.5
312 p2 = p1 * p2_p1
313 T2 = T1 * p2_p1 / rho2_rho1
315 return M2, p2, T2
317 @staticmethod
318 def oblique_p2_p1(M1: float, beta: float, gamma: float = 1.4):
319 """Returns the pressure ratio p2/p1 across an oblique shock.
321 Parameters
322 -----------
323 M1 : float
324 The pre-shock Mach number.
326 beta : float
327 The shock angle specified in radians.
329 Returns
330 --------
331 p2/p1 : float
332 The pressure ratio across the shock.
334 References
335 ----------
336 Peter Jacobs
337 """
338 M1n = M1 * abs(np.sin(beta))
339 p2p1 = 1.0 + 2.0 * gamma / (gamma + 1.0) * (M1n**2 - 1.0)
340 return p2p1
342 @staticmethod
343 def oblique_beta(
344 M1: float, theta: float, gamma: float = 1.4, tolerance: float = 1.0e-6
345 ):
346 """Calculates the oblique shock angle using a recursive algorithm.
348 Parameters
349 ----------
350 M1 : float
351 The pre-shock Mach number.
353 theta : float
354 The flow deflection angle, specified in radians.
356 gamma : float, optional
357 The ratio of specific heats. The default is 1.4.
359 References
360 ----------
361 Peter Jacobs
362 """
363 func = lambda beta: OPM.theta_from_beta(M1, beta, gamma) - theta
365 # Initialise
366 sign_beta = np.sign(theta)
367 theta = abs(theta)
368 b1 = np.arcsin(1.0 / M1)
369 b2 = b1 * 1.05
371 # Check f1
372 f1 = func(b1)
373 if abs(f1) < tolerance:
374 return sign_beta * b1
376 # Check f2
377 f2 = func(b2)
378 if abs(f2) < tolerance:
379 return sign_beta * b2
381 # Finally, solve with secant method
382 # beta = sign_beta * secant(func, b1, b2, tol=tolerance)
383 root_result = root_scalar(f=func, x0=b1, x1=b2, xtol=tolerance, method="secant")
384 if root_result.converged:
385 beta = sign_beta * root_result.root
386 else:
387 raise Exception(
388 f"Cannot converge beta (theta={np.rad2deg(theta):.2f} degrees)."
389 )
391 return beta
393 @staticmethod
394 def theta_from_beta(M1: float, beta: float, gamma: float = 1.4):
395 """Calculates the flow deflection angle from the oblique shock angle.
397 Parameters
398 ----------
399 M1 : float
400 The pre-expansion Mach number.
402 beta : float
403 The shock angle specified in radians.
405 gamma : float, optional
406 The ratio of specific heats. The default is 1.4.
408 Returns
409 --------
410 theta : float
411 The deflection angle, specified in radians.
413 References
414 ----------
415 Peter Jacobs
416 """
417 M1n = M1 * abs(np.sin(beta))
418 t1 = 2.0 / np.tan(beta) * (M1n**2 - 1)
419 t2 = M1**2 * (gamma + np.cos(2 * beta)) + 2
420 theta = np.arctan(t1 / t2)
421 return theta
423 @staticmethod
424 def oblique_M2(M1: float, beta: float, theta: float, gamma: float = 1.4):
425 """Calculates the Mach number following an oblique shock.
427 Parameters
428 ----------
429 M1 : float
430 The pre-expansion Mach number.
432 beta : float
433 The shock angle specified in radians.
435 theta : float
436 The deflection angle, specified in radians.
438 gamma : float, optional
439 The ratio of specific heats. The default is 1.4.
441 Returns
442 --------
443 M2 : float
444 The post-shock Mach number.
446 References
447 ----------
448 Peter Jacobs
449 """
450 M1n = M1 * abs(np.sin(beta))
451 a = 1 + (gamma - 1) * 0.5 * M1n**2
452 b = gamma * M1n**2 - (gamma - 1) * 0.5
453 M2 = (a / b / (np.sin(beta - theta)) ** 2) ** 0.5
454 return M2
456 @staticmethod
457 def oblique_T2_T1(M1: float, beta: float, gamma: float = 1.4):
458 """Returns the temperature ratio T2/T1 across an oblique shock.
460 Parameters
461 -----------
462 M1 : float
463 The pre-shock Mach number.
465 beta : float
466 The shock angle specified in radians.
468 gamma : float, optional
469 The ratio of specific heats. The default is 1.4.
471 Returns
472 --------
473 T2/T1 : float
474 The temperature ratio across the shock.
476 References
477 ----------
478 Peter Jacobs
479 """
480 T2_T1 = OPM.oblique_p2_p1(M1, beta, gamma) / OPM.oblique_rho2_rho1(
481 M1, beta, gamma
482 )
483 return T2_T1
485 @staticmethod
486 def oblique_rho2_rho1(M1: float, beta: float, gamma: float = 1.4):
487 """Returns the density ratio rho2/rho1 across an oblique shock.
489 Parameters
490 -----------
491 M1 : float
492 The pre-shock Mach number.
494 beta : float
495 The shock angle specified in radians.
497 gamma : float, optional
498 The ratio of specific heats. The default is 1.4.
500 Returns
501 --------
502 T2/T1 : float
503 The temperature ratio across the shock.
505 References
506 ----------
507 Peter Jacobs
508 """
509 M1n = M1 * abs(np.sin(beta))
510 rho2_rho1 = (gamma + 1) * M1n**2 / (2 + (gamma - 1) * M1n**2)
511 return rho2_rho1
513 @staticmethod
514 def beta_max(M: float, gamma: float = 1.4):
515 """Returns the maximum shock angle for a given
516 Mach number.
517 """
518 beta_max = np.arcsin(
519 np.sqrt(
520 (1 / (gamma * M**2))
521 * (
522 (gamma + 1) * M**2 / 4
523 - 1
524 + np.sqrt(
525 (gamma + 1)
526 * ((gamma + 1) * M**4 / 16 + (gamma - 1) * M**2 / 2 + 1)
527 )
528 )
529 )
530 )
531 return beta_max
533 def solve_sens(
534 self,
535 sensitivity_filepath: Optional[str] = None,
536 parameters: Optional[list[str]] = None,
537 cells_have_sens_data: Optional[bool] = False,
538 freestream: Optional[FlowState] = None,
539 mach: Optional[float] = None,
540 aoa: Optional[float] = None,
541 cog: Vector = Vector(0, 0, 0),
542 perturbation: float = 1e-3,
543 ) -> SensitivityResults:
544 # TODO - add to cell attributes for visualisation
545 # TODO - return Mach and Temp sensitivities
547 already_run = super().solve_sens(freestream=freestream, mach=mach, aoa=aoa)
548 if already_run:
549 # Already have a result
550 if self.verbosity > 1:
551 print("Attention: this flowstate sensitivity has already been solved.")
552 result = self.flow_sensitivity
554 else:
555 # Get flow
556 flow = self._last_sens_freestream
558 # Check that nominal flow condition has been solved already
559 if self._last_solve_freestream:
560 # Solver has run before
561 if self._last_solve_freestream != flow:
562 # Run flow solver first at the nominal flow conditions
563 if self.verbosity > 0:
564 print("Running flow solver to prepare sensitivities.")
565 self.solve(flow)
566 else:
567 # Solver has not run yet at all, run it now
568 self.solve(flow)
570 if cells_have_sens_data and not parameters:
571 # Need to extract parameters
572 _, parameters = self._load_sens_data(sensitivity_filepath)
573 elif not cells_have_sens_data:
574 # Load sensitivity data and parameters
575 sensdata, parameters = self._load_sens_data(sensitivity_filepath)
577 # Add sensitivity data to cells
578 add_sens_data(cells=self.cells, data=sensdata, verbosity=self.verbosity)
580 # Construct progress bar
581 if self.verbosity > 0:
582 print()
583 desc = "Solving cell flow sensitivities."
584 pbar = tqdm(
585 total=len(self.cells),
586 desc=desc,
587 position=0,
588 leave=True,
589 )
591 # Iterate over cells
592 dFdp = 0
593 dMdp = 0
594 for cell in self.cells:
595 # Initialisation
596 all_directions = [Vector(1, 0, 0), Vector(0, 1, 0), Vector(0, 0, 1)]
597 sensitivities = np.empty(shape=(cell.dndp.shape[1], 3))
598 moment_sensitivities = np.empty(shape=(cell.dndp.shape[1], 3))
600 # Calculate moment dependencies
601 r = cell.c - cog
602 F = cell.flowstate.P * cell.A * cell.n.vec
604 # Calculate orientation of cell to flow
605 theta = np.pi / 2 - np.arccos(
606 np.dot(flow.direction.vec, cell.n.vec)
607 / (cell.n.norm * flow.direction.norm)
608 )
610 # Solve flow for this cell
611 if theta < 0:
612 # Rearward facing cell
613 if theta < np.deg2rad(self.PM_ANGLE_THRESHOLD):
614 _, dP_dtheta, _ = (0.0, 0.0, 0.0)
616 else:
617 # Use Prandtl-Meyer expansion theory
618 _, dP_dtheta, _ = self.dp_dtheta_pm(
619 abs(theta),
620 flow.M,
621 cell.flowstate,
622 flow.P,
623 flow.T,
624 flow.gamma,
625 perturbation,
626 )
628 elif theta > 0:
629 # Use shock theory
630 beta_max = OPM.beta_max(M=flow.M, gamma=flow.gamma)
631 theta_max = OPM.theta_from_beta(
632 M1=flow.M, beta=beta_max, gamma=flow.gamma
633 )
634 if theta > theta_max:
635 # Detached shock
636 _, dP_dtheta, _ = (0.0, 0.0, 0.0)
638 else:
639 # Use oblique shock theory
640 _, dP_dtheta, _ = self.dp_dtheta_obl(
641 theta,
642 flow.M,
643 cell.flowstate,
644 flow.P,
645 flow.T,
646 flow.gamma,
647 perturbation,
648 )
650 else:
651 # Cell is parallel to flow, assume no change
652 _, dP_dtheta, _ = (0.0, 0.0, 0.0)
654 # For each parameter
655 for p_i in range(cell.dndp.shape[1]):
656 # Calculate gradient of pressure with respect to parameter
657 dtheta_dp = (
658 -np.dot(cell.dndp[:, p_i], flow.direction.vec)
659 / (1 - np.dot(cell.n.unit.vec, flow.direction.unit.vec) ** 2)
660 ** 0.5
661 )
663 if np.isnan(dtheta_dp):
664 # Normal vector parallel to flow
665 dtheta_dp = 0
667 dPdp = dP_dtheta * dtheta_dp
669 # Calculate pressure sensitivity for each direction
670 for i, direction in enumerate(all_directions):
671 dF = (
672 dPdp * cell.A * np.dot(cell.n.vec, direction.vec)
673 + cell.flowstate.P
674 * cell.dAdp[p_i]
675 * np.dot(cell.n.vec, direction.vec)
676 + cell.flowstate.P
677 * cell.A
678 * np.dot(-cell.dndp[:, p_i], direction.vec)
679 )
680 sensitivities[p_i, i] = dF
682 # Now evaluate moment sensitivities
683 moment_sensitivities[p_i, :] = np.cross(
684 r.vec, sensitivities[p_i, :]
685 ) + np.cross(cell.dcdp[:, p_i], F)
687 # Append to cell
688 cell.sensitivities = sensitivities
689 cell.moment_sensitivities = moment_sensitivities
691 # Update
692 dFdp += sensitivities
693 dMdp += moment_sensitivities
695 # Update progress bar
696 if self.verbosity > 0:
697 pbar.update(1)
699 if self.verbosity > 0:
700 pbar.close()
701 print("Done.")
703 # Construct dataframes to return
704 df_f = pd.DataFrame(
705 dFdp, columns=["dFx/dp", "dFy/dp", "dFz/dp"], index=parameters
706 )
707 df_m = pd.DataFrame(
708 dMdp, columns=["dMx/dp", "dMy/dp", "dMz/dp"], index=parameters
709 )
711 # Construct results
712 result = SensitivityResults(
713 f_sens=df_f,
714 m_sens=df_m,
715 freestream=flow,
716 )
718 # Save
719 self.flow_sensitivity = result
721 if self.verbosity > 0:
722 print(result)
724 return result
726 @staticmethod
727 def dp_dtheta_obl(
728 theta: float,
729 M1: float,
730 nominal_flowstate: FlowState,
731 p1: float = 1.0,
732 T1: float = 1.0,
733 gamma: float = 1.4,
734 perturbation: float = 1e-3,
735 ):
736 """Solves for the sensitivities at the given point using oblique shock
737 theory.
739 Parameters
740 ----------
741 theta : float
742 The deflection angle, specified in radians.
744 M1 : float
745 The pre-expansion Mach number.
747 nominal_flowstate : FlowState
748 The nominal flowstate at this point.
750 p1 : float, optional
751 The pre-expansion pressure (Pa). The default is 1.0.
753 T1 : float, optional
754 The pre-expansion temperature (K). The default is 1.0.
756 gamma : float, optional
757 The ratio of specific heats. The default is 1.4.
759 perturbation : float, optional
760 The perturbation amount. Perturbations are calculated
761 according to a multiple of (1 +/- perturbation). The
762 default is 1e-3.
764 Returns
765 --------
766 dM_dtheta : float
767 The sensitivity of the post-expansion Mach number with respect
768 to the deflection angle theta.
770 dP_dtheta : float
771 The sensitivity of the post-expansion pressure (Pa) with respect
772 to the deflection angle theta.
774 dT_dtheta : float
775 The sensitivity of the post-expansion temperature (K) with respect
776 to the deflection angle theta.
777 """
778 # Create function handle
779 func = lambda theta: OPM._solve_oblique(theta, M1, p1, T1, gamma)
781 # Calculate derivitive by finite differencing
782 g = OPM._findiff(func, theta, nominal_flowstate, perturbation)
784 return g
786 @staticmethod
787 def dp_dtheta_pm(
788 theta: float,
789 M1: float,
790 nominal_flowstate: FlowState,
791 p1: float = 1.0,
792 T1: float = 1.0,
793 gamma: float = 1.4,
794 perturbation: float = 1e-3,
795 ):
796 """Solves for the sensitivities at the given point using Prandtl-Meyer
797 theory.
799 Parameters
800 ----------
801 theta : float
802 The deflection angle, specified in radians.
804 M1 : float
805 The pre-expansion Mach number.
807 nominal_flowstate : FlowState
808 The nominal flowstate at this point.
810 p1 : float, optional
811 The pre-expansion pressure (Pa). The default is 1.0.
813 T1 : float, optional
814 The pre-expansion temperature (K). The default is 1.0.
816 gamma : float, optional
817 The ratio of specific heats. The default is 1.4.
819 perturbation : float, optional
820 The perturbation amount. Perturbations are calculated
821 according to a multiple of (1 +/- perturbation). The
822 default is 1e-3.
824 Returns
825 --------
826 dM_dtheta : float
827 The sensitivity of the post-expansion Mach number with respect
828 to the deflection angle theta.
830 dP_dtheta : float
831 The sensitivity of the post-expansion pressure (Pa) with respect
832 to the deflection angle theta.
834 dT_dtheta : float
835 The sensitivity of the post-expansion temperature (K) with respect
836 to the deflection angle theta.
837 """
838 # Create function handle
839 func = lambda theta: OPM._solve_pm(theta, M1, p1, T1, gamma)
841 # Calculate derivitive by finite differencing
842 g = OPM._findiff(func, theta, nominal_flowstate, perturbation)
844 return g
846 @staticmethod
847 def _findiff(
848 func: callable,
849 point: any,
850 nominal_flowstate: FlowState,
851 perturbation: float = 1e-3,
852 ):
853 # Generate flow value at points +/- perturbed from nominal point
854 xvals = [point * (1 - perturbation), point, point * (1 + perturbation)]
855 vals = [
856 func(p) for p in [point * (1 - perturbation), point * (1 + perturbation)]
857 ]
858 Mvals = [v[0] for v in vals]
859 pvals = [v[1] for v in vals]
860 Tvals = [v[2] for v in vals]
862 # Inject nominal result
863 Mvals = [Mvals[0], nominal_flowstate.M, Mvals[1]]
864 pvals = [pvals[0], nominal_flowstate.P, pvals[1]]
865 Tvals = [Tvals[0], nominal_flowstate.T, Tvals[1]]
867 dM = np.gradient(Mvals, xvals)
868 dp = np.gradient(pvals, xvals)
869 dT = np.gradient(Tvals, xvals)
871 return dM[1], dp[1], dT[1]
873 def save(self, name: str):
874 # Initialise attributes dictionary
875 attributes = {
876 "pressure": [],
877 "Mach": [],
878 "temperature": [],
879 "method": [],
880 }
882 # Call super method
883 super().save(name, attributes)
885 @staticmethod
886 def _load_sens_data(sensitivity_filepath: str) -> Tuple[pd.DataFrame, list[set]]:
887 """Extracts design parameters from the sensitivity data."""
888 sensdata = pd.read_csv(sensitivity_filepath)
890 parameters = set()
891 for e in sensdata.columns:
892 e: str
893 if e.startswith("dxd") or e.startswith("dyd") or e.startswith("dzd"):
894 # Sensitivity coluns
895 parameters.add(e[3:])
897 return sensdata, list(parameters)