Coverage for /opt/hostedtoolcache/Python/3.11.10/x64/lib/python3.11/site-packages/pysagas/sensitivity/cart3d/cart3d.py: 83%

1import numpy as np 

2import pandas as pd 

3from tqdm import tqdm 

4from typing import List, Optional 

5from pysagas.flow import FlowState 

6from pysagas.sensitivity.calculator import SensitivityCalculator 

7from pysagas.geometry import Vector, Cell, DegenerateCell 

8from pysagas.sensitivity.cart3d.utilities import process_components_file 

9 

10 

11class Cart3DSensitivityCalculator(SensitivityCalculator): 

12 """PySAGAS Cart3D flow sensitivity calculator.""" 

13 

14 solver = "Cart3D" 

15 

16 def __init__( 

17 self, 

18 freestream: FlowState, 

19 sensitivity_filepath: str, 

20 components_filepath: Optional[str] = None, 

21 pointdata: Optional[pd.DataFrame] = None, 

22 celldata: Optional[pd.DataFrame] = None, 

23 write_data: Optional[bool] = False, 

24 verbosity: Optional[int] = 1, 

25 **kwargs, 

26 ) -> None: 

27 """A PySAGAS sensitivity calculator for Cart3D. 

28 

29 Parameters 

30 ---------- 

31 freestream : FlowState 

32 The flow state of the freestream. 

33 

34 sensitivity_filepath : str 

35 The filepath to the geometry sensitivities. 

36 

37 components_filepath : str, optional 

38 The filepath to the Components.i.plt file to be processed. 

39 The default is None. 

40 

41 pointdata : pd.DataFrame, optional 

42 The point data. Must be supplied with celldata. 

43 

44 celldata : pd.DataFrame, optional 

45 The cell data. Must be supplied with pointdata. 

46 

47 write_data : bool, optional 

48 Write the flow data to CSV files. The default is True. 

49 

50 verbosity : int, optional 

51 The verbosity of the code. The defualt is 1. 

52 """ 

53 # Load data 

54 self.sensdata = pd.read_csv(sensitivity_filepath) 

55 if pointdata is not None and celldata is not None: 

56 self.pointdata = pointdata 

57 self.celldata = celldata 

58 else: 

59 self.pointdata, self.celldata = process_components_file( 

60 freestream.a, 

61 freestream.rho, 

62 components_filepath, 

63 write_data=write_data, 

64 verbosity=verbosity, 

65 ) 

66 self.verbosity = verbosity 

67 

68 # Check dimensionality 

69 l1 = len(self.sensdata) 

70 l2 = len(self.pointdata) 

71 if l1 != l2: 

72 raise ValueError( 

73 f"The sensitivity data does not match the point data ({l1} vs. {l2})." 

74 ) 

75 

76 super().__init__(**kwargs) 

77 

78 def _transcribe_cells(self, parameters: List[str]) -> List[Cell]: 

79 """Transcribes the cells from Components.i.plt files into 

80 PySAGAS Cell objects. 

81 

82 Parameters 

83 ----------- 

84 parameters : List[str] 

85 A list of the geometric design parameters. 

86 

87 Returns 

88 -------- 

89 cells : List[Cell] 

90 A list of all transcribed cells. 

91 

92 See Also 

93 -------- 

94 pysagas.geometry.Cell 

95 """ 

96 coordinates = self.pointdata[["Points_0", "Points_1", "Points_2"]] 

97 cell_vertex_ids = self.celldata[ 

98 ["Point Index 0", "Point Index 1", "Point Index 2"] 

99 ] 

100 

101 # Construct cells 

102 if self.verbosity > 0: 

103 print("\nTranscribing cells from Components.i.plt...") 

104 pbar = tqdm( 

105 total=len(self.celldata.index), 

106 position=0, 

107 leave=True, 

108 ) 

109 cells = [] 

110 degen_cells = 0 

111 for cell in self.celldata.index: 

112 vertex_ids = cell_vertex_ids.loc[cell].values 

113 cell_v_coords = [coordinates.loc[v_id].values for v_id in vertex_ids] 

114 vertices = [Vector.from_coordinates(c) for c in cell_v_coords] 

115 

116 # Extract sensitivity information for the cell 

117 vertex_sensitivity_info = [self.sensdata.loc[v_id] for v_id in vertex_ids] 

118 

119 dvdp = np.empty((9, len(parameters))) 

120 for i, p in enumerate(parameters): 

121 r = 0 

122 for v in vertex_sensitivity_info: 

123 for c in ["x", "y", "z"]: 

124 dvdp[r, i] = v[f"d{c}d{p}"] 

125 r += 1 

126 

127 # Create Cell 

128 try: 

129 newcell = Cell.from_points(vertices) 

130 except DegenerateCell: 

131 # Bad cell, skip it 

132 degen_cells += 1 

133 if self.verbosity > 0: 

134 pbar.update(1) 

135 

136 if self.verbosity > 2: 

137 print("\033[1mWarning\033[0m: Degenerate cell.") 

138 continue 

139 

140 # Add geometry sensitivity information 

141 newcell._add_sensitivities(np.array(dvdp)) 

142 

143 # Add flow state at cell 

144 cell_velocity = Vector.from_coordinates( 

145 self.celldata.loc[cell][["U", "V", "W"]].values 

146 ) 

147 cell_flowstate = FlowState( 

148 mach=self.celldata.loc[cell]["M"], 

149 pressure=self.celldata.loc[cell]["p"], 

150 temperature=self.celldata.loc[cell]["T"], 

151 direction=cell_velocity, 

152 ) 

153 newcell.flowstate = cell_flowstate 

154 

155 # Append new Cell 

156 cells.append(newcell) 

157 

158 # Update progress bar 

159 if self.verbosity > 0: 

160 pbar.update(1) 

161 

162 if self.verbosity > 0: 

163 pbar.close() 

164 print("Done.") 

165 

166 if self.verbosity > 1: 

167 if degen_cells > 0: 

168 print(f"{100*degen_cells/len(self.celldata):.2f}% degenerate cells") 

169 

170 return cells