A parallel algorithm for deformable contact problems

In the field of nonlinear computational solid mechanics, contact problems deal with the deformation of separate bodies which interact when they come in touch. Usually, these problems are formulated as constrained minimization problems which may be solved using optimization techniques such as penalty method, Lagrange multipliers, Augmented Lagrangian method, etc. This classical approach is based on node connectivities between the contacting bodies. These connectivities are created through the construction of contact elements introduced for the discretization of the contact interface, which incorporate the contact constraints in the global weak form. These methods are well known and widely used in the resolution of contact problems in engineering and science. As parallel computing platforms are nowadays widely available, solving large engineering problems on high performance computers is a real possibility for any engineer or researcher. Due to the memory and compute power that contact problems require and consume, they are good candidates for parallel computation. Industrial and scientific realistic contact problems involve different physical domains and a large number of degrees of freedom, so algorithms designed to run efficiently in high performance computers are needed. Nevertheless, the parallelization of the numerical solution methods that arises from the classical optimization techniques and discretization approaches presents some drawbacks which must be considered. Mainly, for general contact cases where sliding occurs, the introduction of contact elements requires the update of the mesh graph in a fixed number of time steps. From the point of view of the domain decomposition method for parallel resolution of numerical problems this is a major drawback due to its computational expensiveness, since dynamic repartitioning must be done to redistribute the updated mesh graph to the different processors. On the other hand, some of the optimization techniques modify dynamically the number of degrees of freedom in the problem, by introducing Lagrange multipliers as unknowns. In this work we introduce a Dirichlet-Neumann type parallel algorithm for the numerical solution of nonlinear frictional contact problems, putting a strong focus on its computational implementation. Among its main characteristics it can be highlighted that there is no need to update the mesh graph during the simulation, as no contact elements are used. Also, no additional degrees of freedom are introduced into the system, since no Lagrange multipliers are required. In this algorithm the bodies in contact are treated separately, in a segregated way. The coupling between the contacting bodies is performed through boundary conditions transfer at the contact zone. From a computational point of view, this feature allows to use a multi-code approach. Furthermore, the algorithm can be interpreted as a black-box method as it solves each body separately even with different computational codes. In addition, the contact algorithm proposed in this thesis can also be formulated as a general fixed-point solver for the solution of interface problems. This generalization gives us the theoretical basis to extrapolate and implement numerical techniques that were already developed and widely tested in the field of fluid-structure interaction (FSI) problems, especially those related to convergence ensurance and acceleration. We describe the parallel implementation of the proposed algorithm and analyze its parallel behaviour and performance in both validation and realistic test cases executed in HPC machines using several processors.En el ámbito de la mecánica de contacto computacional, los problemas de contacto tratan con la deformación que sufren cuerpos separados cuando interactúan entre ellos. Comunmente, estos problemas son formulados como problemas de minimización con restricciones, que pueden ser resueltos utilizando técnicas de optimización como la penalización, los multiplicadores de Lagrange, el Lagrangiano Aumentado, etc. Este enfoque clásico está basado en la conectividad de nodos entre los cuerpos, que se realiza a través de la construcción de los elementos de contacto que surgen de la discretización de la interfaz. Estos elementos incorporan las restricciones de contacto en forma débil. Debido al consumo de memoria y a los requerimientos de potencia de cálculo que los problemas de contacto requieren, resultan ser muy buenos candidatos para su paralelización computacional. Sin embargo, tanto la paralelización de los métodos numéricos que surgen de las técnicas clásicas de optimización como los distintos enfoques para su discretización, presentan algunas desventajas que deben ser consideradas. Por un lado, el principal problema aparece ya que en los casos más generales de la mecánica de contacto ocurre un deslizamiento entre cuerpos. Por este motivo, la introducción de los elementos de contacto vuelve necesaria una actualización del grafo de la malla cada cierto número de pasos de tiempo. Desde el punto de vista del método de descomposición de dominios utilizado en la resolución paralela de problemas numéricos, esto es una gran desventaja debidoa su coste computacional. En estos casos, un reparticionamiento dinámico debe ser realizado para redistribuir el grafo actualizado de la malla entre los diferentes procesadores. Por otro lado, algunas técnicas de optimización modifican dinámicamente el número de grados de libertad del problema al introducir multiplicadores de Lagrange como incógnitas. En este trabajo presentamos un algoritmo paralelo del tipo Dirichlet-Neumann para la resolución numérica de problemas de contacto no lineales con fricción, poniendo un especial énfasis en su implementación computacional. Entre sus principales características se puede destacar que no hay necesidad de actualizar el grafo de la malla durante la simulación, ya que en este algoritmo no se utilizan elementos de contacto. Adicionalmente, ningún grado de libertad extra es introducido al sistema, ya que los multiplicadores de Lagrange no son requeridos. En este algoritmo los cuerpos en contacto son tratados de forma separada, de una manera segregada. El acople entre estos cuerpos es realizado a través del intercambio de condiciones de contorno en la interfaz de contacto. Desde un punto de vista computacional, esta característica permite el uso de un enfoque multi-código. Además, este algoritmo puede ser interpretado como un método del tipo black-box ya que permite resolver cada cuerpo por separado, aún utilizando distintos códigos computacionales. Adicionalmente, el algoritmo de contacto propuesto en esta tesis puede ser formulado como un esquema de resolución de punto fijo, empleado de forma general en la solución de problemas de interfaz. Esta generalización permite extrapolar técnicas numéricas ya utilizadas en los problemas de interacción fluido-estructura e implementarlas en la mecánica de contacto, en especial aquellas relacionadas con el aseguramiento y aceleración de la convergencia. En este trabajo describimos la implementación paralela del algoritmo propuesto y analizamos su comportamiento y performance paralela tanto en casos de validación como reales, ejecutados en computadores de alta performance utilizando varios procesadores.Postprint (published version

Recent EUROfusion Achievements in Support of Computationally Demanding Multiscale Fusion Physics Simulations and Integrated Modeling

Integrated modeling (IM) of present experiments and future tokamak reactors requires the provision of computational resources and numerical tools capable of simulating multiscale spatial phenomena as well as fast transient events and relatively slow plasma evolution within a reasonably short computational time. Recent progress in the implementation of the new computational resources for fusion applications in Europe based on modern supercomputer technologies (supercomputer MARCONI-FUSION), in the optimization and speedup of the EU fusion-related first-principle codes, and in the development of a basis for physics codes/modules integration into a centrally maintained suite of IM tools achieved within the EUROfusion Consortium is presented. Physics phenomena that can now be reasonably modelled in various areas (core turbulence and magnetic reconnection, edge and scrape-off layer physics, radio-frequency heating and current drive, magnetohydrodynamic model, reflectometry simulations) following successful code optimizations and parallelization are briefly described. Development activities in support to IM are summarized. They include support to (1) the local deployment of the IM infrastructure and access to experimental data at various host sites, (2) the management of releases for sophisticated IM workflows involving a large number of components, and (3) the performance optimization of complex IM workflows.This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014 to 2018 under grant agreement 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission or ITER.Peer ReviewedPostprint (published version

The cooperative parallel: A discussion about run-time schedulers for nested parallelism

Nested parallelism is a well-known parallelization strategy to exploit irregular parallelism in HPC applications. This strategy also fits in critical real-time embedded systems, composed of a set of concurrent functionalities. In this case, nested parallelism can be used to further exploit the parallelism of each functionality. However, current run-time implementations of nested parallelism can produce inefficiencies and load imbalance. Moreover, in critical real-time embedded systems, it may lead to incorrect executions due to, for instance, a work non-conserving scheduler. In both cases, the reason is that the teams of OpenMP threads are a black-box for the scheduler, i.e., the scheduler that assigns OpenMP threads and tasks to the set of available computing resources is agnostic to the internal execution of each team. This paper proposes a new run-time scheduler that considers dynamic information of the OpenMP threads and tasks running within several concurrent teams, i.e., concurrent parallel regions. This information may include the existence of OpenMP threads waiting in a barrier and the priority of tasks ready to execute. By making the concurrent parallel regions to cooperate, the shared computing resources can be better controlled and a work conserving and priority driven scheduler can be guaranteed.Peer ReviewedPostprint (author's final draft

multiRegionFoam -- A Unified Multiphysics Framework for Multi-Region Coupled Continuum-Physical Problems

This paper presents a unified framework, called multiRegionFoam, for solving multiphysics problems of the multi-region coupling type within OpenFOAM (FOAM-extend). This framework is intended to supersede the existing solver with the same name. The design of the new framework is modular, allowing users to assemble a multiphysics problem region-by-region and coupling conditions interface-by-interface. The present approach allows users to choose between deploying either monolithic or partitioned interface coupling for each individual transport equation. The formulation of boundary conditions is generalised in the sense that their implementation is based on the mathematical jump/transmission conditions in the most general form for tensors of any rank. The present contribution focuses on the underlying mathematical model for these types of multiphysics problems, as well as on the software design and resulting code structure that enable a flexible and modular approach. Finally, deployment for different multi-region coupling cases is demonstrated, including conjugate heat, multiphase flows and fuel-cells

High-Performance Computing: Dos and Don’ts

Computational fluid dynamics (CFD) is the main field of computational mechanics that has historically benefited from advances in high-performance computing. High-performance computing involves several techniques to make a simulation efficient and fast, such as distributed memory parallelism, shared memory parallelism, vectorization, memory access optimizations, etc. As an introduction, we present the anatomy of supercomputers, with special emphasis on HPC aspects relevant to CFD. Then, we develop some of the HPC concepts and numerical techniques applied to the complete CFD simulation framework: from preprocess (meshing) to postprocess (visualization) through the simulation itself (assembly and iterative solvers)

Fluid–structure interaction simulations outperform computational fluid dynamics in the description of thoracic aorta haemodynamics and in the differentiation of progressive dilation in Marfan syndrome patients

Abnormal fluid dynamics at the ascending aorta may be at the origin of aortic aneurysms. This study was aimed at comparing the performance of computational fluid dynamics (CFD) and fluid–structure interaction (FSI) simulations against four-dimensional (4D) flow magnetic resonance imaging (MRI) data; and to assess the capacity of advanced fluid dynamics markers to stratify aneurysm progression risk. Eight Marfan syndrome (MFS) patients, four with stable and four with dilating aneurysms of the proximal aorta, and four healthy controls were studied. FSI and CFD simulations were performed with MRI-derived geometry, inlet velocity field and Young's modulus. Flow displacement, jet angle and maximum velocity evaluated from FSI and CFD simulations were compared to 4D flow MRI data. A dimensionless parameter, the shear stress ratio (SSR), was evaluated from FSI and CFD simulations and assessed as potential correlate of aneurysm progression. FSI simulations successfully matched MRI data regarding descending to ascending aorta flow rates (R2 = 0.92) and pulse wave velocity (R2 = 0.99). Compared to CFD, FSI simulations showed significantly lower percentage errors in ascending and descending aorta in flow displacement (−46% ascending, −41% descending), jet angle (−28% ascending, −50% descending) and maximum velocity (−37% ascending, −34% descending) with respect to 4D flow MRI. FSI- but not CFD-derived SSR differentiated between stable and dilating MFS patients. Fluid dynamic simulations of the thoracic aorta require fluid–solid interaction to properly reproduce complex haemodynamics. FSI- but not CFD-derived SSR could help stratifying MFS patients.This study was funded by Ministerio de Economía y Competitividad (grant no. RTC-2016-5152-1), Fundació la Marató de TV3 (grant no. 20151330), FP7 People: Marie-Curie Actions (grant no. 267128), Instituto de Salud Carlos III (grant nos PI14/0106 and PI17/00381) and ‘la Caixa’ Foundation. M.V. was funded by CompBioMed2, grant agreement ID: 823712, funded under: H2020-EU.1.4.1.3; and SILICOFCM, grant agreement ID: 777204, funded under: H2020-EU.3.1.5.Peer ReviewedPostprint (published version

Fluid–structure interaction analysis of eccentricity and leaflet rigidity on thrombosis biomarkers in bioprosthetic aortic valve replacements

This work intends to study the effect of aortic annulus eccentricity and leaflet rigidity on the performance, thrombogenic risk and calcification risk in bioprosthetic aortic valve replacements (BAVRs). To address these questions, a two-way immersed fluid–structure interaction (FSI) computational model was implemented in a high-performance computing (HPC) multi-physics simulation software, and validated against a well-known FSI benchmark. The aortic valve bioprosthesis model is qualitatively contrasted against experimental data, showing good agreement in closed and open states. Regarding the performance of BAVRs, the model predicts that increasing eccentricities yield lower geometric orifice areas (GOAs) and higher normalized transvalvular pressure gradients (TPGs) for healthy cardiac outputs during systole, agreeing with in vitro experiments. Regions with peak values of residence time are observed to grow with eccentricity in the sinus of Valsalva, indicating an elevated risk of thrombus formation for eccentric configurations. In addition, the computational model is used to analyze the effect of varying leaflet rigidity on both performance, thrombogenic and calcification risks with applications to tissue-engineered prostheses. For more rigid leaflets it predicts an increase in systolic and diastolic TPGs, and decrease in systolic GOA, which translates to decreased valve performance. The peak shear rate and residence time regions increase with leaflet rigidity, but their volume-averaged values were not significantly affected. Peak solid stresses are also analyzed, and observed to increase with rigidity, elevating risk of valve calcification and structural failure. To the authors' knowledge this is the first computational FSI model to study the effect of eccentricity or leaflet rigidity on thrombogenic biomarkers, providing a novel tool to aid device manufacturers and clinical practitioners.This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 713673. The research leading to these results has also received funding from “la Caixa” Foundation, with fellowship ID: LCF/BQ/DI18/11660044, and has been co-funded by the project CompBioMed2 (H2020-EU.1.4.1.3. Grant No. 823712)Peer ReviewedPostprint (published version

Optimization of condensed matter physics application with OpenMP tasking model

The Density Matrix Renormalization Group (DMRG++) is a condensed matter physics application used to study superconductivity properties of materials. It’s main computations consist of calculating hamiltonian matrix which requires sparse matrix-vector multiplications. This paper presents task-based parallelization and optimization strategies of the Hamiltonian algorithm. The algorithm is implemented as a mini-application in C++ and parallelized with OpenMP. The optimization leverages tasking features, such as dependencies or priorities included in the OpenMP standard 4.5. The code refactoring targets performance as much as programmability. The optimized version achieves a speedup of 8.0 × with 8 threads and 20.5 × with 40 threads on a Power9 computing node while reducing the memory consumption to 90 MB with respect to the original code, by adding less than ten OpenMP directives.This work is partially supported by the Spanish Government through Programa Severo Ochoa (SEV2015-0493), by the Spanish Ministry of Science and Technology (project TIN2015-65316-P), by the Generalitat de Catalunya (contract 2017-SGR-1414) and by the BSC-IBM Deep Learning Research Agreement, under JSA “Application porting, analysis and optimization for POWER and POWER AI”. This work was partially supported by the Scientific Discovery through Advanced Computing (SciDAC) program funded by U.S. Department of Energy, Office of Science, Advanced Scientific Computing Research and Basic Energy Sciences, Division of Materials Sciences and Engineering. This research used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725.Peer ReviewedPostprint (author's final draft

On the wake dynamics of an oscillating cylinder via proper orthogonal decomposition

The coherent structures and wake dynamics of a two-degree-of-freedom vibrating cylinder with a low mass ratio at Re=5300 are investigated by means of proper orthogonal decomposition (POD) of a numerical database generated using large-eddy simulations. Two different reduced velocities of U*=3.0 and U*=5.5, which correspond with the initial and super-upper branches, are considered. This is the first time that this kind of analysis is performed in this kind of system in order to understand the role of large coherent motions on the amplification of the forces. In both branches of response, almost 1000 non-correlated in-time velocity fields have been decomposed using the snapshot method. It is seen that a large number of modes is required to represent 95% of the turbulent kinetic energy of the flow, but the first two modes contain a large percentage of the energy as they represent the wake large-scale vortex tubes. The energy dispersion of the high-order modes is attributed to the cylinder movement in the inline and cross-stream directions. Substantially different POD modes have been found in the two branches. While the first six modes resemble those observed in the static cylinder or in the initial branch of a one-degree of freedom cylinder in the initial branch, the modes not only contain information about the wake vortexes in the super-upper branch but also about the formation of the 2T vortex pattern and the Taylor–Görtler structures. It is shown that the 2T vortex pattern is formed by the interplay between the Taylor–Görtler stream-wise vortical structures and the cylinder movement and is responsible for the increase in the lift force and larger elongation in the super-upper branch.This work has been partially financially supported by the Ministerio de Economía, Industria y Competitividad, Secretaría de Estado de Investigación, Desarrollo e Innovación, Spain (Ref. PID2020-116937RB-C21, PID2020-116937RB-C22), and by the European High-Performance Computing Joint Undertaking (JU) under grant agreement No 956104. The JU receives support from the European Union’s Horizon 2020 research and innovation programme and Spain, France, Germany. O. Lehmkuhl work is financed by a Ramón y Cajal postdoctoral contract by the Ministerio de Economía y Competitividad, Secretaría de Estado de Investigación, Desarrollo e Innovación, Spain (RYC2018-025949-I).Peer ReviewedPostprint (published version

Computational analysis of fluid dynamics at the asceding thoracic aorta in Marfan syndrome patients

Els aneurismes aòrtics són una dilatació progressiva i irreversible de la paret aòrtica, que pot causar la ruptura o dissecció dels vasos, el que resulta en una pèrdua catastròfica de sang que condueix a la mort. El tractament farmacològic inicial se centra en aturar el creixement per prevenir la ruptura, però requereix una reparació invasiva oberta o una reparació endovascular en pacients en risc. El maneig del pacient i l'estratificació de risc després del diagnòstic són crítics, especialment en l'aorta ascendent, ja que actualment no hi ha tractaments endovasculars disponibles. Segons les directrius actuals, el diàmetre aòrtic màxim és l'únic criteri geomètric o fluidodinàmic específic del pacient acceptat com a predictor clínic del risc de ruptura. No obstant això, l'anormal fluidodinàmica en l'aorta ascendent s'ha reportat àmpliament com una possible font d'aneurismes aòrtics i la seva comprensió podria millorar l'avaluació del risc del pacient. En aquest estudi, es va avaluar la fluidodinàmica en aortes de controls sans i pacients amb síndrome de Marfan. Per fer això, hem comparat el rendiment de les simulacions de dinàmica de fluids computacional i d'interacció fluid-estructura utilitzant imatges clíniques com a condicions específiques del pacient. També hem dissenyat un sistema in vitro que podria exposar les cèl·lules endotelials aòrtiques humanes a un entorn fluidodinàmic que imita el de les simulacions aòrtiques. L'estudi ha revelat, en pacients Marfan, que considerà l'elasticitat de la paret en les simulacions és essencial per obtenir amb precisió els valors fluidodinàmics que tenen el potencial d'estratificar aquests pacients. En aquest sentit, les simulacions d'interacció fluid-estructura han superat la fluidodinàmica computacional clàssica a un cost computacional moderat. Com a resultat d'aquest estudi, un paràmetre adimensional, la relació d'esforç tallant, ha determinat el seu potencial com a marcador de progressió d'aneurisma en pacients amb Marfan.Los aneurismas aórticos son una dilatación progresiva e irreversible de la pared aórtica, que puede causar la ruptura o disección de los vasos, lo que resulta en una pérdida catastrófica de sangre que conduce a la muerte. El tratamiento farmacológico inicial se centra en detener el crecimiento para prevenir la ruptura, pero se requiere una reparación invasiva abierta o una reparación endovascular en pacientes en riesgo. El manejo del paciente y la estratificación del riesgo después del diagnóstico son críticos, especialmente en la aorta ascendente, ya que actualmente no hay tratamientos endovasculares disponibles. Según las directrices actuales, el diámetro aórtico máximo es el único criterio geométrico o fluidodinámico específico del paciente aceptado como predictor clínico del riesgo de ruptura. Sin embargo, la anormal fluidodinámica en la aorta ascendente se ha reportado ampliamente como una posible fuente de aneurismas aórticos y su comprensión podría mejorar la evaluación del riesgo del paciente. En este estudio, se evaluó la fluidodinámica en aortas de controles sanos y pacientes con síndrome de Marfan. Para hacer esto, hemos comparado el rendimiento de las simulaciones de dinámica de fluidos computacional y de interacción fluido-estructura utilizando imágenes clínicas como condiciones específicas del paciente. También hemos diseñado un sistema in vitro que podría exponer las células endoteliales aórticas humanas a un entorno fluidodinámico que imita el de las simulaciones aórticas. El estudio ha revelado, en pacientes Marfan, que considerar la elasticidad de la pared en las simulaciones es esencial para obtener con precisión los valores dinámicos de los fluidos que tienen el potencial de estratificar a estos pacientes. En este sentido, las simulaciones de interacción fluido-estructura han superado la fluidodinámica computacional clásica a un costo computacional moderado. Como resultado de este estudio, un parámetro adimensional, la relación de esfuerzo cortante, ha demostrado su potencial como marcador de progresión de aneurisma en pacientes con Marfan.Aortic aneurysms are a progressive and irreversible dilation of the aortic wall, which can lead to vessel rupture or dissection, resulting in catastrophic blood loss leading to death. Initial pharmacological treatment is focused on growth arrest to prevent rupture, but invasive open repair or endovascular repair are required in patients at risk. Patient management and risk stratification after diagnosis are critical, especially in the ascending aorta since no endovascular treatments are currently available. According to current guidelines, maximum aortic diameter is the only patient-specific geometrical or fluidodynamic criterion accepted as clinical rupture risk predictor. However, abnormal fluid dynamics at the ascending aorta have been widely reported as potential origin of aortic aneurysms and their understanding could improve the risk assessment of patients. In this study, the fluid dynamics of aortae from healthy controls and patients with Marfan syndrome have been evaluated. To do so, we have compared the performance of computational fluid dynamics and fluid-structure interaction simulations using clinical imaging as patient-specific inputs. We have also designed an in vitro system that could expose human aortic endothelial cells to a fluidodynamic environment that mimics that of aortic simulations. The study has revealed, in Marfan patients, that considering the wall elasticity in simulations is critical to derive precisely fluid dynamic values that hold the potential to stratify such patients. In this sense, fluid-structure interaction simulations have outperformed classic computational fluid dynamics at a moderate computational cost. As a result of this study, a dimensionless parameter, the shear stress ratio, has shown its potential as marker of aneurysm progression in Marfan patients