Παραλληλοποιήσιμες πολυπλεγματικές τεχνικές σε επιλυτές Navier-Stokes για τη μοντελοποίηση ασυμπίεστων ροών

Abstract

The Navier-Stokes equations, that govern the motion of an incompressible or compressible uid, were introduced more than a hundred years ago. Their nonlinearity leads to signicant diculties, sometimes insurmountable, when trying to nd an analytical solution. In order to obtain an exact solution, specic geometries and initial/boundary conditions must be considered. In more complex situations, such as applications in uid structure interaction, low speed aerodynamics, biomechanics etc., the application of numerical methods may providereliablesolutions forthe Navier-Stokesequations. Furthermore, the necessity for a fast high-resolution numerical solver for the incompressible Navier-Stokes equations emerges from real-life simulations, such as ows over hydrofoils, wind-turbine blades, and aircraft wings during takeo and landing. Thus, the objective of this thesis is to develop, explicateanddemonstratetheperformanceofahighlyecientow-solver,whichexploits the computational power of architectures with computing accelerators. In the rst part of this thesis, and after the preliminaries, an elliptic PDE multigrid solver is developed and demonstrated. The solver is capable of handling highly anisotropic 2D Boundary Value Problems (BVPs) and is based on high-order cell-centered Finite Difference Compact schemes and Multigrid techniques. Compact schemes provide a representation of the shorter length scales, when applied to problems with a range of spatial scales, compared to traditional nite dierence approximations. An improvement of the proposed method is the treatment of PDE boundary conditions. Boundary closure formulas for Dirichlet, Neumann, Robin or mixed-type boundary conditions applied to the physical boundary are derived and tested herein for several simulation problems. Moreover, novel multigrid components for cell-centered discretization are being constructed, which could also be generalized to three dimensional problems in a straightforward manner. In particular, the new intergrid operators involve less non-zero entries than common operators in cell-centered grids, preserving at the same time the accuracy of high-order operators. Next, some theoretical convergence results for the multigrid solver using Local Fourier Analysis (LFA)are given in order to improve its multigrid convergence factors. The analysis focuseson both the relaxation method used within the multigrid, as well as in the remaining components of the multigrid method. The proposed multigrid elliptic solver is evaluated on two classical BVPs, so that the fourth-order accuracy of the solver, as well as the boundary treatments, could be validated. The multigrid convergence rates for every acknowledged transfer operator, along with every novel one, are evaluated as to determine the convergence behavior of the corresponding multigrid solver. These results are also compared to the theoretical convergence factors based on the LFA. The concordance of the analytical and numerical results is acceptable. A comparison between the calculated convergence rates to the corresponding values, obtained from the second-order compact scheme, indicates the superiority of the high-order compact scheme. In addition, these convergence rates are as good as the factors one obtains from the vertex-centered case. It is noted, that the cell-centered multigrid numerical analysis is presented for a high-order scheme, opposed to earlier studies. Along with the investigation on the improvement of multigrid techniques for the high-order scheme, iterative methods are also tested for the resulting sparse linear system and compared with the multigrid numerical solver. The comparison results indicate the necessity of a multigrid solver for this kind of problems, as it is proved to be hundreds of times faster than the optimal iterative method when ne discretizations are used. The proposed spatial discretization solver is incorporated in an eective Navier-Stokes solver capable of handling highly anisotropic ow problems. The solver is based on the pressure-velocity coupling and uses fourth-order compact schemes for discretizing each spatial dimension, formulated on a staggered grid arrangement. The temporal discretization is carried out by a fourth-order Runge-Kutta (RK4) method. Incompressibility is enforced at each time step using a global pressure correction method solving a Poisson-type PDE. In this method, the multigrid solver is being used within each stage of the RK4 method to compute the pressure correction. The spatial and temporal fourth-order accuracy of this Navier-Stokes solver are validated for a set of steady and unsteady classical test problems. Further, the performance results indicate that multigrid accelerates the solution procedure more than 10 times, comparing with other solvers in the literature. The numerical study of the sequential Navier-Stokes solver evince that, in cases of grid sizes up to 256×256, the incorporation of the multigrid scheme handles the increasing execution time moderately. However, in case of ner grid sizes, the computational cost becomes intolerable despite the High convergence rates of the Multigrid method. It is noted that the most time-consuming part of the solver is the pressure correction procedure. This time restriction gives motive to redesign and develop an ecient parallel multigrid based Navier-Stokes algorithm, to exploit the benets of modern parallel computer architectureswith accelerators. In order to increase parallelism at each computing phase of the algorithm, the horizontal red-black coloring scheme for grid nodes is chosen. The Block Cyclic Reduction method is also applied for the solution of the arising linear sub-systems, without modifying the multigrid cycling nature in the algorithm. This enables the execution of the entire simulation in the acceleration device, minimizing the communication cost between memory units. In addition, the proposed parallelization exploits the block structure of the coecient matrix, minimizing data storage and increasing again the parallelism. The solver is implemented and examined on three parallel machines with dierent type of accelerator devices. The realization is developed using the OpenACC and OpenMP APIs. The eect of several multigrid components on modern and legacy acceleration architectures is investigated Application’s performance investigation demonstrates that the proposed parallel multigrid solver accomplishes an acceleration of more than ten times over the sequential solver and more than four times over multi-coreCPU-only realizations. In case of highly anisotropic problems, the parallel semi-coarsening multigrid solver is preferred to the full-coarsening one, as the division of labor by the accelerator device provides faster computational rates, in case of non-uniform discretizations. The proposed numerical algorithm can be easily extended for the case of three dimensional ow problems, on curvilinear coordinates, expanding the applicability of the current methodology. Furthermore, the proposed numerical methodology can be applied for solving comparable problems, e.g. Maxwells equations. The design of a parallel algorithm for the utilization of a Heterogeneous Multi-Accelerator architecture using the MPI, OpenMP and OpenACC APIs, is considered to be a promising improvement.

Οι Navier-Stokes ρευστοδυναμικές εξισώσεις για ασυμπίεστες ή συμπιεστές ροές είναι γνωστές τα τελευταία εκατό χρόνια. Ο μη-γραμμικός χαρακτήρας τους καθιστά αδύνατη την αναλυτική επίλυση των προβλημάτων μοντελοποίησης, τα οποία παρουσιάζονται για παράδειγμα στις ερευνητικές περιοχές της Αεροδιαστημικής, Βιοϊατρικής και Εμβιομηχανικής. Οι αναλυτικές λύσεις εξισώσεων κίνησης είναι γνωστές μόνο για ένα πολύ περιορισμένο αριθμό στοιχειωδών ροών και για απλές γεωμετρίες με συγκεκριμένες αρχικές/οριακές συνθήκες. Έτσι είναι αναγκαία η χρήση αριθμητικών μεθόδων επίλυσης για την προσέγγιση των λύσεων τους. Οι αριθμητικές αυτές μέθοδοι υλοποιούνται σε σύγχρονα υπολογιστικά περιβάλλοντα ιδιαίτερα για ρεαλιστικές εφαρμογές. Από τις αρχές του 1980 υπήρξε μια μεγάλη ανάπτυξη αριθμητικών μεθόδων επίλυσης προβλημάτων Ρευστοδυναμικής (Computational Fluid Dynamics-CFD). Παρόλαυτα, ακόμα και στις μέρες μας, η ανάπτυξη υψηλής ακρίβειας αριθμητικών σχημάτων επίλυσης των εξισώσεων Navier-Stokes για ασυμπίεστες ροές είναι αναγκαία για την προσομοίωση μεγάλου ενδιαφέροντος προβλημάτων. Για παράδειγμα, η μοντελοποίηση της ροής γύρω από υδροπτέρυγα, πτέρυγες ανεμογεννητριών και αεροσκαφών, κατά την διαδικασία απογείωσης και προσγείωσης, της ροής αίματος και ροής του αέρα γύρω από πουλιά και έντομα κατά το πέταγμα τους. Τέτοια προβλήματα ροών μπορούν να αντιμετωπιστούν μόνο με την χρήση υψηλής ακρίβειας αριθμητικών μεθόδων. Όμως, η υψηλή ακρίβεια που απαιτείται αυξάνει συνήθως εκθετικά τον υπολογιστικό χρόνο επίλυσης. Ο συνδυασμός υψηλής ποιότητας προσομοιώσεων ασυμπίεστων ροών σε εύλογο χρονικό διάστημα απεικόνισης των ταχυτήτων ροής, αποτέλεσε το κίνητρο για την συγγραφή της παρούσας διδακτορικής διατριβής. Αντικείμενο της αποτελεί η ανάπτυξη, ανάλυση και εφαρμογή ενός αποτελεσματικού επιλυτή ασυμπίεστων μεταβατικών τυρβωδών ροών, ο οποίος θα εκμεταλλεύεται τις σύγχρονες πολυεπεξεργαστικές υπολογιστικές αρχιτεκτονικές με επιταχυντές υπολογισμών. Για την αποδοτική αριθμητική επίλυση προβλημάτων ασυμπίεστων ροών έχουν επικρατήσει μέθοδοι διόρθωσης της πίεσης και της ψευδοσυμπιεστότητας, με την δεύτερη να μην ενδείκνυται για μη μόνιμες ροές. Η αριθμητική μέθοδος που θα αναπτυχθεί θα βασιστεί στην εισαγωγή της πίεσης στην εξίσωση συνέχειας (μέθοδος διόρθωσης πίεσης) ενσωματώνοντας την τεχνική πολλαπλών πλεγμάτων (multigrid technique) στον αριθμητικό επιλυτή διακριτοποίησης. Ειδικότερα, γίνεται χρήση υψηλής τάξης ακρίβειας συμπαγών αριθμητικών σχημάτων πεπερασμένων διαφορών για μετατοπισμένα πλέγματα διακριτοποίησης (staggered grids). Τεχνικές παραλληλοποίησης της επαναληπτικής διαδικασίας επίλυσης του παραγόμενου αραιού και γενικού γραμμικού συστήματος εφαρμόζονται για την αποδοτική υλοποίηση της μεθόδου. Τα συμπαγή σχήματα πεπερασμένωνδιαφορών (compact or implicit schemes) έχουν το πλεονέκτημα, ότι επιτυγχάνουν σε μικρό εύρος πλέγματος μεγαλύτερη τάξη ακρίβειας, σε σχέση με τα κλασικά σχήματα πεπερασμένων διαφορών. Η ιδιότητα αυτή είναι χρήσιμη για την εφαρμογή του σχήματος κοντά στα όρια του υπολογιστικού χωρίου. Μία καινοτομία της προτεινόμενης μεθόδου είναι η επιβολή οποιοδήποτε είδους συνοριακών συνθηκών (Dirichlet, Neumann, Robin and mixed συνοριακές συνθήκες) στον αριθμητικό επιλυτή. Επίσης, σε αντίθεση με τα κλασικά σχήματα πεπερασμένων διαφορών, η χρήση συμπαγών σχημάτων ελαχιστοποιεί το ποσοστό της αριθμητικής διάχυσης, ελαχιστοποιώντας τη διαφορά μεταξύ πραγματικών και αριθμητικών κυματικών αριθμών σε όλες τις διαδιδόμενες στο πλέγμα συχνότητες. Η ελαχιστοποίηση της αριθμητικής διάχυσης είναι επιθυμητή ιδιότητα του αριθμητικού σχήματος, ιδιαίτερα στις περιπτώσεις προσομοίωσης αριθμητικής μετάδοσης κυμάτων (π.χ. αερο-ακουστική και ηλεκτρομαγνητισμός), αλλά και σε προβλήματα με ευρύ συχνοτικό φάσμα ή σε προβλήματα που πραγματοποιούνται ταυτόχρονα σε ευρύ φάρμα χωρικών κλιμάκων.Στην αρχή της παρούσας διατριβής, περιγράφεται η κατασκευή ένος αριθμητικού επιλυτή για ελλειπτικού τύπου ΜΔΕ στις δύο χωρικές διαστάσεις, ο οποίος είναι ικανός να αντιμετωπίσει Προβλήματα Συνοριακών Τιμών με ανισοτροπίες. Το αριθμητικό σχήμα είναι κατάλληλο για γεωμετρικές διαμερίσεις με τους αγνώστους κατανεμημένους στα κέντρα των υπολογιστικών κελιών (cell-centered grids). Επιπλέον, για την εφαρμογή της τεχνικής Πολυπλέγματος κατασκευάζονται πρωτότυποι υψηλής ακρίβειας τελεστές μεταφοράς της πληροφορίας από το πυκνό στο αραιό πλέγμα και αντίστροφα. Ειδικότερα, οι τελεστές αυτοί εμπλέκουν λιγότερες αριθμητικές πράξεις σε σχέση με τους συνηθισμένους για cell-centered διακριτοποιήσεις. Στη συνέχεια, χρησιμοποιώντας ανάλυση Fourier σε τοπικό επίπεδο, δηλαδή σε επίπεδο υπολογιστικού κελίου, παράγονται θεωρητικά αποτελέσματα σύγκλισης της τεχνικής πολυπλέγματος για το προτεινόμενο αριθμητικό σχήμα. Η ανάλυση αυτή επικεντρώνεται τόσο στην επιλογή της μεθόδου χαλάρωσης, όσο και στα υπόλοιπα μέρη που συνθέτουν την τεχνική Πολυπλέγματος, με σκοπό τη βελτιστοποίηση των ρυθμών σύγκλισης του επιλυτή. Ο προτεινόμενος αριθμητικός επιλυτής μελετάται σε δύο προβλήματα δοκιμής, ώστε να επιβεβαιωθεί η τάξη ακρίβειας για κάθε επιλογή συνοριακών συνθηκών. Επίσης, υπολογίζονται οι ρυθμοί σύγκλισης της Πολυπλεγματικής μεθόδου για όλους τους γνωστούς τελεστές μεταφοράς, συμπεριλαμβανομένου και των νέων προτεινόμενων τελεστών. Οι θεωρητικοί ρυθμοί σύγκλισης βρίσκονται σε συμφωνία με τα αντίστοιχα αριθμητικά αποτελέσματα, παρέχοντας σχετική ασφάλεια για την επιλογή των βασικών συστατικών της τεχνικής Πολυπλέγματος. Επιπρόσθετα, οι παρατηρηθείσες ρυθμοί σύγκλισης είναι σημαντικά καλύτεροι σε σχέση με αυτούς των αριθμητικά συμπαγών σχημάτων δεύτερης τάξης ακρίβειας, ενώ ταυτόχρονα μπορούν να συγκριθούν με τους ρυθμούς σύγκλισης που μπορεί κάποιος να επιτύχει σε διακριτοποιήσεις, όταν οι άγνωστοι κατανέμονται στις κορυφές των υπολογιστικών κελιών (vertex-centered grids). Το νέο αυτό σχήμα της Πολυπλεγματικής μεθόδου εφαρμόζεται με επιτυχία σε cell-centered διακριτοποιήσεις για υψηλής ακρίβειας μεθόδους. Ταυτόχρονα με την προσπάθεια βελτιστοποίησης της τεχνικής Πολυπλέγματος του αριθμητικού σχήματος, μελετάται και η διαδικασία επίλυσης του παραγόμενου γραμμικού συστήματος με την συγκριτική εφαρμογή των πιο δημοφιλών επαναληπτικών μεθόδων. Τα συγκριτικά αποτελέσματα που παρουσιάζονται αναδεικνύουν την αναγκαιότητα χρήσης της τεχνικής Πολυπλέγματος για προβλήματα του συγκεκριμένου είδους, αφού παρατηρήθηκε μία αξιόλογη μείωση στο χρόνο επίλυσης του γραμμικού συστήματος, σε σύγκριση με την αποδοτικότερη κλασική επαναληπτική μέθοδο επίλυσης. Ο προτεινόμενος αριθμητικός επιλυτής των εξισώσεων Navier-Stokes μπορεί να αντιμετωπίσει ικανοποιητικά προβλήματα μη-σταθερών ασυμπίεστων ροών, τα οποία παρουσιάζουν διαφορετική κατανομή της φυσικής ποσότητας σε διαφορετικές διευθύνσεις (ανισοτροπίες). Για την διακριτοποίηση των χωρικών μεταβλητών χρησιμοποιούνται συμπαγή σχήματα πεπερασμένων διαφορών τέταρτης τάξης, ενώ για αυτή του χρόνου η κλασική αριθμητική μέθοδος Runge-Kutta τέταρτης τάξης. Η ασυμπιεστότητα του ρευστού εφαρμόζεται σε κάθε χρονικό βήμα του αριθμητικού χρονικού σχήματος μέσω της επίλυσης ενός ελλειπτικού τύπου προβλήματος με σκοπό τη διόρθωση της πίεσης και επακόλουθα του πεδίου ταχυτήτων. Η αποδοτικότητα και η ακρίβεια του αριθμητικού επιλυτή Navier-Stokes επαληθεύεται για προβλήματα δοκιμών σταθερών και μεταβαλλόμενων χρονικά ροών. Επιπρόσθετα, ο προτεινόμενος επιλυτής στα προβλήματα δοκιμής βελτίωσε σημαντικά τους χρόνους εκτέλεσης της επίλυσης των Navier-Stokes. Στις περισσότερες περιπτώσεις, η νέα μέθοδος με τη χρήση της τεχνικής Πολυπλέγματος επιτάχυνε τη διαδικασία επίλυσης πάνω από δέκα φορές. Οι μετρήσεις απόδοσης του αλγορίθμου αναδεικνύουν ότι για προβλήματα μεγέθους διακριτοποίησης μέχρι 256×256 υπολογιστικών κελιών, η ενσωμάτωση της τεχνικής Πολυπλέγματος αντιμετωπίζει με επιτυχία τον αυξανόμενο χρόνο εκτέλεσης που απαιτείται. Ωστόσο, για μεγαλύτερες διακριτοποιήσεις του χωρίου, οι οποίες χρειάζονται μικρότερα χρονικά βήματα διακριτοποίησης, το υπολογιστικό κόστος αποτελεί έναν αποτρεπτικό παράγοντα για την επίλυση ρεαλιστικών προβλημάτων. Επίσης, υπολογίστηκε ότι η πιο χρονοβόρα διαδικασία του αλγορίθμου είναι αυτή της διόρθωσης της πίεσης. Αυτό υπήρξε το κίνητρο για τον επανασχεδιασμό του αλγορίθμου επίλυσης των Navier-Stokes και την ανάπτυξη ενός αποδοτικού παράλληλου αλγορίθμου για σύγχρονες υπολογιστικές αρχιτεκτονικές. Για την αύξηση της παραλληλοποίησης του αλγορίθμου χρησιμοποιήθηκαν τεχνικές όπως: 1) αναδιάταξη των κόμβων του πλέγματος σύμφωνα με το σχήμα ομαδοποίησης κόκκινου-μαύρου(horizontal red-black ordering) με σκοπό την ανεξαρτητοποίηση ομάδων αγνώστων και 2) η μέθοδος Block Cyclic Reduction για την επιτάχυνση της επίλυσης των ομαδοποιημένων γραμμικών υπο-συστημάτων. Οι τεχνικές αυτές εφαρμόστηκαν για να μην αυξηθούν τα βήματα εκτέλεσης του σειριακού αλγορίθμου. Η κατάλληλη οργάνωση των υπολογισμών έκανε εφικτή την εκτέλεση όλης της υπολογιστικής προσομοίωσης σε συσκευή επιτάχυνσης υπολογισμών, με αποτέλεσμα την μείωση στο ελάχιστο του κόστους επικοινωνίας μεταξύ της κύριας μνήμης του υπολογιστικού συστήματος και της μνήμης του επιταχυντή. Επιπλέον, η προτεινόμενη διαδικασία παραλληλοποίησης εκμεταλλεύεται την επαναλαμβανόμενη block δομή του πίνακα των συντελεστών, ελαχιστοποιώντας τις απαιτήσεις σε αποθηκευτικό χώρο. Η υλοποίηση του παράλληλου αλγορίθμου πραγματοποιήθηκε σε περιβάλλον ανάπτυξης του προτύπου OpenACC και δοκιμάστηκε σε τρεις υπολογιστικές αρχιτεκτονικές κοινής μνήμης με διαφορετικών τύπων επιταχυντές. Τα αποτελέσματα απόδοσης του προτεινόμενου αλγορίθμου είναι πολύ ενθαρρυντικά, αφού ο παράλληλος Navier-Stokes επιλυτής επιτυγχάνει μια επιτάχυνση μεγαλύτερη από 10 φορές σε σύγκριση με την σειριακή υλοποίηση του και 4 φορές σε σχέση με τηνυλοποίηση σε κοινής μνήμης πολυεπεξεργαστικές αρχιτεκτονικές χωρίς όμως την χρήση επιταχυντών. Επίσης, η εφαρμογή του σε προβλήματα με μεγάλες ανισοτροπίες ανέδειξε την υπεροχή της Πολυπλεγματικής μεθόδου, που χρησιμοποιεί την τεχνική της ημι-αραίωσης του πλέγματος στην διεύθυνση που παρουσιάζεται η ανισοτροπία σε σχέση με την τεχνική της πλήρους-αραίωσης του και στις δύο διαστάσεις. Αν και η παρούσα ερευνητική διαδικασία επικεντρώνεται σε προβλήματα ορισμένα σε διδιάστατα καρτεσιανά πλέγματα, η επέκτασή της σε τρεις διαστάσεις ή σε καμπυλόγραμμες συντεταγμένες είναι άμεση. Αυτό συμβαίνει διότι είναι εύκολη η επέκταση όλων των διαδικασιών του παραγόμενου αλγορίθμου, καθώς και της επιτάχυνσης μέσω της τεχνικής πολυπλέγματος. Επίσης, η μέθοδος ενδείκνυται για μεταβατικές τυρβώδεις ασυμπίεστες ροές με την μέθοδο προσομοιώσεις τύρβης μεγάλων δινών (large eddy simulation LES), ενώ με τη χρήση έμμεσων αριθμητικών χρονικών μεθόδων και μοντέλων τύρβης μπορούν να δημιουργηθούν προσομοιώσεις τυρβωδών ροών υψηλής πιστότητας.