Hidden Fluid Mechanics: A Navier-Stokes Informed Deep Learning Framework for Assimilating Flow Visualization Data

08/13/2018 ∙ by Maziar Raissi, et al. ∙ 2

We present hidden fluid mechanics (HFM), a physics informed deep learning framework capable of encoding an important class of physical laws governing fluid motions, namely the Navier-Stokes equations. In particular, we seek to leverage the underlying conservation laws (i.e., for mass, momentum, and energy) to infer hidden quantities of interest such as velocity and pressure fields merely from spatio-temporal visualizations of a passive scaler (e.g., dye or smoke), transported in arbitrarily complex domains (e.g., in human arteries or brain aneurysms). Our approach towards solving the aforementioned data assimilation problem is unique as we design an algorithm that is agnostic to the geometry or the initial and boundary conditions. This makes HFM highly flexible in choosing the spatio-temporal domain of interest for data acquisition as well as subsequent training and predictions. Consequently, the predictions made by HFM are among those cases where a pure machine learning strategy or a mere scientific computing approach simply cannot reproduce. The proposed algorithm achieves accurate predictions of the pressure and velocity fields in both two and three dimensional flows for several benchmark problems motivated by real-world applications. Our results demonstrate that this relatively simple methodology can be used in physical and biomedical problems to extract valuable quantitative information (e.g., lift and drag forces or wall shear stresses in arteries) for which direct measurements may not be possible.

READ FULL TEXT VIEW PDF

Authors

page 11

page 12

page 14

page 18

page 20

page 22

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Recent advances in deep learning and computational resources in addition to new data recordings and sensor technologies have the potential to revolutionize our understanding of the physical world in modern application areas such as image recognition krizhevsky2012imagenet, drug discovery gawehn2016deep, and bioinformatics alipanahi2015predicting. Moreover, many concepts from deep learning can be integrated with classical methods in mathematical physics to help us solve sophisticated data assimilation and inverse problems involving complex dynamical systems. This integration of nonlinear dynamics and deep learning opens the door for principled methods for model construction and predictive modeling. Furthermore, recent developments in physics-informed deep learning

have enabled us to train deep neural networks with millions of parameters using very few training examples, for cases where the available data are known to respect a given physical law described by a system of partial differential equations

raissi2017physics_I; raissi2017physics_II. The prior knowledge of an underlying physical law introduces important structure that effectively regularizes the minimization procedure in the training of neural networks, and enables them to generalize well even when only a few training examples are available. An important class of these physical laws govern the fluid motions in nature, which will be the focus of the present work.

Quantifying the flow dynamics essentially requires the knowledge of spatio-temporal fluid velocity and pressure fields, and has been the centerpiece of experimental and theoretical fluid mechanics for centuries. Traditionally, this has been achieved by measuring the instantaneous fluid velocity and pressure experimentally, or by solving the governing equations of fluid mechanics in a precisely defined geometry with proper initial and boundary conditions. Despite significant advances in experimental fluid mechanics (e.g., with the emergence of particle image velocimetry), the use of measurements to reliably infer fluid velocity and pressure/stress fields is not a straightforward task. Furthermore, although experimental measurements of external flows (e.g., flow past a bluff object) are obtained relatively easily, albeit in small subdomains, leading to two or even three dimensional vector fields for velocity, the quantification of velocity fields for internal flows (e.g., blood flow in the vascular networks) could become very difficult or impractical. From the theoretical standpoint, the governing equations of fluid mechanics have been derived from conservation laws (i.e., conservation of mass, momentum and energy) leading to partial differential equations such as the well-known Navier-Stokes equations for Newtonian fluids

batchelor2000introduction. With the increasing computational power, accurate solutions of such equations in a “forward” setting are now available with lower cost even in the turbulence regime using Direct Numerical Simulations (DNS), hence resolving all energetically important scales down to Kolmogorev dissipation scales. However, the assimilation of experimental fluid mechanics data into the mathematical models to infer velocity and pressure fields in an “inverse” setting has not been adequately addressed so far. The main objective of this article is to extend and incorporate the aforementioned physics-informed deep learning algorithms raissi2017physics_I; raissi2017physics_II to leverage the hidden physics of fluid mechanics to infer the latent quantities of interest (e.g., the velocity and pressure fields) through minimalistic usage of data obtained by the “passive scalar” transport (e.g., transport of dye or smoke).

Let us consider the transport of a passive scalar field by a velocity field . Such a problem arises, for example, while studying the spreading of smoke or dye advected by a given velocity field and subject to molecular diffusion. We expect the dynamics of to be governed by an effective transport equation vergassola1997scalar written in the form of

(1)

where is the molecular diffusivity and is the density. Transport of scalar fields in fluid flow has been studied in numerous applications such as aerodynamics, biofluid mechanics, and non-reactive flow mixing to name a few. The use of smoke in wind tunnels or dye in water tunnels for flow visualization and quantification has long been practiced in experimental fluid mechanics barlow2015low. Moreover, recent techniques in planar laser induced fluorescence imaging combined with particle image velocimetry have been developed to assess the relationships between scalar and velocity/vorticity fields koochesfahani1985laser; crimaldi2008planar. The use of scalar transport in conjunction with advanced imaging modalities to quantify blood flow in the vascular networks is now a common practice. For example, coronary computed tomography (CT) angiography is typically performed on multidetector CT systems after the injection of non-diffusible iodine contrast agent, which allows coronary artery visualization and the detection of coronary stenoses voros2011coronary. Another example is the quantification of cerebral blood flow, which is detrimental in the prognostic assessments in stroke patients using a contrast agent and perfusion CT wintermark2002prognostic, and in cognitive neuroscience with the use of functional magnetic resonance imaging that only relies on the blood-oxygen-level dependent contrast to measure brain activity aron2007triangulating.

Inspired by recent developments in physics-informed deep learning raissi2017physics_I; raissi2017physics_II and deep hidden physics models raissi2018deep, we propose to leverage the hidden physics of fluid mechanics (i.e., the Navier-Stokes equations) and infer the latent quantities of interest (e.g., the velocity and pressure fields) by approximating them using deep neural networks. This choice is motivated by modern techniques for solving forward and inverse problems associated with partial differential equations, where the unknown solution is approximated either by a neural network (raissi2017physics_I; raissi2017physics_II; raissi2018forwardbackward; raissi2018multistep) or a Gaussian process (raissi2018numerical; raissi2018hidden; raissi2017inferring; raissi2017machine; raissi2017parametric; perdikaris2017nonlinear; raissi2016deep). Moreover, placing a prior on the solution itself is fully justified by the similar approach pursued in the past century using classical methods of solving partial differential equations such as finite elements and spectral methods, where one would expand the unknown solution in terms of an appropriate set of basis functions. Our focus here is the transport of a passive scalar by incompressible Newtonian flows in both unbounded geometries (i.e., external flows) and confined internal flows. We demonstrate the success of our Navier-Stokes informed deep learning algorithm by recovering the flow velocity and pressure fields solely from time series data collected on the passive scalar in arbitrary domains.

2 Problem Setup and Solution Methodology

Starting with the transport equation (1), we can rewrite it in the following non-dimensional form

(2)

which governs the evolution of the normalized concentration of a passive scalar transported by an incompressible Newtonian fluid whose dynamics are described by the Navier-Stokes and continuity equations (also non-dimensional) given below

(3)

In equations (2) and (3), , , and are the , , and components of the velocity field , respectively, and denotes the pressure. A passive scalar is a diffusive field in the fluid flow that has no dynamical effect (such as the effect of temperature in a buoyancy-driven flow) on the fluid motion itself. Smoke and dye are two typical examples of passive scalars. Here, Re represents the Reynolds number, a dimensionless quantity in fluid mechanics, and Pec denotes the Péclet number, a dimensionless quantity relevant in the study of transport phenomena in flowing fluids. In this work, we assume that the only observables are noisy data on the concentration of the passive scalar. Given such data, scattered in space and time, we are interested in inferring the latent (hidden) quantities , , , and . To solve the aforementioned data assimilation problem, we would like to design an algorithm that is agnostic to the geometry as well as the initial and boundary conditions. This is enabling as it will give us the flexibility to work with data acquired in arbitrarily complex domains such as human arteries or brain aneurysms.

One question that would naturally arise is whether the information on the passive scalar in the training domain and near its boundaries is sufficient to result in a unique velocity field. The answer is that normally there are no guarantees for unique solutions unless some form of boundary conditions are explicitly imposed on the domain boundaries. However, as shown later for the benchmark problems studied in the current work, an informed selection of the training boundaries in the regions where there are sufficient gradients in the concentration of the passive scalar would eliminate the requirement of imposing velocity and pressure boundary conditions. In addition to the proper design of the training domain and the use of the passive scalar , we can improve the model predictions further by introducing an auxiliary variable (essentially the complement of ) that satisfies the transport equation

(4)

This can be clearly seen by a change of variable from to and using equation (2). The complementary nature of the auxiliary variable helps the algorithm better detect the geometry and the corresponding boundary conditions. This makes the training algorithm agnostic to the physical geometry of the problem, hence, keeping the implementation and training significantly simpler. In addition, the region of interest becomes very flexible to be chosen with its boundaries no longer required to be the physical boundaries. Another advantage of the use of is the improvement in the accuracy of the algorithm predictions with no additional cost.

Figure 1: Navier-Stokes informed neural networks:

A plain vanilla densely connected (physics uninformed) neural network, with 10 hidden layers and 50 neurons per hidden layer per output variable (i.e.,

neurons per hidden layer), takes the input variables and outputs , and

. As for the activation functions, we use

. For illustration purposes only, the network depicted in this figure comprises of 2 hidden layers and 7 neurons per hidden layers. We employ automatic differentiation to obtain the required derivatives to compute the residual (physics informed) networks , , , , , and . If a term does not appear in the blue boxes (e.g., or

), its coefficient is assumed to be zero. It is worth emphasizing that unless the coefficient in front of a term is non-zero, that term is not going to appear in the actual “compiled” computational graph and is not going to contribute to the computational cost of a feed forward evaluation of the resulting network. The total loss function is composed of the regression loss of the passive scalar

and its complement on the training data, and the loss imposed by the differential equations . Here, denotes the identity operator and the differential operators , and are computed using automatic differentiation and can be thought of as “activation operators”. Moreover, the gradients of the loss function are backpropogated through the entire network to train the parameters using the Adam optimizer.

Following our earlier work on physics-informed deep learning raissi2017physics_I; raissi2017physics_II and deep hidden physics models raissi2018deep, we approximate the function

by a deep neural network and obtain the following Navier-Stokes informed neural networks (see figure 1) corresponding to equations (2), (3), and (4); i.e.,

(5)

We acquire the required derivatives to compute the residual networks , , , , , and

by applying the chain rule for differentiating compositions of functions using automatic differentiation

baydin2015automatic. A schematic representation of the resulting Navier-Stokes informed neural networks is given in figure 1. It is worth emphasizing that automatic differentiation is different from, and in several respects superior to, numerical or symbolic differentiation – two commonly encountered techniques of computing derivatives. In its most basic description baydin2015automatic

, automatic differentiation relies on the fact that all numerical computations are ultimately compositions of a finite set of elementary operations for which derivatives are known. Combining the derivatives of the constituent operations through the chain rule gives the derivative of the overall composition. This allows accurate evaluation of derivatives at machine precision with ideal asymptotic efficiency and only a small constant factor of overhead. In particular, to compute the required derivatives we rely on Tensorflow

abadi2016tensorflow

, which is a popular and relatively well documented open-source software library for automatic differentiation and deep learning computations. In TensorFlow, before a model is run, its computational graph is defined statically rather than dynamically as for instance in PyTorch

paszke2017automatic. This is an important feature as it allows us to create and compile the computational graph for the Navier-Stokes informed neural networks (5) only once and keep it fixed throughout the training procedure. This leads to significant reduction in the computational cost of the proposed framework.

The shared parameters of the neural networks for , , , , , and can be learned by minimizing the following sum of squared errors loss function

Here, and the first two terms correspond to the training data on the concentration of the passive scalar while the last term enforces the structure imposed by equations (2), (3), and (4) at a finite set of measurement points whose number and locations are taken to be the same as the training data. However, it should be pointed out that the number and locations of the points on which we enforce the set of partial differential equations could be different from the actual training data. Although not pursued in the current work, this could significantly reduce the required number of training data on the concentration of the passive scalar.

3 Results

To generate high-resolution datasets for different benchmark problems studied in the current work, we have employed the spectral/hp-element solver NekTar in which the Navier-Stokes equations (3) along with the transport equation (2) are approximated using high-order semi-orthogonal Jacobi polynomial expansions karniadakis2013spectral. The numerical time integration is performed using a third-order stiffly stable scheme until the system reaches its stationary state. In what follows, a small portion of the resulting dataset corresponding to this stationary solution will be used for model training, while the remaining data will be used to validate our predictions. Our algorithm is agnostic to the choice of initial and boundary conditions as well as the geometry. However, we choose to provide this information for every benchmark problem for the sake of completeness and reproducibility of the numerically generated data. Moreover, all data and codes used in this manuscript are publicly available on GitHub at https://github.com/maziarraissi/HFM.

To obtain the results reported in the remainder of this manuscript, we represent each of the functions , , , , , and by a 10-layer deep neural network with 50 neurons per hidden layer (see figure 1). As for the activation functions, we use . In general, the choice of a neural network’s architecture (e.g., number of layers/neurons and form of activation functions) is crucial and in many cases still remains an art that relies on one’s ability to balance the trade off between expressivity and trainability of the neural network (raghu2016expressive). Our empirical findings so far indicate that deeper and wider networks are usually more expressive (i.e., they can capture a larger class of functions) but are often more costly to train (i.e., a feed-forward evaluation of the neural network takes more time and the optimizer requires more iterations to converge). Moreover, the sinusoid (i.e., ) activation function seems to be numerically more stable than , at least while computing the residual neural networks , (see equation (5)). However, these observations should be interpreted as conjectures rather than as firm results111We encourage the interested reader to check out the codes corresponding to this paper on GitHub at https://github.com/maziarraissi/HFM and experiment with different choices for the neural networks’ architectures.. In this work, we have tried to choose the neural networks’ architectures in a consistent fashion throughout the manuscript. Consequently, there might exist other architectures that could possibly improve some of the results reported in the current work.

As for the training procedure, our experience so far indicates that while training deep neural networks, it is often useful to reduce the learning rate as the training progresses. Specifically, the results reported in the following are obtained after , , and

consecutive epochs of the Adam optimizer

kingma2014adam with learning rates of , , and , respectively. Each epoch corresponds to one pass through the entire dataset. The total number of iterations of the Adam optimizer is therefore given by times the number of data divided by the mini-batch size. The mini-batch size we used is and the number of data points are clearly specified in the following on a case by case basis. Every iterations of the optimizer takes around and seconds, respectively, for two and three dimensional problems on a single NVIDIA Titan X GPU card.

3.1 External flows

As a first example, we consider the prototypical problem of a two dimensional flow past a circular cylinder, known to exhibit rich dynamic behavior and transitions for different regimes of the Reynolds number . Assuming a non-dimensional free stream velocity , cylinder diameter , and kinematic viscosity , the system exhibits a periodic steady state behavior characterized by an asymmetrical vortex shedding pattern in the wake of the cylinder, known as the Kármán vortex street batchelor2000introduction. Importantly, the passive scalar is injected at the inlet within the interval [-2.5, 2.5] using a step function (see figure 2). The choice of this boundary condition only depends on the region of interest in which velocity and pressure fields are inferred using the concentration field. Furthermore, a constant value of diffusivity for the passive scalar is assumed to be given by resulting in . Here, we use zero-slip and zero-concentration boundary conditions on the cylinder wall. Whereas the Péclet number is chosen to be equal to the Reynolds number, there is no restriction on its value as shown for the internal flow benchmark problem. Physically, for the flow of gases such as air with smoke as a passive scalar . This is not, however, the case for liquid flows with dye as the passive scalar since diffusivity of dyes are typically smaller than most fluids, which leads to higher Péclet numbers.

Figure 2: 2D Flow past a circular cylinder: 2D sketch of the simulation domain in which data on the concentration field is generated. We have assumed a uniform free stream velocity profile at the left boundary, a zero pressure outflow condition imposed at the right boundary located downstream of the cylinder, and periodicity for the top and bottom boundaries of the domain. The passive scalar is injected at the inlet from [-2.5, 2.5]. A completely arbitrary training domain in the shape of a flower is depicted in the wake of the cylinder. No information on the velocity is given for this training domain. Another training domain is also shown by a rectangle that includes the cylinder. Information on the velocity is only given on the left boundary of the rectangular domain shown by the dashed line.

A representative snapshot of the input data on the concentration field is shown in figure 2. As illustrated in this figure, the shape and extent of the boundaries of the training domains that we choose for our analysis could be arbitrary and may vary by problem. However, it should be pointed out that there are two important factors that need to be considered when choosing the training domain. First, the concentration field of the passive scalar must be present within the training domain, such that its information can be used to infer other flow variables. Second, to avoid the need for specifying appropriate boundary conditions for velocities, there must exist enough concentration gradient normal to the boundaries (i.e., ) in order for our method to be able to infer a unique solution for the velocity field. As shown in figure 3, excellent agreement can be achieved between the predictions of our algorithm and the exact data within a completely arbitrary training domain downstream of the cylinder, while the inclusion of auxiliary variable can improve the accuracy of predictions further (see figure 4). Here, no input information other than the concentration field for the passive scalar is passed to the algorithm.

Figure 3: 2D Flow past a circular cylinder: A representative snapshot of the input data on the concentration of the passive scalar is shown in the top left panel of this figure, where on the right the same concentration field is reconstructed based on the predictions of our algorithm. The algorithm is capable of accurately reconstructing the velocity and the pressure fields shown in the third row. The exact velocity and pressure fields at the same point in time are plotted for comparison in the second row. Note that the pressure is off by a constant since this is an incompressible flow.
Figure 4: 2D Flow past a circular cylinder: Relative errors between predictions of the model and the corresponding exact concentration, velocity, and pressure fields for the arbitrary training domain downstream of the cylinder. 31 million data points, scattered in space and time, are used both to regress the concentration field and enforce the corresponding partial differential equations. Furthermore, the training is performed with the auxiliary variable (blue line) and without (red dashed line) for comparison.

Next, we choose the region of interest to contain the cylinder (see the rectangular domain in figure 2) so that the fluid forces acting on the cylinder can be inferred. As shown in figure 5, the algorithm is capable of accurately reconstructing the velocity and the pressure fields without having access to sufficient observations of these fields themselves. Note that we use a very small training domain that cannot be used in classical computational fluid dynamics to obtain accurately the Navier-Stokes solutions. It should be emphasized that other than the velocity on the left boundary, no other boundary conditions are given to the algorithm. We need to impose a Dirichlet boundary condition for the velocity at the left boundary simply because the observations of the concentration are not providing sufficient information (i.e., gradients) in front of the cylinder. More notably, there is no need to impose the no-slip condition on the cylinder wall. The presence of the normal concentration gradient naturally allows our algorithm to infer the zero velocity condition on the walls. In regards to the predicted pressure field, we note that due to the nature of the Navier-Stokes equations for incompressible flows, the pressure field is only identifiable up to a constant.

Figure 5: 2D Flow past a circular cylinder: A representative snapshot of the input data on the concentration of the passive scalar is shown in the top left panel of this figure, where on the right the same concentration field is reconstructed based on the predictions of our algorithm. The algorithm is capable of accurately reconstructing the velocity and the pressure fields shown in the third row. The exact velocity and pressure fields at the same point in time are plotted for comparison in the second row. Note that the pressure is off by a constant since this is an incompressible flow.

The learned pressure and velocity fields can be used to obtain the lift and drag forces exerted on the cylinder by the fluid as shown in figure 6. The fluid forces on the cylinder are functions of the pressure and velocity gradients. Consequently, having trained the neural networks, we can use

to obtain the lift and drag forces, respectively. Here, is the outward normal on the cylinder and is the arc length on the surface of the cylinder. We use the trapezoidal rule to approximately compute these integrals. Interestingly, the predicted lift and drag forces are in excellent agreement with the exact ones, both in terms of the frequency of oscillations and the amplitude. The resulting error is within 1% of the test data. Some discrepancy, however, can be observed for the few initial and final time instants, which can be attributed to the lack of data. This can be further clarified by the relative -norm of error results shown in figure 7. Lack of training data on for and leads to weaker neural network predictions for the initial and final time instants. Thus, one should take this into consideration when inference is required within a certain time interval. Note that to compute lift and drag forces, the gradient of the velocity has to be computed on the wall, which can be done analytically using the parametrized surrogate velocity field (i.e., the neural networks). Although no numerical differentiation is needed to compute the gradients, the integration of forces on the surface of the cylinder is approximated by a summation.

Figure 6: Non-dimensional lift and drag forces on 2D cylinder: Comparison between the exact (spectral element solution) and predicted lift (left) and drag (right) forces per unit length of the cylinder.

In addition to the velocity and pressure fields, it is possible to discover other unknown parameters of the flow field such as the Reynolds and Péclet numbers. Although these parameters were prescribed in the 2D flow past the cylinder example, we have tested a case in which both parameters are free to be learned by the algorithm. The results are given in table 1, which shows very good agreement with the exact values. From the practical standpoint, the passive scalar diffusivity and the fluid viscosity (hence, their ratio i.e., the Prandtl number ) may be known in advance. Therefore, discovering the Reynolds number would be sufficient whereas the Péclet number can be computed by Pec = Re Pr.

Figure 7: 2D Flow past a circular cylinder: Relative errors between predictions of the model and the corresponding exact concentration, velocity, and pressure fields for the rectangular training domain that contains the cylinder. 2.8 million data points scattered within the training domain and in time are used for both regressing the concentration field and enforcing the corresponding partial differential equations.
Exact Learned Rel. Error
Pec 100 92.39 7.60%
Re 100 92.47 7.52%
Table 1: 2D flow past a circular cylinder: Learned Reynolds and Péclet numbers considered as free parameters of the model.

The previous benchmark example was a two-dimensional (2D) flow, where we could safely neglect the -coordinate and -component of the velocity from the input and output variables, respectively. For the flow past a cylinder, if we simply increase the Reynolds number beyond a threshold value of ma2002low, the spanwise velocity becomes more prominent due to the effect of the so-called “vortex stretching” batchelor2000introduction; karniadakis1992three. To test the capability of the proposed algorithm in inferring three-dimensional (3D) flow fields, we design another prototype problem of the 3D flow past a finite-size circular cylinder confined between two parallel plates as shown in figure 8, where similar to the previous example, we set . Downstream of the cylinder, the flow exits to an open region, which causes strong 3D effects in the wake of the cylinder. Hence, we set the training domain in the wake of cylinder.

Figure 8: 3D flow past a circular cylinder: 3D sketch for the prototypical problem of flow past a circular cylinder of diameter . Two parallel planes are located apart along the axis. Flow is bounded between the plates for along axis. Periodic boundary conditions are prescribed parallel to the plane at and , and parallel to the plane at and starting from the location where the walls end. All the physical boundaries are shown in black, whereas the gray box located behind the cylinder shows the domain of interest where model training is performed. Zero-Neumann boundary conditions are imposed for velocities and concentration along with zero pressure at the outflow located downstream of the cylinder. A uniform is imposed at the inlet located upstream of the cylinder. Passive scaler is injected at the inlet through a finite region . Zero-velocity and concentration boundary conditions are imposed on each physical boundary. Note that no information on the surfaces of the training (gray box) domain is given.

A representative snapshot of the input data on the concentration field in the wake of the cylinder is plotted in the top left panel of figure 9. The iso-surfaces for all of the fields are also plotted for comparison between the exact data and the predictions of our algorithm. The algorithm is capable of accurately reconstructing the velocity and the pressure fields without having access to any observations of these fields. In particular, no information on the velocity is given on the boundaries of the domain of interest in the training phase. Qualitative comparison between the predictions of our algorithm and the exact data shows good agreement for 3D flows as well, however, the relative -norm errors for velocity and pressure fields are slightly higher than the values for the 2D flow predictions as shown in figure 10. Similar to 2D results, the neural network training lack sufficient data for the initial and final time instants that leads to larger errors in the predictions.

Figure 9: 3D flow past a circular cylinder: The iso-surfaces of exact data and the predictions of our algorithm on the concentration of the passive scalar are shown in top row for a representative time instant within the selected training domain. Using the information on the concentration only, the velocity fields are inferred, are shown in the third row, and are compared with the exact data in the second row. In addition, the exact and predicted pressure fields are plotted in the last row. Note the range of contour levels are set equal between the exact and predicted iso-surfaces of concentration fields and velocity components for better comparison.
Figure 10: 3D flow past a circular cylinder: Relative errors between predictions of the model and the corresponding exact concentration, velocity, and pressure fields. 37 million data points scattered within the training box and in time are used for both regressing the concentration field and enforcing the corresponding partial differential equations.

3.2 Internal flows

We now turn our attention to an important class of flows in confined geometries also known as “internal flows”. While measuring average flow velocity and pressure in ducts, pipes and even blood vessels is now a common practice, quantifying the spatial fields, specifically the shear stresses on the boundaries, is not a trivial task. The first benchmark problem considered here is the transient flow over an obstacle in a 2D channel as shown in figure 11. To make the flow unsteady, we impose a velocity waveform at the inlet of the channel. Passive scalar value on the boundaries is set to zero, where the boundaries are assumed to be impenetrable. It is clear that the auxiliary variable is equal to one on the channel boundaries. In other words, while the passive scalar concentration is convected downstream from the inflow, the complement scalar is “virtually” infused from the boundaries to the flow stream as shown in figure 11

. The presence of the obstacle breaks the symmetry in the flow for which analytical solutions do not exist. Thus, to estimate the velocity field, direct measurements or forward numerical simulations are required.


A representative snapshot of the input data on the concentration field along with its neural network approximation are depicted in the top panel of figure 12. As shown in figure 12, predictions of the algorithm for velocity and pressure fields are in excellent agreement with the exact data. Note that no information for velocity is given on the boundaries of the channel or the obstacle. It should be mentioned that while the concentration gradients normal to the physical boundaries are sufficient to infer the flow variables without the knowledge of the boundary conditions (i.e., no-slip velocity on the wall), the free information from the auxiliary variable in addition to helps the algorithm improve the accuracy of predictions (shown in figure 14) at practically no additional cost.

Figure 11: 2D channel flow over an obstacle: Contours of the concentration field within a 2D channel at a representative time instant. A pulsatile velocity profile () is imposed at the inlet (shown as the inset on the left), and a uniform concentration for the passive scalar is injected to the channel. At the outlet, a zero-Neumann boundary condition for the velocity and concentration is considered, while the pressure is set to zero. Furthermore, the velocity and concentration on the walls are set to zero. The flow is calculated based on the mean inflow velocity and channel hight , whereas we choose a smaller diffusion constant for the passive scalar leading to higher . The training domain (white rectangle) is considered to contain the obstacle with the upper and lower boundaries being the physical wall boundaries. Also shown are contours of the auxiliary variable , which signifies the important role of this variable in helping the algorithm detect the boundaries of the training domain.

One of the real-world examples of the above benchmark problem is encountered in cardiovascular fluid mechanics, where one or multiple coronary arteries are partially blocked by atherosclerosis plaques formed by the lipid accumulation hansson2005inflammation. An accurate prognosis of these “stenotic” vessels is an extremely important step in the treatment planning and decision-making process for patients with cardiovascular disease budoff2008diagnostic. Direct measurements of pressure in the vessel is invasive and bears high risk, whereas computational fluid modeling of blood in the coronaries requires correct reconstruction of the geometry as well as the knowledge of all the inflow/outflow boundary conditions. The proposed inverse approach to infer the velocity and pressure fields with the use of a passive scalar, however, offers a promising alternative to the conventional methods. Here, the passive scalar could be the bolus dye that is typically injected to the blood stream for the purpose of blood flow monitoring and medical imaging.

Figure 12: 2D channel flow over an obstacle: A representative snapshot of the input data on the concentration field within the training domain is plotted in the top left panel alongside the prediction of our algorithm. The algorithm is capable of accurately reconstructing the velocity and the pressure fields without having access to any observations of these fields (shown in the second and third rows). Furthermore, no boundary conditions are specified on the boundaries of the training domain including the physical wall boundaries.

Using the predicted velocity fields, we are able to compute shear stresses everywhere in the training domain. Of particular interest are wall shear stresses, which can be computed using the following equations in 2D;

Here, is the outward normal on the boundary of the domain. Note that to compute wall shear stresses , the gradient of the velocity is required, which can be computed using automatic differentiation. Wall shear stresses are important quantities of interest in many biological processes e.g., in the pathogenesis and progression of vascular diseases that cause aortic aneurysms. Here, we have estimated the temporal wall shear stress magnitudes, , acting on the lower wall using the predictions of neural networks for the 2D channel flow over the obstacle, and plotted them against the results from the spectral/hp element solver in figure 13. The results show excellent agreement between the predictions and numerical estimations. Note that the numerical estimations suffer from a slight aliasing effect (noisy oscillations in the direction) for which dealiasing is required to retrieve a smooth distribution of wall shear stresses. Interestingly, the neural networks prediction is smooth as the velocity gradients of a parameterized velocity field are taken analytically. Furthermore, similar to the previous benchmark problems, slight discrepancies can be observed close to the initial and final times () due to the lack of training data. This is further explained in figure 14, where relative -norm of errors in the velocity fields show a peak close to the initial and final time instants out of the entire prediction time interval. Also shown in the figure is the enhancement in the model prediction accuracy when the auxiliary variable is used as opposed to when is not included during the training phase.

Figure 13: 2D channel flow over an obstacle: Wall shear stress magnitude, , computed on the lower wall () in the training domain of figure 11 as a function of time (). Strong wall shear stress values are observed on the surface of the cylinder (). The plot also shows the periodicity in the flow field.

Similar to the external flow problem past the cylinder, the algorithm is able to learn the Reynolds and Péclet numbers as free parameters. The predicted values of these two parameters are given in table 2, which shows excellent agreement with the exact values. Note that unlike the flow past the cylinder, the physical Re and Pec numbers are not the same. Remarkably, the algorithm is capable of inferring the underlying physics without knowledge of initial and boundary conditions and physical boundaries, where the information is provided only on the transported passive scalar.

Figure 14: 2D channel flow over an obstacle: Relative errors between predictions of the model and the corresponding exact concentration, velocity, and pressure fields. 10 million data points scattered within the training domain and in time are used for both regressing the concentration field and enforcing the corresponding partial differential equations. Furthermore, the training is performed with the auxiliary variable (blue line) and without (red dashed line) for comparison.
Exact Learned Rel. Error
Pec 180 178.95 0.58%
Re 60 59.92 0.12%
Table 2: 2D channel flow over an obstacle: Learned Reynolds and Péclet numbers considered as free parameters of the model.

To further illustrate the implications of the Navier-Stokes informed neural networks in addressing real-world problems, we consider 3D physiologic blood flow in a realistic intracranial aneurysm (ICA) shown in figure 15. The aneurysm is located in the cavernous segment of the right internal carotid artery at the level of the eye and beneath the brain baek2012convergence. Exact concentration fields are generated numerically using realistic boundary conditions, which is a physiologic flow waveform at the inlet along with a uniform concentration for the passive scalar. The strength of the proposed algorithm is its ability to stay agnostic with respect to the geometry as well as initial and boundary conditions. Hence, it is possible to focus only on the regions where velocity and pressure fields are needed, which will significantly reduce the size of data and the cost of training.

We first crop the aneurysm sac out from the rest of geometry, and then use only the scalar data within the ICA sac (right panel of figure 15

) for training where no information is used for the boundary conditions. The exact and predicted concentration, velocity and pressure fields within the ICA sac at a sample time instant are then interpolated on two separate planes perpendicular to

- and -axis, which are shown in figure 16. We observe excellent agreement between the exact and predicted fields given the complexity of the flow field. Furthermore, the relative -norm of errors plotted in figure 17 show a significant drop in the prediction errors away from the initial and final time instants, which justifies an accurate analysis of model predictions in the time interval .

Figure 15: A 3D intracranial aneurysm: The middle panel shows the simulation domain and pressure field at a time instant, whereas on the right the training domain containing only the ICA sac is shown. Two perpendicular planes have been used to interpolate the exact data and the predicted ones for plotting 2D contours in figure 16. A physiologic flow waveform (shown in the inset figure) is prescribed at the inlet along with the uniform concentration for the passive scalar. At the outlet, zero-Neumann boundary conditions are imposed for velocities and concentration, whereas a “Windkessel” type boundary condition is used for the pressure to represent the truncated geometry downstream grinberg2008outflow. The Reynolds and Péclet numbers are estimated based on the mean velocity and lumen diameter at the inlet. Using the kinematic viscosity of blood and assuming a diffusivity constant for the passive scalar equal to the viscosity, we obtain .
Figure 16: A 3D intracranial aneurysm: Contours of the exact fields and model predictions are plotted on two perpendicular planes (shown in figure 15) for concentration , velocity , and pressure fields in each row. The first two columns show the results interpolated on a plane perpendicular to the -axis, and the next two columns are plotted for a plane perpendicular to the -axis. Note that the range of contour levels are set equal for all fields for better comparison.
Figure 17: A 3D intracranial aneurysm: Relative errors between predictions of the model and the corresponding exact concentration, velocity, and pressure fields. 29 million data points scattered within the aneurysm sac and in time are used for both regressing the concentration field and enforcing the corresponding partial differential equations.

4 Discussion and Concluding Remarks

The algorithm developed here is agnostic to the geometry, initial, and boundary conditions, which makes it highly flexible in choosing the domain of interest for data acquisition as well as subsequent training and predictions. Moreover, the current methodology allows us to construct computationally efficient and fully differentiable surrogates for velocity and pressure fields that can be further used to estimate other quantities of interest such as shear stresses and vorticity fields. The predictions presented here are among those cases where a pure machine learning algorithm or a mere scientific computing approach simply cannot reproduce. A pure machine learning strategy has no sense of the physics of the problem to begin with, and a mere scientific computing approach relies heavily on careful specification of the geometry as well as initial and boundary conditions.

Assuming the geometry to be known, to arrive at similar results as the ones presented in the current work, one needs to solve significantly more expensive optimization problems, using conventional computational methods (e.g., finite differences, finite elements, finite volumes, spectral methods, and etc.). The corresponding optimization problems involve some form of “parametrized” initial and boundary conditions, appropriate loss functions, and multiple runs of the conventional computational solvers. In this setting, one could easily end up with very high-dimensional optimization problems that require either backpropagating through the computational solvers

chen2018neural or “Bayesian” optimization techniques shahriari2016taking for the surrogate models (e.g., Gaussian processes). If the geometry is further assumed to be unknown (as is the case in this work), then its parametrization requires grid regeneration, which makes the approach almost impractical.

It must be mentioned that we are avoiding the regimes where the Navier-Stokes equations become chaotic and turbulent (e.g., as the Reynolds number increases). In fact, it should not be difficult for a plain vanilla neural network to approximate the types of complicated functions that naturally appear in turbulence. However, as we compute the derivatives required in equation (5), minimizing the loss function (2) might become a challenge raissi2018deep, where the optimizer may fail to converge to the right values for the parameters of the neural networks. It might be the case that the resulting optimization problem inherits the complicated nature of the turbulent Navier-Stokes equations. Hence, inference of turbulent velocity and pressure fields could be considered in future extensions of this line of research.

The emphasis of this study was to demonstrate, through several prototypical and realistic examples, the ability of the current algorithm to infer the hidden states of the system from partial knowledge of some relevant quantities by leveraging the known underlying dynamics of the system. We have verified that with the use of governing physical laws, our algorithm is able to make accurate predictions for complex 2D/3D flows. One possible limitation of the current study is the use of synthetically generated data on the passive scalar, which are relatively noiseless and clean compared to the realistic measurements. Whereas 3D reconstruction of a passive scalar from a stack of 2D projections is possible, it is not a common practice in the clinical and industrial settings due to the technological complexities as well as computational and processing costs. Thus, extracting information from 2D images directly by taking into account the angle of projections seems to be a more realistic approach. Furthermore, different imaging modalities (e.g., magnetic resonance vs. computed tomography angiography of bolus dye in coronary arteries) have different spatial/temporal resolutions. To address the lack of resolution in time, a viable strategy is to use a high-order time-stepping scheme such as implicit Runge-Kutta for resolving the equations in time using very few snapshots raissi2017physics_I; raissi2017physics_II. These are certainly important factors that have to be taken into account moving forward using the proposed algorithm for real-world fluid mechanics applications, and will be addressed carefully in future implementations.

In this work we have been operating under the assumption of Newtonian and incompressible fluid flow governed by the Navier-Stokes equations. However, the proposed algorithm can also be used when the underlying physics is non-Newtonian, compressible, or partially known. This in fact is one of the advantages of the algorithm in which other unknown parameters such as the Reynolds and Péclet numbers can be inferred in addition to the velocity and pressure fields. When the fluid is non-Newtonian (e.g., blood flow in small vessels), one can encode the momentum equations, where the “constitutive” law for the fluid’s stress-strain relationship is unknown, into a physics-informed deep learning algorithm. Having data on the velocity (or even the passive scalar), it should be possible to learn the constitutive law as well zhao2018active.

As demonstrated in this work, a direct implication of the current method is quantifying the hemodynamics in the vasculature. This could potentially have significant impact on the clinical diagnosis (especially with the non-invasive methods) of vascular diseases associated with important pathologies such as heart attack and stroke. Blood shear stresses acting on the vascular wall are crucial in the prognosis of a vascular disease and their quantification is significantly important clinically shaaban2000wall; zarins1983carotid. Using the proposed method, it is possible to estimate the wall shear stresses at no extra cost. This will simplify the complexities of the state-of-the-art methods in which extracting the exact boundaries of the vessels from the clinical images is required boussel2009phase. In addition to the shear stresses estimation that is of interest, quantification of blood distribution and oxygenation in cerebral arteries using modalities such as perfusion CT wintermark2002prognostic and functional magnetic resonance imaging calamante1999measuring has been practiced and is crucial in the assessment of brain activity for the majority of neurodegenerative diseases such as Alzheimer’s disease greicius2004default and in cognitive neuroscience aron2007triangulating. Furthermore, recent advances in imaging such as the functional photoacoustic microscopy allows 3D blood oxygenation imaging with “capillary-level” resolution yao2015high. Encoding the one-dimensional hemodynamics and oxygen transport equations into physics-informed neural networks, and using the concentration profiles from different imaging modalities will lead to quantifying brain hemodynamics at various length scales.

In conclusion, we have introduced an effective algorithm based on the recently developed physics-informed neural networks framework that is capable of encoding the underlying physical laws that govern a given dataset. Specifically, we have applied the algorithm to an important class of physical laws governing fluid motions. Our inverse approach in this article is unique as we have designed a methodology to infer velocity and pressure fields merely from the knowledge of the time-evolution of a passive scaler.

Acknowledgements

This work received support by the DARPA EQUiPS grant N66001-15-2-4055, the AFOSR grant FA9550-17-1-0013, the NIH grant U01HL116323, and the NSF grant DMS-1736088. Generation of data was performed on XSEDE resources supported by award No. TG-DMS140007. Moreover, all data and codes used in this manuscript will be publicly available on GitHub at https://github.com/maziarraiss/HFM.

References