Jetstream+High-Order+Methods

Click here to return to top level Jetstream page.toc

=**Description of high-order methods**=

Diablo makes use of high-order diagonal norm summation-by-parts (SBP) operators with weakly imposed boundary conditions through simultaneous approximation terms (SAT) to improve stability properties. Boundary conditions imposed strongly (such as the injection method) could destroy the SBP property. One of the most important properties of an operator which satisfies the SBP property is that an energy estimate can be obtained. Various accuracy criteria can be imposed to design application specific operators.

**First Derivative**
Approximating derivatives with high-order finite difference operators has potential for efficiency improvements. In Diablo, derivatives are approximated using diagonal-norm SBP operators. These derivative operators are constructed so that they guarantee stability for linear problems (a necessary but not sufficient condition for stability in non-linear problems).

**Second Derivative**
The second derivative can be formed naively by applying the first derivative operator twice. This is also known as the wide, or non-compact second derivative operator. From an implementation standpoint this is relatively straightforward, but several properties of the operator can be improved by using the compact second derivative operator.

**Artificial Dissipation**
Artificial dissipation is required to stabilize the non-linear solution algorithm, as well as add diagonal dominance to the linear problem at each solution update. The artificial dissipation operator is constructed such that it is a positive semi-definite matrix operator, and incorporates the norm operator used to construct the first derivative operator.

math \mathcal{\hat{R}}_{AD} = H^{-1}D_b B D_b J \hat{Q} math

=**Using High-Order SBP Operators in Diablo**=

There are several input parameters within Diablo, and the default is to use 2nd-order accurate SBP operators. Diablo makes use of a custom derived type to store and apply SBP operators; this can be found in ./common/SBPSAT_Mod.f90. A brief description of the data type can be found in. SBP operators of up to a global order-of-accuracy of 6 are available, but only orders 2,3,4 have been extensively tested with analysis. I've verified the implementation of the 5th and 6th-order operators with a 1D polynomial test, but it seems as though the linear problems are very difficult to solve in a practical analysis problem (might need a better preconditioner).

**Euler Analysis**
Anything within Euler analysis and optimization is programmed to use the SBP operator stored within **blk% Dx**. The order-of-accuracy of the SBP operator is defined by the parameter p. "p" indicates the largest degree of discretized polynomial to which the SBP operator can be applied to obtain the precise derivative. The relevant input parameter here is **diablo% order**, which equals 2*p and is the interior order-of-accuracy of the operator.


 * ~ **Global Order** ||~ **p** ||~ **2p** ||
 * < 2nd-order ||= 1 ||= 2 ||
 * < 3rd-order ||= 2 ||= 4 ||
 * < 4th-order ||= 3 ||= 6 ||
 * < 5th-order ||= 4 ||= 8 ||
 * < 6th-order ||= 5 ||= 10 ||

For example, if I wanted to operator Diablo at 3rd-order I would choose **diablo% order = 4**. This should be the only relevant setting to worry about for Euler analysis directly related to the SBP operator. **diablo% diss_order** is the global order-of-accuracy of the dissipation operator; so for the 3rd-order example above, **diablo% diss_order = 3**. For 4th-order **diablo% order = 6** and **diablo% diss_order = 4**. I've also made a distinction between different boundary schemes in the dissipation operator, controlled by **diablo% diss_bound**. **diablo% diss_bound = 1** is first order at the boundaries and 2 is 2nd-order (neither of these change the global order-of-accuracy of the operator since they represent a high derivative, **diss_order = 3** with **diss_bound = 1** uses a third derivative which is 1st-order). I recommend **diablo% diss_order = 1** as it seems to be more stable.

__Summary:__


 * diablo% order** = {2,4,6,8,10} !-- interior order-of-accuracy, only 2,4,6 have been extensively tested.
 * diablo% SBP_type** = {1,2} !-- only valid for 2nd-order operator; 1=standard 2=extended boundary
 * diablo% diss_order** = {2,3,4,5,6} !-- global order-of-accuracy, only 2,3,4 have been extensively tested.
 * diablo% diss_bound** = {1,2} !-- boundary scheme for dissipation operator

**Viscous Laminar Analysis**
The 2nd-order code written by Michal Osusky is still present in Diablo, but an alternative formulation to most of the functionality provided by his code makes use of the **blk% Dx** operator. To use the **blk% Dx** operator when constructing the viscous contributions to the residual vector or performing any surface integrations involving frictional forces, make sure to set **diablo% HO = .true.**. The default of .false. will use subroutines within Michal's version of the code, which do not make use of the **blk% Dx** operator and instead use hard-coded values (this might actually be a bit faster, its just less general).

To use the non-compact operator **diablo% vis_compact = .false.**, and vice-versa. The compact operator makes use of a new derived type which I have developed with David DRF and is very similar in usage to the SBP derived type and is stored in **blk% Dx2**; more details to come soon.

There is some debate as to what boundary conditions should be applied, and may still be a work in progress. The only boundary condition that provides any option is the wall SAT. This is controlled by **diablo% visc_srf_sat**. Set this to 5 to use Michal's general formulation for the adiabatic wall boundary condition which uses a flux function penalty rather than a penalty on the conservative variables (this flux function version has been upgraded to use **blk%Dx**).

__Summary:__


 * diablo% HO** = {T/F} !-- set to true to use the **blk%Dx** and **blk%Dx2** operators when forming viscous contributions
 * diablo% vis_cmpct** = {T/F} !-- true will use the compact operator stored in **blk%Dx2** and false will apply **blk%Dx** twice.
 * diablo% visc_srf_sat** = 5 !-- use Michal's general flux penalty to enforce adiabatic wall boundary condition

**Spalart-Allmaras Turbulence Model**
Work in progress...

=Derived Types for SBP Operators=

First Derivative Operator
The first derivative operator is stored using the SBP derived type in Diablo. This can be initialized at the beginning of the subroutine/program by type(SBP) :: Dx for example. Dx would be populated by calling call buildFirstDeriv(Dx,order,N) with Dx initialized as above, order representing the desired order-of-accuracy, and N being the number of nodes in the particular direction of the operator. "call buildFirstDeriv(blk%Dx(di), diablo% order, blk%jkmmax(di))" would populate the Dx operator for direction di for the default settings in Diablo.

Refer to the for details of how the various arrays within the SBP derived type are used. Here is a 1-dimensional example of how to apply the first derivative SBP operator to an array y: code subroutine ApplyFirstDeriv(y,dy,Dx,N)

!-- NOTE: I haven't actually used this specific subroutine but it should work, and should be sufficient !-- to get the point across.

integer, parameter :: dp = kind(1.d0) integer, intent(in) :: N real(kind=dp), intent(in) :: y(N) real(kind=dp), intent(out) :: dy(N) type(SBP), intent(in) :: Dx

integer :: j,jp,p real(kind=dp) :: as

dy = 0.d0             !-- initialize derivative to zero

!-- apply SBP operator to all entries of y to form dy do j = 1,N !-- ibeg gives first index for row j, iend gives last index for row j do p = Dx% ibeg(j), Dx% iend(j) jp = j + Dx%ja(p) !-- row offset as = Dx% as(p)    !-- SBP operator coefficient dy(j) = dy(j) + as*y(jp) end do end do

end subroutine ApplyFirstDeriv

code

Second Derivative Operator
The second derivative derived type, SBP2, is closely related to the derived type SBP. More details to come soon on specific array meanings.

Below is a 1-dimensional example of its implementation:

code format="fortran" subroutine ApplySecondDeriv(y,ddy,Dx2,N)

!-- NOTE: Same as above, haven't actually used this subroutine in a program

integer, parameter :: dp = kind(1.d0) integer, intent(in) :: N real(kind=dp), intent(in) :: y(N) real(kind=dp), intent(out) :: ddy(N) type(SBP2), intent(in) :: Dx2

integer :: j,jp,p,jpp real(kind=dp) :: as

ddy = 0.d0   !-- initialize second derivative to zero

!-- apply SBP2 operator to all entries of y to form ddy do j = 1,N do p = Dx2% ibeg(j), Dx2% iend(j) jp = j + Dx2% ja(p) jpp = j + Dx2% ia(p) ddy(jpp) = ddy(jpp) + Dx2% as(p)*y(jp) end do end do

end subroutine ApplySecondDeriv

code