diff --git a/src/GuiSlots.py b/src/GuiSlots.py index ed04306..f816e32 100644 --- a/src/GuiSlots.py +++ b/src/GuiSlots.py @@ -340,6 +340,11 @@ def messageBox(self, message): information(self.parent, 'Information', message, QtWidgets.QMessageBox.Ok) + @QtCore.Slot(str) + def yesnocancelQuestionBox(self, question): + return QtWidgets.QMessageBox.information(self.parent, 'Confirm?', question, + QtWidgets.QMessageBox.Yes | QtWidgets.QMessageBox.No | QtWidgets.QMessageBox.Cancel) + @QtCore.Slot() def onKeyBd(self): # automatically populate shortcuts from PMenu.xml diff --git a/src/MatplotlibWindow.py b/src/MatplotlibWindow.py new file mode 100644 index 0000000..0dda1da --- /dev/null +++ b/src/MatplotlibWindow.py @@ -0,0 +1,56 @@ +import sys +import threading +import random +import time + +import pandas as pd + +from PySide6.QtCore import Qt, QCoreApplication +from PySide6.QtGui import QWheelEvent + +from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg as FigureCanvas +from matplotlib.figure import Figure +from matplotlib.lines import Line2D + + +class MplCanvas(FigureCanvas): + """ + Provides an embeddable matplotlib canvas that can be live-updated + Preserves scrolling ability in GUI + """ + def __init__(self): + self.fig = Figure() + self.ax = self.fig.add_subplot(111) + self.clear() + super().__init__(self.fig) + + def clear(self): + self.ax.clear() + self.ax.grid() + self.ax.set_xlabel("Time iteration") + self.ax.set_ylabel("Value") + self.lines = {} + + def update_plot(self, label: str, data: pd.Series, **kwargs): + if shift_max := kwargs.get(f"{label}_trend", 50): + trend = pd.DataFrame({"_dat": data.values}, index=data.index)["_dat"].ewm(halflife=shift_max).mean() + if not label in self.lines: + line_kwargs = kwargs.get(label, {}) + if "label" in line_kwargs: + print("ERROR: passing label not allowed.") + self.lines[label] = self.ax.plot(list(data.index), list(data.values), **line_kwargs, alpha=0.5, label=label)[0] + # ax.plot returns a list of Line2D objects (hence the [0]) + if shift_max: + line_kwargs["color"] = self.lines[label].get_color() + self.lines[f"{label}_trend"] = self.ax.plot(list(trend.index), list(trend.values), **line_kwargs, label=f"{label}_trend")[0] + self.ax.legend() + else: + self.lines[label].set_data(list(data.index), list(data.values)) + if shift_max: + self.lines[f"{label}_trend"].set_data(list(trend.index), list(trend.values)) + + def wheelEvent(self, event: QWheelEvent): + if event.modifiers() & Qt.KeyboardModifier.ControlModifier: + super().wheelEvent(event) # Let Matplotlib zoom + else: + QCoreApplication.sendEvent(self.parent(), event) diff --git a/src/Meshing.py b/src/Meshing.py index 5a94940..fc718de 100644 --- a/src/Meshing.py +++ b/src/Meshing.py @@ -1238,31 +1238,12 @@ def writeSU2_nolib(wind_tunnel, name=''): # write boundary tags f.write('NMARK= 5\n') - f.write('MARKER_TAG= airfoil\n') - f.write('MARKER_ELEMS= ' + str(num_airfoil_edges) + '\n') - for edge in tags['airfoil']: - f.write(f'3 {edge[0]} {edge[1]}\n') - - f.write('MARKER_TAG= inlet\n') - f.write('MARKER_ELEMS= ' + str(num_inlet_edges) + '\n') - for edge in tags['inlet']: - f.write(f'3 {edge[0]} {edge[1]}\n') - - f.write('MARKER_TAG= outlet\n') - f.write('MARKER_ELEMS= ' + str(num_outlet_edges) + '\n') - for edge in tags['outlet']: - f.write(f'3 {edge[0]} {edge[1]}\n') - - f.write('MARKER_TAG= top\n') - f.write('MARKER_ELEMS= ' + str(num_top_edges) + '\n') - for edge in tags['top']: - f.write(f'3 {edge[0]} {edge[1]}\n') - - f.write('MARKER_TAG= bottom\n') - f.write('MARKER_ELEMS= ' + str(num_bottom_edges) + '\n') - for edge in tags['bottom']: - f.write(f'3 {edge[0]} {edge[1]}\n') - + for tag, edges in tags.items(): + f.write(f'MARKER_TAG= {wind_tunnel.boundary_tags_names[tag]}\n') + f.write(f'MARKER_ELEMS= {str(len(edges))}\n') + for edge in edges: + f.write(f'3 {edge[0]} {edge[1]}\n') + basename = os.path.basename(name) logger.info('SU2 type mesh saved as {}'. format(os.path.join(OUTPUTDATA, basename))) diff --git a/src/PyAero.py b/src/PyAero.py index e1945e4..1fccd84 100644 --- a/src/PyAero.py +++ b/src/PyAero.py @@ -149,6 +149,14 @@ def keyPressEvent(self, event): # progress event super().keyPressEvent(event) + def closeEvent(self, event): + print("custom close event!") + try: + self.toolbox.abortSU2() # FIXME: this doesn't seem to do anything. wrong class's closeEvent?! + except Exception as excpt: + print(f"Aborting simulation processes failed with {excpt}") + self.deleteLater() + class CentralWidget(QtWidgets.QWidget): """ diff --git a/src/Settings.py b/src/Settings.py index fea5452..d8fde42 100644 --- a/src/Settings.py +++ b/src/Settings.py @@ -27,6 +27,12 @@ # path to log files LOGDATA = os.path.join(DATAPATH, 'LOGS') +# path to su2 master-config +SU2CONFIG = os.path.join(DATAPATH, "Config") + +# path to SU2 executable +SU2EXEC = "SU2_CFD" # this should work if it's in PATH, otherwise just put the full path here + # set locale # can be either 'C' or '' # if string is empty then system default locale is used diff --git a/src/ToolBox.py b/src/ToolBox.py index 07fe226..3df0bf8 100644 --- a/src/ToolBox.py +++ b/src/ToolBox.py @@ -1,20 +1,29 @@ # -*- coding: utf-8 -*- import os +import re +import sys +import time import numpy as np +import pandas as pd +import signal from PySide6 import QtGui, QtCore, QtWidgets +from PySide6.QtCore import QTimer import PyAero import Airfoil import FileDialog import FileSystem import SvpMethod +import su2 import SplineRefine import TrailingEdge import Meshing import ContourAnalysis as ca -from Settings import ICONS_L +import MatplotlibWindow +from Settings import ICONS_L, SU2CONFIG, OUTPUTDATA, SU2EXEC +import queue import logging logger = logging.getLogger(__name__) @@ -53,6 +62,7 @@ def __init__(self, parent): # create toolbox items self.itemFileSystem() self.itemAeropython() + self.itemSU2() self.itemBoundaryCondtions() self.itemContourAnalysis() self.itemSplineRefine() @@ -163,6 +173,100 @@ def itemAeropython(self): panelMethodButton.clicked.connect(self.runPanelMethod) + def itemSU2(self): + + form = QtWidgets.QFormLayout() + + label1 = QtWidgets.QLabel(u'Angle of attack (°)') + self.aoaSU2 = QtWidgets.QDoubleSpinBox() + self.aoaSU2.setSingleStep(0.1) + self.aoaSU2.setDecimals(1) + self.aoaSU2.setRange(-10.0, 10.0) + self.aoaSU2.setValue(0.0) + form.addRow(label1, self.aoaSU2) + + label2 = QtWidgets.QLabel('Reynolds number') + self.reynolds_su2 = QtWidgets.QDoubleSpinBox() + self.reynolds_su2.setSingleStep(1e2) + self.reynolds_su2.setDecimals(0) + self.reynolds_su2.setRange(1, 1e8) + self.reynolds_su2.setValue(1e6) + form.addRow(label2, self.reynolds_su2) + + label2 = QtWidgets.QLabel('Freestream Mach number') + self.freestreamSU2 = QtWidgets.QDoubleSpinBox() + self.freestreamSU2.setSingleStep(0.1) + self.freestreamSU2.setDecimals(2) + self.freestreamSU2.setRange(0.0, 100.0) + self.freestreamSU2.setValue(0.4) + form.addRow(label2, self.freestreamSU2) + + label3 = QtWidgets.QLabel('Mesh file') + self.meshfile = QtWidgets.QLineEdit('msh.su2') + form.addRow(label3, self.meshfile) + + label4 = QtWidgets.QLabel(u'Number of cores') + self.n_cores = QtWidgets.QSpinBox() + self.n_cores.setSingleStep(1) + self.n_cores.setRange(1, 64) + self.n_cores.setValue(1) + form.addRow(label4, self.n_cores) + + label4 = QtWidgets.QLabel(u'Restart-file write frequency') + self.n_restartfreq = QtWidgets.QDoubleSpinBox() + self.n_restartfreq.setSingleStep(1) + self.n_restartfreq.setDecimals(0) + self.n_restartfreq.setRange(1, 500) # TODO: allow this to be 0, meaning no restart file + self.n_restartfreq.setValue(100) + form.addRow(label4, self.n_restartfreq) + + label5 = QtWidgets.QLabel(u'Flow-field write frequency') + self.n_flowfreq = QtWidgets.QDoubleSpinBox() + self.n_flowfreq.setSingleStep(1) + self.n_flowfreq.setDecimals(0) + self.n_flowfreq.setRange(1, 500) # TODO: allow this to be 0, meaning no flow file (for whatever reason) + self.n_flowfreq.setValue(5) + form.addRow(label5, self.n_flowfreq) + + n_col = 3 + options = ['"CD"', '"CL"', '"rms[RhoU]"', '"rms[RhoE]"'] # TODO: add others? + self.residuals_cbs = {} + label = QtWidgets.QLabel('Residuals / Quantities to live-plot') + label.setToolTip('Check quantities to plot live') + grid = QtWidgets.QGridLayout() + grid.addWidget(label, 0, 0) + for i, var in enumerate(options): + cb = QtWidgets.QCheckBox(var) + self.residuals_cbs[var] = cb + cb.setChecked(True) + grid.addWidget(cb, i // n_col + 2, i % n_col + 1) + form.addRow(grid) + + label6 = QtWidgets.QLabel(u'Moving-average samples') + self.n_residuals_movingavg = QtWidgets.QDoubleSpinBox() + self.n_residuals_movingavg.setSingleStep(1) + self.n_residuals_movingavg.setDecimals(0) + self.n_residuals_movingavg.setRange(0, 20) + self.n_residuals_movingavg.setValue(5) + form.addRow(label6, self.n_residuals_movingavg) + + runSU2Button = QtWidgets.QPushButton('Run solver') + form.addRow(runSU2Button) + abortSU2Button = QtWidgets.QPushButton('Abort solver') + form.addRow(abortSU2Button) + + self.residuals_canvas = MatplotlibWindow.MplCanvas() + self.residuals_canvas.setMinimumSize(400, 300) + # self.residuals_canvas.setSizePolicy(Q) + form.addRow(self.residuals_canvas) + + self.item_su2 = QtWidgets.QGroupBox('SU2') + self.item_su2.setLayout(form) + + self.simulation_initialised = False + runSU2Button.clicked.connect(self.runSU2) + abortSU2Button.clicked.connect(self.abortSU2) + def itemBoundaryCondtions(self): form = QtWidgets.QFormLayout() @@ -639,30 +743,20 @@ def itemMeshing(self): header_2.setStyleSheet('font-weight: bold;') self.form_bnd.addRow(header_1, header_2) - label = QtWidgets.QLabel('Airfoil') - label.setToolTip('Name of the boundary definition for the airfoil') - self.lineedit_airfoil = QtWidgets.QLineEdit('Airfoil') - self.form_bnd.addRow(label, self.lineedit_airfoil) - - label = QtWidgets.QLabel('Inlet (C-arc)') - label.setToolTip('Name of the boundary definition for the inlet') - self.lineedit_inlet = QtWidgets.QLineEdit('Inlet') - self.form_bnd.addRow(label, self.lineedit_inlet) - - label = QtWidgets.QLabel('Outlet') - label.setToolTip('Name of the boundary definition for the outlet') - self.lineedit_outlet = QtWidgets.QLineEdit('Outlet') - self.form_bnd.addRow(label, self.lineedit_outlet) - - label = QtWidgets.QLabel('Top') - label.setToolTip('Name of the boundary definition for the top of the windtunnel') - self.lineedit_top = QtWidgets.QLineEdit('Top') - self.form_bnd.addRow(label, self.lineedit_top) - - label = QtWidgets.QLabel('Bottom') - label.setToolTip('Name of the boundary definition for the bottom of the windtunnel') - self.lineedit_bottom = QtWidgets.QLineEdit('Bottom') - self.form_bnd.addRow(label, self.lineedit_bottom) + self.lineedits_boundarymarkers = {} + gui_labels = { + "airfoil": "Airfoil", + "inlet": "Inlet (C-arc)", + "outlet": "Outlet", + "top": "Top of Windtunnel", + "bottom": "Bottom of Windtunnel" + } + for boundary, label in gui_labels.items(): + label = QtWidgets.QLabel(label.split(None, 1)[0]) + label.setToolTip(f'Name of the boundary definition for the {label}') + edit = QtWidgets.QLineEdit(boundary) + self.lineedits_boundarymarkers[boundary] = edit + self.form_bnd.addRow(label, edit) self.check_FIRE = QtWidgets.QCheckBox('AVL FIRE') self.check_SU2 = QtWidgets.QCheckBox('SU2') @@ -859,6 +953,7 @@ def makeToolbox(self): self.tb6 = self.addItem(self.item_abc, 'CFD Boundary Conditions') self.tb5 = self.addItem(self.item_ap, 'Aerodynamics (Panel Code)') + self.tb7 = self.addItem(self.item_su2, 'Aerodynamics (SU2)') self.tb3 = self.addItem(self.item_ca, 'Contour Analysis') self.setItemToolTip(0, 'Airfoil database ' + @@ -871,7 +966,8 @@ def makeToolbox(self): ' on Reynolds number and thermodynamics') self.setItemToolTip(4, 'Compute panel based aerodynamic ' + 'coefficients') - self.setItemToolTip(5, 'Analyze the curvature of the ' + + self.setItemToolTip(5, 'Compute aerodynamic coefficients with SU2') + self.setItemToolTip(6, 'Analyze the curvature of the ' + 'selected airfoil') self.setItemIcon(0, QtGui.QIcon(os.path.join(ICONS_L, 'airfoil.png'))) @@ -879,7 +975,8 @@ def makeToolbox(self): self.setItemIcon(2, QtGui.QIcon(os.path.join(ICONS_L, 'mesh.png'))) self.setItemIcon(3, QtGui.QIcon(os.path.join(ICONS_L, 'Fast delivery.png'))) self.setItemIcon(4, QtGui.QIcon(os.path.join(ICONS_L, 'Fast delivery.png'))) - self.setItemIcon(5, QtGui.QIcon(os.path.join(ICONS_L, 'Pixel editor.png'))) + self.setItemIcon(5, QtGui.QIcon(os.path.join(ICONS_L, 'Fast delivery.png'))) + self.setItemIcon(6, QtGui.QIcon(os.path.join(ICONS_L, 'Pixel editor.png'))) # preselect airfoil database box self.setCurrentIndex(self.tb1) @@ -970,6 +1067,128 @@ def runPanelMethod(self): self.parent.slots.messageBox('No airfoil loaded.') return + def runSU2(self): + """Gui callback to run SU2""" + if self.simulation_initialised: + self.parent.slots.messageBox("Simulation already initialised!") + return + + self.timer = QTimer() # this will be the timer for the live-plot + + sample = { + "aoa": self.aoaSU2.value(), + "Re": self.reynolds_su2.value(), + "Ma": self.freestreamSU2.value(), + "y_TE": 0., # FIXME: get this value from airfoil instead! and then add a user message that it will be subtracted + } + + # TODO: allow user to change the working directory? + paths = {"work_subdir": OUTPUTDATA, "template_dir": SU2CONFIG, + "su2_executable": SU2EXEC} + + self.residuals_canvas.clear() + + # Below is code that would find the latest file in the folder. But there is no guarantee, + # that this was run with the CURRENT settings in the GUI! So we would need to store something, e.g. + # a hash of su2 config and su2 mesh, to ensure that the user actually wanted to continue the existing simulation. + # Or just ask the user with a messagebox? + restart_sol = su2.SU2instance.get_restarts(paths['work_subdir']) + if restart_sol: + answer = self.parent.slots.yesnocancelQuestionBox(f'Use restart-flow file at iteration {restart_sol} (Yes), ' + f'or start a new simulation (No)?') + if answer.name == "Cancel": + return + elif answer.name == "No": + restart_sol = 0 + + settings = { + "restart_file_freq": self.n_restartfreq.value(), + "flowfield_freq": self.n_flowfreq.value(), + } + + deviations = { + "MESH_FILENAME": self.meshfile.text() + } + + self.su2_handler = su2.SU2instance(sample, paths, restart_from=restart_sol, deviations=deviations, settings=settings) + + # from this moment on, the abortSU2-method has everything it needs, so we consider it initialised + self.simulation_initialised = True + + self.su2_handler.run(self.n_cores.value()) + + prev_hist = self.su2_handler.get_history() + + while True: # block until at least one iteration has been written + try: + if len(self.su2_handler.get_history()) != len(prev_hist): + break + except (FileNotFoundError, pd.errors.EmptyDataError): + time.sleep(0.1) + + self.restart_lines = [0] + + self.timer.timeout.connect(self.poll_residuals) + self.timer.start(100) # check every 100 ms + + def poll_residuals(self): # apparently this needs to be a method of the main GUI, otherwise plot doesn't update... + # N.B.: this starts to behave strangely if the user does strange things like manually (re)moving certain + # restart- or history files. In short, this method assumes that the various history.csv are consistent! + sim_status = self.su2_handler.is_done() + if sim_status != "Running": + self.timer.stop() # timeout.disconnect(self.poll_residuals) + self.parent.slots.messageBox(f"SU2 exit code: {sim_status}") + return + + try: + df = self.su2_handler.queue.get_nowait() # the is_done() method above puts something into the queue + except queue.Empty: + return + if len(df) == 0: + return + + df = df.set_index('"Time_Iter"') + ymax = -np.inf + ymin = np.inf + for var, cb in self.residuals_cbs.items(): + if not cb.isChecked(): + continue + kwargs = { + f"{var}_trend": self.n_residuals_movingavg.value(), # if 0, it will not plot it + } + _ = self.residuals_canvas.update_plot(var, df[f'{var}'], **kwargs) + ymax = max(ymax, df[f'{var}'].max()) + ymin = min(ymin, df[f'{var}'].min()) + # note: self.residual_canvas is an object of a custom class + self.residuals_canvas.ax.relim() + self.residuals_canvas.ax.autoscale_view() + # now we draw a dashed, grey, vertical line at every restart + freq = self.n_restartfreq.value() + if freq > 0: + current_iter = df.index.max() + next_restart_index = current_iter // freq + 1 + if next_restart_index > max(self.restart_lines): + _ = self.residuals_canvas.ax.plot([(next_restart_index)*freq]*2, [ymin, ymax], "k--", alpha=0.5)[0] + self.restart_lines.append(next_restart_index) + + self.residuals_canvas.fig.tight_layout() + self.residuals_canvas.draw() + + def abortSU2(self): + if not self.simulation_initialised: + return + if self.su2_handler.process.poll() is None: + if sys.platform != 'win32': + self.su2_handler.process.send_signal(signal.SIGKILL) + else: + self.su2_handler.process.kill() + stdout, stderr = self.su2_handler.process.communicate() + self.parent.slots.messageBox(f"Terminated SU2.") + # TODO: use exit_request flags instead; + # https://stackoverflow.com/questions/323972/is-there-any-way-to-kill-a-thread + self.timer.stop() + self.simulation_initialised = False + def spline_and_refine(self): """Spline and refine airfoil""" @@ -1107,12 +1326,9 @@ def exportMesh(self): filename, extension = os.path.splitext(filename) # add boundary definition attributes to mesh object - self.wind_tunnel.boundary_airfoil = self.lineedit_airfoil.text() - self.wind_tunnel.boundary_inlet = self.lineedit_inlet.text() - self.wind_tunnel.boundary_outlet = self.lineedit_outlet.text() - self.wind_tunnel.boundary_top = self.lineedit_top.text() - self.wind_tunnel.boundary_bottom = self.lineedit_bottom.text() + self.wind_tunnel.boundary_tags_names = {boundary: edit.text() for boundary, edit in self.lineedits_boundarymarkers.items()} + name = "msh.su2" if self.check_FIRE.isChecked(): name = filename + '.flma' Meshing.BlockMesh.writeFLMA(self.wind_tunnel, @@ -1127,6 +1343,8 @@ def exportMesh(self): name = filename + '.vtu' Meshing.BlockMesh.writeVTK_nolib(self.wind_tunnel, name=name) + self.meshfile.setText(name) + def exportContour(self): file_dialog = FileDialog.Dialog() diff --git a/src/su2.py b/src/su2.py new file mode 100644 index 0000000..520cd9a --- /dev/null +++ b/src/su2.py @@ -0,0 +1,212 @@ +import copy +import time + +import numpy as np +import pandas as pd +import math +import datetime as dt +import re +import os +import shutil +import subprocess +from multiprocessing import Queue, Process + +class SU2instance: + """ + This class represents an SU2 instance (or, if run with parallel processing, an MPIexec instance) + Functionality includes: + - modifying a "master" configuration and storing the modified copy as "tmp.cfg" for the desired case ("sample") + - executing SU2 (with multiprocessing support) + - reading the history file. (this also puts the history into a queue object which can be accessed from outside) + - checking whether the simulation has finished by itself (i.e., SU2 has decided that it converged) + - checking whether the simulation is "oscillating" (i.e., SU2 keeps iterating over time with a variation too large + to satisfy the convergence criteria) - this has not been tested recently + """ + def __init__(self, sample, paths, restart_from=False, deviations={}, min_iter=300, k=0, dynamic_monitor_var="CD", settings={}): + # initing attributes + self.process = None + self.dyn_sol_k1 = None # iteration number of the first solution of a dynamic simulation (e.g. min cD) + self.dyn_sol_k2 = None # iteration number of the first solution of a dynamic simulation (e.g. max cD) + self.history = None + self.settings = settings + self.restart_from = restart_from + self.shift_max = int(settings.get('shift_max', 50)) + self.dyn_stab_length = int(settings.get('n_iter_convergence', 50)) + self.max_trend_diff = 2e-4 + self.cauchy_iter = int(settings.get('n_iter_convergence', 50)) + self.p_plot = None + self.maxlen = int(settings.get("restart_file_freq", 100)) # this will raise an error if it does not exist, but that's good # 2025-06-28: changed to default + self.flowfield_freq = int(settings.get("flowfield_freq", 100)) + if self.settings.get('live_plot', True): + self.queue = Queue() + self.queue.cancel_join_thread() # avoid blocking of process + # initing object + for location in ['work_subdir', 'template_dir', 'su2_executable']: # out_dir + if location not in list(paths.keys()): + raise ValueError(f"At least one path ({location}) is not defined!") + self.paths = paths + self.cfg = os.path.join(self.paths['work_subdir'], 'tmp.cfg') + self.logfile = os.path.join(self.paths["work_subdir"], "stdout.log") + self.past_history_files = set() + self.past_history = pd.DataFrame() + self.past_history = self.get_history() + self.dyn_monitor_var = f'"{dynamic_monitor_var}"' + self.deviations = deviations + self.deviations.update({'AOA': sample['aoa'] - np.arctan2(sample['y_TE'], 1) * 180 / np.pi, # TODO: check sign!! + 'MACH_NUMBER': sample['Ma'], + 'REYNOLDS_NUMBER': sample['Re'], + 'OUTPUT_WRT_FREQ': f'({self.maxlen}, {self.flowfield_freq})', + 'CONV_WINDOW_CAUCHY_ELEMS': self.cauchy_iter}) + + if self.restart_from > 0: + self.deviations.update({'RESTART_SOL': 'YES', + 'RESTART_ITER': self.restart_from+1, + 'WINDOW_START_ITER': max(min_iter, self.restart_from+1)}) # TODO: check for off-by-one + self.prepare() + + def prepare(self, template=None): # modify config file + if template is None: + template = os.path.join(self.paths['template_dir'], 'su2_master.cfg') + with open(template) as f_src: + with open(self.cfg, "w") as f_dest: + while True: + try: + line = next(f_src) + f_dest.write(self._modify_cfg(line)) + except StopIteration: + break + + def _modify_cfg(self, line): + for parameter, value in self.deviations.items(): + pattern = "^"+parameter+"= " + match = re.search(pattern, line) + if match: + return f'{parameter}= {value}\n' + return line + + def run(self, cores=1): + self.exec(self.paths["su2_executable"], self.cfg, n_cores=cores) + + def exec(self, executable, *args, **kwargs): + n_cores = kwargs.get('n_cores', 1) + if n_cores == 1: + with open(self.logfile, "w") as outfile: + self.process = subprocess.Popen([executable, *args], cwd=self.paths["work_subdir"], stdout=outfile) # Popen: , stdout=subprocess.PIPE, shell=True, preexec_fn=os.setsid) + else: + with open(self.logfile, "w") as outfile: + self.process = subprocess.Popen(['mpiexec', '-n', str(n_cores), executable, *args], cwd=self.paths["work_subdir"], stdout=outfile) # Popen: , stdout=subprocess.PIPE, shell=True, preexec_fn=os.setsid) + # print("SU2 finished.") + + def kill(self): + # os.killpg(os.getpgid(self.process.pid), signal.SIGTERM) + pass + + @staticmethod + def get_restarts(workdir): + restarts = [] + for root, dirs, files in os.walk(workdir): + for file in files: + match = re.search("restart_flow_(\d{5}).dat", file) + if match: + restarts.append(int(match.groups()[0])) + for i in range(len(restarts)-1): + if restarts[-1-i] == restarts[-2-i]+1: + return restarts[-1-i] + return 0 + + def get_history(self): + data_by_files = {} + for root, dirs, files in os.walk(self.paths['work_subdir']): # su2 creates a new history file when continuing from a restart-file + for file in files: + ffull = os.path.join(root, file) + if 'history' in file: + df_temp = pd.read_csv(ffull, sep='\s*,\s*', engine='python') + data_by_files[ffull] = df_temp.groupby('"Time_Iter"').tail(1) # we take the last time iteration + full_history = pd.DataFrame() + fs = sorted(data_by_files.keys()) + + for i in range(len(fs)): + # first, we check the CURRENT file, whether its value is greater than our restart_from; then we ignore it + match = re.search("history_(\d{5}).csv", fs[i]) + if match and int(match.groups()[0]) - 1 > self.restart_from: + continue # the "if match" implies that this never happens for the first history file, which works as intended + + try: + # if there are previous restart files, read the current file only up until where a more recent file starts + # already has data (this can occur due to restarts) + match = re.search("history_(\d{5}).csv", fs[i+1]) + # ^ the number in the NEXT file's name (fs[i+1]) is the number where THIS file (fs[i]) should stop + if match: + maxlen = int(match.groups()[0])-1 + else: + maxlen = -1 # -1 would be the last, but most likely the last line we read is not the last INNER iter of the last time iter, so we skip one + except IndexError: + maxlen = -1 + + df2 = data_by_files[fs[i]] + if maxlen > 0: + df2_trunc = df2[df2['"Time_Iter"'] <= maxlen] + else: + df2_trunc = df2.iloc[:maxlen] + + if full_history.isnull().values.all(): + full_history = df2_trunc + elif not df2_trunc.isnull().values.all(): + full_history = pd.concat([full_history, df2_trunc]) + # else, we just keep full_history (i.e. do nothing) + + self.queue.put(full_history) + return full_history + + def is_done(self): + full_history = self.get_history() + if dt.datetime.now().timestamp() - os.path.getmtime(self.logfile) > 10: # check the file has not been modified in the last ten seconds + with open(self.logfile) as fid: + while True: # check the file for exit messages + try: + line = next(fid) + if "Error Exit" in line: + return "Error Exit" + if "Exit Success" in line: # TODO: find actual name + return "Exit Success" + except StopIteration: + break # in this case we still need to check whether it has reached dynamic stability + period = self._check_dynamically_stable(full_history, column_name=self.dyn_monitor_var) + if period: + val = f"Dynamically Stable. Period: {period}" + else: + val = "Running" + return val + + def _check_dynamically_stable(self, df, column_name='"CD"', ): # TODO: tune + # Extract the specified column + shift_max = self.shift_max + l = self.dyn_stab_length + max_trend_diff = self.max_trend_diff + if len(df) < l+shift_max: + return False + column_values = df[column_name].values # get raw values + periods = [] + trend = pd.DataFrame({'data': column_values})['data'].ewm(halflife=shift_max).mean().values[-shift_max-l:-shift_max] + for val in trend: + if abs(trend[-1] - val)/l > max_trend_diff: + return False # if any value of the drift of the oscillation is outside the limits => not converged + # all last l points are within max_trend_diff => we check whether there is an oscillatory behaviour + # we do this by checking whether the signals "repetition period" (found by autocorrelation) ... + # ... is INDEPENDENT of where we look at the signal, so we shift through shift_max + for shift in range(1, shift_max): + subvals = column_values[-l-shift:-shift] # get last l values + mn = np.mean(subvals) # find mean, so we can offset by it + shifted = subvals - mn # the offset ensures that the values vary around 0, otherwise autocorrel. gives bad values + c = np.correlate(shifted, shifted, 'same')[-l//2:] # autocorrelate the data # for some reason simple division does not work... but it makes sense anyway, just could cause the period to be slightly underestimated + inflection = np.diff(np.sign(np.diff(c))) # find the second-order differences + peaks = (inflection < 0).nonzero()[0] + 1 # find where they are negative + try: + delay = peaks[c[peaks].argmax()] # find the delay where the highest peak occurs (peak of subsequent periods will be lower due to autofill of zeros) + periods.append(delay) + except (ValueError, IndexError): + return False + if len(list(set(periods))) == 1: # period no longer changes + return periods[0] + else: + return False