Navier-Stokes in 2D on DGMulti and TreeMesh

.github/workflows/ci.yml

@@ -69,6 +69,7 @@ jobs:
           - p4est_part1
           - p4est_part2
           - unstructured_dgmulti
+          - parabolic
           - paper_self_gravitating_gas_dynamics
           - misc_part1
           - misc_part2

examples/dgmulti_2d/elixir_advection_diffusion.jl

@@ -0,0 +1,56 @@
+using Trixi, OrdinaryDiffEq
+
+dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial(),
+             surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs),
+             volume_integral = VolumeIntegralWeakForm())
+
+equations = LinearScalarAdvectionEquation2D(1.5, 1.0)
+equations_parabolic = LaplaceDiffusion2D(5.0e-2, equations)
+
+initial_condition_zero(x, t, equations::LinearScalarAdvectionEquation2D) = SVector(0.0)
+initial_condition = initial_condition_zero
+
+# tag different boundary segments
+left(x, tol=50*eps()) = abs(x[1] + 1) < tol
+right(x, tol=50*eps()) = abs(x[1] - 1) < tol
+bottom(x, tol=50*eps()) = abs(x[2] + 1) < tol
+top(x, tol=50*eps()) = abs(x[2] - 1) < tol
+is_on_boundary = Dict(:left => left, :right => right, :top => top, :bottom => bottom)
+mesh = DGMultiMesh(dg, cells_per_dimension=(16, 16); is_on_boundary)
+
+# BC types
+boundary_condition_left = BoundaryConditionDirichlet((x, t, equations) -> SVector(1 + 0.1 * x[2]))
+boundary_condition_zero = BoundaryConditionDirichlet((x, t, equations) -> SVector(0.0))
+boundary_condition_neumann_zero = BoundaryConditionNeumann((x, t, equations) -> SVector(0.0))
+
+# define inviscid boundary conditions
+boundary_conditions = (; :left => boundary_condition_left,
+                         :bottom => boundary_condition_zero,
+                         :top => boundary_condition_do_nothing,
+                         :right => boundary_condition_do_nothing)
+
+# define viscous boundary conditions
+boundary_conditions_parabolic = (; :left => boundary_condition_left,
+                                   :bottom => boundary_condition_zero,
+                                   :top => boundary_condition_zero,
+                                   :right => boundary_condition_neumann_zero)
+
+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-6
+sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol,
+            save_everystep=false, callback=callbacks)
+summary_callback() # print the timer summary

examples/dgmulti_2d/elixir_advection_diffusion_nonperiodic.jl

@@ -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

examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl

@@ -0,0 +1,35 @@
+using Trixi, OrdinaryDiffEq
+
+dg = DGMulti(polydeg = 1, element_type = Tri(), approximation_type = Polynomial(),
+             surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs),
+             volume_integral = VolumeIntegralWeakForm())
+
+equations = LinearScalarAdvectionEquation2D(0.0, 0.0)
+equations_parabolic = LaplaceDiffusion2D(5.0e-1, equations)
+
+function initial_condition_sharp_gaussian(x, t, equations::LinearScalarAdvectionEquation2D)
+  return SVector(exp(-100 * (x[1]^2 + x[2]^2)))
+end
+initial_condition = initial_condition_sharp_gaussian
+
+mesh = DGMultiMesh(dg, cells_per_dimension = (16, 16), periodicity=true)
+semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
+                                             initial_condition, dg)
+
+tspan = (0.0, 0.1)
+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, analysis_callback)
+
+###############################################################################
+# run the simulation
+
+time_int_tol = 1e-6
+sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol,
+            dt = time_int_tol, save_everystep=false, callback=callbacks)
+
+summary_callback() # print the timer summary

examples/dgmulti_2d/elixir_navier_stokes_convergence.jl

@@ -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