Skip to content

Commit

Permalink
Merge branch 'main' into subcell-limiting-minmax
Browse files Browse the repository at this point in the history
  • Loading branch information
bennibolm committed Oct 12, 2023
2 parents 3e6f7f8 + 7fd4503 commit d7dcb75
Show file tree
Hide file tree
Showing 34 changed files with 429 additions and 455 deletions.
2 changes: 1 addition & 1 deletion .github/workflows/CacheNotebooks.yml
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@ jobs:
steps:
# Checks-out your repository under $GITHUB_WORKSPACE, so your job can access it
- name: Checkout caching branch
uses: actions/checkout@v3
uses: actions/checkout@v4
with:
ref: tutorial_notebooks

Expand Down
2 changes: 1 addition & 1 deletion .github/workflows/DocPreviewCleanup.yml
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@ jobs:
runs-on: ubuntu-latest
steps:
- name: Checkout gh-pages branch
uses: actions/checkout@v3
uses: actions/checkout@v4
with:
ref: gh-pages

Expand Down
2 changes: 1 addition & 1 deletion .github/workflows/Documenter.yml
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@ jobs:
build:
runs-on: ubuntu-latest
steps:
- uses: actions/checkout@v3
- uses: actions/checkout@v4
- uses: julia-actions/setup-julia@v1
with:
version: '1.9'
Expand Down
2 changes: 1 addition & 1 deletion .github/workflows/FormatCheck.yml
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ jobs:
with:
version: ${{ matrix.julia-version }}

- uses: actions/checkout@v3
- uses: actions/checkout@v4
- name: Install JuliaFormatter and format
# This will use the latest version by default but you can set the version like so:
#
Expand Down
4 changes: 2 additions & 2 deletions .github/workflows/Invalidations.yml
Original file line number Diff line number Diff line change
Expand Up @@ -19,12 +19,12 @@ jobs:
- uses: julia-actions/setup-julia@v1
with:
version: '1'
- uses: actions/checkout@v3
- uses: actions/checkout@v4
- uses: julia-actions/julia-buildpkg@v1
- uses: julia-actions/julia-invalidations@v1
id: invs_pr

- uses: actions/checkout@v3
- uses: actions/checkout@v4
with:
ref: ${{ github.event.repository.default_branch }}
- uses: julia-actions/julia-buildpkg@v1
Expand Down
2 changes: 1 addition & 1 deletion .github/workflows/ReviewChecklist.yml
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@ jobs:
runs-on: ubuntu-latest
steps:
- name: Check out repository
uses: actions/checkout@v3
uses: actions/checkout@v4
- name: Add review checklist
uses: trixi-framework/add-pr-review-checklist@v1
with:
Expand Down
4 changes: 2 additions & 2 deletions .github/workflows/SpellCheck.yml
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,6 @@ jobs:
runs-on: ubuntu-latest
steps:
- name: Checkout Actions Repository
uses: actions/checkout@v3
uses: actions/checkout@v4
- name: Check spelling
uses: crate-ci/[email protected].9
uses: crate-ci/[email protected].15
2 changes: 1 addition & 1 deletion .github/workflows/benchmark.yml
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ jobs:
arch:
- x64
steps:
- uses: actions/checkout@v3
- uses: actions/checkout@v4
with:
fetch-depth: 0
- run: |
Expand Down
4 changes: 2 additions & 2 deletions .github/workflows/ci.yml
Original file line number Diff line number Diff line change
Expand Up @@ -101,7 +101,7 @@ jobs:
arch: x64
trixi_test: threaded
steps:
- uses: actions/checkout@v3
- uses: actions/checkout@v4
- uses: julia-actions/setup-julia@v1
with:
version: ${{ matrix.version }}
Expand Down Expand Up @@ -175,7 +175,7 @@ jobs:
# Instead, we use the more tedious approach described above.
# At first, we check out the repository and download all artifacts
# (and list files for debugging).
- uses: actions/checkout@v3
- uses: actions/checkout@v4
- uses: actions/download-artifact@v3
- run: ls -R
# Next, we merge the individual coverage files and upload
Expand Down
2 changes: 1 addition & 1 deletion Project.toml
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "Trixi"
uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
authors = ["Michael Schlottke-Lakemper <[email protected]>", "Gregor Gassner <[email protected]>", "Hendrik Ranocha <[email protected]>", "Andrew R. Winters <[email protected]>", "Jesse Chan <[email protected]>"]
version = "0.5.45-pre"
version = "0.5.46-pre"

[deps]
CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
Expand Down
5 changes: 5 additions & 0 deletions docs/src/callbacks.md
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,11 @@ An example elixir using AMR can be found at [`examples/tree_2d_dgsem/elixir_adve
The [`AnalysisCallback`](@ref) can be used to analyze the numerical solution, e.g. calculate
errors or user-specified integrals, and print the results to the screen. The results can also be
saved in a file. An example can be found at [`examples/tree_2d_dgsem/elixir_euler_vortex.jl`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_2d_dgsem/elixir_euler_vortex.jl).
Note that the errors (e.g. `L2 error` or `Linf error`) are computed with respect to the initial condition.
The percentage of the simulation time refers to the ratio of the current time and the final time, i.e. it does
not consider the maximal number of iterations. So the simulation could finish before 100% are reached.
Note that, e.g., due to AMR or smaller time step sizes, the simulation can actually take longer than
the percentage indicates.
In [Performance metrics of the `AnalysisCallback`](@ref performance-metrics) you can find a detailed
description of the different performance metrics the `AnalysisCallback` computes.

Expand Down
5 changes: 3 additions & 2 deletions docs/src/overview.md
Original file line number Diff line number Diff line change
Expand Up @@ -55,14 +55,15 @@ different features on different mesh types.
| Element type | line, square, cube | line, quadᵃ, hexᵃ | quadᵃ | quadᵃ, hexᵃ | simplex, quadᵃ, hexᵃ |
| Adaptive mesh refinement | ✅ | ❌ | ❌ | ✅ | ❌ | [`AMRCallback`](@ref)
| Solver type | [`DGSEM`](@ref) | [`DGSEM`](@ref) | [`DGSEM`](@ref) | [`DGSEM`](@ref) | [`DGMulti`](@ref) |
| Domain | hypercube | mapped hypercube | arbitrary | arbitrary | arbitraryᵇ |
| Domain | hypercube | mapped hypercube | arbitrary | arbitrary | arbitrary |
| Weak form | ✅ | ✅ | ✅ | ✅ | ✅ | [`VolumeIntegralWeakForm`](@ref)
| Flux differencing | ✅ | ✅ | ✅ | ✅ | ✅ | [`VolumeIntegralFluxDifferencing`](@ref)
| Shock capturing | ✅ | ✅ | ✅ | ✅ | ❌ | [`VolumeIntegralShockCapturingHG`](@ref)
| Nonconservative equations | ✅ | ✅ | ✅ | ✅ | ✅ | e.g., GLM MHD or shallow water equations
| Parabolic termsᵇ | ✅ | ✅ | ❌ | ✅ | ✅ | e.g., [`CompressibleNavierStokesDiffusion2D`](@ref)

ᵃ: quad = quadrilateral, hex = hexahedron
ᵇ: curved meshes supported for `SBP` and `GaussSBP` approximation types for `VolumeIntegralFluxDifferencing` solvers on quadrilateral and hexahedral `DGMultiMesh`es (non-conservative terms not yet supported)
ᵇ: Parabolic terms do not currently support adaptivity.

## Time integration methods

Expand Down
5 changes: 5 additions & 0 deletions docs/src/parallelization.md
Original file line number Diff line number Diff line change
Expand Up @@ -186,6 +186,11 @@ julia> set_preferences!(
"libhdf5" => "/path/to/your/libhdf5.so",
"libhdf5_hl" => "/path/to/your/libhdf5_hl.so", force = true)
```
Alternatively, with HDF5.jl v0.17.1 or higher you can use
```julia
julia> using HDF5
julia> HDF5.API.set_libraries!("/path/to/your/libhdf5.so", "/path/to/your/libhdf5_hl.so")
```
For more information see also the
[documentation of HDF5.jl](https://juliaio.github.io/HDF5.jl/stable/mpi/). In total, you should
have a file called LocalPreferences.toml in the project directory that contains a section
Expand Down
215 changes: 215 additions & 0 deletions examples/p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,215 @@
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=(false, 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)

boundary_condition_left_right = BoundaryConditionDirichlet(initial_condition_navier_stokes_convergence_test)

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

# define viscous boundary conditions
boundary_conditions_parabolic = Dict(:x_neg => boundary_condition_left_right,
:x_pos => boundary_condition_left_right,
: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

Original file line number Diff line number Diff line change
Expand Up @@ -66,7 +66,7 @@ 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)
callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback)

###############################################################################
# run the simulation
Expand Down
10 changes: 8 additions & 2 deletions src/callbacks_step/alive.jl
Original file line number Diff line number Diff line change
Expand Up @@ -86,9 +86,15 @@ function (alive_callback::AliveCallback)(integrator)
println(""^100)
println()
elseif mpi_isroot()
t = integrator.t
t_initial = first(integrator.sol.prob.tspan)
t_final = last(integrator.sol.prob.tspan)
sim_time_percentage = (t - t_initial) / (t_final - t_initial) * 100
runtime_absolute = 1.0e-9 * (time_ns() - alive_callback.start_time)
@printf("#timesteps: %6d │ Δt: %.4e │ sim. time: %.4e │ run time: %.4e s\n",
integrator.stats.naccept, integrator.dt, integrator.t, runtime_absolute)
println(rpad(@sprintf("#timesteps: %6d │ Δt: %.4e │ sim. time: %.4e (%5.3f%%)",
integrator.stats.naccept, integrator.dt, t,
sim_time_percentage), 71) *
@sprintf("│ run time: %.4e s", runtime_absolute))
end

# avoid re-evaluating possible FSAL stages
Expand Down
Loading

0 comments on commit d7dcb75

Please sign in to comment.