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.
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 Diablo Developers Guide. 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 Diablo Developers Guide 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:
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
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:
subroutine ApplySecondDeriv(y,ddy,Dx2,N)!-- NOTE: Same as above, haven't actually used this subroutine in a programinteger, parameter::dp=kind(1.d0)integer, intent(in)::Nreal(kind=dp), intent(in)::y(N)real(kind=dp), intent(out)::ddy(N)type(SBP2), intent(in)::Dx2integer::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 ddydo 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)enddoenddoendsubroutine ApplySecondDeriv
Table of Contents
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.
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 Diablo Developers Guide. 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.
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 Diablo Developers Guide 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:
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 ApplyFirstDerivSecond 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: