[WIP] Boundary update with testcases - #27
Conversation
|
flamelet_table_h2 was not working anymore, still has to be updated and checked. |
EvertBunschoten
left a comment
EvertBunschoten
left a comment
There was a problem hiding this comment.
Very nice! Thank you also for doing some clean up in other parts of the code and the test cases. I especially like the extensive visualization functions that have been added to the table generator. Overall, I think some of the functions can be made shorter and there is still one bug where the table plots are not actually visualized when specifying show=True.
| def ClampSourceTerms(self, species_list:list[str], pv_frac:float=0.99, abs_tol:float=1e-3): | ||
| """Clamp source terms of selected species to zero near the burnt (high-PV) boundary. | ||
|
|
||
| For each table level the per-level progress variable maximum is computed. At every node | ||
| where PV >= ``pv_frac * PV_max`` the net, positive, and negative source terms | ||
| (``Y_dot_net-<sp>``, ``Y_dot_pos-<sp>``, ``Y_dot_neg-<sp>``) of each species in | ||
| ``species_list`` are set to zero when their absolute value is below ``abs_tol``. | ||
| In addition, the production terms are clipped to be non-negative and the consumption | ||
| terms non-positive over all nodes, as SU2 evaluates the species source as | ||
| ``source_prod + source_cons * y_aux``. | ||
|
|
||
| Call this method after :meth:`generateTable` and before :meth:`writeSU2Table`. | ||
|
|
||
| :param species_list: species names to clamp (e.g. ``['CO', 'H2']``). | ||
| :type species_list: list[str] | ||
| :param pv_frac: fraction of the PV maximum above which clamping is applied, defaults to 0.99 | ||
| :type pv_frac: float, optional | ||
| :param abs_tol: source terms with an absolute value below this are set to zero, defaults to 1e-3 | ||
| :type abs_tol: float, optional | ||
| :raises Exception: if the table data have not been generated yet. |
There was a problem hiding this comment.
Try to make this functionality independent of the order of operation, similar to how smoothing is applied.
| data_level = self._data_in_table[iLevel] | ||
| pv_nodes = data_level[name_pv].to_numpy(dtype=float) | ||
| pv_max = pv_max_global if pv_is_level_cv else np.max(pv_nodes) | ||
| nodes_near_products = pv_nodes >= pv_frac * pv_max |
There was a problem hiding this comment.
Why only apply this to the burnt side? Would it be possible to apply clamping to the premixed side as well?
There was a problem hiding this comment.
Sure, but for burner stabilized flames the source term does not have to be zero on the burner surface.
There was a problem hiding this comment.
We can use the EquilibriumSolver to calculate the value of the progress variable of the premixed reactants and products. The clamping should not affect the source terms there, since the value of the progress variable at the inflow of the burner-stabilized flamelet solutions is higher than the reactant progress variable.
| def __polylineIntersection(self, polyline_a:np.ndarray[float], polyline_b:np.ndarray[float]): | ||
| """Locate the last point at which two polylines cross. | ||
|
|
||
| :return: the crossing point with the index of the crossed segment in each polyline, or None. | ||
| :rtype: tuple | ||
| """ | ||
| crossing = None | ||
| for iSegment in range(len(polyline_a) - 1): | ||
| start_a = polyline_a[iSegment] | ||
| direction_a = polyline_a[iSegment + 1] - start_a | ||
| for jSegment in range(len(polyline_b) - 1): | ||
| start_b = polyline_b[jSegment] | ||
| direction_b = polyline_b[jSegment + 1] - start_b | ||
|
|
||
| determinant = direction_a[0]*direction_b[1] - direction_a[1]*direction_b[0] | ||
| if abs(determinant) < np.finfo(float).eps: | ||
| continue | ||
|
|
||
| offset = start_b - start_a | ||
| along_a = (offset[0]*direction_b[1] - offset[1]*direction_b[0])/determinant | ||
| along_b = (offset[0]*direction_a[1] - offset[1]*direction_a[0])/determinant | ||
| if 0 <= along_a <= 1 and 0 <= along_b <= 1: | ||
| crossing = (start_a + along_a*direction_a, iSegment, jSegment) | ||
| return crossing |
There was a problem hiding this comment.
This function can be made a static method, for example under Common/Interpolators.py
| if not show: | ||
| plt.close(fig) |
There was a problem hiding this comment.
| if not show: | |
| plt.close(fig) | |
| if show: | |
| plt.show() | |
| if else: | |
| plt.close(fig) | |
| if not table_is_generated: | ||
| self._processTableLevels() | ||
| if self._fluid_data_interpolator is None: | ||
| self._defineFluidDataInterpolator() | ||
| self._checkTableLevelIndex(level_index) | ||
|
|
||
| if table_is_generated: | ||
| table_nodes = self._table_nodes[level_index] | ||
| connectivity = self._table_connectivity[level_index] | ||
| hull_nodes = self._table_hullnodes[level_index] | ||
| data_level = self._data_in_table[level_index] | ||
| else: | ||
| table_nodes, connectivity, hull_nodes, data_level = self.meshTableLevel(level_index) |
There was a problem hiding this comment.
| if not table_is_generated: | |
| self._processTableLevels() | |
| if self._fluid_data_interpolator is None: | |
| self._defineFluidDataInterpolator() | |
| self._checkTableLevelIndex(level_index) | |
| if table_is_generated: | |
| table_nodes = self._table_nodes[level_index] | |
| connectivity = self._table_connectivity[level_index] | |
| hull_nodes = self._table_hullnodes[level_index] | |
| data_level = self._data_in_table[level_index] | |
| else: | |
| table_nodes, connectivity, hull_nodes, data_level = self.meshTableLevel(level_index) | |
| if table_is_generated: | |
| table_nodes = self._table_nodes[level_index] | |
| connectivity = self._table_connectivity[level_index] | |
| hull_nodes = self._table_hullnodes[level_index] | |
| data_level = self._data_in_table[level_index] | |
| else: | |
| self._processTableLevels() | |
| if self._fluid_data_interpolator is None: | |
| self._defineFluidDataInterpolator() | |
| self._checkTableLevelIndex(level_index) | |
| table_nodes, connectivity, hull_nodes, data_level = self.meshTableLevel(level_index) |
| fig, ax = plt.subplots(figsize=(10, 10)) | ||
| ax.triplot(x_nodes, y_nodes, connectivity, linewidth=0.5) | ||
| ax.plot(x_nodes[hull_nodes], y_nodes[hull_nodes], 'ko', ms=3, label="Perimiter nodes") | ||
| ax.set_xlabel(self._table_cv_names[0], fontsize=14) | ||
| ax.set_ylabel(self._table_cv_names[1], fontsize=14) | ||
| ax.set_title(title, fontsize=14) | ||
| ax.legend(fontsize=12) |
There was a problem hiding this comment.
Could these blocks be moved to static functions? The plotting data, labels, etc, could be passed through a dictionary.
| triangulation = mtri.Triangulation(x_nodes, y_nodes, connectivity) | ||
| fig, ax = plt.subplots(figsize=(10, 7), constrained_layout=True) | ||
| colormesh = ax.tripcolor(triangulation, data_level[var_to_plot].to_numpy(dtype=float), \ | ||
| shading='gouraud', cmap='inferno') |
There was a problem hiding this comment.
If we define the colormap through a setter function, we can keep the colormap consistent between plotting functions.
| raise Exception("Boundary polyline should be provided as an array with %i columns" % len(self.__mesh_along_coords)) | ||
| if len(polyline) < 3: | ||
| raise Exception("Boundary polyline should describe at least three points") | ||
| self.__boundary_polyline = polyline |
There was a problem hiding this comment.
There should be one more check here for if the point cloud contains unique points.
2 participants
Proposed Changes
Modify 2D table edges to better follow curved boundaries
add fixed strain rate testcase
add CH4 premixed testcase
To do for this PR:
regression test the new cases
Related Work
PR Checklist
Put an X by all that apply. You can fill this out after submitting the PR. If you have any questions, don't hesitate to ask! We want to help. These are a guide for you to know what the reviewers will be looking for in your contribution.
pre-commit run --allto format old commits.