Solving 1D inviscid Burger's equation with time dependent BC

[ErikQ97]

Hello,

I am trying to solve 1D inviscid Burger's equation with time dependent BC (i.e. u(x=0,t) = sin(2 pi t)), while initially everywhere is zero i.e. u(x,t=0) = 0. However, the solution (t=1) is non-physical as shown below:

The codes is shown as follows:

using Trixi, OrdinaryDiffEq, Plots

function inlet_sine_bc(u_inner, orientation, direction, x, t, surface_flux::FluxLaxFriedrichs, equations::InviscidBurgersEquation1D)
    return SVector(sin(2π * t))
end

equations = InviscidBurgersEquation1D()

solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)
initial_condition(x, t, equations) = 0.0

boundary_conditions = (
    x_neg = inlet_sine_bc,
    x_pos = Trixi.BoundaryConditionDoNothing(),
)

semi = SemidiscretizationHyperbolic(
    TreeMesh(0.0, 1.0, initial_refinement_level=5, n_cells_max = 10^4, periodicity = false),
    equations, initial_condition, solver;
    boundary_conditions = boundary_conditions
)

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

controller = StepsizeCallback(cfl = 0.25)  
save_times = 0:0.01:1.0
sol = solve(ode, SSPRK33(), dt=1e-3, saveat=save_times, callback=controller)

plot(sol)

Any comments are appreciated!

[DanielDoehring]

So what you are actually specifying in this example is the boundary flux $F(u) = \sin(2 \pi t)$, not the state. See for instance the other examples on Burgers which set the fluxes:

Trixi.jl/examples/tree_1d_dgsem/elixir_burgers_rarefaction.jl

Lines 46 to 58 in 6b0f725

	boundary_condition_inflow = BoundaryConditionDirichlet(initial_condition_rarefaction)
	function boundary_condition_outflow(u_inner, orientation, normal_direction, x, t,
	surface_flux_function,
	equations::InviscidBurgersEquation1D)
	# Calculate the boundary flux entirely from the internal solution state
	flux = Trixi.flux(u_inner, normal_direction, equations)
	return flux
	end
	boundary_conditions = (x_neg = boundary_condition_inflow,
	x_pos = boundary_condition_outflow)

Trixi.jl/examples/tree_1d_dgsem/elixir_burgers_shock.jl

Lines 46 to 58 in 6b0f725

	boundary_condition_inflow = BoundaryConditionDirichlet(initial_condition_shock)
	function boundary_condition_outflow(u_inner, orientation, normal_direction, x, t,
	surface_flux_function,
	equations::InviscidBurgersEquation1D)
	# Calculate the boundary flux entirely from the internal solution state
	flux = Trixi.flux(u_inner, normal_direction, equations)
	return flux
	end
	boundary_conditions = (x_neg = boundary_condition_inflow,
	x_pos = boundary_condition_outflow)

In your case, you also probably want to pipe the state through BoundaryConditionDirichlet:

using Trixi, OrdinaryDiffEqSSPRK, Plots

u_init() = 0.5
function inlet_sine_bc(x, t, equation::InviscidBurgersEquation1D)
    return SVector(u_init() + sin(2π * t))
end

equations = InviscidBurgersEquation1D()

solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)
initial_condition(x, t, equations) = u_init()

boundary_conditions = (
    x_neg = BoundaryConditionDirichlet(inlet_sine_bc),
    x_pos = boundary_condition_do_nothing,
)

semi = SemidiscretizationHyperbolic(
    TreeMesh(0.0, 1.0, initial_refinement_level=5, n_cells_max = 10^4, periodicity = false),
    equations, initial_condition, solver;
    boundary_conditions = boundary_conditions
)

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

controller = StepsizeCallback(cfl = 0.25)  
sol = solve(ode, SSPRK33(), dt=1e-3, callback=controller)

plot(sol)

Note that I also set an initial state different from 0 such that there is actually some "flow" from left to right - otherwise stuff will "pile up" quite fast at the left end, and only very slow move through the domain.

[ErikQ97]

Thanks for your reply! So what is the difference between

boundary_conditions = (
    x_neg = inlet_sine_bc,
    x_pos = Trixi.BoundaryConditionDoNothing(),
)

and

boundary_conditions = (
    x_neg = BoundaryConditionDirichlet(inlet_sine_bc),
    x_pos = Trixi.BoundaryConditionDoNothing(),
)

[JoshuaLampert]

You can see the difference by looking at the implementation of BoundaryConditionDirichlet:

Trixi.jl/src/equations/equations.jl

Lines 164 to 182 in 3370de6

	@inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner,
	orientation_or_normal,
	direction,
	x, t,
	surface_flux_function,
	equations)
	u_boundary = boundary_condition.boundary_value_function(x, t, equations)
	# Calculate boundary flux
	if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
	flux = surface_flux_function(u_inner, u_boundary, orientation_or_normal,
	equations)
	else # u_boundary is "left" of boundary, u_inner is "right" of boundary
	flux = surface_flux_function(u_boundary, u_inner, orientation_or_normal,
	equations)
	end
	return flux
	end

As you can see BoundaryConditionDirichlet returns the flux and not the variables itself. Trixi.jl expects the functions you pass to x_neg and x_pos to return the flux on the boundary. In your first version the return value of your inlet_sine_bc is interpreted of the result when computing the numerical flux, which is not what you want. Therefore, you likely want the second version.

[ErikQ97]

You can see the difference by looking at the implementation of BoundaryConditionDirichlet:

Trixi.jl/src/equations/equations.jl

Lines 164 to 182 in 3370de6

@inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner,
orientation_or_normal,
direction,
x, t,
surface_flux_function,
equations)
u_boundary = boundary_condition.boundary_value_function(x, t, equations)

 # Calculate boundary flux 
 if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary 
     flux = surface_flux_function(u_inner, u_boundary, orientation_or_normal, 
                                  equations) 
 else # u_boundary is "left" of boundary, u_inner is "right" of boundary 
     flux = surface_flux_function(u_boundary, u_inner, orientation_or_normal, 
                                  equations) 
 end 

 return flux 

end
As you can see BoundaryConditionDirichlet returns the flux and not the variables itself. Trixi.jl expects the functions you pass to x_neg and x_pos to return the flux on the boundary. In your first version the return value of your inlet_sine_bc is interpreted of the result when computing the numerical flux, which is not what you want. Therefore, you likely want the second version.

Got it! Thank you!

[ranocha]

Please close this issue when your question has been answered.

[ErikQ97]

I tried a simpler case where initially everywhere is zero, i.e. u(x, t=0) = 0, then a constant boundary condition is applied as u(x=0, t) = 1.0. This is a very simple case which should result in a stationary state where u(x, t) = 1.0 at the end of the day (after the shock gone and the inlet constant wave fill in the domain). However, the solver still gives me something like

I think this solution is non-physical. Kindly let me know if I am wrong in codes or/and in analytics. Thank you!

The codes are as follow:

using Trixi, OrdinaryDiffEq, Plots

# Define the sine wave BC
function inlet_sine_bc(x, t, equation::InviscidBurgersEquation1D)
    return ones(length(t))
end

equations = InviscidBurgersEquation1D()

solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)
initial_condition(x, t, equations) = 0.0

boundary_conditions = (
    x_neg = BoundaryConditionDirichlet(inlet_sine_bc),
    x_pos = Trixi.BoundaryConditionDoNothing(),
)

semi = SemidiscretizationHyperbolic(
    TreeMesh(0.0, 1.0, initial_refinement_level=5, n_cells_max = 10^4, periodicity = false),
    equations, initial_condition, solver;
    boundary_conditions = boundary_conditions
)

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

controller = StepsizeCallback(cfl = 0.25)
sol = solve(ode, SSPRK33(), dt=1e-3, callback=controller)

plot(sol)

[DanielDoehring]

Actually you found somewhat like a bug/unexpected behaviour in the code: Since initially you can taken an arbitrary large timestep (as the IC is zero everywhere, the characteristic speed is 0, hence the stepsizeCallback will do only one step (which however becomes unstable over the course of the Runge-Kutta method) with maximum size. So either turn off the stepsize callback here and use a fixed timestep or switch to an adaptive method. Btw, if you would have used AnalysisCallback or AliveCallback you would have seen that only one timestep would be taken, thus hinting at the issue.

[ErikQ97]

Thanks for you reply! It makes great sense!

[ErikQ97]

I tried a sine wave as inlet BC which looks like

function inlet_bc(x, t, equation::InviscidBurgersEquation1D)
    wave = zeros(length(t))
    f = 1.0
    wave .= sin.(2 * π * f * t)
    
    return wave
end

and it give a result as

I think the first peak (or shock) is non-physical because the positive part of the sine wave only has one peak and will result in only one shock, also the characteristics move to the right thus the shock should also go to the right.

Please let me know if I am wrong in the codes or physics. Thank you!

[DanielDoehring]

Hard to tell, does not necessarily look like a bug to me.

[ErikQ97]

I did not see bugs from the codes actually, it seems to be exactly what I want. I am just curious why there are two shocks given a single period sine wave.