Skip to content

Commit

Permalink
Add parabolic BCs for P4estMesh{2} (#1493)
Browse files Browse the repository at this point in the history 
* generalize function signatures to P4estMesh

* add specializations for P4estMesh


d

* add normals

* add surface integrals

* fix type ambiguity

* generalizing `apply_jacobian!` to P4estMesh

* resolving type ambiguity with apply_jacobian!


d

* `apply_jacobian!` -> `apply_jacobian_parabolic!`

* `apply_jacobian!` -> `apply_jacobian_parabolic!`

* switch to `apply_jacobian_parabolic!`

* Update src/solvers/dgsem_tree/dg_1d_parabolic.jl

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

* missed one

* draft of prolong2interfaces and calc_interface_flux

* cache -> cache_parabolic

* adding prolong2boundaries! and calc_boundary_flux_gradients! back

* remove todo

* variable renaming

* extending TreeMesh parabolic functions to P4estMesh

* adding elixir

* comments

* add prolong2boundaries! (untested)

* update test

* initial commit

* fix CI

f

* Update src/solvers/dgsem_p4est/dg_2d_parabolic.jl

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

* Update src/solvers/dgsem_p4est/dg_2d_parabolic.jl

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

* add "no mortars" check

* add curved elixir

* fix gradient bug

* add curved test

* Apply suggestions from code review

Co-authored-by: Erik Faulhaber <[email protected]>
Co-authored-by: Michael Schlottke-Lakemper <[email protected]>

* add comment on mapping

* reuse P4estMesh{2} code

* fix += for muladd

* Update examples/p4est_2d_dgsem/elixir_advection_diffusion_periodic_curved.jl

Co-authored-by: Erik Faulhaber <[email protected]>

* comment

* comments + remove cruft

* add BCs for parabolic P43st

* add tests

* Update examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl

* formatting

* fix CNS convergence elixir and add to tests

* update test values

---------

Co-authored-by: Hendrik Ranocha <[email protected]>
Co-authored-by: Erik Faulhaber <[email protected]>
Co-authored-by: Michael Schlottke-Lakemper <[email protected]>
  • Loading branch information
@jlchan @ranocha
@efaulhaber @sloede
4 people authored Jun 23, 2023
1 parent 1b69182 commit b7be585
Show file tree
Hide file tree
Showing 5 changed files with 544 additions and 8 deletions.
96 changes: 96 additions & 0 deletions examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_curved.jl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,96 @@
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the linear advection-diffusion equation

diffusivity() = 5.0e-2
advection_velocity = (1.0, 0.0)
equations = LinearScalarAdvectionEquation2D(advection_velocity)
equations_parabolic = LaplaceDiffusion2D(diffusivity(), equations)

# Example setup taken from
# - Truman Ellis, Jesse Chan, and Leszek Demkowicz (2016).
# Robust DPG methods for transient convection-diffusion.
# In: Building bridges: connections and challenges in modern approaches
# to numerical partial differential equations.
# [DOI](https://doi.org/10.1007/978-3-319-41640-3_6).
function initial_condition_eriksson_johnson(x, t, equations)
l = 4
epsilon = diffusivity() # TODO: this requires epsilon < .6 due to 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_eriksson_johnson

boundary_conditions = Dict(:x_neg => BoundaryConditionDirichlet(initial_condition),
:y_neg => BoundaryConditionDirichlet(initial_condition),
:y_pos => BoundaryConditionDirichlet(initial_condition),
:x_pos => boundary_condition_do_nothing)

boundary_conditions_parabolic = Dict(:x_neg => BoundaryConditionDirichlet(initial_condition),
:x_pos => BoundaryConditionDirichlet(initial_condition),
:y_neg => BoundaryConditionDirichlet(initial_condition),
:y_pos => BoundaryConditionDirichlet(initial_condition))

# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)

coordinates_min = (-1.0, -0.5)
coordinates_max = ( 0.0, 0.5)

# This maps the domain [-1, 1]^2 to [-1, 0] x [-0.5, 0.5] while also
# introducing a curved warping to interior nodes.
function mapping(xi, eta)
x = xi + 0.1 * sin(pi * xi) * sin(pi * eta)
y = eta + 0.1 * sin(pi * xi) * sin(pi * eta)
return SVector(0.5 * (1 + x) - 1, 0.5 * y)
end

trees_per_dimension = (4, 4)
mesh = P4estMesh(trees_per_dimension,
polydeg=3, initial_refinement_level=2,
mapping=mapping, periodicity=(false, false))

# A semidiscretization collects data structures and functions for the spatial discretization
semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver,
boundary_conditions = (boundary_conditions, boundary_conditions_parabolic))


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

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

# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
# and resets the timers
summary_callback = SummaryCallback()

# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

# The AliveCallback prints short status information in regular intervals
alive_callback = AliveCallback(analysis_interval=analysis_interval)

# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
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-11
sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol,
ode_default_options()..., callback=callbacks)

# Print the timer summary
summary_callback()
209 changes: 209 additions & 0 deletions examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,209 @@
using OrdinaryDiffEq
using Trixi

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

prandtl_number() = 0.72
mu() = 0.01

equations = CompressibleEulerEquations2D(1.4)
equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu=mu(), Prandtl=prandtl_number(),
gradient_variables=GradientVariablesPrimitive())

# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs,
volume_integral=VolumeIntegralWeakForm())

coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y))
coordinates_max = ( 1.0, 1.0) # maximum coordinates (max(x), max(y))

trees_per_dimension = (4, 4)
mesh = P4estMesh(trees_per_dimension,
polydeg=3, initial_refinement_level=2,
coordinates_min=coordinates_min, coordinates_max=coordinates_max,
periodicity=(true, false))

# Note: the initial condition cannot be specialized to `CompressibleNavierStokesDiffusion2D`
# since it is called by both the parabolic solver (which passes in `CompressibleNavierStokesDiffusion2D`)
# and by the initial condition (which passes in `CompressibleEulerEquations2D`).
# This convergence test setup was originally derived by Andrew Winters (@andrewwinters5000)
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

@inline function source_terms_navier_stokes_convergence_test(u, x, t, equations)
y = x[2]

# TODO: parabolic
# 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
inv_gamma_minus_one = inv(equations.gamma - 1)
Pr = prandtl_number()
mu_ = mu()

# 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 / Pr
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 * mu_
+ 2.0 / 3.0 * v2_xy * mu_
- v1_yy * mu_
- v2_xy * mu_ )
# 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 * mu_
- v2_xx * mu_
- 4.0 / 3.0 * v2_yy * mu_
+ 2.0 / 3.0 * v1_xy * mu_ )
# 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 * mu_
+ 2.0 / 3.0 * v2_xy * v1 * mu_
- 4.0 / 3.0 * v1_x * v1_x * mu_
+ 2.0 / 3.0 * v2_y * v1_x * mu_
- v1_xy * v2 * mu_
- v2_xx * v2 * mu_
- v1_y * v2_x * mu_
- v2_x * v2_x * mu_
- 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 ) * mu_
# stress tensor and temperature gradient terms from y-direction
- v1_yy * v1 * mu_
- v2_xy * v1 * mu_
- v1_y * v1_y * mu_
- v2_x * v1_y * mu_
- 4.0 / 3.0 * v2_yy * v2 * mu_
+ 2.0 / 3.0 * v1_xy * v2 * mu_
- 4.0 / 3.0 * v2_y * v2_y * mu_
+ 2.0 / 3.0 * v1_x * v2_y * mu_
- 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 ) * mu_ )

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

initial_condition = initial_condition_navier_stokes_convergence_test

# 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 = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom, heat_bc_top_bottom)

# define inviscid boundary conditions
boundary_conditions = Dict(:y_neg => boundary_condition_slip_wall,
:y_pos => boundary_condition_slip_wall)

# define viscous boundary conditions
boundary_conditions_parabolic = Dict(:y_neg => boundary_condition_top_bottom,
:y_pos => boundary_condition_top_bottom)

semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver;
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)
callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback)

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

time_int_tol = 1e-8
sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol, dt = 1e-5,
ode_default_options()..., callback=callbacks)
summary_callback() # print the timer summary

82 changes: 82 additions & 0 deletions examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,82 @@
using OrdinaryDiffEq
using Trixi

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

# TODO: parabolic; unify names of these accessor functions
prandtl_number() = 0.72
mu() = 0.001

equations = CompressibleEulerEquations2D(1.4)
equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu=mu(),
Prandtl=prandtl_number())

# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)

coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y))
coordinates_max = ( 1.0, 1.0) # maximum coordinates (max(x), max(y))

# Create a uniformly refined mesh
trees_per_dimension = (4, 4)
mesh = P4estMesh(trees_per_dimension,
polydeg=3, initial_refinement_level=2,
coordinates_min=coordinates_min, coordinates_max=coordinates_max,
periodicity=(false, false))

function initial_condition_cavity(x, t, equations::CompressibleEulerEquations2D)
Ma = 0.1
rho = 1.0
u, v = 0.0, 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 = BoundaryConditionNavierStokesWall(velocity_bc_lid, heat_bc)
boundary_condition_cavity = BoundaryConditionNavierStokesWall(velocity_bc_cavity, heat_bc)

# define periodic boundary conditions everywhere
boundary_conditions = Dict( :x_neg => boundary_condition_slip_wall,
:y_neg => boundary_condition_slip_wall,
:y_pos => boundary_condition_slip_wall,
:x_pos => boundary_condition_slip_wall)

boundary_conditions_parabolic = Dict( :x_neg => boundary_condition_cavity,
:y_neg => boundary_condition_cavity,
:y_pos => boundary_condition_lid,
:x_pos => boundary_condition_cavity)

# A semidiscretization collects data structures and functions for the spatial discretization
semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
initial_condition, solver;
boundary_conditions=(boundary_conditions,
boundary_conditions_parabolic))

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

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

summary_callback = SummaryCallback()
alive_callback = AliveCallback(alive_interval=100)
analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
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,
ode_default_options()..., callback=callbacks)
summary_callback() # print the timer summary


Loading

0 comments on commit b7be585

Please sign in to comment.