Notation f (x) Source terms g(x) Boundary values v(x) Velocity field ε Positive parameter Ω s Open domain ∂Ω s Boundary of Ω s R n n-dimensional Euclidean space δ Horizon parameter p(z, µ, Σ) Probability density function L nd δ Nonlocal operator χ 2 α (d) Chi-square distribution B nd δ,A,α (x) Truncated region L δ Nonlocal convection-diffusion operator L nd δ Nonlocal diffusion operator L nc δ Nonlocal convection operator γ nd Diffusion kernel function γ nc Diffusion kernel function I d d-dimensional identity matrix B nd δ,α Truncated influence region of diffusion term B nc δ,α Truncated influence region of convection term γ nd,α Truncated diffusion kernel function γ nc,α Truncated convection kernel function Ω bd Interaction domain about diffusion term Ω bc Interaction domain about convection term Ω c Total interaction domain L δ,α Truncated nonlocal convection-diffusion operator L nd δ,α Truncated nonlocal diffusion operator L nc δ,α Truncated nonlocal convection operator L 2 n0 (Ω) Constrained space B(•, •) Bilinear operator x i Discrete point S i Associated volume T h Grid set B δ (x) Euclidean balls B δ,h (x) Approximate balls △B δ,h (x) Ball difference L 2 n 0 (Ω) Constrained space

Nonlocal models, encompassing peridynamics models [3,19,27,28,30,35], nonlocal diffusion models [9,12,37,41,42], nonlocal advection [8,15,21], nonlocal convection-diffusion models [5,13,33,34], nonlocal Stokes equations [11] and so on, have garnered significant attention in recent decades from both theoretical and computational perspectives [1,16,14,20,25,26,39,43]. For a comprehensive exploration of nonlocal models, we refer to the monograph [17] and the survey article [6]. Following the establishment of the nonlocal vector calculus framework [10], modeling convection-diffusion problems has attracted considerable attention in recent years. A central challenge lies in identifying an appropriate kernel function and integrating convective information into it. A nonlocal, nonlinear advection model was firstly introduced in [8] by extending conventional pointwise concepts to account for nonlocal contributions to the flux. Then, [13] developed a nonlocal convection-diffusion models with volume-constrained boundary condition, while the wellposedness of such models was investigated in [4]. Furthermore, an upwind nonlocal model for convection-diffusion problems with a divergence-free velocity field was proposed in [33], the finite element discritization of which was proved to satisfy the discrete maximum principle. Subsequently, for general velocity fields, a conservative nonlocal convection-diffusion model with a specially constructed upwind convection term was introduced in [34], ensuring both mass conservation and adherence to the maximum principle. To solve nonlocal convection-dominated diffusion problems, Leng et al. [24] introduced an asymptotic compatible Petrov-Galerkin method. Based on the exsiting nonlocal convection-diffusion model, to simulate complicated processes involving chemical reactions,flows and diffusions, a bond-based peridynamic advection-reaction-diffusion model was formulated by Tian et al. [32].

In [36], an innovative Gaussian-type kernel-based nonlocal diffusion model was introduced, revolutionizing the simulation of diffusion processes by incorporating matrix-valued anisotropic coefficients in non-divergence form. This approach allows for the simulation of both isotropic and anisotropic diffusion phenomena by integrating the diffusion matrix into the covariance and establishing the kernel using a multivariate Gaussian distribution. Unlike existing nonlocal models, where the kernel function may be bounded, the Gaussian-type nonlocal model's kernel function is unbounded, enabling the integration of physics information directly into the kernel function. Building upon this approach, we propose and study a novel nonlocal convection-diffusion model corresponding to the following partial differential equation (PDE) model defined on an open domain Ω s ⊂ R d :

Here, ε > 0 serves as a positive parameter, which is relatively more compared to the above velocity field v(x). In this approach, the velocity field is incorporated into the kernel function, as an expectation of a multivariate Gaussian distribution. Considering that the support of the Gaussiantype kernel is unbounded, we strategically truncate the influence region to improve computational efficiency. Thus volume-constrained boundary condition is prescribed. The well-posedness of the proposed nonlocal model is successfully established, and the maximum principle is derived. With direct Riemann quadrature, [31], we further develop a meshfree discretization scheme which satisfies the discrete maximum principle and the asymptotic compatibility. Through extensive numerical experiments, we demonstrate that the proposed meshfree scheme achieves first-order accuracy as the horizon parameter and the grid size simultaneously go to 0 with a fixed ratio (i.e., the so-called δ-convergence).

Compared with existing nonlocal convection-diffusion models, we would like to highlight three key contributions of the proposed model. Firstly, we extend the framework of the Gaussian kernel based nonlocal model from anisotropic diffusion to convection-diffusion. Secondly, our model enables a direct meshfree discretization, which preserves the discrete maximum principle and asymptotic compatible convergence and doesn't need additional modification. Thirdly, our nonlocal convection-diffusion model successfully avoids volume correction issues by introducing an unbounded Gaussian-type kernel function, which eliminates the need for volume correction techniques often required in bounded domain models [18,2,7,29].

The paper is organized as follows. The existing Gaussian-type kerne based nonlocal anisotropic diffusion model in non-divergence form is briefly reviewed in Section 2. In Section 3, the proposed nonlocal convection-diffusion model is presented, along with an analysis of its local limit and the selection of the truncated region. The well-posedness and maximum principle of our nonlocal model are established in Section 4. Section 5 introduces the meshfree discretization scheme for numerically simulating the proposed model. Extensive numerical experiments demonstrating the effectiveness of the proposed model and corresponding numerical scheme are presented in Section 6, followed by concluding remarks in Section 7.

We define a d-dimensional multivariate Gaussian distribution, the probability density function p(z, µ, Σ) is presented as follows:

where µ is the expectation and Σ is the covariance matrix. A Gaussian-type kernel based nonlocal diffusion model in non-divergence form was recently proposed in [36] and given by:

where the kernel function γ(x, x ′ ) is defined by

where δ > 0 is a horizon parameter and A(x) is symmetric positive definite and differentiable. As δ approaches 0, the nonlocal operator L nd δ converges to following form:

Note that the Gaussian-type kernel function γ nd in (2.2) is defined over an unbounded area. To ensure computational efficiency, it is viable to truncate the influence region as follows:

An illustration of the Gaussian function p(x, 0, δ 2 A) with A = [1, 0; 0, 1] and δ = 1/10 and the correspondingly truncated influence region χ 2 α (2) = 36 is presented in Figure 1. Specifically, for any given parameter χ 2 α > 0, A(x) shapes the truncated region, while A(x) and δ jointly determine its size. The truncated kernel function γ d,α is then defined as follows:

A volume constrained boundary is associated with the truncated kernel function (2.3) and the domain Ω c is then defined by

Under the volume constraint Dirichlet boundary condition, the corresponding Gaussian-type kernel based nonlocal diffusion problem is given as follows:

where

The above nonlocal diffusion model showcases its effectiveness in accurately simulating a broad range of diffusion processes of both isotropic and anisotropic types. This paper is mainly to extend the capabilities of this nonlocal modeling approach to the convection process and further develop a nonlocal convection-diffusion model with Gaussian-type kernels, in analog to the classic PDE problem (1.1).

In this section, a Gaussian-type kernel based nonlocal convection-diffusion model is presented, and subsequently, the consistency with its corresponding local counterpart (1.1) is proved.

A nonlocal convection-diffusion operator is first defined by:

Here, L nd δ is a nonlocal diffusion operator and takes the specific form of

where

with I d denoting the d-dimensional identity matrix. We remark that the diffusion kernel function

2). Instead, we define A = I d and move ε to the front of the exponential function to deal with the strong scaling effect caused by the small value of ε in the convection-dominated case (see Remark 1 for detailed explanation).

The nonlocal convection operator L nc is defined by:

where

Note that p(x ′ -x, -v(x)δ, δ 2 I d ) is the probability density of a multivariate normal random variable x ′ -x with expectation µ = -v(x)δ and covariance matrix δ 2 I d . The proposed nonlocal convection term differs from the conventional nonlocal convection term in two significant aspects. Firstly, the parameter δ is intricately linked to the covariance matrix and the expectation. When v(x) is fixed, as δ decreases, the kernel function γ nc (x, x ′ ) exhibits a higher degree of singularity, as illustrated in Figure 2. Furthermore, unlike the diffusion term where the center of its density function p(x ′ -x, 0, δ 2 I d ) remains at the point itself, here for each point x, the center of its density function

is fixed, as δ gets smaller, the center of the density function gets closer to the point x. Secondly, the velocity field v(x) is incorporated into the kernel function γ nc (x, x ′ ). As previously mentioned, the center of the density function for the convection term changes due to the influence of v(x), yet its overall shape remains consistently circular. Consequently, the velocity field exerts a direct and discernible impact on the kernel function γ nc (x, x ′ ), as demonstrated in Figure 3.

The associated nonlocal convection-diffusion problem, subject to the Dirichlet-type boundary condition, is then defined as follows:

where (0, 0)

(0, 0) Remark 1. If we follows exactly γ nd defined (2.2) for the nonlocal diffusion kernel, taking A(x) = εI d would give us

The truncated region then becomes

which forms a circle with the radius r d = χ √ εδ. To numerically guarantee obtaining adequate information within this truncated area, a uniform grid with h = r d /χ = √ εδ is deemed necessary for numerical simulation as discussed in [36], which implies the need for highly dense grids for the convection-dominated problem. To resolve this issue, we propose an alternative approach. For a constant coefficient matrix A, as δ → 0, we obtain with (3.8)

By utilizing the eigenvalue decomposition, we have

where

is the diagonal eigenvalue matrix of A with eigenvalues arranged from smallest to largest and Q is the corresponding eigenvector matrix. Then

where

. The nonlocal operator corresponding to the local limit of the reconstructed differential operator

where

For the nonlocal diffusion term (3.2), according to [36], as δ → 0, L nd δ converges to the local diffusion operator, i.e., L nd δ u → L nd u := ε∆u and the approximation error is O(δ 2 ). Let us consider the nonlocal convection term (3.4) and suppose the velocity field v(x) is differentiable. Under the assumption that the solution u is sufficiently smooth, we get

by applying Taylor expansion at x. According to the definition of first-order moment of multivariate normal distribution, the integral of the first term is

Thus we obtain

Consequently, as δ → 0, L nc δ converges to the local convection operator

and the convergence rate is O(δ). Thus, the PDE problem, which is the local counterpart of the nonlocal convection-diffusion problem of (3.6), is given as follows

where Lu(x) = ε∆u(x) -v(x) • ∇u(x).

In practical computations, truncating the influence regions of γ nd (x, x ′ ) and γ nc (x, x ′ ) for two fundamental reasons is imperative. Firstly, the defined area of the kernel function extends to an unbounded space R n . Secondly, as x ′ moves further away from x, the kernel function experiences rapid decay. By selecting an appropriate cut-off distance wisely, the computational domain can be effectively limited to a finite area without significantly compromising model accuracy. This approach enables the practical implementation of the proposed nonlocal model while faithfully capturing the fundamental physical properties of the system being examined.

For the diffusion term, as established in [36], the truncated influence region for x is

and α selected very close to 0. For the convection term, the random variable x ′ -x is assumed to follow a Gaussian distribution of d-dimensional with a covariance matrix of δ 2 I d . Then, ∥x ′ -x + v(x)δ∥ 2 /δ 2 follows a chi-square distribution χ 2 (d). To define the influence region of the kernel function, let us consider all

α , where 0 < α ≪ 1, and χ 2 α is a parameter of the chi-square distribution χ 2 (d), namely the (1 -α) quantile. Hence, we have:

We then define the truncated influence region for a given point x as

where α is chosen very close to 0. In Figure 4, the projection of the iso-density contour of the truncated influence region B nc δ,α (x) onto the coordinate plane is depicted for different velocity fields in two dimensions, where χ 2 α (2) = 36, implying α ≈ 1.52 × 10 -8 . Across various velocity fields, the iso-density contour projections consistently form circular shapes. It is noteworthy that the center of each circle varies, determined by the offset of -v(x)δ, which is established based on the current point x.

Note that the bounded ranges for the diffusion and convection operators are different from each other, B nd δ,α (x) and B nc δ,α (x) respectively. Thus, we correspondingly define the truncated kernel functions γ nd,α and γ nc,α as follows: for any

Figure 4: Sketch of the truncated region B nc δ,α (x) with χ 2 α = 36. Left: the velocity field v(x) = (1, 1) T , the truncated influence region is the circle with (x-δ, y-δ) as the center and 6δ as the radius. Right: the velocity field v(x) = (3, -3) T , the truncated influence region is the circle with (x -3δ, y + 3δ) as the center and 6δ as the radius.

The corresponding interaction domains about diffusion term and convection term are defined as

For the proposed nonlocal convection-diffusion problem, the volumetric constraints are imposed on the interaction domain Ω bd ∪ Ω bc . In contrast, the nonlocal operator equation is applied in the domain Ω s . These volume constraints are natural extensions of the boundary conditions for problems involving differential equations in the nonlocal context. The interaction domains Ω bd , Ω bc , and the total interaction domains Ω c = Ω bd ∪ Ω bc are visually depicted in Figure 5. Finally, under the volumetric constraint Dirichlet boundary condition, the nonlocal convection-diffusion model with truncated influence regions is given by:

where

with Ω = Ω s ∪ Ω c .

•) is denoted as the L 2 inner product. For x ∈ Ω s , with g(x) = 0 and f ∈ L 2 (Ω s ), the weak form of our models is as follows:

where

and

For γ nd,α and γ nc,α , we first note the existence of a positive constant K 1 (δ), dependent on δ, such that the kernel γ nd,α satisfies

Moreover, there exist positive constants K 2 (δ) and K 3 (δ), dependent on δ, ensuring that the kernel

Based on these conditions, the nonlocal operator L δ,α is bounded in L 2 n0 (Ω), which guarantees that a weak solution is also a strong solution of (3.6) in L 2 n0 (Ω). Consequently, for the nonlocal problem (3.6), the well-posedness with the following result is established.

Theorem 1. (Well-posedness) Suppose δ > 0 is fixed and the kernel γ c,α (x, x ′ ) satisfies

Consequently, a unique solution u ∈ L 2 n0 (Ω) exists for the nonlocal convection-diffusion problem (3.6). Moreover, this solution adheres to the prior estimate:

Proof. First, we can simply obtain

With a proof similar to that of Theorem 2 in [34], we demonstrate the boundedness of the bilinear operator

(Ω), and due to the symmetry of γ nd,α (x, x ′ ) = γ nd,α (x ′ , x), the bilinear operator B nd (u, v) can be expressed as:

Based on (4.4), we obtain For any u, v ∈ L 2 n0 (Ω),

Note that B nc (u, v) can be rewritten as

(4.12)

By utilizing the Cauchy-Schwartz inequality, as well as the inequalities in (4.5) and (4.6), the first term to the right of (4.12) can be qualified as follows:

Similarly, the second term to the right of (4.12) satisfies

Then we obtain

Thus the combination (4.11) and (4.15) gives us

Subsequently we demonstrate that the bilinear operator B(•, •) is coercive on L 2 n0 (Ω). It holds that

Using the inequality in (4.3), we then can get

For the bilinear operator B nc (•, •), we first note that

It is evident that

is always symmetric, even though γ nc,α (x, x ′ ) may not be. Therefore, the first term of the right of (4.18) can be expressed as:

With the assumption (4.7) and the equality

the second term of the right of (4.18) can be rewritten as

With the combination of (4. 19) and (4.20), we obtain

Therefore, using the Lax-Milgram theorem, the problem (4.1) exists as a unique solution u ∈ L 2 n0 (Ω). Furthermore, since

the a priori estimate (4.8) can be obtained.

The condition (4.7) can be considered a nonlocal analog of the following condition

Here is a simple proof: L 2 n 0 (Ω) = {u ∈ L 2 (Ω) | u(x) = 0 on ∂Ω} is defined as a constrained space. For any w ∈ L 2 n 0 (Ω), we have

where n(x) is the normal vector pointing to the outside of the domain. With the Taylor expansion at x, we also have

Hence, as δ → 0, the condition (4.7) can be considered a nonlocal analog of the condition (4.23).

Theorem 2. (Maximum principle) Suppose -L δ,α u < 0 in Ω s , then in the interaction domain Ω c , we can attain a non-negative maximum of u.

Proof. Through proof by contradiction, we aim to demonstrate that a nonnegative maximum cannot be achieved within Ω s . For the sake of contradiction, let's assume that a nonnegative maximum u at x 0 ∈ Ω s can be achieved. Consequently,

Cause u(x ′ ) -u(x 0 ) ≤ 0, we can simply demonstrate that

and

which give us a contradiction with the assumption of -L δ,α u(x 0 ) < 0.

In order to numerically simulate the proposed nonlocal model (3.22), we discretize it by following the meshfree approach proposed in [31]. Suppose that the domain Ω is discretized into nodes {x i }, and each node x i in the reference configuration has a known associated volume S i . The nodes within Ω s are denoted as {x 1 , • • • , x Ns }, while the nodes along the nonlocal boundary Ω c are denoted as {x Ns+1 , • • • , x Ns+Nc }. Collectively, the nodes and the associated volumes constitute a grid set denoted as T h . The method is considered meshfree, signifying the absence of elements or other geometrical connections between the nodes.

For the nonlocal diffusion term in (3.23), at the node x i ∈ Ω s , using the meshfree discretization associated with T h , we have

where S j = B nd δ,α (x i )∩S j . For the nonlocal convection term in (3.23), at the node x i ∈ Ω s , similarly, we approximate it by

where Ŝj = B nc δ,α (x i ) ∩ S j . For each point x i , the geometries of S j and Ŝj may exhibit irregularities. The irregular intersection of these regions holds the potential to introduce quadrature errors, thereby influencing the overall accuracy of the simulation. However, leveraging the unbounded nature of the kernel function in the context of the diffusion term allows us to slightly extend the integral regions from the irregular intersection S j to the regular volume S j as illustrated later in the section of numerical experiments. This strategic extension addresses the issue of irregular intersections while circumventing the need for volume correction. A similar approach is applied to Ŝj for the convection term. This expanded integration technique facilitates a straightforward and precise alignment of the integration area with T h .

In conclusion, the meshfree discretization scheme to solve the nonlocal convection-diffusion problem (3.22) is expressed by: find (u h (x 1 ), u h (x 2 ) . . . , u h (x Ns )) such that

Let us define

)

For the meshfree discretization (5.3), the resulting linear system then can be obtained as

where

Remark 2. In contrast to the proposed Gaussian-type kernel based nonlocal model, many existing nonlocal models restrict nonlocal interactions to bounded neighborhoods, often chosen as Euclidean balls B δ (x). The approximate balls B δ,h (x), typically composed of polygons, impose a challenge when intersecting for meshfree discretization method. An important question arises: to what extent do such approximations impact the nonlocal operators and the corresponding solutions? Recent works have delved into this issue [7,18]. A notable convergence result, presented in Corollary 4.2 of [7], is as follows:

where C e represents a norm-equivalence constant, K represents a positive constant depending on the data f and g but independent of δ and h, and △B δ,h (x) denotes the ball difference.

Theorem 3. The stiffness matrix A h given by (5.4) is an M-matrix which is nonsingular. Thus, the linear system (5.5) is uniquely solvable. Moreover, the discrete maximum principle is satisfied by u h when the boundary values g = 0: if the source terms f ≤ 0 in Ω s , then max 1≤i≤Ns u h (x i ) ≤ 0, and if the source terms f ≥ 0 in Ω s , then

Hence, using (5.4), for any j ̸ = i, we can deduce that a i,j < 0 if

As a result, we establish that the stiffness matrix A h qualifies as an M-matrix. This property ensures the existence of A -1 h and guarantees its absence of negative entries. Consequently, we can infer that the linear system (5.5) has the unique solution, and the discrete maximum principle is satisfied.

Regarding the asymptotic compatibility, while the theoretical proof remains pending, we will provide numerical evidence in Section 6.2 to illustrate that the meshfree discretization(5.3) achieves δ-convergence. Notably, we will demonstrate that the numerical solutions of the nonlocal model under a fixed ratio between δ and h exhibit the first-order convergence towards the corresponding local PDE solution (3.15).

In this section, a series of numerical experiments are carried out in two dimensions. These experiments aim to demonstrate the application of our nonlocal model (3.22) and the effectiveness of the meshfree discretization scheme (5.3). Additionally, the discrete maximum principle is assessed. On the approximation accuracy of L δ,α to L δ , we examine the impact of how to choose the truncation parameter χ 2 α . It is noteworthy that we consistently set χ 2 α = 36 as the default value to truncate the influence region in Examples 1-4, according to the comparative results observed in Example 5.

To test the convergence of the meshfree discretization to solve our nonlocal model (3.22), we initially maintain a constant value for δ.

Example 1. We consider the two-dimensional domain Ω s = (0, 1) × (0, 1), with the diffusion coefficient ε = 1, the vector field v(x, y) = (1, 1) T and the horizon parameter δ = 1/80. Let us choose u(x, y) = sin(x 2 + y 2 ), u(x, y) = e xy , and u(x, y) = xy 5 as the exact solution for distinct scenarios. The boundary values g(x, y) are directly determined rom the exact solutions u(x, y). The source terms f (x, y) are determined by the nonlocal model (3.6).

We utilize a uniform domain partition of Ω s with N × N grid points, where N assumes values of 20, 25, 30, 35, 40, and 50, respectively. Table 1 reports the L 2 errors and convergence rates resulting from the meshfree discretization (5.3). As anticipated, exponential convergence is observed across all cases until the model errors due to the truncation of the influence regions dominate.

We now shift our focus to analyzing the convergence behavior of the proposed meshfree discretization (5.3) as the horizon parameter δ tends towards zero. Specifically, we explore a scenario in which both the grid size h and the horizon parameter δ decrease to zero while maintaining a fixed ratio between them. This so-called δ-convergence test is a widely adopted method to verify the asymptotic compatibility of numerical schemes. We anticipate observing the convergence of the approximate solutions derived from the nonlocal convection-diffusion problem (3.22) towards the solution that corresponds to classical local problem (3.15). This behavior aligns with the recognized continuum limit; that is, greater smoothness and continuity are achieved gradually by the underlying physical system. Through a meticulous examination of the δ-convergence of the proposed meshfree scheme, our objective is to assess its accuracy and reliability. This endeavor also provides valuable insights into the fundamental physics governing the system under scrutiny.

Example 2. Let us consider the two-dimensional domain Ω s = (0, 1) × (0, 1) with the vector field v(x, y) = (1, 1) T , and set u(x, y) = sin(x 2 + y 2 ) as the exact solution for the classical PDE problem (3.15). Two choices of the diffusion coefficients ε = 1 and ε = 10 -7 (strongly convection-dominated) are used respectively. The boundary values g(x, y) are directly derived from the exact solution u(x, y) and the source terms f (x, y) are determined from the nonlocal model (3.6).

According to the analysis of truncation about the influence horizon, the influence region for the diffusion term is given by

and that for the convection term by

At any given point (x, y), the influence region of the diffusion term is represented by a circle centered at (x, y) with a radius of 6δ. However, the influence region of the convection term is a circle centered at (x -δ, y -δ) with a radius of 6δ. The configuration of Ω bd takes the form of an equal-width band encircling the boundary of Ω s , while Ω bc exhibits an unequal-width band surrounding the same boundary. Given that the kernel function is defined over an unbounded region and exponentially decays, volume correction is unnecessary for the irregular tangent part of the boundary region and the mesh. This allows for the appropriate expansion of the boundary scope, facilitating more convenient calculations. It is important to note that enlarging the integration area will not compromise the accuracy of the numerical computations. We divide Ω s into N × N grids and record h = 1/N and set δ = h, δ = 2h, and δ = 4h. Figure 6 provides a visual illustration of the domain Ω = Ω s ∪ Ω c and the influences regions B nd δ,α and B nc δ,α . Table 2 reports the L 2 errors and convergence rates resulting from the meshfree discretization (5.3), and we clearly observe the first order convergence along the grid refinement.

Example 3. The configuration for this example closely mirrors that of Example 2, with one notable distinction: a variable velocity field is introduced, defined as v(x, y) = (sin 2 (πx + πy), cos 2 (πx + πy)) T . This velocity field is visually depicted in Figure 7(a). Similar to Example 2, the influence region for the diffusion term is given by

and that for the convection term by At any given point (x, y), the influence region for the diffusion term is represented by a circle centered at (x, y) with a radius of 6δ. However, the influence region for the convection term is a circle centered at (x-sin 2 (πx+πy)δ, y -cos 2 (πx+πy)δ) with a radius of 6δ. It's worth noting that due to the variability of the velocity field, this leads to the inconsistent displacement of the circle's center point where the convection term is truncated at the boundary. Consequently, this results in an irregular boundary, as illustrated in Figure 7(b). Given that the velocity field is bounded, we extend the boundary by a length of 7δ in both the x and y directions to ensure computational efficiency. This extension yields a regularized boundary denoted as Ω c , as depicted in Figure 7(c). Table 3 reports the L 2 errors and convergence rates resulting from the meshfree discretization (5.3), and we again observe the first-order convergence along the grid refinement.

v(x, y) = (sin 2 (πx + πy), cos 2 (πx + πy)) T N ε = 1 ε = 10 -7 L 2 error CR L 2 error CR 10 8.2979 × 10 -3 -4.9168 × 10 -2 -20 3.3121 × 10 -3 1.33 2.8121 × 10 -2 0.81 40 1.4632 × 10 -3 1.18 1.5941 × 10 -2 0.82 80 6.8512 × 10 -4 1.09 9.6332 × 10 -3 0.73 160 3.3141 × 10 -4 1.05 4.7866 × 10 -3 1.01 320 1.6267 × 10 -4 1.03 2.5196 × 10 -3 0.93

Table 3: Numerical results on L 2 errors and convergence rates for different exact solutions with δ = h in Example 3.

In this subsection, we focus on testing the discrete maximum principle for the meshfree discretization (5.3). We utilize a uniform partition of the domain Ω s with 40×40 grid points and employ the meshfree discretization (5.3) to solve the nonlocal model. Figure 8 showcases the resulting numerical solutions obtained with the two distinct velocity fields, and it is evident that the discrete maximum principle is well maintained in all cases. 6.4. Effect of χ 2 α on the approximation accuracy of L δ,α to L δ In this subsection, we choose different values for χ 2 α to test effectiveness of the truncated nonlocal operator L δ,α defined in (3.22) as an approximation of the nonlocal operator L δ defined in (3.1). Example 5. The experimental settings are similar to Example 2, with the only alteration being the exploration of various values for χ 2 α . Specifically, we set ε = 1 and consider χ 2 α = 9, 16, 25, 36, and 49, respectively. Table 4: Numerical results on L 2 errors with different choice of χ 2 α in Example 5.

exhibit a rapid decrease as χ 2 α increases from 9 to 49. Notably, the disparities in solution errors between χ 2 α = 36 and χ 2 α = 49 are almost negligible. Consequently, we suggest the adoption of χ 2 α = 36 in practical applications. This choice strikes a balance between ensuring the accuracy of the nonlocal convection-diffusion model (3.22) and keeping computational efficiency.

In this paper we introduce a novel nonlocal convection-diffusion model based on the Gaussiantype kernel, expanding upon an existing nonlocal diffusion model. The key innovation lies in the integration of the velocity field into the expectation, utilizing a truncated multivariate Gaussian function as the kernel. The well-posedness and elucidation of certain inherent properties are established to assess the robustness of our proposed model. For the numerical solution, we design a direct meshfree discretization method that adheres to the discrete maximum principle. A series of numerical experiments are also carried out in two dimensions to illustrate the versatility of our model in tackling diverse convection-diffusion problems and the robustness of the proposed numerical scheme.