Some applications of reduced order modelling to slow-fast dynamical systems, turbulence models and elliptic PDEs

Programa de Doctorado "Matemáticas" PhD Dissertation Some applications of reduced order modelling to slow-fast dynamical systems, turbulence models and elliptic PDEs Author Alejandro Bandera Moreno Supervisors Tomás Chacón Rebollo Soledad Fernández García Macarena Gómez Mármol Sevilla, 2023 A mis abuelas. ACKNOWLEDGEMENTS En primer lugar, quiero agradecer su trabajo a mis tutores de tesis: Macarena, Soledad y Tomás. Voy a hacer un breve repaso a nuestra historia conjunta hasta ahora. Todo este viaje comenzó en el curso 2016/17, cuando le propuse a Tomás ser su alumno interno para estudiar sistemas caóticos. Él me aceptó y me presentó a Soledad, con la que estudié bifurcaciones en ecuaciones escalares autónomas. En el curso siguiente, continué como alumno interno de Soledad, a la que se le unió Macarena, para estudiar sistemas de ecuaciones diferenciales ordinarias con múltiples escalas de tiempo. El curso 2018/19 fue mi último año de grado, y defendí mi TFG sobre órbitas canard en sistemas lineales a trozos, tutorizado por Soledad y Macarena. Además, durante ese año disfruté de una Beca de Colaboración en el departamento. En el año del máster ya empecé a meter la cabeza en los modelos reducidos, defendiendo mi TFM sobre métodos de orden reducido para ecuaciones diferenciales. Después han venido cuatro años para la elaboración de esta tesis doctoral. Como veis, son bastantes años los que llevamos trabajando juntos y que espero que sean muchos más en el futuro. Debo hacer también una mención especial a mis profesoras de Matemáticas y Física durante el bachillerato, en particular a Natividad Cortés, mi tutora durante el último año preuniversiario. Gracias a ellas empezó este largo camino que hoy alcanza una nueva meta. Je voudrais également remercier les personnes avec qui j’ai travaillé lors de mes séjours académiques. Simon Brillet et Mejdi Azaïez pour mes séjours à Bordeaux, et Stéphane Menozzi et Alexandre Vidal pour mon séjour à Paris. J’ai hâte de continuer à travailler avec vous à l’avenir. A todos los que me han precedido como estudiantes de doctorado dentro del grupo de investigación M 2 S 2 M, si he podido ver lejos es porque me he subido a vuestros hombros. Espero que mi trabajo le sirva tanto a los que vengan detrás, como a mí me ha servido el vuestro. También a todos los compañeros del departamento, nos hacéis sentir como uno más desde el principio y nos hacéis la vida un poco más fácil. Por supuesto, quiero mencionar a mis compañeros del Colegio Mayor Hernando Colón durante los cinco años de carrera, con vosotros la vida universitaria alcanzó un significado mayor. Me llevo muy buenas experiencias y amistades que no olvidaré v vi nunca. No puedo terminar este párrafo sin agradecer además al equipo de dirección del Colegio Mayor, del que formé parte durante mis dos primeros años de tesis. Gracias a la oportunidad que me disteis, he mejorado en lo personal y en lo profesional. Espero que nuestra relación dure muchos años más. También mencionar a mis compañeros del IMUS durante estos cuatro años de tesis, a todos los que estáis y a todos los que os habéis ido. De no haber sido por vosotros, esto habría sido muy distinto. Han sido muchos los buenos momentos, que espero seguir compartiendo con vosotros durante mucho tiempo. Me gustaría además mencionar a mis compañeros del Doble Grado en Física y Matemáticas, me hicisteis ver que la física y las matemáticas en grupo son más divertidas. Faltan palabras para poder agradecer a mi nebulosa todo lo que ha hecho por mí. Sólo tú sabes lo feliz que me siento de que nos hayamos cruzado, y de que podamos compartir este momento juntos. Pronto vendrán nuevos momentos que celebrar, no me cabe duda. Por último, quiero agradecer todo su esfuerzo a mis padres, Pedro y Conchi, y a mi hermano, Víctor. A los primeros, les traje de cabeza unas semanas antes de Selectividad al decirles que quería estudiar el Doble Grado en Física y Matemáticas, y que este no estaba en Málaga. Gracias por vuestro apoyo incondicional a todas mis decisiones. Esta tesis tiene bastante de vosotros. RESUMEN En las últimas décadas, el progreso tecnológico ha elevado el rol de las simulaciones numéricas a ser una herramienta fundamental en la mayoría de ciencias y aplicaciones tecnológicas. Para ser más exactos, las simulaciones basadas en modelos apoyados en ecuaciones diferenciales adquieren un papel crítico debido a su diversidad de aplicaciones, que cubren desde ingeniería hasta economía, incluyendo áreas tan importantes como la medicina. Estas simulaciones sirven de plataformas virtuales que permiten la realización de experimentos, contribuyendo enormemente a la comprensión de diversas propiedades de los sistemas y sus dinámicas. Sin embargo, trabajar con los sistemas altamente complejos necesarios para aplicaciones realistas y efectivas supone un reto significante. Tales sistemas normalmente tienen asociadas decenas de miles o millones de grados de libertad, haciendo su resolución numérica prohibitivamente cara si se utilizan métodos convencionales. Esta complejidad requiere recursos computacionales extensos, necesitando horas o incluso días de computación, junto a Computación de Alto Rendimiento o arquitecturas informáticas específicas. Esto representa un gran problema, especialmente cuando se necesita una simulación en tiempo real o interactiva, o cuando se consideran múltiples valores de los parámetros, por ejemplo, en la asistencia en los procesos de toma de decisiones o diseño industrial. Esta tesis doctoral tiene como propósito abordar el reto de los altos costes computacionales asociados con los modelos realistas, y es en este punto donde el concepto de Modelado de Orden Reducido (ROM, por sus siglas en inglés) entra en juego. La filosofía subyacente en este campo recae en reemplazar el costoso problema original por uno alternativo, computacionalmente eficiente, que mantenga las propiedades cualitativas y cuantitativas esenciales de la solución original. En general, el uso de las técnicas ROM está motivado por la necesidad de mantener un balance entre la precisión y la eficiencia computacional en varias aplicaciones, permitiendo el análisis y la manipulación de sistemas complejos de una forma más práctica y manejable. Las principales contribuciones de esta tesis doctoral se enmarcan en tres aplicaciones diferentes: sistemas dinámicos lento-rápido, modelos de turbulencia y Ecuaciones en Derivadas Parciales (EDPs) elípticas. vii viii Relativo a la primera, el foco se sitúa en el estudio y simulación de modelos de redes de actividad neuronal que involucran múltiples escalas de tiempo. Estos modelos se tratan por medio de un método inédito basado en la Descomposición Ortogonal Propia, con el objetivo específico de abordar una de las limitaciones que presenta, a saber, la potencial pérdida de estructura del modelo original. En este caso, la pérdida de la separación de escalas temporales. En lo que respecta a los modelos de turbulencia, el énfasis se encuentra en el modelo de Smagorinsky y el método de Bases Reducidas. Aquí, el objetivo es salvar una limitación del método relativa al reto de obtener un estimador de error a posteriori basado en análisis matemático, que depende de la discretización numérica. En particular, desarrollamos un estimador de error basado en la teoría de cascada de energía de Kolmogórov. Por último, relativo a las EDPs elípticas, el objetivo que nos planteamos se centra en la resolución de EDPs elípticas simétricas y en calcular el mejor subespacio que aproxima su solución. La investigación gira en torno al método de la Descomposición Generalizada Propia (PGD), con el propósito de abordar una limitación asociada al cálculo de los modos PGD óptimos. Específicamente, el propósito es explorar la posibilidad de calcular directamente estos modos en una variedad de Grassmann, utilizando el conocido algoritmo del Gradiente Descendente, adaptado a este marco. ABSTRACT In the last decades, the technological progress has elevated the role of numerical simulations to be a fundamental tool in the majority of sciences and technical applications. More precisely, simulations based on models that are supported in differential equations acquire a critical role due to its diversity of applications, that cover from engineering to economy, and even areas so important as medicine. These simulations serve as virtual platforms enabling experimental realization, greatly contributing to the comprehension of diverse system properties and dynamics. Nevertheless, working with highly complex systems necessary for realistic and effective applications poses a significant challenge. Such systems usually have associated ten thousands or millions of degrees of freedom, rendering their numerical solving prohibitively expensive using conventional methods. This complexity demands extensive computational resources, sometimes requiring hours or even days of computation alongside High-Performance Computing or specialized informatics architectures. This fact represents a major problem, especially when a real-time or interactive simulation is needed, or when multiple values of the parameters need consideration, such as assisting decision-making processes or industrial design. This PhD dissertation aim to address the challenge of the high computational costs associated with those realistic models, and it is at this point where the concept of Reduced Order Modelling (ROM) comes into play. The underlying philosophy in this field lies in replacing the original costly problem with an alternative, computationally efficient one, that retains essential qualitative and quantitative properties of the original solution. Overall, the use of ROM techniques is motivated by the need to strike a balance between accuracy and computational efficiency in various applications, enabling the analysis and manipulation of complex systems in a more practical and manageable manner. The main contributions of this PhD dissertation are encompassed in three different applications: slow-fast dynamical systems, turbulence models and elliptic Partial Differential Equations (PDEs). Concerning the first, the focus lies on studying and simulating network models of neuronal activity involving multiple timescales. These models are going to be treated ix 2 introduction (Carrizosa et al., 2017) and control problems (Ito & Ravindran, 1998). Overall, the use of ROM is motivated by the need to strike a balance between accuracy and computational efficiency in various applications, enabling the analysis and manipulation of complex systems in a more practical and manageable manner. There exists several approaches in order to obtain the ROM. They can be broadly classified into different classes (Lassila et al., 2014; Silva et al., 2007): Proper Orthogonal Decomposition (POD) methods (Kerschen et al., 2005), Reduced Basis (RB) methods (Boyaval et al., 2010; Quarteroni et al., 2016), Balancing methods (Gugercin & Antoulas, 2004), simplified physics (Majdabadi et al., 2015), nonlinear manifold methods (Mojgani & Balajewicz, 2021) and projection based reduction, like Proper Generalized Decomposition (PGD) (Chinesta et al., 2011) and Krylov subspace methods (Bai, 2002). In this thesis, we will focus on the study and applications of three of the former approaches, the POD method, the RB method and the PGD method. For the reduction of the nonlinear terms, we highlight the EIM presented in Barrault et al. (2004), and then extended in Grepl et al. (2007). These problems have been selected within the scope of the Spanish national projects "ROM Optimization for Architecture and Design. ROAD" (RTI 2018 - 093521 -B-C 31 ) and "Modelos de orden reducido híbridos aplicados a flujos incompresibles y redes neuronales cerebrales. HYROM" (PID 2021 - 123153 OB-C 21 ), and it also aligns with the development of the European project "Accurate Roms for Insdustrial Applications. ARIA" (H2020-872442). The POD is a numerical method that enables a reduction in the complexity of computer intensive simulations. Its main use is to decompose a physical field, for instance, velocity or pressure, into a set of functions, the so-called POD modes, modulated by some coefficients. This decomposition is performed depending on the different variables that influence its physical behaviors (Kerschen et al., 2005). It is a generalization of the Principal Component Analysis (Pearson, 1901) in the field of statistics, or the Singular Value Decomposition (SVD) in the field of linear algebra, as it refers to singular values and singular vectors of the physical field. In fact, the optimality of the POD is ensured by the celebrated Eckart and Young theorem (Eckart & Young, 1936). Once the dominant structures are extracted, usually a Galerkin projection of the governing differential equations can be employed to obtain a low-order system (Berkooz et al., 1993). The resulting low-order model is usually addressed as standard POD-Galerkin ROM. Addressing RB methods, they were first introduced in the late seventies in Almroth et al. (1978) and Nagy (1979) for structural analysis of beam and arches. In the early introduction 3 eighties, the mathematical basis developed in previous works was extended, which resulted in the presentation of several works on RB methods for structural problems (Almroth et al., 1981; Noor & Peters, 1980). Also in the eighties, the presentation of the first analysis for RB methods was made in Fink and Rheinboldt (1983a,1983b). In those works, the effectiveness of the RB method is provided by a theoretical error estimation for a general non-linear problem. We will be focusing on the application of RB methods to fluid dynamics problems. The first works in this field were presented in the late eighties in Gunzburger (1989) and Peterson (1989). However, the RB method made significant progress at the beginning of the new millennium with the development of a general mathematical background analysis. One of the key advancements was the creation of an a posteriori error bound estimator and an algorithm for the basis functions, which marked a milestone in the expansion of the RB method. The a posteriori error bound estimator is particularly useful in the information selection process for constructing the reduced basis space. This selection is typically performed using a Greedy algorithm. In each iteration of the Greedy algorithm, the snapshot associated with the parameter value that exhibits the highest error between the High Fidelity (HF) and RB solutions is chosen. Determining the actual error precisely can be computationally expensive. Therefore, the weak Greedy algorithm is often employed in RB methods. This algorithm replaces the exact error calculation with an error estimator, which offers a more practicable alternative. The first appearance of the weak Greedy algorithm in the context of RB methods was in Prud’Homme et al. (2002) and in Veroy et al. (2003). Furthermore, these works exploit the parameter dependence of the problem, presenting an offline/online decomposition of it. That is, the problem is divided in two phases, an offline phase, devoted to the computation of the parameter-independent matrices, and an online phase, focused on computing the RB solution depending on the parameter value. Also, in Maday et al. (2002), general a priori convergence properties are studied. The need for an a posteriori error estimator in Reduced Basis (RB) methods is evident in their development. Typically, this estimator relies on the dual norm of the residual, providing a means to certify the method, that is, to yield an error bound for the RB approximation. In the context of linear Partial Differential Equations (PDEs), the process of obtaining an a posteriori error bound estimator is well-established. It involves applying continuity and coercivity results to the actual error equation (Quarteroni et al., 2016). However, when dealing with non-linear PDEs, the development of an a posteriori error bound estimator through mathematical analysis depends strongly 4 introduction on the selected discretization. For instance, the Brezzi-Rappaz-Raviart (BRR) theory (Brezzi et al., 1980) is employed for approximating regular branches of solutions in non-linear PDEs, enabling effective error estimation in these cases. Nevertheless, like any theoretical framework, it has its limitations and challenges. For problems with high-dimensional parameter spaces, the (weak) Greedy Algorithm can be coupled with the POD procedure. This coupled approach, commonly referred as POD+Greedy approach, was introduced in Haasdonk and Ohlberger (2008). Usually, in the POD+Greedy procedure, each greedy step is invoking a temporal compression step by performing a POD. In Haasdonk (2013), the existing convergence results for the Greedy algorithm were extended to the POD-Greedy algorithm case. Unfortunately, the POD methods present a great disadvantage, as it needs knowledge about the original solution before any ROM approximation can be applied. This knowledge has to be large enough to capture the main behaviors of the actual solution over a range of the parameters. This can be extremely expensive in terms of computational time or storage. The PGD methods address this and other issues arising in the POD framework. A PGD method is an iterative numerical technique developed for solving PDEs under specified boundary conditions. Their first appearance was in the mid-2000 in Ammar et al. (2006). This method was designed to alleviate the ’curse of dimensionality’ when solving problems defined in high-dimensional spaces. The main idea behind this technique is to assume that the solution admits a decomposition into a sum of basis functions which do not need to be orthogonal between them, that is, the solution admits a representation as a sum of rank- 1 tensors. Then, the PGD algorithm aims to approximate the solution through successive enrichment by means of a deflation algorithm. In each iteration, a new mode is calculated and incorporated into the approximation. Therefore, the complexity of the reduced problem increases linearly. This is clearly a great improvement and allows us to work efficiently with problems defined in much higher dimensions. Furthermore, thanks to linearity and separability of the representations of the method, only one-dimensional operations are needed. In addition, this approximation is constructed without any a priori knowledge of the solution. There are several ways to approximate the reduced solution. A number of different PGD methods are described in Nouy (2010). In this thesis, we will focus on the PGD method proposed in Azaïez et al. (2018), the so-called intrinsic PGD. This technique is based on the construction of recursive approximations on finite-dimensional optimal subspaces. These approximations are obtained by means of the minimization of the mean parametric error of the residual, introduction 5 similar to the one introduced in Falcó and Nouy (2011). In the latter work Azaïez et al. (2020), this technique is applied along with a deflation algorithm and a Power Iterate (PI) method to compute the PGD modes for parametric elliptic PDEs both, symmetric and non-symmetric. We begin this thesis by dedicating a chapter to introducing the fundamental concepts of model order reduction. Therefore, Chapter 1serves as a foundation for the subsequent sections of the thesis where we delve deeper into the topic. In this chapter, we start in Section 1.1from the concept of differential equations depending on one or various parameters. On the one hand, we introduce a general system of parametric ODEs, mainly used in Part I. On the other hand, we present a general parametrized PDE, that will be the main object of study in Parts II and III. Then, we present in Section 1.2the algebraic equivalent of the previous differential equations, that is, a representation of that problem through operations between matrices and vectors. The interested reader is referred to Quarteroni et al. (2016) for some examples. Section 1.3 is devoted to the argumentation of the need of the development of ROM techniques. Additionally, it provides an overview of the most commonly employed techniques, which will be detailed in the subsequent chapters, namely, the POD method and the Greedy Algorithm, along with a coupled strategy that combines the POD and Greedy procedures. In order to treat possible non-linearities that could be present in the original models, we introduce in Section 1.4a technique that helps us to approximate any non-linear term, the Empirical Interpolation Method (EIM). Furthermore, we dedicate a chapter to the presentation of two models of interest for the rest of the thesis. Hence, Chapter 2acts like a presentation for these two models of interest. On the one hand, we present in Section 2.1a model of neuronal Intracellular Calcium Concentration, presented in Krupa et al. (2013) and studied in Fernández- García and Vidal (2020), among others. This model is used in Part I. On the other hand, in Section 2.2, we present the basic Large Eddy Simulation (LES) turbulence model, the celebrated Smagorinsky turbulence model introduced in Smagorinsky (1963), along with its associated reduced basis problem. This model is used in Part II. The rest of the thesis is then divided into three main parts, each focusing in one particular kind of problem modelled by a system of differential equations: slow-fast dynamical problems, that will be tacked by a Proper Orthogonal Decomposition (POD) approach, turbulence problems, that will be addressed by a Reduced Basis (RB) approach, and symmetric elliptic PDEs, that will be studied under a Proper Generalized Decomposition (PGD) approach. Each part of the thesis explores and 6 introduction analyzes the respective ROM method in detail. It outlines a primary difficulty that we aim to address and introduces the proposed enhancement. These enhancements are designed in a manner, so they serve to improve the study of their corresponding problem. The content of each part is outlined and explained in the following to ensure clarity and coherence throughout the thesis. Every part has an introduction, where we present the context and the main problem. part i In this part, our emphasis is on the study and simulation of network models of neuronal activity that involve different timescales, modelled by a system of ODEs. In particular, we focus on the POD method, with a specific goal of addressing one of the limitations that it presents, namely, the potential loss of structure from the original model. This part aligns with the Spanish national project HYROM. In the last thirty years, extensive research has been groundbreaking in modeling neuronal activity networks. These models use dynamical systems coupled as networks to represent interactions among neuron assemblies, employing various paradigms and scales such as neural mass (Coombes, 2023), neural fields (Faugeras et al., 2008), mean-fields (Byrne et al., 2022), dendritic trees (Traub et al., 2001), and others. However, persistent challenges include difficulties in classifying synchronization features between network cells and the computational cost of reliable simulations for realistic-sized networks. The challenge arises from varying timescales in neuronal activity, reflected in models using multiple timescales systems as cell dynamics within a network. Synchronization mechanisms in these systems differ from harmonic oscillators, particularly in terms of the impact of coupling types and strengths. Despite numerous studies exploring synchronization in the presence of complex phenomena characteristic of multiple timescale dynamics, such as, canard phenomena and bursting oscillations (Coombes & Wedgwood, 2023; Ermentrout & Wechselberger, 2009; Ersöz et al., 2017; Izhikevich, 2001; Semenov & Fradkov, 2021), a comprehensive analysis of synchronization types and associated conditions is needed, especially for large networks, along with their physio-pathological implications. However, simulating large-scale network models can be computationally expensive, and while intrusive ROM approximate accurately the solution of the original system (Bandera Moreno et al., 2022), they may lead to the loss of the problem’s original introduction 7 structure, such as the slow-fast separation of variables common in neuroscience problems. The original problem involves quasi-sparse matrices and vectors, but reductions result in full matrices and vectors. To retain the problem’s structure, various strategies have been proposed, such as a Krylov subspace-based technique for electronic circuits (Freund, 2011), and, more recently, a domain decomposition approach for fluid simulations (Anderson et al., 2022). In this thesis, we present an alternative approach using the POD method and a data-driven procedure based on collinearity analysis. This approach automatically determines the best structure for the reduced model while significantly reducing computational time. Part Iis outlined as follows: • In Chapter 3, we present the theoretical side of this part. Our emphasis is on a system of parametric ODEs, yet it’s worth noting that this technique is adaptable to the realm of PDEs characterized by any qualitative separation in variables. In Section 3.1, we introduce an enhanced version of the original POD technique, designed to maintain the structural integrity of the initial system, that we name the Proper Orthogonal Block Decomposition (POBD). We also present some mathematical results concerning this novel technique in Section 3.1.1. More precisely, we adapt Eckart and Young Theorem to our case, obtaining an error bound for our method, that we present in Proposition 3.1. Furthermore, in Section 3.2, we introduce a preliminary data-driven phase aimed at determining the optimal structure for the reduced model by means of a correlation analysis between variables and group of variables. This step is conducted to retain the maximum possible structural fidelity of the original model while simultaneously minimizing the number of reduced variables. The whole process of determining the structure for the reduced model and performing the POBD is named as Automatic Proper Orthogonal Block Decomposition (APOBD). • Chapter 4is devoted to the validation and numerical study of the methods developed in Chapter 3, namely the POBD and APOBD methods. More precisely, we will apply these reduction techniques to problems arising from the mathematical modelling of biological phenomena exhibited in different cells, like neurons. Section 4.1will be devoted to the study of a model of Intracellular Calcium Concentration, presented in Krupa et al. (2013), studied in Fernández-García and Vidal (2020) and in Bandera Moreno et al. (2022), among others. In Section 4.2.1, we will study the celebrated Hindmarsh-Rose model, introduced in Hindmarsh 8 introduction and Rose (1984) and widely studied (Barrio & Shilnikov, 2011; Bonaventura et al., 2022; Etémé et al., 2017), among others. Finally, Section 4.2.2will be focused on the study of pancreatic β -cells models described by a Fold/Hopf Bursting model (Izhikevich, 2000; Pernarowski et al., 1992; Smolen et al., 1993). Our focus involves network models characterized by a substantial number of interconnected oscillators. To evaluate the efficacy of reduction techniques, we will employ conventional metrics like speed-up, that is, the ratio between the solving computational times of the original and reduced model, and the number of equations in reduced models, in comparison to the original model. Additionally, we will use the sparsity index, an indicator reflecting the reproduction of the original problem’s structure, gauged through the sparsity of the matrices associated with each reduced model. part ii In this part, our emphasis is the study and simulation of turbulent flows, modelled by PDEs. Our focus would be on the RB method, with a specific goal of addressing one of the limitations that it presents, namely, the difficulty of obtaining mathematics-based a posteriori error estimator, and its dependence on the numerical discretization. In particular, we study the development of an a posteriori error estimator based upon the Kolmogórov’s energy cascade theory (Kolmogórov, 1941). This part aligns with the Spanish national projects ROAD and HYROM, and with the European project ARIA. Kolmogorov’s energy cascade theory, also known as K 41 theory, is a foundational concept in the modelling of turbulence (Chacón Rebollo & Lewandowski, 2014; Kolmogórov, 1941; Richardson, 1922). It elucidates the transfer of kinetic energy from large-scale eddies to progressively smaller scales in turbulent flows. This energy cascade, marked by self-similar scaling behavior, culminates in the dissipation of energy into heat at the smallest scales. The theory provides a statistical description of turbulence, serving as a guiding framework for developing turbulence models and conducting numerical simulations of turbulent flows. The energy cascade concept has led to using high-order modes of the POD expansion for modeling small scales in the inertial range. However, achieving precision with this method requires a substantial number of modes, resulting in ROM computation times comparable to a direct Full Order Model (FOM) solver (Hijazi et al., 2020). An alternative approach involves constructing ROMs for established turbulence models, particularly addressing Large Eddy introduction 9 Simulation (LES) models (Germano et al., 1991; Smagorinsky, 1963). In this strategy, the focus is not on modeling turbulence through ROM, but on creating surrogate LES models for significantly faster computations. From the family of LES models, we highlight the Smagorinsky model (Smagorinsky, 1963), in which the eddy viscosity is modeled by a non-linear term, that depends on the mesh size and the modulus of the velocity gradient. Certified RB LES Smagorinsky models for steady flows have demonstrated substantial computing speed-ups for both single and multi-parametric configurations (Ballarin et al., 2020; Caravaca García, 2022; Chacón Rebollo et al., 2023; Chacón Rebollo et al., 2017). Error certification relies on a posteriori error estimation procedures, formulated using the BRR theory for approximating regular branches of non-linear PDEs (Brezzi et al., 1980), leveraging the enhanced regularity of Smagorinsky model solutions. However, this approach has significant limitations. Constructing a posteriori error indicators is closely tied to the specific turbulence model and FOM in use, as well as to the nature of the parameters, whether physical or geometrical, considered in the ROM. For the Smagorinsky turbulence model, there exists a well-established theory (potentially extendable to more general LES models) based on the enhanced regularity compared to the Navier-Stokes equations, facilitating the development of a posteriori error estimation (Caravaca García, 2022). Nevertheless, these estimations are computationally intensive, involving the calculation of stability factors in space-time, that additionally depend on discretization parameters. This part of the thesis focuses on establishing a RB method applicable to a wide range of transient LES parametric turbulence models, irrespective of the FOM employed. The key approach involves an error indicator rooted in the physics of turbulent flows rather than the mathematical analysis of models. This indicator quantifies the deviation between the energy spectrum of a trial solution and the theoretically predicted spectrum by Kolmogórov’s energy cascade theory of turbulence in statistical equilibrium. Part II is outlined as follows: • In Chapter 5, we present the theoretical side of this part. The main purpose behind this chapter is the introduction of an a posteriori error indicator for turbulent flows based only on physical properties of turbulent flows, instead of on a mathematical analysis of the selected model and discretization. More precisely, for the development of this indicator, we have used general results from the statistical turbulence theory, attributed to Kolmogórov (Kolmogórov, 1941). 10 introduction In Section 5.1, we introduce the error indicator based upon Kolmogórov’s energy cascade theory, along with some basics results from this theory. Furthermore, in Section 5.1.2, we develop some error estimates for Kolmogórov’s indicator. Theorem 5.2states that the error between two velocity fields of locally isotropic flows, that is, the statistical properties of the mean flow are invariant under rotations, solutions of the Smagorinsky turbulence model, is fully driven by the indicator. Also, in Theorem 5.3, we state that the error between the FOM and ROM solutions is driven by the estimator plus the projection error on the reduced space. To end this chapter, we describe in Section 5.2a general algorithm based on the a posteriori error indicator for the construction of the reduced basis model of LES turbulent model. Furthermore, as the indicator is not expected to vanish as the reduced dimension grows, but to converge to the indicator for the solution of the FOM, then, the usual criterion for the selection of the best parameter in the Greedy algorithm is no longer useful. Instead, we propose in Algorithm 1new criteria for the selection of the next parameter at each step. Algorithm 1first aims to add to the basis the information of the solutions whose energy spectrum is farther away from the theoretical one. Additionally, if the selected parameter has already been selected in a previous iteration, we are assuming that the reduced error indicator is getting close to the indicator of the FOM solution, and the previous solution plays the role of the full-order one to compare with it. • Chapter 6is devoted to the validation and numerical study of the methods developed in Chapter 5. In particular, we address some academic tests for the Smagorinsky turbulence model. We consider 2 D periodic flows, yet meaningful, as they are designed to present the inertial spectrum predicted by Kolmogórov’s theory. We assume that the Smagorinsky flow already has a part of its spectrum within the inertial range. First, in Section 6.1, we conduct a preliminary study. Our attention in this section is specifically directed towards the handling of the eddy viscosity term, the selection of the POD+Greedy procedure and the study of the time variability of the error indicator. Then, in Section 6.2, we present the implemented test in order to validate the performance of the use of the error indicator. We compare the errors obtained with this strategy versus those obtained with the best error indicator, that is, introduction 11 the exact error. As a result, we obtain errors between the ROM and the FOM solutions quite close to the optimal ones. Furthermore, we obtain a spectral decay of the error as the dimension of the reduced basis space increases, quite close to the one obtained if the exact error is used as indicator. We obtain speed-up rates of computing times around 23 , which is quite satisfying for an evolution turbulence model. part iii In this part, our emphasis is on the resolution of symmetric elliptic PDEs, and the computation of the best subspace that approximates their solution. To that end, we focus on the PGD method, with the specific objective of addressing one of the limitations that it presents, namely, the computation of the optimal PGD modes. In particular, we aim to study the possible calculation of these modes directly in a Grassmann manifold by means of the Gradient Descent (GD) algorithm. This part aligns with the Spanish national projects ROAD and HYROM. Among all the PGD methods, we wish to bring attention to the formulation of the PGD introduced in Azaïez et al. (2018), termed the intrinsic PGD. This approach is founded on constructing recursive approximations within finite-dimensional optimal subspaces. These approximations are derived through minimizing the mean parametric error of the residual, akin to the approach introduced in Falcó and Nouy (2011). In that paper, the authors establish the existence of an optimal subspace, for any given dimension, that achieves the most accurate approximation in terms of the error between the precise solution and the Galerkin solution computed within the subspace, measured by the mean parametric norm associated with the elliptic operator. This mirrors the best approximation characteristic of POD subspaces, with the distinction that in their scenario the norm is parameter-dependent. Then, they employ a deflation technique to construct a sequence of approximate solutions within finite-dimensional optimal subspaces, directly during the online step. Furthermore, they demonstrate that the partial sums of these approximations converge to the continuous solution in the mean parametric elliptic norm. Additionally, they establish a strong connection between the conventional PGD for the specified parametric problem and the deflation algorithm presented in this paper. In fact, both methods coincide if the minimum of the functional is computed exactly. This linkage suggests the feasibility of computing 18 basics of reduced order models 1.1.2Parametric Partial Differential Equations In this section, let us denote by Ω⊂Rd , with d=2,3 , a spatial domain, Γ will denote the boundary of Ω, and QT=If×Ωstands for the spatio-temporal domain. Then, for every parameter vector γ∈ D, we can consider the following PDE ∂tu(γ) + L(γ)u(γ) = f(γ), in QT, (3) where f(γ) = f(t,x;γ) is a given function, and L(γ) = L(x;γ) is a generic elliptic operator acting on the unknown u(γ) = u(t,x;γ). We could also consider a PDE with no time dependence, the so-called elliptic PDEs. In that case, for every parameter vector γ∈ D, we would have the following L(γ)u(γ) = f(γ), in QT, (4) where f(γ) = f(x;γ) is a given function, and L(γ) = L(x;γ) is a generic elliptic operator acting on the unknown u(γ) = u(x;γ) . In the following, we consider that we are working with non-elliptic PDEs, as they are more general. Nevertheless, the similar results can be obtained for the elliptic case. Furthermore, we can consider an initial condition u(0, x;γ) = u0(x,γ),x∈Ω,γ∈ D, where u0is a given function, and we also consider boundary conditions u(t,x,γ) = uD(x,γ),t∈If,x∈ΓD,γ∈ D, ∂u ∂n(t,x,γ) = uN(x,γ),t∈If,x∈ΓN,γ∈ D, with uDand uNgiven functions, and ΓD∪ΓN=Γ, with no-intersecting interiors. Our main objective is to solve the previous PDE problems (4) and (3) numerically, so we introduce its weak formulation. Let us consider X=X(Ω) a suitable Hilbert space endowed with the X -norm, and homogeneous Dirichlet boundary conditions. Then, for every t∈Ifand γ∈ D, we consider the problem (Find u(t;γ)∈Xsuch that (∂tu(t;γ),v)Ω+A(u(t;γ),v;γ) = F(t,v;γ),∀v∈X,(5) 1.1 differential equations depending on parameters 19 where A(·,·;γ):X×X→R , is a bilinear operator associated to the elliptic operator L(γ) and the Neumann boundary condition uN , and F(t,·;γ) is a linear operator associated to function f(γ)for every γ∈ D. The existence and uniqueness of weak solutions of the problems studied in parts II and III, that are like problem (5) or like the variational formulation of problem (4) , is given by the celebrated Lax-Milgram theorem and ore one of its generalizations. Theorem 1.1(Lax-Milgram) Let X be a Hilbert space and A(·,·;γ) a bilinear form in X . We assume that, for every γ∈ D, the bilinear form is •bounded, that is, |A(·,·;γ)| ≤ C∥u∥X∥v∥X, with C >0; and •coercive, that is, A(·,·;γ)≥c∥u∥2 Xwith c >0. Then, for any f∈X′ , with X′ denoting the dual space of X , and for every γ∈ D , there exists a unique solution u ∈X to the equation A(u,v,γ) = f(v),∀v∈X, and the following bound holds ∥u∥X≤1 c∥f∥X′. Space-time discretization As our aim is to solve problem (5) numerically, a discretization is needed. Through this dissertation, we will work with the finite element approach, and we will consider a Galerkin approximation of problem (5). Let us define {Th}h>0 a uniformly regular family of triangulations of Ω , where the lower index h>0 is related to the mesh size. Then, we consider a discrete subspace Xh=Xh(Ω)∈X , of finite dimension Nh , as an inner approximation of a given Hilbert space X. Then, if we set uh(0; γ) as an approximation of u0(γ) in the space Xh , for a given t∈If,γ∈ D, the semi-discretized equation in space is given as (Find uh(t;γ)∈Xhsuch that (∂tuh(t;γ),vh)Ω+A(uh(t;γ),vh;γ) = F(t,vh;γ),∀vh∈Xh. 20 basics of reduced order models For the sake of clarity, we present now an explicit Euler scheme for time discretization. Extension to more general numerical schemes is possible, but the following analysis would be more tedious. Let Nt be a positive integer that denotes the number of time steps that will be considered, ∆t=Tf/Nt the time step and tk=k∆t , for k=1, . . . , Nt . We also denote by uk h(γ)the approximation of uh(tk;γ). Then, the complete space-time discretization of problem (5) can be expressed as          For any k=1, . . . , Nt, and γ∈ D, assuming known uk−1 h(γ)∈Xh, find uk h(t;γ)∈Xhsuch that M(uk h(γ),vh;γ) + A(uk h(t;γ),vh;γ) = Fk(vh;γ) + M(uk−1 h(γ),vh;γ),∀vh∈Xh, (6) where M(·,·;γ):X×X→R is defined by M(u,v;γ) = (u,v)Ω/∆t , and Fk(vh;γ) is a suitable approximation of the second member for any γ∈ D. Time-dependent parametric PDEs will be the main object of study in Part II, devoted to the development of reduced order models of incompressible turbulent flows. While, elliptic parametric PDEs will be the main object of study in Part III, focused on the study of the intrinsic Proper Generalized Decomposition expansion in a Grassmann manifold. 1.2 algebraic representation In order to treat numerically both, ODEs and PDEs, we will be working with the algebraic equivalent of models (1) and (6) , respectively. The algebraic equivalent of a discrete problem, in general terms, is an algebraic equation, that is, a representation of the problem by means of operations between matrices and vectors, see Quarteroni et al. (2016) for some examples. It will also be useful when constructing the reduced models thanks to the techniques that we will present in Section 1.3. Before we begin to develop the algebraic equivalents of models (1) and (6) , we need to introduce some notation that will be useful in the latter. Notation 1.1 With respect to vectors notation: •0N, the null column vector with N components. •1N, the column vector with N components, and all its components equal to 1. 1.2 algebraic representation 21 Notation 1.2 On the other hand, some useful matrix notation will be: •INis the identity matrix of order N. •ON is the null matrix of order N . If it has two subindexes, they will denote the number of rows and the number of columns, respectively. • If d is a column vector, D=diag(d) is a diagonal matrix whose entries are the components of vector d. Bearing these notations in mind, we can focus now on the presentation of the algebraic versions for both, models (1) and (6). 1.2.1Algebraic representation of systems of parametric ordinary differential equations For the construction of the algebraic representation of model (1) , we will consider, like in model (6) , an explicit Euler scheme for time discretization. Again, extensions to more general numeric schemes is possible, but we do not cover it here for the sake of brevity. Let Nt be a positive integer that denotes the number of time steps that will be considered, ∆t=Tf/Nt the time step and tk=k∆t , for k=1, . . . , Nt . We also denote by xk(γ)the approximation of x(tk;γ). The algebraic representation of model (1) , for every k=1, . . . , Nt , will be the following: ∆xk(γ) = f(xk−1,γ) = A(γ)xk−1(γ) + b(xk−1(γ);γ), (7) where ∆xk(γ) denotes the variation of the solution at time step k , A(γ) is the algebraic representative of the linear part of model (1) , and b(·;γ) summarizes the affine part and the possible nonlinearities involved in model (1) . Then, in order to obtain the next value xk(γ), we just need to sum x(γ)k=xk−1(γ) + ∆t∆xk(γ). 1.2.2Algebraic representation of systems of parametric partial differential equations For the construction of the algebraic representation of model (6) , we consider that we are working in a finite element framework. Then, considering the basis functions of the inner approximation Xn of X as the test functions in problem (6) , allows us to 22 basics of reduced order models identify the latter as a system of equations. If we set uk h(γ) the algebraic representative of uk h(γ), and we obtain for k=1, . . . , Nt, a linear system that we have to solve: (M(γ) + A(γ))uk h(γ) = M(γ)uk−1 h(γ) + bk(γ), (8) where M(γ) is the algebraic representative of the time derivative M , A(γ) is the algebraic representative of the elliptic operator and Neumann boundary condition, and bksummarizes the effect of function Fk(·,γ). Here, we have considered, like in model (1) an explicit Euler scheme for time discretization, but again, extensions to more general numeric schemes is possible. If we were to consider a PDE that does not depend on time (4) , we just need to get rid of the algebraic representatives of the time derivative in (8) . In that case, we would obtain the following linear system A(γ)uh(γ) = b(γ), (9) where A(γ) is the algebraic representative of the elliptic operator and Neumann boundary condition, and bsummarizes the effect of function f(γ). Every algebraic representation (7) , (8) and (9) , even in the case when there are no nonlinearities involved, could involve the resolution of matrix-vector multiplications or linear systems of very high dimension. That can be very demanding in terms of computational time and storage. In order to solve this issue, some techniques have been developed. We explain them more in detail in Section 1.3. Algebraic representatives of type (7) will be present in the first part of this dissertation, Part I, when modelling the neural activity by means of networks of coupled oscillators. While a representative of type (8) will be present in Part II, when solving an unsteady incompressible flow. And a representative of type (9) will be present in Part III, when computing the intrinsic Proper Generalized Decomposition expansion in a Grassmann manifold. 1.3 reduced order modelling In this section, we focus on the argumentation of the need of the development of ROM techniques, the main theme of this dissertation, as well as on the presentation of the techniques that we will be using in the following chapters. 1.3 reduced order modelling 23 Although the use of ROM techniques is suitable for both, systems of parametric differential equations and parametric PDEs, the philosophy behind their use is quite different: • In the case of a system of parametric differential equations, like the one presented in (2) , our objective will be to reduce the number of equations considered in the original problem, expressing the solution of the rest as a linear combination of the solved ones. These techniques are some of the so-called hyper-reduction techniques (Farhat et al., 2020; Ryckelynck, 2005). • In the case of parametric PDEs, like the one presented in (3) , the objective behind ROM techniques will be the representation of the solution of the original problem as a linear combination of some basis functions that keep the most of the information of the solution. These techniques include the so-called Reduced Basis Method (Ali et al., 2020; Porsching & Lee, 1987; Quarteroni et al., 2016). In the following, we present some of the most common methods for the construction of the reduced problem, that we will be using along this dissertation. These methods provide, both, the basis and dimension of the reduced problem. 1.3.1Proper Orthogonal Decomposition The POD method is based on the Singular Value Decomposition (SVD) technique, and it is one of the most popular approaches for construction of the reduced model (Kunisch & Volkwein, 2002; Pinnau, 2008; Volkwein, 2013). In the following, we will consider that we have a set of solutions of the original problem, either (1) or (6) , S= (s1, . . . , sNs)∈RN×Ns , where N denotes the number of unknowns of problem (1) and the dimension of the discrete space in problem (6) . Let us denote by snapshot a solution of the problem for a given parameter setting, si as the i-th snapshot of the problem, and Sas the snapshot matrix. It is known that the following property holds for the SVD of the snapshot matrix S (Eckart & Young, 1936; Quarteroni et al., 2016). 24 basics of reduced order models Proposition 1.2 Let VNr={W∈RN×Nr:WTW=INr} be the set of all orthonormal bases of Nr - dimensional subspaces of RN. Then, Ns ∑ i=1 ∥si−UUTsi∥2 2=min W∈VNr Ns ∑ i=1 ∥si−WWTsi∥2 2= r ∑ i=Nr+1 σ2 i, where, σi denotes the i -th singular value, U∈RN×Nr arises from the SVD of S and is a full matrix containing the first Nr left singular vectors in columns corresponding to the Nr greatest singular values, and r is the number of non-zero singular values. Indeed, Proposition 1.2holds for the SVD of any matrix. Therefore, the SVD provides a basis that minimizes the sum of the square of the errors between each snapshot vector wiand its projection onto the reduced subspace, and the error is equal to the sum of the squares of the neglected singular values. Furthermore, Proposition 1.2provides a natural way to select the dimension of the reduced model, as we can select the higher singular values until a certain tolerance criteria is fulfilled. The general POD procedure to obtain the reduced basis and dimension, either for problem (1), or problem (6), can be summarized as follows: 1. First, we need to solve the original problem for several instances of the parameter vector γ∈ Dtrain ∈ D , where the solution is qualitatively similar, to obtain the snapshots in time, and store them in the snapshot matrix, S. 2. Then, we have to compute the SVD of the snapshot matrix S=UΣVT , and sort, in descending order, the singular values and its corresponding left singular vectors. 3. Finally, we can apply Proposition 1.2as a measure of the information lost to select the least amount of basis functions that retain the most information, according to a certain fixed tolerance. For example, if we define ϵtol >0 as the tolerance, then we can use a suitable criterion involving the information retained by each mode, that is, I(Nr) = ∑Nr i=1σ2 i ∑r i=1σ2 i ≥1−ϵ2 tol, 1.3 reduced order modelling 25 or we can choose the first integer Nr such that σ2 Nr≤ϵ2 tol . This will provide the dimension of the reduced model Nr , along with a submatrix of U , the matrix of left singular vectors, that contains the basis functions, {ϕi}Nr i=1. There are two main issues regarding the POD technique. On the one hand, we need to know the solution for several values of the parameter Dtrain , and, as said before, we are trying to avoid the resolution of the computationally expensive original problem. On the other hand, as we have no way to determine the best size for Dtrain , it can be extremely large, and even if we can afford the resolution of the original problem, the computation of the SVD of the snapshot matrix could also be very challenging. The POD technique will be present in the two first parts of the dissertation. To begin with, in Part I, we develop a technique based on the POD that allows to retain the information of the structure of problems like (2) . After that, in Part II, we will use a coupled strategy with the Greedy algorithm (see Section 1.3.2), the so-called coupled POD+Greedy algorithm, that we will introduce in more detail in Section 1.3.3. Now, we proceed to explain another method for the construction of the reduced problem, the Greedy Algorithm. 1.3.2Greedy Algorithm In order to overcome the issue of having to solve the original problem for several instances of the parameter vector γ , we could think of an optimal way to sample D , avoiding repeating information that we have already captured. The Greedy Algorithm is an iterative algorithm that has been proven useful for the construction of the reduced problem (Binev et al., 2011; Buffa et al., 2012; Quarteroni et al., 2016). The Greedy Algorithm is based on iteratively enriching the reduced problem to obtain better approximations of the original solution. In each iteration of the Greedy Algorithm, we add to the reduced problem the information of the original solution of the parameter vector value that provides the larger error between the original solution and the solution of the reduced problem. In the following, we consider that we are working with an elliptic problem (4) or with steady solutions of problem (6) . For the case of unsteady solutions we refer to Section 1.3.3, where a coupled strategy with the POD technique is presented. We also assume that we have a discrete subset of the parameter space Dtrain ⊂ D. The general Greedy Algorithm to obtain the reduced basis and dimension, can be summarized as follows: 26 basics of reduced order models 1. At the first iteration, select an initial parameter vector value γ1∈ Dtrain , define the set of parameter values S1={γ1}, and the reduced space X0={0}. 2. Then, at any iteration i=1, . . . : a) Compute the (steady) solution of the original problem, uh(γi). b) Add the computed solution to the reduced space, Xi=span{Xi−1,uh(γi)}. c) Solve the reduced problem for every parameter value γ∈ D to obtain ui(γ) . d) Compute the errors between the original solution and the solutions of the reduced problem ∥uh(γ)−ui(γ)∥ , for every γ∈ Dtrain , where ∥·∥ is the norm of the corresponding space. e) If maxγ∈Dtrain ∥uh(γ)−ui(γ)∥ is lower than a certain tolerance ϵtol , then stop the algorithm. Otherwise, set γi+1=argmaxγ∈Dtrain ∥uh(γ)−ui(γ)∥ , let Si+1=Si∪ {γi+1}, and repeat from instruction a). However, as it states, the computation of the error ∥uh(γ)−ui(γ)∥ at each step of the Greedy Algorithm still requires the knowledge of the solution of the original problem for every γ∈ Dtrain . So, the general Greedy Algorithm does not prove an improvement against the POD technique. Bearing that in mind, in practice, instead of computing the exact error, we can consider an a posteriori error estimator ∆i(γ) , whose computation is cheaper with respect to the computation of the exact error. The resulting algorithm, changing the exact error by an a posteriori error estimator, is the so-called Weak Greedy Algorithm (Binev et al., 2011; Quarteroni et al., 2016), and it is the usual one for the construction of the reduced problem. A further difference with respect to the POD technique is that the basis obtained for the reduced space is not orthonormal by construction. Therefore, in order to avoid redundancies and possible numerical issues, the reduced space should be orthonormalized with respect to the norm of the original space. For this purpose, we will use the Gram-Schmidt orthonormalization process when needed. The (Weak) Greedy Algorithm will be present in the two latter parts of this dissertation. On the one hand, in Part II, we will use a coupled strategy with the POD technique (see Section 1.3.1) for unsteady problems, that we will explain in more detail in Section 1.3.3. On the other hand, in Part III, a Greedy algorithm based on the Gradient Descent Algorithm in Grassmann manifolds is presented to obtain the modes of the Proper Generalized Decomposition expansion. We will explain this procedure more Chapter 7. 1.3 reduced order modelling 27 1.3.3POD+Greedy algorithm The POD technique (Section 1.3.1) and the Greedy algorithm (Section 1.3.2) can be applied successfully in the case of elliptic or steady problems depending on parameters with no further complexity. However, in the case of unsteady problems, the POD technique can be unaffordable because of the size of the snapshots we will have to consider, and the Greedy Algorithm may converge very slowly due to the variability of the snapshots. For those reasons, we introduce a coupled strategy for unsteady problems combining both, the POD technique and the Greedy Algorithm. We call it the POD+Greedy strategy. The main idea behind this coupled strategy is to use the POD technique considering the time t as the parameter, and the Greedy Algorithm for the parameter vector γ . The whole process for the coupled POD+Greedy strategy can be summarized as follows: 1. At the first iteration, select an initial parameter vector value γ1∈ Dtrain , define the set of parameter values S1={γ1}, and the reduced space X0={0}. 2. Then, at any iteration i=1, . . . : a) Compute the solution of the original problem for an interval of time, {uk h(γi)}Nt k=1, with uk h(γi) = uh(tk,γi), and tk=k∆t. b) Update the reduced space by means of the POD technique, taking into account the previous reduced space Xi−1 and the information about the solution for the new parameter. We obtain the next reduced space Xi of dimension N. c) Solve the reduced problem for every parameter value γ∈ D to obtain uN(γ). d) Compute the errors between the original solution and the solutions of the reduced problem ∥uh(γ)−uN(γ)∥ , for every γ∈ Dtrain , where ∥ · ∥ is the norm of the corresponding space. e) If maxγ∈Dtrain ∥uh(γ)−uN(γ)∥ is lower than a certain tolerance ϵtol , then stop the algorithm. Otherwise, set γi+1=argmaxγ∈Dtrain ∥uh(γ)−uN(γ)∥ , let Si+1=Si∪ {γi+1}, and repeat from instruction a). 34 mathematical models addressed (a) Attractive relaxation limit cycle of the single-cell ICC model together with the S-shaped critical manifold: single and double arrows materialize the slow and fast parts of the limit cycle, respectively. (b) Instance of an orbit of the single-cell ICC model for the case of existence of an attractive equilibrium (black point on the left sheet of the critical manifold). The system does not generate sustained oscillation. (c) MMOs generated by the single-cell ICC model. Left panel: Periodic orbit together with the critical manifold. Right panel: magnified view of the small oscillations near the lower fold of the critical manifold. Figure 2.1: Asymptotic behaviors of system (10) depending on the location of its singular point together with the S-shaped critical manifold. Single and double arrows materialize the slow and fast parts of the limit cycle, respectively. 2.1 models of neuronal intracellular calcium concentration dynamics 35 2.1.2Two coupled-cells model of neuronal Intracellular Calcium Concentration dynamics The two-cells model is built as a linear coupling of two copies of system (10) . In the following, as in Fernández-García and Vidal (2020) and in Bandera Moreno et al. (2022), we will denote (x1,y1,z1) -subsystem as O1 , which depends on x2 , and (x2,y2,z2)-subsystem as O2, which depends on x1. The system then reads as O1           ˙ x1=τ(−y1+f(x1)−ϕf(z1)), ˙ y1=τεk1(x1+a1y1+a2+c(x1−x2)), ˙ z1=τε ϕr(x1)−z1−zb τz, O2           ˙ x2=τ(−y2+f(x2)−ϕf(z2)), ˙ y2=τεk2(x2+a1y2+a2+c(x2−x1)), ˙ z2=τε ϕr(x2)−z2−zb τz,                                (11) with the coupling parameter c∈[−1,1] and where parameters ki,i=1,2 , state the possible heterogeneity between both oscillators, O1and O2. Let us describe the main behavioral patterns obtained in system (11) . We distinguish two cases, on the one hand, we assume that the two cells are identical, that is, ki=1 , i=1,2 , what we call the homogeneous case, that we explain in the following. On the other hand, we let the two cells to be slightly different, that is, ki=1 , for some i=1,2 and k1=k2, what we call the heterogeneous case, and explain afterwards. Homogeneous case The case with ki=1,i=1,2 , that is, the homogeneous case, was already studied in Fernández-García and Vidal (2020). In that work, the authors outline the different synchronization types between O1 and O2 within the attractive limit cycle of system (11) depending on the coupling strength c=0. Here, we summarize the main asymptotic behaviors obtained in Fernández-García and Vidal (2020). We illustrate each case in figures 2.2-2.4. 1. Total Oscillation Death (TOD): O1 and O2 reach asymptotically a stable equilibrium. Figure 2.2ashows the behavior of system (11) in the time interval [0,30] for c=−0.7. 36 mathematical models addressed 2. Relaxation Loss (RL): one oscillator produces relaxation oscillations, i.e. the orbit approaches both Sl and Sr , while the other one produces small oscillations and the orbit remains near Sl . Figure 2.2bshows the behavior of system (11) in the time interval [0, 30]for c=−0.502. 0 5 10 15 20 25 30 -4 -2 0 2 4 x1(t),x2(t) O1O2 0 5 10 15 20 25 30 1 2 3 4 5 z1(t),z2(t) (a) Total oscillation death, with c=−0.7 . The orbit tends to a stable equilibrium. 0 5 10 15 20 25 30 -2 0 2 x1(t), x2(t) O1O2 0 5 10 15 20 25 30 1 2 3 4 5 z1(t), z2(t) (b) Relaxation loss, with c=−0.502 . Only O1 generates spikes in the corresponding ICC time series, O2 does not undergo relaxation. Figure 2.2: The x and z time series (top and bottom panel respectively) of each oscillator O1 and O2 are obtained from the two-cells model (11) for k1=k2=1 and different coupling parameters c<0. 0 5 10 15 20 25 30 -2 0 2 x1(t), x2(t) O1O2 0 5 10 15 20 25 30 1 2 3 4 z1(t),z2(t) (a) Antiphase synchronization, with c=−0.25 . Each spike of each oscillator occurs precisely at the middle time between two successive spikes of the other oscillator. 0 5 10 15 20 25 30 -2 0 2 x1(t),x2(t) O1O2 0 5 10 15 20 25 30 1 2 3 4 z1(t),z2(t) (b) No correlation between the cells, with c=0 . Both oscillators reach asymptotically the exact same orbit with MMOs. Figure 2.3: The x and z time series (top and bottom panel respectively) of each oscillator O1 and O2 are obtained from the two-cells model (11) for k1=k2=1 and different coupling parameters c≤0. 2.1 models of neuronal intracellular calcium concentration dynamics 37 3. Antiphase Synchronization (AS): O1 and O2 reach asymptotically the same orbit, and the phase-shift between them is exactly one half-period. Figure 2.3ashows the behavior of system (11) in the time interval [0,30]for c=−0.25. 4. Uncoupled case (U): when c=0 the oscillators O1 and O2 reach asymptotically the exact same orbit with MMOs. There is a phase shift between the oscillators, which depends on the initial conditions. The phase-shift will vanish if the initial conditions are identical. Figure 2.3bshows the behavior of system (11) in the time interval [0, 30]for c=0. 5. Almost-in-Phase synchronization (AiP): system (11) admits an attractive limit cycle along which O1 and O2 generate MMOs. The two 3 D oscillators are synchronized with a small phase-shift. Figure 2.4ashows the behavior of system (11) in the time interval [15,45] for c=0.1 , as well as a zoom of one of the peaks where we can see that the oscillators are not perfectly in phase. 6. In-phase Locking synchronization (IL): O1 and O2 reach asymptotically the exact same orbit which is precisely the cycle existing for the uncoupled case. Therefore, the attractive limit cycle of system (11) lies on {x1=x2,y1=y2,z1=z2} and is a global attractor of the system. Figure 2.4bshows the behavior of system (11) in the time interval [15, 45]for c=1, as well as a zoom of one of the peaks. In Fernández-García and Vidal (2020), the distribution of these behaviors according to the coupling parameter c has been theoretically studied and numerically shown for symmetric coupled system (11) (with k1=k2). Furthermore, through numerical exploration and heuristic arguments, the authors pinpointed distinctive behaviors emerging within narrow ranges of c values, representing transitions between the aforementioned primary behaviors. To gain deeper insights into the organization of these behaviors and the significance of the most exceptional ones, examining a broader scenario and removing the constraint of symmetry between the two oscillators proves beneficial. Although the theoretical basis relying on this symmetry loses validity, the evanescent behaviors identified in Fernández-García and Vidal (2020) become central behaviors when considering k1=k2 . The role of the heterogeneity parameters and of the coupling parameter in organizing the dynamical structure can be clarified through a numerical investigation. We will discuss this point in more detail in the subsequent section. 38 mathematical models addressed 15 20 25 30 35 40 45 -2 0 2 x1(t),x2(t) O1O2 15 20 25 30 35 40 45 1 2 3 z1(t),z2(t) 20 22 24 26 28 30 -2 0 2 x1(t),x2(t) O1O2 20 22 24 26 28 30 1 2 3 z1(t),z2(t) (a) Almost-in-phase synchronization, with c=0.1 . The two oscillators are synchronized with a slight phase-shift between their spikes. 15 20 25 30 35 40 45 -2 -1 0 1 2 x1(t),x2(t) O1O2 15 20 25 30 35 40 45 1 2 3 z1(t),z2(t) 20 22 24 26 28 30 -2 -1 0 1 2 x1(t),x2(t) O1O2 20 22 24 26 28 30 1 2 3 z1(t),z2(t) (b) In-phase synchronization, with c=1 . The two oscillators spikes are perfectly synchronized. More precisely, the time difference between corresponding spikes of O1and O2tends asymptotically to 0. Figure 2.4: The x and z time series (top and bottom panel respectively) of each oscillator O1 and O2 are obtained from the two-cells model (11)for k1=k2=1and c >0. Heterogeneous case If some heterogeneity is introduced between the two oscillators by considering different values for the recovery variable times-scale parameters k1 and k2 , then, the uncoupled oscillators ( O1 and O2 with c=0 ) admit attractive limit cycles with different intrinsic frequencies. Moreover, the oscillatory time series generated by the ICC variables z1 and z2 along the corresponding limit cycles feature different amplitudes. The theoretical study performed in Fernández-García and Vidal (2020) for the homogeneous case is strongly based on symmetry, therefore it cannot be directly adapted to this case. We present here, the results obtained using numerical simulations for the study of the heterogeneous case, performed in Bandera Moreno et al. (2022). 2.1 models of neuronal intracellular calcium concentration dynamics 39 In the following, we fix k1=1 and consider k2=k>1 , since the results for k<1 can be deduced by exchanging the roles of O1 and O2 . We further explain the synchronization patterns obtained in the heterogeneous case. • For c>0 , while k∈(1,2] , we obtain the same qualitative behavior as in the homogeneous case (k=1). • For c<0 , we find the three behaviors obtained in the homogeneous case, but in the antiphase region c∈[c0+δ,−δ] with c0≈ −0.5 , and 0<δ≪1 , besides the antiphase behavior, we obtain new synchronization patterns for k∈(1,3] that we explain subsequently. First, we need to present a definition that will help us to distinguish between the different behaviors. Definition 2.1 In the following, we say that the signature of a periodic orbit of system (11) is s/1 if, along one period of this orbit, one 3D oscillator spikes s times while the other one spikes once. Hence, with this notation, the antiphase behavior of the homogeneous case, (see Figure 2.3a), would be denoted as 1/1 . Let us focus on case c=−0.25 and describe the main and transitory patterns. • For k∈(1,1.5563] , we obtain a behavior qualitatively similar to the one obtained in the homogeneous case, one oscillation for one oscillation (signature 1/1 ), see Figure 2.5afor a specific instance where k=1.5 . However, as we can observe, the behavior is not in antiphase, as it is the case in the homogeneous case (see Figure 2.3a), because the phase shift is not exactly half of the period. This is an indicator of how the heterogeneity is affecting the system. We refer to the behavior shown in Figure 2.5aas "phase-shifted synchronization". • For k∈[1.6468,2.7350] , one oscillator generates two pulses while the other one generates a single one along a limit cycle period, which corresponds to signature 2/1. See Figure 2.5bfor a specific instance where k=2. • For k∈[2.7582,3] , the behavior of the system changes and one oscillator generates three pulses while the other one generates a single one along a limit cycle period, which corresponds to signature 3/1 . See Figure 2.5cfor a specific instance where k=3. 40 mathematical models addressed 0 5 10 15 20 25 30 -2 0 2 x1(t), x2(t) O1O2 0 5 10 15 20 25 30 1 2 3 4 z1(t), z2(t) (a) Signature 1/1, with k=1.5 . Each oscillator spikes a single time between two consecutive spikes of the other one. 0 5 10 15 20 25 30 -2 0 2 x1(t), x2(t) O1O2 0 5 10 15 20 25 30 1 2 3 4 z1(t), z2(t) (b) Signature 2/1, with k=2 . Oscillator O2 spikes precisely two times between two consecutive spikes of oscillator O1. 0 5 10 15 20 -2 0 2 x1(t), x2(t) O1O2 0 5 10 15 20 1 2 3 4 z1(t), z2(t) (c) Signature 3/1, with k=3 . Oscillator O2 spikes precisely three times between two consecutive spikes of oscillator O1. Figure 2.5: Phase-shifted synchronization with different signatures. The time series are generated with the two-cells model (11)with k1=1and varying k2=k. • For k∈[1.5592,1.6454] , in the transition from signature 1/1 to 2/1 , the system can exhibit a mixture of two signatures: between two successive spikes of O1 , oscillator O2 spikes one and two times, alternatively. Therefore, two signatures 1/1 and 2/1 alternate along the limit cycle. See left panel of Figure 2.6. • For k∈[2.7372,2.7563] , in the transition from signature 2/1 to 3/1 , a similar case occurs: between two successive spikes of O1 , oscillator O2 spikes two and three times, alternatively. Therefore, two signatures 2/1 and 3/1 alternate along the limit cycle. See Figure 2.6, right. 2.1 models of neuronal intracellular calcium concentration dynamics 41 0 5 10 15 20 25 30 -2 0 2 x1(t), x2(t) O1O2 0 5 10 15 20 25 30 1 2 3 4 z1(t), z2(t) 0 5 10 15 20 25 30 -2 0 2 x1(t), x2(t) O1O2 0 5 10 15 20 25 30 1 2 3 4 z1(t), z2(t) Figure 2.6: Transitory signatures of system (11) outputs while varying k . Parameter k1=1 . Left panel: for k=1.6 , the signature alternates between 1/1 and 2/1 . Right panel: for k=2.74 , the signature alternates between 2/1 and 3/1. In Figure 2.7, we represent the signature of system (11) limit cycle for c=−0.05 (left panel), c=−0.25 (central panel) and c=−0.40 (right panel) according to the value of parameter k . The blue lines represent the signature value when system (11) admits a limit cycle, while gray bands materialize the intervals of k value where the signature is not well-defined, corresponding to signature transitions. We note that the representation of the different behaviors in Figure 2.7exhibits a Devil Staircase-like behavior, characterized by different stable behaviors for large interval of one parameter values, separated by transition intervals where the system can exhibit a mixture of these stable behaviors, and even chaotic behaviors, (Bak, 1986). 1 1.5 2 2.5 3 k 1 1.5 2 2.5 3 s 1/1 2/1 3/1 1 1.5 2 2.5 3 k 1 1.5 2 2.5 3 s 1/1 2/1 3/1 1 1.5 2 2.5 3 k 1 1.5 2 2.5 3 3.5 4 4.5 5 s 1/1 2/1 3/1 4/1 5/1 Figure 2.7: Signature of the two-cells model (11) limit cycle according to k values for c=−0.05 (left panel), c=−0.25 (central panel) and c=−0.40 (right panel). The gray bands materialize the intervals of k values for which the signature is not well-defined. Note that there are some ranges of the parameter values inside the transitions where the system does not exhibit one of the behaviors described above, but a more complex combination of the main behaviors. In fact, we also have a canard explosion, 42 mathematical models addressed that is, a very fast transition from small oscillations to relaxation oscillations, with a small variation of the parameters (Benoit et al., 1981; Krupa & Szmolyan, 2001). This transition happens via a family of cycles whose characteristic property is that they contain canard segments, i.e. segments of the orbit that follow the repulsive slow manifold Smfor a time O(1). Now, that we have detailed the observed patterns in the case c=−0.25 , let us come back to the general case. We study the behaviors for k∈(1,3] in the range of the coupling parameter c∈[−0.4, 0), we obtain the subdomains in the (c,k)-plane where the system (11) exhibits each main behavior, as shown in Figure 2.8right, where the region c∈[−0.4, −0.05] and k∈[1,3] is displayed. We can observe that the transition intervals shrink while the coupling parameter values are lowered. -0.5 -0.48 -0.46 -0.44 -0.42 -0.4 c 1 1.5 2 2.5 3 3.5 k RL 2/1 1/1 3/1 4/1 . . . -0.4 -0.35 -0.3 -0.25 -0.2 -0.15 -0.1 c 1 1.5 2 2.5 3 k 1/1 3/1 4/1 2/1 Figure 2.8: Regions of (c,k) -plane associated with fixed signature value of the attractive limit cycle of the two-cells model (11) . The left panel shows that the intervals of k values corresponding to a given signature converge to the singleton {1} , while c decreases to −0.5 . The right panel emphasizes the small strips in (c,k)-plane for which the signature is not well-defined. Finally, we focus now on the range near the transition to the relaxation loss behavior, namely c∈[−0.5, −0.4] . All ranges of the main behaviors converge to the transition point for k=1 , see Figure 2.8left, showing c∈[−0.5, −0.4] and k∈[1,3.5] . We can also observe that the transition to the relaxation loss regime depends on the heterogeneity parameter k. The ICC models will be used in Part Iof this dissertation. Where we present the two clusters model, along with its main behaviors. However, simulating large-scale network models can be computationally expensive, and while the POD method approximate accurately the solution of the original system (Bandera Moreno et al., 2022), they may lead to the loss of the problem’s original structure. Therefore, the main objective of 2.2 reduced order model for the smagorinsky turbulence model 43 Part Iwill be to develop a reduced order model procedure based on the POD, in order to preserve the two clusters structure of model. 2.2 reduced order model for the smagorinsky turbulence model This section will be devoted to the presentation of the basic Large Eddy Simulation (LES) turbulence model, the celebrated Smagorinsky turbulence model introduced in Smagorinsky (1963). In this model, the effect of the subgrid scales on the resolved scales is modeled by eddy diffusion terms (Chacón Rebollo & Lewandowski, 2014; Sagaut, 2005; Smagorinsky, 1963), that is, a characterization of the transport and dissipation of energy in the smaller scales. Although being intrinsically discrete, we present in a first place a continuous model whose discretization can be interpreted as the Smagorinsky turbulence model. The presentation of the continuous model is included to clarify the relationship between the Navier-Stokes model with the Smagorinsky turbulence model. Then, we will derive the respective reduced order model. 2.2.1Smagorinsky turbulence model In this section, we present the Smagorinsky model as a continuous one, although it is intrinsically discrete, associated to a given discretization grid. We will just specify its dependency on the actual grid used in the expression of the eddy viscosity terms. Then, we present the Smagorinsky turbulence model in the discrete framework and some of its properties. We thus consider a polygonal, if d=2 , or polyhedral, if d=3 , bounded domain Ω⊂Rd . We assume that its boundary is split into Γ=ΓD∪ΓN , where ΓD is the boundary relative to the Dirichlet boundary conditions, and ΓN to the Neumann conditions. In order to specify the expression of the eddy viscosity, we will also consider a triangular mesh {Th}h>0 of Ω , that we assume to be regular in the sense of Ciarlet (Ciarlet, 2002). In the following, we will consider that the model depends on just one parameter, the Reynolds number, denoted by γ∈ D . Here, the parametric set D , stands for the positive real numbers. The Reynolds numbers should be large enough to ensure that the flow is in well-developed turbulent regime. For clarity, we will also drop the dependence of uand pon space and time. INTRODUCTION In this part, our emphasis is on the study and simulation of network models of neuronal activity that involve different timescales, modelled by a system of Ordinary Differential Equations (ODEs). In particular, we focus on the Proper Orthogonal Decomposition (POD) method, with a specific goal of addressing one of the limitations that it presents, namely, the potential loss of structure from the original model. This part aligns with the Spanish national project HYROM. The research presented in this part of the thesis is an extended version of the results published in the research paper Bandera Moreno et al. (2022). One of the main challenges of network models of neuronal activity is the difficulty to classify and analyze the synchronization features between network cells, and the cost of reliable simulations for networks of realistic size. Furthermore, these models usually involve different timescales, what is naturally translated into the models by the use of multiple timescales systems as cell dynamics to be coupled as a network. Synchronization mechanisms of multiple timescales systems are quite specific and different from that of harmonic oscillators (especially the impact of the coupling types and strengths on the synchronization). Nevertheless, it is well-known that synchronization of neuron assemblies plays a key role in the neuronal information treatment and diffusion, and overall in the neural behavior. With the aim of reducing the computational complexity of the simulations of this kind of models, in this part, we will focus on the development of reduced order models for network models of neuronal activity. Despite the ability of Reduced Order Modelling (ROM) techniques to capture the main information of these systems (see Bandera Moreno et al. (2022) for an approach that joints POD with the Empirical Interpolation Method (EIM) to a model of intracellular calcium concentration), they fail in preserving the separation of timescales of the original model, what we call the structure of the model. These could play a key role when dealing with systems with more than two timescales or where we want to keep that separation of timescales also in the reduced model. As a motivation, we provide a simple example of the lost of structure. 51 52 introduction In the following, we will consider a dynamical system of N∈N equations, that could be separated into two groups (for example, see Equation (2) ). However, this analysis can be generalized to M<N groups. Without loss of generality, we can consider that the variables can be sorted by group. A dynamical system in this framework could read, ˙x1 ˙x2!= A11 A12 A21 A22 ! x1 x2!+ b1 b2!, (19) where xi and bi,i=1,2, are vectors of dimensions N1 and N2 , respectively, such that N=N1+N2and Aij,i,j=1,2, are real matrices with consistent dimensions. For the sake of clarity, we will consider A11 and A22 as identity matrices multiplied by nonzero real scalars, aand b, respectively, A21 =0 and A12 a full matrix, so A11 A12 A21 A22 != aI A12 ObI!. If we consider system (19) , and we let U be the matrix arising from the POD of some snapshots of the system (see Section 1.3.1for more information), then the reduced order system reads, ˙z=˜ Az+br, (20) where br=UTb= UT 11 UT 21 UT 12 UT 22 ! b1 b2!, and ˜ A=UTAU = UT 11 UT 21 UT 12 UT 22 ! aI A12 ObI! U11 U12 U21 U22 ! = aI+UT 11A11U21 UT 11A12U22 UT 12A12U21 bI+UT 12A12U22 !, and we have lost the structural information of the original problem and the possible separation between variables, for example, when a=1 and b=ε≪1. In order to avoid this lost of structure, we will develop a novel technique that will take into account the structure of the original problem when constructing the reduced order model. That will allow us to preserve the separation of timescales. We will also, introduction 53 provide a data-driven technique in order to determine the structure of the original model, which can be difficult to determine directly. This part is structured as follows: in Chapter 3, we will present a novel technique based on the POD technique that will preserve the structure of the original problem, and we also provide a data-driven technique in order to determine the best structure for the reduced model in order to maintain as most structure of the original model while still providing a low number of reduced variables. Then, in Chapter 4, we will compare these novel techniques with the original POD for three neural network models with multiple timescales. 3 STRUCTURE-PRESERVING METHODS BASED ON THE PROPER ORTHOGONAL DECOMPOSITION In the introduction of this part, we have presented the problematic of the lost of structure when trying to obtain a reduced order model of a dynamical system directly by means of the original POD technique that we presented in Section 1.3.1. To overcome this issue, in this chapter, we develop an extension to the original POD technique, that is able to preserve the structure of the original system. We will focus on the framework of a system of parametric ODEs, although this technique can be also applied in the field of PDEs that exhibit any qualitative separation in the variables. The structure of the chapter is as follows: in Section 3.1, we present the extension of the original POD to systems that exhibit a qualitative separation of the variables, and we also give some mathematical results concerning this novel technique in Section 3.1.1. Then, in Section 3.2, we present a previous data-driven step that can be performed in order to obtain the best structure for the reduced model in order to maintain as most structure of the original model while still providing a low number of reduced variables. 3.1 proper orthogonal block decomposition If we take into account the construction of the reduced model (20) , we can see that, in order to maintain the main properties of the original model when we perform the reduction, we have to take care of the structure of the matrix U arising from the POD technique. Taking into account the separation of variables in the original model (19) and the action of matrix U in the construction of the reduced model (20) , we can see that if matrix Uwere diagonal by blocks, we could keep the structure of the original model. 55 56 structure-preserving methods based on the pod With that in mind, we divide the snapshot matrix S= (w1|. . . |wNs) into two blocks, as follows, ST=ST 1|ST 2 . Now, we perform a Singular Value Decomposition (SVD) for each block, and we obtain a POD basis for each block, S1≈ r1 ∑ i=1 σ(1) iu(1) i(v(1) i)T,S2≈ r2 ∑ i=1 σ(2) iu(2) i(v(2) i)T, where σ(·) i stands for the i -th singular value, u(·) i is the i -th left singular vector, and v(·) i is the i -th right singular vector of the corresponding matrix S· . Proposition 1.2 holds for each block, so as in the original POD method, we have a suitable criterion for the selection of the reduced basis dimension for each variable. We name the reduced dimensions Nr1and Nr2for each block, respectively. Now, we focus on the reconstruction of the original snapshot matrix as a whole. We note that each block matrix has the correct number of columns, so no change in the block right singular vectors v(·) i is needed. However, each block left singular vector u(1) i lays in RN1 , so we need to complete them in order to have a vector in RN as we must, and respectively for vectors u(2) i. A direct concatenation of the block left singular vectors does not work, as we do not know how the product u(k) i(v(l) i)T behaves, and we have no control over it. A successful approach is to add as many zeros as needed to each block left singular vector, taking into account the position of the block matrix in the original matrix, i.e., u(1) i7→ u(1) i 0N2!,u(2) i7→ 0N1 u(2) i!. With this decomposition, we can recover the original matrix as U1ON1,Nr2 ON2,Nr1U2! Σ1ON1,Nr2 ON2,Nr1Σ2! VT 1 VT 2!= U1Σ1VT 1 U2Σ2VT 2!≈ S1 S2!. As we have exploited the separation of the variables by blocks, we name this method the Proper Orthogonal Block Decomposition (POBD). Here, we have considered just two blocks for simplicity, but the method can be extended to any number of blocks Nb. 3.1 proper orthogonal block decomposition 57 3.1.1Properties of the Proper Orthogonal Block Decomposition In this section, we present some properties of the POBD. On the one hand, we will check that this novel technique is able to preserve the structure of the original problem. On the other hand, we will present a result concerning the error bound for the neglected information. We focus first on the preservation of the structure of the original problem, if we select the matrix arising from the POBD U= U1ON1,Nr2 ON2,Nr1U2!, (21) as the new basis for the reduced subspace, the reduced model associated to the original system (19) can be obtained as follows: ˜ A=UTAU = UT 1OT N2,Nr1 OT N1,Nr2UT 2! IN1A12 ON2,N1IN2! U1ON1,Nr2 ON2,Nr1U2! = INr1UT 1A12U2 ONr2,Nr1INr2!= INr1˜ A12 ONr2,Nr1INr2!. We remark that applying the POBD technique allows us to maintain the structure of the original system, that was the aim of this technique. Now, we focus on obtaining an error bound on the neglected information when using the POBD for the construction of the reduced model. As done for the original POD technique, where we were able to obtain an error bound depending on the values of the neglected singular values, thanks to the properties of the SVD of a matrix, see Proposition 1.2, in our case, we are able to obtain a similar error bound applying similar techniques. Proposition 3.1 Let U be an orthonormal matrix built as in (21) for a dynamical system with Nb≥2 blocks. Then, the next error bound holds NS ∑ i=1 ∥wi−UTUwi∥2 2= Nb ∑ j=1 rj ∑ i=Nrj+1 (σ(j) i)2. 58 structure-preserving methods based on the pod Proof : This proof is an adaptation of the proof of Proposition 1.2in (Quarteroni et al., 2016) to our case. We detail it for Nb=2 for the sake of clarity. For a matrix, A∈RN×NS we have ∥A∥2 F=∑NS i=1∥ai∥2 2 , where ai denotes the i -th column of A. Thus, in our case we can write NS ∑ i=1 ∥wi−UUTwi∥2 2=∥S−UUTS∥2 F. (22) Now, as S= S1 S2!= S1 ON2,NS!+ ON1,NS S2!, and UUTS= U1UT 1S1 U2UT 2S2!= U1UT 1S1 ON2,NS!+ ON1,NS U2UT 2S2!, along with the facts that ∥S∥F=∥ST∥F and ∥(A|B)∥2 F=∥A∥2 F+∥B∥2 F , we can decompose the last term of (22) as ∥S−UUTS∥2 F=∥S1−U1UT 1S1∥2 F+∥S2−U2UT 2S2∥2 F. Finally, applying Proposition 1.2to each term on the right-hand side, we obtain the desired result. A similar approach for the decomposition of the domain in computational simulation of fluids was introduced in Anderson et al. (2022). Property 3.1provides a similar error bound to the one obtained in Property 1.2, but, in the former, the action of the blocks is taken into account. We note that the more blocks that we consider the greater the error would be, this will be tackled in the next section, Section 3.2. Now, we give two remarks concerning the properties of the POBD. The first one acknowledging that the original POD technique is a special case of the POBD, and another one about the non-optimality of the error bound for two or more blocks. Remark 3.1 The POBD procedure is suitable for any number of blocks Nb<N and when considering just 1block, we recover the classic POD method. 3.2 automatic proper orthogonal block decomposition 59 Remark 3.2 We remark that despite the similarities between the error bounds of Propositions 1.2and 3.1, for a number of blocks Nb≥2 the latter might not be the optimal for reducing onto a subspace of dimension ∑Nb j=1Nrj. In order to measure the good behavior of this method, taking into account just the number of equations of the reduced system could be not enough, as we could also need to retain some information about the structure of the original model. Bearing that in mind, we propose a new index, that we name the sparsity index, defined as the quotient between the number of zero entries and the number of total elements of the matrix of a system of differential equations. This index provides an indicator of the structure of the models. Original model Non-zero elements ≈N1+N2+N1N2 Total elements (N1+N2)2 Sparsity index ≈N1(N1−1) + N2(N2−1) + N1N2 (N1+N2)2 POD POBD Non-zero elements ≈N2 r≈Nr1+Nr2+Nr1Nr2 Total elements N2 r(Nr1+Nr2)2 Sparsity index ≈0≈Nr1(Nr1−1) + Nr2(Nr2−1) + Nr1Nr2 (Nr1+Nr2)2 Table 3.1: Comparison of the sparsity index for between the original model, and the models arising from the POD and POBD techniques. Table 3.1shows the comparison between the sparsity index of the original model and the ones arising from the POD and POBD techniques for the academic test case (19) . We can easily check that, for this test case, the sparsity index of the POBD systems is much more in concordance with the one of the original model, while the sparsity index of the POD system is always close to zero. 3.2 automatic proper orthogonal block decomposition From Proposition 3.1, we can clearly see that the error bound scales with the number of blocks considered. That gives us an idea that the lesser the number of blocks that we consider, the lesser the error between the solution of the reduced model and the original solution will be. However, we still need to take into account the structure of 66 application to multiple timescales network models three residual behaviors, named D, E and F. We can also identify several boundaries between the behaviors, that we named r,s,t,u,v,w and these changes in behavior can exhibit a transient or not. c c F E D C A r s t u v B w Figure 4.2: Diagram for the different behaviors of model (23) . Zones are represented with capital letters, while the boundary between them are represented with lower case letters. Continuous lines ( u,v,w ) indicate the non-existence of a transient region, while lines with other fashions (r,s,t) indicate the existence of a transient region. Here, we summarize the obtained behaviors of each zone: A. Every cell follows the only one cell behavior. Shown in Figure 4.3a. B. Each cluster follows the two cells antiphase behavior. Shown in Figure 4.3b. C. Each cluster follows the two cells total death behavior. Shown in Figure 4.3c. D. Each cell spikes independent of the other ones. Shown in figures 4.4aand 4.4b. E. This behavior is similar to the total death behavior, but the separation between clusters is no longer clear. Shown in Figure 4.4c. F. The cells inside a cluster spike independent of the other ones, and the cells inside the other cluster remain totally silent. Shown in figures 4.4dand 4.4e. 4.1 intracellular calcium concentration network model 67 (a) Zone A. (cα,cβ) = (1,1). (b) Zone B. (cα,cβ) = (1, −0.25). (c) Zone C. (cα,cβ) = (1, −1). Figure 4.3: Representation of the three of the main behaviors A-C obtained for model (23) . These behaviors represent the behaviors obtained if the intracluster coupling is set to be between 0 and 1, i.e., in-phase synchronization. (a) Zone D. (cα,cβ) = (−0.5, −0.25). (b) Zone D. (cα,cβ) = (−0.5,0.25). (c) Zone E. (cα,cβ) = (−1, −1). (d) Zone F. (cα,cβ) = (−0.1, −0.45). (e) Zone F. (cα,cβ) = (0.01, −0.45). Figure 4.4: Representation of the residual behaviors E-F obtained for model (23) . These behaviors represent the behaviors obtained if the intracluster coupling is set to be lower than 0, i.e., inhibitory synchronization. Now, we present the existing behaviors in the transition between the zones introduced in the previous: u. Cells start to desynchronize, but there are some that remain synchronized. Shown in Figure 4.5a. v. Two clusters. Shown in Figure 4.5b. 68 application to multiple timescales network models w. Relaxation oscillation to total death of the clusters. Shown in Figure 4.5c. r. The cells of a cluster suddenly go from excitation to silence. Shown in Figure 4.5d. s. The cells of a cluster suddenly go from silence to excitation. Shown in Figure 4.5e. t. Cells within each cluster start to desynchronize, the antiphase behavior between cluster starts to vanish. Shown in Figure 4.5f. (a) Boundary u (A-D). (cα,cβ) = (−1, −0.9). (b) Boundary v (A-B). (cα,cβ) = (0.5,0). (c) Boundary w (B-C). (cα,cβ) = (0.5, −0.51). (d) Boundary r (D-E). (cα,cβ) = (−0.7, −0.22). (e) Boundary s (C-E). (cα,cβ) = (−0.6, −1). (f) Boundary t (B-D). (cα,cβ) = (−0.2, −0.32). Figure 4.5: Representation of the evolution of the solution for parameter values in the boundaries of the behaviors obtained for model (23). Heterogeneous case We perform the same study as in the previous section, but we set the heterogeneous parameter randomly in [1,1.5] following a truncated normal distribution N(1.25,0.25) . The results for the Mean Intra Correlation and Mean Inter Correlation are showed in Figure 4.6. 4.1 intracellular calcium concentration network model 69 Mean INTRA 10 20 30 40 50 60 70 80 90 100 cintra 10 20 30 40 50 60 70 80 90 100 cinter -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 (a) Mean Intra Correlation. Mean INTER 10 20 30 40 50 60 70 80 90 100 cintra 10 20 30 40 50 60 70 80 90 100 cinter -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 (b) Mean Inter Correlation. Figure 4.6: Mean Correlations for system (23) for different values of cα and cβ inside the interval [−1,1] . Left panel shows the Mean Intra Correlation, and right panel shows the Mean Inter Correlation. We can see that there are no major changes in the regions of the zones, but that zone C has expanded to the left, reducing the area of behavior E and augmenting the transient r. Zone F also transfers some area to behavior B. 4.1.2Preliminary study In this section, we will build the algebraic version of the original network problem (23) , then, we will present the behaviors that we want to study, and finally, we will determine the a priori number of blocks and their configuration. Algebraic version In order to build the algebraic form of the original network problem (23) , we denote w= (x1, . . . , xN,y1, . . . , yN,z1, . . . , zN) a rearrangement of the vector of unknowns. We note that we have changed the order of the variables in order to help the reader to understand the relationship between them, and that it has no consequences in the following study. We can write problem (23) as ˙w=f(w,γ), with w(0) = w0, (24) where γ= (cα,cβ,k) , the parameter vector, with k= (k1, . . . , kN) the vector of heterogeneity parameters, and cα,cβ the intra-cluster and inter-cluster coupling, respectively. 70 application to multiple timescales network models For the sake of clarity, we define some quantities subsequently: Notation 4.1 K=diag(k)is a diagonal matrix whose entries are the heterogeneity parameters. Notation 4.2 We define the following parameters, Σi=((N−−1)cα+N+cβ,for i =1, . . . , N−, (N+−1)cα+N−cβ,for i =N−+1, . . . , N. Note that, when N is even, Σi=Σ= (N/2 −1)cα+Ncβ/2, ∀i=1, . . . , N. Finally, we can consider the coupling matrix, where top left block is N−×N− , top right is N−×N+, bottom left is N+×N−, and bottom right is N+×N+. M(cα,cβ) =                  Σ1−cα. . . −cα−cβ. . . . . . −cβ −cα....... . .. . .. . . . . .......−cα. . .. . . −cα. . . −cαΣN−−cβ. . . . . . −cβ −cβ. . . . . . −cβΣN−+1−cα. . . −cα . . .. . .−cα....... . . . . .. . .. . .......−cα −cβ. . . . . . −cβ−cα. . . −cαΣN                  We remark here that the matrix M(cα,cβ) admits an affine decomposition depending on cα,cβand Σi. Taking this into account, we can write the algebraic version of system (24) as ˙w=τ((Anp +Ap(k) + B(cα,cβ,k))w+bnp +bp(k) + bnl(w)), (25) where Anp and bnp represent the non-parametric part, Ap(k) and bp(k) the parametric part depending on the heterogeneity parameters, B(cα,cβ) is the parametric part relative to the coupling and bnl(w) denotes the nonlinear part of system (23) . They can be expressed as follows Anp =   4IN−INON ONONON ONON−ε/τzIN   ,Ap(k) = ε   ONONON a0Ka1K ON ONONON   , 4.1 intracellular calcium concentration network model 71 B(cα,cβ) = ε N/2    ONONON KM(cα,cβ)ONON ONONON   , bnp =εzb τz   0N 0N 1N   ,bp(k) = ε   0N a2k 0N   ,bnl(w) =               −x3 1−ϕf(z1) . . . −x3 N−ϕf(zN) 0N εϕr(x1) . . . εϕr(xN)               . With a first look into the model, we can easily see that system (25) depends on many parameters, so the reduction could be unaffordable because of the immense number of parameter vectors that we would have to consider. In order to reduce the number of parameters to consider, we set random heterogeneity parameters ki∈[1,1.5] following a truncated normal distribution N(1.25,0.25) , fixed for all the simulations, that is, every cell is slightly different from each other. Furthermore, we fix cα=1 , that is, we consider that the cells within each cluster are coupled in an excitatory fashion. Therefore, we perform our reduction with respect to the inter-cluster coupling parameter cβ , that is, the coupling between cells of different clusters. This selection is made to be in concordance with previous works (Bandera Moreno et al., 2022), and to avoid chaos and bizarre behaviors, see Figure 4.2. We present in the following section the behaviors in which we will be focusing on. Behaviors to study We focus now on the explanation of the behaviors that will be our object of study for the reduction of problem (25) . We will focus on three main behaviors, for cα=1 , see rightmost border of Figure 4.6, that we explain more in detail in the following. These behaviors differ slightly from the ones presented in Figure 4.2because of the heterogeneity. 1.Almost-in-Phase synchronization (AiP). For cβ∈[0.1,1]. In this range of the parameter, each cell reaches asymptotically an orbit similar to the cycle existing for the relaxation oscillation of the uncoupled case without Mixed Mode Oscillations. All the 3 D oscillators are synchronized with a small 72 application to multiple timescales network models phase-shift or change of frequency depending on its heterogeneity parameter. This behavior is shown in Figure 4.7a. 2.Antiphase Synchronization (AS). For cβ∈[−0.45, −0.05]. In this range of the parameter, each cell within each cluster reach asymptotically the same orbit, and the phase-shift between the clusters is close to half of the period with a small phase-shift or change of frequency depending on its heterogeneity parameter. This behavior is shown in Figure 4.7b. 3.Total Oscillation Death (TOD). For cβ∈[−1, −0.55]. In this range of the parameter, each cell within each cluster reaches asymptotically a stable equilibrium. This behavior is shown in Figure 4.7c. (a) AiP synchronization behavior of system (25) for cβ=1. (b) AS behavior of system (25)for cβ=−0.25. (c) TOD behavior of system (25)for cβ=−1. Figure 4.7: Representation of the different behaviors to study. In blue, the cells belonging to the first cluster and in red, the cells belonging to the second cluster. We only represent x and z , for the sake of clarity. 4.1 intracellular calcium concentration network model 73 Blocks of variables As this test will be the first one performed with the techniques developed in Chapter 3, we need to consider different configurations for them, in order to have a good setup for comparisons with the previous methods. As we have two clusters with three types of variables, we will perform the POBD with 3 and 6blocks, the first POBD will have one block per variable ( x,y,z ), and we will name it POBD 3 , and the second POBD will have one block per cluster variable ( xI,xII,yI,yII,zI,zII ), where subindex I denotes the variables relatives to the first cluster and subindex II denotes the variables relatives to the second cluster. We will name this method POBD6. As APOBD groups the blocks automatically, we start with the maximum number of possible blocks, in our case 6, one per cluster and variable, the same as for POBD6. We have now all the ingredients needed for the study of the techniques developed in Chapter 3for this test case. 4.1.3Numerical Results In this section, we will perform various numerical tests for different number of cells, namely, N=200,500, 1000,2000. In every test case, we will let the original system (25) compute from T=0 to T=50 via a fourth order Runge-Kutta solver with fixed time step of h=10−3 . We perform every method with the data from T=35 to T=50 as we need that the system is almost at an asymptotic behavior. In the construction of reduced systems, the tolerance criteria for the collinearity analysis in the APOBD method is set at 0.2 , that give us two blocks in AiP behavior, one for the fast variables (x) and another one for the slow ones (y,z) , four blocks in AS behavior, (xI) , (xII) , (yI,zI) and (yII,zII) , that is, one block per timescale and cluster, and just one block for TOD behavior, that is, all the variables mixed together. The tolerance criteria of the selection of the reduced dimension in POD methods is set to 10−9 for AiP behavior and 10−6 for AS and TOD behaviors. Nonlinear terms of system (25) are addressed through the Empirical Interpolation Method (EIM, see Section 1.4.1and Barrault et al. (2004) and Grepl et al. (2007)) with a tolerance of 10−6 . In the following, we compare the number of equations, the sparsity of the reduced model via the sparsity index, and the computational time spent to solve the reduced system from T=50 to T=70 , via the speed-ups, defined as the quotient between 74 application to multiple timescales network models the computational time of the original system (25) and the computational time of the reduced system. The objectives of this test case are the following: • First, to validate the new proposed methods. In order to do so, we compare all the methods with N=200 cells in the three behaviors mentioned above. • Then, in order to compare in more detail the methods, we will test them with a more large amount of cells, that is, N=500,1000,2000 for AiP and AS behaviors, as the TOD behavior has not a major interest, as it reproduces the system in a steady state and the reduction is straightforward, as we will see in the following. In Table 4.1, we show the simulation data of the original system (25) for every test case. Test Case N=200 N=500 N=1000 N=2000 Number of equations 600 1500 3000 6000 Sparsity index 0.8867 0.8880 0.8884 0.8887 CPU time ≈6.2 s ≈14 s ≈126 s ≈475 s Table 4.1: Data obtained for system (25) for different test cases depending on the number of cells considered. We show the number of equations, the corresponding sparsity index and the CPU time. Validation of the methods for N =200 The results of the reductions are summarized in Tables 4.2-4.4. In Table 4.2, we see that for this test, POD performs better in terms of speed-up just for AiP behavior, but in terms of sparsity, POBD6provides the best index in the three cases, that is, the best separation of variables. We see that both POBD 3 and APOBD provide a well trade-off between sparsity and speed-up. Method POD POBD3 POBD6APOBD Number of equations 33 26 +12 +10 21 +23 +10 +11 +8+9 26 +16 Sparsity index 0 0.6510 0.7518 0.3685 Speed-up 6.04 4.20 2.70 4.24 Error (·10−3)5.06 4.46 6.95 4.35 Table 4.2: AiP behavior. Comparison of the data for the reduced system for the original system (25) with N=200 . Notation a+...+c indicates the number of equations in each block. In bold fashion, the best results. 4.1 intracellular calcium concentration network model 75 In Table 4.3, we can see the same indexes, but now for AS behavior, POBD6and APOBD are the methods that provide the best speed-ups and sparsity indexes. Method POD POBD3 POBD6APOBD Number of equations 57 53 +14 +8 27 +26 +7+7+4+4 27 +26 +9+8 Sparsity index 0 0.6874 0.7739 0.6937 Speed-up 2.02 1.67 2.20 2.17 Error (·10−3)3.59 7.80 6.24 6.98 Table 4.3: AS behavior. Comparison of the data for the reduced system for the original system (25) with N=200 . Notation a+...+c indicates the number of equations in each block. In bold fashion, the best results. Finally, in Table 4.4, we can see that for TOD behavior, APOBD provides the best speed-up, while POBD6provides the best separation of variables and error. Method POD POBD3 POBD6APOBD Number of equations 1 1 +1+1 1 +1+1+1+1+1 1 Sparsity index 0 0.4444 0.6667 0 Speed-up 34 15.32 9.29 40.06 Error (·10−4)3.09 0.40 0.21 3.09 Table 4.4: TOD behavior. Comparison of the data for the reduced system for the original system (25) with N=200 . Notation a+...+c indicates the number of equations in each block. In bold fashion, the best results. For illustration purposes, in Figure 4.8we show the first 3 basis obtained for each method and each block in AS behavior. Note that the basis for POD (Figure 4.8a) and POBD 3 (Figure 4.8b) are much more erratic, and we cannot identify neither the separation between variables nor clusters, while for POBD6(Figure 4.8c) and APOBD (Figure 4.8d) this separation is clearer, due to the higher sparsity index. This means that, we could not need to go back to the full-dimensional space to analyze the data, allowing us to save more computational time. At this point, we can affirm that this validation test is finished successfully. The reduced methods proposed are able to reproduce the solution of the original problem, while performing a reduction in computational times. It should be noted that while all the methods are comparable in terms of sparsity and error, the novel methods also provide us with a new feature, the perseverance of the structure of the original system. Now, we are able to continue our study. In next sections, we will focus just on the AiP and AS behaviors. As mentioned in the previous, the TOD behavior is 82 application to multiple timescales network models 1.Sustained Oscillation (SO). This behavior was considered in Bonaventura et al. (2022). It is obtained for parameters with classical values a=1 , b=3 , c=1 and d=5 . The separation between timescales is set to ε=10−3 and the external exciting current is set to I=5. For this selection of the parameters, each cell is spiking with no quiescent phase, and all 3 D oscillators are strongly synchronized with a slight phase shift. This behavior is shown in Figure 4.11a. 2.Plateau-like Bursting (PB). This behavior was studied in Barrio and Shilnikov (2011). It is obtained for parameter values a=1 , b=2.52 , c=1 and d=5 . The separation between timescales is set to ε=10−2and the external exciting current is set to I=4. This behavior is characterized by repeated cycles of spiking followed by a prolonged depolarized state known as plateau. All the cells starts their spiking regime at the same time, gradually desynchronize their spikes during the spiking regime, and resynchronize during the quiescence phase. This behavior is shown in Figure 4.11b. (a) SO behavior of system (26)for N =100 cells. (b) PB behavior of system (26)for N =100 cells. Figure 4.11: Representation of the different behaviors to study of system (26) . Left panel shows SO behavior, and right panel shows PB behavior. 4.2 application to bursters models 83 Blocks of variables In this case, we will only consider the POD and APOBD techniques. For the latter, we select one block per variable (x,y,z) for both behaviors, as all the cells belong to the same cluster. Numerical Results In the following study, we will consider various instances for the number of cells N=100,500,1000 . In every test case, we will let the original system (27) compute from T=0 to T=400 via a fourth order Runge-Kutta solver with fixed time step of h=10−2 . We perform both reduction methods with the data from T=200 to T=400 as we need that the system is almost at an asymptotic behavior. This will be enough to capture the main behavior of the system, as shown in Figure 4.11. In the construction of reduced systems, the tolerance criteria for the collinearity analysis in the APOBD method is set at 0.2 , that gives us three blocks for the SO and PB behaviors. The tolerance criteria of the selection of the reduced dimension in POD methods is set to 10−12 for every behavior. Again, nonlinear parts of (27) are addressed by means of the EIM (Barrault et al., 2004; Grepl et al., 2007), now with a tolerance of 10−12. Table 4.7summarizes the simulation data of the original system (27) for both test cases. Test Case N=100 N=500 N=1000 Number of equations 300 1500 3000 Sparsity index 0.8833 0.8878 0.8883 CPU time ≈5.5 s ≈100 s ≈415 s Table 4.7: Simulation data for system (27) for different test cases depending on the number of cells considered. We show the number of equations, the corresponding sparsity index and the CPU time. As done for the previous tests, we compare the number of equations, the sparsity index, and the speed-ups, from T=400 to T=1000 , for the original model and the reduced ones. Furthermore, as an additional way to compare the reduced solutions with the original one, we compute the following a posteriori error indicator, ∆j=max imin k|t(i) j−˜ t(i) k|, (28) 84 application to multiple timescales network models where t(i) j is the j -th time when the original model has a local maximum in the x variable for the cell i , and ˜ t(i) k is the k -th time when the reduced model has a local maximum in the x variable for the cell i . This estimator tells us how shifted are the reduced peaks in comparison to the original ones. Next, we summarize the results obtained for each behavior. Numerical results for SO behavior In this section, we study in detail the SO behavior of system (27) for N=100,500,1000 . Number of cells POD APOBD 100 15 7 +7+3 500 15 7 +7+3 1000 16 7 +7+3 Table 4.8: Comparison of the number of equations of the reduced systems for the original system (27) in SO behavior. The results are summarized in Table 4.8, which shows the number of variables of each cell number case, and figures 4.12a-4.12c, which represent the evolution of sparsity, speed-up and error, respectively. From Table 4.8, we can see that the reduction is stronger when the number of original equations rises, as the number of reduced equations remains almost the same. 0 200 400 600 800 1000 0 0.2 0.4 0.6 0.8 1 Sparsity index SO behaviour POD APOBD (a) Comparison of the sparsity index for each method in SO behavior. 0 200 400 600 800 1000 0 50 100 150 200 250 300 350 Speed-up SO behaviour POD APOBD (b) Comparison of the speed-up for each method in SO behavior. 0 200 400 600 800 1000 10-2 Relative error SO behaviour POD APOBD (c) Comparison of the relative error for each method in SO behavior. Figure 4.12: Comparisons of the results in sparsity index (left), speed-up (center), and relative error (right), for each method in SO behavior of system (27). We can see that although the speed-up is better for the POD method (see Figure 4.12b), APOBD provides an additional order in the relative error, Figure 4.12c. For illustration purposes, in Figure 4.13, we show the first 3 basis obtained for each method and each block in SO behavior. The basis for POD has been represented in 4.2 application to bursters models 85 400 420 440 460 480 500 -100 -80 -60 -40 -20 0 20 40 (a) First basis for the POD method for N=100 cells in sustained oscillation behavior. 400 420 440 460 480 500 -20 -10 0 400 420 440 460 480 500 0 50 100 400 420 440 460 480 500 -40 -20 0 (b) First basis for each block of the APOBD method for N=100 cells in sustained oscillation behavior. Figure 4.13: Comparisons of the basis for each method obtained in the sustained oscillation behavior. In all the cases, the first basis is in blue continuous fashion, the second in red dash-dotted fashion and the third basis in yellow dashed fashion. a) POD. b) APOBD with 3blocks. Figure 4.13a, where we cannot see a clear separation in the timescales. The basis for APOBD has been represented in Figure 4.13b, where this separation is clearer, due to the higher sparsity index. Numerical results for PB behavior In this section, we study in detail the reduction in the PB behavior of system (27) for N=100,500,1000. Number of cells POD APOBD 100 15 7 +6+4 500 15 7 +6+4 1000 15 7 +6+4 Table 4.9: Comparison of the number of equations of the reduced systems for the original system (27) in PB behavior. The results are summarized in Table 4.9, which shows the number of variables of each cell number case, and figures 4.14a-4.14c, which represent the evolution of sparsity, speed-up and error, respectively. From Table 4.9, we can see again that the reduction is stronger when the number of original equations rises, as the number of reduced equations remains almost the same. 86 application to multiple timescales network models It is clear that the APOBD method provides a better speed-up in comparison with the POD for this behavior, maintaining an error of the same order as the one of the original POD method. 0 200 400 600 800 1000 0 0.2 0.4 0.6 0.8 1 Sparsity index PB behaviour POD APOBD (a) Comparison of the sparsity index for each method in PB behavior. 0 200 400 600 800 1000 0 50 100 150 200 Speed-up PB behaviour POD APOBD (b) Comparison of the speed-up for each method in PB behavior. 0 200 400 600 800 1000 10-6 10-5 Relative error PB behaviour POD APOBD (c) Comparison of the relative error for each method in PB behavior. Figure 4.14: Comparisons of the data obtained in sparsity index (left), speed-up (center), and relative error (right), for each method in PB behavior of system (27). For illustration purposes, we show in Figure 4.15 a comparison between the first 3 basis obtained for each method and block in the PB behavior for POD (Figure 4.15a) and APOBD (Figure 4.15b). In the POD case, we cannot observe a separation in the timescales, while in the APOBD case this separation is clearer, due to the higher sparsity index. 400 450 500 550 600 -100 -50 0 50 100 150 200 250 (a) First basis for the POD method for N=100 cells in PB behavior. (b) First basis for each block of the APOBD method for N=100 cells in PB behavior. Figure 4.15: Comparisons of the basis for each method obtained in the PB behavior. In all the cases, the first basis is in blue continuous fashion, the second in red dash-dotted fashion and the third basis in yellow dashed fashion. a) POD. b) APOBD. In Figure 4.16, we show the comparison of the error estimators (28) for both behaviors in the case with N=100 , as it is similar for N=500 and N=1000 . We see that for 4.2 application to bursters models 87 the integrating time, APOBD and POD have a great precision in terms of reproducing the spiking times correctly for the PB behavior as the error is, at most, the integrating time step h , see Figure 4.16b. However, for the SO behavior, as the time increases this precision is lost. This could be due to the transition to a qualitatively different behavior after a stable transient, as it will be discussed in Section 4.3. Nevertheless, according to Figure 4.16a, APOBD performs a more accurate approximation of the spiking times and bursting duration, which is consistent with the relative error shown in Figure 4.12c. 0 20406080100 Spike number (j) 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 j APOBD POD (a) Comparison of the error estimators for SO behavior of system (26)with N =100 cells. 0 102030405060 Spike number (j) 0 0.002 0.004 0.006 0.008 0.01 j APOBD POD (b) Comparison of the errors for PB behavior of system (26)with N =100 cells. Figure 4.16: Comparison of errors of the reduced solutions of system (26) with N=100 , in red fashion, POD, and in blue fashion, APOBD. a) SO. b)PB. It seems that APOBD method performs better when there is a switch between the fast behavior and the slow transient. To test it more in detail, we will compare it with the POD method for another model that exhibits this kind of behavior. 4.2.2Fold/Hopf Bursting Network model In this section, we test the POD and APOBD methods in a case of "Fold/Hopf" (also called "tapered") bursting studied, for instance, in Izhikevich (2000). As reported in Izhikevich (2000), this type of bursting has been found in models of insulin-producing 88 application to multiple timescales network models pancreatic β− cells (Pernarowski, 1994; Smolen et al., 1993) and in models of some enzymatic systems (Holden & Erneux, 1993a,1993b). The model reads:          ˙ x=x−x3/3 −y, ˙ y=ε1(a+x−S(y,z)), ˙ z=ε2kx, where S(y,z) = b/(1+exp((z−y)/d)) , a=1.55 , b=2.5 , d=0.1 , and the timescale separation parameters are set so that the system has three timescales, ε1=0.5 and ε2=0.01 . Note that, then, this model is different from the previous ones in terms of separation timescales as is the only one featuring three timescales 1,ε1and ε2. In this type of burster, the stable equilibrium corresponding to the resting state disappears through a fold bifurcation and the limit cycle attractor corresponding to the spiking state shrinks to a point though a supercritical Hopf bifurcation. Moreover, the fast subsystem undergoes another bifurcation while in the excited state, a fold bifurcation corresponding to the transition from the excited equilibrium to the resting equilibrium, closing the hysteresis loop. The Fold/Hopf Bursting network model can be written as follows, Oi              ˙ xi=xi−x3 i/3 −yi+1 N−1 N ∑ i=1 cij(xj−xi), ˙ yi=ε1(a+xi−S(yi,zi)), ˙ zi=ε2kixi, (29) where i=1, 2, . . . , N . In this test, we consider again a normalized version of the coupling, with cij =1, ∀i=j , and cii =0 . We also add some heterogeneity in the slow equation via ki∈[3.9, 4.1]following a uniform distribution. Preliminary study In this section, we will build the algebraic version of the original network problem (29) , then, we will present the behavior that we want to study, and finally, we will determine the a priori number of blocks and their configuration for the APOBD technique. 4.2 application to bursters models 89 Algebraic version In order to build the algebraic version of model (29) , we begin, as in the previous case, by simplifying the xequations of the model. We can consider D= (dij), with dij =1 N−1 δij N ∑ s=1 csj −cij!, so, the xequations can be rewritten as ˙ xi=xi−x3 i/3 −yi+ N ∑ j=1 dijxj. Now, denoting by γ the parameter vector, the previous simplification allows us to rewrite system (29) as follows w= (Anp +Anp)w+b(γ) + bnl(w;γ), (30) where Anp denotes the non-parametric part, Ap and b the matrix and vector parametric parts, respectively, and bnl stands for the nonlinear part of the system (29) . They can be expressed as follows Ap(γ) =    D ONON ONONON ε2K ONON   ,Anp =   IN−INON ε1INONON ONONON   , b(γ) =    0N ε1a1N 0N   ,bnl(w;γ)    x3/3 1 1+exp (z−y)/d 0N     . Again, we can easily see that model (30) depends on several parameters. In the next epigraph, we will set the behavior that we are trying to model along with its parameters values. Behaviors to study We focus now on the explanation of the behavior that will be our object of study for the reduction of problem (29) . In this case, we will focus on one behavior already described: 90 application to multiple timescales network models 1.Fold/Hopf Bursting (FHB). The resting state disappears through a fold bifurcation, and the spiking state disappears through a supercritical Hopf bifurcation. Moreover, the hysteresis loop is closed through another fold bifurcation. See Figure 4.17. Figure 4.17:FHB behavior of system (30)for N =100 cells. Blocks of variables In this case, we will only consider the POD and APOBD techniques. For the latter, we select one block per variable (x, y, z) for both behaviors, as done in the previous test. Numerical Results In the following study, we will consider various instances for the number of cells N=100,500,1000 . In every test case, we will let the original system (30) compute from T=0 to T=400 via a fourth order Runge-Kutta solver with fixed time step of h=10−2 . We perform both reduction methods with the data from T=200 to T=400 as we need that the system is almost at an asymptotic behavior. This will be enough to capture the main behavior of the system as shown in Figure 4.17. In the construction of reduced systems, the tolerance criteria for the collinearity analysis in the APOBD method is set at 0.2 , that gives us three blocks in the behavior of the study. The tolerance criteria of the selection of the reduced dimension in POD 4.2 application to bursters models 91 methods is set to 10−6 for every behavior, and nonlinear parts of system (30) are addressed by means of the EIM with a tolerance of 10−6. Table 4.10 summarizes the simulation data of the original system (30) for the test case. Test Case N=100 N=500 N=1000 Number of equations 300 1500 3000 Sparsity index 0.8833 0.8867 0.8887 CPU time ≈4 s ≈65 s ≈280 s Table 4.10: Simulation data for system (30) for different test cases depending on the number of cells considered. We show the number of equations, the corresponding sparsity index and the CPU time. Now, we compare the number of equations, the sparsity index, and the speed-ups, from T=400 to T=800 , for the original model and the reduced ones. The results are summarized in Table 4.11, which shows the number of variables of each cell number case, and figures 4.18a-4.18c, which represent the evolution of sparsity, speed-up and error, respectively. Number of cells POD APOBD 100 3 2 +2+2 500 3 2 +2+2 1000 3 2 +2+2 Table 4.11: Comparison of the number of equations of the reduced systems for the original system (30) in FHB behavior. From Table 4.11, we can see again that the reduction is stronger when the number of original equations rises, as the number of reduced equations remains the same. We can again see, that although the speed-up is better for the POD method (see Figure 4.18b), APOBD provides three additional orders in the relative error, Figure 4.18c. So, for almost the same computational time, we obtain a much more precise approximation of the solution of the original problem (30). For illustration purposes, in Figure 4.19, we show the first basis obtained for each method and each block in the studied behavior. The basis for POD are shown in Figure 4.19a, where we cannot see a clear separation in the timescales. The basis for APOBD are shown in Figure 4.19b, where this separation is clearer, thanks to the higher sparsity index. CONCLUSIONS AND OPEN PROBLEMS In this part, we have developed new ROM techniques based on the POD method (Section 1.3.1) for systems that present unknowns with different behaviors, for instance, dynamical systems with multiple timescales. These new techniques have been proven to maintain the structure of the original model, while producing a competitive speedup. Furthermore, we have presented an a priori analysis over the variables of the original model, that could be used in order to group them into blocks. Then, we can apply a data-driven analysis based on the correlation between these group of variables in order to reduce the number of blocks to consider, while maintaining the most part of the structure of the original problem. In the field of applications, we have tested these methods for several biological neural models with multiple timescales that exhibit a rich variety of different qualitative behaviors. These new methods have been proven to produce an accurate approximation of the original solution in much less computational time with respect to solving the original model. The obtained results in speed-up and error are better, or at least, as good as the ones obtained if we build the reduced model with the POD procedure. These new methods also provide a new feature, they keep the structure of the original model. This feature could be useful as we can see the different timescales in the new basis of the reduced model, and could allow the exploitation of this separation of the timescales. Furthermore, we have studied the evolution of the performance of these methods augmenting the number of cells in the original models. The new methods seem to perform better the more cells are considered. We have also noticed that these new techniques seem to perform better when there are more fast variables than slow variables, and where there is a switch between the fast and slow behaviors. Therefore, the proposed techniques based on the POD method and an a priori study of the variables have been proven to be a useful tool for the numerical study of network models. We have shown that they are able to reproduce complex behaviors that these models may exhibit. There are different ways to continue the research presented in this part. We present now some open problems related to them: 99 100 conclusions and open problems • Extension of the presented techniques considering various biological parameters, for example, the coupling and the heterogeneity parameters in ICC model (23). • Use these reduced order models to perform a numerical analysis of the distribution of the behaviors of the models presented with respect to different biological parameters. • Study the reduction of the presented models in the transition between different behaviors. • Application of these techniques to different problems arising from other fields, like physical problems and computational fluid dynamics problems. The research presented in this part of the thesis is an updated version of the research paper Bandera Moreno et al. (2022). And we also foresee the publication of another research paper Bandera Moreno et al. (2024c). Part II REDUCED BASIS MODELLING OF TURBULENCE WITH WELL-DEVELOPED INERTIAL RANGE INTRODUCTION In this part, our emphasis is the study and simulation of turbulent flows, modelled by Partial Differential Equations (PDEs). In particular, we focus on the Reduced Basis (RB) method, with a specific goal of addressing one of the limitations that it presents, namely, the difficulty of obtaining an a posteriori error estimator by means of mathematics-based theories, and its dependence on the numerical discretization. In particular, we study the development of an a posteriori error estimator based upon the Kolmogórov’s energy cascade theory (Kolmogórov, 1941). This part aligns with the Spanish national projects ROAD and HYROM, and the European project ARIA. The research presented in this part of the thesis is an updated version of the results published in the research paper Bandera Moreno et al. (2024a). In order to evaluate the errors between the original solution and the reduced one, the a posteriori error estimates provide a powerful tool to compute estimations of this magnitude in terms of the energy carried by the neglected modes. This technique has been extensively applied to the solution of incompressible flows with success (Busto et al., 2020; Hijazi et al., 2023; Kunisch & Volkwein, 2002; Negri et al., 2015; Stabile & Rozza, 2018). Yet, a joint Proper Orthogonal Decomposition (POD)-Galerkin procedure encounters a significant challenge. To achieve errors below a prescribed level, computing the solution of the original model for multiple parameter values becomes imperative. However, in scenarios concerning incompressible flows, this computation might quickly escalate to an unmanageable scale. Hence, a procedure was developed in order to obtain errors below a certain threshold while minimizing the number of needed solutions of the original model. This issue is addressed by the certified RB method, in which the snapshots are recursively computed by a Greedy Algorithm to progressively reduce the errors. This procedure needs a fast a posteriori error estimation (Grepl et al., 2007; Hesthaven et al., 2016; Quarteroni et al., 2016). It has been applied to incompressible flows with increasing complexity: Stokes, Navier-Stokes and multi-parametric flows (Ali et al., 2020; Ballarin et al., 2015; Deparis, 2008; Deparis & Rozza, 2009; Manzoni, 2014). Moreover, some error indicators based on the dual norm of the discrete time-averaged 103 104 introduction residual have been used, for example in Fick et al. (2018) and in Tsai and Fischer (2022), for the parameter selection in the Greedy Algorithm. Kolmogorov’s energy cascade theory (also known as K 41 theory) is a fundamental concept in the field of fluid mechanics (Chacón Rebollo & Lewandowski, 2014; Kolmogórov, 1941; Richardson, 1922), and it provides a theoretical framework for understanding the energy transfer between different scales in a turbulent flow. In a fully developed turbulent flow, the kinetic energy is transferred from the large-scale eddies to smaller and smaller scales through a series of nonlinear interactions until it reaches the smallest scales, where it is dissipated into heat. This process is known as the energy cascade, and it is characterized by a self-similar scaling behavior in the inertial subrange of scales. Kolmogorov’s theory provides a statistical description of turbulence that has been widely used to guide the development of turbulence models and numerical simulations of turbulent flows. The energy cascade concept has lead to the use of the high-order modes of the POD expansion (see Section 1.3.1) to model the small scales within the inertial range. However, a very large number of modes is needed to achieve a good precision, making the Reduced Order Model (ROM) need computing times similar to a direct numerical solver (Hijazi et al., 2020). An alternative strategy is to build ROMs for well-established turbulence models, actually here we address Large Eddy Simulation (LES) models. In this approach no turbulence modelling procedure by a ROM is intended, instead the purpose is to build surrogate LES models that can be run in highly reduced computing times. Certified RB LES Smagorinsky models for steady flows have been built with high computing speed-ups for single and multi-parametric configurations (Ballarin et al., 2020; Caravaca García, 2022; Chacón Rebollo et al., 2023; Chacón Rebollo et al., 2017). The error certification is based upon a posteriori error estimation procedures, that in their turn are built by the Brezzi-Rappaz-Raviart theory of approximation of regular branches of non-linear PDEs (Brezzi et al., 1980), thanks to the enhanced regularity of the solutions of the Smagorinsky model. However, this procedure faces severe limitations. Indeed, the building of a posteriori error indicators is closely linked to the actual turbulence model and the actual FOM used, as well as to the nature of the parameters (either physical or geometrical) to deal with in the ROM. Also, the extension of this error estimation procedure to transient flows requires very high computational costs, thus limiting the advantage of reduced modeling. introduction 105 In this part, we tackle the construction of a RB method for general transient LES parametric turbulence models, regardless of the FOM used. The basic tool is an error indicator based upon the physics of turbulent flows, rather than on the mathematics used to analyze the models. This error indicator measures the deviation of the energy spectrum of a trial solution versus the theoretical k−5/3 spectrum predicted by the K 41 theory of turbulence in statistical equilibrium. In practice, this indicator leads to a stagnation of the Greedy Algorithm, that selects again a parameter already chosen. When stagnation occurs, we use a correction of this indicator, now measuring the deviation of the spectrum of the trial solution versus the spectrum of the reduced solution obtained in the current reduced basis. Indeed, the ROM solution is intended to approximate the FOM solution, whose spectrum already may have some error with respect to the theoretical K41 spectrum. To validate this indicator, we address some academic tests for the Smagorinsky turbulence model. We consider 2 D periodic flows, yet meaningful, as they are designed to present the inertial spectrum predicted by the K 41 theory. We assume that the Smagorinsky flow already has a part of its spectrum within the inertial range. The Empirical Interpolation Method (EIM) is used to compute reduced approximations of the non-linear eddy viscosity term, see Section 1.4.1. We use various different POD filtering strategies to reduce the time snapshots. The Fast Fourier Transform (FFT) is used to compute the spectrum of the trial solutions. We compare the errors obtained with this strategy (using the error indicator) versus those obtained with the best error indicator, that is, the exact error. As a result, we obtain errors between the ROM and the FOM solutions quite close to the optimal ones. Further, we obtain a spectral decay of the error as the dimension of the reduced basis space increases, quite close to the one obtained if the exact error is used as indicator. We obtain speed-up rates of computing times of nearly 23. This part is structured as follows. In Chapter 5, we introduce error indicator based upon the K 41 theory and state the combined Greedy Algorithm-POD strategy to build the reduced basis from this indicator. We also present in this chapter some error bounds estimates for Kolmogórov’s estimator. The numerical tests are presented in Chapter 6. In particular, we address some academic test for the Smagorinsky turbulence model. 5 REDUCED BASIS MODELLING OF TURBULENT FLOWS BASED ON THE STATISTICAL TURBULENCE THEORY The main purpose behind this chapter is the introduction of an a posteriori error indicator for turbulent flows. This a posteriori error indicator will allow us to select the parameter in the Greedy procedure for the construction of the reduced basis model. For the development of this indicator, we have used general results from the statistical turbulence theory, attributed to Kolmogórov (Kolmogórov, 1941), also known as K 41 theory. This theory presents an expression for the energy cascade (Chacón Rebollo & Lewandowski, 2014; Richardson, 1922) and the main assumption behind this indicator is that the difference between the energy spectrum of a certain solution and the theoretical spectrum predicted by this theory is small for accurate solutions. The structure of the Chapter is as follows. In Section 5.1, we provide error estimates based on the differences of the energy spectrum of two velocity fields, and we introduce an error indicator based upon the K 41 theory. On the other hand, in Section 5.2, we state the combined Greedy Algorithm-POD strategy to build the reduced basis from this indicator. 5.1 a posteriori error indicator based on the statistical turbulence theory Andrei Kolmogórov stated in Kolmogórov (1941) that under suitable similarity and isotropy assumptions for turbulent flow in statistical equilibrium, there exists an inertial range [k1,k2] , in which, the physical magnitude, called energy spectrum E(k) , admits an expression in terms of the turbulent dissipation ε , that is, an indicator of the energy dissipated through viscosity, and the wavenumber k. 107 114 rb modelling of turbulent flows Proof : We start from equation (15) , taking qh=ph , we add and subtract (∂twN,vh)Ω+ ν(∇wN,∇vh). Naming ρN=wN−uh, we obtain the following (∂twN,vh)Ω+a(wN,vh;µ) + aS(uh;uh,vh;µ) + c(uh,uh,vh;µ) =⟨f,vh⟩+ (∂tρN,vh)Ω+a(ρN,vh;µ),∀vh∈Xh. Now, we consider vh=vN∈XN , and we subtract equation (18) considering qN=pN . If we name eN=wN−uN, then, we obtain (∂teN,vN)Ω+a(eN,vN;µ) + (aS(uh;uh,vN;µ)−aS(uN;uN,vN;µ)) +(c(uh,uh,vN;µ)−c(uN,uN,vN;µ)) = (∂tρN,vN)Ω+a(ρN,vN;µ). Considering vN=eN and ν=1/µ , and following the same ideas as in the first part of the proof of Theorem 5.2(see computation of equation (38)), we get to 1 2∂t∥eN∥2+ν∥∇eN∥2 2+aS(uh;eN,eN;µ) = (∂tρN,eN)Ω+a(ρN,eN;µ) + aS(uh;ρN,eN;µ) +c(uh,ρN,eN;µ) + c(ρN,wN,eN;µ)−c(eN,wN,eN;µ) +(aS(uh;uh,eN;µ)−aS(uN;uh,eN;µ)). As done for equation (38) in the proof of Theorem 5.2, we can use the effective viscosity in order to bound the last two terms on the left-hand side. After taking absolute value on the right-hand side and applying the triangular inequality, we obtain the following bound 1 2∂t∥eN∥2+νe∥∇eN∥2 2≤ |(∂tρN,eN)Ω+a(ρN,eN;µ) + aS(uh;ρN,eN;µ)| +|c(uh,ρN,eN;µ) + c(ρN,wN,eN;µ)−c(eN,wN,eN;µ)| +|aS(uh;uh,eN;µ)−aS(uN;uh,eN;µ)|. (42) 5.1 a posteriori error indicator based on statistical turbulence theory 115 The first term on the right-hand side can be upper bounded using Hölder’s and Young’s inequalities as follows |(∂tρN,eN)Ω+a(ρN,eN;µ) + aS(uh;ρN,eN;µ)| ≤ ∥∂tρN∥2∥eN∥2+ν∥∇ρN∥2∥∇eN∥2+∥νT(uh)∇ρN∥2∥∇eN∥2 ≤C∥∂tρN∥2 2+C∥eN∥2 2+ν 2ϵ∗∥∇ρN∥2 2+νϵ∗ 2∥∇eN∥2 2+C 2ϵ∥νT(uh)∇ρN∥2 2+Cϵ 2∥∇eN∥2 2. (43) For the second term on the right-hand side, we need to take into account expression (17) , and Hölder’s and Young’s inequalities and the Sobolev’s embedding L6⊂W1,2 =H1 , as we did on equation (40), to obtain |c(uh,ρN,eN;µ) + c(ρN,wN,eN;µ)−c(eN,wN,eN;µ)| ≤ ∥∇uh∥2∥∇ρN∥2∥∇eN∥2+∥∇wN∥2∥∇ρN∥2∥∇eN∥2 +∥∇wN∥3∥∇eN∥6∥∇eN∥2 2+∥wN∥∞∥eN∥2∥∇eN∥2 2 ≤C 2ϵ∥∇uh∥2 2∥∇ρN∥2 2+Cϵ 2∥∇eN∥2 2+C 2ϵ∥∇wN∥2 2∥∇ρN∥2 2+Cϵ 2∥∇eN∥2 2 +C 2ϵ∥∇wN∥2 3+∥wN∥2 ∞∥eN∥2 2+Cϵ 2∥∇eN∥2 2. (44) In order to bound the last term on the right-hand side, we simply need to apply Hölder’s and Young’s inequalities to get to |(aS(uh;uh,eN;µ)−aS(uN;uh,eN;µ))| ≤ C 2ϵ∥νT(uh)−νT(uN)∥2 ∞∥∇uh∥2 2+Cϵ 2∥∇eN∥2 2. (45) Applying bounds (43) - (45) in equation (42) , and arranging the terms we get to the following bound 1 2∂t∥eN∥2+νe−νϵ∗ 2−5Cϵ 2∥∇eN∥2 2≤C1+C 2ϵ∥∇wN∥2 3+∥wN∥2 ∞∥eN∥2 2 +C∥∂tρN∥2 2+ν 2ϵ∗+C 2ϵ∥∇uh∥2 2+∥∇wN∥2 2∥∇ρN∥2 2+C 2ϵ∥νT(uh)∇ρN∥2 2 +C 2ϵ∥νT(uh)−νT(uN)∥2 ∞∥∇uh∥2 2. 116 rb modelling of turbulent flows Selecting ϵ∗=νe/(2ν) and ϵ=νe/(10C) in the previous equation, and multiplying both terms by 2, we obtain ∂t∥eN∥2+νe∥∇eN∥2 2≤C1+∥∇wN∥2 3+∥wN∥2 ∞ νe∥eN∥2 2 +C∥∂tρN∥2 2+C νeν2+∥∇uh∥2 2+∥∇wN∥2 2∥∇ρN∥2 2+C νe ∥νT(uh)∇ρN∥2 2 +C νe (δh)∥∇uh∥2 2, where we have used Proposition 5.1in the last term, and all the constants are included in C. Applying Grönwall’s inequality, we get to desired estimate. This supports the use of an energy spectrum-based error estimator in the building of a reduced basis model for the Smagorinsky model. 5.1.3Kolmogórov’s a posteriori error indicator In the section 5.1.1, we have seen that the energy spectrum admits a separate expression in terms of the wavenumber k . In this section, we will exploit that property in order to build an a posteriori error indicator for turbulent flows. Let u(x) be a velocity field on a certain domain Ω∈Rd , with d=2,3 . Thanks to the Fourier transform (Fourier, 1888; Stein & Weiss, 1971), we are able to decompose the velocity field in its Fourier modes ˆ u(k) = ZΩu(x)ϕk(x)dx, where ϕk=exp (−2πik·x) . Furthermore, we can recover the velocity field from its Fourier modes thanks to the inverse Fourier transform, u(x) = ZRdˆ u(k)ϕk(x)dk, with ϕkthe complex conjugate of ϕk. 5.1 a posteriori error indicator based on statistical turbulence theory 117 Now, assuming that the velocity field is ergodic and regular enough, that is u(x)∈ L2(Ω) , we can apply Parseval’s identity to obtain an expression of the kinetic energy with respect to its Fourier modes Ekin =ZRd|ˆ u(k)|2dk. Therefore, in the applications we can approximate the total kinetic energy by Ekin ≈ L ∑ k=0 E(k) = L ∑ k=0 ∑ k≤|k|≪k+1 |ˆ u(k)|2dk, where L is related to the maximum wavenumber that can be solved numerically, and ˆ u(k)can be computed by means of the discrete Fourier transform. For a high fidelity solution, we can assume that, in the inertial subrange, this physical property accurately follows the expression E(t,k;γ) = α(t;γ)k−5/3, (46) where α(t;γ)>0 is a positive function depending on time and parameters γ . One of those parameters can be the turbulent dissipation ε of the solution for the Reynolds number, Re . For the sake of clarity, we will suppose that the only parameter acting on α is the Reynolds number, that is, γ=Re , so we can write the energy spectrum as EN(t,k;Re). It is fair to assume that the spectrum of the solution of the original (full-order) model closely follows (46) , as the full-order model is assumed to approximate the continuous problem faithfully in the resolved part of the inertial spectrum. Let us denote EN(t,k;Re) be the energy spectrum associated to the reduced solution that belongs to a reduced space of dimension N for a certain Re and at a certain time t . We define the a posteriori error indicator as follows, ∆N(t;Re) = min α∈RZk2 k1 |EN(t,k;Re)−α(t;Re)k−5/3|dk1/2 . (47) This quantity measures the difference between the spectrum of a reduced solution at a given time and the theoretical spectrum in the inertial range [k1,k2] described by the theory. The main advantage of this error indicator is that it can be applied no matter the numerical discretization and the turbulence model we are working with, as 118 rb modelling of turbulent flows it only involves the spectrum of the trial solution. That means that thanks to the fact that this indicator is build using only physical properties will allow us to avoid the construction of ad hoc a posteriori error estimators for the Reduced Basis discretization, which can be technically demanding, as the ones developed in Chacón Rebollo et al. (2023), Chacón Rebollo et al. (2017), and Deparis (2008), Manzoni (2014), based upon the Brezzi-Rappaz-Raviart theory (Brezzi et al., 1980) and that necessarily are specific for each actual turbulence model and numerical discretization. In practice, we substitute k2 by kc<k2 , where kc is the smallest scale we can solve numerically by means of the turbulence model considered. In most cases, kc is related to the mesh size δ . Since this indicator is valid whenever the turbulence model solves a range of the inertial spectrum, special care should be taken in the construction of the mesh Th , so we are effectively solving in a suitable range, that is, kc∈(k1,k2) , in practice kc≪k2. 5.2 general algorithm to build the reduced basis spaces In this section, we describe a general algorithm based on the a posteriori error indicator ∆N(Re) , presented in (47) , for the construction of the reduced basis model of LES turbulent model. We describe it here in a general framework, as it is independent of the actual LES model considered. We intend to solve parametric turbulent flows depending on a parameter vector γ∈ D ⊂ RP , for any integer P≥1 . In the applications that follow this parameter will be the Reynolds number, although it can also be a different physical or geometrical parameter, o several of these. We shall suppose that γ ranges on a compact set D ⊂ RP . Given a partition, Dtrain, of the parameter space, D, and a given starting parameter vector γ0 , the general algorithm for the Reduced Basis functions selection can be summarized as follows: 1. In iteration i≥1 , solve the original full-order turbulence model for γ∗ i−1∈ Dtrain in some temporal interval, so the turbulence has reached the statistical equilibrium. Store the results of the velocity and pressure, separately. 2. Apply some technique to compress the snapshots for the velocity and pressure (for instance, one of the POD techniques explained in Section 1.3.3), and add them to the previous velocity and pressure Reduced Basis spaces. This will give us the new reduced dimension Ni. 5.2 general algorithm to build the reduced basis spaces 119 3. In order to assure that the inf-sup condition (14) is satisfied for the pair of reduced velocity and pressure spaces, compute the inner pressure supremizer for the basis of the pressure resulting from the previous step (Ballarin et al., 2015), and append it to the basis of the reduced velocity space (see Subsection 2.2.2for a detailed description). If some stabilization techniques are considered (Ali et al., 2020; Chacón Rebollo et al., 2023) this step would be no longer needed to be taken into account. 4. Compute the a posteriori error indicator ∆(i) Ni(γ) for each γ∈ Dtrain associated to the reduced basis of dimension Ni . The error indicator should be computed for some discrete time sets to be determined. 5. Select the new parameter γ∗ i according to some criteria based upon the results of the error indicator ∆(i) Ni. We will focus on this point in Section 5.2.1. 6. Check the stopping criteria. If the criteria are fulfilled, then stop the algorithm, otherwise, repeat from step 1. As said before, this procedure does not depend on the way we obtain the full-order solution, so it is suitable for any turbulence model that satisfies the Kolmogorov’s law and allows us to avoid the construction of a posteriori error estimators specific for each problem, that could involve some technical difficulties. 5.2.1Parameter selection Taking into account that we build the reduced solutions from the numerical approximations of the original solutions, the indicator ∆N(t;γ) is not expected to vanish as N grows. But, it is expected to converge to the indicator for the solution of the full order model ∆h(t;γ), with ∆h(t;γ) = min α∈RZkc k1 |Eh(k;γ)−α(γ)k−5/3|2dk1/2 where Eh(t,k;γ) denotes the spectrum of the full-order solution at time t for k∈(k1,kc) , that is, the resolved inertial range. This spectrum should be close to the theoretical k−5/3 law if the full-order model is accurate, however, some error is to be expected as the full-order model is not exact and real-life flows do not follow strictly the theoretical law. 120 rb modelling of turbulent flows Then, the usual criterion for the selection of the best parameter in the Greedy algorithm, that is γ∗=maxγ∈D ∆N(γ) , is no longer useful. Instead, we propose in Algorithm 1new criteria for the selection of the next parameter at each step. Algorithm 1:Parameter selection criteria Require: Si, the set of previously selected parameters. ∆(i−1) Ni−1and ∆(i) Ni, the reduced energy spectrum of the previous and current iteration, respectively. Goal: Obtain the next parameter γ∗for the Greedy algorithm. Compute γ∗ i=argmaxγ∈Dtrain ∆(i) Ni(γ); if γ∗ i∈Sithen γ∗ i=argmaxγ∈Dtrain\Si|∆(i) Ni(γ)−∆(i−1) Ni−1(γ)| end Algorithm 1first aims to add to the basis the information of the solutions whose energy spectrum is farther away from the theoretical one, then, if the parameter selected that way has already been selected in a previous iteration, we select between the non-selected parameters, the one that provides the largest difference between the energy spectrums of the previous and the current reduced solutions. In the latter case, we are assuming that the reduced error indicator ∆N(γ) is getting close to ∆h(γ) , and the previous solution plays the role of the full-order one to compare with it. 6 APPLICATION TO SMAGORINSKY TURBULENCE MODEL In this chapter, we aim to validate the behavior of the indicator ∆N , equation (47) , developed in the previous chapter. In particular, we address some academic tests for the Smagorinsky turbulence model. We consider 2 D periodic flows, yet meaningful, as they are designed to present the inertial spectrum predicted by Kolmogórov’s theory. We assume that the Smagorinsky flow already has a part of its spectrum within the inertial range. The structure of the chapter is the following. First, we will perform a preliminary study in Section 6.1. In this section, we will focus on the presentation of the space-time discretization and the data of the problem. Then, we will select the POD+Greedy procedure and study the treatment of the eddy viscosity term and the time variability in the error indicator (47) . Section 6.2will be devoted to the presentation of the validation tests and numerical results. 6.1 refinement of the characteristics of the method This section is aimed to performing a preliminary study on the application of the Kolmogórov indicator (47) to the Smagorinsky turbulence model. In particular, we will tune up the method characteristics and play close attention to the selection of the POD+Greedy procedure. 6.1.1Space-time discretization We consider the domain Ω×[0, Tf] , with Ω= [−1/2,1/2]2 and Tf=15 , and we solve the Smagorinsky model (15) for 2 D periodic flows, with periodic boundary conditions and assume no external sources, thus f=0. 121 122 application to smagorinsky turbulence model In order to obtain the approximate solution, we use a semi-discretization in time using the Crank-Nicolson method with dt =0.3125 , so the number of time steps is L=48 , and for the resolution of the resulting problem we use stable finite elements for the pressure velocity, specifically Taylor-Hood i.e., P2−P1 . The mesh on which it is solved is a structured mesh, with m=64 intervals in each of the directions, obtaining a mesh with 8192 triangles and 4225 vertices. 6.1.2Data of the problem and initial Condition The parameter considered is the Reynolds number γ , which we consider moves in the range D= [1000, 16000] . For this range, approximately a time greater than 10 is needed to reach the profile of the Kolmogorov spectrum. In this case there is no inertial spectrum k−3 of the 2 D turbulence, as a consequence of no forcing at large wavenumbers (Kraichnan, 1967). For the initial condition, we need to consider a velocity field whose energy spectrum follows the form (46) . In order to build this velocity field, we will start from its Fourier transform, so we can impose the desired energy spectrum. In particular, we take w0 h= (v,v)where vis such that its Fourier transform is given by ˆ v(k) = (k−(5/3+1)/2 if 0 <k≤m/4, 0 other case. However, this velocity field is completely artificial, so it may not correspond to a realistic physical setting. In order to overcome this issue, we will define our initial condition u0 h as the solution of the Smagorinsky model (15) , with stating point w0 h , for the mean value of the parameter, γ=8500 , at time Tf=15 . Thanks to this, the resulting initial condition u0 h has an energy spectrum in statistical equilibrium, as required by our theory. In Figure 6.1, we can see the initial condition u0 h constructed and its energy spectrum. We can observe that for wavenumbers k between 5 and 32 a good approximation of the initial spectrum is obtained. For values larger than kc=32 , we see a sharp decay of the energy, this is due to the fact that these wavenumbers are outside the circle of the largest radius completely contained in the unit square. 6.1 refinement of the characteristics of the method 123 100101102 k 10-12 10-11 10-10 10-9 10-8 10-7 10-6 10-5 10-4 10-3 10-2 E(k) k-5/3 Figure 6.1: Velocity field of the initial condition u0 h , (left panel) and the corresponding energy spectrum (right panel). 6.1.3Selection of the POD+Greedy procedure In this subsection, we aim to select the method that we will be using in the following numerical study. More precisely, we are comparing the four approaches we presented in Section 1.3.3. For this selection, and in the following numerical results, we will use the following combined velocity-pressure relative error ϵN(t;γ) = ∥uh(t;γ)−uN(t;γ)∥2 H1+∥ph(t;γ)−pN(t;γ)∥2 2 ∥uh(t;γ)∥2 H1+∥ph(t;γ)∥2 2!1/2 . From this definition, we focus on two quantities of interest: the error at the final time step, i.e. ϵN(Tf;γ), and the time-integrated error, i.e. ZTf t0 ϵN(t;γ)dt. In Figure 6.2, we can see the comparison between the methods presented in Section 1.3.3in terms of the l2 and l∞ norms over the parameter, for the same tolerance of the POD procedures. The behavior is similar in both norms, and we can conclude various things from this comparison. On the one hand, we can see that the Reduced POD approach outperforms the others but stagnates in a low dimension of the reduced basis with no further improvement. On the other hand, the complete POD provides a better approximation of the full-order solution for a higher dimension of the reduced basis, achieving almost two more orders of magnitude in the relative error with respect to the Hierarchical POD (Haasdonk, 2017) and Incremental POD (Himpe et al., 2018). For the sake of completeness, we present in Table 6.1a comparison of the evolution of the 130 application to smagorinsky turbulence model It. γN maxγϵN(γ) 1 1000 30 5.54 ·10−1 2 16000 70 6.74 ·10−2 3 2875 96 1.11 ·10−2 4 6625 116 4.43 ·10−3 5 1625 131 1.03 ·10−3 6 11000 143 3.01 ·10−4 7 4125 148 1.28 ·10−4 8 2250 149 8.55 ·10−5 Table 6.2: Convergence story of the POD+Greedy algorithm with the Complete POD approach, using the actual relative error, ϵN(γ), for the parameter selection. It. γN maxγ∆N(γ)maxγ|∆N(γ)−∆N∗(γ)|maxγϵN(γ) 1 1000 30 7.92 ·10−1- 5.54 ·10−1 2 16000 70 5.13 ·10−1- 6.74 ·10−2 3 1625 91 3.53 ·10−12.11 ·10−12.23 ·10−2 4 2250 107 " 4.32 ·10−26.76 ·10−3 5 2875 117 " 1.68 ·10−23.65 ·10−3 6 6000 134 " 1.72 ·10−38.41 ·10−4 7 12875 146 " 1.46 ·10−31.26 ·10−4 8 11625 149 " 5.93 ·10−58.70 ·10−5 Table 6.3: Convergence story of the POD+Greedy algorithm with the Complete POD approach, using Kolmogórov’s indicator, ∆N(γ) , for the parameter selection. Here, ∆N∗(γ) denotes the error indicator in the previous iteration. As we mentioned in Section 5.2.1, the indicator ∆N(γ) is not expected to converge to zero when the dimension Nof the reduced space grows, rather it should converge to ∆h , since the reduced solution is constructed from the finite element approximation. This convergence is shown in the Figure 6.8, where we can easily see that the larger value for the indicator is indeed achieved for γ=16000 , and that after the third iteration, the rest of the values of the indicator are below that threshold. This fact supports the use of Algorithm 1for the selection of the parameter. From Figure 6.8, it can also be seen easily that after the fifth iteration, both, reduced and finite element indicators are very similar. In order to provide a better insight into this matter, we show in Figure 6.9, the absolute difference between the reduced indicators and the finite elements one. After the last iteration, we can see that the major difference between them is of order 10−4. In Table 6.4, we show for several values of the parameter, which some do not correspond to those chosen by the Greedy algorithm, the calculation times of the finite 6.2 numerical results 131 2000 4000 6000 8000 10000 12000 14000 16000 Re 10-4 10-2 100 N(Re) Relative Error Evolution N=30 N=70 N=91 N=107 N=117 N=134 N=146 N=149 Figure 6.7: Relative error ϵN(γ) in the range of the parameter, for each iteration of the POD+Greedy algorithm with the Complete POD procedure. Here, we are using Kolmogórov’s indicator (47)for the parameter selection. 2000 4000 6000 8000 10000 12000 14000 16000 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 N=30 N=70 N=91 N=107 N=117 N=134 N=145 N=149 EF Figure 6.8: Convergence of Kolmogórov’s indicator (47) for the reduced solutions to the finite elements Kolmogórov’s indicator. Indicator from reduced solutions are shown in continuous fashion, while the one from the finite elements methods is shown in dash-dotted black fashion. element method TFE and of the reduced method TRB , as well as the corresponding 132 application to smagorinsky turbulence model 0 5000 10000 15000 10-6 10-4 10-2 100 N=30 N=70 N=91 N=107 N=117 N=134 N=145 N=149 Figure 6.9: Absolute difference between the Kolmogórov’s indicator (47) for the reduced solutions and the finite elements Kolmogórov’s indicator. speed-up. The lower part shows the values corresponding to the indicators ∆N(γ) , ∆h(γ)and the actual relative error ϵN(γ)committed for the different cases. γ3000 6000 9000 12000 15000 TFE(s)62.982 62.579 62.634 62.632 62.535 TRB(s)2.73086 2.74046 2.71848 2.72783 2.71082 Speed-up 23.06 22.84 23.04 22.96 23.07 ∆h(γ)0.327244 0.318601 0.311539 0.309520 0.337642 ∆N(γ)0.326974 0.318606 0.311492 0.309526 0.337656 ϵN(γ)6.27 ·10−52.05 ·10−54.64 ·10−51.28 ·10−51.78 ·10−5 Table 6.4: Validation of the reduced basis model obtained by means of a POD+Greedy approach using Kolmogórov’s indicator (47) for the parameter selection. All the presented data has been obtained from the eight iteration of the POD+Greedy algorithm, with a reduced dimension of N=149 . The speed-up obtained is close to 23 in each of the cases tested, which is satisfying for an evolution turbulence model. The values of the indicators ∆N(γ) , ∆h(γ) , are very close to each other, having at least two coincident decimal places. We can also see that the relative errors for the selected parameters are below 10−4. In Figure 6.10, we show the comparison between the FE (left) and RB (center) velocity field, with the absolute error (right) between them, for a Reynolds value γ=15000 , 6.2 numerical results 133 and times t=20,25,30 seconds (from top to bottom), where we can observe that the difference between both FE and RB velocity field is practically negligible, since the maximum of the absolute error between both fields is less that 10−7. Moreover, in Figure 6.11 we show the temporal evolution of the velocity magnitude for some spatial locations, in (0,0) , (0.25,0.25) and (0.5,0.5) for both the high FE solution and the RB solution for a Reynolds parameter value of γ=15000 . We can observe the regularity of the flow for both the FE and RB solution. Therefore, we can conclude that the error indicator has been successfully extended from parametric linear parabolic problems to LES parametric turbulence models, obtaining relative errors in the order of 10−5. 6.2.2Comparison against equispaced sampling The objective of this second test is to compare the results obtained in the previous test with the ones obtained when considering equispaced parameter values, instead of considering a Greedy algorithm for selecting them. It. N maxγϵN(γ)l2(ϵN(γ)) 1 30 5.54 ·10−15.53 ·10−1 2 70 6.74 ·10−21.95 ·10−1 3 96 2.63 ·10−24.95 ·10−2 4 116 1.27 ·10−21.96 ·10−2 5 126 6.96 ·10−39.18 ·10−3 6 133 3.81 ·10−34.49 ·10−3 7 134 2.85 ·10−33.20 ·10−3 8 138 1.96 ·10−32.12 ·10−3 Table 6.5: Convergence of the algorithm with the Complete POD approach, using equispaced parameter values. Here, we present the maximum relative error and their l2-norm. In Table 6.5, we present the maximum relative error and the l2 -norm of the relative error in the parameter range. We are using the Complete POD approach, presented in Section 1.3.3, for the construction of the reduced basis. We can see that, after the second iteration, these values differ from the ones presented for the use of the actual relative error (Table 6.2) and those for the use of Kolmogórov’s indicator (Table 6.3). While for the use of equispaced parameters after eight iterations the maximum error is above 10−3 , the ones for the use of the actual error and Kolmogórov’s error, after eight iterations the maximum error below 10−4. 134 application to smagorinsky turbulence model Figure 6.10: Comparison between FE (left) and RB (right) velocity field, along with the absolute error between them (right), for γ=15000, at t =20s,25s, 30s from top to bottom. 6.2 numerical results 135 Figure 6.11: Temporal evolution of the FE (left) and RB (right) velocity magnitudes at three different points, x= (0,0),x= (0.25,0.25)and x= (0.5,0.5). We also present in Figure 6.12 the evolution of the relative error in the range of the parameter. We can see that the largest errors appear for lower values of the parameter, decreasing slowly as the reduced basis dimension, N , increases. This is due to the fact that this approach to enrich the reduced basis is blind to the actual behavior of the solution. The approach involving Kolmogórov’s indicator (47) has been proven to identify this kind of behavior, and, therefore, improving the enrichment of the reduced basis. In order to have a better visual of the differences between each approach, we will compare their behaviors in two norms, the l∞ -norm, and l2 -norm. We present, in Figure 6.13a, the comparison of the decay of the l2-norm of the relative error over the parameter range and, in Figure 6.13bthe comparison of the decay of the maximum relative error over the parameter range. We observe a spectral convergence with respect to the dimension of the reduced space N. We also observe that the a posteriori error indicator based on Kolmogorov’s law (47) , provides very similar results compared to the ones obtained using the exact error, in both l2 and maximum relative error decay, when using the Complete POD procedure presented in Section 1.3.3. Also, we observe that it provides better results than just selecting equispaced parameters in the reduced basis construction, specially in the maximum relative error. This is due to the fact that for equispaced sampling the errors are localized in the lower Reynolds numbers, whereas thanks to our approach the information about these parameters are included in the construction of the reduced 136 application to smagorinsky turbulence model 2000 4000 6000 8000 10000 12000 14000 16000 Re 10-4 10-3 10-2 10-1 100 N(Re) Relative error evolution N=30 N=70 N=96 N=116 N=126 N=133 N=134 N=138 Figure 6.12: Relative error ϵN(γ) in the range of the parameter, for each iteration of the POD+Greedy algorithm with Complete POD procedure. Here, we are selecting equispaced parameter values. 20 40 60 80 100 120 140 160 N 10-4 10-3 10-2 10-1 100 101 l2(N(Re)) Comparison error decay Indicator Error Equispaced (a) Comparison of relative l2-norm decay. 20 40 60 80 100 120 140 160 N 10-5 10-4 10-3 10-2 10-1 100 max( N(Re)) Comparison error decay Indicator Error Equispaced (b) Comparison of relative l∞-norm decay. Figure 6.13: Comparison of the relative error decay for the different methods. In both panels, the errors obtained by means of the a posteriori error indicator appear in continuous blue fashion with circles, the errors obtained using the real error appear in dashed red fashion with squares and the errors obtained using equispaced parameters are shown in dash-dotted green fashion with asterisks. basis, therefore reducing the maximum error. Let us recall that the simple equispaced polynomial interpolation may present exponential decay if there is an analytic depen- 6.2 numerical results 137 dence with respect to the parameter in a domain large enough. This fact possibly plays a role here, as the flow is quite smoothly changing as the Reynolds number changes. Therefore, we conclude the validation tests for Kolmogórov’s indicator (47) with very promising results. It has been shown that the application of this indicator provides a similar behavior in error decay as using the exact error, achieving relative errors of order 10−5 after eight iterations. It has also been proven very useful the development of Algorithm 1, as it allows the enrichment of the reduced basis, taking into account that Kolmogórov’s indicator (47) does not tend to 0 as Nincreases. CONCLUSIONS AND OPEN PROBLEMS In this part, we have introduced a Reduced Basis model for turbulence at statistical equilibrium. This Reduced Basis model is based upon an a posteriori error estimation procedure that measures the distance from a trial solution to the theoretical k−5/3 energy spectrum provided by Kolmogórov’s theory. Upon this error estimation procedure, we have carried on an error analysis that clarifies the role of the energy-spectrum error estimator. We have stated that the error between the finite element solution and the reduced one is driven by this estimator, plus the projection error on the reduced space. In order to build the reduced problem, we have presented a coupled POD+Greedy approach (Section 1.3.3). We have adapted the criteria for the selection of the new parameter in the Greedy algorithm, as this indicator is not intended to vanish as the dimensionality of the reduced model increases. To avoid the possible stagnation situations, we compare the energy spectrum of the trial solution with that of the current RB solution as a surrogate of the FOM solution. In the field of the applications, we have tested the RB Smagorinsky model for some academic, yet meaningful, 2D flows, as these present the Kolmogorov’s Theory energy spectrum profile. In the performed tests, the Reynolds number acts as the only parameter. We obtain speed-up ratios of computational time close to 23 , which is very satisfying for an evolution turbulence model. In addition, we have obtained relative errors very similar to those provided by the reduced model built using the actual relative error for the parameter selection. Both errors present a similar exponential decay as the dimension of the reduced spaces increases, in both l2 and l∞ parametric norms. We have also compared these errors against the ones obtained for the reduced model built using equispaced parameters, that present a much slower decay in L∞ parametric norm. One of the main interests of this novel a posteriori error indicator based upon physical properties of the models is that it applies to any kind of parameter (either geometrical or physical) and to any kind of numerical discretization. The only requirement is that the targeted flows have a well-developed inertial spectrum, and that the numerical 139