Stabilization

EulerEquationMod.edgestabilize!

Adds edge-stabilization (see Burman and Hansbo, doi:10.1016/j.cma.2003.12.032) to a given residual. Different methods are available depending on the rank of u:

• For scalar fields, it is assumed that u is a rank-2 array, with the first

dimension for the local-node index, and the second dimension for the element index.

• For vector fields, u is a rank-3 array, with the first dimension for the

index of the vector field, the second dimension for the local-node index, and the third dimension for the element index.

Naturally, the number of entries in the dimension of u (and res) corresponding to the nodes must be equal to the number of nodes in the SBP operator sbp.

Inputs

• sbp: an SBP operator type

• ifaces: list of element interfaces stored as an array of Interfaces

• u: the array of solution data, numDofPerNode x numNodesPerElement x numEl. for scalar fields (ie. numDofPerNode), this can be a 2D array.

• x: Cartesian coordinates stored in (coord,node,element) format

• dxidx: scaled Jacobian of the mapping (as output from mappingjacobian!)

• jac: determinant of the Jacobian array, numNodesPerElement x numEl

• alpha: array of transformation terms (see below)

• stabscale: numfaces x numNodes per face array of scaling parameter

In/Outs

• res: where the result of the integration is stored, same shape as u

Details

The array alpha is used to compute the directional derivative normal to the faces. For a 2-dimensional problem, it can be computed as follows:

  for k = 1:mesh.numelem
    for i = 1:sbp.numnodes
      for di1 = 1:2
        for di2 = 1:2
          alpha[di1,di2,i,k] = (dxidx[di1,1,i,k].*dxidx[di2,1,i,k] + 
                            dxidx[di1,2,i,k].*dxidx[di2,2,i,k])*jac[i,k]
        end
      end
    end
  end

EulerEquationMod.stabscale

This function calculates the edge stabilization scalaing parameter at a node and returns it.

Inputs: * q : vector of conservative variables * dxidx : jacobian of xi wrt x coordinates at the node * nrm : normal vector in xi space * params : ParamType{2}

This is a low level function

EulerEquationMod.stabscale

This function calculate the stabilization scaling parameter across the entire mesh by calling the low level method. This populates eqn.stabscale.

This is a mid level function

EulerEquationMod.GLS

Add Galerkin Least-Squares stabilization to the weak residual. This implementation is only for steady problems and conservative variables

Inputs

• mesh: AbstractMesh type

• sbp : Summation-by-parts operator

• eqn : Equation object used elsewhere

Outputs

• None

EulerEquationMod.SUPG

Add Streamwise Upwind Petrov-Galerkin stabilization to the weak residual. This implementation is only for steady problems

Inputs

• mesh: AbstractMesh type

• sbp : Summation-by-parts operator

• eqn : Equation object used elsewhere

Outputs

• None

EulerEquationMod.calcAxidxi

Calculates AxiDxi + AetaDeta at the element level

Inputs

• Axidxi : Product of flux jacobian with shape function derivative AxiDxi + AetaDeta

• shapefuncderiv: Shape function derivative (Dxi and Deta above)

• Axi : Flux jacobian in the xi direction

• Aeta : Flux jacobian in the eta direction

• ndof : Number of degrees of freedom per node

• nnpe : Number of nodes per element

Outputs

• None

EulerEquationMod.calcEigenFactorization

Computes the eigen value facotrization of the flux jacobian in either xi or eta direction. This is a nodal level function. It is only for 2D conservative variables.

Reference: $Efficient Solution Methods for the Navier–Stokes Equations$, T.H, Pulliam, NASA Ames Research Center

Inputs

• eqn : Euler equation object

• q : Conservative variables at a node

• dxidx: mapping jacobian vector. It can be either [dξ/dx dξ/dy] or [dη/dx dη/dy] depending on what Axi or Aeta being factorized.

• T : Matrix of eigen vectors.

• Tinv : Inverse of T

• Lambda : Diagonal matrix of eigen values

Outputs

• None

EulerEquationMod.calcElementArea

Calculates the element area of a 2D triangular element

Inputs

• coords : Vertex coordinates of the element

Outputs

• area : Area of the triangular element

EulerEquationMod.calcFluxJacobian

It computes the Euler flux Jacobian at the nodal level in te physical space (X & Y axes). It uses conservative variables.

Inputs

• eqn : Equation object

• q : Conservative variables at the node

• Ax : Flux Jacobian at a node in the X-direction

• Ay : Flux Jacobian at a node in the Y-direction

Outputs

• None

EulerEquationMod.calcStabilizationTerm

Calculates the stabilization term tau for all the nodes in the mesh

Inputs

• mesh : Abstract mesh type

• sbp : Summation-by-parts operator

• eqn : Equation object within EulerEquationMod

• tau : Stabilization term being calculated

Outputs

• None

EulerEquationMod.calcTau

Calculates the stabilization term tau for GLS

Inputs

• mesh : Abstract mesh type

• sbp : Summation-by-parts operator

• eqn : Euler equation object

• tau : Stabilization term

Outputs

• None

EulerEquationMod.circumcircleDiameter

Calculates the circumcircle diameter of a triangular element

Inputs

• coords : Vertex coordinates of the element

Outputs

• dia : Diameter of the circumcircle

EulerEquationMod.parametricFluxJacobian

Claculates the flux jacobian in the parametric space for all the nodes in the mesh. It uses conservative variables

Inputs

• mesh : Abstract mesh object

• sbp : Summation-by-parts operator

• eqn : Equation object

Outputs

• None

EulerEquationMod.applyFilter

This function multiplies a filter matrix by the variables at every node in the mesh. The filter matrix is stored in eqn.params.filter_mat

Inputs: mesh sbp eqn opts

Keyword Args: trans=false : transpose the filter matrix or not.

Inputs/Outputs: arr: 3D array to multiply the filter matrix by

EulerEquationMod.calcfilter

This function calculates the filter used by applyFilter(). The filter is calculated as Vfilter_kernelinv(V), where V is a matrix that converts the interpolating SBP basis to a modal basis and filter_kernal is a matrix (typically diagonal) that performs the filtering of the modal basis coefficients.

Inputs: sbp: SBP operator used to compute the solutions filter_name: and String that matches the name of a filter function. This string is used to retrieve the function from a dictionary of all supported filter function. opts: options dictonary

Outputs: F_ret: a numNodesPerElement x numNodesPerElement filter operator

EulerEquationMod.calcLowPassFilter

This function calculates a low pass filter diagonal filter matrix.

Inputs: sbp: SBP operator opts: options dictonary

Outputs: filt: numNodesPerElement x numNodesPerElement diagonal filter matrix

EulerEquationMod.calcModalTransformationOp

This function calculates a matrix operator V that transforms from a modal basis to an interpolating basis for SBP operators. Because the transformation itself is rank deficient, V is augmented with the last n vector of its Q matrix (from the full QR decomposition), where numNodesPerElement - rank(V) = n

Inputs: sbp: SBP operator

Outputs: V_full: a numNodesPerElement x numNodesPerElement generalized Vandermonde matrix

EulerEquationMod.calcRaisedCosineFilter

This function constructs a diagonal matrix filt for a 2 dimensional basis filter, using the raised cosine filter described in Spectral Methods for the Euler Equations: Part I by Hussaini, Kproiva, Salas, Zang

Inputs sbp: SBP operator opts: options dictonary

Outputs: filt: an numNodesPerElement x numNodesPerElement diagonal filter matrix

EulerEquationMod.getPascalLevel

This function returns the level of Pascals Triangle a particular node lies in.

Inputs: node: node index

Outputs: level: integer describing level of Pascals triangle.

EulerEquationMod.applyDissipation

This function multiplies the dissipation matrix operator for each element by the values in a 3D array (typically eqn.q) for each node of the element and stores the results in eqn.res.

The dissipation matrix operators must be stored in eqn.dissipation_mat, a numNodesPerElement x numNodesPerElement x numEl array.

Inputs: mesh sbp eqn opts arr: a 3D array to apply the dissipation matrix to

Aliasing restrictions: no guarantees what happens if arr is eqn.res

EulerEquationMod.calcDissipationOperator

This function calculates the numNodesPerElement x numNodesPerElement dissipation matrix operator for each element and returns them in an array numNodesPerElement x numNodesPerElement x numEl array.

The dissipation matrix operator is calculated as epsilonfilt.'h_jac*filt, where filt is a (usually diagonal) filter matrix, h_jac is a scaling term that incorporates information about the shape of the element.

Inputs: mesh sbp eqn opts dissipation_name: an String of the function name used to retrieve the function that generates the matrix filt from a dictonary

Outputs: dissipation_mat: a numNodesPerElement x numNodesPerElement x numEl array

Aliasing restrictions: none

EulerEquatoinMod.getDissipatoinFilterOperator

This function gets the dissipation filter operator used to construction the artificial dissipation matrix.

The filter is calculated as Vfilter_kernelinv(V), where V is a matrix that converts the interpolating SBP basis to a modal basis and filter_kernal is a matrix (typically diagonal) that performs the filtering of the modal basis coefficients.

Inputs: sbp: SBP operator filter: function to call to calculate the filter kernal

Outputs: F: a numNodesPerElement x numNodesPerElement filter matrix