diff src/SbpOperators/volumeops/derivatives/dissipation.jl @ 1355:102ebdaf7c11 feature/variable_derivatives

Merge default
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 08 Feb 2023 21:21:28 +0100
parents d60a10ad6579
children 356ec6a72974
line wrap: on
line diff
--- a/src/SbpOperators/volumeops/derivatives/dissipation.jl	Fri Feb 03 22:50:42 2023 +0100
+++ b/src/SbpOperators/volumeops/derivatives/dissipation.jl	Wed Feb 08 21:21:28 2023 +0100
@@ -1,4 +1,16 @@
-function undivided_dissipation(g::EquidistantGrid, p, direction)
+"""
+    undivided_skewed04(g::EquidistantGrid, p, direction)
+
+Create undivided difference operators approximating the `p`th derivative. The
+operators do not satisfy any SBP-property and are meant to be used for
+building artificial dissipation terms.
+
+The operators and how they are used to create accurate artifical dissipation
+is described in "K. Mattsson, M. Svärd, and J. Nordström, “Stable and Accurate
+Artificial Dissipation,” Journal of Scientific Computing, vol. 21, no. 1, pp.
+57–79, Aug. 2004"
+"""
+function undivided_skewed04(g::EquidistantGrid, p, direction)
     T = eltype(g)
     interior_weights = T.(dissipation_interior_weights(p))
 
@@ -20,7 +32,7 @@
     return D, Dᵀ
 end
 
-undivided_dissipation(g::EquidistantGrid{1}, p) = undivided_dissipation(g, p, 1)
+undivided_skewed04(g::EquidistantGrid{1}, p) = undivided_skewed04(g, p, 1)
 
 function dissipation_interior_weights(p)
    if p == 0
@@ -47,8 +59,16 @@
 dissipation_lower_closure_size(weights) = midpoint(weights) - 1
 dissipation_upper_closure_size(weights) = length(weights) - midpoint(weights)
 
-dissipation_lower_closure_stencils(interior_weights) = ntuple(i->Stencil(interior_weights..., center=i                       ), dissipation_lower_closure_size(interior_weights))
-dissipation_upper_closure_stencils(interior_weights) = ntuple(i->Stencil(interior_weights..., center=length(interior_weights)-dissipation_upper_closure_size(interior_weights)+i), dissipation_upper_closure_size(interior_weights))
+function dissipation_lower_closure_stencils(interior_weights)
+    stencil(i) = Stencil(interior_weights..., center=i)
+    return ntuple(i->stencil(i), dissipation_lower_closure_size(interior_weights))
+end
+
+function dissipation_upper_closure_stencils(interior_weights)
+    center(i) = length(interior_weights) - dissipation_upper_closure_size(interior_weights) + i
+    stencil(i) = Stencil(interior_weights..., center=center(i))
+    return ntuple(i->stencil(i), dissipation_upper_closure_size(interior_weights))
+end
 
 function dissipation_transpose_lower_closure_stencils(interior_weights)
     closure = ntuple(i->dissipation_transpose_lower_closure_stencil(interior_weights, i), length(interior_weights))