search

trixi-framework / Trixi.jl

Trixi.jl: Adaptive high-order numerical simulations of conservation laws in Julia
https://trixi-framework.github.io/Trixi.jl
MIT License
541 stars 110 forks source link

Specifying partially periodic boundary conditions with `P4estMesh` is inconsistent with `TreeMesh` #1554

Open apey236 opened 1 year ago

apey236 commented 1 year ago

The main difference is that we cannot specify a boundary_condition_periodic for P4estMesh. This example illustrates this:

using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations

equations = CompressibleEulerEquations3D(1.4)

initial_condition = initial_condition_constant
dirichlet = BoundaryConditionDirichlet(initial_condition)
boundary_conditions = Dict( :x_neg => boundary_condition_periodic,
                         :x_pos => boundary_condition_periodic,
                         :y_neg => dirichlet,
                         :y_pos => dirichlet,
                         :z_neg => boundary_condition_periodic,
                         :z_pos => boundary_condition_periodic)

solver = DGSEM(polydeg=4, surface_flux=flux_lax_friedrichs,
               volume_integral=VolumeIntegralWeakForm())

coordinates_min = (-1.0, -1.0, -1.0) # minimum coordinates (min(x), min(y), min(z))
coordinates_max = ( 1.0,  1.0,  1.0) # maximum coordinates (max(x), max(y), max(z))

trees_per_dimension = (2, 2, 2)

mesh = P4estMesh(trees_per_dimension, polydeg=3,
                 coordinates_min=coordinates_min, coordinates_max=coordinates_max,
                 periodicity=(true, false, true), initial_refinement_level=2)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions=boundary_conditions)

###############################################################################
# ODE solvers, callbacks etc.

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

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

stepsize_callback = StepsizeCallback(cfl=1.2)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        stepsize_callback)

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

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary

The stack trace that I get is: 

julia> include("../examples/p4est_3d_dgsem/elixir_euler_free_stream_nonperiodic.jl")
ERROR: Key :x_pos is not a valid boundary name
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] initialize!(boundary_types_container::Trixi.UnstructuredSortedBoundaryTypes{2, Tuple{Trixi.BoundaryConditionPeriodic, BoundaryConditionDirichlet{typeof(initial_condition_constant)}}}, cache::NamedTuple{(:elements, :interfaces, :boundaries, :mortars, :fstar_threaded, :fstar_tmp_threaded, :u_threaded), Tuple{Trixi.P4estElementContainer{3, Float64, Float64, 4, 5, 6}, Trixi.P4estInterfaceContainer{3, Float64, 5}, Trixi.P4estBoundaryContainer{3, Float64, 4}, Trixi.P4estMortarContainer{3, Float64, 4, 6}, Vector{Array{Float64, 4}}, Vector{Array{Float64, 3}}, Vector{Array{Float64, 3}}}})
   @ Trixi ~/.julia/dev/Trixi/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl:71
 [3] Trixi.UnstructuredSortedBoundaryTypes(boundary_conditions::Dict{Symbol, Any}, cache::NamedTuple{(:elements, :interfaces, :boundaries, :mortars, :fstar_threaded, :fstar_tmp_threaded, :u_threaded), Tuple{Trixi.P4estElementContainer{3, Float64, Float64, 4, 5, 6}, Trixi.P4estInterfaceContainer{3, Float64, 5}, Trixi.P4estBoundaryContainer{3, Float64, 4}, Trixi.P4estMortarContainer{3, Float64, 4, 6}, Vector{Array{Float64, 4}}, Vector{Array{Float64, 3}}, Vector{Array{Float64, 3}}}})
   @ Trixi ~/.julia/dev/Trixi/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl:36
 [4] digest_boundary_conditions(boundary_conditions::Dict{Symbol, Any}, mesh::P4estMesh{3, Float64, Static.False, P4est.PointerWrappers.PointerWrapper{P4est.LibP4est.p8est}, P4est.PointerWrappers.PointerWrapper{P4est.LibP4est.p8est_ghost_t}, 5, 4}, solver::DGSEM{LobattoLegendreBasis{Float64, 5, SVector{5, Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}, Trixi.LobattoLegendreMortarL2{Float64, 5, Matrix{Float64}, Matrix{Float64}}, SurfaceIntegralWeakForm{FluxLaxFriedrichs{typeof(max_abs_speed_naive)}}, VolumeIntegralWeakForm}, cache::NamedTuple{(:elements, :interfaces, :boundaries, :mortars, :fstar_threaded, :fstar_tmp_threaded, :u_threaded), Tuple{Trixi.P4estElementContainer{3, Float64, Float64, 4, 5, 6}, Trixi.P4estInterfaceContainer{3, Float64, 5}, Trixi.P4estBoundaryContainer{3, Float64, 4}, Trixi.P4estMortarContainer{3, Float64, 4, 6}, Vector{Array{Float64, 4}}, Vector{Array{Float64, 3}}, Vector{Array{Float64, 3}}}})
   @ Trixi ~/.julia/dev/Trixi/src/semidiscretization/semidiscretization_hyperbolic.jl:204
 [5] SemidiscretizationHyperbolic(mesh::P4estMesh{3, Float64, Static.False, P4est.PointerWrappers.PointerWrapper{P4est.LibP4est.p8est}, P4est.PointerWrappers.PointerWrapper{P4est.LibP4est.p8est_ghost_t}, 5, 4}, equations::CompressibleEulerEquations3D{Float64}, initial_condition::typeof(initial_condition_constant), solver::DGSEM{LobattoLegendreBasis{Float64, 5, SVector{5, Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}, Trixi.LobattoLegendreMortarL2{Float64, 5, Matrix{Float64}, Matrix{Float64}}, SurfaceIntegralWeakForm{FluxLaxFriedrichs{typeof(max_abs_speed_naive)}}, VolumeIntegralWeakForm}; source_terms::Nothing, boundary_conditions::Dict{Symbol, Any}, RealT::Type, uEltype::Type, initial_cache::NamedTuple{(), Tuple{}})
   @ Trixi ~/.julia/dev/Trixi/src/semidiscretization/semidiscretization_hyperbolic.jl:71
 [6] top-level scope
   @ ~/.julia/dev/Trixi/examples/p4est_3d_dgsem/elixir_euler_free_stream_nonperiodic.jl:33
jlchan commented 1 year ago

It looks like not specifying boundary_condition_periodic for :xneg, :xpos, :zneg, :zpos fixes things. This is also consistent with https://github.com/jlchan/Trixi.jl/blob/b7be5856eba029a3d91257166be4ee62514ae0dd/examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl - I don't know why we didn't see that before...

ranocha commented 1 year ago

Good to know taht there is a way to make it work right now. But it is definitely confusing and should be fixed

jlchan commented 1 year ago

Related: https://github.com/trixi-framework/Trixi.jl/issues/1510