Second order finite volume integral in 2D

[sbairos]

I mostly just copied what you had in trixi-framework#2022 but made it run in 2D.

Convergence seems decent.

l2
rho                 rho_v1              rho_v2              rho_e               
error     EOC       error     EOC       error     EOC       error     EOC       
6.56e-04  -         4.27e-04  -         4.27e-04  -         1.32e-03  -         
1.60e-04  2.04      1.02e-04  2.06      1.02e-04  2.06      3.17e-04  2.06      
3.91e-05  2.03      2.50e-05  2.03      2.50e-05  2.03      7.79e-05  2.03      

mean      2.03      mean      2.05      mean      2.05      mean      2.04      
----------------------------------------------------------------------------------------------------
linf
rho                 rho_v1              rho_v2              rho_e               
error     EOC       error     EOC       error     EOC       error     EOC       
1.93e-03  -         1.39e-03  -         1.39e-03  -         3.94e-03  -         
5.85e-04  1.72      4.20e-04  1.73      4.20e-04  1.73      1.07e-03  1.88      
1.62e-04  1.86      1.18e-04  1.83      1.18e-04  1.83      2.87e-04  1.90      

mean      1.79      mean      1.78      mean      1.78      mean      1.89      
----------------------------------------------------------------------------------------------------
Dict{Symbol, Any} with 3 entries:
  :variables => ("rho", "rho_v1", "rho_v2", "rho_e")
  :l2        => [2.03433, 2.04642, 2.04642, 2.04116]
  :linf      => [1.78902, 1.78095, 1.78095, 1.8899]

[github-actions]

Review checklist

This checklist is meant to assist creators of PRs (to let them know what reviewers will typically look for) and reviewers (to guide them in a structured review process). Items do not need to be checked explicitly for a PR to be eligible for merging.

Purpose and scope

• The PR has a single goal that is clear from the PR title and/or description.
• All code changes represent a single set of modifications that logically belong together.
• No more than 500 lines of code are changed or there is no obvious way to split the PR into multiple PRs.

Code quality

• The code can be understood easily.
• Newly introduced names for variables etc. are self-descriptive and consistent with existing naming conventions.
• There are no redundancies that can be removed by simple modularization/refactoring.
• There are no leftover debug statements or commented code sections.
• The code adheres to our conventions and style guide, and to the Julia guidelines.

Documentation

• New functions and types are documented with a docstring or top-level comment.
• Relevant publications are referenced in docstrings (see example for formatting).
• Inline comments are used to document longer or unusual code sections.
• Comments describe intent ("why?") and not just functionality ("what?").
• If the PR introduces a significant change or new feature, it is documented in NEWS.md with its PR number.

Testing

• The PR passes all tests.
• New or modified lines of code are covered by tests.
• New or modified tests run in less then 10 seconds.

Performance

• There are no type instabilities or memory allocations in performance-critical parts.
• If the PR intent is to improve performance, before/after time measurements are posted in the PR.

Verification

• The correctness of the code was verified using appropriate tests.
• If new equations/methods are added, a convergence test has been run and the results
  are posted in the PR.

Created with ❤️ by the Trixi.jl community.

[DanielDoehring]

Could you add one more iteration of the convergence test? If that yields the expected order you can open this PR at the main repo

[sbairos]

So I actually wasn't able to do an additional level of convergence_test because I was running out of RAM. I ended up decreasing the initial refinement level so that I can run 4 levels of the convergence test without OOMing. Here's the new results:

####################################################################################################
l2
rho                 rho_v1              rho_v2              rho_e               
error     EOC       error     EOC       error     EOC       error     EOC       
2.50e-03  -         1.64e-03  -         1.64e-03  -         5.33e-03  -         
6.56e-04  1.93      4.27e-04  1.95      4.27e-04  1.95      1.32e-03  2.02      
1.60e-04  2.04      1.02e-04  2.06      1.02e-04  2.06      3.17e-04  2.06      
3.91e-05  2.03      2.50e-05  2.03      2.50e-05  2.03      7.79e-05  2.03      

mean      2.00      mean      2.01      mean      2.01      mean      2.03      
----------------------------------------------------------------------------------------------------
linf
rho                 rho_v1              rho_v2              rho_e               
error     EOC       error     EOC       error     EOC       error     EOC       
5.80e-03  -         4.56e-03  -         4.56e-03  -         1.49e-02  -         
1.93e-03  1.58      1.39e-03  1.71      1.39e-03  1.71      3.94e-03  1.92      
5.85e-04  1.72      4.20e-04  1.73      4.20e-04  1.73      1.07e-03  1.88      
1.62e-04  1.86      1.18e-04  1.83      1.18e-04  1.83      2.87e-04  1.90      

mean      1.72      mean      1.76      mean      1.76      mean      1.90      
----------------------------------------------------------------------------------------------------
Dict{Symbol, Any} with 3 entries:
  :variables => ("rho", "rho_v1", "rho_v2", "rho_e")
  :l2        => [2.0002, 2.0131, 2.0131, 2.03272]
  :linf      => [1.72094, 1.75754, 1.75754, 1.90058]

Without this change, I was seeing:

ERROR: LoadError: AssertionError: New length would exceed capacity
Stacktrace:
 [1] resize!
   @ ~/trixi/Trixi.jl/src/auxiliary/containers.jl:51 [inlined]
 [2] refine_uniformly!(t::Trixi.SerialTree{2, Float64}, max_level::Int64)
   @ Trixi ~/trixi/Trixi.jl/src/meshes/abstract_tree.jl:212
 [3] macro expansion
   @ ~/.julia/packages/TrixiBase/38UFo/src/trixi_timeit.jl:64 [inlined]
 [4] initialize!(mesh::TreeMesh{…}, initial_refinement_level::Int64, refinement_patches::Tuple{}, coarsening_patches::Tuple{})
   @ Trixi ~/trixi/Trixi.jl/src/meshes/tree_mesh.jl:173
 [5] TreeMesh(coordinates_min::Tuple{…}, coordinates_max::Tuple{…}; n_cells_max::Int64, periodicity::Bool, initial_refinement_level::Int64, refinement_patches::Tuple{}, coarsening_patches::Tuple{}, RealT::Type)
   @ Trixi ~/trixi/Trixi.jl/src/meshes/tree_mesh.jl:165
 [6] top-level scope
   @ ~/trixi/Trixi.jl/examples/tree_2d_dgsem/elixir_euler_convergence_pure_fvO2.jl:24
in expression starting at /home/purple/trixi/Trixi.jl/examples/tree_2d_dgsem/elixir_euler_convergence_pure_fvO2.jl:24
Some type information was truncated. Use `show(err)` to see complete types.

[DanielDoehring]

Oh okay, yeah looks good!

Can you also check one of the StructuredMesh testcases with a non-trivial mesh if this works as well?

[sbairos]

Ahh good idea with trying StructuredMesh. There seems to be something wrong for them, I seem to consistently end up with negative values for pressure and density in an Euler example, which crashes the simulation after the first timestep. Will need to investigate further

[DanielDoehring]

okay, crashing in the first timestep indeed seems to hint that something is wrong. You might need to check if there are anny differences between the already implemented plain first-order finite volume integral between 1D and 2D

[sbairos]

Alright, finally figured it out and you were correct. Turns out that in 2D, StructuredMesh does not share calcflux_fv() with TreeMesh since StructuredMesh has to deal with cell orientations. Can't believe how much debugging it took to finally arrive at that 😅

Anyways, here's the convergence_test for StructuredMesh (which looks a bit suspiciously good in my opinion):

####################################################################################################
l2
rho                 rho_v1              rho_v2              rho_e               
error     EOC       error     EOC       error     EOC       error     EOC       
2.75e-03  -         1.78e-03  -         1.78e-03  -         5.59e-03  -         
7.22e-04  1.93      4.33e-04  2.04      4.33e-04  2.04      1.34e-03  2.06      
1.78e-04  2.02      1.04e-04  2.07      1.04e-04  2.07      3.23e-04  2.06      
4.41e-05  2.02      2.53e-05  2.03      2.53e-05  2.03      7.92e-05  2.03      

mean      1.99      mean      2.05      mean      2.05      mean      2.05      
----------------------------------------------------------------------------------------------------
linf
rho                 rho_v1              rho_v2              rho_e               
error     EOC       error     EOC       error     EOC       error     EOC       
1.08e-02  -         7.87e-03  -         7.87e-03  -         2.42e-02  -         
3.19e-03  1.76      2.28e-03  1.79      2.28e-03  1.79      7.39e-03  1.71      
9.69e-04  1.72      5.80e-04  1.98      5.80e-04  1.98      1.93e-03  1.93      
3.39e-04  1.52      1.81e-04  1.68      1.81e-04  1.68      5.01e-04  1.95      

mean      1.66      mean      1.81      mean      1.81      mean      1.86      
----------------------------------------------------------------------------------------------------

[DanielDoehring]

Good job!

I have a couple minor requests, but overall this looks quite good!

[DanielDoehring]

You can also take a look at trixi-framework#2678 where I revisited some comments to the FV subcell method.

[sbairos]

I had to do reconstruction_O2_full instead of reconstruction_O2_inner here because inner was crashing after around 70 iterations with:

  LoadError: DomainError with -1.7433242556262465:
  sqrt was called with a negative real argument but will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).
  DomainError detected in the user `f` function. This occurs when the domain of a function is violated.
  For example, `log(-1.0)` is undefined because `log` of a real number is defined to only output real
  numbers, but `log` of a negative number is complex valued and therefore Julia throws a DomainError
  by default. Cases to be aware of include:
  
  * `log(x)`, `sqrt(x)`, `cbrt(x)`, etc. where `x<0`
  * `x^y` for `x<0` floating point `y` (example: `(-1.0)^(1/2) == im`)
  
  Within the context of SciML, this error can occur within the solver process even if the domain constraint
  would not be violated in the solution due to adaptivity. For example, an ODE solver or optimization
  routine may check a step at `new_u` which violates the domain constraint, and if violated reject the
  step and use a smaller `dt`. However, the throwing of this error will have halted the solving process.
  
  Thus the recommended fix is to replace this function with the equivalent ones from NaNMath.jl
  (https://github.com/JuliaMath/NaNMath.jl) which returns a NaN instead of an error. The solver will then
  effectively use the NaN within the error control routines to reject the out of bounds step. Additionally,
  one could perform a domain transformation on the variables so that such an issue does not occur in the
  definition of `f`.
  
  For more information, check out the following FAQ page:
  https://docs.sciml.ai/Optimization/stable/API/FAQ/#The-Solver-Seems-to-Violate-Constraints-During-the-Optimization,-Causing-DomainErrors,-What-Can-I-Do-About-That?
  
  Stacktrace:
    [1] throw_complex_domainerror(f::Symbol, x::Float64)
      @ Base.Math ./math.jl:33
    [2] sqrt
      @ ./math.jl:686 [inlined]
    [3] flux_hllc(u_ll::SVector{4, Float64}, u_rr::SVector{4, Float64}, normal_direction::SVector{2, Float64}, equations::CompressibleEulerEquations2D{Float64})
      @ Trixi ~/trixi/Trixi.jl/src/equations/compressible_euler_2d.jl:1810
    [4] calcflux_fvO2!
      @ ~/trixi/Trixi.jl/src/solvers/dgsem_structured/dg_2d.jl:325 [inlined]
    [5] fvO2_kernel!
      @ ~/trixi/Trixi.jl/src/solvers/dgsem_tree/dg_2d.jl:314 [inlined]
    [6] macro expansion
      @ ~/trixi/Trixi.jl/src/solvers/dgsem_tree/dg_2d.jl:341 [inlined]
    [7] macro expansion
      @ ~/.julia/packages/Polyester/almvr/src/closure.jl:525 [inlined]
    [8] macro expansion
      @ ~/trixi/Trixi.jl/src/auxiliary/auxiliary.jl:213 [inlined]
(stack trace goes on but is full of excessively huge calls after this that won't be displayed here well)

[DanielDoehring]

Ok, if that works fine - have you tried reducing the CFL? Since the spatial discretization changes you might need to adapt CFL as well.

[DanielDoehring]

Looks good! You can open the PR at the main repo soon!