Skip to content

Commit

Permalink
Navier-Stokes in 2D on DGMulti and TreeMesh (#1165)
Browse files Browse the repository at this point in the history
* first draft of container for Navier-Stokes constants and fluxes

* remove unneeded temperature computation

* draft of elixir with possible boundary condition structure

* added manufactured solution and source term

* fix typo in function name for MMS

* update variable names for consistency. improve comments

* fix dumb typos in new equation system name

* actually export new equations

* add comment near variable_mapping declaration.

* Apply suggestions from code review

Co-authored-by: Hendrik Ranocha <[email protected]>

* parabolic equations now exported separately

* change name to CompressibleNavierStokesEquations2D

* export NS with proper name

* explicitly require compressible Euler in the type parameter

* name kinematic viscosity nu for consistency with Lattice-Boltzmann

* add promotion in constructor

* make Reynolds, Prandtl, Mach, and kappa keyword arguments

* update constructor call in elixir

* reduce computation by exploiting stress tensor symmetry

* fix unpacking of flux

* modifying parabolic cache creation

in cache, we assume we take the gradient of all hyperbolic variables. since the number of parabolic variables can differ from the number of hyperbolic variables, we pass in the hyperbolic equations to `create_cache_parabolic` now

* comments

* comments

* formatting and renaming equations to equations_hyperbolic


formatting comments

* fix unpacking of gradients in flux

* adding CNS BCs

* adding lid-driven cavity elixir

* adding variable transform, editing cons2prim for CNS

* add prim2cons for NS (inconsistent right now)

* add draft of DGMulti Navier Stokes convergence elixir

* converging solution using elixir for TreeMesh with BCs

* fixing DGMulti advection diffusion elixir convergence

* naming equations as parabolic/hyperbolic

* generalizing transform_variables

* add TODO


more todos

* additional checks on get_unsigned_normal

* adding isothermal BC

* commenting out unused CNS functions for now

* fix call to transform_variables

* comments and cleanup

* changing default solver and Re for cavity

* adding more advection diffusion tests

* label tests

* add gradient_variables field to SemidiscretizationHyperbolicParabolic

* Revert "add gradient_variables field to SemidiscretizationHyperbolicParabolic"

This reverts commit 063b602.

* allowing for specialization in transform_variables

adding a function `gradient_variable_transformation` which should get specialized if the gradient variables are not conservative variables

* formatting and comments

* reverting elixir

* comments

* standardizing time tol

* minor fix to CNS boundary flux for convenience

make it so that the density state is computed correctly, even though it's not used

* formatting + comments

* using primitive variables in viscous flux instead of conservative

* minor formatting

* add CNS tests

* fix test

* testing if scoping issues fix TreeMesh tests

* decreasing timestep tol for TreeMesh parabolic test

* enabling periodic BCs for TreeMesh + more tests

* fix density BC flux (unused, but could be useful)

* adding non-working TreeMesh elixirs

* adding AD checks

* standardizing parameters in convergence elixirs

* minor cleanup

* revert DGMulti CNS elixir setup back to the one used in tests

* adding TreeMesh CNS convergence test

* removing redundant elixir

* add more tests

* add more test

* Apply suggestions from code review

Co-authored-by: Hendrik Ranocha <[email protected]>

* add docstrings

Co-authored-by: Hendrik Ranocha <[email protected]>
Co-authored-by: Hendrik Ranocha <[email protected]>
Co-authored-by: Jesse Chan <[email protected]>
  • Loading branch information
4 people authored Jul 29, 2022
1 parent e8a24a7 commit 6002ba6
Show file tree
Hide file tree
Showing 15 changed files with 1,326 additions and 66 deletions.
69 changes: 69 additions & 0 deletions examples/dgmulti_2d/elixir_advection_diffusion_nonperiodic.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,69 @@
using Trixi, OrdinaryDiffEq

dg = DGMulti(polydeg = 3, element_type = Quad(), approximation_type = Polynomial(),
surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs),
volume_integral = VolumeIntegralWeakForm())

get_diffusivity() = 5.0e-2

equations = LinearScalarAdvectionEquation2D(1.0, 0.0)
equations_parabolic = LaplaceDiffusion2D(get_diffusivity(), equations)

# from "Robust DPG methods for transient convection-diffusion."
# Building bridges: connections and challenges in modern approaches to numerical partial differential equations.
# Springer, Cham, 2016. 179-203. Ellis, Truman, Jesse Chan, and Leszek Demkowicz."
function initial_condition_erikkson_johnson(x, t, equations)
l = 4
epsilon = get_diffusivity() # TODO: this requires epsilon < .6 due to the sqrt
lambda_1 = (-1 + sqrt(1 - 4 * epsilon * l)) / (-2 * epsilon)
lambda_2 = (-1 - sqrt(1 - 4 * epsilon * l)) / (-2 * epsilon)
r1 = (1 + sqrt(1 + 4 * pi^2 * epsilon^2)) / (2 * epsilon)
s1 = (1 - sqrt(1 + 4 * pi^2 * epsilon^2)) / (2 * epsilon)
u = exp(-l * t) * (exp(lambda_1 * x[1]) - exp(lambda_2 * x[1])) +
cos(pi * x[2]) * (exp(s1 * x[1]) - exp(r1 * x[1])) / (exp(-s1) - exp(-r1))
return SVector{1}(u)
end
initial_condition = initial_condition_erikkson_johnson

# tag different boundary segments
left(x, tol=50*eps()) = abs(x[1] + 1) < tol
right(x, tol=50*eps()) = abs(x[1]) < tol
bottom(x, tol=50*eps()) = abs(x[2] + .5) < tol
top(x, tol=50*eps()) = abs(x[2] - .5) < tol
entire_boundary(x, tol=50*eps()) = true
is_on_boundary = Dict(:left => left, :right => right, :top => top, :bottom => bottom,
:entire_boundary => entire_boundary)
mesh = DGMultiMesh(dg; coordinates_min=(-1.0, -0.5), coordinates_max=(0.0, 0.5),
cells_per_dimension=(16, 16), is_on_boundary)

# BC types
boundary_condition = BoundaryConditionDirichlet(initial_condition)

# define inviscid boundary conditions
boundary_conditions = (; :left => boundary_condition,
:top => boundary_condition,
:bottom => boundary_condition,
:right => boundary_condition_do_nothing)

# define viscous boundary conditions
boundary_conditions_parabolic = (; :entire_boundary => boundary_condition)

semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, dg;
boundary_conditions=(boundary_conditions, boundary_conditions_parabolic))

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

summary_callback = SummaryCallback()
alive_callback = AliveCallback(alive_interval=10)
analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg))
callbacks = CallbackSet(summary_callback, alive_callback)

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

time_int_tol = 1e-8
sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol,
save_everystep=false, callback=callbacks)
summary_callback() # print the timer summary
206 changes: 206 additions & 0 deletions examples/dgmulti_2d/elixir_navier_stokes_convergence.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,206 @@
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the ideal compressible Navier-Stokes equations

get_Re() = 100
get_Pr() = 0.72

equations = CompressibleEulerEquations2D(1.4)
# Note: If you change the Navier-Stokes parameters here, also change them in the initial condition
# I really do not like this structure but it should work for now
equations_parabolic = CompressibleNavierStokesEquations2D(equations, Reynolds=get_Re(), Prandtl=get_Pr(),
Mach_freestream=0.5, kappa=1.0)

# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial(),
surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs),
volume_integral = VolumeIntegralWeakForm())

top_bottom(x, tol=50*eps()) = abs(abs(x[2]) - 1) < tol
is_on_boundary = Dict(:top_bottom => top_bottom)
mesh = DGMultiMesh(dg, cells_per_dimension=(16, 16); periodicity=(true, false), is_on_boundary)

# Define initial condition
# Note: If you change the parameters here, also change it in the corresponding source terms
function initial_condition_navier_stokes_convergence_test(x, t, equations)

# Amplitude and shift
A = 0.5
c = 2.0

# convenience values for trig. functions
pi_x = pi * x[1]
pi_y = pi * x[2]
pi_t = pi * t

rho = c + A * sin(pi_x) * cos(pi_y) * cos(pi_t)
v1 = sin(pi_x) * log(x[2] + 2.0) * (1.0 - exp(-A * (x[2] - 1.0)) ) * cos(pi_t)
v2 = v1
p = rho^2

return prim2cons(SVector(rho, v1, v2, p), equations)
end

initial_condition = initial_condition_navier_stokes_convergence_test

@inline function source_terms_navier_stokes_convergence_test(u, x, t, equations)

y = x[2]

# we currently need to hardcode these parameters until we fix the "combined equation" issue
# see also https://github.com/trixi-framework/Trixi.jl/pull/1160
kappa = 1
inv_gamma_minus_one = inv(equations.gamma - 1)
Pr = get_Pr()
Re = get_Re()

# Same settings as in `initial_condition`
# Amplitude and shift
A = 0.5
c = 2.0

# convenience values for trig. functions
pi_x = pi * x[1]
pi_y = pi * x[2]
pi_t = pi * t

# compute the manufactured solution and all necessary derivatives
rho = c + A * sin(pi_x) * cos(pi_y) * cos(pi_t)
rho_t = -pi * A * sin(pi_x) * cos(pi_y) * sin(pi_t)
rho_x = pi * A * cos(pi_x) * cos(pi_y) * cos(pi_t)
rho_y = -pi * A * sin(pi_x) * sin(pi_y) * cos(pi_t)
rho_xx = -pi * pi * A * sin(pi_x) * cos(pi_y) * cos(pi_t)
rho_yy = -pi * pi * A * sin(pi_x) * cos(pi_y) * cos(pi_t)

v1 = sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t)
v1_t = -pi * sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * sin(pi_t)
v1_x = pi * cos(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t)
v1_y = sin(pi_x) * (A * log(y + 2.0) * exp(-A * (y - 1.0)) + (1.0 - exp(-A * (y - 1.0))) / (y + 2.0)) * cos(pi_t)
v1_xx = -pi * pi * sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t)
v1_xy = pi * cos(pi_x) * (A * log(y + 2.0) * exp(-A * (y - 1.0)) + (1.0 - exp(-A * (y - 1.0))) / (y + 2.0)) * cos(pi_t)
v1_yy = (sin(pi_x) * ( 2.0 * A * exp(-A * (y - 1.0)) / (y + 2.0)
- A * A * log(y + 2.0) * exp(-A * (y - 1.0))
- (1.0 - exp(-A * (y - 1.0))) / ((y + 2.0) * (y + 2.0))) * cos(pi_t))
v2 = v1
v2_t = v1_t
v2_x = v1_x
v2_y = v1_y
v2_xx = v1_xx
v2_xy = v1_xy
v2_yy = v1_yy

p = rho * rho
p_t = 2.0 * rho * rho_t
p_x = 2.0 * rho * rho_x
p_y = 2.0 * rho * rho_y
p_xx = 2.0 * rho * rho_xx + 2.0 * rho_x * rho_x
p_yy = 2.0 * rho * rho_yy + 2.0 * rho_y * rho_y

# Note this simplifies slightly because the ansatz assumes that v1 = v2
E = p * inv_gamma_minus_one + 0.5 * rho * (v1^2 + v2^2)
E_t = p_t * inv_gamma_minus_one + rho_t * v1^2 + 2.0 * rho * v1 * v1_t
E_x = p_x * inv_gamma_minus_one + rho_x * v1^2 + 2.0 * rho * v1 * v1_x
E_y = p_y * inv_gamma_minus_one + rho_y * v1^2 + 2.0 * rho * v1 * v1_y

# Some convenience constants
T_const = equations.gamma * inv_gamma_minus_one * kappa / Pr
inv_Re = 1.0 / Re
inv_rho_cubed = 1.0 / (rho^3)

# compute the source terms
# density equation
du1 = rho_t + rho_x * v1 + rho * v1_x + rho_y * v2 + rho * v2_y

# x-momentum equation
du2 = ( rho_t * v1 + rho * v1_t + p_x + rho_x * v1^2
+ 2.0 * rho * v1 * v1_x
+ rho_y * v1 * v2
+ rho * v1_y * v2
+ rho * v1 * v2_y
# stress tensor from x-direction
- 4.0 / 3.0 * v1_xx * inv_Re
+ 2.0 / 3.0 * v2_xy * inv_Re
- v1_yy * inv_Re
- v2_xy * inv_Re )
# y-momentum equation
du3 = ( rho_t * v2 + rho * v2_t + p_y + rho_x * v1 * v2
+ rho * v1_x * v2
+ rho * v1 * v2_x
+ rho_y * v2^2
+ 2.0 * rho * v2 * v2_y
# stress tensor from y-direction
- v1_xy * inv_Re
- v2_xx * inv_Re
- 4.0 / 3.0 * v2_yy * inv_Re
+ 2.0 / 3.0 * v1_xy * inv_Re )
# total energy equation
du4 = ( E_t + v1_x * (E + p) + v1 * (E_x + p_x)
+ v2_y * (E + p) + v2 * (E_y + p_y)
# stress tensor and temperature gradient terms from x-direction
- 4.0 / 3.0 * v1_xx * v1 * inv_Re
+ 2.0 / 3.0 * v2_xy * v1 * inv_Re
- 4.0 / 3.0 * v1_x * v1_x * inv_Re
+ 2.0 / 3.0 * v2_y * v1_x * inv_Re
- v1_xy * v2 * inv_Re
- v2_xx * v2 * inv_Re
- v1_y * v2_x * inv_Re
- v2_x * v2_x * inv_Re
- T_const * inv_rho_cubed * ( p_xx * rho * rho
- 2.0 * p_x * rho * rho_x
+ 2.0 * p * rho_x * rho_x
- p * rho * rho_xx ) * inv_Re
# stress tensor and temperature gradient terms from y-direction
- v1_yy * v1 * inv_Re
- v2_xy * v1 * inv_Re
- v1_y * v1_y * inv_Re
- v2_x * v1_y * inv_Re
- 4.0 / 3.0 * v2_yy * v2 * inv_Re
+ 2.0 / 3.0 * v1_xy * v2 * inv_Re
- 4.0 / 3.0 * v2_y * v2_y * inv_Re
+ 2.0 / 3.0 * v1_x * v2_y * inv_Re
- T_const * inv_rho_cubed * ( p_yy * rho * rho
- 2.0 * p_y * rho * rho_y
+ 2.0 * p * rho_y * rho_y
- p * rho * rho_yy ) * inv_Re )

return SVector(du1, du2, du3, du4)
end

# BC types
velocity_bc_top_bottom = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x, t, equations)[2:3])
heat_bc_top_bottom = Adiabatic((x, t, equations) -> 0.0)
boundary_condition_top_bottom = BoundaryConditionViscousWall(velocity_bc_top_bottom, heat_bc_top_bottom)

# define inviscid boundary conditions
boundary_conditions = (; :top_bottom => boundary_condition_slip_wall)

# define viscous boundary conditions
boundary_conditions_parabolic = (; :top_bottom => boundary_condition_top_bottom)

semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, dg;
boundary_conditions=(boundary_conditions, boundary_conditions_parabolic),
source_terms=source_terms_navier_stokes_convergence_test)


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

# Create ODE problem with time span `tspan`
tspan = (0.0, 0.5)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()
alive_callback = AliveCallback(alive_interval=10)
analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg))
callbacks = CallbackSet(summary_callback, alive_callback)

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

time_int_tol = 1e-8
sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol,
save_everystep=false, callback=callbacks)
summary_callback() # print the timer summary
70 changes: 70 additions & 0 deletions examples/dgmulti_2d/elixir_navier_stokes_lid_driven_cavity.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,70 @@
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the ideal compressible Navier-Stokes equations

equations = CompressibleEulerEquations2D(1.4)
# Note: If you change the Navier-Stokes parameters here, also change them in the initial condition
# I really do not like this structure but it should work for now
equations_parabolic = CompressibleNavierStokesEquations2D(equations, Reynolds=1000, Prandtl=0.72,
Mach_freestream=0.1, kappa=1.0)

# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
dg = DGMulti(polydeg = 3, element_type = Quad(), approximation_type = GaussSBP(),
surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs),
volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha))

top(x, tol=50*eps()) = abs(x[2] - 1) < tol
rest_of_boundary(x, tol=50*eps()) = !top(x, tol)
is_on_boundary = Dict(:top => top, :rest_of_boundary => rest_of_boundary)
mesh = DGMultiMesh(dg, cells_per_dimension=(16, 16); is_on_boundary)

function initial_condition_cavity(x, t, equations::CompressibleEulerEquations2D)
Ma = 0.1
rho = 1.0
u, v = 0, 0
p = 1.0 / (Ma^2 * equations.gamma)
return prim2cons(SVector(rho, u, v, p), equations)
end
initial_condition = initial_condition_cavity

# BC types
velocity_bc_lid = NoSlip((x, t, equations) -> SVector(1.0, 0.0))
velocity_bc_cavity = NoSlip((x, t, equations) -> SVector(0.0, 0.0))
heat_bc = Adiabatic((x, t, equations) -> 0.0)
boundary_condition_lid = BoundaryConditionViscousWall(velocity_bc_lid, heat_bc)
boundary_condition_cavity = BoundaryConditionViscousWall(velocity_bc_cavity, heat_bc)

# define inviscid boundary conditions
boundary_conditions = (; :top => boundary_condition_slip_wall,
:rest_of_boundary => boundary_condition_slip_wall)

# define viscous boundary conditions
boundary_conditions_parabolic = (; :top => boundary_condition_lid,
:rest_of_boundary => boundary_condition_cavity)

semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, dg;
boundary_conditions=(boundary_conditions, boundary_conditions_parabolic))


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

# Create ODE problem with time span `tspan`
tspan = (0.0, 10.0)
ode = semidiscretize(semi, tspan);

summary_callback = SummaryCallback()
alive_callback = AliveCallback(alive_interval=10)
analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg))
callbacks = CallbackSet(summary_callback, alive_callback)

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

time_int_tol = 1e-8
sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol,
save_everystep=false, callback=callbacks)
summary_callback() # print the timer summary
2 changes: 1 addition & 1 deletion examples/tree_2d_dgsem/elixir_advection_diffusion.jl
Original file line number Diff line number Diff line change
Expand Up @@ -77,7 +77,7 @@ callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback)
# run the simulation

# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
time_int_tol = 1.0e-8
time_int_tol = 1.0e-11
sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol,
save_everystep=false, callback=callbacks)

Expand Down
Loading

0 comments on commit 6002ba6

Please sign in to comment.