Initial waveform optimization based on the adjoint variable and the finite element methods
Takahiko Kurahashi*
Department of Mechanical Engineering, Nagaoka University of Technology
Nagaoka University of Technology, 1603-1 Kamitomioka, Nagaoka, Niigata 940-2188, JAPAN Email: [email protected]
Shigetoshi Tai
Department of Mechanical Engineering, Nagaoka University of Technology
Nagaoka University of Technology, 1603-1 Kamitomioka, Nagaoka, Niigata 940-2188, JAPAN Email: [email protected]
Summary
In this study, we present the initial waveform optimization analysis based on the adjoint variable and the finite element methods. As the governing equation, a shallow-water equation is employed to simulate flow behavior, and the finite element method using the stabilized bubble function element and the Crank-Nicolson method were used for spatial and temporal discretization. The results of the investigation for the adjoint boundary condition are shown in this paper.
Keywords: Initial waveform optimization, Adjoint variable method, Finite element method.
1 Introduction
In this study, the initial waveform optimization analysis is carried out based on the adjoint variable and the finite element method. To obtain the optimized initial waveform using the time history of observed water elevation, the performance function, i.e., the cost functional, should be defined[1],[2]. The performance function is defined as the sum of the squares of the residual between the computed and observed water elevations. The performance function, is expressed as
J = 1 2
∫ t
ft
0{ η − η obs. } T [Q] { η − η obs. } dt, (1) where [Q] is the weighting diagonal matrix, η is the computed water elevation and η obs. is the observed water elevation.
The purpose of this study is to find the distribution of the initial water elevation η (t 0 ) to minimize the performance function under the constraints of the shallow-water equation, the initial and the boundary conditions. Minimizing the performance function means that the computed water elevation should be as close as possible to the observed water elevation at the observation points. The results of the investigation of the adjoint boundary condition in the inverse analysis based on the adjoint variable method are shown in this study. In addition, the performance function extended by the governing equation, i.e., a shallow-water equation, and the adjoint
variable is expressed by J ∗ in this paper. J ∗ is referred to as the Lagrange function and the extended performance function.
2 Numerical experiment
The purpose of this study is to determine the distribution of water elevation to bring the computed water elevation close to the observed water elevation at the observation points based on the adjoint variable method. The finite element mesh and the location of observation points are shown in Figure 1. The total number of nodes and elements are respectively 1681 and 3200. The computational conditions for this analysis are listed in Table 1. The weighting constant Q at the observation point is 1.0 and at the other points are 0.0. In this study, the artificial observed water elevation obtained by using the exact initial water elevation shown in Figure 2 is employed. The mean water depth h is set at 10.0 m in the computational domain. The distribution of the initial water elevation shown in Figure 2, is expressed in the form of the following equations.
η (x,y) = 0.5 × ( cos
( π r 30 )
+1.0 )
r ≤ 30.0, (2) η (x, y) = 0.0 r > 30.0. (3) Here, the radius r is calculated as
r =
√
(x − 50) 2 + (y − 50) 2 . (4)
Table 1: Computational conditions
Number of time steps 4000
Time increment ∆t (sec) 0.005
Manning coefficient of roughness n (sec/m 1/3 ) 0.025 Kinematic eddy viscosity ν (m 2 /sec) 10.0
Convergence criteria ε 10 −6
Figure 1: Finite element mesh and definition of computational domain and boundaries
10 10.5 11
E+H(m)
0 20 40 60 80 100
X (m) 0
20 40
60 80
100 Y (m)
X Y
Z
E+H (m) 11 10.9 10.8 10.7 10.6 10.5 10.4 10.3 10.2 10.1 10 9.9
Figure 2: Distribution of initial water depth ( η + h) (Exact solution for estimation of initial water elevation)
If the adjoint boundary condition on the open boundary Γ D is defined as the adjoint flux s x = 0.0 the distribution of the gradient ∂η ∂ J (t
∗0
) at the first iteration is obtained, as shown in Figure 3. The adjoint flux s x is related variable with respect to the flow quantity flux for x direction. In this case, the iterative computation for the inverse analysis is not carried out, since the gradient distribution is not appropriately expressed. The distribution of the initial water elevation is obtained as the linear combination for the gradient vector with respect to the initial water elevation.
Therefore, if the gradient vector at the first iteration is not obtained precisely, the inverse analysis computation cannot be continued. In this case, it appears that if the gradient value is concentrically distributed, the analysis can be carried out. Here, let us consider the case where the adjoint boundary condition on the open boundary Γ D is defined as λ x = 0.0. The adjoint variable λ x is related variable with respect to the flow velocity for x direction.
Consequently, the gradient distribution is obtained as shown in Figure 4. In this case, the iterative computation for the inverse analysis is carried out, and the distribution of the initial water elevation is appropriately obtained.
X (m)
Y(m)
0 50 100
0 10 20 30 40 50 60 70 80 90 100
15 14 13 12 11 10
Figure 3: Distribution of gradient of extended
performance function with respect to initial water
elevation ( conventional method )
X (m)
Y(m)
0 50 100
0 10 20 30 40 50 60 70 80 90 100
10.2 10.1 10 9.9 9.8 9.7
Figure 4: Distribution of gradient of extended performance function with respect to initial water elevation ( present method )
The computational results are shown as follows. Figure 5 shows the variation in the performance function. It can be seen that the performance function decreases and converges. This means the computed water elevation is close to the observed water elevation. The distribution of the initial water elevation, shown in Figure 6, can thus be obtained. Figures 7 and 8 show a comparison of the initial water level between the computed and the exact values on the lines X = 50.0 m and Y = 50.0 m. It can be seen that the computed initial water elevation is lower than the exact initial water elevation.
1e-05 1e-04 0.001 0.01 0.1 1
0 20 40 60 80 100 120
1e-07 1e-06 1e-05 1e-04 0.001 0.01 0.1 1
Relative performance function Difference of relative performance functions
Number of iterations Convergence histories
Performance func. "J"
Difference of "J"
Figure 5: Variation of performance function and absolute value of difference of performance function between (l) and (l + 1) iterations
10 10.5 11
E+H(m)
0 20 40 60 80 100
X (m) 0
20 40
60 80
100 Y (m)
X Y
Z
E+H (m) 11 10.9 10.8 10.7 10.6 10.5 10.4 10.3 10.2 10.1 10 9.9
Figure 6: Distribution of optimized initial water depth ( η + h)
-0.2 0 0.2 0.4 0.6 0.8 1
0 20 40 60 80 100
Water elevation (m)
Distance (m)
Comparison of exact and assimilated initial water elevation Assimilated W.L.
Exact W.L.
Figure 7: Comparison of optimized and exact water elevation (X=50m)
-0.2 0 0.2 0.4 0.6 0.8 1
0 20 40 60 80 100
Water elevation (m)
Distance (m)
Comparison of exact and assimilated initial water elevation Assimilated W.L.
Exact W.L.
Figure 8: Comparison of optimized and exact
water elevation (Y=50m)
Figures 9 - 12 show the time histories of water elevation at the observation points Nos. 1 - 4. The solid line indicates the assimilated water elevation and the broken line indicates the observed water elevation. At all the observation points, the assimilated water elevation is in good agreement with the observed water elevation. It is seen in Figures 9 and 10 that there is difference early on. It appears that this is not a type of the error in the computation of the shallow water flow. This is because the observed data in Figures 9 and 10 are also obtained by using the shallow water code in the same way as the code used in the computation of the estimation for the initial water elevation. There is thus the possibility of a kind of numerical error in the computation of the inverse analysis, since the assimilated data at No. 3 was the same as the data at No. 4.
-0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5
0 2 4 6 8 10 12 14 16 18 20
Water elevation (m)
Time (sec)
Comparison of water elevation at observation point No.1 Assimilated data
Observed data
Figure 9: Comparison of assimilated and observed water elevation at observation point No.1
-0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5
0 2 4 6 8 10 12 14 16 18 20
Water elevation (m)
Time (sec)
Comparison of water elevation at observation point No.2 Assimilated data
Observed data
Figure 10: Comparison of assimilated and observed water elevation at observation point No.2
-0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5
0 2 4 6 8 10 12 14 16 18 20
Water elevation (m)
Time (sec)
Comparison of water elevation at observation point No.3 Assimilated data
Observed data
Figure 11: Comparison of assimilated and observed water elevation at observation point No.3
-0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5
0 2 4 6 8 10 12 14 16 18 20
Water elevation (m)
Time (sec)
Comparison of water elevation at observation point No.4 Assimilated data
Observed data
Figure 12: Comparison of assimilated and observed water elevation at observation point No.4
In this computation, the adjoint boundary condition on
the open boundary Γ D is defined as λ x = 0.0. The adjoint
flux s x integrated in time duration [t 0 ,t f ] and on the open
boundary Γ D must therefore be zero to satisfy the stationary
condition δ J ∗ = 0, under conditions where the performance
function converges. This is because, in the equation of
the first variation of the extended performance function,
the term for the adjoint flux integrated in time and on the
open boundary should be zero. In general, the adjoint
flux on the open boundary is defined as zero to satisfy the
stationary condition. However, in this study, the initial
water elevation can be obtained by computation in which
the adjoint boundary condition on the open boundary is
defined as λ x = 0.0. Therefore, it is necessary to investigate
whether the value of the adjoint flux integrated in time
and on the open boundary is headed for zero if the type
of adjoint boundary condition is changed, i.e., the adjoint
boundary condition on the open boundary is defined as λ x = 0.0. Figure 13 shows the history of the adjoint flux
∫ t
ft
0∫
Γ
Ds x d Γ dt. The value of the adjoint flux is reduced by 92.5% and can be close to zero. The adjoint flux on the open boundary has to theoretically converge to zero to satisfy the stationary condition. However, the adjoint flux on the open boundary does not completely converge to zero, since this computation is carried out numerically and iteratively. In this computation, the adjoint flux on the open boundary can consequently be reduced sufficiently compared with the value at the first iteration, and can approach zero, as shown in Figure 13. It therefore appears that if the data assimilation problem with open boundary can’t be solved based on the traditional method, the present method can be useful in solving for the problem.
0.001 0.01 0.1 1
0 20 40 60 80 100 120
Relative value (adjoint flux)
Number of iterations
Variation of adjoint flux integrated on open boundary Adjoint flux
