Tree Mesh 2D Dirichlet

[maxbertrand1996]

src/solvers/dgsem_tree/dg_2d.jl did not yet consider non periodic boundary conditions, if the underlying equation has non conservative terms.
This has been added using a similar logic as in src/solvers/dgsem_tree/dg_1d.jl

[maxbertrand1996]

Yes, it only calls on the surface flux. I thought this is intended, as the inital slip wall consition also only calls on the surface flux?
But it does call on the non conxervative surface flux term. (Or at least tries to.)

I also dir a simulation with a discontinuous bottom topography.

This is the underlying code

using LinearAlgebra
using Plots
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the shallow water equations with a discontinuous
# bottom topography function

equations = ShallowWaterEquations2D(gravity_constant=1.0, H0=7.0)

# An initial condition with constant total water height and zero velocities to test well-balancedness.
# Note, this routine is used to compute errors in the analysis callback but the initialization is
# overwritten by `initial_condition_discontinuous_well_balancedness` below.
function initial_condition_well_balancedness(x, t, equations::ShallowWaterEquations2D)

  inicenter = SVector(0.0, 0.0)
  x_norm = x - inicenter
  r = norm(x_norm)

  # Calculate primitive variables
  H =  r < 0.1 ? 9.0 : 7.0
  
  v1 = 0.0
  v2 = 0.0
  # bottom topography taken from Pond.control in [HOHQMesh](https://github.com/trixi-framework/HOHQMesh)
  x1, x2 = x
  b = (  1.5 / exp( 0.5 * ((x1 - 1.0)^2 + (x2 - 1.0)^2) )
       + 0.75 / exp( 0.5 * ((x1 + 1.0)^2 + (x2 + 1.0)^2) ) )
  return prim2cons(SVector(H, v1, v2, b), equations)
end

initial_condition = initial_condition_well_balancedness
boundary_condition = boundary_condition_slip_wall

###############################################################################
# Get the DG approximation space

volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
surface_flux = (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal)
solver = DGSEM(polydeg=4, surface_flux=surface_flux,
               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

###############################################################################
# Get the TreeMesh and setup a periodic mesh

coordinates_min = (-1.0, -1.0)
coordinates_max = ( 1.0,  1.0)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=2,
                n_cells_max=10_000,
                periodicity = false)

# Create the semi discretization object
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver;
                                    boundary_conditions = boundary_condition)

###############################################################################
# ODE solver

tspan = (0.0, 1.0)
ode = semidiscretize(semi, tspan)

###############################################################################
# Workaround to set a discontinuous bottom topography and initial condition for debugging and testing.

# alternative version of the initial conditinon used to setup a truly discontinuous
# bottom topography function for this academic testcase of well-balancedness.
# The errors from the analysis callback are not important but the error for this lake at rest test case
# `∑|H0-(h+b)|` should be around machine roundoff
# In contrast to the usual signature of initial conditions, this one get passed the
# `element_id` explicitly. In particular, this initial conditions works as intended
# only for the TreeMesh2D with initial_refinement_level=2.
function initial_condition_discontinuous_well_balancedness(x, t, element_id, equations::ShallowWaterEquations2D)

  inicenter = SVector(0.0, 0.0)
  x_norm = x - inicenter
  r = norm(x_norm)

  # Calculate primitive variables
  H =  r < 0.1 ? 9.0 : 7.0
  
  v1 = 0.0
  v2 = 0.0
  x1, x2 = x
  b = (  1.5 / exp( 0.5 * ((x1 - 1.0)^2 + (x2 - 1.0)^2) )
  + 0.75 / exp( 0.5 * ((x1 + 1.0)^2 + (x2 + 1.0)^2) ) )

  # Setup a discontinuous bottom topography using the element id number
  if element_id == 7
    b = 2.0 + 0.5 * sin(2.0 * pi * x[1]) + 0.5 * cos(2.0 * pi * x[2])
  end

  return prim2cons(SVector(H, v1, v2, b), equations)
end

# point to the data we want to augment
u = Trixi.wrap_array(ode.u0, semi)
# reset the initial condition
for element in eachelement(semi.solver, semi.cache)
  for j in eachnode(semi.solver), i in eachnode(semi.solver)
    x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations, semi.solver, i, j, element)
    u_node = initial_condition_discontinuous_well_balancedness(x_node, first(tspan), element, equations)
    Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, j, element)
  end
end

###############################################################################
# run the simulation

sol = solve(ode, RDPK3SpFSAL49(), abstol=1.0e-8, reltol=1.0e-8, save_everystep=true);

pyplot()

animation = @animate for k= 1:5:700
  pd = PlotData2D(sol.u[k], semi)
  wireframe(pd["H"])
  surface!(pd["b"], zlim=(0,10), camera = (30,20), title="t=$(k-1)")
end

gif(animation, "dam_break_2d.gif", fps=10)

[andrewwinters5000]

Yes, it only calls on the surface flux. I thought this is intended, as the inital slip wall consition also only calls on the surface flux?

Yeah, this is consistent with how the wall BC is set for the other mesh types with the shallow water equations. The last example run you provided looks good.

[andrewwinters5000]

The implementation of the wall BC for TreeMesh looks good now. Is there a corresponding test that was added?

[maxbertrand1996]

Not yet.

What are you thinking of?
Like a well balanced test?

[maxbertrand1996]

I added a well balanced test which is analogue to examples/tree_2d_dgsem/elixir_shallowwater_well_balanced.jl called examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wall.jl which uses the slip wall condition as a boundary condition.

The original elixir had the lake at rest error at the end of the simulation of $9.90775192e-15$ and the new elixir of $1.11096018e-14$

[maxbertrand1996]

I also added another elixir examples/tree_2d_dgsem/elixir_shallowwater_source_terms_dirichlet.jl analogous to examples\\tree_2d_dgsem\\elixir_shallowwater_source_terms.jl which uses the inital condition as its boundary condition.
But for some reason when execute

julia> convergence_test("examples\\tree_2d_dgsem\\elixir_shallowwater_source_terms_dirichlet.jl", 4)

I get the following error message

ERROR: UndefVarError: result not defined
Stacktrace:
 [1] find_assignment(expr::Expr, destination::Symbol)
   @ Trixi C:\Users\maxbe\Repos\maxbertrand1996\Trixi.jl\src\auxiliary\special_elixirs.jl:219
 [2] extract_initial_resolution(elixir::String, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ Trixi C:\Users\maxbe\Repos\maxbertrand1996\Trixi.jl\src\auxiliary\special_elixirs.jl:239
 [3] convergence_test(mod::Module, elixir::String, iterations::Int64; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ Trixi C:\Users\maxbe\Repos\maxbertrand1996\Trixi.jl\src\auxiliary\special_elixirs.jl:62
 [4] convergence_test
   @ C:\Users\maxbe\Repos\maxbertrand1996\Trixi.jl\src\auxiliary\special_elixirs.jl:56 [inlined]
 [5] #convergence_test#1157
   @ C:\Users\maxbe\Repos\maxbertrand1996\Trixi.jl\src\auxiliary\special_elixirs.jl:142 [inlined]
 [6] convergence_test(elixir::String, iterations::Int64)
   @ Trixi C:\Users\maxbe\Repos\maxbertrand1996\Trixi.jl\src\auxiliary\special_elixirs.jl:142
 [7] top-level scope
   @ REPL[14]:1

caused by: UndefVarError: result not defined
Stacktrace:
 [1] find_assignment(expr::Expr, destination::Symbol)
   @ Trixi C:\Users\maxbe\Repos\maxbertrand1996\Trixi.jl\src\auxiliary\special_elixirs.jl:219
 [2] extract_initial_resolution(elixir::String, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ Trixi C:\Users\maxbe\Repos\maxbertrand1996\Trixi.jl\src\auxiliary\special_elixirs.jl:229
 [3] convergence_test(mod::Module, elixir::String, iterations::Int64; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ Trixi C:\Users\maxbe\Repos\maxbertrand1996\Trixi.jl\src\auxiliary\special_elixirs.jl:62
 [4] convergence_test
   @ C:\Users\maxbe\Repos\maxbertrand1996\Trixi.jl\src\auxiliary\special_elixirs.jl:56 [inlined]
 [5] #convergence_test#1157
   @ C:\Users\maxbe\Repos\maxbertrand1996\Trixi.jl\src\auxiliary\special_elixirs.jl:142 [inlined]
 [6] convergence_test(elixir::String, iterations::Int64)
   @ Trixi C:\Users\maxbe\Repos\maxbertrand1996\Trixi.jl\src\auxiliary\special_elixirs.jl:142
 [7] top-level scope
   @ REPL[14]:1

Does anybody of you know what might be causing this?

[andrewwinters5000]

Does anybody of you know what might be causing this?

You need to add a carriage return to the end of the new elixir_shallowwater_source_terms_dirichlet.jl such that the convergence_test parses the code properly, then it should work

[sloede]

You need to add a carriage return to the end of the new elixir_shallowwater_source_terms_dirichlet.jl such that the convergence_test parses the code properly, then it should work

😳 Uh. Is this a known (and fixable) issue in the implementation of convergence_test? This seems like an inconvenient "feature"...

[andrewwinters5000]

You need to add a carriage return to the end of the new elixir_shallowwater_source_terms_dirichlet.jl such that the convergence_test parses the code properly, then it should work

😳 Uh. Is this a known (and fixable) issue in the implementation of convergence_test? This seems like an inconvenient "feature"...

Not sure how easy it is to fix. I played around with the code for about twenty minutes until I noticed that without this carriage return the command expr = Meta.parse("begin $code end") on line 225 of src/auxiliary/special_elixirs.jl returns

$(Expr(:incomplete, "incomplete: \"begin\" at none:1 requires end"))

because the end somehow gets "lost"

[ranocha]

Not sure how easy it is to fix. I played around with the code for about twenty minutes until I noticed that without this carriage return the command expr = Meta.parse("begin $code end") on line 225 of src/auxiliary/special_elixirs.jl returns

$(Expr(:incomplete, "incomplete: \"begin\" at none:1 requires end"))

because the end somehow gets "lost"

Should be fixed by #1280

[andrewwinters5000]

Not sure how easy it is to fix. I played around with the code for about twenty minutes until I noticed that without this carriage return the command expr = Meta.parse("begin $code end") on line 225 of src/auxiliary/special_elixirs.jl returns

$(Expr(:incomplete, "incomplete: \"begin\" at none:1 requires end"))

because the end somehow gets "lost"

Should be fixed by #1280

Yes, the change from #1280 fixes this newline issue

[maxbertrand1996]

Thanks a lot, the convergence_test works now

####################################################################################################
l2
h                   h_v1                h_v2                b
error     EOC       error     EOC       error     EOC       error     EOC
1.66e-03  -         2.29e-02  -         3.45e-02  -         6.27e-05  -
1.10e-04  3.91      1.03e-03  4.47      2.50e-03  3.78      3.94e-06  3.99
6.91e-06  4.00      4.68e-05  4.46      6.77e-05  5.21      2.47e-07  4.00
4.34e-07  3.99      2.73e-06  4.10      2.98e-06  4.51      1.54e-08  4.00

mean      3.97      mean      4.34      mean      4.50      mean      4.00
----------------------------------------------------------------------------------------------------
linf
h                   h_v1                h_v2                b
error     EOC       error     EOC       error     EOC       error     EOC
1.37e-02  -         1.02e-01  -         1.98e-01  -         1.82e-04  -
1.22e-03  3.50      6.27e-03  4.03      1.08e-02  4.20      1.21e-05  3.91
7.54e-05  4.01      3.54e-04  4.15      4.47e-04  4.60      7.70e-07  3.98
4.75e-06  3.99      2.32e-05  3.93      2.24e-05  4.32      4.83e-08  3.99

mean      3.83      mean      4.04      mean      4.37      mean      3.96
----------------------------------------------------------------------------------------------------
Dict{Symbol, Any} with 3 entries:
  :variables => ("h", "h_v1", "h_v2", "b")
  :l2        => [3.96622, 4.34438, 4.49952, 3.99681]
  :linf      => [3.83252, 4.03565, 4.37067, 3.95958]

The results are even better than for elixir_shallowwater_source_terms.jl which uses periodic boundaries which might be expected.

####################################################################################################
l2
h                   h_v1                h_v2                b
error     EOC       error     EOC       error     EOC       error     EOC
1.65e-03  -         2.14e-02  -         2.70e-02  -         6.27e-05  -
1.11e-04  3.90      1.02e-03  4.39      2.93e-03  3.20      3.94e-06  3.99
6.91e-06  4.00      4.44e-05  4.52      6.82e-05  5.43      2.47e-07  4.00
4.34e-07  4.00      2.62e-06  4.08      2.92e-06  4.55      1.54e-08  4.00

mean      3.96      mean      4.33      mean      4.39      mean      4.00
----------------------------------------------------------------------------------------------------
linf
h                   h_v1                h_v2                b
error     EOC       error     EOC       error     EOC       error     EOC
1.34e-02  -         9.76e-02  -         1.61e-01  -         1.82e-04  -
1.21e-03  3.46      6.40e-03  3.93      1.29e-02  3.64      1.21e-05  3.91
7.55e-05  4.01      3.19e-04  4.33      4.76e-04  4.76      7.70e-07  3.98
4.74e-06  3.99      1.87e-05  4.09      2.41e-05  4.30      4.83e-08  3.99

mean      3.82      mean      4.12      mean      4.23      mean      3.96
----------------------------------------------------------------------------------------------------
Dict{Symbol, Any} with 3 entries:
  :variables => ("h", "h_v1", "h_v2", "b")
  :l2        => [3.9641, 4.33063, 4.39123, 3.99681]
  :linf      => [3.82165, 4.11574, 4.23446, 3.95958]

I also added the example elixirs to test_tree_2d_shallowwater.jl

[maxbertrand1996]

I get shown that there are still changes requested for this PR, but I do not see what needs to be addressed.

[ranocha]

I get shown that there are still changes requested for this PR, but I do not see what needs to be addressed.

@andrewwinters5000 requested changes before, so he needs to approve the new version to clear the old status and allow merging this PR.

[andrewwinters5000]

I get shown that there are still changes requested for this PR, but I do not see what needs to be addressed.

@andrewwinters5000 requested changes before, so he needs to approve the new version to clear the old status and allow merging this PR.

Yes, I will review everything today.

[andrewwinters5000]

The new tests and feature look good. I just left one last request to see if we can simplify the code.

[andrewwinters5000]

Sorry for the confusion about the nonconservative flux. The suggested simplifications should work and shrink the code.

[maxbertrand1996]

I hope now everything is fine

[andrewwinters5000]

Yes, everything looks fine now. This can merge once the tests all pass.