Skip to content

Commit fb1c625

Browse files
DanielDoehringranocha
DanielDoehring
and
ranocha
authored
Couette + Poiseuille flow (#2419)
* add poiseuille * Couette * comment * rename * Update test/test_parabolic_2d.jl * Update test/test_parabolic_2d.jl * prim errors --------- Co-authored-by: Hendrik Ranocha <[email protected]>
1 parent 3872868 commit fb1c625

File tree

3 files changed

+333
-0
lines changed

3 files changed

+333
-0
lines changed

‎examples/p4est_2d_dgsem/elixir_navierstokes_couette_flow.jl‎

Lines changed: 141 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,141 @@
1+
using OrdinaryDiffEqLowStorageRK
2+
using Trixi
3+
4+
###############################################################################
5+
# semidiscretization of the ideal compressible Navier-Stokes equations
6+
7+
# Fluid parameters
8+
gamma() = 1.4
9+
prandtl_number() = 0.71
10+
11+
# Parameters for compressible, yet viscous set up
12+
Re() = 100
13+
Ma() = 1.2
14+
15+
# Parameters that can be freely chosen
16+
v_top() = 1
17+
rho_in() = 1
18+
height() = 1.0
19+
20+
# Parameters that follow from Reynolds and Mach number + adiabatic index gamma
21+
nu() = v_top() * height() / Re()
22+
23+
c() = v_top() / Ma()
24+
p_over_rho() = c()^2 / gamma()
25+
p_in() = rho_in() * p_over_rho()
26+
mu() = rho_in() * nu()
27+
28+
equations = CompressibleEulerEquations2D(gamma())
29+
equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu(),
30+
Prandtl = prandtl_number())
31+
32+
# Set inflow, impose known/expected velocity profile
33+
@inline function initial_condition(x, t, equations)
34+
v1 = x[2] / height() * v_top() # Linear profile from 0 to v_top
35+
v2 = 0.0
36+
37+
prim = SVector(rho_in(), v1, v2, p_in())
38+
return prim2cons(prim, equations)
39+
end
40+
41+
len() = 5 * height() # Roughly constant at this len of the channel
42+
coordinates_min = (0.0, 0.0) # minimum coordinates (min(x), min(y))
43+
coordinates_max = (len(), height()) # maximum coordinates (max(x), max(y))
44+
45+
trees_per_dimension = (36, 12)
46+
mesh = P4estMesh(trees_per_dimension, polydeg = 1,
47+
coordinates_min = coordinates_min, coordinates_max = coordinates_max,
48+
periodicity = (false, false))
49+
50+
# Free outflow
51+
function boundary_condition_outflow(u_inner, normal_direction::AbstractVector, x, t,
52+
surface_flux_function,
53+
equations::CompressibleEulerEquations2D)
54+
# Calculate the boundary flux entirely from the internal solution state
55+
return Trixi.flux(u_inner, normal_direction, equations)
56+
end
57+
58+
### Hyperbolic boundary conditions ###
59+
bs_hyperbolic = Dict(:x_neg => BoundaryConditionDirichlet(initial_condition), # Weakly enforced inflow BC
60+
:x_pos => boundary_condition_outflow, # Free outflow/extended domain
61+
# Top/Bottom of channel: Walls
62+
:y_neg => boundary_condition_slip_wall,
63+
:y_pos => boundary_condition_slip_wall)
64+
65+
### Parabolic boundary conditions ###
66+
67+
velocity_bc_top_left = NoSlip((x, t, equations) -> SVector(x[2] / height() * v_top(), 0))
68+
# Use isothermal for inflow - adiabatic should also work
69+
heat_bc_top_left = Isothermal() do x, t, equations_parabolic
70+
Trixi.temperature(initial_condition(x, t,
71+
equations_parabolic),
72+
equations_parabolic)
73+
end
74+
bc_parabolic_top_left = BoundaryConditionNavierStokesWall(velocity_bc_top_left,
75+
heat_bc_top_left)
76+
77+
velocity_bc_bottom = NoSlip((x, t, equations) -> SVector(0, 0))
78+
heat_bc_bottom = Adiabatic((x, t, equations) -> 0)
79+
boundary_condition_bottom = BoundaryConditionNavierStokesWall(velocity_bc_bottom,
80+
heat_bc_bottom)
81+
82+
# On right end: Just copy the state/gradients
83+
@inline function boundary_condition_copy(flux_inner,
84+
u_inner,
85+
normal::AbstractVector,
86+
x, t,
87+
operator_type::Trixi.Gradient,
88+
equations::CompressibleNavierStokesDiffusion2D{GradientVariablesPrimitive})
89+
return u_inner
90+
end
91+
@inline function boundary_condition_copy(flux_inner,
92+
u_inner,
93+
normal::AbstractVector,
94+
x, t,
95+
operator_type::Trixi.Divergence,
96+
equations::CompressibleNavierStokesDiffusion2D{GradientVariablesPrimitive})
97+
return flux_inner
98+
end
99+
100+
bcs_parabolic = Dict(:x_neg => bc_parabolic_top_left,
101+
:x_pos => boundary_condition_copy,
102+
:y_neg => boundary_condition_bottom,
103+
:y_pos => bc_parabolic_top_left)
104+
105+
solver = DGSEM(polydeg = 3, surface_flux = flux_hll)
106+
107+
semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
108+
initial_condition, solver,
109+
boundary_conditions = (bs_hyperbolic,
110+
bcs_parabolic))
111+
112+
###############################################################################
113+
114+
tspan = (0, 5)
115+
ode = semidiscretize(semi, tspan)
116+
117+
summary_callback = SummaryCallback()
118+
119+
analysis_interval = 1000
120+
analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
121+
extra_analysis_errors = (:l2_error_primitive,
122+
:linf_error_primitive))
123+
124+
alive_callback = AliveCallback(analysis_interval = analysis_interval)
125+
126+
save_solution = SaveSolutionCallback(dt = 1.0,
127+
save_initial_solution = true,
128+
save_final_solution = true,
129+
solution_variables = cons2prim)
130+
131+
callbacks = CallbackSet(summary_callback,
132+
analysis_callback,
133+
alive_callback,
134+
save_solution)
135+
136+
###############################################################################
137+
138+
time_int_tol = 1e-7
139+
sol = solve(ode, RDPK3SpFSAL49();
140+
abstol = time_int_tol, reltol = time_int_tol,
141+
ode_default_options()..., callback = callbacks)

‎examples/p4est_2d_dgsem/elixir_navierstokes_poiseuille_flow.jl‎

Lines changed: 142 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,142 @@
1+
using OrdinaryDiffEqLowStorageRK
2+
using Trixi
3+
4+
###############################################################################
5+
# semidiscretization of the ideal compressible Navier-Stokes equations
6+
7+
# Fluid parameters
8+
const gamma = 1.4
9+
const prandtl_number = 0.72
10+
11+
# Parameters for compressible, yet viscous set up
12+
const Re = 1000
13+
const Ma = 0.8
14+
15+
# Parameters that can be freely chosen
16+
const v_in = 1
17+
const rho_in = 1
18+
const height = 1.0
19+
20+
# Parameters that follow from Reynolds and Mach number + adiabatic index gamma
21+
const nu = v_in * height / Re
22+
23+
const c = v_in / Ma
24+
const p_over_rho = c^2 / gamma
25+
const p_in = rho_in * p_over_rho
26+
const mu = rho_in * nu
27+
28+
equations = CompressibleEulerEquations2D(gamma)
29+
equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu,
30+
Prandtl = prandtl_number)
31+
32+
# Set naive inflow
33+
@inline function initial_condition(x, t, equations)
34+
rho = rho_in
35+
v1 = v_in
36+
v2 = 0.0
37+
p = p_in
38+
39+
prim = SVector(rho, v1, v2, p)
40+
return prim2cons(prim, equations)
41+
end
42+
43+
const len = 10 * height # Roughly constant at this len of the channel
44+
coordinates_min = (0.0, 0.0) # minimum coordinates (min(x), min(y))
45+
coordinates_max = (len, height) # maximum coordinates (max(x), max(y))
46+
47+
trees_per_dimension = (36, 12)
48+
mesh = P4estMesh(trees_per_dimension, polydeg = 1,
49+
coordinates_min = coordinates_min, coordinates_max = coordinates_max,
50+
periodicity = (false, false))
51+
52+
# Free outflow
53+
function boundary_condition_outflow(u_inner, normal_direction::AbstractVector, x, t,
54+
surface_flux_function,
55+
equations::CompressibleEulerEquations2D)
56+
# Calculate the boundary flux entirely from the internal solution state
57+
return Trixi.flux(u_inner, normal_direction, equations)
58+
end
59+
60+
### Hyperbolic boundary conditions ###
61+
bs_hyperbolic = Dict(:x_neg => BoundaryConditionDirichlet(initial_condition), # Weakly enforced inflow BC
62+
:x_pos => boundary_condition_outflow, # Free outflow/extended domain
63+
# Top/Bottom of channel: Walls
64+
:y_neg => boundary_condition_slip_wall,
65+
:y_pos => boundary_condition_slip_wall)
66+
67+
### Parabolic boundary conditions ###
68+
69+
velocity_bc_inflow = NoSlip((x, t, equations) -> SVector(v_in, 0))
70+
# Use isothermal for inflow - adiabatic should also work
71+
heat_bc_inflow = Isothermal() do x, t, equations_parabolic
72+
Trixi.temperature(initial_condition(x, t,
73+
equations_parabolic),
74+
equations_parabolic)
75+
end
76+
bc_parabolic_inflow = BoundaryConditionNavierStokesWall(velocity_bc_inflow, heat_bc_inflow)
77+
78+
velocity_bc_wall = NoSlip((x, t, equations) -> SVector(0, 0))
79+
heat_bc_wall = Adiabatic((x, t, equations) -> 0)
80+
boundary_condition_wall = BoundaryConditionNavierStokesWall(velocity_bc_wall, heat_bc_wall)
81+
82+
# On right end: Just copy the state/gradients
83+
@inline function boundary_condition_copy(flux_inner,
84+
u_inner,
85+
normal::AbstractVector,
86+
x, t,
87+
operator_type::Trixi.Gradient,
88+
equations::CompressibleNavierStokesDiffusion2D{GradientVariablesPrimitive})
89+
return u_inner
90+
end
91+
@inline function boundary_condition_copy(flux_inner,
92+
u_inner,
93+
normal::AbstractVector,
94+
x, t,
95+
operator_type::Trixi.Divergence,
96+
equations::CompressibleNavierStokesDiffusion2D{GradientVariablesPrimitive})
97+
return flux_inner
98+
end
99+
100+
bcs_parabolic = Dict(:x_neg => bc_parabolic_inflow,
101+
:x_pos => boundary_condition_copy,
102+
# Top/Bottom of channel: Walls
103+
:y_neg => boundary_condition_wall,
104+
:y_pos => boundary_condition_wall)
105+
106+
solver = DGSEM(polydeg = 3, surface_flux = flux_hll)
107+
108+
semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
109+
initial_condition, solver,
110+
boundary_conditions = (bs_hyperbolic,
111+
bcs_parabolic))
112+
113+
###############################################################################
114+
115+
tspan = (0, 10)
116+
ode = semidiscretize(semi, tspan)
117+
118+
summary_callback = SummaryCallback()
119+
120+
analysis_interval = 1000
121+
analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
122+
extra_analysis_errors = (:l2_error_primitive,
123+
:linf_error_primitive))
124+
125+
alive_callback = AliveCallback(analysis_interval = analysis_interval)
126+
127+
save_solution = SaveSolutionCallback(dt = 1.0,
128+
save_initial_solution = true,
129+
save_final_solution = true,
130+
solution_variables = cons2prim)
131+
132+
callbacks = CallbackSet(summary_callback,
133+
analysis_callback,
134+
alive_callback,
135+
save_solution)
136+
137+
###############################################################################
138+
139+
time_int_tol = 1e-7
140+
sol = solve(ode, RDPK3SpFSAL49();
141+
abstol = time_int_tol, reltol = time_int_tol,
142+
ode_default_options()..., callback = callbacks)

‎test/test_parabolic_2d.jl‎

Lines changed: 50 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -945,6 +945,31 @@ end
945945
end
946946
end
947947

948+
@trixi_testset "elixir_navierstokes_poiseuille_flow.jl" begin
949+
@test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem",
950+
"elixir_navierstokes_poiseuille_flow.jl"),
951+
l2=[
952+
0.028671228188785286,
953+
0.2136420195921885,
954+
0.009953689550858224,
955+
0.13216036594768157
956+
],
957+
linf=[
958+
0.30901218409540543,
959+
1.3488655161645846,
960+
0.1304661713119874,
961+
1.2094591729756736],
962+
tspan=(0.0, 1.0))
963+
# Ensure that we do not have excessive memory allocations
964+
# (e.g., from type instabilities)
965+
let
966+
t = sol.t[end]
967+
u_ode = sol.u[end]
968+
du_ode = similar(u_ode)
969+
@test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
970+
end
971+
end
972+
948973
@trixi_testset "elixir_navierstokes_kelvin_helmholtz_instability_sc_subcell.jl" begin
949974
@test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem",
950975
"elixir_navierstokes_kelvin_helmholtz_instability_sc_subcell.jl"),
@@ -1000,6 +1025,31 @@ end
10001025
end
10011026
end
10021027

1028+
@trixi_testset "elixir_navierstokes_couette_flow.jl" begin
1029+
@test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem",
1030+
"elixir_navierstokes_couette_flow.jl"),
1031+
l2=[
1032+
0.009585252225488753,
1033+
0.007939233099864973,
1034+
0.0007617512688442657,
1035+
0.027229870237669436
1036+
],
1037+
linf=[
1038+
0.027230029149270862,
1039+
0.027230451118692933,
1040+
0.0038642959675975713,
1041+
0.04738248734987671
1042+
])
1043+
# Ensure that we do not have excessive memory allocations
1044+
# (e.g., from type instabilities)
1045+
let
1046+
t = sol.t[end]
1047+
u_ode = sol.u[end]
1048+
du_ode = similar(u_ode)
1049+
@test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
1050+
end
1051+
end
1052+
10031053
@trixi_testset "elixir_navierstokes_blast_reflective.jl" begin
10041054
@test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem",
10051055
"elixir_navierstokes_blast_reflective.jl"),

0 commit comments

Comments
 (0)