Mercurial > repos > public > sbplib_julia

comparison src/SbpOperators/boundaryops/normal_derivative.jl @ 1657:9c84e97de895 feature/sbp_operators/laplace_curvilinear

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix bug in normal derivative
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 27 Jun 2024 16:40:50 +0200
parents 89456aa6fa80
children cc9d18a5ff2d
comparison
equal deleted inserted replaced
1656:89456aa6fa80 1657:9c84e97de895
28 scaled_stencil = scale(closure_stencil,h_inv) 28 scaled_stencil = scale(closure_stencil,h_inv)
29 return BoundaryOperator(g, scaled_stencil, boundary) 29 return BoundaryOperator(g, scaled_stencil, boundary)
30 end 30 end
31 31
32 function normal_derivative(g::MappedGrid, stencil_set::StencilSet, boundary) 32 function normal_derivative(g::MappedGrid, stencil_set::StencilSet, boundary)
33 k = grid_id(boundary)
33 b_indices = boundary_indices(logicalgrid(g), boundary) 34 b_indices = boundary_indices(logicalgrid(g), boundary)
34 35
35 k = grid_id(boundary) 36 # Compute the weights for the logical derivatives
36
37
38 # Compute the weights for the logival derivatives
39 g⁻¹ = geometric_tensor_inverse(g) 37 g⁻¹ = geometric_tensor_inverse(g)
40 α = map(CartesianIndices(g⁻¹)[b_indices...]) do I # TODO: Fix iterator here 38 α = map(CartesianIndices(g⁻¹)[b_indices...]) do I # TODO: Fix iterator here
41 gᵏⁱ = g⁻¹[I][k,:] 39 gᵏⁱ = g⁻¹[I][k,:]
42 gᵏᵏ = g⁻¹[I][k,k] 40 gᵏᵏ = g⁻¹[I][k,k]
43 41
44 gᵏⁱ./sqrt(gᵏᵏ) 42 gᵏⁱ./sqrt(gᵏᵏ)
45 end 43 end
46 44
47 45 # Assemble difference operator
48 mapreduce(+,1:ndims(g)) do i 46 mapreduce(+,1:ndims(g)) do i
49 if i == k 47 if i == k
50 ∂ξᵢ = normal_derivative(logicalgrid(g), stencil_set, boundary) 48 ∂_ξᵢ = normal_derivative(logicalgrid(g), stencil_set, boundary)
51 else 49 else
52 ∂ξᵢ = first_derivative(logicalgrid(g), stencil_set, i) 50 e = boundary_restriction(logicalgrid(g), stencil_set, boundary)
51 ∂_ξᵢ = e ∘ first_derivative(logicalgrid(g), stencil_set, i)
53 end 52 end
54 53
55 αᵢ = componentview(α,i) 54 αᵢ = componentview(α,i)
56 DiagonalTensor(αᵢ) ∘ ∂ξᵢ 55 DiagonalTensor(αᵢ) ∘ ∂_ξᵢ
57 end 56 end
58 end 57 end