{'fraction_non_alphanumeric': 0.03713582626353014, 'fraction_numerical': 0.009206518559308909, 'mean_word_length': 4.739456790123457, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Learning nonlocal constitutive models with neural networks

Constitutive and closure models play important roles in computational mechanics and computational physics in general. Classical constitutive models for solid and fluid materials are typically local, algebraic equations or flow rules describing the dependence of stress on the local strain and/or strain-rate. Closure models such as those describing Reynolds stress in turbulent flows and laminar--turbulent transition can involve transport PDEs (partial differential equations). Such models play similar roles to constitutive relation, but they are often more challenging to develop and calibrate as they describe nonlocal mappings and often contain many submodels. Inspired by the structure of the exact solutions to linear transport PDEs, we propose a neural network representing a region-to-point mapping to describe such nonlocal constitutive models. The range of nonlocal dependence and the convolution structure are derived from the formal solution to transport equations. The neural network-based nonlocal constitutive model is trained with data. Numerical experiments demonstrate the predictive capability of the proposed method. Moreover, the proposed network learned the embedded submodel without using data from that level, thanks to its interpretable mathematical structure, which makes it a promising alternative to traditional nonlocal constitutive models.


Introduction
Computational mechanics, and more generally, computational physics predominantly rely on solving partial differential equations (PDEs) describing the conservation laws for mass, momentum, energy, and species.
However, a complementary yet equally critical component to the problem specification is a closure model. Such closure models come in various forms, ranging from local, algebraic constitutive models (e.g., stressstrain relations for linear elastic solids, Newtonian fluids, and viscoelastic fluids) to empirical correlations (e.g., drag force laws in particle-laden flows) and partial differential equations (PDEs, e.g., Reynolds stresses transport in turbulent flows [1], intermittency transport in laminar-turbulence transition [2,3]). All these closure models play similar roles in computational mechanics in that they express the unknown quantities in the PDEs with quantities that are known or solved for. However, an important conceptual difference between them is that classical constitutive models describe the macroscopic property of the materials, while the PDE-based models describe the property of the flow and thus can be considered generalized constitutive relations. In particular, the PDE-based models in generalized constitutive relations tend to be nonlocal because the unresolved physics they aim to represent is nonlocal. This is in stark contrast to the classical local constitutive models for material behaviors, which are usually local. As a result, such PDE-based constitutive models are much more difficult to calibrate due to their high complexities, the presence of submodels with interactions, and the lack of data. In the present work we explore a new approach for this class of nonlocal, PDE-based constitutive models in the context of data-driven modeling.

Examples and challenges in nonlocal, PDE-based constitutive models
An important class of nonlocal constitutive models in transport problems takes the form of convectiondiffusion-reaction PDEs. This is in contrast to the classical constitutive models, which are typically linear or nonlinear algebraic relations between fluxes (e.g., heat flux and stress) and the gradient of primary variables (e.g., displacement, velocity, and temperature). Here we use the Reynolds stress transport models for turbulent flows as an example to motivate the current work, but it is worth noting that such nonlocal, PDE-based constitutive models are representative in a broad range of applications. Turbulent flows are described by the Navier-Stokes equations, which can be written for incompressible flows of constant density ρ and kinematic viscosity ν as: where v and p denote velocity and pressure, respectively, while t indicates temporal coordinate, and ∇ indicates spatial derivatives. Since the flows involve turbulent structures of a wide range of scales for high Reynolds number flows in practice, resolving those scales are prohibitively expensive. Consequently, engineering modeling often relies on descriptions of the mean flow as opposed to the instantaneous flow above. The mean flow equations are obtained by decomposing the flow into mean component and fluctuation components, i.e., v =v + v and p =p + p (wherev denotes mean velocity, written as u for brevity hereafter, andp is the mean pressure). Taking the ensemble average of the equation above leads to the Reynolds-averaged Navier-Stokes (RANS) equations describing the mean flow: where the mean velocity u and mean pressurep are the primary variables to be solved for. The RANS equations have a similar form to the Navier-Stokes equation except for the velocity covariance term τ ≡ −v i v j , referred to as Reynolds stress tensor [1]. In this paper we call the Reynolds stress τ as the closure variable, for which a closure model is needed for the RANS equations (2), i.e., the primary equations, to be complete. Such closure models are called turbulence models, which include many categories such as the commonly used eddy viscosity models and the more advanced Reynolds stress transport models. An important class among them are Reynolds stress transport models, which are based on transport PDEs derived from the original Navier-Stokes equations and take the following form [4]: where T (τ ) , P (τ ) , E (τ ) , and Π (τ ) indicate flux (due to diffusion and triple velocity correlation), production (energy extraction from mean flow), turbulence dissipation (destruction due to viscosity), and pressurestrain-rate (energy exchange among different components of the Reynolds stress tensor), respectively. We refer to the Reynolds stress transport equation as the constitutive (closure) PDE for the RANS equation (the primary equation). Unfortunately, the Reynolds stress transport equation itself needs to be closed by using covariance terms of even higher order. The unclosed terms include velocity triple correlation, pressurestrain-rate, and dissipation, among others. In particular, the pressure-strain-rate term plays a critical role and has a large impact on the predictive performances of the Reynolds stress models, but it is notoriously difficult to model. Consequently, this difficulty impairs the robustness of Reynolds stress models and restricts their usage to a small fraction of practical turbulent flows despite their theoretical superiority [5,6].
The Reynolds stress transport equation (3) implies a nonlocal mapping from the primary variable (mean velocity u) to the closure variable (Reynolds stress τ ). Following Lumley [7,8], such a nonlocal mapping in the turbulent constitutive relation is illustrated in Fig. 1a (where the generic closure variable c should be interpreted as τ ), and the mapping can be written as follows: where x indicates spatial coordinates, ξ is the spatial offset from x, and R x is the region of dependence (influence region) for the closure variable τ (x). In contrast, the widely used eddy viscosity models can be loosely considered local since they are derived based on weak equilibrium assumption [1,9], where the anisotropy b of the Reynolds stress (i.e., its orientation and aspect ratio) depends only on the local strainrate, i.e., b(x) = F (u(x)), where b = − τ 2k + I 3 and I is the identity tensor. This is illustrated in Fig. 1b.
Note that the magnitude of the Reynolds stress (i.e., turbulent kinetic energy k) in such eddy viscosity models does contain upstream information, and as such the generic closure variable c in Fig. 1b should be interpreted as the Reynolds stress anisotropy b. This fundamental difference contributes to the theoretical superiority of the Reynolds stress transport models over eddy viscosity models mentioned above. The example of Reynolds stress transport models above highlights a number of intrinsic difficulties with the traditional PDE-based constitutive models. In particular, they typically contain layers of submodels embedded in the PDEs. Moreover, the data needed to develop and individually calibrate such submodels are rarely available. Finally, interactions and coupling among the closure terms are difficult to quantify. Some of these difficulties are shared by other turbulence models such as the two-equation turbulence models [10] and the laminar-turbulent transition models based on transport equations, as well as continuum-based modeling of granular flows [11]. For example, both classes of the commonly used laminar-turbulent transition models, the intermittency transport model [2,12] and the amplification factor transport (AFT) model [3,13], are based on PDEs of the following form: where φ indicates the closure variable (intermittency γ or amplification factorñ) in the respective models; , and E (φ) denote diffusion transport, production, and dissipation/destruction (for intermittency transport model only) for the corresponding primary variable φ. In particular, the production term in both models are closed by layers of submodels with empirical calibrations.

Neural networks for nonlocal constitutive modeling
In view of the common difficulty of such PDE-based constitutive models in transport-dominated physics, we propose a neural network that can represent the nonlocal, region-to-point mapping. The motivations are twofold: (1) to resolve the conflict between the need for calibrating many submodels and the lack of data at such levels, and (2) to improve the robustness of the resulted constitutive models. Our data-driven approach to nonlocal constitutive modeling draws inspiration from the seminal work of Lumley [7]. He pointed out that turbulence in homogeneous shear flow at high Reynolds number behaves like a viscoelastic solid, based on which he argued that a universal constitutive relation exists for turbulence far removed from boundaries (e.g., wakes, jets, and the outer boundary layer in wall-bounded flows). A turbulence constitutive model is defined as a functional mapping from mean velocities to the Reynolds stresses [7]. The dependence can be nonlocal, with the nonlocality (size of the influence region) determined from the spatial scales of the physics [see e.g. , 14]. However, by definition, it is only determined by the local structure of the turbulence and should not depend on initial or boundary conditions (except for a level of time and length scales being set therefrom).
In the present work, we study a broader class of nonlocal constitutive models mapping the mean velocity to the closure variable c in the context of data-driven modeling. The details of the nonlocal models and the role of the closure variable c are explained in Section 2.1. Such closure models are typically governed by the underlying transport processes. In some cases the exact physics may be written in explicit transport PDEs (e.g., the Reynolds stress transport equations) while in other cases an exact PDE is not available (e.g., intermittency or amplification factor in transition) and is heuristically motivated instead. For brevity, hereafter the nonlocal mapping in Eq. (4) above is written in a more compact form for a general closure variable c as where u indicates the velocities at all points in the region R x (the shaded region in Fig. 1a). In these models, the closure variable c at location x depends on u in the region R x , which may be spatially varying depending on the local physics. Typically, the influence diminishes for points further away from x (i.e., as the offset distance ξ increases), which is indicated by gradually fading shade towards the edge of the shaded region. In contrast, classical constitutive models mostly involve local, point-to-point mapping as shown in We propose representing the region-to-point mapping in nonlocal constitutive models (as depicted in Fig. 1 and Eq. (4)) by using a neural network. The structure of the neural network is chosen to mimic that of the underlying transport PDEs in the classical constitutive models, e.g., Eqs. (3) and (5) for Reynolds stress models and transition models, respectively. Specifically, the neural network consists of a convolution between a kernel and a production term. The kernel is a nonlocal function of the velocities on the influence region, while the production term is a local function of the velocity field u and the closure variable c. The development and application of such a neural network-based constitutive model thus consists of the following steps (taking turbulence modeling as an example) as illustrated in Fig. 2: (a) perform direct simulation, experiments, or inferences to generate training data by extracting the regionto-point mapping from representative flow configurations and regions therein (Fig. 2a), (b) train the neural network of region-to-point mapping u → c by using the generated data (Fig. 2b); once trained, the neural network model replaces the PDEs traditionally used to describe nonlocal constitutive relations (e.g., the Reynolds stress transport equation (3), where closure variable c is Reynolds stress τ ); (c) use the trained constitutive neural network c = F NN [ u ] as the closure to the primary PDEs, e.g., the RANS momentum equation in turbulence modeling to solve for mean velocities and pressure (Fig. 2c).
The constitutive neural network proposed here is distinct from the local, point-to-point data-driven constitutive models [e.g., 15-18, in turbulence modeling]. Clearly, the purely local models are not able to reflect the nonlocal physics implied in many constitutive models as outlined above. Recently, field-tofield mappings or global operators parameterized by machine learning models for approximating solutions of the PDEs [see, e.g., [19][20][21][22][23][24][25][26][27][28] have emerged as a new computation tool. Such tools can speed up outerloop applications such as uncertainty quantification [29][30][31] and design optimization [32][33][34], and also hold the promise to serve as nonlocal constitutive modelling. However, some trained global mapping therein has initial conditions embedded in the network and thus cannot serve as constitutive models. In contrast to most work cited above, this work focuses on a region-to-point mapping better suiting the needs for representing nonlocal constitutive models instead of finding global PDE solutions with machine learning models.
In this work we will focus on steps (a) and (b) outlined in Figs. 2a and 2b, respectively, i.e., generating data and training a nonlocal constitutive model, while the integration with the primary solver is left for future work. Specifically, we demonstrate the concept of the constitutive neural network on a convectiondiffusion-reaction equation, which is a hypothetical constitutive model yet a prototype for many actual constitutive models in computational mechanics. Moreover, the velocity field in the convection-diffusionreaction equation is specified and frozen, i.e., it is not affected by the closure variable c. The objective is to demonstrate that the proposed network architecture representing a nonlocal, region-to-point mapping, is a viable approach to describe nonlocal constitutive relations typically depicted by PDEs. Compared to such ad hoc and often heuristically justified closure PDEs, the neural network-based constitutive model is easier to train and has more flexibility and expressive power. We expect this preliminary study shall pave way for full-fledged data-driven turbulence models [35] and data-driven constitutive models in other related fields, such as computational mechanics [36][37][38][39][40] and non-Newtonian fluid dynamics [41]. Finally, we point out that the framework lacks rotational invariance, and this deficiency is addressed in our follow-on work [42].

Contributions of present work
The present work aims to demonstrate the merits of utilizing the underlying mathematical structure and physics to guide the design of neural network models. Even though the motivation example (linear differential equation) for the design of the network structure is generally invalid for intended practical applications with nonlinearity, the machine learning framework is much more disciplined than a purely black-box model, The neural network architecture in panel (b) is inspired by the structure of the formal solution to transport PDEs (using Reynolds stress transport equation (3) as an example), which is written as a convolution G P (u) and represents nonlocal constitutive relations. A convolutional network is thus constructed to represent the nonlocal mapping from velocity field to the kernel G, and a fully connected network is used to model the production term P as a local, pointwise function of u (or, more precisely, the strain rate S and rotation rate Ω). The framework lacks rotational invariance, and this deficiency is addressed in our follow-on work [42].
and it performs very well in the nonlinear problems. The same spirit has been used in Refs. [43,44] in inverse problems like tomographic reconstruction and wave scattering inversion. The present work focuses on transport constitutive relations that have a convection-diffusion-reaction PDE structure. Such problems are commonly found in computational physics and mechanics, particularly in fluid dynamics.
Neural networks have been used to solve PDEs [45][46][47][48][49] and to build surrogate models for the solutions of a class of parameterized PDEs [22,50]. In such works the intent was to solve clearly specified PDEs by using the solutions of similar PDEs as training data. In contrast, the present work aims to learn from data a nonlocal mapping implied by an unknown PDE, albeit with its structure assumed to be known (convectiondiffusion-reaction type). In physical modeling it is common to know only the structure and not the specific form of a PDE, particularly when modeling unresolved processes. For example, the widely used Spalart-Allmaras turbulence model is based on a transport PDE of a viscosity-like variableν without rigorous derivation [51], so are the intermittency and amplification factor transport equations in laminar-turbulent transition modeling [2,3]. Moreover, even though the turbulence dissipation needed for two-equation and Reynolds stress models can be derived exactly, an ad hoc form is often used [1,52]. As such, in this work we aim to learn a region-to-point mapping that embodies the physics described by the underlying (possibly unknown) PDE without requiring data at the submodel levels (e.g., those for the pressure-strain-rate, triple correlation in Reynolds stress model). Given the context of the present work in constitutive modeling, which by definition shall be free from the boundary and initial conditions, we consider only the physics governed by PDEs. This is consistent with the comment made in Lumley [7], who, through analogy to St Venant's principle in solid mechanics, stated that it is unlikely for a turbulent constitutive relation to exist near the boundaries. This insight applies not only to turbulence models but also to many other nonlocal constitutive relations. In practice, turbulence modeling within RANS or LES (large eddy simulation) solvers, the boundary regions near solid walls (fluid-solid interfaces) or free surfaces (liquid-gas interfaces) are handled by separate yet complementary closure models, i.e., wall models [53] and free surface models [54], respectively.
The rest of the paper is organized as follows. Section 2 describes the problem the neural network-based constitutive modeling aims to solve and presents the technical details of the proposed framework, including the PDE-inspired neural network architecture, the choice of region size in the nonlocal mapping, and the details of data generation and model training. Section 3 presents detailed verification of the proposed method, followed by multiple test cases to demonstrate the predictive capability of the constitutive neural network.
Finally, Section 4 concludes the paper.

Problem statement
We consider a steady-state convection-diffusion-reaction PDE describing the transport of an active tracer u · ∇c − ∇ · ν∇c + E = P(c; u) where c is the concentration of the tracer, which is the closure variable (analogous to Reynolds stress τ in RANS equations), u is the given flow field (whose governing equation contains unclosed term involving c), ν is the diffusion coefficient, P is a production term depending only locally on the flow velocity and concentration, and E contains the dissipation term and any sources that are nonlocal functions of u and c. The active tracer concentration variable c corresponds to unclosed terms such as Reynolds stresses, intermittency, or amplification factor in the examples given in Section 1. In this work we consider the scenario where ν is a specified constant, E = ζc with ζ > 0, and P is a function of the flow field and concentration: where S =S/S max , Ω =Ω/Ω 0 ,S = ∇u+(∇u) andΩ = ∇u−(∇u) are the magnitudes, of the strain rate and rotation rate tensors, respectively, with · denoting tensor norm;S max andΩ 0 are normalization constants. The functions are defined as: The behaviors of function f (S) and h(c) are presented in Fig. 3. The function h(c) has three distinct regimes, which is representative of the closure submodels in fluid dynamics. In particular, h(c) has a transition regime (shaded in Fig. 3b) with a steep gradient, which can pose challenges for constitutive modeling. A similar submodel can be found in the particle-laden flow simulations, where the drag on a spherical particle as a function of Reynolds number can be primarily divided into laminar and turbulent flow regimes, with a transition regime in between [55].
The overall problem statement above is reminiscent of the closure models in fluid dynamics, particularly the laminar-turbulence transition equation in Eq. (5). Specifically, while the whole ad hoc transport PDE serves as a constitutive model in a mean flow solver, the equation itself contains production and destruction terms that require closure (with several layers of embedded, regime-dependent models as amplified in the setup of Eq. (9c)). The Reynolds transport model, on the other hand, has additional difficulties involving inter-component coupling and energy transfer through pressure-strain-rate term Π and wall models, which are not considered in present work. As such, the present work can be seen as a first attempt towards a data-driven Reynolds stress transport model.

Proposed network architecture
The objective of this work is to use a neural network as a surrogate model that represents the nonlocal mapping described by the convection-diffusion-reaction PDEs. The mapping not only embodies the nonlocal dependence due to convection, diffusion, and other nonlocal production terms, but also contains the effects of the submodels, i.e., when the production is a nonlinear function of the flow field u and active tracer concentration c.
Our neural network architecture for nonlocal constitutive modeling (Fig. 2b) is inspired by the linear simplification of the convection-diffusion-reaction equation. To see the motivation, we write the transport equations in the constitutive models in an abstract form: where u is the flow field (determined from the primary solver for which the constitutive model serve as closure), c is the primary variable in the constitutive model of concern, N denotes the differential operators including convective and diffusive transport terms along with any nonlocal production terms, P is the production function containing local, pointwise dependence on c and u. Notice that the turbulence models Eq. (3), the transport model Eq. (5), and the prototype model Eq. (7) can all fit into the form of Eq. (10) above. The separation of local and nonlocal operators in Eq. (10) is often justified by the relevant physics.
For example, in Reynolds stress transport equation (3), the production, dissipation, and velocity triple correlation terms are all local, while the pressure-strain-rate is nonlocal due to the elliptic nature of the pressure equation.
If Eq. (10) is linear in c, we can denote the linear differential operator associated with the flow field u by L u and assume the production only depends on u, i.e., : which is defined on an infinite domain R d with vanishing boundary condition, then the solution can be written as [see, e.g., 56]: Here G u is the Green's function corresponding to, specifically the inverse operator of, the linear differential operator L u , P = P(u(x )), and denotes the convolution operator. Motivated by the structure of Eq. (12), we propose the following network architecture to represent the general nonlocal constitutive relation.
Specifically, we intend to learn a nonlocal mapping from a patch of flow field u to the concentration variable c at the point of interest x. In this work, we focus on two-dimensional (2D) problems and assume the field is discretized on a regular grid. Suppose the point of interest x corresponds to the grid index (i, j) and we use a stencil of size (2k + 1) × (2k + 1) (k ∈ Z + ) to approximate the velocity field in the considered influence region: where is the offset index set with stride n 0 ∈ Z + (positive integer). Here k (along with the discretization grid) implies the size of stencil, and the parameter n 0 determines how frequent the data is sampled in the region (patch). If n 0 > 1, the velocity field in the considered region is downsampled. A few examples of combinations of k and n 0 are shown in Fig. 4 to illustrate how they determine the region size and the spacing between stencil points.
The proposed architecture consists of two distinct modules as shown in Fig. 2b: (1) a convolutional network to learn the mapping u i,j → G i,j from the patch of the velocity u i,j surrounding a point to the Green's kernel; (2) a fully connected neural network to learn the local mapping u i,j → P i,j from the velocity at point (i, j) to the production term at the same point. In order to guarantee the universality of the constitutive relation, the mapping represented by the fully connected network is shared by each point.
Finally, with these two modules above, the final tracer concentration c is computed as a discrete convolution of G and P as follows: This is a discrete counterpart of the formal solution in Eq. (12) for the linear convection-diffusion-reaction equation. We remark that the structure in Eqs. (12) and (15) can be generalized to three dimensions (3D) straightforwardly such that our learning approach is capable of dealing with 3D problems in the same way.
Note that no data at the submodel level, e.g., the Green's kernel G or the production term P will be used, since such data are rarely available in our application scenarios. In a linear problem (Eq. (11)), the ground truth of both G and P are known with clear interpretation. It is expected that both mappings, u i,j → G i,j and u i,j → P i,j , can be learned approximately (up to a multiplication constant among G and P ) given the data. This will be verified in two preliminary prediction cases in Section 3.1 (Case I ) and Section 3.2 (Case II ). In Case I, we consider a steady-state convection-diffusion-reaction PDE describing the transport of an active tracer u · ∇c − ∇ · ν∇c + E = P where P is a random production term generated by Gaussian process. We will verify that the Green's kernel can be learned approximately. In Case II, we consider a similar PDE where P is a production term depending only locally on the flow velocity. We will verify that the submodel on P can be learned approximately in two different situations: For the nonlinear problem (Eq. (7)) of our interest, the exact structure in Eq. (12) no longer holds, but we expect that the network structure inspired by the linear problem, together with the nonlinearity of the networks, provide sufficient approximation capability for the resulted nonlocal mapping. This will be verified in Section 3.3. We remarked that while the production term is chosen to be a product of two or three univariate functions for convenience, this fact is never used in the data-driven modeling, and thus all the conclusions should be equally valid for general production models.

Form of training data and region size
As introduced above, the nonlocal, region-to-point mapping is trained by data consisting of pairs of ( u , c) (or ( u , P , c) in the special Case I ). When preparing the training data u , we need to specify k and n 0 in Eq. (14), which essentially dictates the size of the region considered and the downsampling level therein. The necessary size of the region should be determined by the physical influence (depicted in Fig. 2), which is in turn determined by the underlying PDE, while the downsampling level should be determined by the budget of points to be used. The question of determining the patch size amounts to finding the minimum region size of a truncated domain R x such that the true solution can be approximated well by the where | · | denote absolute value. We have derived the region size in a simplified setting with details provided in Appendix A and only a summary is provided below.
We consider a steady-state convection-diffusion-reaction equation with constant coefficients. The diffusion and dissipation constants are denoted by ν and ζ. Suppose we use a square region as the influence region with error tolerance . The half-length of its side used for both the upstream and downstream direction should be approximately = 2ν log It can be seen that the region size 2 increases with larger diffusion coefficient ν and velocity magnitude u but decreases with larger dissipation coefficient ζ. In fact, the diffusion and dissipation coefficients control the local equilibrium. The larger ζ is or the smaller ν is, the more local the problem becomes. We remark that in practice an accurate size of the influence region can be spatially inhomogeneous and difficult to obtain. Our derivation provides a simple formula to apply in general situations and it can serve as a rule-ofthumb guideline for practice. In the derivation, we have assumed the truncated region is square in order to employ the analytical solution in 1D. Moreover, we have calculated the size conservatively by treating the upstream and downstream regions equally. The difference between the streamwise and lateral velocities is ignored as well. Although the derivation is based on these simplifications, we will see in Section 3.1 that the best learning performance is achieved when stencil size parameter k and sampling density parameter n 0 are properly chosen such that the underlying region size matches that suggested by the theoretical derivation.

Implementation
The proposed neural network is implemented in Python by using the machine learning library Py-Torch [57]. The training data are generated by solving the underlying convection-diffusion-reaction PDEs with FiPy [58], a finite volume PDE solver. The code for training data generation and neural network training are both made available in a public GitHub repository [59], and thus the results presented here can be straightforwardly reproduced and further developed.

Results
In this section we evaluate the performance of the proposed method in three cases, with underlying flow fields being uniform flows (constant velocity in the entire domain), plane Poiseuille flows, and Taylor-Green vortexes. First, in Section 3.1 we use the uniform flow with randomly sampled production field P(x) to verify that the proposed network structure is consistent with the solution structure for linear PDEs and that the proposed algorithm can learn the Green's kernel well. We will also demonstrate that the theoretically derived region size is indeed near-optimal numerically. As a further verification, in Section 3.2 we learn from data generated from Eq. (16) with flow-dependent production P(u(x)). We demonstrate that the proposed architecture is capable of learning the underlying submodel for the production term, i.e., its dependence on the flow field. Finally, the setup in Section 3.3 showcases the predictive capability of the proposed framework for nonlinear PDE in which the production term is a function of both the flow field and concentration as described in Eq. (8).
In each numerical experiment, we draw 100 random samples of physical coefficients to construct 100 PDEs and solve the PDEs to generate training data. The diffusion coefficient ν and dissipation coefficient ζ are constant in different PDEs belonging to the same case. The domain of the flow field is [0, 1] × [0, 1] with periodic boundary condition and discretized with grid size 0.01. In total we use 500,000 pairs of ( u , c) (or ( u , P , c) in Case I ) as our training data and reserve 50,000 data points for testing. The details of network architectures and training parameters are provided in Appendix B. We evaluate the prediction performance of the trained neural network by the normalized prediction error, which is defined as: where the summation is performed over all of the N s training or prediction data points.

Case I: Verification of learning Green's kernel
In this case, we generate training data by solving the transport PDE (16) with uniform velocity fields and random production fields P(x). The diffusion coefficient and dissipation coefficient in the PDE are set as ν = 0.05 and ζ = 12, respectively. The velocities are in the x-direction and uniformly distributed between 0 and 2.0. The random production fields P(x) are generated by sampling from a zero-mean Gaussian process (also called Gaussian random field) with the periodic covariance kernel [60,61], i.e., P ∼ GP(0, The objective of this verification is to demonstrate the capability of the proposed neural architecture in recovering the mathematical structure of the analytical solution that is used to motivate it. Specifically, the training data imply a region-to-point mapping ( u i,j , P i,j ) → c i,j , while the underlying exact mapping is c = G u P . In general, the Green's kernel depends on the flow field u. In this case, since the velocity field is uniform, the analytical kernel G becomes constant that can be easily compared with the learning result.
We remark that in this special case, the discretized constant kernel can be interpreted as the inverse of the finite difference approximation of the operator L[c] ≡ (u · ∇ + ∇ · ν∇ + ζ)[c] (see Appendix A for more details). For the region-to-point mapping in this case, by default we use 11 × 11 stencil points and distance between the adjacent stencil points is 0.06 (downsampling parameter n 0 = 6) such that the size of region covered by these points is approximately 0.6 × 0.6, which is close to the suggested region size 2 determined by ν, ζ, u and error tolerance = 0.2 according to Eq. (21). A comparison of the predicted tracer concentration field and the ground truth is presented in Fig. 6. We can see in Figs. 6a and 6b that the predicted concentration field and the ground truth have good overall agreement. This can be seen more clearly in the cross section plots in Fig. 6c and 6d where the predicted and the ground truth concentrations have rather good overall agreement except for some minor discrepancy.
Note that the dramatic difference between the original production and the resulting concentration reveals that they are related via a non-trivial mapping due to the convection and diffusion effect of the flow. The learning of this mapping from data is the subject of our discussion below.
It is noteworthy that the neural network not only predicts the closure variable c very well, but also achieves it for the right reason: it learned the Green's kernel correctly, as demonstrated in Fig. 7. In this case, we know the analytic form of the Green's function and we take its discretization with the corresponding grid size as the analytical kernel. The learned kernel G is of size 11×11, and we compare it with the analytical kernel. Since the velocity u is in the x-direction, the kernel is essentially 1D, depending only on the relative coordinate in the x-direction. Accordingly, we compute the average of the learned kernel along the y-direction and obtain the weight in the x-direction. The 1D weights of both the learned kernel and theoretical kernel are normalized such that the maximum is 1.0. We can see that the trend of learned G is quite similar to that of the analytical G in the x-direction, explaining the good performance of the concentration prediction.
Notice that G determines the weight of the production at different points relative to the center point. Since the productions in the upstream have a greater influence on the center point than those in the downstream, the weights on the left (upstream) side are larger than the weights of the symmetric positions on the right (downstream) side, and our learned kernel reveals such physics correctly. Slight differences of the weights are observed in the left-most (upstream) and right-most (downstream) positions, which may be caused by the padding used in our convolutional neural networks. However, the results are not sensitive to the padding mode that is used (zeros, reflect or replicate padding). It is also possible that the slight difference is caused by the cutoff of the influence region and downsampling procedure.  Fig. 4. A total number of 20 cases with various combinations of region size and stencil size are investigated, and the prediction error of concentration in each case are presented in Fig. 8.
From Fig. 8 we can see that, with each fixed stencil size (2k + 1) × (2k + 1), the smallest prediction error is achieved when the region size is close to the theoretical optimum, i.e., D/(2 ) ≈ 1. This observation suggests that, given the same number of data points as the computational budget on the input size, the learning is most effective when the information is drawn from the proper physical domain of influence. Otherwise, when the region size D/(2 ) departs significantly from the theoretical optimum, e.g., when D/(2 ) approach the lower limit 0.5 or the higher limit 1.5, the prediction error increases dramatically. This is because when D/(2 ) is much smaller than 1, the stencil fails to account for the points whose velocity u has significant influences on the concentration c. On the other hand, when the region size D/(2 ) is much larger than 1, the information is diluted by using many points whose velocity has little influence on the concentration c.

Case II: Verification of learning submodel for PDEs
In this case, we consider two families of flow fields shown in Fig. 9. The first is plane Poiseuille flows, whose velocity field u = (u, v) is where α is the pressure gradient parameter (i.e., α = − 1 2µ dp dx , with µ being the dynamic viscosity). Another family of flows is Taylor-Green vortexes [62], whose velocity filed u = (u, v) is where U 0 is the maximum velocity. For the flows above, we consider three scenarios of specified production fields: Case IIa The production is only a function of the normalized strain rate magnitude S, i.e., P(u) = f (S), with f (S) defined in Eq. (9a) and depicted in Fig. 3a,
This is a more challenging case that requires a more comprehensive investigation. It is omitted here for brevity (see Appendix C for details).
Given any combination of a velocity field and a production field generated as detailed above, we then generate the training data ( u, c ) by solving the transport PDE (17). The diffusion coefficient in the PDE is set to ν = 0.02, and dissipation coefficient is ζ = 20. The stencil size for the region-to-point mapping in both families of flows is set to be 11 × 11, and the downsampling parameter n 0 is 2 such that the region size is approximately 0.2 × 0.2, which is very close to the suggested region size 2 determined by ν, ζ, u and error tolerance = 0.2 according to Eq. (21).
The training data define a region-to-point mapping ( u i,j ) → c i,j , while the underlying exact mapping is c = G u P(S) or c = G u P(S, Ω). The production is a function of the strain rate (and rotation rate), so both the Green's kernel and production depend on the flow field u. Apart from learning the Green's kernel, the constitutive neural network also correctly learns the production term P . In fact, here we will demonstrate that the proposed network can even learn the production submodel P(S, Ω) itself, i.e., the functional dependence of production P on the strain rate magnitude S and rotation rate Ω.
In Case IIb, i.e., with production P(u) = f (S)g(Ω), the prediction error plots (omitted here for brevity; see Appendix D) show excellent agreement between prediction and ground truth of the tracer concentration c. The normalized prediction errors for the Poiseuille flow and Taylor-Green vortex are 0.62% and 1.05%, respectively. A comparison of the predicted tracer concentration field and the ground truth for a representative test case (out of a total of five) is shown in Fig. 10, which suggests that the concentration fields are predicted very well. Compared to Case I with a uniform velocity field, the predicted concentrations for both plane Poiseuille flow and Taylor-Green vortex are closer to the ground truths and much smoother, even though the flows are more complex. This is likely because the production field is smoother (compared to the production field drawn from the Gaussian processes). The network not only predicts the closure variable c accurately, but also learns the underlying production. This is demonstrated in Fig. 11. As the solution c(u) = G(u) P (u) is a convolution of the Green's kernel and the production term, we can only expect to learn both G and P up to a scaling constant. This is because the solution remains the same with respect to an arbitrary scaling β to P accompanied by a scaling of 1/β to G. Therefore, before the comparison we first multiplied the learned production fieldP by a scaling constantβ, i.e., P =βP . The scaling factorβ is chosen such that the squared discrepancy between the learned productionP and the ground truth P * is minimized, i.e.,β = arg min β βP − P * 2 . With such a scaling, it can be seen that the predicted production fields in Fig. 11 show similar patterns to that of the ground truths for both the plane Poiseuille flow and the Taylor-Green vortex. Given that the neural network predicts well the pattern of the production field, we further verify that the network is able to learn the underlying submodel of the production. In Case IIa, where P(u) = f (S), we can see from Even in the more challenging Case IIb, where P(u) = f (S)g(Ω), the neural network is still able to learn the submodel of the production term, which is demonstrated in Fig. 13. However, for the plane Poiseuille flow, a caveat is that the strain rate and the rotation rate are co-linear. That is,S =Ω = | ∂u ∂y | as ∂v ∂x = 0, and thus S and Ω only differ by a factor of the normalization constant. This coincidental co-linearity is due to the simplicity of the flow field (homogeneity in the stream direction). The co-linearity can also be confirmed from the scatter plot in Fig. 14 (for Case III ), as the projection of the points for the Poiseuille flow onto the S-Ω plane would form a straight line (not shown). As a result, the bivariate production function P(S, Ω) degenerates to a univariate function of S only. In light of this observation, the neural network successfully learned the functional relationship between the production P(S, Ω) = f (S)g(Ω) and the strain rate S (see Fig. 13a). For Taylor-Green vortex, the correlation between the strain rate and the rotation rate is much weaker, and thus the production P(S, Ω) is determined jointly by both strain rate S and rotation rate Ω (both normalized properly). Consequently, the neural network is able to learn both f (S) and g(Ω) from the data, as shown in Fig. 13b and 13c, respectively. We note that, in all cases, there are discrepancies between the learned productions near the upper and lower bounds of the strain rate S and rotation rate Ω. This is mainly caused by the uneven distribution of strain rate S and rotation rate Ω in the training data. Specifically, in plane Poiseuille flows, S = 1 4 | ∂u ∂y | = α 2 |y − 0.5|, where α is uniformly sampled in the interval between 1 and 4. In this way, the probability of S being picked close to is 1 is relatively small only when α ≈ 4 and y ≈ 0 or 1, which leads to the uneven data distribution. A similar situation happens to Taylor-Green vortexes in terms of both strain rate S and rotation rate Ω.
In Appendix C, we study further a scenario (Case IIc) showing the effect of the data distribution on the learning of submodel, especially when the data in the regime of interest is not well represented.   and 2.69%, respectively. A comparison of the predicted tracer concentration field and the ground truth for a representative case is shown in Fig. 15. For both flows, the predicted concentration field is similar to its corresponding ground truth, although there are slight differences seen: for plane Poiseuille flows, the region with low concentrations of prediction is larger than that of the ground truth; for Taylor-Green vortexes, truth for a representative extrapolation case is shown in Fig. 16. For both flows, we can still see a similar pattern of concentrations between the prediction and ground truth, albeit the distinctions near the boundary and the centerline.

Conclusion
Many constitutive relations in computational mechanics, particularly those involving fluid flows and transport phenomena, are formulated as partial differential equations (PDEs) that describe the nonlocal mapping from primary variables to closure variables. These PDEs are often heuristically motivated and contain many submodels, which pose challenges for calibration or even for obtaining necessary data for their development. In this work we propose using neural networks to represent the region-to-point mappings in nonlocal constitutive relations. Such a representation leads to a constitutive neural network with a simple structure, which is inspired by the solution to linear convection-diffusion-reaction equations and is much easier to train than traditional models. We verify that the proposed network can learn the Green's kernel from the data generated by underlying linear PDEs. We further demonstrate that the data can be used to learn the submodels of a wide range of complexities without explicitly using the data from the latter level.
Finally, we show the model has good predictive capability for general constitutive relations in data generated from nonlinear PDEs.
Despite the preliminary successes demonstrated here, we point out several directions of future study following the present work. First and foremost, the performance of the learned closure model in solving the primary PDEs (e.g., the RANS momentum equation in turbulence modeling to solve for mean velocities and pressure) should be systematically examined in future work. To this end, how to integrate the learned constitutive relation with other boundary conditions when their effects are nonnegligible is certainly worthy of future exploration. Moreover, we assume a regular mesh to present the main idea of our data-driven approach to nonlocal constitutive modeling. Future work concerns extending the approach to unstructured meshes to facilitate the usage on more practical data coming from large-scale industrial CFD simulations and experiments. Our follow-on work [42] has shown a lot of promise in this direction. Finally, more efficient and adaptive sampling of data should be studied to boost the learning of more complicated flow fields.
The methodology developed here is demonstrated on 2D problems only. It can be extended to 3D quite straightforwardly and its numerical performance should be validated as well.
Appendix A. Analysis of region size in 1D transport equation

. Analytic solution
The steady-state convection-diffusion-reaction equation in 1D with constant coefficients is written as: with boundary condition c(a) = c(b) = 0. Here ν, u, ζ denote diffusion, flow velocity and dissipation constants, respectively. ν, ζ are both positive, and without loss of generality, we assume u ≥ 0. If we define the corresponding differential operator L as the associated Green's function G(x; x ), with a fixed x ∈ [a, b] satisfies LG(x; with G(a; x ) = G(b; x ) = 0. Here δ(x) denotes the Dirac delta function. With the Green's function, the solution to (A.1) can be written as in the same form as the first equality of (12).
To eliminate the effect of boundary conditions and facilitate the analysis in the meantime, we further assume a = −∞, b = ∞. In this case, it can be straightforwardly verified that where λ 1 , λ 2 are two roots of the associated characteristic polynomial −νy 2 + uy + ζ = 0: The solution to (A.1) becomes Notice that it is exactly the convolution form presented in Eq. (12), indicating the nonlocal dependence on the source. Specifically, recalling that we have assumed u ≥ 0, so we can interpret the two terms in can see that the influence region increases with larger diffusion coefficient ν and velocity magnitude u but decreases with larger dissipation coefficient ζ. The quantitative relation between the influence region and these parameters will be derived below. In the discrete sense, the Green's function can also be recovered from the finite difference method.
Suppose we have a uniform mesh grid in 1D: x 0 < x 1 < · · · < x N and denote the concentration and source at these grid points by c i and P i , respectively. The discretization (with upwind and central scheme for the convection and diffusion terms and zero boundary condition) of Eq. (A.1) leads to the following system of linear equations: where coefficients a = u/∆x and d = ν/∆x 2 come from the discretization of the advection and the diffusion terms, respectively. Denote the coefficient matrix in (A.8) with , then the Green's kernel is closely related to L −1 . In particular, the solution indicates how the production terms at different location points contribute to the concentration through the weights L −1 . Indexing these weights according to the relative coordinates, we can see that they approximately equal to the discrete value of the Green's function, as shown in Fig To get more qualitative insight of the relationship between the size of the truncated region and the corresponding accuracy, we consider a simplified setting. Suppose the production term P (x) varies a little such that it is roughly a constant. In this situation, in order to achieve a relative error tolerance for the contribution in the upstream region (the first term in Eq. (A.10)), we require Similarly, the analysis of the downstream region gives us Noticing that |λ 2 | < |λ 1 |, we see the influence region of the upstream is larger than that of the downstream.
However, in practical 2D and 3D applications, we may not know the direction of flow in advance, so for simplicity we shall choose a conservative patch size 2 − to guarantee the small error. Furthermore, in the 2D case, we can assume the truncated domain is square such that we can reuse the result in the 1D case, in which u is replaced by the magnitude u as a conservative choice. Taking the above simplifications into consideration, we arrive at the final half-length of the patch = 2ν log (A.14) Appendix B. Details of network architecture and training parameters The details of the employed network architecture are listed in Table B.1. The entire neural network consists of two sub-networks: (1) a convolutional neural network and (2) a fully connected network.
The convolutional neural network learns a nonlinear region-to-region mapping u i,j → G i,j . It consists of six convolutional layers without pooling layers. The number of channels in each layer is chosen to be the powers of two empirically. The convolutional kernels in all layers are of size 3 × 3, which proves to be very effective in image recognition [63,64]. The amount of zero padding is set as one to ensure the output and input of each layer is of the same size. Other padding modes including reflect and replicate padding modes are also tested when learning the Green's kernel in Cass II, with only slight differences. Batch normalization [65] is used in each layer except the last layer to stabilize the learning process and accelerate the training.
The fully connected network learns a pointwise mapping u i,j → P i,j , applicable for every stencil point in the influence region u i,j . In fact, this mapping consists of two steps: u i,j → (S i,j , Ω i,j ) and (S i,j , Ω i,j ) → P i,j . The first step calculates the physical properties (strain rate and rotation rate) of the flow without any trainable parameters. The second step corresponds to this fully connected network. We set the number of layers as five by default. When the source function is P = P(c; u) (Case III ) or P = h(S) (Appendix C), we add one more layer (32 × 32) between layer 1 and 2 to ensure enough complexity of the fully connected network. We select the rectified linear activation function (ReLU) as the activation function because it is easier to train and often achieves good performance.
We adopt the Adam optimizer [66] to optimize the parameters in neural networks. In all the cases the training takes 500 epochs with constant learning rate 0.001. The only exception is Case I where the random production generated by the Gaussian process makes it harder to train the neural network. In Case I, the training takes 2000 epochs, and the learning rate is 0.001 at the beginning and then multiplied by 0.7 every 600 epochs. Learning submodels (e.g., P(S)) from closure variable data (e.g., { u , c}) is of great practical importance as closure models in computational physics often have separate embedded models to account for different unresolved physical processes. While we have demonstrated that the proposed framework has successfully learned univariate production submodel P(S) = 4 sin (2πS) + 6S 2 + 5e S , we further evaluate our framework in learning a more challenging production model. Specifically, in this assessment the submodel P(S) has three distinct regimes with steep gradient in the middle regime as shown in Eq. (9c) (and depicted in Fig. 12): The overall trend of the submodel function h(S) can be learned correctly, as shown in Fig. C.19a.
However, the learned model failed to capture the sharp transition in the region S ∈ [0.5, 0.55]. Rather, the learned function has a relatively mild transition connecting the descending and ascending regimes. This is caused by the distribution of the training data used to learn the function. Specifically, we generated the training data under the plane Poiseuille flows, in which the velocity along the central axis is fixed and the pressure gradient parameter α is sampled from a uniform distribution on the interval [1,4]. Therefore, the distribution of velocity gradient (and thus strain rate magnitude S) is not explicitly controlled but implicitly determined by the flow field. The probability distribution of S is presented as overlaying shades in Fig. C.19.
It can be seen that the percentage of the data in the transition regime of h(S) is rather small (31,380 or 6.3% of the 500,000 data points in total) in this regime. As a result, they are not sufficiently weighted in the training process for learning the steep gradient in the middle regime. The learning of submodels can be improved by augmenting the data in the region of interest (the middle regime S ∈ [0.5, 0.55] in this example). In some applications this is well justified. For example, consider the submodel C d (Re) of drag coefficient C d on a spherical particle as a function of the particle Reynolds number Re. The range of Re in the transition regime (which connects the laminar regime at low Re and the turbulent regime at high Re) is known a priori from experimental evidence and the physics of transition. Hence, the augmentation of data in this regime can be performed. In fact, such augmentation shares a similar spirit with stratified sampling and importance sampling, which are commonly used techniques in statistics and machine learning [67][68][69][70][71]. After we augment the training data by increasing the number of points in the transition regime (which could have been achieved equally well by increasing sample weights in this region instead), the learned function agrees much better in the transition regime, as shown in Fig. C.19b. Unfortunately, the agreement in the high strain-rate regime (S ∈ [0.7, 1.0]) deteriorated as compared to Fig. C.19a, as the relative weights in this regime are too small. When we increase the percentage of samples in this regime in the training data, the submodel learning can be improved. The best result in this experiment in terms of the learned submodel is shown in Fig. C.19c, which is obtained by the data distribution shown in the corresponding background shade.
In summary, the results here demonstrate that the learning of the submodel has achieved overall success, but the performance is affected by the data distribution of the submodel variable S. If the overall behavior of the submodel is known, we can use such prior knowledge to guide the sampling in the training data and thus to improve the learning of the submodel. This idea is consistent with the well-known importance sampling techniques in statistics. Otherwise, some adaptive sampling techniques [67][68][69][70][71] may also help to improve the learning performance.

Appendix D. Additional plots in results
The prediction error plots for Cases IIb and III are presented in Figs. D.20 and D.21, respectively.