Easy-to-implement hp-adaptivity for non-elliptic goal-oriented problems

Easy-to-implement hp-adaptivity for non-elliptic goal-oriented problems Felipe Vinicio Caro Guti´errez Supervised by David Pardo and Elisabete Alberdi December 2023 (cc)2023 FELIPE VINICIO CARO GUTIERREZ (cc by 4.0) Easy-to-implement hp-adaptivity for non-elliptic goal-oriented problems Felipe Vinicio Caro Guti´errez Supervised by David Pardo and Elisabete Alberdi December 2023 This dissertation has been possible with the support of the University of the Basque Country (UPV/EHU) grant No. PRE2018-084258; the BCAM “Severo Ochoa” accreditation of excellence (SEV-2017-0718); the Basque Government through the BERC 2018-2021 program; and the Consolidated Research Group MATHMODE (IT1294-19; IT1456-22) given by the Department of Education. i Acknowledgements I am grateful for the past four years, my life’s most remarkable and transformative period. Following two challenging years pursuing a master’s degree in physics applied to the ocean, fate guided me back to the enchanting realm of mathematics. First, I want to express my sincere gratitude to my supervisor, David Pardo, for graciously accepting me as his student. Joining David’s group has been incredible, and I am profoundly grateful for the opportunity. With his guidance and support, I am writing my PhD Dissertation. I aspire to emulate David’s many professional and personal qualities. He exemplifies the ideal supervisor, and I am privileged to have him as my mentor. I want to express my deep gratitude to Elisabete Alberdi, my co-supervisor, for her endless patience and persistent support throughout my journey. Our enlightening discussions in her office, where she introduced me to the intriguing realm of Finite Element Methods, will forever be embedded in my memory. I appreciate her tenacious confidence in my capabilities and unwavering trust in my work. I want to express my deepest gratitude to Vincent Darrigrand for his invaluable contributions throughout these years and during my stay in Toulouse. His unwavering support in assisting me with code-related challenges has been truly remarkable. Vincent’s guidance and endless advice on “doing things in the right way” have been instrumental in shaping my personal and professional growth. He has provided technical expertise and imparted invaluable knowledge on discerning the distinctions between the correct and incorrect approaches. Thanks to his guidance, I have developed my own set of principles and criteria in my personal and professional initiatives. I want to express my genuine gratitude to Julen Alvarez-Aramberri for his invaluable contributions to enriching my research activities. Julen has shared his profound professional skills throughout our three years as colleagues, impacting my academic journey indelibly. More importantly, he has imparted one of the most vital life lessons: the importance of being honest with oneself. Julen’s guidance has not only shaped me into a better scientist but also a better individual. I am genuinely grateful for his support and mentorship. I want to express my fervent appreciation to my wonderful colleagues at BCAM, including Jon Ander, Oscar, Carlos, and Ana. I am also grateful for the professional connections I have developed with them. In particular, I owe a debt of ii Acknowledgements gratitude to Judit, with whom I shared my first scientific congress. That experience was truly unforgettable, and I will always treasure the memories. I look forward to meeting any of them, whether at BCAM or elsewhere, as it will always be a joyful occasion. I want to thank Prof. Maciej Paszy´nski for his warm and invaluable support in welcoming me to the Institute of Computer Science at AGH University in Krakow. My time in Krakow was genuinely enriching, and the professional experiences I gained were immeasurable. I am deeply grateful for the guidance and mentorship provided by Prof. Paszy´nski during my stay. Additionally, I would like to express my sincere appreciation to Albert Oliver-Serra, Micha l Jungiewicz, Maciej Wo´zniak, Eirik Valseth, and Maciej Smo lka. Our shared professional experiences at AGH were precious, and I am grateful for the collaboration and insights we exchanged. Their contributions have undoubtedly played a significant role in shaping my academic and professional development. I sincerely tribute Dariel Hern´andez, who transcends the definition of a mere colleague at BCAM. Dariel’s tenacious dedication to his craft and infectious passion for dancing has enriched our professional relationship and fostered a genuine friendship. Words alone are insufficient to express the depth of my gratitude towards him. I am profoundly thankful for his constant support, camaraderie, and the joy he brings to our shared experiences. I am honored to call you both a colleague and a dear friend. Lastly, I would like to express my sincere appreciation to Cambr´e Dancing School and all those associated with it. You have given me unforgettable moments, travels, and shared dancing congresses that I will cherish forever. I want to give a special mention to Sergio and Andrea, whose guidance and support were instrumental in helping me discover my passion for dance. They provided me with the knowledge and encouragement to push beyond my limits and embrace the joy of dancing. Dancing has become an integral part of my life, thanks to the countless hours we spent honing our skills and creating unforgettable memories. I am forever grateful to Sergio, Andrea, and everyone at Cambr´e for enriching my life through the art of dance. iii Abstract The Finite Element Method (FEM) has become a foundational numerical technique in computational mechanics and civil engineering since its inception by Courant in 1943 [56]. Originating from the Ritz method and variational calculus, the FEM was primarily employed to derive solutions for vibrational systems. A distinctive strength of the FEM is its capability to represent mathematical models through the weak variational formulation of Partial Differential Equations (PDEs), facilitating computational feasibility even in intricate geometries. However, the search for accuracy often imposes a significant computational task. In the FEM, adaptive methods have emerged to balance the accuracy of solutions with computational costs. The h-adaptive FEM designs more efficient meshes by reducing the mesh size hlocally while keeping the polynomial order of approximation pfixed (usually p= 1,2). An alternative approach to the h-adaptive FEM is the p-adaptive FEM, which locally enriches the polynomial space pwhile keeping the mesh size hconstant. By dynamically adapting hand p, the hp-adaptive FEM achieves exponential convergence rates. Adaptivity is crucial for obtaining accurate solutions. However, the traditional focus on global norms, such as L2or H1, might only sometimes serve the requirements of specific applications. In engineering, controlling errors in specific domains related to a particular Quantity of Interest (QoI) is often more critical than focusing on the overall solution. That motivated the development of Goal-Oriented Adaptive (GOA) strategies. In this dissertation, we develop automatic Goal-Oriented (GO) hp-adaptive algorithms tailored for non-elliptic problems. These algorithms shine in terms of robustness and simplicity in their implementation, attributes that make them especially suitable for industrial applications. A key advantage of our methodologies is that they do not require computing reference solutions on globally refined grids. Nevertheless, our approach is limited to anisotropic pand isotropic h refinements. We conduct multiple tests to validate our algorithms. We probe the convergence behavior of our GO hand p-adaptive algorithms using Helmholtz and convection-diffusion equations in one-dimensional scenarios. We test our GO hpadaptive algorithms on Poisson, Helmholtz, and convection-diffusion equations in two dimensions. We use a Helmholtz-like scenario for three-dimensional cases to highlight the adaptability of our GO algorithms. iv Abstract We also create efficient ways to build large databases ideal for training Deep Neural Networks (DNNs) using hp Multi-Adaptive Goal-Oriented (MAGO) FEM. As a result, we efficiently generate large databases, possibly containing hundreds of thousands of synthetic datasets or measurements. v Resumen El m´etodo de elementos finitos (MEF) aproxima soluciones a ecuaciones diferenciales parciales (EDPs). Bas´andose en el m´etodo de Ritz y el c´alculo variacional, Courant desarroll´o el MEF en 1943 [56]. Desde entonces se ha convertido en una t´ecnica fundamental en la mec´anica computacional y la ingenier´ıa civil que se ha utilizado para resolver una amplia gama de problemas, incluyendo an´alisis estructural, mec´anica de fluidos y sistemas vibratorios. Una fortaleza distintiva del MEF es su capacidad para representar modelos matem´aticos a trav´es de la formulaci´on variacional d´ebil de las EDPs, facilitando la viabilidad computacional incluso en geometr´ıas intrincadas. Sin embargo, la b´usqueda de precisi´on a menudo impone una tarea computacional significativa. Debido a los altos costos computacionales de ciertos problemas, han surgido m´etodos adaptativos para equilibrar la precisi´on de las soluciones con los costos computacionales. El MEF adaptativo es un m´etodo num´erico que permite aproximar soluciones de forma m´as precisa con menor costo computacional. El MEF adaptativo hdise˜na mallas m´as eficientes reduciendo el tama˜no de malla localmente mientras mantiene el orden del polinomio de aproximaci´on pfijo (generalmente p= 1,2). Una alternativa al MEF adaptativo hes el MEF adaptativo p, que enriquece localmente el espacio de polinomios pmanteniendo constante el tama˜no de malla h. Al combinar din´amicamente ambos m´etodos, el MEF adaptativo hp logra tasas de convergencia exponenciales. El enfoque tradicional de la adaptatividad en normas globales (L2oH1) s´olo sirve para ciertas aplicaciones. En ingenier´ıa, controlar errores en dominios espec´ıficos relacionados con una cantidad de inter´es es a menudo m´as cr´ıtico que controlar errores globales. Debido a esta necesidad, surge la adaptatividad orientada a un objetivo espec´ıfico. En este trabajo, desarrollamos algoritmos autom´aticos orientados a un objetivo hp dise˜nados para problemas no el´ıpticos. Estos algoritmos se destacan en t´erminos de robustez y simplicidad en su implementaci´on, atributos que los hacen especialmente adecuados para aplicaciones industriales. Una ventaja clave de nuestras metodolog´ıas es que no requieren calcular soluciones de referencia en mallas globalmente refinadas. Sin embargo, nuestro enfoque se limita a refinamientos anisotr´opicos pe isotr´opicos h. Los resultados num´ericos 1D muestran la convergencia de nuestros algoritmos orientados a un objetivo, tanto hcomo p, usando las ecuaciones de Helmholtz y vi Resumen convecci´on-difusi´on. Adem´as, los resultados num´ericos en 2D muestran la convergencia de los algoritmos hp usando las ecuaciones de Poisson, Helmholtz y convecci´on-difusi´on. Tambi´en, probamos estos algoritmos hp en casos 3D con la ecuaci´on de Helmholtz para demostrar la versatilidad de nuestros algoritmos. Finalmente, extendemos nuestros algoritmos orientados a un objetivo hp para generar grandes bases de datos confiables e ideales para entrenar redes neuronales. Como resultado, mostramos la generaci´on eficiente de grandes bases de datos potencialmente con cientos de miles de datos sint´eticos. vii LIST OF FIGURES 7.9. Box plots for different adaptive grids with a threshold maximum relative error set at 10−5. ................... 79 7.10. Computational domain Ω, where homogeneous Dirichlet boundary conditions are imposed on ∂Ω. Additionally, we define Ωlas the support of the QoI l(ϕ), and Ωfas the support of the source function. 81 7.11. Absolute value of the solutions of our cross-shaped domain Poisson example................................. 82 7.12. hp-adapted meshes for our 1-sample cross-shaped domain example. 84 7.13. hp-adapted meshes for our 5-sample cross-shaped domain example. 85 7.14. hp-adapted meshes for our 10-sample cross-shaped domain example. 86 7.15. hp-adapted meshes for our 50-sample cross-shaped domain example. 87 7.16. hp-adapted meshes for our 100-sample cross-shaped domain example. 88 7.17. Box plots for different adaptive grids with a threshold maximum relative error set at 1.0%. ................... 89 7.18. Absolute value of the solutions of our Poisson example. . . . . . . 91 7.19. hp-adapted meshes for our 1-sample grid-based domain example. . 92 7.20. hp-adapted meshes for our 5-sample grid-based domain example. . 93 7.21. hp-adapted meshes for our 10-sample grid-based domain example. 94 7.22. hp-adapted meshes for our 50-sample grid-based domain example. 95 7.23. hp-adapted meshes for our 100-sample grid-based domain example. 96 7.24. Box plots for different adaptive grids with a threshold maximum relative error set at 1.0%.................... 97 xiv List of Tables 7.1. The computational cost based on the factorization cost of generating the database using the SAGO strategy. . . . . . . . . . . . . 80 7.2. The computational cost based on the factorization cost of generating the database using the MAGO strategy. . . . . . . . . . . . 80 7.3. The computational cost based on the factorization cost of generating the database using the SAGO strategy. . . . . . . . . . . . . 90 7.4. The computational cost based on the factorization cost of generating the database using the MAGO strategy. . . . . . . . . . . . 90 7.5. The computational cost based on the factorization cost of generating the database using the SAGO strategy. . . . . . . . . . . . . 98 7.6. The computational cost based on the factorization cost of generating the database using the MAGO strategy. . . . . . . . . . . . 98 xv 1. Introduction 1.1. Motivation In recent years, the Finite Element Method (FEM) has gained significant popularity as one of the most extensively utilized numerical techniques in computational mechanics and civil engineering. The beginnings of the FEM can be traced back to Courant’s pioneering work in 1943 [56], where he employed the Ritz method of numerical analysis and the minimization of variational calculus to derive approximate solutions for vibration systems. However, the computational success and widespread of the FEM can be attributed to the contributions of Turner et al. in 1956 [187] and Clough in 1960 [49]. The FEM has revolutionized various knowledge areas, driven by its primary application in structural mechanics [25, 27, 94, 165, 211]. Its significant impact extends to disciplines such as earthquake engineering, transforming the understanding and practices in these fields [48, 51]. Furthermore, through continuous research, the FEM application has expanded beyond structural mechanics. It has applications in various disciplines, including fluid mechanics, thermal analysis, and electrical engineering. The popularity of the FEM can be attributed to its capability to represent mathematical models through the weak variational formulation of Partial Differential Equations (PDEs). This formulation enables decomposing the problem domain into finite elements, with a corresponding number of unknowns called Degrees of Freedom (DoF). This decomposition makes it computationally feasible to obtain accurate solutions even in complex geometries (see, e.g., [69, 98, 109, 212, 213, 214] among others). We refer to the interested reader seeking a comprehensive mathematical foundation of the FEM to [33, 47, 127, 164]. Despite the significant advancements made in FEM over the past century [125], the computational cost of achieving highly accurate solutions remains a challenge. As the desired level of solution accuracy increases, the number of unknowns and computational resources required also escalate, potentially resulting in computationally expensive calculations that may be prohibitive in practice. The h-adaptive FEM addresses the computational costs of increasing solution accuracy. The method designs more efficient meshes by locally reducing the mesh size hwhile keeping the polynomial order of approximation pfixed (typically 1 1. Introduction p= 1,2). This dynamic adjustment of the mesh resolution through h-adaptivity facilitates the acquisition of accurate solutions while mitigating computational costs. The classical h-adaptive FEM involves locally refining elements by reducing their size h[18]. It has successfully achieved convergence rates regarding DoF through mesh adaptation [133]. Pioneering works by Babuˇska and Rheinboldt [15, 16, 17] have laid the foundation for this approach. However, it is essential to note that this method has limitations in overcoming algebraic convergence rates, resulting in slow convergence. Furthermore, the practical implementation of this method may be constrained by limited computer resources, as the computational demands can present significant challenges. An alternative approach to the h-adaptive FEM is the p-adaptive FEM [19, 42, 71, 183], which locally enriches the polynomial space pwhile keeping the mesh size hconstant. This method proves to be more practical for problems with smooth solutions, as it can achieve the same level of accuracy with a slightly refined mesh. One of the key advantages of the p-adaptive FEM is that by increasing the polynomial order of approximation p, it attains exponential convergence rates while simultaneously reducing the number of Degrees of Freedom (nDoF) required to achieve a desired level of accuracy. Non-smooth problems are prevalent in computational mechanics, especially in regions characterized by e.g. re-entrant corners and material interfaces, demanding precise simulations for accurate results. To address this, a combined approach of both adaptive techniques, namely the hp-adaptive FEM [87, 88], has emerged as an efficient alternative. This approach enables a more precise mesh refinement by adjusting the element size hnear singularities and the polynomial approximation order pin regions with smooth solutions. By dynamically adapting both hand p, the hp-adaptive FEM achieves exponential convergence rates, even in the presence of singularities [14], thereby offering higher accuracy for the same nDoF. To gain insight into the historical development of the FEM, it is valuable to refer to the works of Babuˇska [13] and Oden [126]. 1.2. Literature review 1.2.1. Advances in hp-adaptivity Adaptivity entails the selective modification of specific subdomains approximations within the computational domain rather than uniformly altering the approximation over the entire domain. By focusing on relevant subdomains, adaptivity aims to optimize the accuracy and efficiency of the solution while minimizing computational costs. This iterative process concentrates computational resources 2 1. Introduction on regions where accuracy improvements are most crucial, resulting in improved overall efficiency and accuracy of the solution. Adaptivity is critical in optimizing computational resources, particularly when they are limited. The primary objective is to achieve the highest level of accuracy while minimizing the nDoF required. The critical components for successful mesh adaptation include a posteriori error estimates [2, 3, 4] based on the computed solution, local error indicators, and a strategy that utilizes these indicators to adapt the mesh automatically [22]. Clough’s work [50] stands out as a pioneering contribution to developing a fully automated computer program for FEM analysis. Additionally, we shall mention Bank et al. [24] for their pioneering work in developing a global mesh adaptive algorithm. A wide range of h-adaptive algorithms are available, and here are a few notable examples. Deuflhard et al. [67] introduced the KASKADE code [75, 167], which utilizes hierarchical finite element bases as proposed by Yserentant [204]. In addition to KASKADE, other notable codes for addressing nonlinear problems include PLTMG, developed by Bank [23], and NFEARS [115, 116], developed by Mesztenyi and collaborators, among them. We also encounter, the work of Karniadakis et al. [101, 202, 206] in spectral/hp elements applied to incompressible and compressible flow problems. This approach combines the h-adaptive FEM with the desirable numerical properties of spectral methods. One of the complexities of this method is the requirement of two compatible meshes, which adds a challenge to the computational process. In addition to h-adaptive algorithms, B. A. Szab´o et al. [1, 73, 181, 182] employed a p-adaptive process and rely on a priori assumptions to design a mesh that is adequately adapted to the exact solution. Moreover, in hp-adaptive algorithms, G. W. Zumbusch [215] introduced an hp-adaptive algorithm based on the adaptive multilevel code, KASKADE. Additionally, J. Sch¨oberl [173] developed a mesh generator capable of generating new meshes (re-meshing) to support the hp-adaptive process. The work of Demkowicz et al. [62, 64, 66], and its applications [7, 8, 37, 80, 81, 84, 139, 140, 141, 145, 147, 149], proposed a method that produces optimal hp-meshes by minimizing the local projection error based on a reference solution. However, this approach requires implementing a Projection-Based Interpolation (PBI) and involves computationally expensive computations on a globally refined (h 2, p+1)-grid. In addition, ensuring continuity via the 1-irregularity rule leads to complex implementations. Other hp strategies in the field include the Texas three-step approach [128], which involves alternating between hand p-refinements. However, this method often produces suboptimal results. Another strategy, proposed in [5], is based on the local regularity of the exact solution. Its suitability for industrial applica3 1. Introduction tions remains uncertain, and it shares this limitation with specific Discontinuous Galerkin (DG) methods [10, 38, 39, 54, 63, 82, 83, 93, 154]. For a comprehensive review and comparison of existing hp-adaptive strategies up to 2014, please refer to [118]. Implementing high-order hp-meshes presents several challenges, particularly regarding the occurrence of hanging nodes during local h-refinements [68, 179]. These nodes must be constrained to ensure solution continuity. However, managing the data structures necessary to handle hanging nodes is complex and involves numerous technical difficulties. To simplify implementation, especially in higher dimensions, researchers [62, 184], among others, limit their algorithms to the 1-irregularity rule, which allows for a maximum of one level of hanging nodes. To address these challenges and reduce implementation complexity, Zander et al. introduced a novel data structure in their work [207, 208, 210] that supports hp-discretizations and inherently eliminates hanging nodes. Their approach utilizes hierarchical basis functions in hand pon a multi-level grid, employing uniform refinements with many Dirichlet nodes to ensure continuity and enable local refinements. Replacing global uniform refinements with isotropic refinements over selected elements eliminates hanging nodes while simplifying existing data structures for hp-refinements. Kopp et al. [104, 105] have extended these data structures to arbitrary dimensions [105] and space-time discretizations [104], expanding the approach’s applicability. In 2020, Darrigrand et al. [59] proposed a new automatic hp-adaptive meshrefinement strategy for elliptic problems that build upon Zander’s data structures [207, 208, 210]. Their approach not only eliminates mesh irregularities caused by hanging nodes but also avoids implementations of local projections (e.g., PBI [66]) that require the maintenance of multiple grids in the data structures. This easy-to-implement hp-strategy consists of a general (user-defined) refinement step followed by a specific mesh coarsening step. The method uses quadrilateral elements and alternates between global hor p-refinements with local and quasi-optimal hp-unrefinements (similarly to [29, 40]). In particular, the method eliminates basis functions with the lowest contributions to the solution energy at each hp-unrefinement step. The coarsening-based strategy described earlier provides a significant benefit. It can address and rectify inevitable “mistakes” that may have occurred due to undesired basis functions introduced during global refinements or in the preasymptotic regime. Moreover, subsequent unrefinement iterations can further enhance the results, improving upon any potential non-optimal results that may have arisen due to the approximate quasi-orthogonality assumption of the basis functions. Due to the inherent complexity of the hp-adaptive algorithms, both convergence 4 1. Introduction [59, 62, 161], and optimality are desirable properties. Optimality is obtaining the best solution using limited DoF. Canuto et al. [40] have provided proof of optimality in 1D and 2D problems, demonstrating the ability to achieve optimal solutions within the given DoF constraints. On the other hand, convergence measures how closely the computed solution approaches the exact solution of the problem. For a comprehensive understanding of convergence theory in the context of FEM, Ciarlet [46, 47] provides a valuable reference. We shall mention specific algorithms that have provided convergence proofs [30, 36, 41, 58]. 1.2.2. Advances in Goal-Oriented adaptivity Adaptivity aims to maximize the efficiency of computational resources while achieving the desired level of accuracy in the solution. The conventional approach to adaptivity, which estimates the error in a global norm (e.g., L2or H1), may sometimes fail to align with the specific requirements of applications. The need to control errors in specific Quantities of Interest (QoIs), rather than the overall energy of the solution, is common in many engineering applications. These requirements have driven the development of Goal-Oriented Adaptive (GOA) strategies. The development of Goal-Oriented (GO) adaptivity, aimed at efficiently approximating specific Quantity of Interest (QoI) with reduced computational cost, can be attributed to the pioneering works of Rannacher et al. [26, 162, 163]. Peraire and Patera [114, 137, 138, 152, 153, 172] further expanded upon these foundational studies. These researchers focused on deriving a posteriori error estimates that explicitly target the error in the QoIs. Traditional approaches for representing the error in the QoI involve utilizing the direct and adjoint solutions and the global bilinear form of the problem. This representation is then partitioned into local and computable quantities, which are used to guide local refinements (see, for example, [134]). In the context of goaloriented error estimation, Prudhomme and Oden [129, 130, 156, 157] developed a procedure that employs global functions defined over the entire computational domain to represent the error in the QoI. They also proposed a method to estimate lower and upper bounds on the QoI error using global energy error estimates, with the bounds determined by the sum of local indicators. The convergence analysis of adaptive algorithms can be attributed to the early works of D¨orfler and Morin [70, 121]. Before 2006, most goal-oriented methods were not proven to converge, although there were two exceptions [57, 120]. However, significant progress has been made since then, with the development of algorithms that exhibit exponential convergence rates for specific solution properties. For instance, Mommer [119] proposed an adaptive finite element method for approximating functionals of the solution of symmetric elliptic second-order 5 1. Introduction boundary value problems. In 2012, Pollock’s dissertation [155] presented a convergence theory for a class of goal-oriented adaptive finite element algorithms, including works on second-order non-symmetric [92] and semilinear [91] elliptic equations. Moreover, Feischl [76] performed an abstract analysis of optimal GO adaptivity. Numerical results demonstrating convergence have been provided by Darrigrand et al. [60, 61] and Valseth et al. [189], offering insightful examples. GO adaptivity has gained significant importance in various engineering applications, such as electromagnetics [142, 143, 144]. A noteworthy example is the work in [178], where the authors devised a GOA strategy that eschews explicit error estimates for guiding hp-refinements. Instead, they employ a suitable reference solution to recover an approximate error function, which provides a substantially more accurate approximation than the one obtained on the coarse mesh. In a related study [148], the authors further investigated the effectiveness of the GO hp-adaptive strategy by employing analytical techniques such as the Fourier transform and Bessel functions. Specifically, they focused on a problem involving the radiation of a loop antenna wrapped around a metallic cylinder into a conductive medium. The application of GO adaptivity in structural problems has its roots in the seminal works by Oden et al. [132] and Vemaganti et al. [197]. These pioneering studies laid the foundation for the theory and methodologies of GO adaptivity in modeling heterogeneous materials. Subsequent advancements in the field were made by Oden et al. [131], who explored GO adaptivity in discrete lattice models, and Romkes et al. [169], who investigated elastostatic problems of heterogeneous materials with material properties expressed as functions of random variables. In 2012, Jhurani et al. [96, 97] introduced a framework for numerical homogenization and GO adaptivity for non-linear lattice elasticity problems based on the MoorePenrose pseudo-inverse of element stiffness matrices. Furthermore, Panetier et al. [135], Verdugo et al. [198], and Waeytens et al. [200] made notable contributions to the application of GO adaptivity in the field of viscoelasticity. In the context of linear viscoelasticity, the works of Chamoin et al. [45] and Ladev`eze et al. [107, 108] are worth mentioning, as they developed error bounds for outputs of interest. The application of GO adaptivity in the context of fluid-structure interactions can be traced back to the pioneering research of Th. Dune [72]. Dune’s work introduced an innovative Eulerian framework for modeling fluid-structure interactions, which incorporated a posteriori GO error estimation as a fundamental component of the methodology. In [86], authors developed a nonlinear GO error estimation procedure tailored explicitly to analyze Navier-Stokes incompressible fluid flows with structural interactions. During his Ph.D. dissertation, K. G. van der Zee made significant contributions to fluid-structure interactions [190], fur6 1. Introduction ther extending the understanding and application of GO adaptivity in this area. One example of his contributions can be seen in [193], where authors developed a GO error estimator tailored for finite-element discretizations of fluid-structureinteraction problems. Their study focused on a model problem involving steady Stokes flow in a 2D channel with a flexible section of the channel wall. Moreover, valuable contributions were made in free-boundary problems, as demonstrated in [194, 195]. Additionally, in [196], GO error estimation in the context of free-boundary problems, where GO error estimation was applied using isogeometric analysis, was explored. K. G. van der Zee and colleagues also made other noteworthy contributions. In [192], they presented rigorous derivations of exact linearized adjoints for a coupled fluid-structure problem. At the same time, in [191], they developed a posteriori estimate of errors in the QoI for the nonlinear system of evolution equations embodied in the Cahn-Hilliard model of binary phase transition. 1.3. Main contributions of the dissertation The present dissertation summarizes the main contributions as follows. First, we extend the energy-based approach proposed by Darrigrand et al. [59] to the context of hand p-GOA algorithms. To achieve this, we combine the energy-based approach with an alternative pseudo-dual operator for representing the error in the QoI [60]. Our proposed approach is based on defining a new representation for the residual error of the adjoint problem, which exhibits better properties than the original bilinear form (e.g., positive definiteness). This new representation has been successfully used in previous studies [61, 123] and allows us to compute the error in the QoI in a way similar to classical approaches. As a result, we obtain automatic GO hand p-adaptive algorithms for non-elliptic problems. Second, we extend the energy-based-adaptive hp-strategy proposed by Darrigrand [59] to non-elliptic equations. To achieve this, we provide an alternative estimation of the energy contribution in terms of an inner product that depends on the bilinear form of the problem. As a result, we obtain an automatic hpadaptive algorithm for non-elliptic problems. Third, we extend Darrigrand’s strategy [59] to GOA approaches for both elliptic and non-elliptic problems. To achieve this, we use the adjoint problem to construct an upper bound of the error representation expressed in terms of an inner product that depends on the bilinear form of the problem. As a result, we obtain an automatic GO hp-adaptive algorithm for elliptic and non-elliptic problems. Our algorithms exhibit robustness and straightforward implementation, making them suitable for industrial applications. Notably, our approaches do not 7 1. Introduction require the computation of reference solutions on very fine grids, unlike other methods such as [66]. Our approach is limited to anisotropic pand isotropic hrefinements. However, recent work by Zander et al. [209] has extended multi-level data structures to support anisotropic h-refinements. To showcase the effectiveness of our algorithms, we demonstrate the convergence of our hand p-adaptive algorithms in 1D Helmholtz and convection-diffusion equations. Additionally, we test and analyze our hp-adaptive algorithm in three different 2D problems based on Poisson, Helmholtz, and convection-diffusion equations. Furthermore, we provide numerical results for a 3D Helmholtz-like problem. Although it is possible to construct suitable a posteriori error estimators [4, 26, 163] to enhance the refinement step of the algorithm, this possibility is outside the scope of this dissertation. Lastly, we extend Caro et al.’s [43] work to parametric PDEs. We develop an efficient way to generate reliable databases containing hundreds of thousands of synthetic data or measurements while minimizing computational costs for training Deep Neural Networks (DNNs). Due to the limited capabilities of Deep Learning (DL) techniques in solving PDEs, we approximate the forward operator. We adopt a modified version of the GO hp-adaptive FEM strategy [43, 44], unlike Hashemian et al.’s [90] study, which used a refined Isogeometric Analysis (IGA) approach to create databases of up to 100,000 Earth models. 1.4. Outline In this dissertation, we discuss the data structures presented by Zander et al. [207, 208, 210] in Section 2.1 of Chapter 2. We also introduce the concept of removable basis functions in Section 2.1.1, an essential idea in this dissertation. In Chapter 3, we present the adaptive strategy and element-wise error indicators. Our coarsening policy is introduced in Section 3.1, and we define the concept of projectors in Section 3.2, which applies to a single finite element mesh. We derive error indicators in Section 3.3, which guide the adaptivity for energy-norm and GO adaptivity. The methodology is applied to both elliptic and non-elliptic problems. Chapter 4 provides numerical results for 1D problems using the hand p-GOA algorithms proposed in this dissertation. We detail the proposed algorithms in Section 4.1 and outline the error indicators used in our hand padaptive algorithms in Section 4.2. We present numerical results demonstrating the convergence of the proposed hand p-GOA algorithms for 1D Helmholtz and convection-diffusion equations in Section 4.3. Finally, Section 4.4 summarizes the numerical results presented in this chapter. Chapter 5 illustrates the performance of our hp-adaptive algorithm numerically. We demonstrate the exponential convergence behavior of the approach for various 2D problems. Specifically, Section 8 3. Goal-Oriented coarsening strategy Algorithm 2: hp-unrefinement policy Input: A given mesh Output: An hp-unrefined mesh do Compute the solution on the current mesh; Compute the element-wise error indicators; Unrefine the mesh by eliminating the removable basis functions with low error indicators; When no contributions are below a given tolerance, exit; end ; 3.2. Projectors For dimension d∈ {1,2,3}, let Ω ⊂Rdbe an open bounded domain with a Lipschitz-continuous boundary ∂Ω, and let H(Ω) be a Hilbert functional space on Ω (simply denoted as Hin the following). For a given continuous bilinear form bdefined on H×H, let us define our problem with the following abstract variational formulation: Find u∈Hsuch that b(u, ϕ) = f(ϕ),∀ϕ∈H,(3.1) where fis a linear form. The discrete counterpart of this abstract variational formulation reads as follows: Find uF∈HFsuch that b(uF, ϕF) = f(ϕF),∀ϕF∈HF,(3.2) where HF:= span {ϕ1, . . . , ϕnF}is a finite element discretization Tof H, such that HF⊂H,F={ϕi}nF i=1 is a set of basis functions ϕi, and nF= dim (HF). Besides, uFcorresponds to the Galerkin approximation of uin HF. Some hp techniques handle a fine and a coarse mesh at the same time (see, e.g., [62, 64]). In addition to the coding difficulties derived from this fact, they typically need to define and implement projection operators (such as the Projection-Based Interpolation (PBI)) to link both grids. One of the main characteristics of our “painless” approach is continuously operating on a single mesh. While it simplifies the implementation, it requires defining a simple projector that simulates the presence of a coarse mesh without the trouble of handling one. 15 3. Goal-Oriented coarsening strategy For a given subset of basis functions S ⊂ F that generates the space HS⊂HF, we define our projection operator ΠS F:HF−→ HSas ΠS FuF:=X ϕi∈S uiϕi,(3.3) that is, we extract the coefficients of uFcorresponding to the basis functions in S, and we set the others to zero. For any element K, we denote by RKthe set of removable basis functions (see Section 2.1.1) associated to K, by |RK|its cardinality, and by HRKits associated space. Additionally, we define the subset of essential basis functions EKas EK:=F \ RK, while its associated space is denoted by HEK. These spaces satisfy that HEK⊂HF,HRK⊂HF, and HF=HEK∪HRK, with HEK∩HRK=∅. As a consequence, we can express any uF∈HF, as: uF= ΠEK FuF+ ΠRK FuF.(3.4) Since we consider a single mesh at a time, the solution uEKin EKassociated to eq. (3.2) is, in fact, never computed. Instead, we employ the projection of uF into EKto approximate it when necessary. 3.3. Error indicators Let ∥·∥ebe the energy norm associated with the Hilbert space H. For elliptic problems (given by symmetric and positive-definite bilinear forms), we define this energy from the bilinear form of the problem b, that is, ∥·∥2 e=b(·,·). For each non-elliptic problem, we shall define an alternative operator anot necessarily the original bilinear form such that |b(ϕ, ψ)|≤|a(ϕ, ψ)| ∀ϕ, ψ ∈Hand ∥·∥2 e=a(·,·) is the energy norm of the problem (i.e., adefines an inner product). We emphasize that the choice of these operators might highly influence the results of the adaptive process, which is usually an essential ingredient of adaptive strategies. With this in mind, our objective is to provide representative element-wise error indicators that drive the hp-coarsening steps (see Algorithm 2). For that, we consider isotropic and anisotropic indicators that are problem-dependent. In the following subsections, we derive only the isotropic error estimators ηK,∀K∈ T for a wide range of problems (see [59], for anisotropic indicators). To select what basis functions to unrefine, we compute the error indicators’ average (per degree of freedom) for the removable basis functions. We subsequently eliminate the removable basis function whose contribution is smaller than a percentage of this average. For further details and implementation technicalities, see [59]. 16 3. Goal-Oriented coarsening strategy In the following, we summarize the results from Darrigrand et al. [59] for elliptic energy-norm-based adaptive problems from the energy-norm perspective. After that, we extend these results to non-elliptic equations, and finally, we consider GO adaptivity applied to elliptic and non-elliptic problems. We can obtain all the proposed results by assuming (quasi)-b-orthogonality of the basis functions. However, this assumption is strong and unneeded for the energy-based adaptivity, and, therefore, we only employ it for GO adaptivity. To do so, let us denote by “≲” the inequality that holds up to a constant; that is, we represent a≤Cb by a≲b, with a, b, C ∈R, and let us define the L2-inner product of two possible complex and possibly vector-valued functions g1and g2 as: ⟨g1, g2⟩L2(Ω) =ZΩ (g∗ 1)Tg2dΩ,(3.5) where gTis the transpose of g, while g∗ 1represents the complex conjugate of g1. 3.3.1. Energy-norm based elliptic problems For a given element K∈ T , the objective is to quantify how much energy we lose in the solution when removing a subset of basis functions of the set of removable basis functions RK. Specifically, we want to compute ∥uF−uEK∥2 e. If this quantity is small, we guarantee that the energy of the removed set of basis functions is insignificant. Therefore, the fine and the unrefined meshes would provide comparable results. Analogously to Cea’s lemma proof, we derive: ∥uF−uEK∥2 e=b(uF−uEK, uF−uEK) (3.6) =buF−uEK, uF−ΠEK FuF+buF−uEK,ΠEK FuF−uEK(3.7) ≤ ∥uF−uEK∥e uF−ΠEK FuF e,(3.8) where we have used the b-orthogonality of uF−uEKwith HEKand the CauchySchwarz inequality. Therefore, ∥uF−uEK∥e≤ uF−ΠEK FuF e= ΠRK FuF e.(3.9) It is then natural to define the following element-wise error indicator: ηK:= ΠRK FuF  2 e,∀K∈ T .(3.10) 3.3.2. Extension to energy-based non-elliptic problems Again, our purpose is to compute ∥uF−uEK∥2 eto eliminate the removable basis functions with a low contribution to the solution. For that, let us start with the 17 3. Goal-Oriented coarsening strategy triangular inequality, which provides that ∥uF−uEK∥e≤ uF−ΠEK FuF e+ ΠEK FuF−uEK e.(3.11) Let us assume now that bsatisfies the discrete inf-sup condition: ∃γ > 0,inf ϕ∈HEK sup ψ∈HEK b(ϕ, ψ) ∥ϕ∥e∥ψ∥e ≥γ. (3.12) Then, using this inequality and the b-orthogonality of uF−uEKwith respect to HEK, we control the second term of eq. (3.11): γ ΠEK FuF−uEK e≤sup ψ∈HEK bΠEK FuF−uEK, ψ ∥ψ∥e (3.13) ≤sup ψ∈HEK bΠEK FuF−uF, ψ+b(uF−uEK, ψ) ∥ψ∥e (3.14) ≤sup ψ∈HEK Mb ΠEK FuF−uF e∥ψ∥e ∥ψ∥e (3.15) ≤Mb uF−ΠEK FuF e,(3.16) where Mbis the continuity constant of b. Therefore, ∥uF−uEK∥2 e≲ uF−ΠEK FuF  2 e= ΠRK FuF  2 e.(3.17) Accordingly, we define the element-wise indicator as follows: ηK:= ΠRK FuF  2 e,∀K∈ T .(3.18) The coarsening step will unrefine the elements that exhibit small ηK. Therefore, eq. (3.17) ensures that the problem’s energy loss will be negligible when removing these basis functions. 3.3.3. Extension to Goal-Oriented adaptivity GOA techniques aim to approximate specific quantities of finite element solutions rather than the global energy of the problem. These quantities with particular engineering applications are often called influence functions or Quantities of Interest (QoIs). Thus, the objective is to produce a space HFwith a minimum dimension such that the error in the Quantity of Interest (QoI) is below a userprescribed tolerance. To control the error in the QoI, we introduce the following adjoint problem [130, 156] associated to eq. (3.1): 18 3. Goal-Oriented coarsening strategy Find v∈Hsuch that b(ϕ, v) = l(ϕ),∀ϕ∈H,(3.19) where l:H−→ Ris a linear continuous form. Hence, the QoI of the solution uF is denoted by l(uF). The discrete equivalent of this problem is given by: Find vF∈HFsuch that b(ϕF, vF) = l(ϕF),∀ϕF∈HF,(3.20) where vFstands for the Galerkin approximation of the solution vto the adjoint problem associated with the space HF. For the mathematical analysis, we also consider the solution vEKin EKassociated with eq. (3.20), although we never compute it in practice. For a given element K∈ T , we want to quantify how much the QoI changes when removing some basis functions from the set of removable basis functions RKassociated with K. That is, we need to control |l(uF)−l(uEK)|,∀K∈ T . Since HEK⊂HF, Galerkin orthogonality ensures that b(uF−uEK, ϕ) = 0,∀ϕ∈HEK.(3.21) Then, l(uF)−l(uEK) = b(uF−uEK, vF) = b(uF−uEK, vF−vEK).(3.22) Using eq. (3.4) on vF, we have that: l(uF)−l(uEK) = buF−uEK,ΠRK FvF+ ΠEK FvF−vEK(3.23) =buF−uEK,ΠRK FvF+buF−uEK,ΠEK FvF−vEK.(3.24) Again, thanks to Galerkin orthogonality the second term vanishes. Then, applying eq. (3.4) on uFto the remaining term, we have that l(uF)−l(uEK) = bΠRK FuF+ ΠEK FuF−uEK,ΠRK FvF(3.25) =bΠRK FuF,ΠRK FvF+bΠEK FuF−uEK,ΠRK FvF.(3.26) Additionally, if we assume that EKis (quasi) b-orthogonal to RKdue to the (quasi)-orthogonality assumption of the basis functions, then bΠEK FuF−uEK,ΠRK FvF≃0,(3.27) 19 3. Goal-Oriented coarsening strategy and consequently, |l(uF)−l(uEK)| ≃ bΠRK FuF,ΠRK FvF≤aΠRK FuF,ΠRK FvF.(3.28) Then, we define the element-wise indicators as ηK:=aΠRK FuF,ΠRK FvF,∀K∈ T .(3.29) Here again, eq. (3.28) ensures that eliminating the basis functions associated with small indicators during the coarsening process should have a limited effect on the error of the QoI. Remark: Since bis continuous on Hwith respect to the energy norm, we also have |l(uF)−l(uEK)| ≃ bΠRK FuF,ΠRK FvF≲ ΠRK FuF e ΠRK FvF e,(3.30) and we could also define the element-wise indicators based on the above equation. Notice that if we select lto be the source term in the adjoint problem defined by eq. (3.19), with eq. (3.30) we recover the element-wise indicators derived previously in eqs. (3.10) and (3.18). However, in the forthcoming numerical results, we employ the estimators based on eq. (3.29). 3.3.4. Error indicators using a pseudo-dual operator The adjoint problem is often employed in the literature to guide GO refinements (see, e.g., [130, 156]). In addition, for the case of indefinite or non-symmetric problems, we further need to introduce an inner product (symmetric and positive definite form) to guide the refinements. To overcome this issue, we first define ΠEK FvFas a projection of the dual solution vFinto a given subset of essential basis functions EK. Such projections can be trivially implemented in the context of the multi-level data structures proposed in Zander et al. [207, 208, 210]; but not when using traditional data structures like those described in [62, 64, 65]. Then, we introduce a pseudo-dual bilinear form ˆ b, in this case, defined by the 1D Laplace operator (although it is possible to select other symmetric positive definite bilinear forms) to solve the following residual-based pseudo-dual problem: Find ˜εsuch that ˆ b(ϕF,˜ε) = l(ϕF)−bϕF,ΠEK FvF,∀ϕ∈H.(3.31) 20 3. Goal-Oriented coarsening strategy In previous work, Romkes et al. [168] introduced an elliptic error representation. Later, Darrigrand et al. [60] utilized this concept in traditional data structures. However, their approach required dealing with two grids (fine and coarse) and PBI operators [62, 64, 66], which made implementation and mathematical analysis highly complex. In contrast, we define problem (3.31) using a simpler approach. We use the projection of vFinto EK, denoted as ΠEK FvF. Thus, we define ηKas the error indicator associated with the element Kas follows ηK:= ˆ bΠRK FuF,˜ε,∀K∈ T ,(3.32) i.e., we define the operator a(·,·) simply as a(·,·) = ˆ b(·,·). 21 4. 1D Numerical results for Goal-Oriented hand p-adaptivity This chapter describes our hand p-adaptive strategies tailored to address 1D problems governed by Helmholtz and convection-diffusion equations. These adaptive algorithms offer a distinctive approach, focusing on minimizing the error in a specific Quantity of Interest (QoI) rather than the global error. We will comprehensively describe our adaptive algorithms, elaborating on the error indicators utilized throughout this chapter. Our approach incorporates a pseudo-dual operator given by eq. (3.31), which proves advantageous for non-elliptic GoalOriented (GO) problems. The numerical results were published in Caro et al. [44]. The hand p-adaptive algorithms proposed in this chapter follow the next refinement pattern: first, we perform a global and uniform hor p-refinement (for the hand p-adaptive versions, respectively). Then, we perform a coarsening step, removing some basis functions. This procedure is illustrated in Algorithm 1, and it was already introduced in [59] in the context of energy-norm adaptivity. The critical part is the coarsening step we depict in Algorithm 2. The critical step here is the computation of the element-wise error indicators described in section 3.3. In particular, we employ the eq. (3.32) to compute the error indicators utilized throughout this chapter. To illustrate the performance of our adaptive strategies, we consider two problems governed by Helmholtz and convection-diffusion equations. We provide the evolution of the relative error in the QoI for hand p-adaptivity and different values of the Partial Differential Equation (PDE) parameters. To define the relative error in the QoI, we compute l(u) on a globally refined mesh. Then, we define the relative error in a QoI in percentage as follows: eQoI rel :=|l(u)−l(uTc)| |l(u)|·100,(4.1) where uis the solution in a fine grid, while uTcis the solution associated with a coarser unrefined mesh. In some cases where the exact solution is available, we 22 4. 1D Numerical results for Goal-Oriented hand p-adaptivity will replace the fine grid solution uwith the exact solution, and we will directly compute eQoI rel . 4.1. Helmholtz Goal-Oriented problem Let us consider the following wave propagation problem: Find usuch that, −u′′ −k2u=1(0,2 5)in (0,1) ,(4.2) u(0) = 0,(4.3) u′(1) = 0.(4.4) We define the QoI as l(u) = 5·R4 5 3 5 u dx. Figures 4.1 and 4.2 show the evolution of eQoI rel by using hand p-adaptivity, respectively. Note that the larger the number of Degrees of Freedom (nDoF) per wavelength, the faster eQoI rel decreases. For example, in Figure 4.1, for k= 7·2π, 10 Degrees of Freedom (DoF) per wavelength are sufficient to enter into the so-called asymptotic regime. In contrast, for k= 28·2π, we need to consider at least 40 DoF per wavelength. In Figure 4.2, we select the initial mesh size such that the nDoF per wavelength is at least 3. This way, we satisfy the Nyquist rate. Both Figures 4.1 and 4.2 show optimal convergence rates in both hand p-adaptivity. As a curiosity, we observe that the curves in Figure 4.1 are parallel, while the ones in Figure 4.2 coincide. That occurs due to dispersion (pollution) error, which quickly disappears with the p-method. Figure 4.3 shows the solutions for the case k= 7 ·2π. We also provide the corresponding hand p-adaptive meshes. For the p-adaptive mesh, we show the mesh obtained in the 6th iteration, containing high approximation orders. To visualize the h-adaptive mesh, we display the mesh obtained in the 5th iteration. During this iteration, the refinements are denser in areas where the solution changes rapidly or exhibits sharp gradients. As a result, the element sizes in these regions are smaller than in other areas. 23 4. 1D Numerical results for Goal-Oriented hand p-adaptivity 3 10 40 200 600 10−2 10−1 100 101 102 103 DoF per wavelength Relative error in the QoI (%) k= 7 ·2π k = 14 ·2π k= 28 ·2π Figure 4.1.: Evolution of eQoI rel using h-adaptivity. Initial mesh size h=1 30 and uniform p= 1. 4 10 20 30 10−11 10−8 10−5 10−2 101 DoF per wavelength Relative error in the QoI (%) k= 7 ·2π k = 14 ·2π k= 28 ·2π Figure 4.2.: Evolution of eQoI rel using p-adaptivity. Uniform mesh size h=1 30. 24 5. 2D Numerical results for hp-adaptivity Ωf Ωl Ω ΓD 0 1 1 Figure 5.1.: Our singular Poisson example is defined over the domain Ω. The Dirichlet boundary is denoted by ΓD. The source function is supported on Ωf, and the QoI l(ϕ) is supported on Ωl. 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x y 0 0.71 1.41 ·10−2 Direct solution (a) Solution to the direct problem. 0 0.2 0.4 0.6 0.8 10 0.2 0.4 0.6 0.8 1 x y 0 0.11 0.23 Adjoint solution (b) Solution to the adjoint problem. Figure 5.2.: Direct and adjoint solutions of our singular Poisson example. 31 5. 2D Numerical results for hp-adaptivity 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x y 1234567891011 Order of approximation (a) Final hp-adapted mesh with polynomial orders in the x-direction. 0 0.2 0.4 0.6 0.8 10 0.2 0.4 0.6 0.8 1 x y 1234567891011 Order of approximation (b) Final hp-adapted mesh with polynomial orders in the y-direction. 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x y (c) Final h-adapted mesh, p= 1. 101102103104105106 10−10 10−8 10−6 10−4 10−2 100 Number of DoFs, N(log scale) Relative error in % (log scale) hp (p+ 2) h(p= 1) h(p= 2) (d) Evolution of eQoI rel in the process. Figure 5.3.: Final hand hp-adapted meshes for our singular Poisson example. 32 5. 2D Numerical results for hp-adaptivity 5.2. Wave propagation problem We consider the following non-elliptic problem based on Helmholtz’s equation. Find usuch that, −∆u−k2u=1Ωfin Ω,(5.7) u= 0 on ΓD,(5.8) ∇u·n = 0 on ΓN,(5.9) where Ω = (0,1)2\1 4,3 42⊂R2, Ωf=0,1 42⊂Ω, and k= (8 ·2π, 2π). The complex-valued kindicates the medium is lossy. ΓDand ΓNstand for the parts of the boundary ∂Ω where we impose homogeneous Dirichlet and Neumann boundary conditions, respectively. From eq. (5.3), we define Ωl=3 4,12⊂Ω. Figure 5.4 shows the domain of this hyperbolic (non-elliptic) problem. Ωf Ωl ΓD ΓN Ω Figure 5.4.: Our wave propagation example is defined over the domain Ω with a hole in the middle (marked in gray). The Dirichlet boundary is denoted by ΓD, while the Neumann boundary condition is denoted by ΓN. The source function is supported on Ωf, and the QoI l(ϕ) is supported on Ωl. 5.2.1. Energy-norm adaptivity For GO adaptivity, Figures 5.5a and 5.5b show the solutions to the direct and adjoint problems, respectively. Figure 5.6 shows the final hand hp-adapted meshes and Figure 5.7 shows the evolution of ˜eenergy rel and eQoI rel . The initial uniform 33 5. 2D Numerical results for hp-adaptivity mesh is composed of twelve root elements. We perform a double h-hierarchical refinement on the initial mesh to obtain a fine mesh to start the adaptivity. 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x y 0 4.99 9.98 ·10−4 Direct solution (a) Solution to the direct problem. 0 0.2 0.4 0.6 0.8 10 0.2 0.4 0.6 0.8 1 x y 0 0.54 1.08 ·10−2 Adjoint solution (b) Solution to the adjoint problem. Figure 5.5.: Absolute value of the direct and adjoint solutions of our wave propagation example in a lossy medium. For the h-adapted case, we observe heavy refinements around the source; however, almost no refinement occurs near the QoI. That happens due to the lossy nature of the problem. As a result, we observe a proper energy-norm convergence, as shown in Figure 5.7a, but a poor convergence behavior in the QoI, as demonstrated in Figure 5.7b. When implementing the hp-adaptive strategy, the refinements tend to be denser around the source than in the vicinity of the QoI. However, some non-trivial refinements still occur around the QoI. Despite this, the relative error in the QoI, denoted as eQoI rel , still converges to a level of 10−3% with just 20k unknowns. We define the operators b(·,·) and a(·,·) associated with the above problem as follows: b(·,·):=⟨∇· ,∇·⟩L2(Ω) −k2⟨· ,·⟩L2(Ω) , a (·,·):=⟨∇· ,∇·⟩L2(Ω)+k2⟨· ,·⟩L2(Ω). (5.10) Once more, ∥·∥2 e=a(·,·) defines our energy norm and |b(ϕ, ψ)| ≤ |a(ϕ, ψ)|,∀ϕ, ψ ∈ H. 34 5. 2D Numerical results for hp-adaptivity 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x y 1234567891011 Order of approximation (a) Final hp-adapted mesh with polynomial orders in the x-direction. 0 0.2 0.4 0.6 0.8 10 0.2 0.4 0.6 0.8 1 x y 1234567891011 Order of approximation (b) Final hp-adapted mesh with polynomial orders in the y-direction. 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x y (c) Final h-adapted mesh, p= 1. Figure 5.6.: Final hand hp-adapted meshes for our wave propagation example in a lossy medium. 35 5. 2D Numerical results for hp-adaptivity 102103104105 10−12 10−9 10−6 10−3 100 Number of DoFs, N(log scale) Relative error in % (log scale) hp (p+ 2) h(p= 1) h(p= 2) (a) Evolution of ˜eenergy rel in the process. 102103104105 10−3 10−2 10−1 100 101 102 103 Number of DoFs, N(log scale) hp (p+ 2) h(p= 1) h(p= 2) (b) Evolution of eQoI rel in the process. Figure 5.7.: Energy-norm adaptivity. Evolution of ˜eenergy rel and eQoI rel in our wave propagation example in a lossy medium. 103104 10−7 10−5 10−3 10−1 101 103 Number of DoFs, N Relative error in % (log scale) GOA energy-norm (a) Evolution of goal-oriented adaptivity. 103104 10−12 10−9 10−6 10−3 100 Number of DoFs, N GOA energy-norm (b) Evolution of energy-norm adaptivity. Figure 5.8.: Convergence history of eQoI rel and ˜eenergy rel for the energy-norm and GO hp-adaptive strategies. 36 5. 2D Numerical results for hp-adaptivity 5.2.2. Goal-Oriented adaptivity Figure 5.9 shows the final hand hp-adapted meshes and the evolution of eQoI rel . The initial mesh is uniform and composed of twelve root elements. As in the energy-norm adaptivity, we perform a double h-hierarchical refinement on the initial mesh to obtain a fine mesh to start the adaptivity. We observe heavy h-refinements around four localized singularities at the interior corners of the domain. In addition, we recover exponential convergence rates for the hand for the hp-adaptive versions. As a result, we construct a hp-adapted mesh with 20k unknowns that delivers a relative error in the QoI of 10−6% (three orders of magnitude better than in Figure 5.7b). To better illustrate this idea, Figure 5.8 compares the evolution of eQoI rel and ˜eenergy rel when executing the energy-norm and the GO hp-adaptive strategies in our wave propagation example in a lossy medium. Figure 6.4a shows a relative error in the QoI three orders of magnitude better when performing GO adaptivity than considering energy-norm adaptivity. Figure 6.4b shows that the ˜eenergy rel rapidly converges when employing energy-norm adaptivity, while with the hp-adaptive GO strategy, the rapid initial convergence stagnates at the level of 10−6%. As expected, this situation is also noticeable in terms of h-adaptivity (see Figures 5.7 and 5.9d). 37 5. 2D Numerical results for hp-adaptivity 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x y 1234567891011 Order of approximation (a) Final hp-adapted mesh with polynomial orders in the x-direction. 0 0.2 0.4 0.6 0.8 10 0.2 0.4 0.6 0.8 1 x y 1234567891011 Order of approximation (b) Final hp-adapted mesh with polynomial orders in the y-direction. 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x y (c) Final h-adapted mesh, p= 1. 102103104105 10−7 10−5 10−3 10−1 101 103 Number of DoFs, N(log scale) Relative error in % (log scale) hp (p+ 2) h(p= 1) h(p= 2) (d) Evolution of eQoI rel in the process. Figure 5.9.: Final hand hp-adapted meshes for our singular GO wave propagation example in a lossy medium and the evolution of eQoI rel . 38 5. 2D Numerical results for hp-adaptivity 5.3. Convection-dominated diffusion problem 5.3.1. Convection-dominated diffusion: example 1 We consider the following non-elliptic problem based on the convection-dominated diffusion equation. Find usuch that, −ε∆u+σ· ∇u=fin Ω,(5.11) u= 0 on ∂Ω. The selection of a suitable norm to measure the error in problems based on eq. (5.11) is an open research subject. For instance, authors of [77, 78] use the standard energy norm, in [79] a balanced norm, and in [180, 199] different norms from the previous ones. In here, we define the operators b(·,·) and a(·,·) associated with the above problem as follows: b(·,·):=ε⟨∇· ,∇·⟩L2(Ω) +⟨σ∇· ,·⟩L2(Ω) , a (·,·):= (ε+C)⟨∇· ,∇·⟩L2(Ω) , (5.12) where ∥·∥2 e=a(·,·) is our energy norm and C∈R+. We select this definition of a(·,·) by bounding from above the convective term of b(·,·) using a meshindependent constant Cfor the Poincar´e inequality that also includes the effect of σ1 2. 5.3.1.1. Energy-norm adaptivity For this example, we consider ε= 10−3as the diffusive coefficient, σ= (3,1)T, and Ω = (0,1)2. The load function fis a linear continuos form on Hand it is selected so that the solution uis of the form: u(x, y) = eε x(x−1) cosh 5001 2+σ−1(x, −y)−2 .(5.13) Figure 5.10 shows the solution of this convection-dominated diffusion example. The initial uniform mesh is composed of thirty-six root elements. Figure 5.11 1It is essential to consider a mesh-independent norm a(·,·)1/2since we approximate some errors by computing the difference of the norm of two approximated solutions evaluated on different grids. 2The actual value of the constant C is unneeded in practice since we compute relative error indicators; in our case, we select (C+ε) = 1. 39 5. 2D Numerical results for hp-adaptivity shows the final energy-norm hand hp-adapted meshes and the evolution of erel. As expected, we observe heavy h-refinements around the line that characterizes the solution. In the hp-adapted case, we also observe an increase in the polynomial order in some of the elements near this characteristic line. We also observe exponential convergence rates (see Figure 5.11d). 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x y −1.4·10−40.5 1 Value of the solution Figure 5.10.: Solution of the convection-dominated diffusion example 1. 5.3.2. Convection-dominated diffusion: example 2 We now consider a more challenging setting with advection skew to the mesh. We solve a similar problem to the one depicted in Figure 9.3 of [55] (see Figure 5.12). Our convection-dominated diffusion problem is governed by eq. (5.11) on the domain Ω = (0,1)2, with ε= 10−4,σ= (cos θ, sin θ)T,θ= arctan(2), and zero Dirichlet boundary conditions, as depicted in Figure 5.12a. We define our source term f(with support in Ωfand illustrated in Figure 5.12b) 40 5. 2D Numerical results for hp-adaptivity 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x y −8·10−40.1 0.21 Solution (a) Solution at iteration 17. 0 0.2 0.4 0.6 0.8 10 0.2 0.4 0.6 0.8 1 x y 1234567891011 Order of approximation (b) hp-adapted mesh at iteration 17. 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x y −5.2·10−40.1 0.21 Solution (c) Solution at iteration 21. 0 0.2 0.4 0.6 0.8 10 0.2 0.4 0.6 0.8 1 x y 1234567891011 Order of approximation (d) hp-adapted mesh at iteration 21. Figure 5.16.: Numerical solutions and hp-adapted meshes (polynomial orders in the x-direction) at iterations 17 and 21. 47 5. 2D Numerical results for hp-adaptivity 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x y −1.86 ·10−30.13 0.27 Direct solution (a) Solution to the direct problem. 0 0.2 0.4 0.6 0.8 10 0.2 0.4 0.6 0.8 1 x y −5·10−52.4 4.8 Adjoint solution (b) Solution to the adjoint problem. Figure 5.17.: Direct and adjoint numerical solutions of the convection-dominated diffusion problem for GO adaptivity. 48 5. 2D Numerical results for hp-adaptivity 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x y 1234567891011 Order of approximation (a) Final hp-adapted mesh with polynomial orders in the x-direction. 0 0.2 0.4 0.6 0.8 10 0.2 0.4 0.6 0.8 1 x y 1234567891011 Order of approximation (b) Final hp-adapted mesh with polynomial orders in the y-direction. 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x y (c) Final h-adapted mesh, p= 1. 101102103104105 10−7 10−5 10−3 10−1 101 103 Number of DoFs, N(log scale) Relative error in % (log scale) hp (p+ 2) h(p= 1) h(p= 2) (d) Evolution of eQoI rel in the process. Figure 5.18.: Final hand hp-adapted meshes for our second convectiondominated diffusion example and the evolution of eQoI rel . 49 6. 3D Numerical results for hp-adaptivity 6.1. Wave propagation problem Let us consider the following non-elliptic problem based on heterogeneous Helmholtz’s equation. Find usuch that, −∇ · (σ∇u)−k2u=1Ωfin Ω,(6.1) u= 0 on ΓD,(6.2) ∇u·n = 0 on ΓN,(6.3) where Ω = (0,1)3⊂R3, Ωf=0,1 43⊂Ω, and k= (4 ·2π, 2π). ΓDand ΓN stand for the parts of the boundary ∂Ω where we impose homogeneous Dirichlet and Neumann boundary conditions, respectively. We impose Dirichlet boundary conditions on the 3 faces whose intersection is (0,0,0) and Neumann boundary on the 3 faces whose intersection is (1,1,1). ΓD:= ([0,1] ×[0,1] × {0})∪([0,1] × {0} × [0,1]) ∪({0} × [0,1] ×[0,1]) ,(6.4) ΓN:= ((0,1) ×(0,1) × {1})∪((0,1) × {1} × (0,1)) ∪({1} × (0,1) ×(0,1)) . (6.5) Here, σ(x) =          1 if x∈Ω1=0<x<1,0< y < 1 2,0< z < 1, 103if x∈Ω2=1 2<x<1,1 2< y < 1,0< z < 1 2, 10 if x∈Ω3=1 2<x<1,1 2< y < 1,1 2< z < 1, 10−2if x∈Ω4=0<x<1 2,1 2< y < 1,0< z < 1. We define the operators b(·,·) and a(·,·) associated with the above problem as follows: b(·,·):=⟨∇· , σ∇·⟩L2(Ω) −k2⟨· ,·⟩L2(Ω) , a (·,·):=⟨∇· , σ∇·⟩L2(Ω)+k2⟨· ,·⟩L2(Ω). (6.6) 50 6. 3D Numerical results for hp-adaptivity Once again, ∥·∥2 e=a(·,·) is our energy norm and |b(ϕ, ψ)| ≤ |a(ϕ, ψ)|,∀ϕ, ψ ∈ H. Figure 6.1 displays the different materials in the domain. Following the definition of eq. (5.3), we select Ωl=3 4,13⊂Ω. For Goal-Oriented (GO) adaptivity, Figures 6.2a and 6.2b show the solutions of the direct and adjoint problems, respectively. 0.01 1 10 1000 Diffusivity of the materials Figure 6.1.: Diffusive coefficient values for the different materials in the domain. 6.1.1. Energy-norm adaptivity Figure 6.3 displays the final hp-adapted meshes for our 3D wave propagation example in a lossy medium using energy-norm adaptivity. The initial uniform mesh is composed of sixty-four root elements. As expected, we observe heavy h-refinements near different materials’ interfaces. Figure 6.4 shows the corresponding convergence curves. As in the 2D case, the energy-norm hp-adaptivity provides proper convergence results in terms of energy. However, the convergence of the energy-norm adaptivity in terms of the error in the Quantity of Interest (QoI) is slow, especially in the pre-asymptotic regime. 6.1.2. Goal-Oriented adaptivity Figure 6.5 displays the final hp-adapted meshes for our 3D wave propagation example in a lossy medium using GO adaptivity. The initial uniform mesh is com51 6. 3D Numerical results for hp-adaptivity 0 2.1 4.19 ·10−3 Direct solution (a) Solution to the direct problem. 0 8.65 ·10−20.17 Adjoint solution (b) Solution to the adjoint problem. Figure 6.2.: Absolute value of the direct and adjoint solutions of our 3D wave propagation example in a lossy medium. posed of sixty-four root elements. As expected, we observe heavy h-refinements near different materials’ interfaces. When using GO adaptivity, the evolution of the error in the QoI exhibits much better behavior, while the energy convergence becomes suboptimal, as expected. As computational problems grow in complexity and scale, they pose significant challenges to our computational capabilities. Developing parallel computational strategies for finite element discretization schemes [21, 150] solves these challenges. By distributing tasks and computations across multiple processors or computational nodes, these strategies can address complex engineering problems, increasing computational capacity and improving efficiency. This dissertation follows an algorithm [99] that distributes the computational domain among participating processes. It subdivides the domain into sub-domains of relatively equal computational cost and assigns them to different processes, thus optimizing resource utilization. As the algorithm progresses, dynamic rebalancing techniques are employed to redistribute tasks, ensuring optimal load distribution across processes. Introducing adaptivity creates the challenge of balancing computational workload—techniques such as limiting refinements and process aggregation address this. Additionally, communication efficiency is supported by aggregating data into larger sets, minimizing the frequency and latency of data transfers. 52 6. 3D Numerical results for hp-adaptivity 1234567891011 Order of approximation (a) Final hp-adapted mesh with polynomial orders pin the x-direction. 1234567891011 Order of approximation (b) Final hp-adapted mesh with polynomial orders pin the y-direction. 1234567891011 Order of approximation (c) Final hp-adapted mesh with polynomial orders pin the z-direction. Figure 6.3.: Energy-norm adaptivity. Final hp-adapted meshes for our 3D wave propagation example in a lossy medium. 53 6. 3D Numerical results for hp-adaptivity 102103104105 10−1 100 101 102 Number of DoFs, N Relative error in % (log scale) GOA energy-norm (a) Evolution of goal-oriented adaptivity. 102103104105 10−7 10−5 10−3 10−1 101 Number of DoFs, N GOA energy-norm (b) Evolution of energy-norm adaptivity. Figure 6.4.: Convergence history of eQoI rel and ˜eenergy rel for the energy-norm and GO hp-adaptive strategies. 54 6. 3D Numerical results for hp-adaptivity 1 2 3 4 5 6 7 8 9 Order of approximation (a) Final hp-adapted mesh with polynomial orders pin the x-direction. 123456789 Order of approximation (b) Final hp-adapted mesh with polynomial orders pin the y-direction. 123456789 Order of approximation (c) Final hp-adapted mesh with polynomial orders pin the z-direction. Figure 6.5.: GO adaptivity. Final hp-adapted meshes for our 3D wave propagation example in a lossy medium. 55 Part II. Goal-Oriented hp-adaptivity for parametric PDEs. 56 7. Database generation for DL inversion Algorithm 4: Multi-Adaptive Goal-Oriented adaptive process Input: PDE, initial finite element mesh, Ssamples of model parameters, definition of the QoI Output: A final hp-adapted mesh while eq. (7.12) is not satisfied do Perform a global refinement (user-defined); while error indicators above threshold do for each sample miwhere i= 1 to Sdo Solve the forward problem for sample miusing eq. (7.8); Solve the adjoint problem for sample miusing eq. (7.9); Calculate error indicators for the i-th sample using eq. (7.7); end Compute the error indicators using eq. (7.11), which combines those from all samples into a single measure; Remove basis functions with low error indicators to unrefine the mesh; end Update error in the QoI; end 63 7. Database generation for DL inversion 7.2. Generation of databases We refer to the two main stages of the database production process as the Adaptive and Generation processes. In the Adaptive process, which is the initial stage, we construct a highly accurate hp-grid, referred to as the adapted mesh, capable of accommodating an arbitrary number of samples (SA). Subsequently, in the Generation process, we employ this adapted mesh to solve multiple FEM problems, generating the required data. This approach’s key aspect is using a single hp-FEM, i.e., the adapted mesh, throughout the entire Generation process. To elaborate further, we begin by generating a reduced number of samples of model parameters (SA) that parametrize the PDE of the problem. We then construct the adapted mesh, ensuring it satisfies eq. (7.12), where the maximum error for all the SAsamples falls below a user-prescribed tolerance. The central concept behind this approach is the anticipation of achieving low errors when employing the adapted mesh with samples different from those used in the adaptive process. In the subsequent Generation process, we solve one FEM problem for each sample in SG, which represents a set of additional samples we consider. This process allows us to obtain accurate synthetic data or measurements database. The overall process can be summarized as follows: 1. Adaptive process: a) Generate SAsamples of model parameters to be used in the Adaptive process. b) Construct the adapted grid by employing the hp-FEM following the guidelines described in Section 7.1. 2. Generation process: a) Generate SGadditional samples of model parameters specifically for the Generation process. b) For each sample in SG, solve a FEM problem using the adapted mesh, which was specifically designed during the Adaptive process to deliver highly accurate solutions for a wide range of model parameters. By following this approach, we can efficiently generate a reliable database of accurate synthetic data or measurements using a single adapted mesh for multiple samples, thereby reducing computational expenses and maintaining high accuracy across various scenarios. 7.2.1. Computational costs of MAGO The GOA strategy, elaborated upon in Part I of this dissertation and presented in Algorithm 3, comprises a series of refining and coarsening steps. We use a 64 7. Database generation for DL inversion direct solver to solve each FEM problem, contributing to the computational cost of building the hp-mesh. The total cost is given by: CGOA = NrIter X i=1 NrCoarse X j=1 Cas Nij +Can Nij +Cfa Nij + 2 Cso Nij +Ces Nij.(7.13) Each component of the cost corresponds to specific operations: Cas for assembling the matrix, Can for the analysis part of the direct solver, Cfa for factorization, Cso for solving the linear system of equations after factorization (i.e., backward elimination), and Ces for computing the error estimators. Additionally, Nij represents the number of Degrees of Freedom (nDoF) of the meshes at each iteration iof the adaptive process and each coarsening step jassociated with each iteration i. Notably, the factor of 2 in the Cso term accounts for solving both the forward and adjoint problems. Since we use a direct solver, the extra cost for solving the adjoint problem associated with the forward problem reduces to only backward and forward substitutions. The costs related to the finest grid with NfDegrees of Freedom (DoF) dominate those associated with the coarser grids. Consequently, we can approximate Equation (7.13) by: CGOA (Nf)≈Cas (Nf) + Can (Nf) + Cfa (Nf)+2Cso (Nf) + Ces (Nf).(7.14) Thus, the approximate costs of generating a database of SGsamples with the Single-Adaptive Goal-Oriented (SAGO) and the MAGO approaches are as follows: SAGO approach: We approximate the computational cost CSAGO of generating one GOA mesh for each of the SGsamples by: CSAGO = SG X i=1 CGOA N(i) f.(7.15) MAGO approach: The cost CMAGO of generating the database with the MAGO strategy is the sum of the costs of constructing the adapted mesh CAplus the cost of actually generating the data CG, that is, CMAGO =CA+CG. 65 7. Database generation for DL inversion We approximate the cost of the Adaptive process with SAsamples by: CA=SACGOA (Nmago f),(7.16) where Nmago frepresents the nDoF in the fine mesh adapted using the MAGO strategy. After generating the adapted coarse mesh of size Nmago c, we proceed to generate the data. The approximate cost of the Generation process is then given by: CG=SGCas (Nmago c) + Can (Nmago c) + Cfa (Nmago c)+2Cso (Nmago c).(7.17) In our MAGO approach, all samples share the same discretization, which conforms with the material parameters of the PDE. This design allows us to precompute and reuse certain information across different samples. As detailed in the following subsection, we perform precomputations for the integrals of the global matrices and the analysis part of the direct solver of equations, thereby reducing the assembling and analysis processes to a single occurrence. Considering that the cost of computing the estimators is comparable to the cost of assembling [146], we obtain the following approximations: CA= 2 Cas (Nmago f) + Can (Nmago f) + SACfa (Nmago f)+2Cso (Nmago f),(7.18) and CG=Cas (Nmago c) + Can (Nmago c) + SGCfa (Nmago c)+2Cso (Nmago c).(7.19) Compared to the SAGO approach, the costs associated with assembling, analysis, and estimation occur only once due to the precomputations, contributing to improved efficiency and reduced computational overhead. Factorization costs dominate other aspects in traditional C0-continuous FEM problems when using a direct solver. It scales as (see, e.g., [52, 53]): ON(1+(d−1)/2),(7.20) where d= 1,2,3 represents the dimension of the problem and Ndenotes the nDoF. The final approximate costs for the SAGO and MAGO approaches are as follows: CSAGO ≈ SG X i=1 Cfa Ni f≈SGCfa (Nsago f),(7.21) 66 7. Database generation for DL inversion where Nsago frepresents an average value (using eq. (7.20)) of the nDoF of the SG fine grids, and CMAGO =SACfa (Nmago f)+SGCfa (Nmago c).(7.22) While SA≪SG, the relationship between Nmago f,Nmago c, and Nsago fis not generalizable. This relationship depends on various factors, including the initial mesh configuration, the refinement and coarsening criteria, the convergence behavior of the solution, and the problem’s complexity. In some instances, investing resources in building a sufficiently good adapted mesh with only a fraction of samples is reasonable. The gains in the Generation part of the strategy will likely compensate for this computational effort, considering the specific problem. Consequently, we expect that CMAGO < CSAGO. However, it is essential to note that the number of samples SAand the mesh size significantly impact the accuracy of the generated data. Thus, a tradeoff exists between accuracy and the cost of obtaining the adapted mesh. 7.2.2. Precomputations of the global matrices In many adaptive FEM implementations, the integrals associated with the bilinear form and error indicators are calculated element by element. However, when dealing with many material samples S(possibly reaching hundreds of thousands), computing all these integrals for each sample becomes computationally expensive. To overcome this challenge, we take advantage of our materials being piecewise constant and conforming to discretization. We perform a clever optimization by precomputing and storing the integrals for an arbitrary unitary sample, where material properties are assumed to equal one. Once these integrals are precomputed, we can reuse this information across all samples instead of recalculating the integrals for each sample. This technique significantly accelerates the integration process and reduces the computational cost of handling many material samples. The bilinear form in the problem may consist of multiple terms, denoted as: bmi(·,·) = Mb X j=1 bmi j(·,·),(7.23) where Mbrepresents the number of terms in the bilinear form. For each material sample mi, we compute the contributions associated with each element Kas follows: bmi(·,·)K= Mb X j=1 mj i(K)b1 j(·,·)K,(7.24) 67 7. Database generation for DL inversion Algorithm 5: Precomputation of element-wise matrices Input: Given variational formulation Output: Pre-computed element-wise unitary matrices for Each term in the bilinear form (j= 1, . . . , Mb)do for Each element in the finite element discretization (K∈ T )do Compute and store the element-wise unitary matrix b1 j(·,·)K; end end Algorithm 6: Precomputation: assembling the global matrices Input: Variational formulation of the problem, Ssamples, pre-computed matrices Output: Assembled global matrices for all samples for Each sample mi(i= 1, . . . , S)do for Each term in the bilinear form (j= 1, . . . , Mb)do for Each element in the finite element discretization (K∈ T )do Initialize the element-wise matrix to zero: [bmi(·,·)]K= 0; Load the pre-computed element-wise unitary matrix b1 j(·,·)K(see Algorithm 5); Load the material property of the element (mj i(K)); Update the value of the element matrix [bmi(·,·)]K= [bmi(·,·)]K+mj i(K)b1 j(·,·)K; Assemble the global matrix by inserting [bmi(·,·)]Kinto bmi(·,·); end end We obtain the fully assembled matrix associated with the i-th sample. end 68 7. Database generation for DL inversion where mj i(K) represents the scalar material property associated with the Kth element of the i-th material sample and the j-th term of the bilinear form. Additionally, b1 j(·,·)Kcorresponds to the j-th unitary element-wise bilinear term (sub-matrix) associated with the element K. Specifically, mj i(K) = 1 for the bilinear terms independent of the material properties. To optimize the computation process, we pre-compute and store all the unitary integrals in b1 j(·,·)Kfor all elements in the discretization, K∈ T , and for j= 1, . . . , Mb. This way, we only need to compute these integrals once and then load the pre-computed unitary sub-matrices for each of the Ssamples. By multiplying them with the corresponding material property for each element, we can efficiently assemble the global bilinear matrices and compute the error indicators as scalar products according to eq. (7.10). The overall process is summarized using the following algorithms: Algorithm 5 shows how we compute and store the unitary element-wise matrices; Algorithms 6 and 7 explain how we construct the global matrices and compute the error indicators, respectively, utilizing the pre-computed information. Algorithm 7: Precomputation: computation of the error indicators Input: Error indicator, variational formulation of the error, Ssamples, pre-computed matrices, forward and adjoint solutions Output: Error indicator for all samples at the same time for Each sample mi(i= 1, . . . , S)do for Each term in the bilinear form (j= 1, . . . , Mb)do for Each element in the finite element discretization (K∈ T )do Initialize the element-wise matrix to zero: [ami(·,·)]K= 0; Load the pre-computed element-wise unitary matrix a1 j(·,·)Kfrom Algorithm 5; Load the material property of the element as mj i(K); Update the element matrix as [ami(·,·)]K= [ami(·,·)]K+mj i(K)a1 j(·,·)K; end end Compute the error indicators for the i-th sample using eq. (7.7); end Compute a single error indicator considering all the samples according to eq. (7.11); Notably, the precomputation of global matrices can be utilized in both the Adaptive and Generation parts, enhancing the efficiency of the adaptive FEM. 69 7. Database generation for DL inversion 7.3. Numerical results This section demonstrates the performance of the MAGO approach in generating large databases for various problems. It highlights the method’s capability to adaptively construct meshes and compute the QoI across a broad spectrum of sample configurations. The primary objective is to design an optimized single hp-mesh, whose size is as small as possible, to efficiently determine the QoI for all samples. These QoI are denoted as l(umi) for i= 1, . . . , S. 7.3.1. Definitions We categorize the numerical results into three categories based on the purpose they serve: 1. To demonstrate the quasi-exponential convergence of the MAGO strategy, showcasing how the MAGO approach achieves rapid convergence in adaptive mesh generation. 2. To verify the accuracy of the produced measurements, assess the precision of the computed QoI values obtained using the MAGO approach. 3. To highlight the numerical advantages of the MAGO approach, quantifying its benefits in terms of computational efficiency and mesh size reduction. The problems considered include 2D scenarios involving Poisson and Helmholtz equations. The Hilbert space Hchosen for all problems is defined as H= {u∈H1(Ω) |u= 0 on ΓD}, where ΓDdenotes the boundary with Dirichlet boundary conditions. The mesh is designed specifically for all scenarios’ materials, sources, and desired QoI. 7.3.1.1. Convergence of the MAGO adaptivity We investigate the convergence behavior of the MAGO approach by varying the number of samples SAused to construct adapted meshes. This follows the procedure detailed in Item 1 of Section 7.2. We present visualizations of the final adapted hp-meshes for various cases and introduce quasi-exponential convergence curves, showcasing the effectiveness of the MAGO strategy. To calculate convergence, we calculate two relative errors in the QoI concerning the nDoF during the adaptive processes: the maximum relative error emax rel and the mean relative error emean rel . 70 7. Database generation for DL inversion These errors are computed among the SAsamples, where the value of SAvaries based on specific examples. The maximum relative error is given by: emax rel = max i=1,...,S  l(umi)−l(umi T) l(umi) ·100.(7.25) The mean relative error is: emean rel =1 S S X i=1  l(umi)−l(umi T) l(umi) ·100.(7.26) In these equations, umiand umi Trepresent the solutions linked to a given model mion a fine and a coarser mesh, respectively. The maximum relative error highlights the worst-case error among samples, whereas the mean relative error offers an average accuracy overview. The maximum, mean, and traditional (erel) relative errors are identical for a single sample. 7.3.1.2. Computational costs of generating the database Equations (7.21) and (7.22) provide an approximate estimation of the computational costs for the SAGO and MAGO approaches, respectively, where the cost is influenced by the nDoF, as articulated in Equation (7.20). Consequently, we approximate the computational cost CSAGO of generating one GOA mesh for each of the SGsamples by: CSAGO ≈ SG X i=1 Nsago f(i)(1+(d−1)/2) .(7.27) In addition, we provide an approximate estimation of the computational costs for the MAGO by: CMAGO ≈ SA X i=1 Nmago f(i)(1+(d−1)/2) + SG X i=1 Nsago c(i)(1+(d−1)/2) .(7.28) While we may omit certain additional costs for clarity, it is essential to highlight that generating the database using the MAGO approach yields significant savings, as detailed in Section 7.2. 7.3.1.3. Generating model parameter samples We generate all model parameter samples randomly. Specifically, we use a uniform distribution over the interval [−1,3] to derive the values of log10 (mi). Consequently, the values of σican vary by up to four orders of magnitude, ranging from 10−1to 103. 71 7. Database generation for DL inversion 7.3.2. Wave propagation example We consider the following non-elliptic problem based on Helmholtz’s equation. Find usuch that −∇ · (∇u)−jσ(x)u= 1 in Ω,(7.29) u= 0 on ∂Ω.(7.30) 7.3.2.1. Example: wave propagation problem We employ a 5 ×5 grid for numerical computations, encompassing a square computational domain, Ω = [0,1]2. Within this domain, we distinguish two regions, Ωfand Ωl, which symbolize the source and the QoI, respectively. These regions are defined within Ω as Ωf=1 20,3 202and Ωl=17 20,19 202, with the origin situated at the bottom-left corner of the domain. Figure 7.2 shows the computational domain Ω, its boundary ∂Ω (subject to Dirichlet conditions), and the locations of Ωland Ωf. In this depiction, Ωldefines the region of the QoI function l(ϕ), whereas Ωfis the region for the source function. Ω Ωf Ωl ∂Ω Figure 7.2.: Our grid-based domain example is defined over the domain Ω. The Dirichlet boundary condition is denoted by ∂Ω. The source function is supported on Ωf, and the QoI l(ϕ) is supported on Ωl. Figures 7.3a and 7.3b show the absolute values of forward and adjoint numerical solutions on a logarithmic scale. 72 7. Database generation for DL inversion 7.3.2.3. Wave propagation: accuracy To showcase the efficiency of the MAGO approach, we present statistical properties of errors via standard box plots. These plots offer a detailed depiction of the error distribution. Tukey [186] introduced these box plots in 1977 to provide a robust data representation. The box plots visually represent how the maximum relative error in the QoI varies with different numbers of training samples for the adaptive process. As shown in Figure 7.9, the box plots represent various values of SAin the MAGO process. We consider the adaptive grids when they reach a maximum relative error, emax rel that drops below 10−5. Every number over each upper whisker represents the nDoF in each hp grid with the maximum relative error reduced to under 10−5. The trend suggests that as the number of training samples for the adaptive process increases, the maximum relative error in the QoI tends to decrease or remain stable. The variation in relative errors becomes more confined with increasing training samples, as indicated by the tightening spread of the box plots. 1 5 10 100 1000 10−8 10−7 10−6 10−5 10−4 10−37,779 15,440 16,86715,87017,477 Number of training samples for the adaptive process Max relative error in the QoI (%) Figure 7.9.: Box plots for different adaptive grids with a threshold maximum relative error set at 10−5. 7.3.2.4. Wave propagation: computational costs We estimate the computational cost based on the factorization cost, which constitutes the most resource-intensive part of data generation and a significant expense in many FEM codes. In Tables 7.1 and 7.2, we approximate the computational 79 7. Database generation for DL inversion expenses in terms of Floating Point Operations (FLOPs) associated with database generation using the SAGO and MAGO approaches. We compute CSAGO using Equation (7.27). The maximum relative error is constrained to be under 10−5. We observe that the MAGO approach is more cost-effective for solving problems in this case, as CMAGO < CSAGO. The median values of the maximum relative error are below 10−6, indicating a significantly higher accuracy by one order of magnitude in the results, and demonstrated in Figure 7.9. This increased accuracy is achieved at a reduced cost, making the MAGO approach particularly suitable for solving challenging problems. Number of DoF CSAGO SG Nsago f1051071091011 41259 8.4238 ·1011 8.4238 ·1013 8.4238 ·1015 8.4238 ·1017 Table 7.1.: The computational cost based on the factorization cost of generating the database using the SAGO strategy. Number of DoF CMAGO SG SANmago cNmago f1051071091011 5 15440 52351 1.9191 ·1011 1.9185 ·1013 1.9185 ·1015 1.9185 ·1017 10 16867 57171 2.1919 ·1011 2.1906 ·1013 2.1906 ·1015 2.1906 ·1017 100 15870 53661 2.0117 ·1011 1.9994 ·1013 1.9992 ·1015 1.9992 ·1017 1000 17477 60381 2.4588 ·1011 2.3120 ·1013 2.3105 ·1015 2.3105 ·1017 Table 7.2.: The computational cost based on the factorization cost of generating the database using the MAGO strategy. 80 7. Database generation for DL inversion 7.3.3. Poisson example We consider the following elliptic problem based on the Poisson equation. Find usuch that −∇ · (σ(x)∇u) = 1 in Ω,(7.31) u= 0 on ∂Ω.(7.32) 7.3.3.1. Example: cross-shaped domain Poisson problem We address a Poisson problem over a domain Ω in a two-dimensional space, represented on a 5 ×5 grid. The domain Ω resembles a cross and is defined by Ω = [0,1] ×1 5,4 5∪1 5,4 5×[0,1]. Please refer to Figure 7.10 to visualize the domain. Within this domain, there are two notable regions: Ωf, the source area, and Ωl, the QoI area. Both Ωfand Ωlare subregions within Ω. Specifically, Ωfis the square defined by x∈1 5,2 5and y∈1 5,2 5, and Ωlis the square with x∈3 5,4 5and y∈3 5,4 5. The origin of the coordinate system is the bottom-left corner of the domain. Figures 7.11a and 7.11b showcase the absolute values of forward and adjoint numerical solutions, respectively, on a logarithmic scale. Ω Ωf Ωl ∂Ω Figure 7.10.: Computational domain Ω, where homogeneous Dirichlet boundary conditions are imposed on ∂Ω. Additionally, we define Ωlas the support of the QoI l(ϕ), and Ωfas the support of the source function. 81 7. Database generation for DL inversion 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x y 0 1.73 3.46 ·10−3 Direct solution (a) Solution to the direct problem. 0 0.2 0.4 0.6 0.8 10 0.2 0.4 0.6 0.8 1 x y 0 0.73 1.46 ·10−3 Adjoint solution (b) Solution to the adjoint problem. Figure 7.11.: Absolute value of the solutions of our cross-shaped domain Poisson example. 82 7. Database generation for DL inversion 7.3.3.2. Cross-shaped domain Poisson: convergence We present the numerical results for our MAGO process with one sample in Figure 7.12, and for SAequal to 5, 10, 50, and 100 in Figures 7.13 to 7.16, respectively. The convergence is quasi-optimal, and the three norms (l1,l2, and l∞) yield similar results as in previous examples. Predictably, regions with significant material coefficient fluctuations witness more substantial mesh refinement, resulting in more compact mesh elements (denoted as h) near intersections of multiple materials. In addition, Figure 7.12c represents the material properties. The quasi-optimal exponential convergence graph for this scenario, along with the progress of emax rel and emean rel , can be found in Figure 7.12d. 83 7. Database generation for DL inversion 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x y 1234567891011 Order of approximation (a) Final hp-adapted mesh with polynomial orders pin the x-direction. 0 0.2 0.4 0.6 0.8 10 0.2 0.4 0.6 0.8 1 x y 1234567891011 Order of approximation (b) Final hp-adapted mesh with polynomial orders pin the y-direction. 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x y 0.1 500.05 1,000 Value of σ(x) (log scale) (c) Values for the materials in the domain. 101102103104105 10−7 10−5 10−3 10−1 101 Number of DoFs, N(log scale) Max relative error in % (log scale) l∞ (d) Evolution of emax rel in the process. Figure 7.12.: hp-adapted meshes for our 1-sample cross-shaped domain example. 84 7. Database generation for DL inversion 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x y 1234567891011 Order of approximation (a) Final hp-adapted mesh with polynomial orders pin the x-direction. 0 0.2 0.4 0.6 0.8 10 0.2 0.4 0.6 0.8 1 x y 1234567891011 Order of approximation (b) Final hp-adapted mesh with polynomial orders pin the y-direction. 101102103104105 10−1 100 101 102 Number of DoFs, N(log scale) Max relative error in % (log scale) l∞l1l2 (c) Evolution of emax rel in the process. 101102103104105 10−1 100 101 Number of DoFs, N(log scale) Mean relative error in % (log scale) l∞l1l2 (d) Evolution of emean rel in the process. Figure 7.13.: hp-adapted meshes for our 5-sample cross-shaped domain example. 85 7. Database generation for DL inversion 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x y 1234567891011 Order of approximation (a) Final hp-adapted mesh with polynomial orders pin the x-direction. 0 0.2 0.4 0.6 0.8 10 0.2 0.4 0.6 0.8 1 x y 1234567891011 Order of approximation (b) Final hp-adapted mesh with polynomial orders pin the y-direction. 101102103104105 10−1 100 101 102 Number of DoFs, N(log scale) Max relative error in % (log scale) l∞l1l2 (c) Evolution of emax rel in the process. 101102103104105 10−1 100 101 Number of DoFs, N(log scale) Mean relative error in % (log scale) l∞l1l2 (d) Evolution of emean rel in the process. Figure 7.14.: hp-adapted meshes for our 10-sample cross-shaped domain example. 86 7. Database generation for DL inversion 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x y 1234567891011 Order of approximation (a) Final hp-adapted mesh with polynomial orders pin the x-direction. 0 0.2 0.4 0.6 0.8 10 0.2 0.4 0.6 0.8 1 x y 1234567891011 Order of approximation (b) Final hp-adapted mesh with polynomial orders pin the y-direction. 101102103104105 100 101 102 Number of DoFs, N(log scale) Max relative error in % (log scale) l∞l1l2 (c) Evolution of emax rel in the process. 101102103104105 10−1 100 101 Number of DoFs, N(log scale) Mean relative error in % (log scale) l∞l1l2 (d) Evolution of emean rel in the process. Figure 7.15.: hp-adapted meshes for our 50-sample cross-shaped domain example. 87 7. Database generation for DL inversion 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x y 1234567891011 Order of approximation (a) Final hp-adapted mesh with polynomial orders pin the x-direction. 0 0.2 0.4 0.6 0.8 10 0.2 0.4 0.6 0.8 1 x y 1234567891011 Order of approximation (b) Final hp-adapted mesh with polynomial orders pin the y-direction. 101102103104105 100 101 102 Number of DoFs, N(log scale) Max relative error in % (log scale) l∞l1l2 (c) Evolution of emax rel in the process. 101102103104105 10−1 100 101 Number of DoFs, N(log scale) Mean relative error in % (log scale) l∞l1l2 (d) Evolution of emean rel in the process. Figure 7.16.: hp-adapted meshes for our 100-sample cross-shaped domain example. 88 7. Database generation for DL inversion 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x y 1234567891011 Order of approximation (a) Final hp-adapted mesh with polynomial orders pin the x-direction. 0 0.2 0.4 0.6 0.8 10 0.2 0.4 0.6 0.8 1 x y 1234567891011 Order of approximation (b) Final hp-adapted mesh with polynomial orders pin the y-direction. 102103104105 100 101 102 Number of DoFs, N(log scale) Max relative error in % (log scale) l∞l1l2 (c) Evolution of emax rel in the process. 102103104105 10−1 100 101 Number of DoFs, N(log scale) Mean relative error in % (log scale) l∞l1l2 (d) Evolution of emean rel in the process. Figure 7.22.: hp-adapted meshes for our 50-sample grid-based domain example. 95 7. Database generation for DL inversion 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x y 1234567891011 Order of approximation (a) Final hp-adapted mesh with polynomial orders pin the x-direction. 0 0.2 0.4 0.6 0.8 10 0.2 0.4 0.6 0.8 1 x y 1234567891011 Order of approximation (b) Final hp-adapted mesh with polynomial orders pin the y-direction. 102103104105 100 101 102 Number of DoFs, N(log scale) Max relative error in % (log scale) l∞l1l2 (c) Evolution of emax rel in the process. 102103104105 10−1 100 101 Number of DoFs, N(log scale) Mean relative error in % (log scale) l∞l1l2 (d) Evolution of emean rel in the process. Figure 7.23.: hp-adapted meshes for our 100-sample grid-based domain example. 96 7. Database generation for DL inversion 7.3.3.7. Grid-based domain Poisson: accuracy Figure 7.24 displays a series of box plots corresponding to different adaptive grids. The x-axis denotes the number of samples (SA) used during the MAGO process to formulate the final hp grid. We stop every MAGO adaptation once the maximum relative error (emax rel ) is reduced to under 1.0%. Every number over each upper whisker represents the nDoF in each hp grid with the maximum relative error reduced to under 1.0%. The general trend suggests that as the number of training samples for the adaptive process increases, the maximum relative error in the QoI tends to decrease. The spread of relative errors becomes greater with more training samples. However, the median error values do not change dramatically after fifty training samples, suggesting decreasing returns in error reduction with additional samples. 1 5 10 50 100 10−5 10−4 10−3 10−2 10−1 100 101 102832 18,77738,407 73,48781,573 Number of training samples for the adaptive process Max relative error in the QoI (%) Figure 7.24.: Box plots for different adaptive grids with a threshold maximum relative error set at 1.0% 7.3.3.8. Grid-based domain Poisson: computational costs In Tables 7.5 and 7.6, we approximate the computational expenses in terms of FLOPs associated with database generation using the SAGO and MAGO approaches. We compute CSAGO using Equation (7.27). The maximum relative error is constrained to be under 1.0%. While solving problems using the SAGO approach appears to be cost-effective, it is essential to acknowledge a limitation, as CSAGO < CMAGO. 97 7. Database generation for DL inversion Number of DoF CSAGO SG Nsago f1051071091011 28428 6.4888 ·1011 6.4888 ·1013 6.4888 ·1015 6.4888 ·1017 Table 7.5.: The computational cost based on the factorization cost of generating the database using the SAGO strategy. Number of DoF CMAGO SG ANmago cNmago f1051071091011 5 18777 186161 2.5770 ·1011 2.5730 ·1013 2.5730 ·1015 2.5730 ·1017 10 38407 428635 7.5550 ·1011 7.5272 ·1013 7.5269 ·1015 7.5269 ·1017 50 73487 769965 2.0259 ·1012 1.9925 ·1014 1.9921 ·1016 1.9921 ·1018 100 81573 831627 2.4056 ·1012 2.3306 ·1014 2.3298 ·1016 2.3298 ·1018 Table 7.6.: The computational cost based on the factorization cost of generating the database using the MAGO strategy. 98 Part III. Main achievements, conclusions and future work 99 8. Main Achievements 8.1. Peer-reviewed Publications 2023 F. V. Caro, V. Darrigrand, J. Alvarez-Aramberri, and D. Pardo. A MultiAdaptive-Goal-Oriented Strategy to Generate Massive Databases of Parametric PDEs. To be submitted to Computer Methods in Applied Mechanics and Engineering in October 2023. 2022 F. V. Caro, V. Darrigrand, J. Alvarez-Aramberri, E. Alberdi, and D. Pardo. A painless multi-level automatic goal-oriented hp-adaptive coarsening strategy for elliptic and non-elliptic problems.Computer Methods in Applied Mechanics and Engineering, 401:115641, 2022. Impact Factor: 7.2, Quartile: Q1, Scimago Ranking. https://doi.org/10.1016/j.cma.2022.115641 2022 F. V. Caro, V. Darrigrand, J. Alvarez-Aramberri, E. A. Celaya, and D. Pardo. 1D Painless Multi-level Automatic Goal-Oriented hand pAdaptive Strategies Using a Pseudo-Dual Operator. In Computational Science – ICCS 2022, pages 347–357, 2022. https://doi.org/10.1007/978-3-031-08754-7_43 8.2. International Conferences 2023 F. V. Caro, V. Darrigrand, J. Alvarez-Aramberri, and D. Pardo. Generation of Massive Databases for Deep Learning Inversion Using A GoalOriented hp-Adaptive Strategy. XI International Conference on Adaptive Modeling and Simulation, Gothenburg, Sweden, [June 19-21, 2023]. 2022 F. V. Caro, V. Darrigrand, J. Alvarez-Aramberri, E. Alberdi, and D. Pardo. A Painless Automatic hp-Adaptive Coarsening Strategy For NonSPD problems: A Goal-Oriented Approach. 15th World Congress on Computational Mechanics and 8th Asian Pacific Congress on Computational Mechanics, Yokohama, Japan, [July 31 - August 5, 2022]. 100 8. Main Achievements 2022 F. V. Caro, V. Darrigrand, J. Alvarez-Aramberri, E. Alberdi, and D. Pardo. 1D Painless Multi-Level Automatic Goal-Oriented h and p Adaptive Strategies using a Pseudo-Dual Operator. 22nd International Conference on Computational Science, London, United Kingdom, [June 21-23, 2022]. 2022 F. V. Caro, V. Darrigrand, J. Alvarez-Aramberri, E. Alberdi, and D. Pardo. Goal-Oriented hp-Adaptive Finite Element Methods: A Painless Multilevel Automatic Coarsening Strategy For Non-SPD Problems. 8th European Congress on Computational Methods in Applied Sciences and Engineering, Oslo, Norway, [June 5-9, 2022]. 2021 F. V. Caro, V. Darrigrand, E. Alberdi, and D. Pardo. A Painless GoalOriented hp-Adaptive Strategy for Indefinite Problems. 16th U.S. National Congress on Computational Mechanics, Chicago, U.S.A, [July 25-29, 2021]. 2021 F. V. Caro, V. Darrigrand, E. Alberdi, and D. Pardo. Goal-Oriented hpAdaptive Finite Element Methods: A Painless Multi-level Automatic Coarsening Strategy. 10th International Conference on Adaptive Modeling and Simulation, Gothenburg, Sweden, [June 21-23, 2021]. https://doi.org/10.23967/admos. 2021.044. 2021 F. V. Caro, V. Darrigrand, E. Alberdi, and D. Pardo. Painless Multi-level Automatic Goal-Oriented hp-Adaptive Coarsening Strategy. XVI Congreso de Matem´atica Aplicada, Gij´on, Spain, [June 14-18, 2021]. 8.3. Seminars 2022 F. V. Caro, V. Darrigrand, J. Alvarez-Aramberri, E. Alberdi, and D. Pardo. A Boundary Value Problem: A Painless Multi-Level hp-Adaptive Case. Centro Universitario de Ciencias Exactas e Ingenier´ıas, Universidad de Guadalajara, Guadalajara, M´exico, [March 9, 2022]. 2021 F. V. Caro, V. Darrigrand, E. Alberdi, and D. Pardo. A Painless Multilevel Automatic Goal-Oriented hp-Adaptive Coarsening Strategy. Red interuniversitaria de Ciencias-RIdeC, Lima, Per´u, [July 13, 2021]. 101 8. Main Achievements 8.4. Research Stays 2023 AGH University of Science and Technology, Krakow (Poland) Supervisor: Maciej Paszynski Date: 2 February 2023 - 31 March 2023 (58 days) 2021 CNRS-IRIT-ENSEEIHT, Toulouse (France) Supervisor: Vincent Darrigrand Date: 24 September 2021 - 25 November 2021 (61 days) 2020 CNRS-IRIT-ENSEEIHT, Toulouse (France) Supervisor: Vincent Darrigrand Date: 1 November 2020 - 4 December 2020 (34 days) 8.5. Implemented software In this dissertation, I used the FEM library from the MathMode group1. The group initially designed this library to address elliptic problems using an energybased-adaptive hp-strategy with H1-conforming discretizations. The library, written in Fortran90, supports solving problems in 1D, 2D (using quadrilateral elements), and 3D (using hexahedral elements). I contributed to the software in two significant ways. First, I expanded the energy-based-adaptive hp-strategy to a Goal-Oriented (GO) hp-adaptive algorithm that now handles both elliptic and non-elliptic problems. In this effort, I collaborated with Dr. Vincent Darrigrand and Dr. Julen Alvarez-Aramberri to introduce an upper bound of the error representation expressed through an inner product depending on the problem’s bilinear form. Furthermore, we collaborated to enhance the adaptive hp-strategy to fit the Multi-Adaptive GoalOriented (MAGO) framework for solving parametric Partial Differential Equations (PDEs). Our method seeks to produce reliable synthetic data or measurements, which experts can utilize for solving Inverse Problems (IPs) or training Neural Networks (NNs). We implemented this using piecewise-constant materials that align with the discretization. Instead of computing integrals for each sample, which would be time-consuming, we precomputed and saved the integrals for a unitary sample only once, optimizing the stiffness matrix computation process. 1https://www.mathmode.science/home 102 9. Conclusions and Future Work 9.1. Conclusions This dissertation mainly focuses on expanding an energy-based hp-adaptive algorithm previously limited to elliptic problems to both elliptic and non-elliptic problems under a Goal-Oriented (GO) framework. Chapter 4 proposes hand p-Goal-Oriented Adaptive (GOA) strategies suitable for both elliptic and potentially non-elliptic problems. These strategies use hierarchical basis functions to handle the hanging nodes, first performing a global and uniform refinement and then a coarsening step to remove certain basis functions. To determine which basis functions to remove, we employ an unconventional symmetric and positive definite bilinear form that quantifies the error in the Quantity of Interest (QoI). We test these algorithms on 1D Helmholtz and convectiondiffusion problems by applying the Laplace operator’s pseudo-dual problem. The numerical results show a linear convergence rate for the hscenarios and a quasiexponential rate for the pscenarios. Chapter 5 introduces an automatic adaptive mesh-generation strategy alternating between refinement and quasi-optimal hp-unrefinement procedures. Identifying which basis functions to remove efficiently presents a challenge. To address this, we extend a coarsening strategy previously tailored for energy-norm adaptivity to address non-elliptic problems and GOA strategies. Precisely, we consider the relevance of each basis function to the solution using an inner product associated with the problem’s bilinear form. Based on these evaluations, each coarsening step eliminates certain basis functions. The algorithm’s design simplifies implementation using hierarchical data structures that avoid the conventional 1-irregularity rule, which usually deals with hanging nodes. Our numerical results, which include 2D problems such as Poisson, Helmholtz, and convectiondominated equations, validate the algorithm’s robustness and fast convergence. The resulting algorithm is easy to implement, and due to its robustness and rapid convergence, it shows potential for easy adaptation to industrial scenarios. Chapter 6 highlights the strengths of our algorithm, showcasing its performance on a 3D heterogeneous Helmholtz equation-based problem. In Chapter 7, we introduce the Multi-Adaptive Goal-Oriented (MAGO) strategy to address the computational costs and dataset requirements associated with 103 9. Conclusions and Future Work accurately training a Deep Neural Network (DNN) to mimic the forward solver. Building upon our previously developed GOA approach for non-parametric Partial Differential Equations (PDEs), MAGO demonstrates promising results in efficiently generating a single hp-mesh. This optimized mesh ensures accurate computation of the QoI for multiple samples within a single GOA process. By combining the individual errors from all samples using l1,l2, and l∞norms, the MAGO approach provides sufficiently accurate solutions for all scenarios, including wave propagation examples with both h-refinements towards material discontinuities and strong p-refinements. The accuracy assessment of MAGO ’s adaptivity through box plots indicates that a more significant number of samples involved in the adaptive process (SA) leads to improved hp-grid results. Consequently, the statistical properties associated with the maximum relative error decrease as SA increases. The results underscore the robustness, speed, and computational efficiency of MAGO as an alternative for generating reliable databases while ensuring high accuracy. Furthermore, the computational costs in terms of Floating Point Operations (FLOPs) of the Single-Adaptive Goal-Oriented (SAGO) and MAGO strategies, based on factorization, are detailed in Tables 7.1 and 7.2, respectively. Notably, the MAGO approach demonstrates its effectiveness in problem-solving within this context, as evidenced by CMAGO < CSAGO. This observation underscores that the MAGO approach attains a higher level of accuracy while simultaneously reducing costs, rendering it a particularly suitable choice for addressing challenging problems. 9.2. Future Work In this dissertation, we identify potential paths for future research. One significant avenue is the extension of algorithms to address multi-physics problems, notably H(curl) and H(div). Enhanced parallelization and factorization techniques can reduce future computational resource requirements. Moreover, it is crucial to validate the efficacy of our algorithms in real-world scenarios such as Magnetotellurics, Controlled Sources, and Logging While Drilling. Furthermore, our approach to generating expansive databases, explicitly designed for DNN training, can be improved. By integrating our strategy with Machine Learning (ML) methodologies, we can expedite and improve the DNN training processes. An in-depth analysis of the impact of the nature and distribution of various random samples on Deep Learning (DL) inversion could provide critical insights for optimization. 104 BIBLIOGRAPHY [61] V. Darrigrand, A. Rodr´ıguez-Rozas, I. Muga, D. Pardo, A. Romkes, and S. Prudhomme. Goal-oriented adaptivity using unconventional error representations for the multidimensional Helmholtz equation. International Journal for Numerical Methods in Engineering, 113(1):22–42, 2018. (cited in page(s) 6, 7) [62] L. Demkowicz. Computing with hp-ADAPTIVE FINITE ELEMENTS. Vol. 1. One and Two Dimensional Elliptic and Maxwell Problems. Chapman and Hall/CRC., 2006. (cited in page(s) 3, 4, 5, 15, 20, 21) [63] L. Demkowicz, J. Gopalakrishnan, and A. H. Niemi. A class of discontinuous Petrov–Galerkin methods. Part III: Adaptivity. Applied Numerical Mathematics, 62(4):396–427, 2012. Third Chilean Workshop on Numerical Analysis of Partial Differential Equations (WONAPDE 2010). (cited in page(s) 4) [64] L. Demkowicz, J. Kurtz, D. Pardo, M. Paszy´nski, W. Rachowicz, and A. Zdunek. Computing with hp-ADAPTIVE FINITE ELEMENTS: Volume II Frontiers: Three Dimensional Elliptic and Maxwell Problems with Applications (1st ed.). Chapman and Hall/CRC., 2008. (cited in page(s) 3, 11, 15, 20, 21) [65] L. Demkowicz, J. Oden, W. Rachowicz, and O. Hardy. Toward a universal h-p adaptive finite element strategy, part 1. Constrained approximation and data structure. Computer Methods in Applied Mechanics and Engineering, 77(1):79–112, 1989. (cited in page(s) 20) [66] L. Demkowicz, W. Rachowicz, and P. Devloo. A Fully Automatic hpAdaptivity. Journal of Scientific Computing, 17:117–142, 2002. (cited in page(s) 3, 4, 8, 11, 21) [67] P. Deuflhard, P. Leinen, and H. Yserentant. Concepts of an adaptive hierarchical finite element code. IMPACT of Computing in Science and Engineering, 1(1):3–35, 1989. (cited in page(s) 3) [68] P. Di Stolfo, A. Schr¨oder, N. Zander, and S. Kollmannsberger. An easy treatment of hanging nodes in hp-finite elements. Finite Elements in Analysis and Design, 121:101–117, 2016. (cited in page(s) 4) [69] J. Donea and A. Huerta. Finite Element Methods for Flow Problems. John Wiley & Sons, 2003. (cited in page(s) 1) [70] W. D¨orfler. A convergent adaptive algorithm for poisson’s equation. SIAM Journal on Numerical Analysis, 33(3):1106–1124, 1996. (cited in page(s) 5) 111 BIBLIOGRAPHY [71] M. R. Dorr. The approximation theory for the p-version of the finite element method. SIAM Journal on Numerical Analysis, 21(6):1180–1207, 1984. (cited in page(s) 2) [72] T. Dunne. An Eulerian approach to fluid–structure interaction and goaloriented mesh adaptation. International Journal for Numerical Methods in Fluids, 51(9-10):1017–1039, 2006. (cited in page(s) 6) [73] A. D¨uster, E. Rank, and B. A. Szab´o. The p-Version of the Finite Element and Finite Cell Methods. Encyclopedia of Computational Mechanics Second Edition, pages 1–35, 2017. (cited in page(s) 3) [74] W. E, J. Han, and A. Jentzen. Deep Learning-Based Numerical Methods for High-Dimensional Parabolic Partial Differential Equations and Backward Stochastic Differential Equations. Communications in Mathematics and Statistics, 5(4):349–380, 2017. (cited in page(s) 57) [75] B. Erdmann, R. Roitzsch, and F. A. Bornemann. KASKADE: numerical experiments. Number TR 91-1 in Technical Report. Konrad-Zuse-Zentrum f¨ur Informationstechnik Berlin, 1991. (cited in page(s) 3) [76] M. Feischl, D. Praetorius, and K. G. van der Zee. An Abstract Analysis of Optimal Goal-Oriented Adaptivity. SIAM Journal on Numerical Analysis, 54(3):1423–1448, 2016. (cited in page(s) 6) [77] S. Franz and T. Linß. Superconvergence analysis of the Galerkin FEM for a singularly perturbed convection–diffusion problem with characteristic layers. Numerical Methods for Partial Differential Equations, 24(1):144– 164, 2008. (cited in page(s) 39) [78] S. Franz and G. Matthies. Local projection stabilisation on S-type meshes for convection–diffusion problems with characteristic layers. Computing, 87(3):135–167, 2010. (cited in page(s) 39) [79] S. Franz and H.-G. Roos. Error estimation in a balanced norm for a convection-diffusion problem with two different boundary layers. Calcolo, 51(3):423–440, 2014. (cited in page(s) 39) [80] L. E. Garcia-Castillo, D. Pardo, and L. F. Demkowicz. Energy-norm-based and goal-oriented automatic hp adaptivity for electromagnetics: Application to waveguide discontinuities. IEEE Transactions on Microwave Theory and Techniques, 56(12):3039–3049, 2008. (cited in page(s) 3) 112 BIBLIOGRAPHY [81] L. E. Garc´ıa-Castillo, D. Pardo, I. G´omez-Revuelto, and L. F. Demkowicz. A two-dimensional self-adaptive hp finite element method for the characterization of waveguide discontinuities. Part I: Energy-norm based automatic hp-adaptivity. Computer Methods in Applied Mechanics and Engineering, 196(49):4823–4852, 2007. (cited in page(s) 3) [82] S. Giani and E. Hall. An a-posteriori error estimate for hp-adaptive DG methods for elliptic eigenvalue problems on anisotropically refined meshes. Computing, 95(1):319–341, 05 2013. (cited in page(s) 4) [83] S. Giani, D. Sch¨otzau, and L. Zhu. An a-posteriori error estimate for hpadaptive DG methods for convection–diffusion problems on anisotropically refined meshes. Computers & Mathematics with Applications, 67(4):869– 887, 2014. High-order Finite Element Approximation for Partial Differential Equations. (cited in page(s) 4) [84] I. Gomez-Revuelto, L. E. Garc´ıa-Castillo, S. Llorente-Romano, and D. Pardo. A three-dimensional self-adaptive hp finite element method for the characterization of waveguide discontinuities. Computer Methods in Applied Mechanics and Engineering, 249-252:62–74, 2012. Higher Order Finite Element and Isogeometric Methods. (cited in page(s) 3) [85] I. Goodfellow, Y. Bengio, and A. Courville. Deep Learning. MIT Press, 2016. (cited in page(s) 57) [86] T. Gr¨atsch and K.-J. Bathe. Goal-oriented error estimation in the analysis of fluid flows with structural interactions. Computer Methods in Applied Mechanics and Engineering, 195(41):5673–5684, 2006. John H. Argyris Memorial Issue. Part II. (cited in page(s) 6) [87] B. Guo and I. Babuˇska. The h−pversion of the finite element method. Part 1: The basic approximation results. Computational Mechanics, 1(1):21–41, 1986. (cited in page(s) 2) [88] B. Guo and I. Babuˇska. The h−pversion of the finite element method. Part 2: General results and applications. Computational Mechanics, 1(3):203– 220, 1986. (cited in page(s) 2) [89] J. Han, A. Jentzen, and W. E. Solving high-dimensional partial differential equations using deep learning. Proceedings of the National Academy of Sciences, 115(34):8505–8510, 2018. (cited in page(s) 57) [90] A. Hashemian, D. Garcia, J. A. Rivera, and D. Pardo. Massive database generation for 2.5D borehole electromagnetic measurements using refined 113 BIBLIOGRAPHY isogeometric analysis. Computers & Geosciences, 155:104808, 2021. (cited in page(s) 8) [91] M. Holst and S. Pollock. Convergence of Goal-Oriented Adaptive Finite Element Methods for Nonsymmetric Problems. Numerical Methods for Partial Differential Equations, 32(2):479–509, 2016. (cited in page(s) 6) [92] M. Holst, S. Pollock, and Y. Zhu. Convergence of goal-oriented adaptive finite element methods for semilinear problems. Computing and Visualization in Science, 17(1):43–63, 2015. ID: Holst2015. (cited in page(s) 6) [93] P. Houston, C. Schwab, and E. S¨uli. Discontinuous hp-Finite Element Methods for Advection-Diffusion-Reaction Problems. SIAM Journal on Numerical Analysis, 39(6):2133–2163, 2002. (cited in page(s) 4) [94] T. Hughes. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Dover Publications, Inc., 2000. (cited in page(s) 1) [95] U. Iturrar´an-Viveros, A. M. Mu˜noz-Garc´ıa, O. Castillo-Reyes, and K. Shukla. Machine Learning as a Seismic Prior Velocity Model Building Method for Full-Waveform Inversion: A Case Study from Colombia. Pure and Applied Geophysics, 178(2):423–448, 2021. (cited in page(s) 57) [96] C. Jhurani and L. Demkowicz. Multiscale modeling using goal-oriented adaptivity and numerical homogenization. Part I: Mathematical formulation and numerical results. Computer Methods in Applied Mechanics and Engineering, 213-216:399–417, 2012. (cited in page(s) 6) [97] C. Jhurani and L. Demkowicz. Multiscale modeling using goal-oriented adaptivity and numerical homogenization. Part II: Algorithms for the Moore–Penrose pseudoinverse. Computer Methods in Applied Mechanics and Engineering, 213-216:418–426, 2012. (cited in page(s) 6) [98] C. Johnson. Numerical Solution of Partial Differential Equations by the Finite Element Method. Dover Publications, 2012. (cited in page(s) 1) [99] J. N. Jomo, N. Zander, M. Elhaddad, A. ¨ Ozcan, S. Kollmannsberger, R.-P. Mundani, and E. Rank. Parallelization of the multi-level hp-adaptive finite cell method. Computers & Mathematics with Applications, 74(1):126–142, 2017. 5th European Seminar on Computing ESCO 2016. (cited in page(s) 52) [100] M. Jungiewicz and A. Smywi´nski-Pohl. Towards textual data augmentation for neural networks: synonyms and maximum loss. Computer Science, 20(1):57–83, 2019. (cited in page(s) 57) 114 BIBLIOGRAPHY [101] G. E. Karniadakis and S. J. Sherwin. Spectral/hp Element Methods for Computational Fluid Dynamics. Numerical Mathematics and Scientific Computation. Oxford University Press, 2nd edition, 2013. (cited in page(s) 3) [102] E. Kharazmi and G. E. Karniadakis. Variational Physics-Informed Neural Networks for Solving Partial Differential Equations. arXiv preprint arXiv:1912.00873, 2019. (cited in page(s) 57) [103] Y. Khoo, J. Lu, and L. Ying. Solving parametric PDE problems with artificial neural networks. European Journal of Applied Mathematics, 32(3):421– 435, 2021. (cited in page(s) 57) [104] P. Kopp, V. Calo, E. Rank, and S. Kollmannsberger. Space-time hp-finite elements for heat evolution in laser powder bed fusion additive manufacturing. Engineering with Computers, 38:4879–4893, 2022. (cited in page(s) 4) [105] P. Kopp, E. Rank, V. M. Calo, and S. Kollmannsberger. Efficient multi-level hp-finite elements in arbitrary dimensions. Computer Methods in Applied Mechanics and Engineering, 401:115575, 2022. (cited in page(s) 4) [106] D. La Torre, H. Kunze, F. Mendivil, M. Ruiz Galan, and R. Zaki. Inverse Problems: Theory and Application to Science and Engineering 2015. Mathematical Problems in Engineering, 2015:796094, 2015. (cited in page(s) 57) [107] P. Ladev`eze. Upper error bounds on calculated outputs of interest for linear and nonlinear structural problems. Comptes Rendus M´ecanique, 334(7):399–407, 2006. (cited in page(s) 6) [108] P. Ladev`eze. Strict upper error bounds on computed outputs of interest in computational structural mechanics. Computational Mechanics, 42(2):271– 286, 2008. TY - JOUR. (cited in page(s) 6) [109] S. Larsson and V. Thom´ee. Partial Differential Equations with Numerical Methods, volume 45. Springer Berlin, Heidelberg, 2003. (cited in page(s) 1) [110] D. J. Lary, A. H. Alavi, A. H. Gandomi, and A. L. Walker. Machine learning in geosciences and remote sensing. Geoscience Frontiers, 7(1):3– 10, 2016. Special Issue: Progress of Machine Learning in Geosciences. (cited in page(s) 57) [111] Y. LeCun, Y. Bengio, and G. Hinton. Deep learning. Nature, 521(7553):436–444, 2015. (cited in page(s) 57) 115 BIBLIOGRAPHY [112] W. K. Liu, S. Li, and H. S. Park. Eighty Years of the Finite Element Method: Birth, Evolution, and Future. Archives of Computational Methods in Engineering, 29(6):4431–4453, 2022. (cited in page(s) 57) [113] L. Lu, X. Meng, Z. Mao, and G. E. Karniadakis. DeepXDE: A Deep Learning Library for Solving Differential Equations. SIAM Review, 63(1):208– 228, 2021. (cited in page(s) 57) [114] Y. Maday, A. T. Patera, and J. Peraire. A general formulation for a posteriori bounds for output functionals of partial differential equations; application to the eigenvalue problem. C. R. Acad. Sci. Paris S´er. I Math., 328(9):823–828, 1999. (cited in page(s) 5) [115] C. Mesztenyi, A. Miller, and W. Szymczak. FEARS: Details of Mathematical Formulation Univac 1100. Technical Report BN-994, University of Maryland, College Park, 1982. (cited in page(s) 3) [116] C. Mesztenyi and W. C. Rheinboldt. NFEARS: A Nonlinear Adaptive Finite Element Solver. Technical Report ICMA-87-113, Department of Mathematics and Statistics, University of Pittsburgh, 1987. (cited in page(s) 3) [117] T. Mitchell. Machine Learning. McGraw-Hill International Edit, McGraw Hill Higher Education, 1st edition, 1997. (cited in page(s) 57) [118] W. F. Mitchell and M. A. McClain. A Comparison of hp-adaptive Strategies for Elliptic Partial Differential Equations. ACM Trans. Math. Softw., 41(1), oct 2014. (cited in page(s) 4) [119] M. S. Mommer and R. Stevenson. A Goal-Oriented Adaptive Finite Element Method with Convergence Rates. SIAM Journal on Numerical Analysis, 47(2):861–886, 2009. (cited in page(s) 5) [120] K.-S. Moon, E. von Schwerin, A. Szepessy, and R. Tempone. Convergence Rates for an Adaptive Dual Weighted Residual Finite Element Algorithm. BIT Numerical Mathematics, 46(2):367–407, 2006. (cited in page(s) 5) [121] P. Morin, R. H. Nochetto, and K. G. Siebert. Convergence of adaptive finite element methods. SIAM Review, 44(4):631–658, 2002. (cited in page(s) 5) [122] C. D. Mote Jr. Global-local finite element. International Journal for Numerical Methods in Engineering, 3(4):565–574, 1971. (cited in page(s) 11) [123] J. Mu˜noz-Matute, E. Alberdi, D. Pardo, and V. M. Calo. Time-domain goal-oriented adaptivity using pseudo-dual error representations. Computer 116 BIBLIOGRAPHY Methods in Applied Mechanics and Engineering, 325:395–415, 2017. (cited in page(s) 7) [124] A. K. Noor. Global-local methodologies and their application to nonlinear analysis. Finite Elements in Analysis and Design, 2(4):333–346, 1986. (cited in page(s) 11) [125] J. T. Oden. Historical Comments on Finite Elements, page 152–166. Association for Computing Machinery, New York, NY, USA, 1990. (cited in page(s) 1) [126] J. T. Oden. Finite elements: An introduction. In P. G. Ciarlet and J. L. Lions, editors, Handbook of Numerical Analysis, volume 2, pages 3–15. North– Holland, Amsterdam, 1991. (cited in page(s) 2) [127] J. T. Oden and L. Demkowicz. Applied Functional Analysis. CRC press, 2017. (cited in page(s) 1) [128] J. T. Oden and A. Patra. A parallel adaptive strategy for hp finite element computations. Computer Methods in Applied Mechanics and Engineering, 121(1):449–470, 1995. (cited in page(s) 3) [129] J. T. Oden and S. Prudhomme. New approaches to error estimation and adaptivity for the stokes and oseen equations. International Journal for Numerical Methods in Fluids, 31(1):3–15, 1999. (cited in page(s) 5) [130] J. T. Oden and S. Prudhomme. Goal-oriented error estimation and adaptivity for the finite element method. Computers & Mathematics with Applications, 41(5):735–756, 2001. (cited in page(s) 5, 18, 20) [131] J. T. Oden, S. Prudhomme, and P. Bauman. On the extension of goal-oriented error estimation and hierarchical modeling to discrete lattice models. Computer Methods in Applied Mechanics and Engineering, 194(34):3668–3688, 2005. (cited in page(s) 6) [132] J. T. Oden and K. S. Vemaganti. Estimation of Local Modeling Error and Goal-Oriented Adaptive Modeling of Heterogeneous Materials: I. Error Estimates and Adaptive Algorithms. Journal of Computational Physics, 164(1):22–47, 2000. (cited in page(s) 6) [133] P. Oswald. On function spaces related to finite element approximation theory. Z. Anal. Anwend., 9(1):43–64, 1990. (cited in page(s) 2) 117 BIBLIOGRAPHY [134] J. S. Ovall. Asymptotically exact functional error estimators based on superconvergent gradient recovery. Numerische Mathematik, 102(3):543– 558, 2006. (cited in page(s) 5) [135] J. Panetier, P. Ladev`eze, and L. Chamoin. Strict and effective bounds in goal-oriented error estimation applied to fracture mechanics problems solved with XFEM. International Journal for Numerical Methods in Engineering, 81(6):671–700, 2010. (cited in page(s) 6) [136] G. Pang, L. Lu, and G. E. Karniadakis. fPINNs: Fractional PhysicsInformed Neural Networks. SIAM Journal on Scientific Computing, 41(4):A2603–A2626, 2019. (cited in page(s) 57) [137] M. Paraschivoiu and A. T. Patera. A hierarchical duality approach to bounds for the outputs of partial differential equations. Computer Methods in Applied Mechanics and Engineering, 158(3):389–407, 1998. (cited in page(s) 5) [138] M. Paraschivoiu, J. Peraire, and A. T. Patera. A posteriori finite element bounds for linear-functional outputs of elliptic partial differential equations. Computer Methods in Applied Mechanics and Engineering, 150(1):289–312, 1997. Symposium on Advances in Computational Mechanics. (cited in page(s) 5) [139] D. Pardo. Integration of HP-Adaptivity with a Two Grid Solver: Applications to Electromagnetics. PhD thesis, The University of Texas at Austin, 2004. (cited in page(s) 3) [140] D. Pardo. Multigoal-oriented adaptivity for hp-finite element methods. Procedia Computer Science, 1(1):1953–1961, 2010. ICCS 2010. (cited in page(s) 3) [141] D. Pardo, L. Demkowicz, C. Torres-Verd´ın, and C. Michler. PML Enhanced with a Self-Adaptive Goal-Oriented hp-Finite Element Method: Simulation of Through-Casing Borehole Resistivity Measurements. SIAM Journal on Scientific Computing, 30(6):2948–2964, 2008. (cited in page(s) 3) [142] D. Pardo, L. Demkowicz, C. Torres-Verd´ın, and M. Paszynski. TwoDimensional High-Accuracy Simulation of Resistivity Logging-WhileDrilling (LWD) Measurements Using a Self-Adaptive Goal-Oriented hp Finite Element Method. SIAM Journal on Applied Mathematics, 66(6):2085– 2106, 2006. (cited in page(s) 6) 118 BIBLIOGRAPHY [143] D. Pardo, L. Demkowicz, C. Torres-Verd´ın, and M. Paszynski. A selfadaptive goal-oriented hp-finite element method with electromagnetic applications. Part II: Electrodynamics. Computer Methods in Applied Mechanics and Engineering, 196(37):3585–3597, 2007. Special Issue Honoring the 80th Birthday of Professor Ivo Babuˇska. (cited in page(s) 6) [144] D. Pardo, L. Demkowicz, C. Torres-Verd´ın, and L. Tabarovsky. A goaloriented hp-adaptive finite element method with electromagnetic applications. Part I: electrostatics. International Journal for Numerical Methods in Engineering, 65(8):1269–1309, 2006. (cited in page(s) 6) [145] D. Pardo, L. E. Garc´ıa-Castillo, L. F. Demkowicz, and C. Torres-Verd´ın. A two-dimensional self-adaptive hp finite element method for the characterization of waveguide discontinuities. Part II: Goal-oriented hp-adaptivity. Computer Methods in Applied Mechanics and Engineering, 196(49):4811– 4822, 2007. (cited in page(s) 3) [146] D. Pardo, M. Paszynski, N. Collier, J. Alvarez, L. Dalcin, and V. M. Calo. A survey on direct solvers for Galerkin methods. SeMA Journal, 57(1):107– 134, 2012. (cited in page(s) 66) [147] M. Paszy´nski, L. Demkowicz, and D. Pardo. Verification of goal-oriented HP-adaptivity. Computers & Mathematics with Applications, 50(8):1395– 1404, 2005. (cited in page(s) 3) [148] M. Paszy´nski, L. Demkowicz, and D. Pardo. Verification of goal-oriented hp-adaptivity. Computers & Mathematics with Applications, 50(8):1395– 1404, 2005. (cited in page(s) 6) [149] M. Paszy´nski, D. Pardo, and V. Calo. Parallel simulations of 3d DC borehole resistivity measurements with goal-oriented self-adaptive hp finite element method. Journal of the Serbian Society for Computational Mechanics/Vol, 6(2):1–18, 2012. (cited in page(s) 3) [150] M. Paszy´nski and L. Demkowicz. Parallel, fully automatic hp-adaptive 3D finite element package. Engineering with Computers, 22(3):255–276, 12 2006. (cited in page(s) 52) [151] M. Paszy´nski, R. Grzeszczuk, D. Pardo, and L. Demkowicz. Deep Learning Driven Self-adaptive Hp Finite Element Method. In International Conference on Computational Science, pages 114–121. Springer, 2021. (cited in page(s) 57) 119 BIBLIOGRAPHY [152] J. Peraire and A. T. Patera. Bounds for linear-functional outputs of coercive partial differential equations: local indicators and adaptive refinement. In M. Papadrakakis, A. Samartin, L. Eca, and Y. Tsompanakis, editors, Advances in Adaptive Computational Methods in Mechanics, volume 47 of Studies in Applied Mechanics, pages 199–216. Elsevier Science BV, Amsterdam, 1998. (cited in page(s) 5) [153] J. Peraire and A. T. Patera. Asymptotic a posteriori finite element bounds for the outputs of noncoercive problems: the helmholtz and burgers equations. Computer Methods in Applied Mechanics and Engineering, 171(1):77– 86, 1999. (cited in page(s) 5) [154] S. Petrides and L. F. Demkowicz. An adaptive DPG method for high frequency time-harmonic wave propagation problems. Computers & Mathematics with Applications, 74(8):1999–2017, 2017. (cited in page(s) 4) [155] S. Pollock. Convergence of Goal-Oriented Adaptive Finite Element Methods. PhD thesis, University of California, San Diego, 2012. (cited in page(s) 6) [156] S. Prudhomme and J. T. Oden. On goal-oriented error estimation for elliptic problems: application to the control of pointwise errors. Computer Methods in Applied Mechanics and Engineering, 176(1):313–331, 1999. (cited in page(s) 5, 18, 20) [157] S. Prudhomme and J. T. Oden. Computable error estimators and adaptive techniques for fluid flow problems. In Error Estimation and Adaptive Discretization Methods in Computational Fluid Dynamics, volume 25 of Lecture Notes in Computational Science and Engineering, page 207. Springer, Berlin, 2003. (cited in page(s) 5) [158] V. Puzyrev. Deep learning electromagnetic inversion with convolutional neural networks. Geophysical Journal International, 218(2):817–832, 05 2019. (cited in page(s) 57) [159] V. Puzyrev and A. Swidinsky. Inversion of 1D frequencyand time-domain electromagnetic data with convolutional neural networks. Computers & Geosciences, 149:104681, 2021. (cited in page(s) 57) [160] W. Rachowicz and L. Demkowicz. An hp-adaptive finite element method for electromagnetics: Part 1: Data structure and constrained approximation. Computer Methods in Applied Mechanics and Engineering, 187(1):307–335, 2000. (cited in page(s) 11) 120