How to extract solution when using Trixi.jl

Specific Domains
1

I have trouble extracting the solution in a format of points where each point stores (x, u) where x is the coordinate (position) and u is the value at that position.

Even though I can extract the coordinates x and values u, I could not find a way to turn them into one solution. The issue is that coordinates x are duplicated (due to the stitching of nodes in Trixi meshes) and values u corresponding to the same coordinates have different value.

To my understanding one of the core features of Trixi is that it does the stitching of nodes automatically with precision guarantees.

So, how do I obtain stitched final solution with precision guarantees?

I am aware that I could obtain some solution using some interpolation, but I do not know which one should I choose and also how to implement it neatly in code (i.e. what is the recommended way of achieving the solution).

Below is the code illustrating the issue, where overlapping points have different values.

using Trixi
using OrdinaryDiffEq
using DifferentiationInterface

ode_solver = QNDF(; autodiff = AutoFiniteDiff())

N = 3
K = 2^4

L = 2.0 * π  # Length of domain
tf = 2.0 * π
a = 2.0 * π # Advection velocity
p = [a]

equations = LinearScalarAdvectionEquation1D(a)
coordinates_min = 0.0
coordinates_max = L
mesh = TreeMesh(coordinates_min, coordinates_max, initial_refinement_level=Int(log2(K)), n_cells_max=10_000)
solver = DGSEM(polydeg=N)

function initial_condition(x, t, equations::LinearScalarAdvectionEquation1D)
    u = sin(x[1])
    return SVector(u)
end

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
tspan = (0.0, tf)
ode_trixi = semidiscretize(semi, tspan)
summary_callback = SummaryCallback()
callbacks = CallbackSet(summary_callback)

sol_trixi = solve(ode_trixi, ode_solver; abstol=1e-12, reltol=1e-12, callback=callbacks)

print("u in a form of 1x4x16 array")
display(collect(Trixi.wrap_array(sol_trixi.u[end], semi)))

print("x in a form of 1x4x16 array")
display(semi.cache.elements.node_coordinates)

print("u in a form of 64 element array")
display(sol_trixi.u[end])

print("Check sol.u at overlapping coordinates x")
print("sol_trixi.u[end][4]: ", sol_trixi.u[end][4])
print("sol_trixi.u[end][5]: ", sol_trixi.u[end][5])

output:

u in a form of 1x4x16 array
1×4×16 Array{Float64, 3}:
[:, :, 1] =
 -0.978384  -0.994985  -0.997163  -0.983036

[:, :, 2] =
 -0.983132  -0.957436  -0.892337  -0.838117

[:, :, 3] =
 -0.838207  -0.774127  -0.65166  -0.565602

;;; … 

[:, :, 14] =
 -0.18315  -0.288565  -0.451422  -0.545403

[:, :, 15] =
 -0.545437  -0.632994  -0.758542  -0.82462

[:, :, 16] =
 -0.824686  -0.881056  -0.950181  -0.978296

x in a form of 1x4x16 array
1×4×16 Array{Float64, 3}:
[:, :, 1] =
 0.0  0.108539  0.28416  0.392699

[:, :, 2] =
 0.392699  0.501238  0.676859  0.785398

[:, :, 3] =
 0.785398  0.893938  1.06956  1.1781

;;; …

[:, :, 14] =
 5.10509  5.21363  5.38925  5.49779

[:, :, 15] =
 5.49779  5.60633  5.78195  5.89049

[:, :, 16] =
 5.89049  5.99903  6.17465  6.28319

u in a form of 64 element array
64-element Vector{Float64}:
 -0.9783838741096459
 -0.9949847972300191
 -0.9971633658442055
 -0.983035815133657
 -0.9831319428423163
 -0.9574364376157815
 -0.8923369221222619
 -0.8381168581916705
 -0.8382070853905123
 -0.7741270595574815
 -0.6516602710616084
 -0.5656022071382435
 -0.5656727975540656
 -0.47296385416080833
 -0.31177425104684203
 -0.20697974724482202
 -0.20701995412685925
 -0.09979618939617992
  0.07557657244916079
  0.18315350289048568
  0.1831498006757271
  ⋮
  0.5656727975541287
  0.47296385416080733
  0.31177425104686446
  0.20697974724481225
  0.20701995412687982
  0.09979618939617542
 -0.07557657244914776
 -0.18315350289062346
 -0.18314980067567863
 -0.28856454054929537
 -0.4514215478930889
 -0.5454032925014901
 -0.5454366585824441
 -0.632993935040229
 -0.7585416848168363
 -0.8246203749222623
 -0.8246857296156699
 -0.8810557410255903
 -0.9501807264247821
 -0.9782964804633923

Check sol.u at overlapping coordinates x
sol_trixi.u[end][4]: -0.983035815133657
sol_trixi.u[end][5]: -0.9831319428423163
2

You are right that Trixi.jl has two values of the solution at interfaces. This comes from the discontinuous Galerkin method, which is used in Trixi.jl.

What do you mean by “turn them into one solution” or “stiched final solution”? Do you want to obtain an interpolation, i.e., a function you can evaluate at any x to obtain a solution u?

3

Yes, obtaining an interpolation is the final goal. By “one solution”, I mean obtaining arrays of x and u where there are no duplicates in x. I know that I could obtain such “one solution” e.g. by dropping duplicates, but that is definitely not the way to achieve the best precision and lowest error.

4

Yes, probably you should not try to drop values. If you want to have an interpolation and you use DGSEM, you should use a Gauss-Lobatto-Legendre (GLL) interpolation as the inherent interpolation method of the solver. We have code in Trixi.jl, which can be helpful, but is not part of the public API. For example, we have gauss_lobatto_nodes_weights and polynomial_interpolation_matrix, which can be used to interpolate to other arbitrary points by using GLL interpolation.
Edit: This code block used in the visualization could also be helpful.