Accessible Published by De Gruyter March 16, 2016

High-order relaxation approaches for adjoint-based optimal control problems governed by nonlinear hyperbolic systems of conservation laws

Elimboto M. Yohana and Mapundi K. Banda

Abstract

A computational investigation of optimal control problems which are constrained by hyperbolic systems of conservation laws is presented. The general framework is to employ the adjoint-based optimization to minimize the cost functional of matching-type between the optimal and the target solution. Extension of the numerical schemes to second-order accuracy for systems for the forward and backward problem are applied. In addition a comparative study of two relaxation approaches as solvers for hyperbolic systems is undertaken. In particular optimal control of the 1-D Riemann problem of Euler equations of gas dynamics is studied. The initial values are used as control parameters. The numerical flow obtained by optimal initial conditions matches accurately with observations.

MSC 2010: 49Mxx; 65Mxx

1 Introduction

Conservation laws are time-dependent systems of hyperbolic partial differential equations which describe, for example, the conservation of quantities such as mass, momentum and energy, and are usually nonlinear. In one dimensional space, the equations take the form:

(1.1)tu+xF(u)=0,u(x,t=0)=u0(x)t[0,),x(,)

where u: ℝ× ℝ+ → ℝm is a vector with m conserved quantities uj, and F: ℝm → ℝm is a smooth vector valued flux function, in which each jth component fj(u) is a function of components uj of u. Equations of type (1.1) are popularly referred to as Cauchy problems.

In this paper the computational solutions of optimization problems governed by a system of hyperbolic conservation laws of the form of equation (1.1) will be investigated. As a prototype, the following problem will be considered:

(1.2)minu0J(u(,T),ud);subjecttotu+xF(u)=0u(x,t=0)=u0(x)t[0,T),x(,).

Contrary to the existing results in [7, 34, 42], the main contribution of this paper is the extension of the adjoint method to second-order relaxing schemes, in general, for numerical optimization. In addition the second-order approach will be applied to nonlinear systems using the Euler systems of gas dynamics as a prototype for the first time. In particular, a detailed computational study of the discrete velocity kinetic system formulation [3, 32] of equation (1.1) will be undertaken. The performance of this formulation will be compared with the relaxation approach in [24]. In this approach the hyperbolic systems of conservation laws (HCL) are rewritten as relaxation systems with a stiff source term. For more discussion on relaxation approaches, the reader is referred to [4, 8, 11]. Analysis on the existence and uniqueness of the solution for the relaxation approach described in [24] was given in [43]. Further study on the relaxation approach in [3, 4] is found in [29, 3133]. The relaxation method or the discrete velocity formulation was chosen due to its promising features of simplicity which can lead to generalization to both higher orders and high dimensional systems of HCLs without further modification. The semilinearity structure of the relaxation system allows for Riemannsolvers free treatment and avoids the computation of Jacobians. The relaxation approximation preserves the hyperbolic nature of the system at the expense of additional stiff source terms and additional equations.

HCLs are commonly applied in the optimal control problems [10, 13, 15, 27, 41] – finding some geometry that optimizes performance subject to a set of constraints [20]. Several authors including [23] have employed the nonsmooth optimization in combination with the adjoint methods for subgradient computation. They studied the optimal control of flows with discontinuities and tested the approach using one-dimensional (1D) Riemann problem for Euler equations. In the cases where gradient-based methods were employed, either discontinuities were ignored or means to circumvent their effects were employed. In many situations shocks were smoothed using numerical dissipation. It has been shown that smoothing is sometimes equivalent to modifying the cost function [28]. It is generally known that the semi-group generated by a HCL is not differentiable in L1 even in the scalar, one-dimensional (1D) case [7, 10, 1214] and the solution of HCL is needed in the optimization cycle, it is therefore important to pay close attention to its solution. To circumvent this challenge a notion of shift differentiability is introduced in [10, 12, 13].

In this paper an adjoint-based approach will be considered to solve the optimal control problems numerically. Adjoint and sensitivity calculus based on shift differentiability in the optimal control of entropy solutions of scalar conservation laws with a source term is discussed in [41]. First-order necessary optimality conditions for systems of conservation laws are given in [16]. Discussion on adjoint-based optimization of problems governed by partial differential equations (PDEs) is presented in [20, 21]. Many adjoint-based softwares for CFD have been developed by different pioneers. These include: adjoint-based optimal designs with an application to designing business jets [20], adjoint approach to aerodynamic designs [18, 19], adjoint approach to shape and airfoil designs [1, 17] and continuous adjoint formulation [2, 35]. For trade-o_ between continuous and discrete adjoint approach to automatic aerodynamic optimization, consult [30]. The adjoint approach is robust in the sense that all sensitivities are calculated only once via the adjoint equation in each iteration cycle regardless of the number of control parameters [21].

In combination with the adjoint approach, the relaxation method becomes more appealing due to the simplicity of the relaxation approaches. Other attempts to apply this approach to systems of HCL can be found in [34, 42] where formally first-order schemes were applied. In the current presentation a framework for higher-order extensions is presented. In Section 2 the mathematical framework is summarised. Thereafter the numerical approach is presented in Section 3 which is followed by test cases in Section 4. Conclusions can be found in Section 5.

2 The mathematical formulation

This section presents the mathematical framework and develops the adjoint systems which will be solved using numerical approaches in the next section. The relaxation systems under consideration will be presented. Thereafter the adjoint formulation will be derived.

2.1 Relaxation approaches

The relaxation approach [4, 24] transforms a nonlinear conservation law into a system of semi-linear equations with a nonlinear source term. A good approximation to the original conservation law is achieved by solving relaxation systems for a positive parameter ε ≪ 1. Such relaxation systems are stiff. These systems can be solved numerically by avoiding computationally costly Riemann solvers. Here, two classes of relaxation approaches are considered, namely the relaxation approximation [24] and the discrete kinetic model [3] for the purpose of a comparative study in optimization.

The first relaxation fomulation for the Cauchy problem in equation (1.1) is obtained by replacing (1.1) by a semi-linear system with a stiff lower order term as follows:

(2.1)tuε+xvε=0tvε+A2xuε=1ε(vεF(uε))

where ε is the small positive parameter called the relaxation rate, vε ∈ ℝm is the artificial variable and A = diag{a1, a2,…, am} is a positive diagonal matrix. As, the relaxation system (2.1) formally relaxes to the original HCL. For the solution of equation (2.1) to approach solutions of equation (1.1) as ε → 0, the following sub-characteristic condition needs to be satisfied [24, 25, 43]:

A2F(uε)20uε

where F′(uε) is the Jacobian matrix of the flux function F. In addition, the linear hyperbolic part of equation (2.1) has two characteristic variables:

(2.2)vε  ±  Auε.

The above relaxation system will be referred to as the JXM in the sequel and it is solved subject to the following initial conditions:

uε(x,  t  =  0)  =  u0(x),        vε(x,  t  =  0)  =  F(u0(x))

with u0(x) given in equation (1.1).

Next a different relaxation system presented in [3, 29, 32] namely the discrete kinetic system will be introduced. In this relaxation system, which is simply a BGK model [9], the system of conservation equations in equation (1.1) is approximated by the following semi-linear systems:

(2.3)tfkε  +  λkxfkε  =  1ε(Mk(Pfε)fkε),  k    {1,  .  .  .  ,  L}

where, ε > 0 each fkε    ℝm,  λk are fixed velocities, P is defined by Pf = Σkfk. In addition the following initial conditions will be imposed:

fkε(x,t=0)=Mk(u0(x)).

The Maxwellian functions Mkm → ℝm depend on the quantities uε = Pfε, the flux F and the velocities λk. These functions are continuous and piecewise C1. To link the system in equation (1.1) and (2.3), the following compatibility conditions must be satisfied

(2.4)ΣK=1NMk(u)=u,            ΣK=1NλkMk(u)=F(u)

for all u. A model of equation (2.3) satisfying condition (2.4) can have the following Maxwellian: N = 2, λ1 = −λ, λ2 = λ

M1(u)=12(uF(u)λ),        M2(u)=12(u+F(u)λ).

In the sequel the relaxation system in equation (2.3) will be referred to as the DKS. For more analytical properties of the DKS, the reader may refer to [4, 5]. In general, the Maxwellian functions take the form [5]:

(2.5)Mk(u)=αku+βkF(u),                                      αk,βk    .              

In the next section, the optimization problem to be considered will be introduced in conjuction with the adjoint-based method, which is also discussed in the subsequent sections.

2.2 Adjoint approach to optimization

To solve the problem in equation (1.2) the matching-type objective functional of the form:

(2.6)J(u(x,T;u0);u0ud)=12Ω|u(x,T;u0)ud(x)|2dΩ.

where uis the solution of (1.1) at a terminal time T, u0 is the initial condition and ud is the target (desired) solution. The idea is to solve the equivalent problems which minimise (2.6) subject to either JXM (equation (2.1)) or DKS (equation (2.3)). Hence the initial value u0(x) that will generate an optimal solution u(∙, T)which matches the given target solution ud at terminal time T will be identified. The optimization problem will be solved using an adjoint-type approach.

Firstly, the problem with constraints that are the JXM system will be discussed. Considering optimization of the objective function in equation (2.6) above constrained by JXM (2.1), where u0, the initial condition, is the control variable;

u(x, T; u0) is the solution at time T and ud is the desired profile. This problem can thus be re-written as an unconstrained optimal control problem,

L=L(u(x,T),u0,μ;ud)=J(u(x,T),u0;ud)+0TΩμt[tu+xvtv+A2xu+1ε(vf(u))]dxdt

where μ = (p, q), p, qϵ ℝm is the co-state variable which is assumed to be a smooth function with compact support in Ω and μ = 0on the boundaries of Ω.

In [7] the first-order optimality system was given as: Consider the optimization of the objective function in equation (2.6) above. Note that the superscript ∊ can be dropped for simplicity in notation while the superscript t stands for transpose. To derive the first-order optimality system, the first variations of with respect to each of the variables in μ, u, vand u0 are set equal to zero. Hence one obtains:

(2.7)tuε+xvε=0,uε(x,t=0)=u0ε(x)tvε+A2xuε=1ε(vε-F(uε)),vε(x,t=0)=F(u0ε(x))tpA2xq=F(u)qε,p(x,t=T)=pT(x)tqxq=qε,q(x,t=T)=qT(x).

Setting the partial derivative of .∆u0 to zero, gives the optimality condition

u0J=Ω[p0F(u0)q0ε]dx.

This simplifies to the gradient

(2.8)u0J=p0  +  F(u0)Tq0.

The system (2.7) which comprises the JXM system (2.1) with initial conditions, the adjoint equations with terminal conditions and the gradient (2.8) together form the first-order optimality system [7].

A similar approach is undertaken for the DKS. Consider the Maxwellians of the form in equation (2.5):

(2.9)Mk(u)=αku+βkF(u)=αki=1Nfi+βki=1Nλifi.

Hence for the discrete kinetic model one obtains the Lagrangian of the form:

L=L(u(x,T)(u0),μk;ud)=J(u(x,T)(u0);ud)+k=1N(0TΩμk[tfk+λkxfk1ε(Mk(u)fk)]dxdt).

Introducing the Maxwellians in equation (2.9) as above and re-writing the Lagrangian functional gives

L=J(u(x,T)(u0);ud)+k=1N(0TΩμk[tfk+λkxfk1ε(αki=1Nfi+βki=1Nλififk)]dxdt).

The first-order optimality system for the above Lagrangian takes the form:

(2.10)tfkε+λkxfkε=1ε(Mk(pfε)fkε),k{1,,L}tpkλkxpk=1ε[i=1N(αipi+βiλipi)pi]

with

pi(x,T)=Mi(uT,ud),i=1,,L.

The framework introduced in this section serves as a foundation for the numerical solution of hyperbolic partial differential equation as well as for the optimization process. In the next section, the discretization methods for the relaxations systems (2.1) and (2.3) will be considered, as well as for the derived adjoint systems (2.7) and (2.10).

3 Discretizations of the relaxation systems

In this section discretization approaches for the problems discussed in the previous section will be discussed. In the solution approach, the semi-discrete method in combination with the Implicit–Explicit (IMEX) [7, 37] Runge–Kutta schemes will be applied. The numerical schemes will be tested on one-dimensional nonlinear systems of hyperbolic partial differential equations in the next section.

3.1 Discretization of the systems based on JXM

Let hj denote a grid cell width, with grid spacing hj = xj+1/2xj−1/2, where xj+1/2 = (j+1/2)hj, and a uniform discrete time step, ∆t = tn+1tn for n = 0, 1, 2, . . . . Next approximate Wj+1/2n=W(χj+1/2,tn), and define

Dxwj=wj+1/2wj1/2hj

to obtain the following discretization of the forward system in equation (2.1):

tuj+1hj(vj+1/2vj1/2)=0tvj+1hjA2(uj+1/2uj1/2)=1ε(vjFj)

where

Fj=1hjxj1/2xj+1/2F(u)dx=F(1hjxj1/2xj+1/2udx)+Ο(h2)=F(uj)+O(h2)

see [24] for details. Applying the first-order upwind scheme to the characteristic variables in equation (2.2) gives the following, for each component of u and v:

uj+1/2=12(uj+uj+1)12a1(vj+1vj)vj+1/2=12(vj+vj+1)12a(uj+1uj)

where a can be chosen as max{|a1|, |a2|, …, |am|}, also see Section 4 for specific examples. Similarly, second-order schemes can be derived based on the MUSCL formulation as presented in [24]. Note that here the piecewise linear interpolation to the components of v±Au are considered [24].

An Implicit–Explicit (IMEX) algorithm presented in [7] for the time discretization of the relaxation system (2.1) will be considered. See also, for example, more recent work on Total Variation diminishing (TVD) Runge– Kutta time discretizations construction for relaxation systems in [38, 39].

Following the same discretization as in [7, 24], starting with initial conditions ujn,  vjn=F(uj*) one proceeds as follows:

uj*=ujn,uj(1)=uj*ΔtDxvj*,ujn+1=uj1,vj*=vinΔtε(vj*f(uj*))vj(1)=vj*Δta2Dxuj*vjn+1=vj1.

Following a similar Runge–Kutta time discretization algorithm as for the first-order scheme, one obtains the following second-order time-discretization:

(3.1)u*=un,vn=f(u*)v*=vn+Δtε(v*f(u*))u(1)=u*ΔtDxv*v(1)=v*Δta2Dxu*.

Hence

un+1=12(un+u(2)),              vn+1=12(vn+v(2)).

3.2 Discretization of systems based on DKS

A discretization of a BGK-like model is considered:

tfiε+λixfiε=1ε(Mi(pfε)fiε),i{1,,N}

which is the model based on a kinetic approximation of the problem (1.1) with initial conditions

fε(x,0)=f0ε(x)=Mi(u0)

where Mi are Lipschitz (piecewise C1) continuous functions called Maxwellians, defined on Rm, and other conditions are satisfied as discussed under Section 2.

The first-order spatial discretization of the discrete kinetic model is simple upwinding [3]. For the Runge– Kutta time discretization scheme, the operator splitting approach is applied and split the relaxation system (3.1) into stiff ODE and an advection system. These will be presented below for completeness:

Considering Maxwellian functions Mi(u) = αiu+ βiF(u), Mi(u) ϵ ℝm (i = 1, 2, . . .), one obtains, for the stiff ordinary differential equation (ODE):

uj*=ujn,fj,in=Mj,i(uj*)tfi=1ε(Mi(u)fi)    fj,i*=fj,in+Δtε(Mj,i(uj*)fj,i*).

And for the advection system,

tfi+λixfi=0fj,in+1/2=fj,inλiΔt2hj(fj+1,i*fj1,i*)+|λi|Δt2hj(fj+1,i*2fj,i*+fj1,i*)ujn+1=jfj,in+1/2,fj,in=Mj,i(ujn+1).

For practical reasons, for example, for two-velocities (N = 2) discrete kinetic model, the Maxwellians take the form [3],

Mi(u)=12(u±F(u)λi),              λ1=λ,λ2=λ.

The second-order spatial discretisation is applied as follows [3]: consider Xi = |λ|△t/hj. For positive velocity λi, and for all j ≥ 0,

fj,in+1/2=(1χi)fj,in+χifj1,inhjχi(1χi)2(σj,iσj1,i).

For all j ≥ 0, and for non-positive velocities, λi,

fj,in+1/2=(1χi)fj,in+χifj+1,inhjχi(1Xi)2(σj+1,iσj,i).

A minmod slope limiter is used for all j ≥ 0,

σi,jn=minmod(X1,i,jΔfj+1/2,inhj,X2,i,jΔfj1/2,inhj)

where

minmod(a,b)=sgn(a)+sgn(b)2min(|a|,|b|)

and

Δfj+1/2,in=fj+1,infj,in,              Δfj1/2,in=fj,infj1,in

where X1,i,j and X2,i,j can be chosen according to [3].

To demonstrate the implementation of the temporal discretization a 3-velocities model relaxation system, corresponding to λ1 = − λ2 = λ > 0 and λ0 = 0 will be applied. This was proposed in [5], where f+ corresponds to λ1, f0 corresponds to λ0 and f corresponds to λ2

tf++λxf+=1ε(M+(u)f+)tf0=1ε(M0(u)f0)tfλxf=1ε(M(u)f).

For Runge–Kutta time discretization, the system above is split into two parts: the collision/stiff ODE part

tf+=1ε(M+(u)f+)tf0=1ε(M0(u)f0)tf=1ε(M(u)f)

and the transport system

tf++λxf+=0tf0=0tfλxf=0.

As for the first-order scheme, an IMEX scheme is applied, which takes two stages, one for stiff part and another for advection system, thus, for ujn=uj*,

tf+=1ε(M+(u)f+)tf0=1ε(M0(u)f0)tf=1ε(M(u)f)
fj,+*=fj,+n1ε(Mj,+(uj*)fj,+*)fj,0*=fj,0n1ε(Mj,0(uj*)fj,0*)fj,*=fj,n1ε(Mj,(uj*)fj,*)fj,+(1)=fj,+*ΔtλDxfj,+*=fj,+*χi(fj,+*fj1,+*)hjχi(1χi)2(σj,+*σj1,+*)fj,(1)=fj,*+ΔtλDxfj,*=fj,*+χi(fj+1,*fj,*)hjχi(1χi)2(σj+1,*σj,*)fj,+**=fj,+(1)+Δtε(Mj,+(uj**)fj,+**)+2Δtε(Mj,+(uj*)fj,+*)fj,0**=fj,0(1)+Δtε(Mj,0(uj**)fj,0**)+2Δtε(Mj,0(uj*)fj,0*)fj,**=fj,(1)+Δtε(Mj,(uj**)fj,+**)+2Δtε(Mj,(uj*)fj,*)fj,+(2)=fj,+**ΔtλDxfj,+**=fj,+**χi(fj,+**fj1,+**)hjχi(1χi)2(σj,+**σj1,+**)fj,(2)=fj,**+ΔtλDxfj,**=fj,**+χi(fj+1**fj,**)hjχi(1χi)2(σj+1,**σj,**)fj,+n+1=12(fj,+(2)+fj,+n),fj,n+1=12(fj,(2)+fj,n),fj,0n+1=fj,0**

where

ujn+1=12i(fj,i(2)+fj,in).

In the following section, the discrete version of the adjoint relaxation system in (2.7) will be derived. Analogously, the discretization process for both space and time is achieved separately, and then, the two semi-discrete schemes are merged to obtain a fully discrete scheme.

3.3 Discretization of the adjoint systems

A similar approach for the discretization of the adjoint system as above will be applied, following a similar approach to [7]. The adjoint system is solved backwards in time.

The set of adjoint equations (2.7) is considered.

As shown in [7], the characteristic variables p± Aq satisfy

t(p±Aq)±A2x(p±Aq)=0t(p±A(q))±A2χ(p±A(q))=0.

The adjoint equation is solved backwards in time, thus an upwind discretization for each component of the linear system advects p ± aq and −p ± a(−q) with velocity ∓a. Therefore,

pj+1/2=pj+pj+12aqj+1qj2    qj+1/2=qj+qj+12pj+1pj2a.

The time discretization of the adjoint system takes the form [7]:

pj(1)=pjn+1,pj*=pj(1)+Δtεqj*fj(uj*)Δta2Dxqj1,pjn=pj*,qj(1)=qjn+1qj*=qj(1)Δtεqj*ΔtDxpj(1)qjn=qj*.

For the second-order approach an MUSCL method is applied. This was first developed by Sweby [40]. Schemes developed using this approach satisfy the TVD property and are utmost second-order accurate.

Similarly, the second-order in time and space discretization is similar to the first order scheme, see [7].

A second-order accurate discretization of the adjoint equations in space can be summarised as follows:

pj+1/2=12(pj+pj+1+a(qj+1qj)+12(σjσj+1+)).

Similarly,

qj+1/2=12(qj+1+qj)+12a(pj+pj+1)14a(σj+1++σj).

Replacing p by −p and q by −q, one obtains:

(3.2)pj+1/2=12(pjpj+1+a(qj+1+qj)+12(σjσj+1+))
(3.3)qj+1/2=12(qj+1qj)+12a(pjpj+1)14a(σj+1++σj)

where

σj=(pj+1aqj+1pj+aqj)φ(ϑj)σj+1+=(pj+2+aqj+2pj+1aqj+1)φ(ϑj+1+).

Futhermore,

φ(ϑj±)=φ(pj±aqj(pj1±aqj1)pj+1±aqj+1(pj±aqj)).

To conform with the format of the adjoint equations, equations (3.2) and (3.3) are written as

pj+1/2=[12(pj+pj+1)a2(qjqj+1)14(σjσj+1+)]qj+1/2=[12(qj+1+qj)12a(pjpj+1)+14a(σj+1++σj)].

The second-order TVD Runge–Kutta time discretization takes the form (also see [7]):

pj(2)=12pjn+1,qj(2)=12qjn+1qj**=(εε+Δt)qj(2)+(εε+Δt)Δt2hj[(pj+1(2)pj+1(2))+a(qj+1(2)2qj(2)+qj1(2))]+(εε+Δt)Δt4hj[σj1(σj++σj)+σj1+]pj**=pj(2)+qj**Δtεf(uj**)+a2Δt2hj[(qj+1(2)qj1(2))+1a(pj+1(2)2pj(2)+pj1(2))]aΔt4hj[σj1+(σj+σj)σj+1+]pj(1)=pj**,qj(1)=qj**qj*=(εε+Δt)qj(1)(εε+Δt)2Δtεqj**+(εε+Δt)Δt2hj[(pj+1(1)pj1(1))+a(qj+1(1)2qj(1)+qj1(1))]+(εε+Δt)Δt4hj[σj1+(σj++σj)+σj+1+]pj*=pj(1)Δtεf(uj**)(qj*2qj**)+a2Δt2hj[(qj+1(1)qj1(1))+1a(pj+1(1)2pj(1)+pj1(1))]aΔt4hj[σj1+(σj+σj)σj+1+]pjn=12pjn+1+pj*,qjn=12qjn+1+qj*.

Without loss of generality, the discretization of the adjoint system (2.10) is considered for a 3-velocity system:

tpiλixpi=1ε(Mi(u)fi),    i  =  1,  .  .  .  ,  3.  

If the system is written in an extended form

(3.4)tp++λxp+=1ε(M+(u)f+)tp0=1ε(M0(u)f0)tpλxp=1ε(M(u)f).

Thus, second-order Runge–Kutta time discretization for the adjoint system (3.4) would be

pj,+(2)=12pj,+n+1,pj,(2)=12pj,n+1pj,+**=pj,+(2)Δtε(Mj,+(uj)fj,+)λΔtDxpj,+(2)pj,**=pj,(2)+Δtε(Mj,(uj)fj,)+λΔtDxpj,(2)pj,+(1)=pj,+**,pj,(1)=pj,**pj,+*=pj,+(1)1ε(Mj,+(uj)fj,+)+2Δtε(Mj,+(uj)fj,+)λΔtDxpj,+(1)pj,*=pj,(1)Δtε(Mj,(uj)fj,)+2Δtε(Mj,(uj)fj,)+λΔtDxpj,(1)pj,+n=12(pj,+n+1+pj,+*),pj,n12(pj,n+1+pj,*).

3.4 Algorithm for gradient computing

The optimal control method under consideration involves the computation of the cost functional (2.6) gradient. This gradient is used to modify the design variable (in this case initial condition) which adjusts in such a way to produce an optimal solution that matches the given target. To compute the gradient of (2.6) and incorporate it in the optimal control method, the following algorithm can be applied.

Algorithm 3.1. The following is the algorithm to minimize the functional:

  1. Consider the constraint in relaxation form (2.1) and (2.3);

  2. Derive its adjoint system;

  3. Solve for the flow variable u(·, T) forward in time;

  4. Use the solution of the flow variable above at the terminal time T to solve for the adjoint variable backward in time;

  5. Use the adjoint variable and the control variable u0 to evaluate the gradient of the cost functional (2.6), ∇u0J;

  6. Update the control variable by using the gradient obtained above and the chosen step size by making a step in the negative gradient direction,

    αnew=αoldu0J;
  7. Repeat steps 3 to 6 until minimization is reached.

This algorithm was applied to optimize the Euler equations of gas dynamics. To ensure convergence of the optimization process the Armijo step-size rule described in [6] was considered to automatically choose the optimal value of α while making steps towards the steepest descent direction.

4 Numerical results

This Section presents a comparative study of the second-order numerical results based on implementation of the two relaxing schemes derived under Section 3. The main focus was to investigate systems of Euler equations by first computing their solutions and then employing adjoint-based control summarized by Algorithm 3.1 as an optimal control strategy. Two relaxing schemes: the JXM scheme described in [24] and the discrete kinetic scheme (DKS) presented in [3, 32] are compared. Schemes applied are based on the method of lines associated with Runge–Kutta type time discretizations which are total variation diminishing (TVD). Second-order accuracy in solution is achieved through second-order MUSCL type space discretization coupled with a second-order Runge–Kutta time splitting scheme.

All numerical results that will be presented in the sequel were performed with Intel core i5 processor with 2.67 GHz, and 4 GB RAM, and the programs are developed using the Python scripting language.

4.1 Numerical discretizations of spatial and temporal domains

To discretize both the time and spatial domains, consider a bounded domain of ℝ, [xL, xR]. For the sake of simplicity, the domain [xL, xR] is divided uniformly into a sequence of M + 1 gridpoints, m = 0, ..., M + 1, such that x0 = xL, xM = xR; with a mesh size ∆x = 1/(M − 1) each and gridcell width ∆x = Xj+1/2xj−1/2. The temporal domain is considered to be [0, T], discretized into N time levels, tn, with time step, ∆t = tn+1tn, and the time horizon is given by T = Nt. In the spirit of finite volumes, the approximate solution based on cell averages u(xj, tn) will be denoted by ujn.

The approximate solutions were computed on a uniform mesh, with 200, 400, 800, and 1600 gridpoints and present results for different values of time. The time step is fixed at ∆t selected based on the CFL condition as ∆t = CFL∆x/amax for stability, where amax is the maximum of the characteristic speeds. Without loss of generality, solutions are computed on the spatial computational domain [0, 1], and a mesh with M + 1 grid points. Where possible, numerical solutions are compared with exact solutions or a reference solution solved with the number of spatial grid points M = 1600.

For the relaxing scheme DKS, the macroscopic variable u is linked to microscopic fi by the Maxwellians given in equation (2.5) where Maxwellians must be monotone preserving according to [3].

Next, the following variables related to DKS, that will be used throughout the course of this section are defined: for a system with two velocity models λ1 = −λ2, while for three velocities discrete kinetic scheme, λ3 = −λ1, where λ2 is set to 0. Furthermore, α = α1 = α3, -α2 = 1 − 2α, β = β2, hence

                β1=12(1λ3β),        β3=12(1λ3β).

For two-velocity models, obviously the diagonal relation

M1(u)=12(u+F(u)λ1),        M2(u)=12(u+F(u)λ2)

between macroscopic u and Maxwellians is obtained. Here, numerical results for discrete kinetic model are restricted to two or three velocities schemes. A relaxation rate ε = 10−8 is considered for both schemes.

In the sequel, firstly numerical results obtained with the schemes as derived in Section 3 will be presented for various gas dynamic problems. The aim is to demonstrate the performance of these accurate numerical schemes to be incorporated in adjoint-based optimal control process.

4.2 Numerical solutions for systems of equations

Here, the solution of one-dimensional system of equations will be considered, and in particular, these solutions will act as an illustration of the optimization process under consideration.

Consider the approximate numerical solution of the one-dimensional Euler equations of gas dynamics in a conserved form

(4.1)          tρ+xm=0tm+x(ρu2+p)=0tE+x(u(E+p))=0

where ρ, u, m = ρu, p, and E are, respectively, the density, velocity, momentum, pressure, and total energy of the gas. For a perfect gas, E is related to other quantities by

E=pγ1+12ρu2

where γ= 1.4 constitutes the thermodynamic property of the gas and is the ratio of specific heat constants.

The Jacobian of the fluxes of (4.1), therefore, has three eigenvalues given by

λ=uc,λ0=u,λ+=u+c

where c=γp/p=γRT is the local sound speed, R the universal gas constant and T is the temperature of the gas. For computational tests, the choices according to [24], i.e., a1, a2, a3 = sup|uc|, sup|u|, sup|u + c|, respectively, or a1 = a2 = a3 = max(sup|uc|, sup|u|, sup|u + c|) can be applied. The two relaxation systems (2.1) and (2.3) with appropriate parameters will be solved.

4.3 Sod shock tube problem

The computational tests involved experimentation with the Sod’s data for the Shock Tube problem. Detailed descriptions of the shock tube problem is found in [22, 26].

The schemes derived in Section 3 were implemented with the Riemann data defined in this way: uL = (ρL, vL, ρL) corresponding to the density, velocity and pressure on the left part of the domain, i.e., 0 ≤ x < 0.5 and uR = (ρR, vR, ρR) is the data corresponding to density, velocity and pressure on the right part of the domain, i.e., 0.5 ≤ x ≤ 1. All tests (unless stated otherwise) for the system of the Euler equations are carried out for the following set of values: T = 0.17, the computational space domain is [0, 1], M = 400 grid-points, CFL = 0.75 and ε = 10−8. In addition, for the second-order relaxing schemes the minmod slope limiter was used.

The solution of the Euler equations was computed with the following pair of Riemann data:

uL=  (1.25,  0.0,  1.2),      uR=  (0.25,  0.0,  0.25).

For this a1 = 1.0, a2 = 2.5, a3 = 5.2 is considered. Results obtained with the first and second order relaxing schemes are presented in Fig. 1. In this figure it can be observed that the two schemes agree very well with each other.

Figure 1 First and second order numerical solutions for density, velocity and pressure profiles with JXM and DKS, for 1D Euler equations at time, T = 0.17, uL = (1.25, 0.0, 1.2), uR = (0.25, 0.0, 0.25), for a1 = 1.0, a2 = 2.5, a3 = 5.2, M = 400 and ε = 10−8.

Figure 1

First and second order numerical solutions for density, velocity and pressure profiles with JXM and DKS, for 1D Euler equations at time, T = 0.17, uL = (1.25, 0.0, 1.2), uR = (0.25, 0.0, 0.25), for a1 = 1.0, a2 = 2.5, a3 = 5.2, M = 400 and ε = 10−8.

The second example is solved with initial condition uL = (1.45, 0.0, 1.5), uR = (0.45, 0.0, 0.5). This time a1 = 2.2, a2 = 2.5, a3 = 5.0 is chosen. Solutions computed over a time length T = 0.17 are displayed in Fig. 2. Density, velocity and pressure profiles for both schemes are reasonably equivalent.

Figure 2 First and second order numerical solutions for density, velocity and pressure profiles with JXM and DKS, for 1D Euler equations at time, T = 0.17, uL = (1.45, 0.0, 1.5), uR = (0.45, 0.0, 0.5); a1 = 2.2, a2 = 2.5, a3 = 5.0, M = 400 and ε = 10−8.

Figure 2

First and second order numerical solutions for density, velocity and pressure profiles with JXM and DKS, for 1D Euler equations at time, T = 0.17, uL = (1.45, 0.0, 1.5), uR = (0.45, 0.0, 0.5); a1 = 2.2, a2 = 2.5, a3 = 5.0, M = 400 and ε = 10−8.

Finally, the last example evolved from initial data uL = (2.5, 0.0, 2.0), uR = (0.5, 0.0, 0.6) computed over same time horizon T = 0.17 as for previous examples, and the present simulated results are displayed in Fig. 3. Simulations show that, solutions obtained with two different relaxing schemes are comparable.

Figure 3 First and second order numerical solutions for density, velocity and pressure profiles with JXM and DKS, for 1D Euler equations at time, T = 0.17, uL = (2.5, 0.0, 2.0), uR = (0.5, 0.0, 0.6); a1 = 3.5, a2= 4.5, a3 = 5.5, M = 400 and ε = 10−8.

Figure 3

First and second order numerical solutions for density, velocity and pressure profiles with JXM and DKS, for 1D Euler equations at time, T = 0.17, uL = (2.5, 0.0, 2.0), uR = (0.5, 0.0, 0.6); a1 = 3.5, a2= 4.5, a3 = 5.5, M = 400 and ε = 10−8.

These examples demonstrate that the two relaxing schemes under consideration are basically equivalent and give similar results in one-dimension as expected.

All numerical results for the two relaxing schemes presented for different meshes show that the two schemes, the JXM and the DKS give similar results. The meshes of 200, 400, 800, 1600 grid points were considered for the sake of grid convergence. Clearly, results show schemes convergence with mesh refinement, but beyond 400 grid points, numerical experiments reveal that refinement is no longer necessary as solutions are already sharp enough at this point (see Figs. 4 and 5). Thus, only numerical results for 400 grid points are presented, and this grid will be used for optimization as well.

Figure 4 Comparing first-order and second-order solutions for different grids with JXM scheme, T = 0.17, uL = (2.5, 0.0, 2.0), uR = (0.5, 0.0, 0.6); a1 = 2.0, a2 = 3.5, a3 = 5.5, ε = 10−8.

Figure 4

Comparing first-order and second-order solutions for different grids with JXM scheme, T = 0.17, uL = (2.5, 0.0, 2.0), uR = (0.5, 0.0, 0.6); a1 = 2.0, a2 = 3.5, a3 = 5.5, ε = 10−8.

Figure 5 Comparing first-order and second-order solutions for different grids with DKS. Same parameters as in Fig. 4.

Figure 5

Comparing first-order and second-order solutions for different grids with DKS. Same parameters as in Fig. 4.

It needs to be noted that some sets of data introduce more viscosity, see for example in Fig. 3, considerably smearing contact discontinuities while others produce sharper solutions even for first-order schemes especially with increased grid refinement. In addition, these numerical results further verify the applicability and the robustness of the methods considered.

4.4 Adjoint-optimization tests

The optimal control is carried out by matching the numerical solution to the target for a given time length as discussed in Algorithm 3.1. The minimization process, therefore, involves the cost functional of L2 norm as defined in (2.6). In addition, the Wolfe conditions [6, 36] to restrict the choice of the step-size which is used to modify the functional gradient (2.8) that perturb the flow solution is employed. The control parameters are chosen to be initial values of the primitive variables: density, velocity and pressure. Existence of the optimal control solution can be conceptualized from numerical results. A good match between the numerical solution and the desired one for both first and second order schemes was obtained.

The process is initiated by choosing some initial guess u0, and then using the relaxing schemes derived under Section 3 to solve the hyperbolic conservation laws in order to obtain the solution u(x, T), a function of u0 at some terminal time, T. The functional gradient (2.8) is calculated and used at every optimization cycle to modify the design parameter u0, and each modified u0 is evolved by the relaxing scheme until the optimal solution u(x, T) that matches the target is attained. The initial solution, denoted V(x, T) (u0), is obtained by solving the flow equation once forward in time, where V is any primitive quantity (density, velocity and pressure for the Euler system of equations). Optimal denotes the numerical solution that matches the desired profile (solution) and the desired profile is called Target.

The stopping criterion used for the functional (2.6) is defined using absolute values, |J(u0, ud)| < 10−4.

In the first example, the initial guess for the target solution comprises the Sod shock tube data for the system of Euler equations (4.1),

uL=  (1.0,  0.0,  1.0)  ,    uR=  (0.125,  0.0,  0.1).

The target solution for the initial data (4.2) is obtained using a1 = 1.0, a2 = 1.68, a3 = 5.045 over T = 0.17. The optimal control results presented are first-order and second-order in time and space for the two relaxing schemes presented in Section 3.

The first example considered matches the target to the optimal solution computed from the set of data

(4.2)uL=  (1.2,  0.078,  1.2),    uR=  (0.325,  0.285,  0.295).

a1 = 2.183, a2 = 3.004, a3 = 4.286 over the time T = 0.17. The optimal solution for density, pressure and velocity were found to be in a very good agreement with the target for both schemes. First and second order optimal control results for JXM scheme and DKS are, respectively, displayed in Figs. 6 and 7.

Figure 6 First-order (Left) and second-order (Right) 1D Euler equations optimal control results. The optimal solution is shown by the red solid line and target by black solid line with squares at data points, T = 0.17, uL = (1.2, 0.078, 1.2), uR = (0.325, 0.285, 0.295), a1 = 2.183, a2 = 3.004, a3 = 4.286,M= 400 and ε = 10−8.

Figure 6

First-order (Left) and second-order (Right) 1D Euler equations optimal control results. The optimal solution is shown by the red solid line and target by black solid line with squares at data points, T = 0.17, uL = (1.2, 0.078, 1.2), uR = (0.325, 0.285, 0.295), a1 = 2.183, a2 = 3.004, a3 = 4.286,M= 400 and ε = 10−8.

Figure 7 First-order (Left) and second-order (Right) 1D Euler equations optimal control results obtained with DKS. Optimal solution is depicted by the red solid line and target by black solid line with squares at data points, same parameters as used for Fig. 6.

Figure 7

First-order (Left) and second-order (Right) 1D Euler equations optimal control results obtained with DKS. Optimal solution is depicted by the red solid line and target by black solid line with squares at data points, same parameters as used for Fig. 6.

The second example considered in which also the optimal solutions for the density, velocity and pressure match accurately with the target solutions. Optimal solutions are solved for values of initial condition,

uL=  (1.2,  0.2,  1.25),      uR=  (0.32,  0.73,  0.32)

a1 = 2.47, a2 = 3.5, a3 = 4.36. The computations were carried out over usual time T = 0.17. First and second order results for this set of data are displayed in Fig. 8 for JXM scheme and in Fig. 9 for DKS.

Thirdly, results for initial Riemann data

uL=  (1.24,  0.1852,  1.25),    uR=  (0.366,  0.629,  0.33)
Figure 8 First-order (Left) and second-order (Right) 1D Euler equations optimal control results. Optimal (red solid line) and target (black solid line with squares at data points) solutions, for density, velocity and pressure profiles with JXM scheme at time, T = 0.17, uL = (1.2, 0.2, 1.25), uR = (0.32, 0.73, 0.32), a1 = 2.47, a2 = 3.5, a3 = 4.36,M= 400 and ε = 10−8.

Figure 8

First-order (Left) and second-order (Right) 1D Euler equations optimal control results. Optimal (red solid line) and target (black solid line with squares at data points) solutions, for density, velocity and pressure profiles with JXM scheme at time, T = 0.17, uL = (1.2, 0.2, 1.25), uR = (0.32, 0.73, 0.32), a1 = 2.47, a2 = 3.5, a3 = 4.36,M= 400 and ε = 10−8.

Figure 9 First-order (Left) and second-order (Right) 1D Euler equations optimal control results. Optimal (red solid line) and target (black solid line with squares at data points) solutions, for density, velocity and pressure profiles with DKS. Same parameters as in Fig. 8.

Figure 9

First-order (Left) and second-order (Right) 1D Euler equations optimal control results. Optimal (red solid line) and target (black solid line with squares at data points) solutions, for density, velocity and pressure profiles with DKS. Same parameters as in Fig. 8.

are presented. The following parameters are chosen: a1 = 1.96, a2 = 2.9, a3 = 4.33. As for the previous examples, a good match between the optimal and target solutions is obtained. Results are visualized in Figs. 10 and 11.

Finally, an example presented in [23] is considered. In this case, initial optimal solutions are solved for values of Sod data initial condition,

uL=  (1.0,  0.0,  1.0)  ,    uR=  (0.125,  0.0,  0.1).
Figure 10 First-order (Left) and second-order (Right) 1D Euler equations optimal control results. Optimal (red solid line) and target (black solid line with squares at data points) solutions, for density, velocity and pressure profiles with JXM scheme at time, T = 0.17, uL = (1.24, 0.1852, 1.25), uR = (0.366, 0.629, 0.33), a1 = 1.96, a2 = 2.9, a3 = 4.33,M= 400 and ε = 10−8.

Figure 10

First-order (Left) and second-order (Right) 1D Euler equations optimal control results. Optimal (red solid line) and target (black solid line with squares at data points) solutions, for density, velocity and pressure profiles with JXM scheme at time, T = 0.17, uL = (1.24, 0.1852, 1.25), uR = (0.366, 0.629, 0.33), a1 = 1.96, a2 = 2.9, a3 = 4.33,M= 400 and ε = 10−8.

Figure 11 First-order (Left) and second-order (Right) 1D Euler equations optimal control results. Optimal (red solid line) and target (black solid line with squares at data points) solutions, for density, velocity and pressure profiles obtained with the DKS. Same parameters as in Fig. 10.

Figure 11

First-order (Left) and second-order (Right) 1D Euler equations optimal control results. Optimal (red solid line) and target (black solid line with squares at data points) solutions, for density, velocity and pressure profiles obtained with the DKS. Same parameters as in Fig. 10.

and the target is the solution of the initial conditions

uL=  (1.1,  0.0,  1.1),    uR=  (0.2,  0.0,  0.2)

for a1 = 1.6, a2 = 2.82, a3 = 4, 25. Computations are carried out over usual time T = 0.17. First and second order results for this example are given in Fig. 12 for JXM scheme and in Fig. 13 for DKS. A good match between optimal and target solutions was observed.

Figure 12 First-order (Left) and second-order (Right) 1D Euler equations optimal control results. Optimal (red solid line) and target (black solid line with squares at data points) solutions, for density, velocity and pressure profiles with JXM scheme at time, T = 0.17, uL = (1.0, 0.0, 1.0), uR = (0.1, 0.0, 0.125) for optimal, and uL = (1.1, 0.0, 1.1), uR = (0.2, 0.0, 0.2) for target. a1 = 1.96, a2 = 2.9, a3 = 4.33,M= 400 and ε = 10−8.

Figure 12

First-order (Left) and second-order (Right) 1D Euler equations optimal control results. Optimal (red solid line) and target (black solid line with squares at data points) solutions, for density, velocity and pressure profiles with JXM scheme at time, T = 0.17, uL = (1.0, 0.0, 1.0), uR = (0.1, 0.0, 0.125) for optimal, and uL = (1.1, 0.0, 1.1), uR = (0.2, 0.0, 0.2) for target. a1 = 1.96, a2 = 2.9, a3 = 4.33,M= 400 and ε = 10−8.

Figure 13 First-order (Left) and second-order (Right) 1D Euler equations optimal control results. Optimal (red solid line) and target (black solid line with squares at data points) solutions, for density, velocity and pressure profiles with DKS at time, T = 0.17. Same parameters as in Fig. 12.

Figure 13

First-order (Left) and second-order (Right) 1D Euler equations optimal control results. Optimal (red solid line) and target (black solid line with squares at data points) solutions, for density, velocity and pressure profiles with DKS at time, T = 0.17. Same parameters as in Fig. 12.

Figure 14 Convergence history for the solution of the optimization problem computed with the first order (left) and second-order (right) relaxing schemes for both JXM and discrete kinetic schemes DKS.

Figure 14

Convergence history for the solution of the optimization problem computed with the first order (left) and second-order (right) relaxing schemes for both JXM and discrete kinetic schemes DKS.

In conclusion it can be pointed out that the numerical approach based on the two relaxation formulation give very accurate results for the systems of conservation laws. This is achieved while the details of the flow are under-resolved. There is no need for following the shocks along the characteristics or such technical computations.

4.4.1 Functional Convergence

The optimal control results presented above will be further analysed by giving a brief convergence history for the two relaxing schemes. Under this analysis, the representative optimal control example associated with a set of data

uL=  (1.2,  0.078,  1.2),    uR=  (0.325,  0.285,  0.295).

whose computations were carried out over time T = 0.17 is chosen. The two graphs below summarize the progressive minimization of the cost functional with the number of iterations for both the first-order and second order relaxing schemes. Results show that second-order schemes need fewer iterations to converge than first-order ones. However, it must also be pointed out that the number of optimization iterations is independent of the grid size.

4.4.2 Comparison of computation time

Besides qualitatively and physically comparing optimal control results for the JXM scheme and the DKS discussed in the previous sections, a brief discussion of the computation time taken for simulation of these results will follow. The time needed for the JXM scheme to converge is larger than that needed for the DKS. Obviously, time taken for the algorithm to converge for both schemes increases with the number of discretization points, M. The computation time for a representative example is reported in Tables 1 and 2 for the space discretization with M= 200, 400, 800 grid points against the number of iterations, NI. All computations are performed on a 2.67 GHz Intel Core dual i5 processor using Python 2.5.6.

Table 1

Computational time for first-order JXM and DKS.

MNIJXMDKS
200221.8245× 1021.26 × 102
400229.07168 × 1024.9572 × 102
800221.9837 × 1031.2961 × 103
Table 2

Computational time for second-order JXM and DKS.

MNIJXMDKS
200154.3234 × 1022.97 × 102
400151.7576 × 1031.1453 × 103
800153.2761 × 1032.370 × 103

5 Conclusions

In this work we extended the adjoint based framework to second-order accurate solution approaches for the hyperbolic systems of conservation laws. Some theoretical perspectives from previous works have been substantiated by the numerical results obtained.

Research focused on the adjoint approach to optimize a problem constrained by nonlinear systems of Hyperbolic Conservation Laws. The optimality conditions were derived giving the optimality systems for the two relaxation approaches discussed. Thereafter their corresponding relaxing schemes of both first and second order in time and space were derived. Contrary to the existing results, the adjoint scheme derived in [7] was extended to systems and up to second-order optimal control results were obtained. Through generalization new adjoint relaxing schemes were developed also for the discrete-velocity kinetic model. The coupling information obtained by solving flow and adjoint systems during optimization cycle for both relaxation approaches, demonstrate that the first and second order numerical results obtained are promising and comparable.

From the point of view of efficiency, computer runtimes for both schemes, of the flow equations and that of the adjoint equations are almost equal. But in general the discrete-velocity kinetic scheme tends to be more efficient in terms of runtime. The storage requirements are almost equal since each scheme fundamentally handles and processes almost the same amount of data during computations.

The hyperbolic systems of equations develop discontinuities due to the interaction of characteristics even if their initial conditions are smooth. As a result solving the backward problem seems to be more sensitive. Information is hardly reversible for interesting problems. Thus stability conditions need to be accurately satisfied at every time step especially for the adjoint problem.

For systems of hyperbolic conservation laws, the Euler equations of gas dynamics have been considered. It is expected that these results can be extended to other nonlinear systems of hyperbolic conservation laws. This study could have contributions in these areas of applications and also presents challenges to researchers and scientists for further investigations.

It can, therefore, be concluded that the study lays down a foundation for possible application of adjoint optimization schemes to other application areas. Areas of further study may include incorporation of schemes higher than second-order in optimization processes. It is also possible to extend this work to deal with multi-dimensional problems. More importantly, higher than first-order optimality conditions for adjoint based optimization are still elusive.

Acknowledgment

EMY acknowledges financial support from the African Institute of Mathematical Sciences (AIMS), University of the Witwatersrand (Wits) and University of Tanzania for granting study leave to pursue this work.

Funding

This work is based on the research supported in part by the National Research Foundation (NRF) of South Africa UID: 74260, UID: 81299 and UID: 85566.

References

[1] W. K. Anderson and D. L. Bonhaus, Airfoil design on unstructured grids for turbulent flows, AIAA37(1999), 185–191. Search in Google Scholar

[2] W. K. Anderson and V. Venkatakrishnan, Aerodynamic design optimization on unstructured grids with a continuous adjoint formulation, AIAA97(1997) 0643. Search in Google Scholar

[3] D. Aregba-Driollet and V. Milisic, Kinetic approximation of a boundary value problem for conservation laws, J. Num. Math. 97(2004), 595–633. Search in Google Scholar

[4] D. Aregba-Driollet and R. Natalini, Convergence of relaxation schemes for conservation laws, Applicable Analysis61(1996), 163–193. Search in Google Scholar

[5] D. Aregba-Driollet and R. Natalini, Discrete kinetic schemes for multidimensional conservation laws, SIAM J. Numer. Anal. 37(2000), 1973–2004. Search in Google Scholar

[6] L. Armijo, Minimization of functions having Lipschitz continuous first partial derivatives, Pacific J. Math16(1966), 1–3. Search in Google Scholar

[7] M. K. Banda and M. Herty, Adjoint IMEX-based schemes for control problems governed by hyperbolic conservation laws, Comput Optim Appl. (2010). Search in Google Scholar

[8] M. K. Banda and M. Seaïd, Higher-order relaxation schemes for hyperbolic systems of conservation laws, J. Numer. Math. 13(2005) 171–196. Search in Google Scholar

[9] P. Bhatnagar, E. Gross, and M. Krook, A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems, Phys. Rev. 94(1954), 511–525. Search in Google Scholar

[10] S. Bianchini, On the shift differentiability of the flow generated by a hyperbolic system of conservation laws, Discrete Contin. Dynam. Systems6(2000), 329–350. Search in Google Scholar

[11] F. Bouchut and T. Morales de Luna, Semi-discrete entropy satisfying approximate Riemann solvers. The case of the Suliciu relaxation approximation, J. Sci. Comput. 41(2009), 483–509. Search in Google Scholar

[12] A. Bressan and G. Guerra, Shift-differentiability of the flow generated by a conservation law, Discrete Contin. Dynam. Systems3(1997), 35–58. Search in Google Scholar

[13] A. Bressan and M. Lewicka, Shift differentials of maps in BV spaces, In: Nonlinear Theory of Generalized Functions, Vienna, 1997. Chapman and Hall/CRC Res. Notes Math. 401, Chapman and Hall/CRC, Boca Raton, FL (1999), 47–61. Search in Google Scholar

[14] A. Bressan and A. Marson, A maximum principle for optimally controlled systems of conservation laws, Rend. Sem. Univ. Podova94(1995), 79–94. Search in Google Scholar

[15] A. Bressan and A. Marson, A variational calculus for discontinuous solutions to conservation laws, Commun. Partial Di_er. Eqn. 20(1995), 1491–1552. Search in Google Scholar

[16] A. Bressan and W. Shen, Optimality conditions for solutions to hyperbolic balance laws. Control methods in PDE dynamical systems, Contemp. Math. 426(2007), 129–152. Search in Google Scholar

[17] H. Cabuk, C.H. Sung and V. Modi, Adjoint operator approach to shape design for internal incompressible flow, In: Pennyslvania State Univ., 3rd International Conference on Inverse Design Concepts and Optimization in Engineering Sciences (1991). Search in Google Scholar

[18] J. Elliot, Aerodynamic optimization based on the Euler and Navier-Stokes equations using unstructured grids, PhD Thesis, MIT, Department of Aeronautics and Astronomy, 1998. Search in Google Scholar

[19] J. Elliot and J. Peraire, Practical 3D aerodynamic design and optimization using unstructured meshes, AIAA J. 35(1997), 1479–1485. Search in Google Scholar

[20] M. B. Giles and N. A. Pierce, An introduction to the Adjoint Approach to design, Flow, Turbul. Combust. 65(2000), 393– 415. Search in Google Scholar

[21] M. D. Gunzburger, Perspectives in Flow Control and Optimization, Advances in design and control, SIAM, PA. USA, 2003. Search in Google Scholar

[22] C. Hirsch, Numerical Computation of Internal and External Flows, John Wiley and Sons, 1988. Search in Google Scholar

[23] C. Homescu and I. M. Navon, Optimal control of flow with discontinuities, J. Comp. Phys. 187(2003), 660–682. Search in Google Scholar

[24] S. Jin and Z. Xin, The relaxation schemes for systems of conservation laws in arbitrary space dimensions, Comm. Pure Appl. Math48(1995), 235–276. Search in Google Scholar

[25] C. Lattanzio and D. Serre, Convergence of a relaxation scheme for hyperbolic systems of conservation laws, Numer. Math. 88(2001), 121–134. Search in Google Scholar

[26] R. J. LeVeque, Numerical Methods for Conservation Laws, Birkhäuser Verlag, Berlin, 1992. Search in Google Scholar

[27] Z. Liu and A. Sandu, On the properties of discrete adjoints of numerical methods for the advection equation, Int. J. for Num. Meth. in Fluids56(2008), 769–803. Search in Google Scholar

[28] T. Matsuzawa and M. Hafez, Optimum shape design using adjoint equations for compressible flows with shock waves, International Journal for Computational Fluid Dynamics7(1998), 343–365. E. M. Yohana and M. K. Banda, Relaxation approaches for adjoint control Ë 71 Search in Google Scholar

[29] V. Milisic, Stability and convergence of discrete kinetic approximations to an initial-boundary value problem for conservation laws, Proceedings of the Amer. Math. Soc. 97(2004), 595–633. Search in Google Scholar

[30] S. K. Nadarajah and A. Jameson, A Comparison of the Continuous and Discrete Adjoint Approach to Automatic Aerodynamic Optimization, AIAA 2000-0667 (2000). Search in Google Scholar

[31] R. Natalini, Convergence to equilibrium for the relaxation approximations of conservation laws, Comm. Pure Appl. Math. 49 (1996), 795–823. Search in Google Scholar

[32] R. Natalini, A discrete kinetic approximation of entropy solutions to multidimensional scalar conservation laws, J. Differential Equations148(1998), 292–317. Search in Google Scholar

[33] R. Natalini and A. Terracina, Convergence of a relaxation approximation to a boundary value problem for conservation laws, Comm. Partial Differential Equations26(2001), 1235–1252. Search in Google Scholar

[34] J.-M. T. Ngnotchouye, M. Herty, M. K. Banda and S. Veelken, Relaxation Approaches to the Optimal Control of the Euler Equations, Comp. and Appl. Math. 30(2011), 399–425. Search in Google Scholar

[35] N. T. Nguyen, Continuous Adjoint-Based Optimization of Hyperbolic Equations with Nonlinear Differential Equation Constraints on Periodic Boundary Conditions, Preprint. Search in Google Scholar

[36] J. Nocedal and S. J. Wright, Numerical optimization, Springer-Verlag, New York, 1999. Search in Google Scholar

[37] L. Pareschi and G. Russo, Implicit-Explicit Runge-Kutta schemes for stiff systems of differential equations, In: Recent Trends in Numerical Analysis (Eds. Brugano, Trigiante), 3, 2000, 269–284. Search in Google Scholar

[38] L. Pareschi and G. Russo, Asymptotically Strong Stability Preserving schemes for hyperbolic systems with sti_ relaxation. In: Hyperbolic problems: Theory, Numerics, Applications: Proceedings of the Ninth International Conference on Hyperbolic Problems held in Caltech, Pasadena, 2003, 241–255. Search in Google Scholar

[39] L., Pareschi and G. Russo, Implicit-Explicit Runge-Kutta schemes and applications to hyperbolic systems with relaxation, J. Sci. Comput. 252005, 129–155. Search in Google Scholar

[40] P. K. Sweby, High resolution schemes using flux limiters for hyperbolic conservation laws, SIAM J. Num. Anal. 211984, 995–1011. Search in Google Scholar

[41] S. Ulbrich, A sensitivity and adjoint calculus for discontinuous solutions of hyperbolic conservation laws with source terms, SIAM J. Control Optim. 412002, 740. Search in Google Scholar

[42] S. Veelken, M. Herty, J.-M. T. Ngnotchouye and M. K. Banda, Optimal Control of the Euler Equations via Relaxation Approaches, PAMM102010, 595–596. Search in Google Scholar

[43] W. C. Wang and Z. Xin, Asymptotic limit of initial boundary-value problems for conservation laws with relaxation extensions, Comm. Pure Appl. Math. 51(1998), 505–535. Search in Google Scholar

Received: 2013-3-26
Accepted: 2014-8-9
Published Online: 2016-3-16
Published in Print: 2016-3-1

© 2016 by Walter de Gruyter Berlin/Boston