Curvilinear FDSBP #1295

andrewwinters5000 · 2022-12-12T19:00:18Z

Extends the capability of the finite difference SBP solver to work on UnstructuredMesh2D. In principle, can be extended to P4estMesh and/or StructuredMesh if desired. Basically boils down to creating new constructors for the element containers and computing the metric terms.

Central FDSBP convergence
Central FDSBP free-stream
Improve performance; use mul! when possible. Avoid allocations and the Matrix command
Upwind FDSBP convergence - done in Upwind SBP on curved meshes #1857
Upwind FDSBP free-stream - done in Upwind SBP on curved meshes #1857

ranocha

This looks already great! I pushed some changes. What's left to do in terms of performance is to rewrite the volume terms. I will look at them later.

src/solvers/fdsbp_unstructured/containers_2d.jl

ranocha · 2022-12-13T07:56:55Z

src/solvers/fdsbp_unstructured/containers_2d.jl

+# construct the metric terms for a FDSBP element "block". Directly use the derivative matrix
+# applied to the node coordinates.
+# TODO: FD performance; this is slow and allocates. Rewrite to use `mul!`
+# TODO: FD; How to make this work for the upwind solver because basis has three available derivative matrices


This is a very good question that will probably require some analytical investigations. A first attempt may be to use the central operator associated to upwind SBP operators - you should be able to access it as basis.central.

I am thinking about a helper function like get_derivative_matrix(basis) that would return either D_SBP or basis.central. I want to try this soon

src/solvers/fdsbp_unstructured/containers_2d.jl

src/solvers/fdsbp_unstructured/fdsbp_2d.jl

ranocha · 2022-12-13T09:08:20Z

src/solvers/fdsbp_unstructured/fdsbp_2d.jl

+function calc_error_norms(func, u, t, analyzer,
+                          mesh::UnstructuredMesh2D, equations, initial_condition,
+                          dg::FDSBP, cache, cache_analysis)
+  # TODO: FD. This is rather inefficient right now and allocates...


I added this already in the TreeMesh version and hope to fix it with a new interface in SummationByPartsOperators.jl, see ranocha/SummationByPartsOperators.jl#161.

examples/unstructured_2d_fdsbp/elixir_euler_free_stream.jl

Co-authored-by: Hendrik Ranocha <[email protected]>

src/solvers/fdsbp_unstructured/fdsbp_2d.jl

ranocha · 2022-12-14T03:16:33Z

src/solvers/fdsbp_unstructured/fdsbp_2d.jl

+    # x direction
+    for j in eachnode(dg)
+      for i in eachnode(dg)
+        Ja1 = get_contravariant_vector(1, contravariant_vectors, i, j, element)
+        f_element[i, j] = flux(u_element[i, j], Ja1, equations)
+      end
+      mul!(view(du_vectors, :, j, element), D, view(f_element, :, j),
+           one(eltype(du)), one(eltype(du)))
+    end
+
+    # y direction
+    for i in eachnode(dg)
+      for j in eachnode(dg)
+        Ja2 = get_contravariant_vector(2, contravariant_vectors, i, j, element)
+        f_element[i, j] = flux(u_element[i, j], Ja2, equations)
+      end
+      mul!(view(du_vectors, i, :, element), D, view(f_element, i, :),
+           one(eltype(du)), one(eltype(du)))
+    end


This is not the de-aliased form, so it will not be stable for linear advection, will it?

That is true. This is the standard form, so it is not provably stable for variable coefficients from the metric terms. This should be easy to put in though.

We should think about an interface. Right now, it is still experimental, so we are free to do what we think is best.

I think we should be somewhat consistent with the VolumeIntegralWeakForm - which is also using the plain form without any guarantees. We could add a specialization (at first to the VolumeIntegralStrongForm, later also to the VolumeIntegralWeakForm) that controls whether we use the plain form or the "correct" de-aliased form. What do you think?

That is a good idea. It would be equivalent to a particular version of the flux differencing volume integral. After all, the split forms were originally an FD concept so it should be possible to get all that functionality working here. It is mostly an issue of making it fast because the DG version does everything loop based due to dense D matrices, whereas the FD should exploit sparsity.

Yeah, that is definitely a bigger issue - requiring new interfaces in SummationByPartsOperators.jl or the conversion to sparse matrices Jesse does for DGMulti.
I was thinking about the special case of strong/weak form volume terms with de-aliasing, which should be basically equivalent to flux differencing with the central flux but can hopefully be implemented more easily using mul!.

I added this to the list in #1284

…nto curvi-fdsbp

ranocha · 2022-12-20T07:54:15Z

@andrewwinters5000 What about cleaning up this PR containing only changes for the classical central SBP operators and postponing the upwind stuff to another PR later?

… duplicate code

…th biased metric terms

ranocha · 2023-01-17T05:55:54Z

src/equations/compressible_euler_2d.jl

@@ -970,6 +970,46 @@ end
  return SVector(f1m, f2m, f3m, f4m)
 end

+@inline function splitting_lax_friedrichs(u, ::Val{:plus}, normal_direction::AbstractVector,


We need to remember to update the docstring of splitting_lax_friedrichs

ranocha · 2023-01-17T06:07:22Z

src/equations/compressible_euler_2d.jl

+                                          equations::CompressibleEulerEquations2D)
+  rho, rho_v1, rho_v2, rho_e = u
+  v1 = rho_v1 / rho
+  v2 = rho_v2 / rho
+  v_dot_n = normal_direction[1] * v1 + normal_direction[2] * v2
+  p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2))
+
+  a = sqrt(equations.gamma * p / rho)
+  H = (rho_e + p) / rho
+  lambda = 0.5 * (v_dot_n + a * norm(normal_direction))
+
+  f1p = 0.5 * rho * v_dot_n + lambda * u[1]
+  f2p = 0.5 * rho * v_dot_n * v1 + 0.5 * p * normal_direction[1] + lambda * u[2]
+  f3p = 0.5 * rho * v_dot_n * v2 + 0.5 * p * normal_direction[2] + lambda * u[3]
+  f4p = 0.5 * rho * v_dot_n * H + lambda * u[4]
+
+  return SVector(f1p, f2p, f3p, f4p)
+end


This is not the same as rotating to x, applying the Cartesian splitting, and rotating back:

julia> using Revise; using Trixi julia> equations = CompressibleEulerEquations2D(1.4) ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ │ CompressibleEulerEquations2D │ │ ════════════════════════════ │ │ #variables: ………………………………………………… 4 │ │ │ variable 1: …………………………………………… rho │ │ │ variable 2: …………………………………………… rho_v1 │ │ │ variable 3: …………………………………………… rho_v2 │ │ │ variable 4: …………………………………………… rho_e │ └──────────────────────────────────────────────────────────────────────────────────────────────────┘ julia> u = SVector(1.0, 0.2, -0.3, 2.0) 4-element SVector{4, Float64} with indices SOneTo(4): 1.0 0.2 -0.3 2.0 julia> normal_direction = SVector(0.3, -0.5) 2-element SVector{2, Float64} with indices SOneTo(2): 0.3 -0.5 julia> let norm_ = Trixi.norm(normal_direction) normal_vector = normal_direction / norm_ u_rotated = Trixi.rotate_to_x(u, normal_vector, equations) f = splitting_lax_friedrichs(u_rotated, Val{:plus}(), 1, equations) f_rotated = Trixi.rotate_from_x(f, normal_vector, equations) * norm_ f_rotated - splitting_lax_friedrichs(u, Val{:plus}(), normal_direction, equations) end 4-element SVector{4, Float64} with indices SOneTo(4): 0.00011898020814316013 2.3796041628665332e-5 -3.5694062443025754e-5 0.0002379604162865423

You mentioned some unusual behavior required for the splittings in general. Is this behavior what we want/need?

ranocha · 2023-01-17T06:10:51Z

src/equations/numerical_fluxes.jl

+# TODO: FD. Works properly for `splitting_lax_friedrichs` but needs further
+#       testing for other splittings.
+@inline function (numflux::FluxUpwind)(u_ll, u_rr, normal_direction::AbstractVector,
+                                       equations::AbstractEquations{2})
+  @unpack splitting = numflux
+  # Compute splitting in generic normal direction with specialized
+  # eigenvalues estimates calculated inside the `splitting` function
+  f_tilde_m = splitting(u_rr, Val{:minus}(), normal_direction, equations)
+  f_tilde_p = splitting(u_ll, Val{:plus}() , normal_direction, equations)
+  return f_tilde_m + f_tilde_p
+end


Continuing the comment above, this is not consistent with the approach to compute numerical fluxes by rotation:

julia> using Revise; using Trixi julia> equations = CompressibleEulerEquations2D(1.4) ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ │ CompressibleEulerEquations2D │ │ ════════════════════════════ │ │ #variables: ………………………………………………… 4 │ │ │ variable 1: …………………………………………… rho │ │ │ variable 2: …………………………………………… rho_v1 │ │ │ variable 3: …………………………………………… rho_v2 │ │ │ variable 4: …………………………………………… rho_e │ └──────────────────────────────────────────────────────────────────────────────────────────────────┘ julia> u_ll = SVector(1.0, 0.2, -0.3, 2.0) 4-element SVector{4, Float64} with indices SOneTo(4): 1.0 0.2 -0.3 2.0 julia> u_rr = SVector(1.1, 0.3, -0.2, 2.1) 4-element SVector{4, Float64} with indices SOneTo(4): 1.1 0.3 -0.2 2.1 julia> normal_direction = SVector(0.3, -0.5) 2-element SVector{2, Float64} with indices SOneTo(2): 0.3 -0.5 julia> numflux = FluxUpwind(splitting_lax_friedrichs) FluxUpwind(splitting_lax_friedrichs) julia> numflux(u_ll, u_rr, normal_direction, equations) - FluxRotated(numflux)(u_ll, u_rr, normal_direction, equations) 4-element SVector{4, Float64} with indices SOneTo(4): 0.00999999999999987 0.002735925833319497 -0.001804120520855823 0.019080092708350427

This appears to be strange, doesn't it?

ranocha · 2023-01-17T06:14:14Z

src/solvers/fdsbp_unstructured/containers_2d.jl

+# function init_element!(elements, element, basis::AbstractDerivativeOperator, corners_or_surface_curves)
+function init_element!(elements, element, basis::DerivativeOperator, corners_or_surface_curves)


Why didn't you choose the version with basis::AbstractDerivativeOperator?

andrewwinters5000 · 2023-01-19T09:19:38Z

test/test_trixi.jl

+          # Ignore the message about ODE Symbols nonsense. DO NOT COMMIT! This is local to my machine
+          "┌ Warning: Backwards compatability support of the new return codes to Symbols will be deprecated with the Julia v1.9 release. Please see https://docs.sciml.ai/SciMLBase/stable/interfaces/Solutions/#retcodes for more information
+└ @ SciMLBase ~/.julia/packages/SciMLBase/0IYdl/src/retcodes.jl:360\n"


This is a reminder to myself to remove this extra ignore information. I accidentally committed it.

andrewwinters5000 and others added 5 commits December 12, 2022 19:50

add FDSBP for curvilinear meshes

77cd96e

new elixirs for free-stream and convergence test

00b0eff

add tests

6229c85

avoid allocations in init_element!

ef47e92

avoid allocations in calc_surface_integral!

f512d03

ranocha requested changes Dec 13, 2022

View reviewed changes

andrewwinters5000 and others added 6 commits December 13, 2022 12:01

Apply suggestions from code review

5f968f7

Co-authored-by: Hendrik Ranocha <[email protected]>

add helper function to get the derivative matrix from basis

3cf1886

increased tolerance for free stream test

14b1d96

WIP: improve performance of volume terms

46a9866

remove allocations in volume terms

e3b2772

improve performance of volume terms

6e6aedc

ranocha reviewed Dec 13, 2022

View reviewed changes

src/solvers/fdsbp_unstructured/fdsbp_2d.jl Show resolved Hide resolved

ranocha reviewed Dec 14, 2022

View reviewed changes

ranocha and others added 3 commits December 14, 2022 08:46

Merge branch 'main' into curvi-fdsbp

674bddd

Merge branch 'curvi-fdsbp' of github.com:andrewwinters5000/Trixi.jl i…

89e97b3

…nto curvi-fdsbp

fix summary output when discretization uses VolumeIntegralStrongForm

c2d1877

ranocha mentioned this pull request Dec 20, 2022

Finish upwind SBP implementation #1284

Open

8 tasks

andrewwinters5000 added 5 commits December 29, 2022 19:36

unify UnstructuredElement2D constructors for DGSEM and FDSBP. Removes…

101aa57

… duplicate code

Merge branch 'main' into curvi-fdsbp

9230503

new container type for biased metric terms. Upwind volume integral wi…

3c2b1e6

…th biased metric terms

add two tests for the upwind curvilinear solver

588453e

Merge branch 'main' into curvi-fdsbp

c545a6b

ranocha reviewed Jan 17, 2023

View reviewed changes

andrewwinters5000 commented Jan 19, 2023

View reviewed changes

andrewwinters5000 mentioned this pull request Dec 13, 2023

Central SBP finite difference solver for UnstructuredMesh2D #1773

Merged

andrewwinters5000 mentioned this pull request Mar 4, 2024

Upwind SBP on curved meshes #1857

Merged

andrewwinters5000 closed this Mar 6, 2024

andrewwinters5000 deleted the curvi-fdsbp branch March 6, 2024 10:31

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Curvilinear FDSBP #1295

Curvilinear FDSBP #1295

andrewwinters5000 commented Dec 12, 2022 •

edited

Loading

ranocha left a comment

ranocha Dec 13, 2022

andrewwinters5000 Dec 13, 2022

ranocha Dec 13, 2022

ranocha Dec 14, 2022

andrewwinters5000 Dec 14, 2022

ranocha Dec 14, 2022 •

edited

Loading

andrewwinters5000 Dec 14, 2022

ranocha Dec 14, 2022

ranocha Dec 20, 2022

ranocha commented Dec 20, 2022

ranocha Jan 17, 2023

ranocha Jan 17, 2023

ranocha Jan 17, 2023

ranocha Jan 17, 2023

andrewwinters5000 Jan 19, 2023

		# function init_element!(elements, element, basis::AbstractDerivativeOperator, corners_or_surface_curves)
		function init_element!(elements, element, basis::DerivativeOperator, corners_or_surface_curves)

Curvilinear FDSBP #1295

Curvilinear FDSBP #1295

Conversation

andrewwinters5000 commented Dec 12, 2022 • edited Loading

ranocha left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ranocha Dec 14, 2022 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ranocha commented Dec 20, 2022

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

andrewwinters5000 commented Dec 12, 2022 •

edited

Loading

ranocha Dec 14, 2022 •

edited

Loading