Skip to content

Commit

Permalink
set true discontinuities in 2D
Browse files Browse the repository at this point in the history
  • Loading branch information
patrickersing committed Apr 22, 2024
1 parent 2eec19d commit db5f660
Show file tree
Hide file tree
Showing 4 changed files with 168 additions and 82 deletions.
41 changes: 41 additions & 0 deletions examples/tree_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,47 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,

tspan = (0.0, 2.0)
ode = semidiscretize(semi, tspan)
###############################################################################
#=
Workaround for TreeMesh2D to set true discontinuities for debugging and testing.
Essentially, this is a slight augmentation of the `compute_coefficients` where the `x` node values
passed here are slightly perturbed in order to set a true discontinuity that avoids the doubled
value of the LGL nodes at a particular element interface.
=#

# 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 i in eachnode(semi.solver), j in eachnode(semi.solver)
x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations,
semi.solver, i, j, element)
# Changing the node positions passed to the initial condition by the minimum
# amount possible with the current type of floating point numbers allows setting
# discontinuous initial data in a simple way. In particular, a check like `if x < x_jump`
# works if the jump location `x_jump` is at the position of an interface.
if i == 1 && j == 1 # bottom left corner
x_node = SVector(nextfloat(x_node[1]), nextfloat(x_node[2]))
elseif i == 1 && j == nnodes(semi.solver) # top left corner
x_node = SVector(nextfloat(x_node[1]), prevfloat(x_node[2]))
elseif i == nnodes(semi.solver) && j == 1 # bottom right corner
x_node = SVector(prevfloat(x_node[1]), nextfloat(x_node[2]))
elseif i == nnodes(semi.solver) && j == nnodes(semi.solver) # top right corner
x_node = SVector(prevfloat(x_node[1]), prevfloat(x_node[2]))
elseif i == 1 # left boundary
x_node = SVector(nextfloat(x_node[1]), x_node[2])
elseif j == 1 # bottom boundary
x_node = SVector(x_node[1], nextfloat(x_node[2]))
elseif i == nnodes(semi.solver) # right boundary
x_node = SVector(prevfloat(x_node[1]), x_node[2])
elseif j == nnodes(semi.solver) # top boundary
x_node = SVector(x_node[1], prevfloat(x_node[2]))
end

u_node = initial_condition_dam_break(x_node, first(tspan), equations)
Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, j, element)
end
end

###############################################################################
# Callbacks
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -10,8 +10,10 @@ using TrixiShallowWater
equations = ShallowWaterMultiLayerEquations2D(gravity_constant = 1.0,
rhos = (0.9, 0.95, 1.0))

# This academic testcase sets up a discontinuous bottom topography
# function and initial condition to test entropy conservation.
# This test case uses a special work around to setup a truly discontinuous bottom topography
# function and initial condition for this academic testcase of entropy conservation. First, a
# dummy initial_condition_dam_break is introduced to create the semidiscretization. Then the initial
# condition is reset with the true discontinuous values from initial_condition_discontinuous_dam_break.

function initial_condition_dam_break(x, t, equations::ShallowWaterMultiLayerEquations2D)
# Bottom topography
Expand Down Expand Up @@ -65,6 +67,49 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition,
tspan = (0.0, 0.5)
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
# test case and initial condition.
# 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 specific mesh loaded above!

function initial_condition_discontinuous_dam_break(x, t, element_id,
equations::ShallowWaterMultiLayerEquations2D)
# Bottom topography
b = 0.3 * exp(-0.5 * ((x[1] - sqrt(2) / 2)^2 + (x[2] - sqrt(2) / 2)^2))

# Left side of discontinuity
IDs = [1, 2, 5, 6, 9, 10, 13, 14]
if element_id in IDs
H = SVector(1.0, 0.8, 0.6)
# Right side of discontinuity
else
H = SVector(0.9, 0.7, 0.5)
b += 0.1
end

v1 = zero(H)
v2 = zero(H)
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_dam_break(x_node, first(tspan), element,
equations)
Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, j, element)
end
end

###############################################################################
# Callbacks

Expand Down
80 changes: 40 additions & 40 deletions test/test_tree_2d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -617,28 +617,28 @@ end # 2LSWE
@test_trixi_include(joinpath(EXAMPLES_DIR,
"elixir_shallowwater_multilayer_dam_break.jl"),
l2=[
0.006943714920464853,
0.007232726823994391,
0.011985810753673342,
0.005778397278079132,
0.005736810674427555,
0.006946616973320985,
6.143372810603547e-5,
6.116048140831271e-5,
7.85179348120683e-5,
0.005455457843589704,
0.006867196239306527,
0.007212095601190101,
0.014309127811972761,
0.005801387118671292,
0.0057555558162007536,
0.006951246408270715,
6.084817535430214e-5,
6.068733435363927e-5,
7.807791370377477e-5,
0.003857583749542185,
],
linf=[
0.0264647222938314,
0.02644739838686727,
0.09689908398594305,
0.021624857501702666,
0.02083541825119268,
0.029889274751729353,
0.0003438446930977006,
0.00033898478545986335,
0.0006204362806211161,
0.05596747200337038,
0.028792046975677804,
0.036503344705121676,
0.19267105917517297,
0.028685238833612656,
0.02761561086221269,
0.029892140357308587,
0.0004158435510054702,
0.0004119057098728111,
0.0007709896609887244,
0.10000011323773067,
],
tspan=(0.0, 0.25))
# Ensure that we do not have excessive memory allocations
Expand All @@ -655,28 +655,28 @@ end # 2LSWE
@test_trixi_include(joinpath(EXAMPLES_DIR,
"elixir_shallowwater_multilayer_dam_break.jl"),
l2=[
0.006846433642482368,
0.007157521853802785,
0.012242149312850049,
0.00562742302611423,
0.005619979904755309,
0.006853160451226712,
5.8363153191008005e-5,
5.848194044876578e-5,
7.584883097507706e-5,
0.005455457843589704,
0.007162130568658371,
0.007507487998680027,
0.01500636533827742,
0.005593279150854283,
0.005603720137202845,
0.0068519457928639385,
5.748086634734491e-5,
5.7787474070077486e-5,
7.492163249697633e-5,
0.003857583749542185,
],
linf=[
0.026631073359625668,
0.025615528815443434,
0.0817112552789852,
0.020564308208119463,
0.020040028098637117,
0.024332298183399693,
0.0003206729530360014,
0.00031809818924231675,
0.00042369787238329603,
0.05596747200337038,
0.037711432402825124,
0.0427565001934887,
0.17275386567546575,
0.025348821041845597,
0.02434895938442919,
0.03396954814003092,
0.0003026537474088657,
0.00029947595638190504,
0.00037146307629417057,
0.10000011323773067,
],
surface_flux=(flux_lax_friedrichs,
flux_nonconservative_ersing_etal),
Expand Down
80 changes: 40 additions & 40 deletions test/test_unstructured_2d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -600,28 +600,28 @@ end # 2LSWE
@test_trixi_include(joinpath(EXAMPLES_DIR,
"elixir_shallowwater_multilayer_dam_break.jl"),
l2=[
0.008264086240458098,
0.008904024684630463,
0.015040884852261934,
0.0068380853981447965,
0.006803190764735809,
0.007681935438533932,
0.00014126104190197007,
0.00013800738645477302,
0.00015518676788954237,
0.0062500008375575705,
0.008133096076522545,
0.009059379390936938,
0.0156472431836693,
0.006787446774526842,
0.0067495775839911685,
0.007737394569507473,
0.0001835435993951076,
0.00018347939156198144,
0.00023500246849234892,
0.004003203849568449,
],
linf=[
0.0271281524526652,
0.038446397522241216,
0.10189873392468171,
0.02279784901518067,
0.02344053175169211,
0.0384823318069141,
0.0015613916350439184,
0.0015293762865072726,
0.0020984165518842324,
0.05808926980412543,
0.0268742067416933,
0.05243441379266703,
0.1590063807340898,
0.028581575095277572,
0.024486565205936787,
0.03953111184090555,
0.002625120329243771,
0.0026387938388929785,
0.0038007822188701103,
0.10000000026183736,
],
tspan=(0.0, 0.25))
# Ensure that we do not have excessive memory allocations
Expand All @@ -638,28 +638,28 @@ end # 2LSWE
@test_trixi_include(joinpath(EXAMPLES_DIR,
"elixir_shallowwater_multilayer_dam_break.jl"),
l2=[
0.008289485483638842,
0.008744140233824497,
0.014162544835509572,
0.0067108677626928565,
0.006699962946263383,
0.0075441143283640315,
0.00012772255623776208,
0.0001243555119128215,
0.00013480820924093037,
0.0062500008375575705,
0.0089213866827554,
0.009502261925352482,
0.017690193299797062,
0.006845800468099525,
0.006849466485572439,
0.00761550757687342,
0.0004491107702429648,
0.00043570105077582146,
0.0005832036370132536,
0.004003203849568449,
],
linf=[
0.024574998802901732,
0.024386474982493633,
0.09088941363545519,
0.01828534534858913,
0.017626892238781267,
0.02546243728251761,
0.0007231987972067275,
0.0007131387402676157,
0.0007959906366350036,
0.05808926980412543,
0.04276512135712615,
0.054061877851283885,
0.18896205511083025,
0.018902398329971263,
0.017901194463979406,
0.025135202096835247,
0.004065137704454743,
0.0038290081477974037,
0.005930216968236107,
0.10000000026183736,
],
surface_flux=(flux_lax_friedrichs,
flux_nonconservative_ersing_etal),
Expand Down

0 comments on commit db5f660

Please sign in to comment.