Predictions of the future state of a flow, here internal 2D or 3D flow of an incompressible viscous fluid in a domain , + ç ∇ -+ ∇ = ( ), and ∇ ç = 0, in , 0 < d , (
= 0 on and ( , 0) = 0 ( ),
are improved, Kalnay [1], by incorporating / assimilating observations of the flow ( , ) for 0 < d .
The goal of data assimilation is to combine incomplete, sparse, noisy, and possibly time-delayed observations, , with incomplete and approximate knowledge of a flow's dynamic laws to produce a more accurate prediction of flow statistics. Nudging-based assimilation, from 1964 Luenberger [2], is amenable to both analysis, e.g. Biswas and Price [3], Azouani, Olson and Titi [4], Cao, Giorgini, Jolly and Pakzad [5], Larios, Rebholz, and Zerfas [6] (among hundreds of papers), and straightforward implementation. However, nudging is non-optimal since no criterion is minimized, Lakshmivarahan and Lewis [7]. Uniform in-time convergence has even been proven [6] under parameter conditions reviewed in Section 2. A difficulty remaining is that users must select effective nudging parameters and the analytical theory places parameter restrictions, see Section 2.1, far beyond ranges found effective in computational experience, summarized in Kalnay [1]. Since á priori analysis develops conditions from a series of worst-case estimates, these restrictions may be pessimistic, and self-adaptive parameter selection, developed herein, is natural.
To develop adaptive-methods, we adopt the simplest interesting setting. Let the true velocity be denoted ( , ) and its observed values = be an 2 projection of the true values on a finite-dimensional subspace (associated with a length-scale denoted ). The continuum, nudged approximation ( , ) satisfies select parameter , set ( , 0) = 0 ( ) and solve + ç ∇ --( -) + ∇ = ( ), ∇ ç = 0, in .
(1.
3)
The initial and boundary conditions are = 0 on , = 0 on , ( , 0) = 0 ( ), and ( , 0) = 0 ( ).
To support the necessity for adapting , in Section 2 we summarized error analysis (inspired by the analysis in [6]) showing that for small enough and large enough nudging is accurate uniformly in time. Specifically, the conditions from this analysis are
2D -condition: -2 1
) 3D -condition: 2 [ -2048 19683 -3 1 + 0 ‖∇ ‖ 4 ] e 0 > 0. (1.6) With explicit time discretization, the nudging term alters the CFL stability condition and a third parameter constraint emerges, detailed in Farhat, Larios, Martinez, and Whitehead [8]. For high Reynolds number flows in 2D and 3D, in Section 2.1 we investigate what turbulence phenomenology suggests large enough and small enough in (1.4)-(1.6) mean. There results in
in 3D: ≳ ý(Ā 5 ) and ≲ ý(Ā -3 ).
(
While sufficient (in theory), these are impractical (in practice). Extensive practical experience, Kalnay [1], shows that the nudging parameter is best chosen for intermediate size, balancing data errors with model errors and numerical errors. Thus the parameter requirements (1.4)-(1.6) are far beyond parameter values in problems where assimilation is used. The gap is likely because of conditions (1.4)-(1.6) emerge from analysis through a series of worst-case estimates which may not occur simultaneously and (1.7), (1.8) emerge from the assumption that the flow is fully developed turbulence. To address this issue, in Section 3 we present 2 low complexity algorithms for self-adaptive parameter selection. These act to find a smaller effective value in response to the local flow behavior and are only slightly more expensive in operations and memory than á priori selection.
Adaptivity exploits -sensitivity to improve approximations. It is based on 3 principles: localization (of the global criteria), estimation (of the deviation from a desired state), and decision (of how to choose the next parameter value). The condition that ‖ (ç, ) -(ç, )‖ is decreasing is localized to decreasing every time step. The estimator is ‖ ( -)‖ and the decision tree is to increase for large estimates and decrease for very small estimators. We compare two decision strategies based on derived analytical criteria for decrease and simply asking if a reduction is observed.
Section 4 presents tests of the adaptive nudging methods for problems with model errors with both á priori and self-adaptive parameter selections. We find both adaptive-methods produce effective values smaller than (1.7)-(1.8) and accurate velocities. We also observe that an H-condition seems to be necessary for accurate, long-time approximations. Addressing the H-condition is an open problem. In many applications, depends on the density of sensors.
Nudging for data assimilation stems from the 1964 work of Luenberger [2] and has extensive tests in geophysical flow simulations, e.g., [1,7,[9][10][11][12]. Its form also arises under the terms various types of damping and time relaxation, e.g., [13], with large numbers of papers for each realization. The analytical work on nudging is similarly enormous; we mention three that we (and most others) build upon. Azouani, Olson, and Titi [4] gave an early and possibly even the first complete convergence analysis of a nudging algorithm. Their analysis produced a similar -condition as herein in Theorem 2.1 p. 254. Biswas and Price [3] gave a remarkable analysis in 3D without assuming á priori strong solutions. Their analysis recovered similar conditions as herein on (their equation (3.27)) and (their (3.8)). Biswas and Price [3] also provided an early algorithm for adaptation of which retained a Reynolds number dependent lower bound on in their equation (4.4). In 2019 Larios, Rebholz, and Zerfas [6] proved uniform in time error bounds in a way that narrowed the gap between analysis and computational tests. In Zerfas, Rebholz, Schneier, and Iliescu [14], the nudging parameter (here ) is adapted when the nudged system's kinetic energy is out of the bounds that the NSE energy inequality imposes. See also García-Archilla and Novo [15] for discrete-time error analysis.
None of these papers sought an optimal . Nudging emerged from Luenberger's work on control theory was then developed in independent directions. Along its path, several reconnected nudging to control and optimization, determining methods for optimization of the nudging parameter. For work on this approach (not considered herein) involving solving the adjoint problem, see Hoke and Anthes [9], Stauffer and Bao [11], Zou, Navon, LeDimet [12] and Navon [10]. Carlson, Hudson, and Larios [16] and Martinez [17] have studied recovering the unknown viscosity parameter for the 2D Navier-Stokes equations (NSE). Parameter estimates for the Lorenz equation are studied in Carlson, Hudson, Larios, Martinez, Ng, and Whitehead [18]. Finally, after this report's submission, two interesting reports exploring other approaches to adapting have become available on ArXiv, Carlson, Farhat, Martinez, and Victor [19,20].
The notation for the 2 ( ) norm and inner product is ‖ ç ‖ and (ç, ç), respectively. Using ‖ ç ‖ , we indicate the ( ) norm. The numerical tests use a standard finite element discretization, see [21] for details. We also use the inequalities derived by Ladyzhenskaya [22].
Theorem 1.1 (The Ladyzhenskaya Inequalities, See Ladyzhenskaya [22]). For any vector function ∶ R ³ R with compact support and with the indicated norms finite,
, ( = 2),
‖∇ ‖, ( = 3).
To motivate the necessity of adaptivity, for completeness, we review the (now standard since [3,4,6]) proof that the error ³ 0 uniformly in as ³ @ and uniformly in as ³ @. The (direct) proof shows how conditions small enough and large enough emerge organically for the structure of the NSE and it motivates one adaptive algorithm of Section 3. This conditions small enough and large enough are interpreted under 2D and 3D turbulence phenomenology in Section 2.1 below, leading to (1.7). We assume here that is an 2 projection and satisfies
For the 3D case, we assume 1 that ‖∇ ‖ 4 * 1 (0, ). This assumption is not necessary (since the work of Biswas and Price [3]) but shortens the analysis.
Proposition 2.1. Let = -. Suppose is an 2 projection with the approximation property (2.1) above. In 2D let the parameter conditions (1.4), (1.5) hold and in 3D (1.4), (1.6). Then, the error ³ 0 uniformly in as ³ @ and uniformly in as ³ @. In particular,
Proof. The 2D case: By subtracting (1.3) from (1.1) and taking the inner product with , there follows
By the 2D Ladyzhenskaya and arithmetic-geometric mean inequalities
Thus,
Provided the term in the first bracket is non-negative, the condition (1.4), we proceed as follows. Denote
. Thus, provided (1.5), i.e., 1 One can assume a stronger condition, such as ‖∇ ‖ * @ (0, ), and obtain an estimate that superficially appears better behaved concerning the Reynolds number or a weaker condition, such as Prodi-Serrin, and obtain an estimate that has a worse explicit dependence. This is because some dependence is subsumed within the assumption.
Computer Methods in Applied Mechanics and Engineering 433 (2025) 117526
] and the error ³ 0 uniformly in as ³ @ and uniformly in as ³ @.
The 3D case: Starting at
By the 3D Ladyzhenskaya inequality and an arithmetic-geometric mean inequality(with exponents 4 and 4/3) 2
Thus,
Following the 2D case, the error satisfies
] e 0 > 0. ¦
In the analysis, conditions on and arise naturally for uniform in convergence. The question we now investigate for higher Reynolds number, turbulent flows is How severe are these conditions for practical problems? Let the large scale length and velocity be denoted , , e.g.,
has units of length and of -1 . In 2D, the conditions are -2 2 1 2 e 0 and -2 1
For fully developed turbulent flows in 2D with periodic boundary conditions, Alexakis and Doering [23] have proven the following, consistent with phenomenology,
where Ā = and = average energy-input wave number,
In 2D, phenomenology interprets the -condition to mean
The -condition, -2 2
1 2 e 0, can be rearranged to read
Since ≳ ý(Ā 3∕2 ) the -condition becomes ( ∕ ) 2 d ý(Ā -5∕2 ). These conditions on and are severe. In 3D, the scaling becomes worse. For fully developed turbulent flows in 3D with periodic boundary conditions, the following has been proved by Doering and Foias [24], again consistent with turbulent phenomenology in the large limit
Under the (often called a spherical cow assumption) condition that 2 The equality ( ç ∇ç ∇ , ) = ( ç ∇ , ) requires skew symmetry of the trilinear form ( ç ∇ , ). After space discretization, an analogous equality holds with explicitly skew symmetrized trilinear forms or when divergence-free elements are used.
Computer Methods in Applied Mechanics and Engineering 433 (2025) 117526
we then estimate
The condition (1.6), big enough, has the interpretation
The small enough condition then becomes
1
The 3D condition * ≳ ý(Ā 5 ) and ∕ ≲ ý(Ā -3 ) are severe restrictions. For problems with smooth solutions (e.g., test problems constructed by the method of manufactured solutions) the -condition can be altered as follows. For a vector function ( ) define ( ) as
Here represents an average length-scale of the solution, analogous to the Taylor micro-scale (which would also involve time averaging). Then ‖∇ ‖ d -1 ‖ ‖. The term that gives rise to the -condition is now estimated by
It can now be subsumed into the -condition (in 3D here). This gives long time estimates when is small with respect to ( ) (rather than Ā ) from the inequality.
The severe constraints, in 2D ≳ ý(Ā 3∕2 ) and ≳ ý(Ā 5 ) in 3D can possibly be improved since analysis gives a sufficient (but possibly not sharp), worst-case lower bound. It is therefore interesting to develop methods for self-adaptive parameter selection that respond to the local state of the flow. The condition on is linked to the condition on ; improving the latter improves the former. Further, changing requires changing observation locations rather than just picking a new user-defined parameter. We thus focus on adapting . Reformulating nudging to improve the H-condition is an open problem.
The first algorithm is not grounded in the á priori analysis other than the result that a decrease of ‖ ‖ is possible and ‖ ( )‖ d ‖ ‖. It uses the fact that ‖ ‖ is computable in a reasonable but heuristic way as follows.
Choose 0 > 0. Set upper safety factor, Factor e 1, and lower tolerance, The second algorithm is grounded in the á priori analysis that led to conditions (1.5), (1.6). We first develop it in 2D. Changing the parameter values at each time step means the nudging term is no longer autonomous. Thus we set = ( ). Since we are not addressing the -condition the analysis is shortened by setting = . Alternatively, the analysis below could proceed (with more terms) with the -condition present and assumed satisfied. For the tests in Section 4, . We begin with
For the nonlinear term we rearrange and decompose on ( ç ∇ç ∇ , ) = ( ç ∇ , ) rather than ( ç ∇ , ) since is computed. Following the subsequent steps (making this one change) we obtain
Denote ( ) ∶= 2 ( ) --1 ‖∇ ‖ 2 , 0 d ( ) = ‖ ‖ 2 then 2 + ( ) d 0. Thus, the error at time satisfies and the contribution to the evolution of the error of one-time step of size satisfies
This suggests the following algorithm with piecewise constant .
Choose 0 > 0. Given such that
‖∇ ( )‖ 2 e 0 > 0.
Proceed to the next step after resetting
[Too small value] If
Repeat step after resetting
In 3D the only change is to replace
In both cases, the integrals over < < +1 are evaluated by quadrature using a rule with accuracy compatible with the time discretization employed. In our tests the trapezoidal approximation is used with the approximate velocities at , +1 . The 3D multiplier -3 reflects the greater complexity of 3D turbulence, Section 2.1. Users can define a maximum number of retries, denoted as , for both algorithms. In our tests in Section 4, we set = 5.
In this section, we conduct three tests to validate the adaptive algorithms. In the first test, we calculate the rate of convergence with a test with exact solutions. Then we test the adaptive algorithms with a complex flow between offset cylinders from [25] with a higher Reynolds number. In the third test, we use the test case from [26], which involves flows past a flat plate with a of 50.
We note that the two adaptive strategies have similar but different aims. The first aims to enforce ‖ ( ( + ))‖ d ‖ ( ( ))‖. The second aims to enforce ‖ ( + )‖ d ‖ ( )‖. Since ‖ ( )‖ d ‖ ‖ we expect the first to yield smaller -values than the second.
In the numerical example, we verify the theoretical temporal convergence results with the BDF2 (second-order backward differentiation formula) time discretization described and analyzed in Section 3 of Larios, Rebholz, and Zerfas, [6]. The following exact solution for the NSE is considered in the domain = (0, 1) × (0, 1). The velocity and the pressure are as follows:
( , , ) = (cos , sin ), and ( , , ) = ( -)(1 + ).
We insert these in the NSE to calculate the body force ( , ). We choose = 1. We create a fine mesh resolution that consists of 55554 degrees of freedom (dof) to isolate the convergence rates for of the errors by making the spacial errors small. We run the code in the time interval [0, 2]. We use a Scott-Vogelius finite element pair with a barycenter refined mesh. We observed similar results with Taylor-Hood elements and a skew-symmetrized nonlinearity. The measurements were made by 2 ( ) norm. In Algorithm 1, we set a safety factor and a lower Tolerance as follows. The factor is user-defined to avoid growing too big when the projection error grows slightly.
Doubling when ‖ ( ( + ))‖ > Factor ‖ ( ( ))‖.
We set Factor = 1.3.
where Tol is also defined by the user and Tol = 0.2 in this test.
In Table 1, we observe second-order convergence in for velocity. This is an optimal rate and is evidence that the implementation is correct. We also present the values of the maximum value of through adaptation to understand how large the generated values become. The small values seem to be in response to the solution's smaller . According to the results of this accuracy test, theoretical expectations are compatible with numerical results, except that, oddly, the -values of Algorithm 1 are larger than Algorithm 2. This effect reverses over a longer time. Due to this reversal of the expected values, we monitor the behavior of the relative errors and values for a longer time, final time = 10, see Figs. 1 and 2. The relative errors are ý(10 -6 ) for a long time. Then they start growing, but the growth speed is slowing down and the errors appear to be saturating around 10 -6 similar to logistic growth. Logistic growth is a long-accepted model of error growth in fluid dynamics, e.g., Lorenz [27]. For Algorithm 1, the adapted value is less than 50, and for Algorithm 2, is growing until it reaches the maximum value 10 6 . Algorithm 1 provided smaller -values is consistent with expectations, based on ‖ ( )‖ d ‖ ‖.
In this example, we are testing the adaptive nudging algorithms with a higher flow with complex features and without exact solutions. The domain is a disk with a smaller off-center obstacle inside. Let the outer circle radius 1 = 1 and the inner circle radius The computational domain with its triangulation is presented in Fig. 3. A counterclockwise rotation drives the flow with the body force given below:
with no-slip boundary conditions. The outer circle remains stationary. The Delaunay algorithm generates the mesh with 75 mesh points on the outer circle and 60 mesh points on the inner circle. The final time = 10, time step size = 0.01, = 10 -3 , = 1, = 1, and = . Initial condition ( , , 0) = 0 and the Dirichlet boundary condition = 0 on . The relative error is defined as ‖ -‖ ‖ ‖ starting at * = 1. In Fig. 6, we can see the adaptive algorithms and the constant nudging effectively decrease the relative error from = 1 to = 2. After = 2, the relative error grows fast until it reaches ý(1) as in the previous test, again suggesting -condition violation. In Fig. 7, the values of Algorithms 1 and 2 grow with time until they reach the maximum value. We observe that the of Algorithm 1 grows fast at first, then the growth rate gets smaller, and for Algorithm 2 at the beginning barely changes and sharply grows after = 2. These two tests suggest that both a -condition and an H-condition are essential for long-time convergence. This test also suggests that data from a turbulence model does not alleviate non-satisfaction of nudging's -condition.
We adopt this test from [26]. The computational domain is a [-7, 20] × [0, 20] rectangular channel with 0.125 × 1 flat plate obstacle so that | | H 540. The statistics are relative to solution size which also makes them area independent. The inflow velocity is set with = ï1, 0ð and no force is applied. No-slip boundary conditions are applied on walls and the plate whereas a weak zero-traction boundary condition is enforced on the outflow boundary as in [26]. We choose Ā = 50, = 1 Ā , = 0.02 and the final time = 81, three times turnover time. The flow is at rest at = 0. We approximate the exact velocity via finer mesh as seen in Fig. 8 and solve the equation with coarser mesh. The total dof for the fine mesh is 27373 and for coarse mesh is 15037. Since this test is a through-flow problem, it interrogates errors over 1 through-flow time, not over longer times. Still, it is a complex flow with many interesting features and an accepted test problem.
We compare Algorithms 1, 2 and the constant with different initial values, where = 1, 10, 100, and 1000. We set Factor = 1 and = 0.2, and the maximum of = 10 6 . We define the relative error ‖ -‖ ‖ ‖ , and relative projection error
. In Figs. 10-16, both the relative error and the relative projection error decrease with larger values for Algorithm 1 and the constant . Algorithm 1 performs better when the initial value is small because it effectively adapts to reduce the projection error. On the other hand, Algorithm 2 is robust to initial values, adjusting to obtain the optimal value for smaller errors at each time step. We observe that the projection error is smaller than the error. Furthermore, the value for Algorithm 1 is smaller than Algorithm 2 in the first turnover time, and it is bigger than Algorithm 2 in long times (see Fig. 9). We obtain smaller relative error and projection errors for Algorithm 2 than for Algorithm 1 and the constant . We would like to point out here that if and are solved on the same mesh (as in, e.g., [26]) instead of fine mesh coarse mesh convention as has been done here, the error is negligible, ý(machine epsilon).
Two adaptive-strategies were derived and tested: enforcing a decrease in the observed error, ‖ ( -)‖, and enforcing a computable sufficient condition for the full error, ‖ -‖, decrease. The second, based on a stronger condition, generally yielded larger values. The tests at high Ā have the adaptive algorithm steadily increasing . This means the large error analysis in Diegel, Li, and Rebholz [32] is of special interest. In one test (see Fig. 6) the second responded more quickly to transient behavior and resulted in significantly smaller errors. Adaptive selection of the nudging parameter improves performance. It does not, however, eliminate the necessity of the -condition related to the density of observations. This means revising the standard formulation of nudging to improve the -condition is an important open problem. Other critical open problems include the effect of time delays in the nudging control term and the effectiveness of nudging for correction of model errors. The precise realization of Algorithm 2 and its supporting discrete-time analysis is also an open problem for higher-order linear multi-step methods.
The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: William Layton reports financial support was provided by National Science Foundation. Rui Fang reports financial support was provided by National Science Foundation. Farjana Siddiqua reports financial support was provided by National Science Foundation. Aytekin Cibik reports financial support was provided by TUBITAK. If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Computer Methods in Applied Mechanics and Engineering 433 (2025) 117526