Skip to content

4. Using SBP operators

Jump to bottom Edit New page
Jason Hicken edited this page May 2, 2015 · 14 revisions

The following are defined for both 2D and 3D elements.

calcnodes

This function returns the node coordinates for an SBP operator. The nodes are ordered first by vertex, then by edge (with nodes ordered in sequence along the directed edge), then by face (if appropriate), and then finally by volumn nodes. 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 (see also test/test_useoperators.jl)

x[:,:,1] = calcnodes(sbp, vtx)

weakdifferentiate!

Applies the SBP stiffness matrix to data in u and adds the result to 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
  • di: direction index of the operator that is desired (di=1 for Qx, etc)
  • u: the array that the operator is applied to

In/Outs

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

Example (see also test/test_useoperators.jl)

weakdifferentiate!(sbp, di, u, res)

differentiate!

Applies the SBP differentiation matrix operator, D, to data in u and adds the result to 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
  • di: direction index of the operator that is desired (di=1 for Dx, etc)
  • u: the array that the operator is applied to

In/Outs

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

Example (see also test/test_useoperators.jl)

differentiate!(sbp, di, u, res)

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

In/Outs

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

Example (see also test/test_useoperators.jl)

volumeintegrate!(sbp, u, res)

mappingjacobian!

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 dxidx 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

  • dxidx: 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

Example (see also test/test_useoperators.jl)

mappingjacobian!(sbp, x, dxidx, jac)

boundaryintegrate!

Integrates a numerical flux over a boundary using appropriate mass matrices defined on the element faces. 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
  • bndryfaces: list of boundary faces stored as an array of Boundarys
  • u: the array of data that the boundary integral depends on
  • dξdx: Jacobian of the element mapping (as output from mappingjacobian!)
  • bndryflux: function to compute the numerical flux over the boundary

In/Outs

  • res: where the result of the integration is stored

Example (see also test/test_useoperators.jl)

# 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)
# compute the Jacobian of the mapping
dξdx = zeros(Float64, (2,2,sbp.numnodes,2))
jac = zeros(Float64, (sbp.numnodes,2))
mappingjacobian!(sbp, x, dξdx, jac)
# identify the boundary faces
bndryfaces = Array(Boundary, 4)
bndryfaces[1] = Boundary(1,1) 
bndryfaces[2] = Boundary(1,3) # face 3 of element 1
bndryfaces[3] = Boundary(2,1)
bndryfaces[4] = Boundary(2,2)
# define a field to integrate, and a storage location
u = rand(Float64, (sbp.numnodes,4))
res = zeros(u)
# define a numerical flux function
function bndryflux(u, dξdx, nrm)
  return u*sum(nrm.'*dξdx)
end
# integrate function bndryflux over the boundary
boundaryintegrate!(sbp, bndryfaces, u, dξdx, bndryflux, res)

SummationByParts.edgestabilize!

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

  • For scalar fields, it is assumed that u is a rank-2 array, with the first
    dimension for the local-node index, and the second dimension for the element
    index.
  • For vector fields, u is a rank-3 array, with the first dimension for the
    index of the vector field, the second dimension for the local-node index, and
    the third dimension for the element index.

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

Inputs

  • sbp: an SBP operator type
  • ifaces: list of element interfaces stored as an array of Interfaces
  • u: the array of solution data
  • x: Cartesian coordinates stored in (coord,node,element) format
  • dξdx: scaled Jacobian of the mapping (as output from mappingjacobian!)
  • jac: determinant of the Jacobian
  • α: array of transformation terms (see below)
  • stabscale: function to compute the edge-stabilization scaling (see below)

In/Outs

  • res: where the result of the integration is stored

Details

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

for k = 1:mesh.numelem                                                                                
  for i = 1:sbp.numnodes                                                                              
    for di1 = 1:2                                                                                     
      for di2 = 1:2                                                                                   
        α[di1,di2,i,k] = (dξdx[di1,1,i,k].*dξdx[di2,1,i,k] +                                          
                          dξdx[di1,2,i,k].*dξdx[di2,2,i,k])*jac[i,k]                                  
      end                                                                                             
    end                                                                                               
  end                                                                                                 
end

The function stabscale has the signature

function stabscale(u, dξdx, nrm)

where u is the solution at a node, dξdx is the (scaled) Jacobian at the same
node, and nrm is the normal to the edge in reference space. stabscale should
return the necessary scaling to ensure the edge-stabilization has the desired
asymptotic convergence rate.

Add a custom footer
Add a custom sidebar
Clone this wiki locally