Reference

SummationByParts.bndrynodalexpansion

Computes the transformation matrix that maps the Proriol orthogonal polynomials to polynomials that are nodal on the boundary nodes, i.e. if E is the transformation matrix and P is the matrix of Proriol polys, with the polynomials listed by row, then P*E = I, when restricted to the boundary nodes.

Inputs

• cub: symmetric cubature rule
• vtx: vertices of the right simplex
• d: maximum total degree for the Proriol polynomials

Outputs

• E: transformation matrix

SummationByParts.nodalexpansion

Computes the transformation matrix that maps the Proriol orthogonal polynomials to polynomials that are nodal on the cubature nodes, i.e. if C is the transformation matrix and P is the matrix of Proriol polys, with the polynomials listed by row, then P*C = I, when restricted to the nodes. This function can be used to construct the operators in the spectral-element method of Giraldo and Tayler.

Inputs

• cub: symmetric cubature rule
• vtx: vertices of the right simplex
• d: maximum total degree on the edges
• e: maximum degree of bubble functions on the interior

Outputs

• C: transformation matrix

SummationByParts.boundaryoperators

Finds the symmetric part of the SBP operators, Ex, Ey (Ez). These operators coorespond to boundary integrals in the divergence theorem, and are related to the mass matrices of the boundary faces.

Inputs

• cub: symmetric cubature rule
• vtx: vertices of the right simplex
• d: maximum total degree for the Proriol polynomials

Outputs

• Ex, Ey (Ez): symmetric parts of the SBP first derivative operators

SummationByParts.boundarymassmatrix

Returns the (dense) mass matrix for a set of nodes on a reference boundary. This mass matrix can be used for boundary integration or to impose boundary conditions weakly. The array bndryindices is also returned, which is a list of element-node indices for each boundary (see also SymCubatures.getfacenodeindices).

Inputs

• cub: symmetric cubature rule
• vtx: vertices of the right simplex
• d: maximum total degree for the polynomials

Outputs

• Hbndry: reference boundary mass matrix
• bndryindices: list of nodes that lie on each boundary

SummationByParts.accuracyconstraints

Returns the accuracy constraints on the asymmetric part of the SBP stiffness matrices. These constraints are linear, and for each coordinate-direction operator (i.e. Qx, Qy,...) the system matrix A is the same; only the right-hand side changes.

The columns in A are ordered assuming only the strictly lower triangular part of the operators are the unknowns. These unknowns are ordered by row and then column. For example, entry Q21 = -Q12 is the number 1 variable, and Q32 = -Q23 is the number 3 variable.

Inputs

• cub: symmetric cubature rule
• vtx: vertices of the right simplex
• d: maximum total degree for the Proriol polynomials
• E: the symmetric part of the SBP stiffness matrices
• dl: (optional) the minimum total degree for the Proriol polynomials

Outputs

• A: the system matrix for the linear accuracy constraints
• bx,by (bz): the right-hand-sides of the accuracy constraints

SummationByParts.commuteerror

Returns the commute-error objective value. For 2D SBP operators, this is defined as ||H(DxDy - DyDx)||^2, where the norm is the Frobenius norm. For 3D operators, the error is defined as ||H(DxDy - DyDx)||^2 + ||H(DxDz - DzDx)||^2 + ||H(DyDz - DzDx||^2.

Inputs

• w: cubature rule weights
• Qxpart,Qypart (Qzpart): Q operators that satisfy the accuracy conditions
• Z: basis for the null space of the accuracy constraints (may be empty)
• reducedsol: the weights for Z; the first [1:numnodes] elements are for Qx

Outputs

• f: commute-error objective value

SummationByParts.buildoperators

Construct and return the SBP matrix operators, specifically the diagonal norm matrix and the stiffness matrices.

Inputs

• cub: symmetric cubature rule
• vtx: vertices of the right simplex
• d: maximum total degree for the Proriol polynomials

Outputs

• w: the diagonal norm stored as a 1D array
• Q: the stiffness matrices

Note: when a fifth input is included, it is interpreted as the degree of the interior bubble functions, and the Q operators returned are the spectral-element differentiation operators of Giraldo and Taylor

SummationByParts.buildsparseoperators

Construct and return sparse SBP matrix operators, specifically the diagonal norm matrix and the stiffness matrices. These are sparse in the sense that there are zeros in the S matrices, but they are not returned as sparse matrices.

Inputs

• cub: symmetric cubature rule
• vtx: vertices of the right simplex
• d: maximum total degree for the Proriol polynomials

Outputs

• w: the diagonal norm stored as a 1D array
• Q: the stiffness matrices

SummationByParts.buildMinConditionOperators

Construct and return SBP matrix operators that minimize the condition number of a model advection problem. These operators have diagonal mass and boundary operators.

Inputs

• cub: symmetric cubature rule
• vtx: vertices of the right simplex
• d: maximum total degree for the Proriol polynomials
• tol: (optional) gradient norm tolerance for the optimization problem
• vertices: (optional) if true, include vertices in the operator
• opthist: (optional) if true, show the optimization history

Outputs

• w: the diagonal norm stored as a 1D array
• Q: the stiffness matrices

SummationByParts.getnodepermutation

The node ordering produced by SymCubature is not convenient for mapping local to global node indices in the global residual assembly procedure. This function returns a reordering that is more suited for local-to-global mapping.

Note: the edge parameters of cub are assumed to be less than T(0.5).

WARNING: the faceperm array has not been thoroughly tested

Inputs

• cub: symmetric cubature rule
• d: maximum polynomial degree for which the SBP operator is exact

Outputs

• perm: a permutation vector of indices
• faceperm: a permutation vector for the face indices

SummationByParts.buildoperators_pocs

Construct SBP operators using the Projection Onto Convex Sets (POCS) algorithm

Inputs

• cub: symmetric cubature rule
• vtx: vertices of the right simplex
• d: maximum total degree for the Proriol polynomials
• vertices: (optional) if true, include vertices in the operator

Outputs

• D: the derivative matrices
• H: the norm/mass matrix
• E: the boundary integration matrices
• Q: the stiffness matrices
• S: the skew-symmetric matrices

SummationByParts.buildfacereconstruction

Builds a matrix operator that can reconstruct a field from a set of volume nodes to a set of face nodes. The reconstruction operator is only constructed for one face, but a permutation array is returned that allows the same operator to be applied on all the faces.

Inputs

• facecub: symmetric cubature rule for the face
• cub: symmetric cubature rule for the volume
• vtx: vertices of the right simplex
• d: maximum total degree for the Proriol polynomials

Returns

• R: the volume- to face-node reconstruction operator
• perm: a permutation array that allows R to be applied on all faces

SummationByParts.buildfacederivatives

Builds matrix operators that can differentiate a polynomial field of degree d from a set of volume nodes to a set of face nodes. The derivative operators are only constructed for one face, but a permutation array is returned that allows the same operators to be applied on all the faces.

Note: the derivative operators are for the tangential and normal directions, and do not necessarily correspond to the directions ξ and η (and ζ) used for the volume derivatives. These face derivatives are intended for edge stabilization.

Inputs

• facecub: symmetric cubature rule for the face
• cub: symmetric cubature rule for the volume
• vtx: vertices of the right simplex
• d: maximum total degree for the Proriol polynomials

Returns

• D: derivative operators in [face node, vol node, direction] format
• perm: a permutation array that allows D to be applied on all faces

Cubature.quadrature{T}

This high-level function computes and returns a symmetric cubature of requested accuracy on the interval [-1,1]

Inputs

• q: maximum degree of polynomial for which the cubature is exact
• T: the data type used to represent the cubature
• internal: if true, all nodes are strictly internal (default false)
• tol: tolerance with which to solve the cubature

Outputs

• cub: a symmetric cubature for the interval [-1,1]
• vtx: vertices, [-1,1]

Cubature.cubatureresidual

This method computes the residuals, F, between a cubature, defined by cub, and the true value of an integral. Each residual corresponds with an orthogonal polynomial on the simplex up to degree q. The Jacobian, dF, of the residual, with respect to the quadrature nodes and weights, is also returned.

Inputs

• cub: defines the nodes and weights of the cubature via symmetry orbits
• q: maximum degree of the othogonal polynomials used in the conditions
• calc_grad: indicates whether to calculate the gradients the orthogonal polynomials

Outputs

• F: the accuracy conditions for orthogonal polynomials up to degree q
• dF: derivative of F with respect to x, y, (z,) w, in that order

Cubature.equivalenceconstant{T}

Computes the equivalence constant for a given cubature; that is, it finds the maximum eigenvalue for the matrix pk^T H pm, where H = diag(weights) and pk denotes the orthogonal polynomial evaluated at the cubature points.

Inputs

• cub: symmetric cubature rule
• vtx: vertices of the right simplex
• q: maximum degree of polynomial for which the cubature is to be tested

Outputs

• λmax: maximum eigenvalue, which is the equivalence constant

Cubature.getTetCubatureDiagE{T}

Returns a cubature rule and vertices for the SBP DiagE operators on tetrahedra; these are cubatures that have nodes on the boundary that are analogous to LG or LGL quadrature rules, which then leads to diagonal E.

Inputs

• q: maximum degree of polynomial for which the cubature is exact
• T: the data type used to represent the cubature
• faceopertype: the operator type on the facets of the tetrahedron
• tol: tolerance with which to solve the cubature

Outputs

• cub: a symmetric cubature for the right tetrahedron
• vtx: vertices for the right tetrahedron

Cubature.getTetCubatureGamma{T}

Returns a cubature rule and vertices for the SBP Gamma operators on tetrahedra; these are operators with (p+1)(p+2)/2 nodes on each face, where, typically, p = (q+1)/2.

Inputs

• q: maximum degree of polynomial for which the cubature is exact
• T: the data type used to represent the cubature
• tol: tolerance with which to solve the cubature

Outputs

• cub: a symmetric cubature for the right tetrahedron
• vtx: vertices for the right tetrahedron

Cubature.getTetCubatureOmega{T}

Returns a cubature rule and vertices for the SBP Omega operators on tetrahedra; these are cubatures that are analogous to Gauss-Legendre in 1D, and they are strictly internal to the tet.

Inputs

• q: maximum degree of polynomial for which the cubature is exact
• T: the data type used to represent the cubature
• tol: tolerance with which to solve the cubature

Outputs

• cub: a symmetric cubature for the right tetrahedron
• vtx: vertices for the right tetrahedron

Cubature.getTriCubatureDiagE{T}

Returns a cubature rule and vertices for the SBP DiagE operators on triangles; these are cubatures that have nodes on the boundary that correspond with LG or LGL quadrature rules, which then leads to diagonal E.

Inputs

• q: maximum degree of polynomial for which the cubature is exact
• T: the data type used to represent the cubature
• vertices: if true then vertices are included
• tol: tolerance with which to solve the cubature

Outputs

• cub: a symmetric cubature for the right triangle
• vtx: vertices for the right triangle

Cubature.getTriCubatureDiagELineLGL{T}

Returns a cubature rule and vertices for the Line-LGL SBP DiagE operators on triangles; these are cubatures that have nodes on the boundary that correspond with LGL quadrature rules, which then leads to diagonal E.

Inputs

• q: maximum degree of polynomial for which the cubature is exact
• T: the data type used to represent the cubature
• tol: tolerance with which to solve the cubature

Outputs

• cub: a symmetric cubature for the right triangle
• vtx: vertices for the right triangle

Cubature.getTriCubatureDiagE{T}

Returns a cubature rule and vertices for facets of the SBP DiagE operators on tetrahedra; these should not be used for 2D problems as they do not satisfy the accruacy requirements along their edge.

Inputs

• q: maximum degree of polynomial for which the cubature is exact
• T: the data type used to represent the cubature
• faceopertype: the operator type on the facets of the tetrahedron
• tol: tolerance with which to solve the cubature

Outputs

• cub: a symmetric cubature for the right triangle
• vtx: vertices for the right triangle

Cubature.getTriCubatureGamma{T}

Returns a cubature rule and vertices for the SBP Gamma operators on triangles; these are operators with p+1 nodes on each face, where, typically, p = (q+1)/2.

Inputs

• q: maximum degree of polynomial for which the cubature is exact
• T: the data type used to represent the cubature
• tol: tolerance with which to solve the cubature

Outputs

• cub: a symmetric cubature for the right triangle
• vtx: vertices for the right triangle

Cubature.getTriCubatureOmega{T}

Returns a cubature rule and vertices for the SBP Omega operators on triangles; these are cubatures that are analogous to Gauss-Legendre in 1D, and they are strictly internal to the triangle.

Inputs

• q: maximum degree of polynomial for which the cubature is exact
• T: the data type used to represent the cubature
• tol: tolerance with which to solve the cubature

Outputs

• cub: a symmetric cubature for the right triangle
• vtx: vertices for the right triangle

Cubature.getTriCubatureOmegaLG{T}

Returns a cubature rule and vertices for the SBP Omega operators on triangles; these are cubatures derived by eliminating nodes in the line-LG cubature.

Inputs

• q: maximum degree of polynomial for which the cubature is exact
• T: the data type used to represent the cubature
• tol: tolerance with which to solve the cubature

Outputs

• cub: a symmetric cubature for the right triangle
• vtx: vertices for the right triangle

Cubature.getTriCubatureOmegaLineLG{T}

Returns a cubature rule and vertices for the SBP Omega operators on triangles; these are cubatures derived by using tensor-product of half of the nodes in the LG quadrature rules.

Inputs

• q: maximum degree of polynomial for which the cubature is exact
• T: the data type used to represent the cubature
• tol: tolerance with which to solve the cubature

Outputs

• cub: a symmetric cubature for the right triangle
• vtx: vertices for the right triangle

Cubature.pointCubature

This returns a (trivial) point cubature and default vertex -1

Inputs

• T: the data type used to represent the cubature

Outputs

• cub: a symmetric cubature for point
• vtx: vertex, [-1]

Cubature.quadratureUniform{T}

This high-level function computes and returns a uniform cubature of requested accuracy on the interval [-1,1]

Inputs

• q: maximum degree of polynomial for which the cubature is exact
• N: number of nodes (N >= q+1)
• T: the data type used to represent the cubature

Outputs

• cub: a symmetric cubature for the interval [-1,1]
• vtx: vertices, [-1,1]

Cubature.solvecubature!{SymCub{T}}

Attempts to solve for the nodes and weights of a cubature that is exact for polynomials of degree r <= q. The nodes and weights of the cubature are defined by cub, which is a parametric abstract type (see symcubatures.jl).

Inputs

• cub: symmetric cubature rule
• q: maximum (desired) degree for which the cubature is exact
• mask: array of indicies of parameters and weights that are free
• tol: tolerance with which to solve the accuracy conditions
• hist: if true, print the residual-norm convergence history
• xinit: initial parameter guess
• delta1: purturbation constant that determines how much the parameters are perturbed in the case of stagnation with fmin <= 0.1
• delta2: purturbation constant that determines how much the parameters are perturbed in the case of stagnation with fmin > 0.1

In/Outs

• cub: on entry, defines the initial guess for the cubature nodes and weights. on exit, defines the nodes and weights that satisfy the desired accuracy.

Cubature.solvecubatureweights!{SymCub{T}}

Attempts to solve for the weights of a cubature that is exact for polynomials of degree r <= q. The weights (and nodes) of the cubature are defined by cub, which is a parametric abstract type (see symcubatures.jl).

Inputs

• q: maximum (desired) degree for which the cubature is exact
• tol: tolerance with which to solve the accuracy conditions
• hist: if true, print the residual-norm convergence history

In/Outs

• cub: on entry, defines the initial guess for the cubature nodes and weights. on exit, defines the nodes and weights that satisfy the desired accuracy.

Cubature.tetcubature{T}

This high-level function computes and returns a symmetric cubature of requested accuracy on the right tetrahedron.

Inputs

• q: maximum degree of polynomial for which the cubature is exact
• T: the data type used to represent the cubature
• internal: if true, all nodes are strictly internal (default false)
• facequad: if true, the cubatures' face nodes coincide with a quadrature
• tol: tolerance with which to solve the cubature

Outputs

• cub: a symmetric cubature for the right tetrahedron
• vtx: vertices for the right tetrahedron

Cubature.tricubature{T}

Deprecated; this function will be removed in the future

Inputs

• q: maximum degree of polynomial for which the cubature is exact
• T: the data type used to represent the cubature

Outputs

• cub: a symmetric cubature for the right triangle
• vtx: vertices for the right triangle

SummationByParts.deriveTriCubatureOmega

This function derives quadrature rules for SBP-Omega operators on the triangle.

Inputs

• q: the degree of the operator
• vertices : if true, vertices are present in the set of nodes
• midedges : if true, edge midpoints are present in set of nodes
• centroid : if true, centroid is present in set of nodes
• numedge : number of unique edge parameters
• numS21 : number of S21 orbits (vertex to opposite face)
• numS111 : number of S111 orbits
• delta1: purturbation constant that determines how much the parameters are perturbed in the case of stagnation with fmin <= 0.1
• delta2: purturbation constant that determines how much the parameters are perturbed in the case of stagnation with fmin > 0.1
• verbose: print out iteration results
• xinit: initial guess of the prameter and/or the parameter and weights
• xinit_sym_group: list of the symmetry group ordering provided in xinit

Outputs

• cub: a symmetric cubature for the right triangle
• vtx: vertices for the right triangle

SummationByParts.deriveTriCubatureGamma

This function derives quadrature rules for SBP-Gamma operators on the triangle.

Inputs

• q: the degree of the operator
• vertices : if true, vertices are present in the set of nodes
• midedges : if true, edge midpoints are present in set of nodes
• centroid : if true, centroid is present in set of nodes
• numedge : number of unique edge parameters
• numS21 : number of S21 orbits (vertex to opposite face)
• numS111 : number of S111 orbits
• delta1: purturbation constant that determines how much the parameters are perturbed in the case of stagnation with fmin <= 0.1
• delta2: purturbation constant that determines how much the parameters are perturbed in the case of stagnation with fmin > 0.1
• verbose: print out iteration results
• xinit: initial guess of the prameter and/or the parameter and weights
• xinit_sym_group: list of the symmetry group ordering provided in xinit

Outputs

• cub: a symmetric cubature for the right triangle
• vtx: vertices for the right triangle

SummationByParts.deriveTriCubatureDiagE

This function derives quadrature rules for SBP-Omega operators on the triangle.

Inputs

• q: the degree of the operator
• vertices : if true, vertices are present in the set of nodes
• midedges : if true, edge midpoints are present in set of nodes
• centroid : if true, centroid is present in set of nodes
• numedge : number of unique edge parameters
• numS21 : number of S21 orbits (vertex to opposite face)
• numS111 : number of S111 orbits
• delta1: purturbation constant that determines how much the parameters are perturbed in the case of stagnation with fmin <= 0.1
• delta2: purturbation constant that determines how much the parameters are perturbed in the case of stagnation with fmin > 0.1
• verbose: print out iteration results
• xinit: initial guess of the prameter and/or the parameter and weights
• xinit_sym_group: list of the symmetry group ordering provided in xinit
• xedge: parameters for the edge symmetry groups
• xedge_sym_group: list of the symmetry group ordering provided in xedge

Outputs

• cub: a symmetric cubature for the right triangle
• vtx: vertices for the right triangle

SummationByParts.deriveTetCubatureOmega

This function derives quadrature rules for SBP-Omega operators on the tetrahedron.

Inputs

• q: the degree of the operator
• vertices : if true, vertices are present in the set of nodes
• midedges : if true, edge midpoints are present in set of nodes
• centroid : if true, centroid is present in set of nodes
• facecentroid : if true, face centroids are present in the set of nodes
• numedge : number of unique edge parameters
• numfaceS21 : number of S21 face orbits (same tri orbit on face)
• numfaceS111 : number of S111 face orbits (same tri orbit on face)
• numS31 : number of S31 orbits (vertex to opposite face)
• numS22 : number of S22 orbits
• numS211: number of S211 orbits
• numS1111: number of S1111 orbits
• delta1: purturbation constant that determines how much the parameters are perturbed in the case of stagnation with fmin <= 0.1
• delta2: purturbation constant that determines how much the parameters are perturbed in the case of stagnation with fmin > 0.1
• verbose: print out iteration results
• xinit: initial guess of the prameter and/or the parameter and weights
• xinit_sym_group: list of the symmetry group ordering provided in xinit

Outputs

• cub: a symmetric cubature for the right triangle
• vtx: vertices for the right triangle

SummationByParts.deriveTetCubatureGamma

This function derives quadrature rules for SBP-Gamma operators on the tetrahedron.

Inputs

• q: the degree of the operator
• vertices : if true, vertices are present in the set of nodes
• midedges : if true, edge midpoints are present in set of nodes
• centroid : if true, centroid is present in set of nodes
• facecentroid : if true, face centroids are present in the set of nodes
• numedge : number of unique edge parameters
• numfaceS21 : number of S21 face orbits (same tri orbit on face)
• numfaceS111 : number of S111 face orbits (same tri orbit on face)
• numS31 : number of S31 orbits (vertex to opposite face)
• numS22 : number of S22 orbits
• numS211: number of S211 orbits
• numS1111: number of S1111 orbits
• delta1: purturbation constant that determines how much the parameters are perturbed in the case of stagnation with fmin <= 0.1
• delta2: purturbation constant that determines how much the parameters are perturbed in the case of stagnation with fmin > 0.1
• verbose: print out iteration results
• xinit: initial guess of the prameter and/or the parameter and weights
• xinit_sym_group: list of the symmetry group ordering provided in xinit

Outputs

• cub: a symmetric cubature for the right triangle
• vtx: vertices for the right triangle

SummationByParts.deriveTetCubatureDiagE

This function derives quadrature rules for SBP diagonal-E operators on the tetrahedron.

Inputs

• q: the degree of the operator
• vertices : if true, vertices are present in the set of nodes
• midedges : if true, edge midpoints are present in set of nodes
• centroid : if true, centroid is present in set of nodes
• facecentroid : if true, face centroids are present in the set of nodes
• numedge : number of unique edge parameters
• numfaceS21 : number of S21 face orbits (same tri orbit on face)
• numfaceS111 : number of S111 face orbits (same tri orbit on face)
• numS31 : number of S31 orbits (vertex to opposite face)
• numS22 : number of S22 orbits
• numS211: number of S211 orbits
• numS1111: number of S1111 orbits
• delta1: purturbation constant that determines how much the parameters are perturbed in the case of stagnation with fmin <= 0.1
• delta2: purturbation constant that determines how much the parameters are perturbed in the case of stagnation with fmin > 0.1
• verbose: print out iteration results
• xinit: initial guess of the prameter and/or the parameter and weights
• xinit_sym_group: list of the symmetry group ordering provided in xinit
• xedge: parameters for the edge symmetry groups
• xedge_sym_group: list of the symmetry group ordering provided in xedge

Outputs

• cub: a symmetric cubature for the right triangle
• vtx: vertices for the right triangle

SummationByParts.differentiate!

Applies the SBP differentiation matrix operator, D, to data in flux and stores the result in res. Different methods are available depending on the rank of flux:

• For scalar fields, it is assumed that flux 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, flux 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 flux (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
• di: direction index of the operator that is desired (di=1 for Dx, etc)
• flux: the array that the operator is applied to
• ±: PlusFunctor to add to res, MinusFunctor to subract

In/Outs

• res: where the result of applying inv(H)*Q[:,:,di] to u is stored

SummationByParts.differentiateElement!

This is the single-element variant of differentiate! Applies the SBP differentiation matrix operator, D, to data in flux and stores the result in res. Different methods are available depending on the rank of flux:

• For scalar fields, it is assumed that flux is a rank-1 array, with first

and only dimension for the local-node index.

• For vector fields, flux is a rank-2 array, with the first dimension for

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

Naturally, the number of entries in the dimension of flux (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
• di: direction index of the operator that is desired (di=1 for Dx, etc)
• flux: the array that the operator is applied to
• ±: PlusFunctor to add to res, MinusFunctor to subract
• trans (optional): if true, the transpose operation is applied

In/Outs

• res: where the result of applying inv(H)*Q[:,:,di] to u is stored

SummationByParts.differentiate_rev!

This is the reverse differentiated version of differentiate!. See differentiate.jl for further details of the primal method. This function is differentiated with respect to the primal version's flux variable.

Inputs

• sbp: an SBP operator type
• di: direction index of the operator that is desired (di=1 for Dx, etc)
• res_bar: vector applied to the left of the D operator
• ±: PlusFunctor to add to res_bar, MinusFunctor to subract

In/Outs

• flux_bar: the result of the vector matrix product between D and res_bar

SummationByParts.differentiateElement_rev!

This is the reverse differentiated version of differentiateElement!. See differentiate.jl for further details of the primal method. This function is differentiated with respect to the primal version's flux variable.

Inputs

• sbp: an SBP operator type
• di: direction index of the operator that is desired (di=1 for Dx, etc)
• res_bar: vector applied to the left of the D operator
• ±: PlusFunctor to add to res_bar, MinusFunctor to subract

In/Outs

• flux_bar: the result of the vector matrix product between D and res_bar

SummationByParts.directionalDifferentiateElement!

Performs a directional derivative (in reference space) at a given node. The input field u is for a single element, not a collection of elements.

WARNING: In the case of a vector field u, the directional derivative is added to the output Ddir; the user must zero this before.

Inputs

• sbp: an SBP operator type
• dir: a direction vector for the directional derivative
• u: the field that is being differentiated (either a scalar or vector)
• i: index of the node at which the derivative is desired

Returns or In/Outs

• Ddir: derivative of u in direction dir

SummationByParts.weakdifferentiate!

Applies the SBP stiffness matrix (or its transpose) to data in flux and adds to or subtracts from res. Different methods are available depending on the rank of flux:

• For scalar fields, it is assumed that flux 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, flux 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 flux (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
• di: direction index of the operator that is desired (di=1 for Qx, etc)
• flux: the array that the operator is applied to
• ± : PlusFunctor to add to res, MinusFunctor to subract
• trans (optional): if true, the transpose operation is applied

In/Outs

• res: where the result of applying Q[:,:,di] to u is stored

SummationByParts.weakDifferentiateElement!

This is the single-element variant of weakdifferentiate!. Applies the SBP stiffness matrix (or its transpose) to data in flux and adds to or subtracts from res. Different methods are available depending on the rank of flux:

• For scalar fields, it is assumed that flux is a rank-1 array, with the only

dimension for the local-node index.

• For vector fields, flux is a rank-2 array, with the first dimension for

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

Naturally, the number of entries in the dimension of flux (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
• di: direction index of the operator that is desired (di=1 for Qx, etc)
• flux: the array that the operator is applied to
• ± : PlusFunctor to add to res, MinusFunctor to subract
• trans (optional): if true, the transpose operation is applied

In/Outs

• res: where the result of applying Q[:,:,di] to u is stored

SummationByParts.weakdifferentiate_rev!

This is the reverse differentiated version of weakdifferentiate!. See weakdifferentiate.jl for further details of the primal method. This function is differentiated with respect to the primal version's flux variable.

Inputs

• sbp: an SBP operator type
• di: direction index of the operator that is desired (di=1 for Qx, etc)
• res_bar: vector applied to the left of the Q operator
• ± : PlusFunctor to add to res_bar, MinusFunctor to subract
• trans (optional): if true, the transpose operation is applied

In/Outs

• flux_bar: the result of the vector matrix product between Q and res_bar

SummationByParts.weakDifferentiateElement_rev!

This is the reverse differentiated version of weakDifferentiateElement!. See weakdifferentiate.jl for further details of the primal method. This function is differentiated with respect to the primal version's flux variable.

Inputs

• sbp: an SBP operator type
• di: direction index of the operator that is desired (di=1 for Qx, etc)
• res_bar: vector applied to the left of the Q operator
• ± : PlusFunctor to add to res_bar, MinusFunctor to subract
• trans (optional): if true, the transpose operation is applied

In/Outs

• flux_bar: the result of the vector matrix product between Q and res_bar

SummationByParts.weakDifferentiateElement_jac!

This function applies the SBP stiffness matrix (or its transpose) to a diagonal, or block diagonal, matrix of flux jacobians dfluxdu. The resulting matrix is added to or subtracted from dresdu. Different methods are available depending on the rank of dfluxdu:

• For scalar fields, it is assumed that dfluxdu is a rank-1 array, with the

only dimension for the local-node index. dresdu is a rank-2 array, with size sbp.numnodes x sbp.numnodes.

• For vector fields, dfluxdu is a rank-3 array, with the first and second

dimensions for indices of the vector field (at a particular node) and the second dimension for the local-node index. dresdu is a rank-4 array; the first 2 dimensions are of size nvar = size(dfluxdu,1), while the third and fourth dimensions are of size sbp.numnodes.

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

Inputs

• sbp: an SBP operator type
• di: direction index of the operator that is desired (di=1 for Qx, etc)
• dfluxdu: array of the derivative of the flux w.r.t. the state
• ± : PlusFunctor to add to res, MinusFunctor to subract
• trans (optional): if true, the transpose operation is applied

In/Outs

• dresdu: stores Q[:,:,di]diag(dfluxdu) or Q[:,:,di]'diag(dfluxdu)

SummationByParts.edgestabilize!

Applies edge stabilization to a given field, differentiating in the direction specified by dirvec, and scaling by the tau field.

Inputs

• sbpface: an SBP face operator type
• ifaces: list of element interfaces stored as an array of Interfaces
• dirvec: direction to differentiate in [xi coord, face node, L/R, face] format
• tau: scaling term in [face node, face] format
• u: field being stablized in [vol node, element] format
• ±: PlusFunctor to add to res, MinusFunctor to subract

In/Outs

• res: where the result is stored in [vol node, element] format

SummationByParts.integratefunctional!

Integrates a given scalar (or vector) field over the boundary faces.

• For scalar fields, the dimensions of uface correspond to [face-node index, boundary index] and the scalar functional is a return value
• For vector fields, the dimensions of uface correspond to [field index, face-node index, boundary index] and the dimensions of fun correspond to [field index].

Inputs

• sbpface: an SBP face operator type
• bndryfaces: list of boundary faces stored as an array of Boundarys
• flux: array of field data that is being integrated
• ±: PlusFunctor to add to fun, MinusFunctor to subract

In/Outs

• fun: functional value (or vector) being contributed to by the integration

Returns

• fun: in the case of the scalar version, the functional value is returned

SummationByParts.integrateBoundaryFunctional!

This is the single-face variant of integratefunctional!. Integrates a given scalar (or vector) field over the boundary faces.

• For scalar fields, the dimensions of uface correspond to [face-node index] and the scalar functional is a return value

• For vector fields, the dimensions of uface correspond to [field index, face-node index] and the dimensions of fun correspond to [field index].

Inputs

• sbpface: an SBP face operator type
• face: the face of the element to integrate
• flux: array of field data that is being integrated
• ±: PlusFunctor to add to fun, MinusFunctor to subract

In/Outs

• fun: functional value (or vector) being contributed to by the integration

Returns

• fun: in the case of the scalar version, the functional value is returned

SummationByParts.boundaryintegrate!

Scales flux values at boundary cubature points by cubature weights, and then performs transposed interpolation/extrapolation back to volume nodes. Different methods are available depending on the rank of flux:

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

first dimension for the face-node index, and the second dimension for the boundary index.

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

the index of the vector field, the second dimension for the face-node index, and the third dimension for the boundary index.

The dimensions of res are still based on elements; the last dimension is for the element index and the second-last dimension is for the element-local node index.

Inputs

• sbpface: an SBP AbstractFace type
• bndryfaces: list of boundary faces stored as an array of Boundarys
• flux: array of flux data that is being integrated
• ±: PlusFunctor to add to res, MinusFunctor to subract

In/Outs

• res: where the result of the integration is stored

WARNING: the order of the boundaries in bndryfaces and flux must be consistent.

SummationByParts.boundaryFaceIntegrate!

This is the single-face variant of boundaryintegrate!. Scales flux values at boundary cubature points by cubature weights, and then performs transposed interpolation/extrapolation back to volume nodes. Different methods are available depending on the rank of flux:

• For scalar fields, it is assumed that flux is a rank-1 array, with the

first and only dimension for the face-node index.

• For vector fields, flux is a rank-2 array, with the first dimension for

the index of the vector field, and the second dimension for the face-node index.

The dimensions of res are still based on elements; the last dimension (in the scalar case, the only dimension) is for the element-local node index.

Inputs

• sbpface: an SBP AbstractFace type
• face: the face of the element to integrate and project back to the element
• flux: array of flux data that is being integrated
• ±: PlusFunctor to add to res, MinusFunctor to subract

In/Outs

• res: where the result of the integration is stored

SummationByParts.interiorfaceintegrate!

Scales flux values at element-interface cubature points by cubature weights, and then performs transposed interpolation/extrapolation back to volume nodes. Different methods are available depending on the rank of flux:

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

first dimension for the face-node index, and the second dimension for the interface index.

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

the index of the vector field, the second dimension for the face-node index, and the third dimension for the interface index.

The dimensions of res are still based on elements; the last dimension is for the element index and the second-last dimension is for the element-local node index.

Inputs

• sbpface: an SBP AbstractFace type
• ifaces: list of element interfaces stored as an array of Interfaces
• flux: array of flux data that is being integrated
• ±: PlusFunctor to add to res, MinusFunctor to subract

In/Outs

• res: where the result of the integration is stored

WARNING: the order of the interfaces in ifaces and flux must be consistent.

SummationByParts.interiorFaceIntegrate!

This is the single-face variant of interiorfaceintegrate!. Scales flux values at element-interface cubature points by cubature weights, and then performs transposed interpolation/extrapolation back to volume nodes. Different methods are available depending on the rank of flux:

• For scalar fields, it is assumed that flux is a rank-1 array, with the

first and only dimension for the face-node index

• For vector fields, flux is a rank-2 array, with the first dimension for

the index of the vector field, and the second dimension for the face-node index

The dimensions of resL and resR are still based on elements; the last dimension (in the scalar case, the only dimension) is for the element-local node index.

Inputs

• sbpface: an SBP AbstractFace type
• ifaces: list of element interfaces stored as an array of Interfaces
• flux: array of flux data that is being integrated
• ±: PlusFunctor to add to res, MinusFunctor to subract

In/Outs

• resL: where the result of the integration is stored for the left element
• resR: where the result of the integration is stored for the right element

SummationByParts.integratefunctional_rev!

This is the reverse differentiated version of integratefunctional!. See faceintegrate.jl for further details of the primal method. This function is differentiated with respect to the primal version's flux variable.

Inputs

• sbpface: an SBP face operator type
• bndryfaces: list of boundary faces stored as an array of Boundarys
• fun_bar: incoming quantity that multiplies the functional from the left
• ±: PlusFunctor to add to flux_bar, MinusFunctor to subract

In/Outs

• flux_bar: result of the vector-Jacobian product.

SummationByParts.integrateBoundaryFunctional_rev!

This is the reverse differentiated version of integrateBoundaryFunctional!. See faceintegrate.jl for further details of the primal method. This function is differentiated with respect to the primal version's flux variable.

Inputs

• sbpface: an SBP face operator type
• face: the face of the element to integrate
• fun_bar: incoming quantity that multiplies the functional from the left
• ±: PlusFunctor to add to flux_bar, MinusFunctor to subract

In/Outs

• flux_bar: result of the vector-Jacobian product.

SummationByParts.boundaryintegrate_rev!

This is the reverse differentiated version of boundaryintegrate!. See faceintegrate.jl for further details of the primal method. This function is differentiated with respect to the primal version's flux variable.

Inputs

• sbpface: an SBP AbstractFace type
• bndryfaces: list of boundary faces stored as an array of Boundarys
• res_bar: vector applied to the left of the (R^T*B) operator
• ±: PlusFunctor to add to flux_bar, MinusFunctor to subract

In/Outs

• flux_bar: result of the vector matrix product between (R^T*B) and res_bar

SummationByParts.boundaryFaceIntegrate_rev!

This is the reverse differentiated version of boundaryFaceIntegrate!. See faceintegrate.jl for further details of the primal method. This function is differentiated with respect to the primal version's flux variable.

Inputs

• sbpface: an SBP AbstractFace type
• face: the face of the element to integrate and project back to the element
• res_bar: vector applied to the left of the (R^T*B) operator
• ±: PlusFunctor to add to flux_bar, MinusFunctor to subract

In/Outs

• flux_bar: result of the vector matrix product between (R^T*B) and res_bar

SummationByParts.interiorfaceintegrate_rev!

This is the reverse differentiated version of interiorfaceintegrate!. See faceintegrate.jl for further details of the primal method. This function is differentiated with respect to the primal version's flux variable.

Inputs

• sbpface: an SBP AbstractFace type
• ifaces: list of element interfaces stored as an array of Interfaces
• res_bar: vector applied to the left of the (R^T*B) operator
• ±: PlusFunctor to add to flux_bar, MinusFunctor to subract

In/Outs

• flux_bar: result of the vector matrix product between (R^T*B) and res_bar

SummationByParts.interiorFaceIntegrate_rev!

This is the reverse differentiated version of interiorFaceIntegrate!. See faceintegrate.jl for further details of the primal method. This function is differentiated with respect to the primal version's flux variable.

Inputs

• sbpface: an SBP AbstractFace type
• ifaces: list of element interfaces stored as an array of Interfaces
• resL_bar: vector (left element) applied to the left of the (R^T*B) operator
• resR_bar: vector (right element) applied to the left of the (R^T*B) operator
• ±: PlusFunctor to add to flux_bar, MinusFunctor to subract

In/Outs

• flux_bar: result of the vector-matrix product

SummationByParts.boundaryFaceIntegrate_jac!

Given the face-node flux Jacobians, this method computes the Jacobian of the chain boundaryFaceInterpolate! –> face-flux evaluation –> boundaryFaceIntegrate! and adds the contributions to the relevant element's Jacobian matrix. Different methods are available depending on the rank of dfluxdu:

• For scalar fields, it is assumed that dfluxdu is a rank-1 array, with the

only dimension for the face local-node index. dresdu is the Jacobian of the element residual with respect to the state.

• For vector fields, dfluxdu is a rank-3 array, with the first and second

dimensions for indices of the vector field (at a particular face node) and the second dimension for the face local-node index. dresdu is a rank-4 array that is the Jacobian of the element residual with respect to the state; the first 2 dimensions are of size nvar = size(dfluxdu,1), while the third and fourth dimensions are of size sbp.numnodes.

Inputs

• sbpface: an SBP AbstractFace type
• face: the relevant face (index) of the element
• dfluxdu: array of the derivative of the flux w.r.t. the state trace
• ±: PlusFunctor to add to res, MinusFunctor to subract

In/Outs

• dresdu: Jacobian of element residual w.r.t. the state

SummationByParts.interiorFaceIntegrate_jac!

Given the face-node flux Jacobians, this method computes the Jacobian of the chain interiorFaceInterpolate! –> face-flux evaluation –> interiorFaceIntegrate! and adds the contributions to the adjacent elements' Jacobians matrices. Different methods are available depending on the rank of dfluxduL and dfluxduR:

• For scalar fields, it is assumed that dfluxduL and dfluxduR are rank-1

arrays, with the only dimension for the face local-node index. dresLduL is the Jacobian of the left-element residual with respect to the left-element solution, and similarly for dresLduR, dresRduL, dresRduR.

• For vector fields, dfluxduL and dfluxduR are rank-3 arrays, with the

first and second dimensions for indices of the vector field (at a particular face node) and the second dimension for the face local-node index. dresLduL is a rank-4 array that is the Jacobian of the left-element residual with respect to the left-element solution; the first 2 dimensions are of size nvar = size(dfluxduL,1), while the third and fourth dimensions are of size sbp.numnodes. Similarly for dresLduR, dresRduL, dresRduR.

Inputs

• sbpface: an SBP AbstractFace type
• ifaces: list of element interfaces stored as an array of Interfaces
• dfluxduL: array of the derivative of the flux w.r.t. the left trace
• dfluxduR: array of the derivative of the flux w.r.t. the right trace
• ±: PlusFunctor to add to res, MinusFunctor to subract

In/Outs

• dresLduL: Jacobian of left residual w.r.t. left state
• dresLduR: Jacobian of left residual w.r.t. right state
• dresRduL: Jacobian of right residual w.r.t. left state
• dresRduR: Jacobian of right residual w.r.t. right state

SummationByParts.volumeintegrate!

Applies the SBP mass matrix operator, H, to data in u and stores the result in res. 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
• u: the array that the operator is applied to
• ±: PlusFunctor to add to res, MinusFunctor to subract

In/Outs

• res: where the result of applying H to u is stored

SummationByParts.volumeIntegrateElement!

This is the single-element variant of volumeIntegrate!. Applies the SBP mass matrix operator, H, to data in u and stores the result in res. Different methods are available depending on the rank of u:

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

and only dimension for the local-node index.

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

index of the vector field, and the second dimension for the local-node 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
• u: the array that the operator is applied to
• ±: PlusFunctor to add to res, MinusFunctor to subract

In/Outs

• res: where the result of applying H to u is stored

SummationByParts.volumeintegrate_rev!

This is the reverse differentiated version of volumeintegrate!. See volumeintegrate.jl for further details of the primal method. This function is differentiated with respect to the primal version's u variable.

Inputs

• sbp: an SBP operator type
• res_bar: vector applied to the left of the H operator
• ±: PlusFunctor to add to res, MinusFunctor to subract

In/Outs

• u_bar: the result of the vector matrix product between H and res_bar

SummationByParts.volumeIntegrateElement_rev!

This is the reverse differentiated version of volumeIntegrateElement!. See volumeintegrate.jl for further details of the primal method. This function is differentiated with respect to the primal version's u variable.

Inputs

• sbp: an SBP operator type
• res_bar: vector applied to the left of the H operator
• ±: PlusFunctor to add to res, MinusFunctor to subract

In/Outs

• u_bar: the result of the vector matrix product between H and res_bar

SummationByParts.calcMappingJacobian!

Uses a given set of Lagrangian element nodes to determine an analytical (polynomial) mapping, and then uses this mapping to determine the Jacobian of the mapping. The approach varies depending on the dimension of the problem:

• For 2-dimensional problems the exact Jacobian of the mapping is used;
• For 3-dimensional problems an optimization problem is solved

Inputs

• sbp: an SBP operator type
• mapdegree: the polynomial degree of the mapping
• xref: Lagrangian nodes in reference space; [coord, Lagrangian node]
• xlag: Lagrangian nodes in physical space; [coord, Lagrangian node, element]
• Eone: Exone, Eyone (Ez*one); [sbp node, coord, element] see notes below

In/Outs

• xsbp: location of the SBP nodes in physical space; [coord, sbp node, element]
• dξdx: scaled Jacobian of mapping; [ref coord, phys coord, sbp node, element]
• jac: the determinant of the Jacobian; [sbp node, element]

Notes

The array Eone is the product of the boundary operators, in physical space, with the vector of ones (see Crean et al., Entropy-Conservative, Multidimensional Summation-By-Parts Discretization of the Euler Equations, as well as the test in test/test_mappingjacobian.jl). These products are used to define the metric invariants. They are not needed by the 2-dimensional code, and so this array can be passed empty in that case.

SummationByParts.calcMappingJacobianElement!

Single element variant of calcMappingJacobian!. Uses a given set of Lagrangian element nodes to determine an analytical (polynomial) mapping, and then uses this mapping to determine the Jacobian of the mapping. The approach varies depending on the dimension of the problem:

• For 2-dimensional problems the exact Jacobian of the mapping is used;
• For 3-dimensional problems an optimization problem is solved

Inputs

• sbp: an SBP operator type
• mapdegree: the polynomial degree of the mapping
• xref: Lagrangian nodes in reference space; [coord, Lagrangian node]
• xlag: Lagrangian nodes in physical space; [coord, Lagrangian node]
• Eone: Exone, Eyone (Ez*one); [sbp node, coord] see notes below

In/Outs

• xsbp: location of the SBP nodes in physical space; [coord, sbp node]
• dξdx: the scaled Jacobian of the mapping; [ref coord, phys coord, sbp node]
• jac: the determinant of the Jacobian; [sbp node]

Notes

The array Eone is the product of the boundary operators, in physical space, with the vector of ones (see Crean et al., Entropy-Conservative, Multidimensional Summation-By-Parts Discretization of the Euler Equations, as well as the test in test/test_mappingjacobian.jl). These products are used to define the metric invariants. They are not needed by the 2-dimensional code, and so this array can be passed empty in that case.

SummationByParts.mappingjacobian!

Evaluates the Jacobian of the mapping from face-reference coordinates to physical coordinates, as well as the determinant of the Jacobian.

Inputs

• sbpface: an SBP face operator type
• ifaces: list of element interfaces stored as an array of Interfaces
• x: the physical coordinates in [coord, node, elem] format

In/Outs

• dξdx: the Jacobian in [ξ coord, x coord, face node, L/R, face] format
• jac: the determinant in [face node, L/R, face] format

SummationByParts.mappingjacobian!

Deprecated:

Evaluates the (scaled) Jacobian of the mapping from reference coordinates to physical coordinates, as well as the determinant of the Jacobian. The values returned in dξdx are scaled by the determinant, so they have the same units as the boundary measure (i.e. length in 2D, or length^2 in 3D). This scaling is adopted, because conservation laws written in conservative form in the reference frame use the scaled Jacobian.

Inputs

• sbp: an SBP operator type
• x: the physical coordinates; 1st dim = coord, 2nd dim = node, 3rd dim = elem

In/Outs

• dξdx: the scaled Jacobian of the mapping; 1st dim = ref coord, 2nd dim = phys coord, 3rd dim = node, 3rd dim = elem
• jac: the determinant of the Jacobian; 1st dim = node, 2nd dim = elem

SummationByParts.calcMappingJacobian_rev!

Forms the reverse-mode of algorithmic differentiation product for the method calcMappingJacobian!. Specifically, it computes xlag_bar = xsbp_bar^T * ∂xsbp/∂xlag + dξdx_bar^T * ∂dξdx/∂xlag + jac_bar^T * ∂jac/∂xlag and Eone_bar = dξdx_bar^T * ∂dξdx/∂xlag where the various quantities are defined below (the _bar denotes the reverse mode). Note that the input and output order of the arguments follows that used in calcMappingJacobian!.

Inputs

• sbp: an SBP operator type
• mapdegree: the polynomial degree of the mapping
• xref: Lagrangian nodes in reference space; [coord, Lagrangian node]
• xsbp_bar: gradient w.r.t. SBP nodes; [coord, sbp node]
• dξdx: scaled Jacobian of mapping (as output from calcMappingJacobian!)
• dξdx_bar: gradient w.r.t. Jacobian; [ref coord, phys coord, sbp node, element]
• jac: the determinant of the Jacobian (as output from calcMappingJacobian!)
• jac_bar: gradient w.r.t. determinant of the Jacobian; [sbp node, element]

In/Outs

• xlag_bar: gradient w.r.t. Lagrangian nodes; [coord, Lagrangian node, element]
• Eone_bar: gradient w.r.t. Exone, Eyone (Ez*one); [sbp node, coord, element]

Notes

See calcMappingJacobian! for an explanation of Eone; it is only needed in the 3D case, but Eone_bar needs to be supplied in both 2D and 3D.

SummationByParts.mappingjacobian_rev!

Deprecated:

Forms the reverse-mode of algorithmic differentiation product for the method mappingjacobian!. Specifically, it computes x_bar = dξdx_bar^T * ∂dξdx/∂x + jac_bar^T * ∂jac/∂x where the various quantities are defined below (the _bar denotes the reverse mode). Note that the input and output order of the arguments follows that used in mappingjacobian!.

Inputs

• sbp: an SBP operator type
• x: the physical coordinates; 1st dim = coord, 2nd dim = node, 3rd dim = elem
• dξdx_bar: gradient w.r.t. Jacobian; [ref coord, phys coord, sbp node, element]
• jac_bar: gradient w.r.t. determinant of the Jacobian; [sbp node, element]

In/Outs

• x_bar: gradient w.r.t. SBP nodes; [coord, Lagrangian node, element]

SummationByParts.calcFaceNormals!

Uses a given set of Lagrangian face nodes to determine an analytical (polynomial) mapping, and then uses this mapping to determine the scaled face-normal vector.

Inputs

• sbpface: an SBP face operator type
• mapdegree: the polynomial degree of the mapping
• xref: Lagrangian nodes in reference space; [coord, Lag node]
• xlag: Lagrangian nodes in physical space; [coord, Lag node, face]

In/Outs

• xsbp: SBP-face nodes in physical space; [coord, sbp node, face]
• nrm: scaled face-normal at the sbpface nodes; [component, sbp node, face]

SummationByParts.calcFaceNormals!

Uses a given set of Lagrangian face nodes to determine an analytical (polynomial) mapping, and then uses this mapping to determine the scaled face-normal vector.

Inputs

• sbpface: an SBP face operator type
• mapdegree: the polynomial degree of the mapping
• xref: Lagrangian nodes in reference space; [coord, Lag node]
• xlag: Lagrangian nodes in physical space; [coord, Lag node, face]

In/Outs

• xsbp: SBP-face nodes in physical space; [coord, sbp node, face]
• nrm: scaled face-normal at the sbpface nodes; [component, sbp node, face]

SummationByParts.facenormal!

This is the single-face variant of calcFaceNormals!. Uses a given set of Lagrangian face nodes to determine an analytical (polynomial) mapping, and then uses this mapping to determine the scaled face-normal vector.

Inputs

• sbpface: an SBP face operator type
• mapdegree: the polynomial degree of the mapping
• xref: Lagrangian nodes in reference space; [coord, Lagrangian node]
• xlag: Lagrangian nodes in physical space; [coord, Lagrangian node]

In/Outs

• xsbp: location of the SBP-face nodes in physical space; [coord, sbp node]
• nrm: scaled face-normal at the sbpface nodes

SummationByParts.facenormal!

This is the single-face variant of calcFaceNormals!. Uses a given set of Lagrangian face nodes to determine an analytical (polynomial) mapping, and then uses this mapping to determine the scaled face-normal vector.

Inputs

• sbpface: an SBP face operator type
• mapdegree: the polynomial degree of the mapping
• xref: Lagrangian nodes in reference space; [coord, Lagrangian node]
• xlag: Lagrangian nodes in physical space; [coord, Lagrangian node]

In/Outs

• xsbp: location of the SBP-face nodes in physical space; [coord, sbp node]
• nrm: scaled face-normal at the sbpface nodes

SummationByParts.calcFaceNormals_rev!

This is the reverse differentiated version of calcFaceNormals!. See facenormal.jl for further details of the primal method. This function is differentiated with respect to the primal version's xsbp and nrm variables.

Note: xlag must be provided, but is only needed for the TetFace method.

Inputs

• sbpface: an SBP face operator type
• mapdegree: the polynomial degree of the mapping
• xref: Lagrangian nodes in reference space; [coord, Lag node]
• xlag: Lagrangian nodes in physical space; [coord, Lag node, face]
• xsbp_bar: multiplies d(xsbp)/d(xlag) from left; [coord, sbp node, face]
• nrm_bar: multiplies d(nrm)/d(xlag) from the left; [component, sbp node, face]

In/Outs

• xlag_bar: result of vector Jacobian product; [coord, Lag node, face]

SummationByParts.calcFaceNormals_rev!

This is the reverse differentiated version of calcFaceNormals!. See facenormal.jl for further details of the primal method. This function is differentiated with respect to the primal version's xsbp and nrm variables.

Note: xlag must be provided, but is only needed for the TetFace method.

Inputs

• sbpface: an SBP face operator type
• mapdegree: the polynomial degree of the mapping
• xref: Lagrangian nodes in reference space; [coord, Lag node]
• xlag: Lagrangian nodes in physical space; [coord, Lag node, face]
• xsbp_bar: multiplies d(xsbp)/d(xlag) from left; [coord, sbp node, face]
• nrm_bar: multiplies d(nrm)/d(xlag) from the left; [component, sbp node, face]

In/Outs

• xlag_bar: result of vector Jacobian product; [coord, Lag node, face]

SummationByParts.facenormal_rev!

This is the reverse differentiated version of facenormal!. See facenormal.jl for further details of the primal method. This function is differentiated with respect to the primal version's xsbp and nrm variables.

Note: xlag must be provided, but is only needed for the TetFace method.

Inputs

• sbpface: an SBP face operator type
• mapdegree: the polynomial degree of the mapping
• xref: Lagrangian nodes in reference space; [coord, Lagrangian node]
• xlag: Lagrangian nodes in physical space; [coord, Lagragnian node]
• xsbp_bar: multiplies d(xsbp)/d(xlag) from left; [coord, sbp node]
• nrm_bar: multiplies d(nrm)/d(xlag) from the left; [component, sbp node]

In/Outs

• xlag_bar: result of vector Jacobian product; [coord, Lagrangian node]

SummationByParts.facenormal_rev!

This is the reverse differentiated version of facenormal!. See facenormal.jl for further details of the primal method. This function is differentiated with respect to the primal version's xsbp and nrm variables.

Note: xlag must be provided, but is only needed for the TetFace method.

Inputs

• sbpface: an SBP face operator type
• mapdegree: the polynomial degree of the mapping
• xref: Lagrangian nodes in reference space; [coord, Lagrangian node]
• xlag: Lagrangian nodes in physical space; [coord, Lagragnian node]
• xsbp_bar: multiplies d(xsbp)/d(xlag) from left; [coord, sbp node]
• nrm_bar: multiplies d(nrm)/d(xlag) from the left; [component, sbp node]

In/Outs

• xlag_bar: result of vector Jacobian product; [coord, Lagrangian node]

SummationByParts.levenberg_marquardt

Levenberg-Marquardt Algorithm (LMA)

Inputs

• fun: function to be optimized
• cub: cubature data
• q: the degree of the quadrature rule
• mask: a vector of the index of the parameters to be determined
• np: number of particles
• xinit: initial parameter guess
• xL: lower bound on the parameters
• xR: upper bound on the parameters
• nu: parameter controling exploration (spliting between Newton's and steepest decent methods)
• maxiter: maximum number of iterations
• tol: tolerance to stop iteration
• verbose: boolean to print results of iteration periodically

Outputs

• fmin: the optimized function value
• v: the minimizer (solution)
• iter: number of itrations

SummationByParts.pso

Particle Swarm Optimization (PSO) Algorithm

Inputs

• fun: function to be optimized
• ne: number of parameters to be found
• cub: cubature data
• q: the degree of the quadrature rule
• mask: a vector of the index of the parameters to be determined
• np: number of particles
• xinit: initial parameter guess
• xL: lower bound on the parameters
• xR: upper bound on the parameters
• delta1: purturbation constant that determines how much the parameters are perturbed in the case of stagnation with fmin <= 0.1
• delta2: purturbation constant that determines how much the parameters are perturbed in the case of stagnation with fmin > 0.1
• maxiter: maximum number of iterations (default is 1000)
• tol: tolerance to stop iteration (default is 1e-14)
• save_iter: boolean to save results of each iteration
• verbose: boolean to print results of iteration periodically

Outputs

• fmin: the optimized function value
• xmin: the minimizer (solution)
• f_all: all function evaluations

OrthoPoly.diffjacobipoly{T}

Evaluate the first derivative of a Jacobi Polynomial at some points.

Inputs

• x: points at which to evaluate polynomial derivative
• alpha,beta: define the type of Jacobi Polynomial (alpha + beta != 1)
• N: polynomial degree

Outputs

• dP - derivative of polynomial evaluated at x

OrthoPoly.diffproriolpoly{T}

Evaluate the derivatives of a Proriol orthogonal polynomial basis function on the right tetrahedron.

Notes: the derivatives are computed using the complex-step method (since there are many outputs and only 3 inputs); therefore, a different method should be used for verification of this method.

Inputs

• x,y,z: locations at which to evaluate the derivative
• i,j,k: index triple that defines the basis function to differentiate; see Hesthaven and Warburton's Nodal DG book, for example, for a reference.

Outputs

• dPdx,dPdy,dPdz: derivatives of basis function at (x,y,z)

OrthoPoly.diffproriolpoly

Evaluate the derivatives of a Proriol orthogonal polynomial basis function on the right triangle.

Inputs

• x,y: locations at which to evaluate the derivative
• i,j: index pair that defines the basis function to differentiate; see Hesthaven and Warburton's Nodal DG book, for example, for a reference.

Outputs

• dPdx,dPdy: derivative of basis function at (x,y)

OrthoPoly.jacobipoly{T}

Evaluate a Jacobi polynomial at some points. Based on JacobiP in Hesthaven and Warburton's nodal DG book.

Inputs

• x: points at which to evaluate polynomial
• alpha,beta: define the type of Jacobi Polynomial (alpha + beta != 1)
• N: polynomial degree

Outputs

• P: the polynomial evaluated at x

OrthoPoly.lglnodes

Computes the Legendre-Gauss-Lobatto (LGL) quadrature nodes and weights on the interval [-1,1]. The LGL nodes are the zeros of (1-x^2)*P'N(x), where PN(x) denotes the Nth Legendre polynomial.

Reference: C. Canuto, M. Y. Hussaini, A. Quarteroni, T. A. Tang, "Spectral Methods in Fluid Dynamics," Section 2.3. Springer-Verlag 1987

Inputs

• N: highest degree (number of nodes = N+1)
• T: number type

Outputs

• x: the LGL nodes
• w: the LGL weights

Julia version adapted from Matlab code written by Greg von Winckel - 04/17/2004 Contact: gregvw@chtm.unm.edu

OrthoPoly.lgnodes

Computes the Legendre-Gauss (LG) quadrature nodes and weights on the interval [-1,1]. The LG nodes are the zeros of PN(x), where PN(x) denotes the Nth Legendre polynomial.

Inputs

• N: number of nodes
• T: number type

Outputs

• x: the LG nodes
• w: the LG weights

Julia version adapted from Matlab code written by Greg von Winckel - 02/25/2004 Contact: gregvw@chtm.unm.edu

OrthoPoly.proriolpoly{T}: method for a right tetrahedron

Evaluate Proriol orthogonal polynomial basis function on the right tetrahedron.

Inputs

• x,y,z: locations at which to evaluate the polynomial
• i,j,k: index triple that defines the basis function to evaluate; see Hesthaven and Warburton's Nodal DG book, for example, for a reference.

Outputs

• P: basis function at (x,y,z)

OrthoPoly.proriolpoly{T}: method for right triangle

Evaluate Proriol orthogonal polynomial basis function on the right triangle.

Inputs

• x,y: locations at which to evaluate the polynomial
• i,j: index pair that defines the basis function to evaluate; see Hesthaven and Warburton's Nodal DG book, for example, for a reference.

Outputs

• P: basis function at (x,y)

OrthoPoly.vandermonde{T}

Evaluate the Vandermonde matrix using the Proriol polynomials on the right tetrahedron.

Inputs

• p: the maximum total degree of the polynomial
• x,y,z: locations at which to evaluate the derivative
• compute_grad: indicates whether to compute the gradient of the Vandermonde matrix

Outputs

• V: the Vandermonde matrix
• Vdx,Vdy,Vdz: derivatives of the Vandermonde matrix at (x,y,z)

OrthoPoly.vandermonde{T}

Evaluate the Vandermonde matrix using the Proriol polynomials on the right triangle.

Inputs

• p: the maximum total degree of the polynomial
• x,y: locations at which to evaluate the derivative
• compute_grad: indicates whether to compute the gradient of the Vandermonde matrix

Outputs

• V: the Vandermonde matrix
• Vdx,Vdy: derivatives of the Vandermonde matrix at (x,y)

OrthoPoly.vandermonde_arnoldi{T}

Evaluate the Vandermonde with Arnoldi on the right tetrahedron.

Inputs

• p: the maximum total degree of the polynomial
• x,y,z: locations at which to evaluate the derivative
• compute_grad: indicates whether to compute the gradient of the Vandermonde matrix

Outputs

• V: the Vandermonde matrix
• Vdx,Vdy,Vdz: derivatives of the Vandermonde matrix at (x,y,z)
• Hes: the Hessenberg matrix

OrthoPoly.vandermonde_arnoldi{T}

Evaluate the Vandermonde with Arnoldi on the right triangle.

Inputs

• p: the maximum total degree of the polynomial
• x,y: locations at which to evaluate the derivative
• compute_grad: indicates whether to compute the gradient of the Vandermonde matrix

Outputs

• V: the Vandermonde matrix
• Vdx, Vdy: derivatives of the Vandermonde matrix at (x,y)
• Hes: the Hessenberg matrix

OrthoPoly.vandermonde_arnoldi{T}

Evaluate the Vandermonde with Arnoldi on the line.

Inputs

• p: the maximum total degree of the polynomial
• x: locations at which to evaluate the derivative
• compute_grad: indicates whether to compute the gradient of the Vandermonde matrix

Outputs

• V: the Vandermonde matrix
• Vdx: derivatives of the Vandermonde matrix at (x)
• Hes: the Hessenberg matrix

OrthoPoly.vandermonde_full{T}

Evaluate the Vandermonde and complements it with the nullspace to create a full rank matrix.

Inputs

• p: the maximum total degree of the polynomial
• x: locations at which to evaluate the derivative

Outputs

• Vfull: the full rank Vandermonde matrix

OrthoPoly.vandermonde_monomial{T}

Evaluate the Vandermonde matrix using monomials on the right tetrahedron.

Inputs

• p: the maximum total degree of the polynomial
• x,y,z: locations at which to evaluate the derivative
• compute_grad: indicates whether to compute the gradient of the Vandermonde matrix
• compute_integ: indicates whether to compute the integral of the monomial basis functions

Outputs

• V: the Vandermonde matrix
• Vdx,Vdy, Vdz: derivatives of the Vandermonde matrix at (x,y,z)
• Vinteg: the integral of each basis function in the Vandermonde matrix

OrthoPoly.vandermonde_monomial{T}

Evaluate the Vandermonde matrix using monomials on the right triangle.

Inputs

• p: the maximum total degree of the polynomial
• x,y: locations at which to evaluate the derivative
• compute_grad: indicates whether to compute the gradient of the Vandermonde matrix
• compute_integ: indicates whether to compute the integral of the monomial basis functions

Outputs

• V: the Vandermonde matrix
• Vdx,Vdy: derivatives of the Vandermonde matrix at (x,y)
• Vinteg: the integral of each basis function in the Vandermonde matrix

OrthoPoly.vandermonde_monomial{T}

Evaluate the Vandermonde matrix using monomials on the line.

Inputs

• p: the maximum total degree of the polynomial
• x: locations at which to evaluate the derivative
• compute_grad: indicates whether to compute the gradient of the Vandermonde matrix
• compute_integ: indicates whether to compute the integral of the monomial basis functions

Outputs

• V: the Vandermonde matrix
• Vdx: derivatives of the Vandermonde matrix at (x)
• Vinteg: the integral of each basis function in the Vandermonde matrix

SBP.getLineSegSBPLobbato

Returns Gauss-Lobbato type elements, that have nodes on the element boundary

Inputs

• degree: maximum polynomial degree for which the derivatives are exact
• Tsbp: floating point type used for the operators

Returns

• sbp: a Gauss-Lobbato operator of the appropriate degree

SBP.getLineSegSBPLegendre

Returns Gauss-Legendre type elements, that do not have nodes on the element boundary

Inputs

• degree: maximum polynomial degree for which the derivatives are exact
• Tsbp: floating point type used for the operators

Returns

• sbp: a Legendre-Gauss operator of the appropriate degree

SBP.getTriSBPGamma

Returns SBP-Gamma type elements, that have nodes on the element boundary

Inputs

• degree: maximum polynomial degree for which the derivatives are exact
• Tsbp: floating point type used for the operators

Returns

• sbp: an SBP-Gamma operator of the appropriate degree

SBP.getTriSBPOmega

Returns SBP-Omega type elements, that have no nodes on the element boundary

Inputs

• degree: maximum polynomial degree for which the derivatives are exact
• Tsbp: floating point type used for the operators

Returns

• sbp: an SBP-Omega operator of the appropriate degree

SBP.getTriSBPDiagE

Returns SBP-DiagE type elements, whose boundary nodes are positioned at cubature points

Inputs

• degree: maximum polynomial degree for which the derivatives are exact
• Tsbp: floating point type used for the operators

Returns

• sbp: an SBP-DiagE operator of the appropriate degree

SBP.getTetSBPGamma

Returns SBP-Gamma type elements, that have nodes on the element boundary

Inputs

• degree: maximum polynomial degree for which the derivatives are exact
• Tsbp: floating point type used for the operators

Returns

• sbp: an SBP-Gamma operator of the appropriate degree

SBP.getTetSBPOmega

Returns SBP-Omega type elements, that have no nodes on the element boundary

Inputs

• degree: maximum polynomial degree for which the derivatives are exact
• Tsbp: floating point type used for the operators

Returns

• sbp: an SBP-Omega operator of the appropriate degree

SBP.getTetSBPDiagE

Returns SBP-DiagE type elements, whose boundary nodes are positioned at cubature points

Inputs

• degree: maximum polynomial degree for which the derivatives are exact
• Tsbp: floating point type used for the operators

Returns

• sbp: an SBP-DiagE operator of the appropriate degree

SBP.getLineSegFace

Returns a trival face for line-segment elements.

Inputs

• degree: face integration is exact for polys of degree 2*degree
• volcub: cubature rule for the associated "volume"
• vtx: vertices of the line-segment

Returns

• sbpface: an SBP face type for line-segment elements

SBP.TriFace

Defines a face between two TriSBP operators with the same cubature nodes

Fields

• degree : face integration is exact for polys of degree 2*degree
• numnodes : number of cubature nodes
• stencilsize : number of nodes in the reconstruction stencil
• dstencilsize : number of nodes in the derivative operator stencils
• cub : a symmetric cubature type for triangle faces (i.e. edges)
• vtx : the vertices of the face in reference space, [-1,1]
• wface : mass matrix (quadrature) for the face
• interp[:,:] : volume-to-face-nodes reconstruction operator
• perm[:,:] : permutation for volume nodes so interp can be used on all sides
• deriv[:,:] : derivative operators for face-based coordinate system
• dperm[:,:] : permutation for volume nodes so deriv can be used on all sides
• nbrperm[:,:] : permutation for face nodes on neighbour element

SBP.getTriFaceForDiagE

Returns a quadrature that can be used to construct SBP operators with diagonal boundary operators, E.

Inputs

• degree: face integration is exact for polys of degree 2*degree
• volcub: cubature rule for the associated volume
• vtx: vertices of the triangle
• vertices: (optional) if true, the face cubature includes the vertices

Returns

• sbpface: an SBP face type for triangle elements

SBP.TetFace

Defines a face between two TetSBP operators with the same cubature nodes

Fields

• degree : face integration is exact for polys of degree 2*degree
• numnodes : number of cubature nodes
• stencilsize : number of nodes in the reconstruction stencil
• dstencilsize : number of nodes in the derivative operator stencils
• cub : a symmetric cubature type for tetrahedral faces (i.e. triangles)
• vtx : the vertices of the face in the reference space of the face
• wface : mass matrix (quadrature) for the face
• interp[:,:] : volume-to-face-nodes reconstruction operator
• perm[:,:] : permutation for volume nodes so interp can be used on all sides
• deriv[:,:] : derivative operators for face-based coordinate system
• dperm[:,:] : permutation for volume nodes so deriv can be used on all sides
• nbrperm[:,:] : permutation for face nodes on neighbour element

SBP.getTetFaceForDiagE

Returns a quadrature that can be used to construct SBP operators with diagonal boundary operators, E.

Inputs

• degree: face integration is exact for polys of degree 2*degree+1
• volcub: cubature rule for the associated volume
• vtx: vertices of the triangle

Returns

• sbpface: an SBP face type for triangle elements

SymCubatures.LineSymCub

Defines a symmetric quadrature rule on the interval [-1,1]. Current choices are Legendre-Gauss-Lobatto (LGL) or Legendre-Gauss (LG) rules.

Fields

• numparams : total number of nodal degrees of freedom
• numweights : total number of unique weights
• numnodes : total number of nodes
• vertices : if true, vertices (ends of interval) are in the set of nodes
• centroid : if true, centroid is present in set of nodes
• numedge : number of unique edge parameters
• numsym : number of node sets in each symmetry group (2 for [-1,1])
• params : the actual values of the orbit nodal parameters
• weights : values of the unique weights

SymCubatures.OperType

OperType is an abtract type that defines the type of SBP operators It has the subtypes Omega, Gamma, and DiagE.

SymCubatures.PointSymCub

Defines the trivial point cubature for uniformity across methods.

Fields

• numparams : total number of nodal degrees of freedom
• numweights : total number of unique weights
• numnodes : total number of nodes
• vertices : if true, vertices (ends of interval) are in the set of nodes
• centroid : if true, centroid is present in set of nodes
• numedge : number of unique edge parameters
• numsym : number of node sets in each symmetry group (0 for [1])
• params : the actual values of the orbit nodal parameters
• weights : values of the unique weights

SymCubatures.SymCub

SymCub is an parametric abstract type that defines cubatures for symmetric nodal distributions. It is parameterized on T in order to allow for future implementations of arbitrary precision types. The parameterization also permits the use of the complex-step method for verification.

SymCubatures.TetSymCub

Used to define symmetric cubature rules on the tetrahedron. The params array determines the position of the parameterized nodes, and the weights array determines the value of the weight for each symmetric orbit. Note that boolean fields are used to activate some degenerate symmetry orbits. For example, vertices are a special case of several of the orbits, and should be activated by setting vertices=true rather than relying on a specific value of a parameter.

Fields

• numparams : total number of nodal degrees of freedom
• numweights : total number of unique weights
• numnodes : total number of nodes
• vertices : if true, vertices are present in the set of nodes
• midedges : if true, edge midpoints are present in set of nodes
• centroid : if true, centroid is present in set of nodes
• facecentroid : if true, face centroids are present in the set of nodes
• numedge : number of unique edge parameters
• numfaceS21 : number of S21 face orbits (same tri orbit on face)
• numfaceS111 : number of S111 face orbits (same tri orbit on face)
• numS31 : number of S31 orbits (vertex to opposite face)
• numS22 : number of S22 orbits
• numS211: number of S211 orbits
• numS1111: number of S1111 orbits
• numsym : number of node sets in each symmetry group (5 groups for Tet)
• params : the actual values of the orbit nodal parameters
• weights : values of the unique weights

SymCubatures.TriSymCub

Used to define symmetric cubature rules on the triangle. The params array determines the position of the parameterized nodes, and the weights array determines the value of the weight for each symmetric orbit. Note that boolean fields are used to activate some degenerate orbits. For example, vertices are a special case of several of the orbits, and should be activated by setting vertices=true rather than relying on a specific value of a parameter.

Fields

• numparams : total number of nodal degrees of freedom
• numweights : total number of unique weights
• numnodes : total number of nodes
• vertices : if true, vertices are present in the set of nodes
• midedges : if true, edge midpoints are present in set of nodes
• centroid : if true, centroid is present in set of nodes
• numedge : number of unique edge parameters
• numS21 : number of S21 orbits (vertex to opposite face)
• numS111 : number of S111 orbits
• numsym : number of node sets in each symmetry group (3 groups for Tri)
• params : the actual values of the orbit nodal parameters
• weights : values of the unique weights

SymCubatures.calcjacobian

Returns the Jacobian of the nodal coordinates and weights with respect to their parameters. In other words, returns the block-rectangular matrix [Jcoords, 0; 0, Jweights], where Jcoords is the Jacobian of the coordinates, Jweights is the Jacobian of the weights, and the 0s indicate zero blocks of the appropriate size.

Inputs

• cub: symmetric cubature rule
• vtx: vertices of the cubature domain

Outputs

• Jac: Jacobian of the nodal coordinates and weights

SymCubatures.calcjacobianofnodes

Returns the Jacobian of the nodes with respect to the orbit parameters.

Notes: Jac stores all the x-coordinate Jacobians first, then y (then z)

Inputs

• cub: symmetric cubature rule
• vtx: vertices of the cubature domain

Outputs

• Jac: Jacobian of the mapping from node parameters to nodes

SymCubatures.calcjacobianofweights

Returns the Jacobian of the nodal weights with respect to the unique weights. The resulting Jacobian is a rectangular matrix of ones and zeros that indicates the mapping from the unique weights to the nodal weights.

Inputs

• cub: symmetric cubature rule

Outputs

• Jac: Jacobian of the mapping from (unique) weights to nodal weights

SymCubatures.calcnodes

Use the orbital parameter values to compute a cubature's nodal coordinates. The second dimension of the vtx array of vertices does not need to match the dimension as the cubature; for example, a line quadrature can be over a line in 1D, 2D, 3D or ND, and a triangle cubature can be in 2D, 3D, or ND, etc.

Inputs

• cub: symmetric cubature rule
• vtx: vertices that define the domain

Outputs

• x: cubature's nodal coordinates (potentially in a subspace)

SymCubatures.calcweights

Map the unique cubature weights to the weights of all nodes.

Inputs

• cub: symmetric cubature rule

Outputs

• w: cubature's weights at all nodes

SymCubatures.calcweights_symgroup

Map the cubature weights of all nodes to the unique weights of their symmetry groups.

Inputs

• cub: symmetric cubature rule

Outputs

• w: cubature's weights of the symmetry group

SymCubatures.findleftperm!

For a matrix A, we are given a right permutation of the columns, A[:,permR]. This function attempts to find the left permultation of rows such that A[permL,:] = A[:,permR]

Inputs

• A: a rectangular matrix for which the left permutation is sought
• permR: the given right permutation of the columns

Outputs

• permL: the left permutation of the rows, if it exists

Returns

• true if the permutation exists, false otherwise

SymCubatures.getInternalParamMask

Returns the set of parameter indices corresponding to internal nodes; this is useful when finding cubature rules for which we wish to fix the boundary nodes and only allow the internal nodes to move.

Inputs

• cub: symmetric cubature rule

Returns

• mask: integer array of parameter indices associated with internal nodes

SymCubatures.getbndrynodeindices

Returns the indices of the nodes that lie on the boundary, in their natural order. See getfacenodeindices for a method returns node indices for each face.

Inputs

• cub: a symmetric cubature rule whose boundary-node indices are sought

Outputs

• bndryindices: indicies of nodes that lie on boundary

SymCubatures.getfacebasedpermutation

Returns a permutation of the volume nodes (or a subset of them) for each face, such that the same face operator can be applied to all faces. This is useful for volume-to-face interpolation or differentiation.

Inputs

• cub: a symmetric cubature rule for which a face-based permutation is sought
• faceonly: if true, only face nodes are used in the permutation.

Outputs

• perm: permutation of the volume nodes for each face

SymCubatures.getfacenodeindices

Returns the indices of the nodes that lie on each face. See getbndrynodeindices for a method that returns a single array of boundary nodes.

Inputs

• cub: a symmetric cubature rule whose boundary-node indices are sought

Outputs

• bndryindices: indicies of nodes that lie on boundary; there is a separate column of indices for each edge/face.

SymCubatures.getfacevertexindices

Returns the indices of the vertices that make up each face. This is useful when building nodes on a given face using Barycentric coordinates

Inputs

• cub: a symmetric cubature rule whose face vertices are sought

Outputs

• facevtx: subarray facevtx[:,f] lists the vertices of face f

SymCubatures.getinteriornodeindices

Returns the indices of the nodes that are strictly interior.

Inputs

• cub: a symmetric cubature rule whose interior-node indices are sought

Outputs

• indices: indicies of nodes that are strictly interior.

SymCubatures.getneighbourpermutation

At element interfaces, the cubature nodes of the common face will not match when natural ordering is provided. This routine produces the permutation that makes the 'right' element's nodes match the 'left' element's nodes. The permutation depends on the face dimension:

• For line segment faces, i.e. points, a trival permutation is returned
• For triangle faces, i.e. line segments, there is only one possible orientation.
• For tetrahedral faces, i.e. triangles, there are three orientations:

1. The "1" vertex from each element's face is coincident;
2. The "1" vertex from face 1 coincides with "2" vertex from face 2;
3. The "1" vertex from face 1 coincides with "3" vertex from face 2.

Inputs

• cub: a symmetric cubature rule for which the permutation is sought

Outputs

• perm: permutation of the interface nodes for each possible orientation

SymCubatures.getnumboundarynodes

Returns the number of (explicit) boundary nodes

Notes: if the parameter value for an internal orbit is such that the corresponding node lies on the boundary, this node is NOT included in the boundary-node count returned.

Inputs

• cub: symmetric cubature rule

Outputs

• numboundary: number of boundary nodes

SymCubatures.getnumfacenodes

Returns the number of nodes on an individual face of the element.

Inputs

• cub: symmetric cubature rule

Outputs

• numfacenodes: number of nodes on a face

SymCubatures.getpermutation

Returns a permutation of the cubature nodes based on a reordering of the vertices from their canonical orientation. That is, finds the node ordering such that the new nodes have the same order as the old nodes, but relative to vtxperm. This is useful for face-based operations.

Inputs

• cub: a symmetric cubature rule for which the permutation is sought
• vtxperm: the permutation that is applied to the vertices

Outputs

• perm: permutation of the cubature nodes

SymCubatures.setparams!

Sets the nodal parameters for any parameterized symmetry orbits in the cubature.

Inputs

• params: parameter values

In/Outs

• cub: symmetric cubature rule whose nodal parameters are being updated

SymCubatures.setweights!

Sets a cubature's (unique) weights.

Inputs

• weights: cubature weights grouped by orbit

In/Outs

• cub: symmetric cubature rule whose weights are being updated

SBP.AbstractSBP

AbstractSBP is a parametric abstract type that defines summation-by-parts finite-difference operators.

SBP.LineSegSBP

Defines diagonal-norm SBP first-derivative operators on a line segment.

Fields

• degree : maximum polynomial degree for which the derivatives are exact
• numnodes : number of nodes on line segment required for these operators
• cub : a symmetric cubature type for line segments (usually LG or LGL)
• vtx : vertices of the reference element in computational space
• w : cubature weights, i.e. the diagonal SBP norm, stored as an array
• Q[:,:,1] : discrete stiffness matrix operator

SBP.TriSBP

Defines diagonal-norm SBP first-derivative operators on a right-triangle.

Fields

• degree : maximum polynomial degree for which the derivatives are exact
• numnodes : number of nodes in the triangle required for these operators
• cub : a symmetric cubature type for triangles
• vtx : vertices of the reference element in computational space
• w : cubature weights, i.e. the diagonal SBP norm, stored as an array
• Q[:,:,i] : discrete stiffness matrix operator in ith coordinate direction

SBP.SparseTriSBP

Defines diagonal-norm SBP first-derivative operators on a right-triangle using a cubature rule that is greater than 2*p-1. This provides additional flexiblity in the SBP operator that is used to make a sparse S.

Fields

• degree : maximum polynomial degree for which the derivatives are exact
• numnodes : number of nodes in the triangle required for these operators
• cub : a symmetric cubature type for triangles
• vtx : vertices of the reference element in computational space
• w : cubature weights, i.e. the diagonal SBP norm, stored as an array
• Q[:,:,i] : discrete stiffness matrix operator in ith coordinate direction

SBP.TetSBP

Defines diagonal-norm SBP first-derivative operators on a right-tetrahedron.

Fields

• degree : maximum polynomial degree for which the derivatives are exact
• numnodes : number of nodes in the tetrahedron required for these operators
• cub : a symmetric cubature type for tetrahedra
• vtx : vertices of the reference element in computational space
• w : cubature weights, i.e. the diagonal SBP norm, stored as an array
• Q[:,:,i] : discrete stiffness matrix operator in ith coordinate direction

SBP.SparseTetSBP

Defines diagonal-norm SBP first-derivative operators on a right-tetrahedron using a cubature rule that is greater than 2*p-1. This provides additional flexiblity in the SBP operator that is used to make a sparse S.

Fields

• degree : maximum polynomial degree for which the derivatives are exact
• numnodes : number of nodes in the tetrahedron required for these operators
• cub : a symmetric cubature type for tetrahedra
• vtx : vertices of the reference element in computational space
• w : cubature weights, i.e. the diagonal SBP norm, stored as an array
• Q[:,:,i] : discrete stiffness matrix operator in ith coordinate direction

SBP.AbstractFace

AbstractFace is a parametric abstract type that defines face-based data and operations (e.g. volume-to-face reconstruction, face integration, etc) for summation-by-parts finite-difference operators.

SBP.DenseFace

DenseFace is a parametric abstract type that defines face-based data and operations (e.g. volume-to-face reconstruction, face integration, etc) for summation-by-parts finite-difference operators. This is a subtype for which interpolation is a dense matrix.

SBP.LineSegFace

Defines a "face" between two LineSegSBP operators with the same cubature nodes.

Fields

• degree : face integration is exact for polys of degree 2*degree
• numnodes : number of cubature nodes (always 1)
• stencilsize : number of nodes in the reconstruction stencil
• dstencilsize : number of nodes in the derivative operator stencils
• cub : a symmetric cubature type for line-segment faces (i.e. points)
• vtx : the vertices of the face in reference space, [-1]
• wface : mass matrix (quadrature) for the face (always 1.0)
• interp[:,:] : volume-to-face-nodes reconstruction operator
• perm[:,:] : permutation for volume nodes so interp can be used on both sides
• deriv[:,:] : derivative operators for face-based coordinate system
• dperm[:,:] : permutation for volume nodes so deriv can be used on both sides
• nbrperm[:,:] : permutation for face nodes on neighbour element

SBP.SparseFace

SparseFace is a parametric abstract type that defines face-based data and operations (e.g. volume-to-face reconstruction, face integration, etc) for summation-by-parts finite-difference operators in the case where the face-cubature nodes and volume nodes coincide (i.e. diagonal E operators).

SBP.TriSparseFace

Defines a face between two TriSBP operators with the same cubature nodes, in which the face-cubature nodes and volume nodes coincide (i.e. diagonal E operators).

Fields

• degree : face integration is exact for polys of degree 2*degree
• numnodes : number of cubature nodes
• dstencilsize : number of nodes in the derivative operator stencils
• cub : a symmetric cubature type for triangle faces (i.e. edges)
• vtx : the vertices of the face in reference space, [-1,1]
• wface : mass matrix (quadrature) for the face
• perm[:,:] : maps volume nodes to face nodes on each side
• deriv[:,:] : derivative operators for face-based coordinate system
• dperm[:,:] : permutation for volume nodes so deriv can be used on all sides
• nbrperm[:,:] : permutation for face nodes on neighbour element

SBP.TetSparseFace

Defines a face between two TetSBP operators with the same cubature nodes, in which the face-cubature nodes and volume nodes coincide (i.e. diagonal E operators).

Fields

• degree : face integration is exact for polys of degree 2*degree
• numnodes : number of cubature nodes
• dstencilsize : number of nodes in the derivative operator stencils
• cub : a symmetric cubature type for tetrahedral faces (i.e. triangles)
• vtx : the vertices of the face in the reference space of the face
• wface : mass matrix (quadrature) for the face
• perm[:,:] : permutation for volume nodes to face nodes on each side
• deriv[:,:] : derivative operators for face-based coordinate system
• dperm[:,:] : permutation for volume nodes so deriv can be used on all sides
• nbrperm[:,:] : permutation for face nodes on neighbour element

ODLCommonTools.Boundary

Used to identify boundary faces in a finite-element grid.

Fields

• element : index of the element to which the boundary face belongs
• face : the face index of the boundary (local index to the element)

Example

To mark face 2 of element 7 to be a boundary face, use Boundary(7,2)

ODLCommonTools.Interface

Used to identify interfaces between elements in a finite-element grid.

Fields

• elementL : index of the so-called left element in the pair
• elementR : index of the so-called right element in the pair
• faceL : the face index of the interface with respect to the left element
• faceR : the face index of the interface with respect to the right element
• orient : orientation of the 'right' element relative to the 'left'

Example

Consider an interface between elements 2 and 5. Suppose the interface is on face 1 of element 2 and face 3 of element 5. Furthermore, suppose element 5 has orientation 1 relative to element 1 (defintion of orientation TBD). This can be indicated as Interface(2,5,1,3,1)

SummationByParts.getNumFaceNodes

Returns the number of SBP element nodes on a face.

Inputs

• sbp: an SBP operator

Returns

• numfacenodes: number of nodes on (one) face of the element

SummationByParts.getnbrnodeindex

Returns the face-node index on face.faceR equivalent to index i on face.faceL.

Inputs

• sbp: an SBP operator
• face: an element interface
• i: face-node index on face.faceL

Returns

• j: face-node index on face.faceR

SummationByParts.calcnodes

This function returns the node coordinates for an SBP operator. It basically calls calcnodes for the underlying SymCubature. This function assumes the element mapping is linear, i.e. edges are lines.

Inputs

• sbp: an SBP operator
• vtx: the vertices that define the element

Outputs

• x: the node coordinates; 1st dimension is the coordinate, the second the node

Example

  # define a third-order accurate SBP on triangles
  sbp = TriSBP{Float64}(degree=2)
  # build a simple 2-element grid on a square domain
  x = zeros(Float64, (2,sbp.numnodes,2))
  vtx = [0. 0.; 1. 0.; 0. 1.]
  x[:,:,1] = calcnodes(sbp, vtx)
  vtx = [1. 0.; 1. 1.; 0. 1.]
  x[:,:,2] = calcnodes(sbp, vtx)

SummationByParts.calcminnodedistance

Returns the minimum distance between distinct nodes on an element with straight sides

Inputs

• sbp: an SBP operator
• vtx: the vertices that define the element

Returns

• mindist: the minimum distance between distinct nodes

SummationByParts.buildinterpolation

Builds a matrix operator that can reconstruct a field located at the sbp nodes to an auxlliary set of nodes.

Inputs

• sbp: an SBP operator
• xinterp: points to interpolate to in ref coords, size = [ndim,numpoints]
• d=sbp.degree: (optional) interpolation is exact for degree d polys

Returns

• R: the interpolation operator, size = [numpoints, sbp.numnodes]

SummationByParts.permuteinterface!

For a given array of values on a faces, permutes the node values (in place) to be in the orientation specified by the orient field of the corresponding Interface. Methods are available for scalar and vector fields.

Inputs

• sbp: an SBPFace operator
• ifaces: an array of Interfaces

In/Outs

• uface: the array of face values, must have dimensions [n x sbpface.numnodes x length(ifaces] if 3D array (where n is arbitrary) or [sbpface.numnodes x length(ifaces] if 2D. The permutation is applied to the second array dimension for the 3D case or the first dimension in the 2D case.

SummationByParts.permuteface!

This function applys a permutation to the data on a particular face.

Inputs

• permvec: vector specifying the permutation to apply
• workarr: a temporary array, same size as face_data, that is overwritten during the computation

In/Outs

• face_data: an N-D array containing the data to be pemuted, where the permutation is applied to the second dimension of the array for N=2 and the first dimension if N=1. N > 2 is not currently supported

SummationByParts.basispursuit!

Finds an approximate solution to the underdetermined problem Ax = b that is sparse using the alternating direction method of multipliers (ADMM).

Inputs

• A: matrix in the linear equation that must be satisfied
• b: vector in the linear equation that must be satisfied
• rho (optional) : augmented Lagrangian parameter
• alpha (optional) : over-relaxation parameter (usually between 1.0 and 1.8)
• hist (optional): output the convergence history to the terminal
• abstol (optional): absolute tolerance
• reltol (optional): relative tolerance

In/Outs

• x: sparse solution of the problem

Notes

This is a direct translation of Boyd et al's Matlab implementation of ADMM for basis pursuit. This method is not well suited to high-accuracy solutions, so, for our purposes, it is best used as a means of identifying the sparsity pattern.

SummationByParts.calcSparseSolution!

Finds a solution to the underdetermined problem Ax = b that is sparse. Uses basispursit! to find an approximate solution, which is used to eliminate columns from A, from which an accurate solution is found.

Inputs

• A: matrix in the linear equation that must be satisfied
• b: vector in the linear equation that must be satisfied

In/Outs

• x: sparse solution of the problem

SummationByParts.absMatrix!

Computes the absolulte value of a symmetric matrix A; that is, using the eigenvalue factorization A = EΛE^T, this function returns Aabs = E|Λ|E^T, where |Λ| is the elementwise absolute value applied to the diagonal eigenvalue matrix.

Inputs

• A: symmetric matrix whose absolute value is sought

In/Outs

• Aabs: matrix where the absolute value is stored

SummationByParts.calcMatrixEigs!

Finds the eigenvalues of the given matrix in order of increasing modulus; that is, this works for symmetric and non-symmetric square matrices. This method is basically a front end for eigfact.

Inputs

• A: matrix whose eigenvalues are desired

In/Outs

• λ: eigenvalues of A sorted in increasing modulus

SummationByParts.calcMatrixEigs_rev!

The reverse-mode differentiated version of calcMatrixEigs!. The math behind this is from Mike Giles report "An extended collection of matrix derivative results for forward and reverse mode algorithmic differentiation."

Inputs

• A: matrix whose eigenvalues are desired
• λ_bar: ∂f/∂λ vector that multiplies derivatives from left

In/Outs

• λ: eigenvalues of A sorted in increasing modulus
• A_bar: derivatives ∂f/∂A

SummationByParts.conditionObj

Let x = Znull*xred + xperp be the (unique) entries in a skew symmetric matrix S, and let E be a symmetric positive-definite matrix. This routine computes an approximation of the condition number of the matrix A = (S + E). The matrix A corresponds to a weak-form discretization of linear advection. The condition number is approximated using KS aggregation to ensure the objective is differentiable.

Inputs

• xred: a reduced-space for the entries in the skew-symmetric matrix
• p: defines the KS parameter; as p tends to infinity, we get obj = kappa(A)
• xperp: a particular solution that satisfies the SBP accuracy conditions
• Znull: matrix that defines the null-space of the SBP accuracy conditions
• E: symmetric positive-definite matrix (boundary operator for an SBP matrix)

Returns

• obj: approximate condition number of A as defined above.

SummationByParts.conditionObjGrad!

Computes the gradient of the function conditionObj with respect to xred, and returns it in the array g.

Inputs

• xred: a reduced-space for the entries in the skew-symmetric matrix
• p: defines the KS parameter; as p tends to infinity, we get obj = kappa(A)
• xperp: a particular solution that satisfies the SBP accuracy conditions
• Znull: matrix that defines the null-space of the SBP accuracy conditions
• E: symmetric positive-definite matrix (boundary operator for an SBP matrix)

In/Outs

• g: gradient of the objective conditionObj with respect to xred

SummationByParts.eigenvalueObj

Let x = Znullxred + xperp be the (unique) entries in a skew symmetric matrix S, and let E be a symmetric matrix. This routine returns the spectral radius of the matrix A = diagm(1./w)(S + |E|), where |⋅| is the elementwise absolute value. The matrix A corresponds to a strong-form discretization of linear advection.

Inputs

• xred: a reduced-space for the entries in the skew-symmetric matrix
• p: not used presently
• xperp: a particular solution that satisfies the SBP accuracy conditions
• Znull: matrix that defines the null-space of the SBP accuracy conditions
• w: diagonal norm entries in an SBP operator
• E: symmetric matrix, usually the boundary operator for an SBP matrix

Returns

• obj: the 2p-norm of the moduli of the eigenvalues of A as defined above.

SummationByParts.eigenvalueObjGrad!

Computes the gradient of the function eigenvalueObj with respect to xred, and returns it in the array g.

Inputs

• xred: a reduced-space for the entries in the skew-symmetric matrix
• p: not used at present
• xperp: a particular solution that satisfies the SBP accuracy conditions
• Znull: matrix that defines the null-space of the SBP accuracy conditions
• w: diagonal norm entries in an SBP operator
• E: symmetric matrix, usually the boundary operator for an SBP matrix

In/Outs

• g: gradient of the objective eigenvalueObj with respect to xred

SummationByParts.truncErr

This function computes the truncation error from the derivative operator.

Inputs

• d: the degree of the operator
• x: the coordinates of the quadrature points
• w: the weights of the quadrature rule
• Q: the weak differentiation matrices

In/Outs

• trunc_err: the truncation error

SummationByParts.computeConditionNumber

This function computes the condition number of A=S+E for operators on triangular elements.

Inputs

• p: the degree of the operator
• opertype: the SBP operator family type
• vtx: the vertices of the
• vertices: indicates whether vertices are included

In/Outs

• condA: the condition number of A=S+E

SummationByParts.pocssparses

This function uses the Projection Onto Convex Sets (POCS) algorithm to construct a sparse skew-symmetric S matrix.

Inputs

• S: the skew-symmetric matrix
• H: the norm/mass matrix
• E: the boundary integration matrix
• V: the Vandermonde matrix
• Vdx: the gradient of the Vandermonde matrix

In/Outs

• S: the sparsified skew-symmetric matrix

SummationByParts.quadTruncErr

This function computes the integration truncation error of the quadrature rule on the right triangle or tetrahedron.

Inputs

• cub: a symmetric cubature for the right triangle
• q: the quadarature degree

In/Outs

• quadTruncErr: the quadature truncation error

SummationByParts.checkInteriorNodeLocaton

This function checks if all interior nodes are within the interior of the right triangle or tetrahedron; i.e., not outside or on the boundaries.

Inputs

• cub: a symmetric cubature for the right triangle
• vtx: the vertices of the triangle

Note: The vertex is assumed to be on [-1,1] interval in each dimension