diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml index bc324c689b..bb5a32f72e 100644 --- a/.github/workflows/SpellCheck.yml +++ b/.github/workflows/SpellCheck.yml @@ -10,4 +10,4 @@ jobs: - name: Checkout Actions Repository uses: actions/checkout@v3 - name: Check spelling - uses: crate-ci/typos@v1.15.0 + uses: crate-ci/typos@v1.16.0 diff --git a/NEWS.md b/NEWS.md index 9b46ba565f..35c7039b2e 100644 --- a/NEWS.md +++ b/NEWS.md @@ -9,6 +9,7 @@ for human readability. #### Added - Experimental support for 3D parabolic diffusion terms has been added. +- Capability to set truly discontinuous initial conditions in 1D. #### Changed diff --git a/Project.toml b/Project.toml index 9d51e4dcff..828f4778f7 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.5.29-pre" +version = "0.5.32-pre" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" diff --git a/docs/literate/src/files/adding_nonconservative_equation.jl b/docs/literate/src/files/adding_nonconservative_equation.jl index 08dd631058..110fa48607 100644 --- a/docs/literate/src/files/adding_nonconservative_equation.jl +++ b/docs/literate/src/files/adding_nonconservative_equation.jl @@ -147,7 +147,7 @@ plot(sol) # above. error_1 = analysis_callback(sol).l2 |> first -@test isapprox(error_1, 0.0002961027497) #src +@test isapprox(error_1, 0.00029609575838969394) #src # Next, we increase the grid resolution by one refinement level and run the # simulation again. diff --git a/docs/literate/src/files/shock_capturing.jl b/docs/literate/src/files/shock_capturing.jl index b165f7ec8b..afa34cbf06 100644 --- a/docs/literate/src/files/shock_capturing.jl +++ b/docs/literate/src/files/shock_capturing.jl @@ -48,7 +48,7 @@ # with the total energy $\mathbb{E}=\max\big(\frac{m_N^2}{\sum_{j=0}^N m_j^2}, \frac{m_{N-1}^2}{\sum_{j=0}^{N-1} m_j^2}\big)$, # threshold $\mathbb{T}= 0.5 * 10^{-1.8*(N+1)^{1/4}}$ and parameter $s=ln\big(\frac{1-0.0001}{0.0001}\big)\approx 9.21024$. -# For computational efficiency, $\alpha_{min}$ is introduced und used for +# For computational efficiency, $\alpha_{min}$ is introduced and used for # ```math # \tilde{\alpha} = \begin{cases} # 0, & \text{if } \alpha<\alpha_{min}\\ diff --git a/docs/src/github-git.md b/docs/src/github-git.md index ad5991d87a..57b63073e7 100644 --- a/docs/src/github-git.md +++ b/docs/src/github-git.md @@ -112,7 +112,7 @@ branch, and the corresponding pull request will be updated automatically. Please note that a review has nothing to do with the lack of experience of the person developing changes: We try to review all code before it gets added to `main`, even from the most experienced developers. This is good practice and -helps to keep the error rate low while ensuring the the code is developed in a +helps to keep the error rate low while ensuring that the code is developed in a consistent fashion. Furthermore, do not take criticism of your code personally - we just try to keep Trixi.jl as accessible and easy to use for everyone. @@ -121,7 +121,7 @@ Once your branch is reviewed and declared ready for merging by the reviewer, make sure that all the latest changes have been pushed. Then, one of the developers will merge your PR. If you are one of the developers, you can also go to the pull request page on GitHub and and click on **Merge pull request**. -Voilá, you are done! Your branch will have been merged to +Voilà, you are done! Your branch will have been merged to `main` and the source branch will have been deleted in the GitHub repository (if you are not working in your own fork). diff --git a/docs/src/styleguide.md b/docs/src/styleguide.md index de367c086c..60e227204c 100644 --- a/docs/src/styleguide.md +++ b/docs/src/styleguide.md @@ -53,9 +53,9 @@ of your PR), you need to install JuliaFormatter.jl first by running ```shell julia -e 'using Pkg; Pkg.add("JuliaFormatter")' ``` -You can then recursively format all Julia files in the Trixi.jl repo by executing +You can then recursively format the core Julia files in the Trixi.jl repo by executing ```shell -julia -e 'using JuliaFormatter; format(".") +julia -e 'using JuliaFormatter; format(["benchmark", "ext", "src", "utils"])' ``` from inside the Trixi.jl repository. For convenience, there is also a script you can directly run from your terminal shell, which will automatically install JuliaFormatter in a @@ -65,3 +65,14 @@ utils/trixi-format.jl ``` You can get more information about using the convenience script by running it with the `--help`/`-h` flag. + +### Checking formatting before committing +It can be convenient to check the formatting of source code automatically before each commit. +We use git-hooks for it and provide a `pre-commit` script in the `utils` folder. The script uses +[JuliaFormatter.jl](https://github.com/domluna/JuliaFormatter.jl) just like formatting script that +runs over the whole Trixi.jl directory. +You can copy the `pre-commit`-script into `.git/hooks/pre-commit` and it will check your formatting +before each commit. If errors are found the commit is aborted and you can add the corrections via +```shell +git add -p +``` diff --git a/examples/dgmulti_3d/elixir_advection_tensor_wedge.jl b/examples/dgmulti_3d/elixir_advection_tensor_wedge.jl new file mode 100644 index 0000000000..4f43f2571a --- /dev/null +++ b/examples/dgmulti_3d/elixir_advection_tensor_wedge.jl @@ -0,0 +1,56 @@ +using OrdinaryDiffEq +using Trixi +using LinearAlgebra + +############################################################################### +equations = LinearScalarAdvectionEquation3D(1.0, 1.0, 1.0) + +initial_condition = initial_condition_convergence_test + +# Define the polynomial degrees for the polynoms of the triangular base and the line +# of the tensor-prism +tensor_polydeg = (3, 4) + +dg = DGMulti(element_type = Wedge(), + approximation_type = Polynomial(), + surface_flux = flux_lax_friedrichs, + polydeg = tensor_polydeg) + + +cells_per_dimension = (8, 8, 8) +mesh = DGMultiMesh(dg, + cells_per_dimension, + coordinates_min = (-1.0, -1.0, -1.0), + coordinates_max = (1.0, 1.0, 1.0), + periodicity = true) + + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg, + boundary_conditions=boundary_condition_periodic) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 5.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg)) + +alive_callback = AliveCallback(analysis_interval=analysis_interval) + +# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step +stepsize_callback = StepsizeCallback(cfl=1.0) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, stepsize_callback) + + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), dt = 1.0, + save_everystep=false, callback=callbacks); + +summary_callback() # print the timer summary \ No newline at end of file diff --git a/examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_curved.jl b/examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_curved.jl new file mode 100644 index 0000000000..55682f73fc --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_curved.jl @@ -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() diff --git a/examples/p4est_2d_dgsem/elixir_euler_sedov.jl b/examples/p4est_2d_dgsem/elixir_euler_sedov.jl index 9f5247e8c4..d5d8e0c78b 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_sedov.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_sedov.jl @@ -11,7 +11,7 @@ equations = CompressibleEulerEquations2D(1.4) initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) The Sedov blast wave setup based on Flash -- http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 """ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) # Set up polar coordinates @@ -20,7 +20,7 @@ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEq y_norm = x[2] - inicenter[2] r = sqrt(x_norm^2 + y_norm^2) - # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 + # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) E = 1.0 p0_inner = 3 * (equations.gamma - 1) * E / (3 * pi * r0^2) diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl new file mode 100644 index 0000000000..8111df8251 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl @@ -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 + diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl new file mode 100644 index 0000000000..051f4defe5 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl @@ -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 + + diff --git a/examples/p4est_3d_dgsem/elixir_euler_sedov.jl b/examples/p4est_3d_dgsem/elixir_euler_sedov.jl index 00da413285..6fa285b556 100644 --- a/examples/p4est_3d_dgsem/elixir_euler_sedov.jl +++ b/examples/p4est_3d_dgsem/elixir_euler_sedov.jl @@ -11,7 +11,7 @@ equations = CompressibleEulerEquations3D(1.4) initial_condition_medium_sedov_blast_wave(x, t, equations::CompressibleEulerEquations3D) The Sedov blast wave setup based on Flash -- http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 with smaller strength of the initial discontinuity. """ function initial_condition_medium_sedov_blast_wave(x, t, equations::CompressibleEulerEquations3D) @@ -22,7 +22,7 @@ function initial_condition_medium_sedov_blast_wave(x, t, equations::Compressible z_norm = x[3] - inicenter[3] r = sqrt(x_norm^2 + y_norm^2 + z_norm^2) - # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 + # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) E = 1.0 p0_inner = 3 * (equations.gamma - 1) * E / (4 * pi * r0^2) diff --git a/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_jeans_instability.jl b/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_jeans_instability.jl index 1774e39513..fb445616cd 100644 --- a/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_jeans_instability.jl +++ b/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_jeans_instability.jl @@ -15,7 +15,7 @@ The classical Jeans instability taken from - Dominik Derigs, Andrew R. Winters, Gregor J. Gassner, Stefanie Walch (2016) A Novel High-Order, Entropy Stable, 3D AMR MHD Solver with Guaranteed Positive Pressure [arXiv: 1605.03572](https://arxiv.org/abs/1605.03572) -- Flash manual https://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel.pdf +- Flash manual https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node189.html#SECTION010131000000000000000 in CGS (centimeter, gram, second) units. """ function initial_condition_jeans_instability(x, t, @@ -32,7 +32,7 @@ function initial_condition_jeans_instability(x, t, pres0 = 1.5e7 # dyn/cm^2 delta0 = 1e-3 # set wave vector values for perturbation (units 1/cm) - # see FLASH manual: https://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel.pdf + # see FLASH manual: https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node189.html#SECTION010131000000000000000 kx = 2.0*pi/0.5 # 2π/λ_x, λ_x = 0.5 ky = 0.0 # 2π/λ_y, λ_y = 1e10 k_dot_x = kx*x[1] + ky*x[2] @@ -49,7 +49,7 @@ function initial_condition_jeans_instability(x, t, equations::HyperbolicDiffusionEquations2D) # gravity equation: -Δϕ = -4πGρ # Constants taken from the FLASH manual - # https://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel.pdf + # https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node189.html#SECTION010131000000000000000 rho0 = 1.5e7 delta0 = 1e-3 diff --git a/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_sedov_blast_wave.jl b/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_sedov_blast_wave.jl index f7bb5bbb01..8933224a2c 100644 --- a/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_sedov_blast_wave.jl +++ b/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_sedov_blast_wave.jl @@ -15,14 +15,14 @@ Adaptation of the Sedov blast wave with self-gravity taken from A purely hyperbolic discontinuous Galerkin approach for self-gravitating gas dynamics [arXiv: 2008.10593](https://arxiv.org/abs/2008.10593) based on -- http://flash.uchicago.edu/site/flashcode/user_support/flash4_ug_4p62/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114100000000000000 Should be used together with [`boundary_condition_sedov_self_gravity`](@ref). """ function initial_condition_sedov_self_gravity(x, t, equations::CompressibleEulerEquations2D) # Set up polar coordinates r = sqrt(x[1]^2 + x[2]^2) - # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash4_ug_4p62/node184.html#SECTION010114000000000000000 + # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114100000000000000 r0 = 0.125 # = 4.0 * smallest dx (for domain length=8 and max-ref=8) E = 1.0 p_inner = (equations.gamma - 1) * E / (pi * r0^2) @@ -59,7 +59,7 @@ Adaptation of the Sedov blast wave with self-gravity taken from A purely hyperbolic discontinuous Galerkin approach for self-gravitating gas dynamics [arXiv: 2008.10593](https://arxiv.org/abs/2008.10593) based on -- http://flash.uchicago.edu/site/flashcode/user_support/flash4_ug_4p62/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114100000000000000 Should be used together with [`initial_condition_sedov_self_gravity`](@ref). """ function boundary_condition_sedov_self_gravity(u_inner, orientation, direction, x, t, @@ -122,7 +122,7 @@ Adaptation of the Sedov blast wave with self-gravity taken from A purely hyperbolic discontinuous Galerkin approach for self-gravitating gas dynamics [arXiv: 2008.10593](https://arxiv.org/abs/2008.10593) based on -- http://flash.uchicago.edu/site/flashcode/user_support/flash4_ug_4p62/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114100000000000000 Should be used together with [`boundary_condition_sedov_self_gravity`](@ref). """ function initial_condition_sedov_self_gravity(x, t, equations::HyperbolicDiffusionEquations2D) @@ -143,7 +143,7 @@ Adaptation of the Sedov blast wave with self-gravity taken from A purely hyperbolic discontinuous Galerkin approach for self-gravitating gas dynamics [arXiv: 2008.10593](https://arxiv.org/abs/2008.10593) based on -- http://flash.uchicago.edu/site/flashcode/user_support/flash4_ug_4p62/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114100000000000000 Should be used together with [`initial_condition_sedov_self_gravity`](@ref). """ function boundary_condition_sedov_self_gravity(u_inner, orientation, direction, x, t, diff --git a/examples/structured_1d_dgsem/elixir_advection_shockcapturing.jl b/examples/structured_1d_dgsem/elixir_advection_shockcapturing.jl index 9a81acfe51..313812fe08 100644 --- a/examples/structured_1d_dgsem/elixir_advection_shockcapturing.jl +++ b/examples/structured_1d_dgsem/elixir_advection_shockcapturing.jl @@ -9,7 +9,9 @@ advection_velocity = 1.0 """ initial_condition_composite(x, t, equations::LinearScalarAdvectionEquation1D) -Slight simplification of +Wave form that is a combination of a Gaussian pulse, a square wave, a triangle wave, +and half an ellipse with periodic boundary conditions. +Slight simplification from - Jiang, Shu (1996) Efficient Implementation of Weighted ENO Schemes [DOI: 10.1006/jcph.1996.0130](https://doi.org/10.1006/jcph.1996.0130) @@ -60,7 +62,7 @@ volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; solver = DGSEM(basis, surface_flux, volume_integral) # Create curved mesh -cells_per_dimension = (125,) +cells_per_dimension = (120,) coordinates_min = (-1.0,) # minimum coordinate coordinates_max = (1.0,) # maximum coordinate mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max, diff --git a/examples/structured_1d_dgsem/elixir_euler_sedov.jl b/examples/structured_1d_dgsem/elixir_euler_sedov.jl index ee950b3aaa..9d7be21a5c 100644 --- a/examples/structured_1d_dgsem/elixir_euler_sedov.jl +++ b/examples/structured_1d_dgsem/elixir_euler_sedov.jl @@ -11,7 +11,7 @@ equations = CompressibleEulerEquations1D(1.4) initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D) The Sedov blast wave setup based on Flash -- http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 """ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D) # Set up polar coordinates @@ -19,7 +19,7 @@ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEq x_norm = x[1] - inicenter[1] r = abs(x_norm) - # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 + # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) # r0 = 0.5 # = more reasonable setup E = 1.0 @@ -78,8 +78,8 @@ save_solution = SaveSolutionCallback(interval=100, stepsize_callback = StepsizeCallback(cfl=0.5) -callbacks = CallbackSet(summary_callback, - analysis_callback, +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, save_solution, stepsize_callback) diff --git a/examples/structured_2d_dgsem/elixir_euler_sedov.jl b/examples/structured_2d_dgsem/elixir_euler_sedov.jl index ed1bfab3be..efc3b6627c 100644 --- a/examples/structured_2d_dgsem/elixir_euler_sedov.jl +++ b/examples/structured_2d_dgsem/elixir_euler_sedov.jl @@ -11,7 +11,7 @@ equations = CompressibleEulerEquations2D(1.4) initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) The Sedov blast wave setup based on Flash -- http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 """ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) # Set up polar coordinates @@ -20,7 +20,7 @@ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEq y_norm = x[2] - inicenter[2] r = sqrt(x_norm^2 + y_norm^2) - # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 + # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) E = 1.0 p0_inner = 3 * (equations.gamma - 1) * E / (3 * pi * r0^2) @@ -59,12 +59,12 @@ function mapping(xi, eta) y = eta + 0.125 * (cos(1.5 * pi * xi) * cos(0.5 * pi * eta)) x = xi + 0.125 * (cos(0.5 * pi * xi) * cos(2 * pi * y)) - + return SVector(x, y) end - + cells_per_dimension = (16, 16) - + mesh = StructuredMesh(cells_per_dimension, mapping, periodicity=true) # create the semidiscretization diff --git a/examples/structured_3d_dgsem/elixir_euler_sedov.jl b/examples/structured_3d_dgsem/elixir_euler_sedov.jl index 8f428495b4..e0595437c9 100644 --- a/examples/structured_3d_dgsem/elixir_euler_sedov.jl +++ b/examples/structured_3d_dgsem/elixir_euler_sedov.jl @@ -11,7 +11,7 @@ equations = CompressibleEulerEquations3D(1.4) initial_condition_medium_sedov_blast_wave(x, t, equations::CompressibleEulerEquations3D) The Sedov blast wave setup based on Flash -- http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 with smaller strength of the initial discontinuity. """ function initial_condition_medium_sedov_blast_wave(x, t, equations::CompressibleEulerEquations3D) @@ -22,11 +22,11 @@ function initial_condition_medium_sedov_blast_wave(x, t, equations::Compressible z_norm = x[3] - inicenter[3] r = sqrt(x_norm^2 + y_norm^2 + z_norm^2) - # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 + # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) E = 1.0 p0_inner = 3 * (equations.gamma - 1) * E / (4 * pi * r0^2) - p0_outer = 1.0e-3 + p0_outer = 1.0e-3 # Calculate primitive variables rho = 1.0 @@ -52,8 +52,8 @@ indicator_sc = IndicatorHennemannGassner(equations, basis, volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; volume_flux_dg=volume_flux, volume_flux_fv=surface_flux) - -solver = DGSEM(polydeg=polydeg, surface_flux=surface_flux, volume_integral=volume_integral) + +solver = DGSEM(polydeg=polydeg, surface_flux=surface_flux, volume_integral=volume_integral) # Mapping as described in https://arxiv.org/abs/2012.12040 function mapping(xi, eta, zeta) diff --git a/examples/tree_1d_dgsem/elixir_burgers_shock.jl b/examples/tree_1d_dgsem/elixir_burgers_shock.jl index 987fb320ad..00b5314e19 100644 --- a/examples/tree_1d_dgsem/elixir_burgers_shock.jl +++ b/examples/tree_1d_dgsem/elixir_burgers_shock.jl @@ -21,7 +21,7 @@ surface_flux = flux_lax_friedrichs volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; volume_flux_dg=surface_flux, volume_flux_fv=surface_flux) - + solver = DGSEM(basis, surface_flux, volume_integral) coordinate_min = 0.0 @@ -59,7 +59,7 @@ end boundary_conditions = (x_neg=boundary_condition_inflow, x_pos=boundary_condition_outflow) - + initial_condition = initial_condition_shock semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, @@ -79,7 +79,7 @@ analysis_callback = AnalysisCallback(semi, interval=analysis_interval) alive_callback = AliveCallback(analysis_interval=analysis_interval) -stepsize_callback = StepsizeCallback(cfl=0.9) +stepsize_callback = StepsizeCallback(cfl=0.85) callbacks = CallbackSet(summary_callback, diff --git a/examples/tree_1d_dgsem/elixir_euler_positivity.jl b/examples/tree_1d_dgsem/elixir_euler_positivity.jl index 7942937151..966661e889 100644 --- a/examples/tree_1d_dgsem/elixir_euler_positivity.jl +++ b/examples/tree_1d_dgsem/elixir_euler_positivity.jl @@ -11,7 +11,7 @@ equations = CompressibleEulerEquations1D(1.4) initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D) The Sedov blast wave setup based on Flash -- http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 """ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D) # Set up polar coordinates @@ -19,7 +19,7 @@ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEq x_norm = x[1] - inicenter[1] r = abs(x_norm) - # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 + # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) # r0 = 0.5 # = more reasonable setup E = 1.0 diff --git a/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave.jl b/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave.jl index 746a7cf1ba..106ccacf4f 100644 --- a/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave.jl +++ b/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave.jl @@ -11,7 +11,7 @@ equations = CompressibleEulerEquations1D(1.4) initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D) The Sedov blast wave setup based on Flash -- http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 """ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D) # Set up polar coordinates @@ -19,7 +19,7 @@ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEq x_norm = x[1] - inicenter[1] r = abs(x_norm) - # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 + # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) # r0 = 0.5 # = more reasonable setup E = 1.0 diff --git a/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave_pure_fv.jl b/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave_pure_fv.jl index 00b80dbae9..ebe8fa7ceb 100644 --- a/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave_pure_fv.jl +++ b/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave_pure_fv.jl @@ -11,7 +11,7 @@ equations = CompressibleEulerEquations1D(1.4) initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D) The Sedov blast wave setup based on Flash -- http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 """ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D) # Set up polar coordinates @@ -19,7 +19,7 @@ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEq x_norm = x[1] - inicenter[1] r = abs(x_norm) - # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 + # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) # r0 = 0.5 # = more reasonable setup E = 1.0 diff --git a/examples/tree_1d_dgsem/elixir_eulermulti_two_interacting_blast_waves.jl b/examples/tree_1d_dgsem/elixir_eulermulti_two_interacting_blast_waves.jl index 353093e5f7..8196619418 100644 --- a/examples/tree_1d_dgsem/elixir_eulermulti_two_interacting_blast_waves.jl +++ b/examples/tree_1d_dgsem/elixir_eulermulti_two_interacting_blast_waves.jl @@ -88,7 +88,7 @@ ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() -analysis_interval = 100 +analysis_interval = 1000 analysis_callback = AnalysisCallback(semi, interval=analysis_interval) diff --git a/examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube.jl b/examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube.jl index 1c07fc4fdd..c5727109d9 100644 --- a/examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube.jl +++ b/examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube.jl @@ -94,7 +94,7 @@ amr_callback = AMRCallback(semi, amr_controller, adapt_initial_condition=true, adapt_initial_condition_only_refine=true) -stepsize_callback = StepsizeCallback(cfl=0.8) +stepsize_callback = StepsizeCallback(cfl=0.65) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, diff --git a/examples/tree_1d_dgsem/elixir_mhdmulti_ec.jl b/examples/tree_1d_dgsem/elixir_mhdmulti_ec.jl index 34fdce6634..69ea0551be 100644 --- a/examples/tree_1d_dgsem/elixir_mhdmulti_ec.jl +++ b/examples/tree_1d_dgsem/elixir_mhdmulti_ec.jl @@ -4,8 +4,8 @@ using Trixi ############################################################################### # semidiscretization of the ideal MHD equations -equations = IdealGlmMhdMulticomponentEquations1D(gammas = (2.0, 2.0, 2.0), - gas_constants = (2.0, 2.0, 2.0)) +equations = IdealGlmMhdMulticomponentEquations1D(gammas = (2.0, 2.0, 2.0), + gas_constants = (2.0, 2.0, 2.0)) initial_condition = initial_condition_weak_blast_wave diff --git a/examples/tree_1d_dgsem/elixir_mhdmulti_es.jl b/examples/tree_1d_dgsem/elixir_mhdmulti_es.jl index 8ca32194b9..93cf3e0fdb 100644 --- a/examples/tree_1d_dgsem/elixir_mhdmulti_es.jl +++ b/examples/tree_1d_dgsem/elixir_mhdmulti_es.jl @@ -4,8 +4,8 @@ using Trixi ############################################################################### # semidiscretization of the ideal MHD equations -equations = IdealGlmMhdMulticomponentEquations1D(gammas = (2.0, 2.0, 2.0), - gas_constants = (2.0, 2.0, 2.0)) +equations = IdealGlmMhdMulticomponentEquations1D(gammas = (2.0, 2.0, 2.0), + gas_constants = (2.0, 2.0, 2.0)) initial_condition = initial_condition_weak_blast_wave diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl b/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl index be6a2cb166..1469afec1c 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl @@ -8,9 +8,34 @@ using Trixi equations = ShallowWaterEquations1D(gravity_constant=9.81) -# Note, this initial condition is used to compute errors in the analysis callback but the initialization is -# overwritten by `initial_condition_ec_discontinuous_bottom` below. -initial_condition = initial_condition_weak_blast_wave +# Initial condition with a truly discontinuous water height, velocity, and bottom +# topography function as an academic testcase for entropy conservation. +# The errors from the analysis callback are not important but `∑∂S/∂U ⋅ Uₜ` should +# be around machine roundoff. +# Works as intended for TreeMesh1D with `initial_refinement_level=4`. If the mesh +# refinement level is changed the initial condition below may need changed as well to +# ensure that the discontinuities lie on an element interface. +function initial_condition_ec_discontinuous_bottom(x, t, equations::ShallowWaterEquations1D) + # Set the background values + H = 4.25 + v = 0.0 + b = sin(x[1]) # arbitrary continuous function + + # Setup the discontinuous water height and velocity + if x[1] >= 0.125 && x[1] <= 0.25 + H = 5.0 + v = 0.1882 + end + + # Setup a discontinuous bottom topography + if x[1] >= -0.25 && x[1] <= -0.125 + b = 2.0 + 0.5 * sin(2.0 * pi * x[1]) + end + + return prim2cons(SVector(H, v, b), equations) +end + +initial_condition = initial_condition_ec_discontinuous_bottom ############################################################################### # Get the DG approximation space @@ -37,46 +62,6 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) tspan = (0.0, 2.0) ode = semidiscretize(semi, tspan) -############################################################################### -# Workaround to set a discontinuous bottom topography and initial condition for debugging and testing. - -# alternative version of the initial conditinon used to setup a truly discontinuous -# bottom topography function and initial condition for this academic testcase of entropy conservation. -# The errors from the analysis callback are not important but `∑∂S/∂U ⋅ Uₜ` should be around machine roundoff -# In contrast to the usual signature of initial conditions, this one get passed the -# `element_id` explicitly. In particular, this initial conditions works as intended -# only for the TreeMesh1D with `initial_refinement_level=4`. -function initial_condition_ec_discontinuous_bottom(x, t, element_id, equations::ShallowWaterEquations1D) - # Set the background values - H = 4.25 - v = 0.0 - b = sin(x[1]) # arbitrary continuous function - - # setup the discontinuous water height and velocity - if element_id == 10 - H = 5.0 - v = 0.1882 - end - - # Setup a discontinuous bottom topography using the element id number - if element_id == 7 - b = 2.0 + 0.5 * sin(2.0 * pi * x[1]) - end - - return prim2cons(SVector(H, v, b), equations) -end - -# point to the data we want to augment -u = Trixi.wrap_array(ode.u0, semi) -# reset the initial condition -for element in eachelement(semi.solver, semi.cache) - for i in eachnode(semi.solver) - x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations, semi.solver, i, element) - u_node = initial_condition_ec_discontinuous_bottom(x_node, first(tspan), element, equations) - Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, element) - end -end - ############################################################################### # Callbacks diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_shock_capturing.jl b/examples/tree_1d_dgsem/elixir_shallowwater_shock_capturing.jl index 50241126a2..62346d7b5a 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_shock_capturing.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_shock_capturing.jl @@ -7,24 +7,37 @@ using Trixi equations = ShallowWaterEquations1D(gravity_constant=9.812, H0=1.75) -function initial_condition_stone_throw(x, t, equations::ShallowWaterEquations1D) - # Set up polar coordinates - inicenter = 0.15 - x_norm = x[1] - inicenter[1] - r = abs(x_norm) +# Initial condition with a truly discontinuous velocity and bottom topography. +# Works as intended for TreeMesh1D with `initial_refinement_level=3`. If the mesh +# refinement level is changed the initial condition below may need changed as well to +# ensure that the discontinuities lie on an element interface. +function initial_condition_stone_throw_discontinuous_bottom(x, t, equations::ShallowWaterEquations1D) + + # Calculate primitive variables + + # flat lake + H = equations.H0 + + # Discontinuous velocity + v = 0.0 + if x[1] >= -0.75 && x[1] <= 0.0 + v = -1.0 + elseif x[1] >= 0.0 && x[1] <= 0.75 + v = 1.0 + end - # Calculate primitive variables - H = equations.H0 - # v = 0.0 # for well-balanced test - v = r < 0.6 ? 1.75 : 0.0 # for stone throw + b = ( 1.5 / exp( 0.5 * ((x[1] - 1.0)^2 ) ) + + 0.75 / exp( 0.5 * ((x[1] + 1.0)^2 ) ) ) - b = ( 1.5 / exp( 0.5 * ((x[1] - 1.0)^2 ) ) - + 0.75 / exp( 0.5 * ((x[1] + 1.0)^2 ) ) ) + # Force a discontinuous bottom topography + if x[1] >= -1.5 && x[1] <= 0.0 + b = 0.5 + end - return prim2cons(SVector(H, v, b), equations) + return prim2cons(SVector(H, v, b), equations) end -initial_condition = initial_condition_stone_throw +initial_condition = initial_condition_stone_throw_discontinuous_bottom boundary_condition = boundary_condition_slip_wall @@ -62,49 +75,13 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions = boundary_condition) ############################################################################### -# ODE solvers, callbacks etc. +# ODE solver tspan = (0.0, 3.0) ode = semidiscretize(semi, tspan) -# Hack in a discontinuous bottom for a more interesting test -function initial_condition_stone_throw_discontinuous_bottom(x, t, element_id, equations::ShallowWaterEquations1D) - - inicenter = 0.15 - x_norm = x[1] - inicenter[1] - r = abs(x_norm) - - # Calculate primitive variables - H = equations.H0 # flat lake - # Discontinuous velocity set via element id number - v = 0.0 - if element_id == 4 - v = -1.0 - elseif element_id == 5 - v = 1.0 - end - - b = ( 1.5 / exp( 0.5 * ((x[1] - 1.0)^2 ) ) - + 0.75 / exp( 0.5 * ((x[1] + 1.0)^2 ) ) ) - - # Setup a discontinuous bottom topography using the element id number - if element_id == 3 || element_id == 4 - b = 0.5 - end - - return prim2cons(SVector(H, v, b), equations) -end - -# point to the data we want to augment -u = Trixi.wrap_array(ode.u0, semi) -# reset the initial condition -for element in eachelement(semi.solver, semi.cache) - for i in eachnode(semi.solver) - x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations, semi.solver, i, element) - u_node = initial_condition_stone_throw_discontinuous_bottom(x_node, first(tspan), element, equations) - Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, element) - end -end +############################################################################### +# Callbacks summary_callback = SummaryCallback() diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl index 67c1098b17..60770d158f 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl @@ -3,20 +3,34 @@ using OrdinaryDiffEq using Trixi ############################################################################### -# Semidiscretization of the two-layer shallow water equations for a dam break test with a -# discontinuous bottom topography function to test entropy conservation +# Semidiscretization of the two-layer shallow water equations for a dam break +# test with a discontinuous bottom topography function to test entropy conservation equations = ShallowWaterTwoLayerEquations1D(gravity_constant=9.81, H0=2.0, rho_upper=0.9, rho_lower=1.0) -############################################################################### -# Workaround to set a discontinuous bottom topography and initial condition. +# Initial condition of a dam break with a discontinuous water heights and bottom topography. +# Works as intended for TreeMesh1D with `initial_refinement_level=5`. If the mesh +# refinement level is changed the initial condition below may need changed as well to +# ensure that the discontinuities lie on an element interface. +function initial_condition_dam_break(x, t, equations::ShallowWaterTwoLayerEquations1D) + v1_upper = 0.0 + v1_lower = 0.0 + + # Set the discontinuity + if x[1] <= 10.0 + H_lower = 2.0 + H_upper = 4.0 + b = 0.0 + else + H_lower = 1.5 + H_upper = 3.0 + b = 0.5 + end + + return prim2cons(SVector(H_upper, v1_upper, H_lower, v1_lower, b), equations) +end -# This test case uses a special work around to setup a truly discontinuous bottom topography -# function and initial condition for this academic testcase of entropy conservation. First, a -# dummy initial condition is introduced to create the semidiscretization. Then the initial condition -# is reset with the true discontinuous values from initial_condition_dam_break. Note, that this -# initial condition only works for TreeMesh1D with `initial_refinement_level=5`. -initial_condition = initial_condition_convergence_test +initial_condition = initial_condition_dam_break ############################################################################### # Get the DG approximation space @@ -25,7 +39,6 @@ volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) solver = DGSEM(polydeg=3, surface_flux=(flux_fjordholm_etal, flux_nonconservative_fjordholm_etal), volume_integral=VolumeIntegralFluxDifferencing(volume_flux)) - ############################################################################### # Get the TreeMesh and setup a non-periodic mesh @@ -39,54 +52,22 @@ mesh = TreeMesh(coordinates_min, coordinates_max, boundary_condition = boundary_condition_slip_wall # create the semidiscretization object -semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions=boundary_condition) ############################################################################### -# ODE solvers, callbacks etc. +# ODE solvers tspan = (0.0,0.4) ode = semidiscretize(semi, tspan) -# Initial conditions dam break test case -function initial_condition_dam_break(x, t, element_id, equations::ShallowWaterTwoLayerEquations1D) - v1_upper = 0.0 - v1_lower = 0.0 - - # Set the discontinuity - if element_id <= 16 - H_lower = 2.0 - H_upper = 4.0 - b = 0.0 - else - H_lower = 1.5 - H_upper = 3.0 - b = 0.5 - end - - return prim2cons(SVector(H_upper, v1_upper, H_lower, v1_lower, b), equations) -end - - -# point to the data we want to augment -u = Trixi.wrap_array(ode.u0, semi) -# reset the initial condition -for element in eachelement(semi.solver, semi.cache) - for i in eachnode(semi.solver) - x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, - equations, semi.solver, i, element) - u_node = initial_condition_dam_break(x_node, first(tspan), element, equations) - Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, element) - end -end - - - +############################################################################### +# Callbacks summary_callback = SummaryCallback() analysis_interval = 500 -analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=false, +analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=false, extra_analysis_integrals=(energy_total, energy_kinetic, energy_internal,)) stepsize_callback = StepsizeCallback(cfl=1.0) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl index 8383565683..e07bc04d76 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl @@ -8,21 +8,28 @@ using Trixi equations = ShallowWaterEquations1D(gravity_constant=9.81, H0=3.25) -# An initial condition with constant total water height and zero velocities to test well-balancedness. -# Note, this routine is used to compute errors in the analysis callback but the initialization is -# overwritten by `initial_condition_discontinuous_well_balancedness` below. -function initial_condition_well_balancedness(x, t, equations::ShallowWaterEquations1D) +# Setup a truly discontinuous bottom topography function for this academic +# testcase of well-balancedness. The errors from the analysis callback are +# not important but the error for this lake-at-rest test case +# `∑|H0-(h+b)|` should be around machine roundoff. +# Works as intended for TreeMesh1D with `initial_refinement_level=3`. If the mesh +# refinement level is changed the initial condition below may need changed as well to +# ensure that the discontinuities lie on an element interface. +function initial_condition_discontinuous_well_balancedness(x, t, equations::ShallowWaterEquations1D) # Set the background values H = equations.H0 v = 0.0 + b = 0.0 - # bottom topography inspired by from Pond.control in [HOHQMesh](https://github.com/trixi-framework/HOHQMesh) - b = (1.5 / exp( 0.5 * ((x[1] - 1.0)^2) )+ 0.75 / exp(0.5 * ((x[1] + 1.0)^2))) + # Setup a discontinuous bottom topography + if x[1] >= 0.5 && x[1] <= 0.75 + b = 2.0 + 0.5 * sin(2.0 * pi * x[1]) + end return prim2cons(SVector(H, v, b), equations) end -initial_condition = initial_condition_well_balancedness +initial_condition = initial_condition_discontinuous_well_balancedness ############################################################################### # Get the DG approximation space @@ -50,41 +57,6 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) tspan = (0.0, 100.0) ode = semidiscretize(semi, tspan) -############################################################################### -# Workaround to set a discontinuous bottom topography and initial condition for debugging and testing. - -# alternative version of the initial conditinon used to setup a truly discontinuous -# bottom topography function for this academic testcase of well-balancedness. -# The errors from the analysis callback are not important but the error for this lake at rest test case -# `∑|H0-(h+b)|` should be around machine roundoff. -# In contrast to the usual signature of initial conditions, this one get passed the -# `element_id` explicitly. In particular, this initial conditions works as intended -# only for the TreeMesh1D with `initial_refinement_level=3`. -function initial_condition_discontinuous_well_balancedness(x, t, element_id, equations::ShallowWaterEquations1D) - # Set the background values - H = equations.H0 - v = 0.0 - b = 0.0 - - # Setup a discontinuous bottom topography using the element id number - if element_id == 7 - b = 2.0 + 0.5 * sin(2.0 * pi * x[1]) - end - - return prim2cons(SVector(H, v, b), equations) -end - -# point to the data we want to augment -u = Trixi.wrap_array(ode.u0, semi) -# reset the initial condition -for element in eachelement(semi.solver, semi.cache) - for i in eachnode(semi.solver) - x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations, semi.solver, i, element) - u_node = initial_condition_discontinuous_well_balancedness(x_node, first(tspan), element, equations) - Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, element) - end -end - ############################################################################### # Callbacks diff --git a/examples/tree_1d_fdsbp/elixir_advection_upwind.jl b/examples/tree_1d_fdsbp/elixir_advection_upwind.jl new file mode 100644 index 0000000000..5c50e1a6c6 --- /dev/null +++ b/examples/tree_1d_fdsbp/elixir_advection_upwind.jl @@ -0,0 +1,57 @@ +# !!! warning "Experimental implementation (upwind SBP)" +# This is an experimental feature and may change in future releases. + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear scalar advection equation equation + +equations = LinearScalarAdvectionEquation1D(1.0) + +function initial_condition_sin(x, t, equation::LinearScalarAdvectionEquation1D) + return SVector(sinpi(x[1] - equations.advection_velocity[1] * t)) +end + +D_upw = upwind_operators(SummationByPartsOperators.Mattsson2017, + derivative_order = 1, + accuracy_order = 4, + xmin = -1.0, xmax = 1.0, + N = 16) +flux_splitting = splitting_lax_friedrichs +solver = FDSBP(D_upw, + surface_integral = SurfaceIntegralUpwind(flux_splitting), + volume_integral = VolumeIntegralUpwind(flux_splitting)) + +coordinates_min = -1.0 +coordinates_max = 1.0 +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + n_cells_max = 10_000) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_sin, solver) + + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 2.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval) + +alive_callback = AliveCallback(analysis_interval=analysis_interval) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback) + + +############################################################################### +# run the simulation + +sol = solve(ode, RDPK3SpFSAL49(); abstol=1.0e-6, reltol=1.0e-6, + ode_default_options()..., callback=callbacks); +summary_callback() # print the timer summary diff --git a/examples/tree_2d_dgsem/elixir_euler_positivity.jl b/examples/tree_2d_dgsem/elixir_euler_positivity.jl index e40dc3b47a..4c7dd7eb6c 100644 --- a/examples/tree_2d_dgsem/elixir_euler_positivity.jl +++ b/examples/tree_2d_dgsem/elixir_euler_positivity.jl @@ -11,7 +11,7 @@ equations = CompressibleEulerEquations2D(gamma) initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) The Sedov blast wave setup based on Flash -- http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 """ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) # Set up polar coordinates @@ -20,7 +20,7 @@ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEq y_norm = x[2] - inicenter[2] r = sqrt(x_norm^2 + y_norm^2) - # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 + # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) # r0 = 0.5 # = more reasonable setup E = 1.0 diff --git a/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave.jl b/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave.jl index da7e1d55c9..512e582237 100644 --- a/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave.jl +++ b/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave.jl @@ -11,7 +11,7 @@ equations = CompressibleEulerEquations2D(gamma) initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) The Sedov blast wave setup based on Flash -- http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 """ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) # Set up polar coordinates @@ -20,7 +20,7 @@ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEq y_norm = x[2] - inicenter[2] r = sqrt(x_norm^2 + y_norm^2) - # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 + # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) # r0 = 0.5 # = more reasonable setup E = 1.0 diff --git a/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_neuralnetwork_perssonperaire.jl b/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_neuralnetwork_perssonperaire.jl index 5671578937..5fd32da2e5 100644 --- a/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_neuralnetwork_perssonperaire.jl +++ b/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_neuralnetwork_perssonperaire.jl @@ -23,7 +23,7 @@ equations = CompressibleEulerEquations2D(gamma) initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) The Sedov blast wave setup based on Flash -- http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 """ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) # Set up polar coordinates @@ -32,7 +32,7 @@ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEq y_norm = x[2] - inicenter[2] r = sqrt(x_norm^2 + y_norm^2) - # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 + # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) # r0 = 0.5 # = more reasonable setup E = 1.0 diff --git a/examples/tree_2d_fdsbp/elixir_euler_convergence.jl b/examples/tree_2d_fdsbp/elixir_euler_convergence.jl index 0843cece67..2a6c291f0b 100644 --- a/examples/tree_2d_fdsbp/elixir_euler_convergence.jl +++ b/examples/tree_2d_fdsbp/elixir_euler_convergence.jl @@ -5,7 +5,8 @@ using OrdinaryDiffEq using Trixi ############################################################################### -# semidiscretization of the linear advection equation +# semidiscretization of the compressible Euler equations + equations = CompressibleEulerEquations2D(1.4) initial_condition = initial_condition_convergence_test diff --git a/examples/tree_2d_fdsbp/elixir_euler_kelvin_helmholtz_instability.jl b/examples/tree_2d_fdsbp/elixir_euler_kelvin_helmholtz_instability.jl index 1e58badf47..e63343852a 100644 --- a/examples/tree_2d_fdsbp/elixir_euler_kelvin_helmholtz_instability.jl +++ b/examples/tree_2d_fdsbp/elixir_euler_kelvin_helmholtz_instability.jl @@ -5,7 +5,8 @@ using OrdinaryDiffEq using Trixi ############################################################################### -# semidiscretization of the linear advection equation +# semidiscretization of the compressible Euler equations + equations = CompressibleEulerEquations2D(1.4) function initial_condition_kelvin_helmholtz_instability(x, t, equations::CompressibleEulerEquations2D) diff --git a/examples/tree_2d_fdsbp/elixir_euler_vortex.jl b/examples/tree_2d_fdsbp/elixir_euler_vortex.jl index abaf3d494d..c1bee8f9c4 100644 --- a/examples/tree_2d_fdsbp/elixir_euler_vortex.jl +++ b/examples/tree_2d_fdsbp/elixir_euler_vortex.jl @@ -5,7 +5,8 @@ using OrdinaryDiffEq using Trixi ############################################################################### -# semidiscretization of the linear advection equation +# semidiscretization of the compressible Euler equations + equations = CompressibleEulerEquations2D(1.4) """ diff --git a/examples/tree_3d_dgsem/elixir_euler_sedov_blast_wave.jl b/examples/tree_3d_dgsem/elixir_euler_sedov_blast_wave.jl index 336c09e921..3641878149 100644 --- a/examples/tree_3d_dgsem/elixir_euler_sedov_blast_wave.jl +++ b/examples/tree_3d_dgsem/elixir_euler_sedov_blast_wave.jl @@ -15,14 +15,14 @@ Adaptation of the Sedov blast wave with self-gravity taken from A purely hyperbolic discontinuous Galerkin approach for self-gravitating gas dynamics [arXiv: 2008.10593](https://arxiv.org/abs/2008.10593) based on -- http://flash.uchicago.edu/site/flashcode/user_support/flash4_ug_4p62/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 Should be used together with [`boundary_condition_sedov_self_gravity`](@ref). """ function initial_condition_sedov_self_gravity(x, t, equations::CompressibleEulerEquations3D) # Calculate radius as distance from origin r = sqrt(x[1]^2 + x[2]^2 + x[3]^2) - # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash4_ug_4p62/node184.html#SECTION010114000000000000000 + # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 r0 = 0.25 # = 4.0 * smallest dx (for domain length=8 and max-ref=7) E = 1.0 p_inner = (equations.gamma - 1) * E / (4/3 * pi * r0^3) @@ -60,7 +60,7 @@ Adaptation of the Sedov blast wave with self-gravity taken from A purely hyperbolic discontinuous Galerkin approach for self-gravitating gas dynamics [arXiv: 2008.10593](https://arxiv.org/abs/2008.10593) based on -- http://flash.uchicago.edu/site/flashcode/user_support/flash4_ug_4p62/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 Should be used together with [`initial_condition_sedov_self_gravity`](@ref). """ function boundary_condition_sedov_self_gravity(u_inner, orientation, direction, x, t, diff --git a/examples/tree_3d_fdsbp/elixir_euler_convergence.jl b/examples/tree_3d_fdsbp/elixir_euler_convergence.jl index 576a07e6ab..6aafa1b5cc 100644 --- a/examples/tree_3d_fdsbp/elixir_euler_convergence.jl +++ b/examples/tree_3d_fdsbp/elixir_euler_convergence.jl @@ -6,7 +6,6 @@ using Trixi ############################################################################### # semidiscretization of the compressible Euler equations - equations = CompressibleEulerEquations3D(1.4) initial_condition = initial_condition_convergence_test diff --git a/examples/unstructured_2d_dgsem/elixir_euler_sedov.jl b/examples/unstructured_2d_dgsem/elixir_euler_sedov.jl index 3d5a391bd9..570a208469 100644 --- a/examples/unstructured_2d_dgsem/elixir_euler_sedov.jl +++ b/examples/unstructured_2d_dgsem/elixir_euler_sedov.jl @@ -11,7 +11,7 @@ equations = CompressibleEulerEquations2D(1.4) initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) The Sedov blast wave setup based on Flash -- http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 """ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) # Set up polar coordinates @@ -20,7 +20,7 @@ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEq y_norm = x[2] - inicenter[2] r = sqrt(x_norm^2 + y_norm^2) - # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 + # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) E = 1.0 p0_inner = 3 * (equations.gamma - 1) * E / (3 * pi * r0^2) diff --git a/src/Trixi.jl b/src/Trixi.jl index 4fcb859a3c..6fc62f5052 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -83,7 +83,8 @@ import SummationByPartsOperators: integrate, semidiscretize, upwind_operators # DGMulti solvers -@reexport using StartUpDG: StartUpDG, Polynomial, Gauss, SBP, Line, Tri, Quad, Hex, Tet +@reexport using StartUpDG: StartUpDG, Polynomial, Gauss, TensorProductWedge, SBP, Line, Tri, + Quad, Hex, Tet, Wedge using StartUpDG: RefElemData, MeshData, AbstractElemShape # TODO: include_optimized diff --git a/src/auxiliary/auxiliary.jl b/src/auxiliary/auxiliary.jl index 115d055c0c..1f7d30d6aa 100644 --- a/src/auxiliary/auxiliary.jl +++ b/src/auxiliary/auxiliary.jl @@ -132,7 +132,7 @@ end default_example() Return the path to an example elixir that can be used to quickly see Trixi.jl in action on a -[`TreeMesh`]@(ref). See also [`examples_dir`](@ref) and [`get_examples`](@ref). +[`TreeMesh`](@ref). See also [`examples_dir`](@ref) and [`get_examples`](@ref). """ function default_example() joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_basic.jl") @@ -142,7 +142,7 @@ end default_example_unstructured() Return the path to an example elixir that can be used to quickly see Trixi.jl in action on an -[`UnstructuredMesh2D`]@(ref). This simulation is run on the example curved, unstructured mesh +[`UnstructuredMesh2D`](@ref). This simulation is run on the example curved, unstructured mesh given in the Trixi.jl documentation regarding unstructured meshes. """ function default_example_unstructured() @@ -155,7 +155,7 @@ end Return the default options for OrdinaryDiffEq's `solve`. Pass `ode_default_options()...` to `solve` to only return the solution at the final time and enable **MPI aware** error-based step size control, whenever MPI is used. -For example, use `solve(ode, alg; ode_default_options()...)` +For example, use `solve(ode, alg; ode_default_options()...)`. """ function ode_default_options() if mpi_isparallel() @@ -213,8 +213,8 @@ might be provided by other packages such as [Polyester.jl](https://github.com/Ju This macro does not necessarily work for general `for` loops. For example, it does not necessarily support general iterables such as `eachline(filename)`. -Some discussion can be found at https://discourse.julialang.org/t/overhead-of-threads-threads/53964 -and https://discourse.julialang.org/t/threads-threads-with-one-thread-how-to-remove-the-overhead/58435. +Some discussion can be found at [https://discourse.julialang.org/t/overhead-of-threads-threads/53964](https://discourse.julialang.org/t/overhead-of-threads-threads/53964) +and [https://discourse.julialang.org/t/threads-threads-with-one-thread-how-to-remove-the-overhead/58435](https://discourse.julialang.org/t/threads-threads-with-one-thread-how-to-remove-the-overhead/58435). """ macro threaded(expr) # Use `esc(quote ... end)` for nested macro calls as suggested in diff --git a/src/auxiliary/math.jl b/src/auxiliary/math.jl index 27c1bed5ca..4ecf7dd3fc 100644 --- a/src/auxiliary/math.jl +++ b/src/auxiliary/math.jl @@ -51,7 +51,7 @@ Given ε = 1.0e-4, we use the following algorithm. - Agner Fog. Lists of instruction latencies, throughputs and micro-operation breakdowns for Intel, AMD, and VIA CPUs. - https://www.agner.org/optimize/instruction_tables.pdf + [https://www.agner.org/optimize/instruction_tables.pdf](https://www.agner.org/optimize/instruction_tables.pdf) """ @inline function ln_mean(x, y) epsilon_f2 = 1.0e-4 @@ -166,8 +166,10 @@ checks necessary in the presence of `NaN`s (or signed zeros). # Examples +```jldoctest julia> max(2, 5, 1) 5 +``` """ @inline max(args...) = @fastmath max(args...) @@ -183,8 +185,10 @@ checks necessary in the presence of `NaN`s (or signed zeros). # Examples +```jldoctest julia> min(2, 5, 1) 1 +``` """ @inline min(args...) = @fastmath min(args...) diff --git a/src/auxiliary/mpi.jl b/src/auxiliary/mpi.jl index 2c485b4832..c85c23670b 100644 --- a/src/auxiliary/mpi.jl +++ b/src/auxiliary/mpi.jl @@ -72,7 +72,7 @@ You must pass this function as a keyword argument to OrdinaryDiffEq.jl's `solve` when using error-based step size control with MPI parallel execution of Trixi.jl. -See the "Advanced Adaptive Stepsize Control" section of the [documentation](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/) +See the "Advanced Adaptive Stepsize Control" section of the [documentation](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/). """ ode_norm(u::Number, t) = @fastmath abs(u) function ode_norm(u::AbstractArray, t) @@ -125,6 +125,6 @@ You should pass this function as a keyword argument to OrdinaryDiffEq.jl's `solve` when using error-based step size control with MPI parallel execution of Trixi.jl. -See the "Miscellaneous" section of the [documentation](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/) +See the "Miscellaneous" section of the [documentation](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/). """ ode_unstable_check(dt, u, semi, t) = isnan(dt) diff --git a/src/auxiliary/p4est.jl b/src/auxiliary/p4est.jl index 93b5166cd8..968af339cb 100644 --- a/src/auxiliary/p4est.jl +++ b/src/auxiliary/p4est.jl @@ -24,35 +24,38 @@ function init_p4est() return nothing end +# for convenience to either pass a Ptr or a PointerWrapper +const PointerOrWrapper = Union{Ptr{T}, PointerWrapper{T}} where {T} + # Convert sc_array of type T to Julia array -function unsafe_wrap_sc(::Type{T}, sc_array::Ptr{sc_array}) where {T} - sc_array_obj = unsafe_load(sc_array) +function unsafe_wrap_sc(::Type{T}, sc_array_ptr::Ptr{sc_array}) where {T} + sc_array_obj = unsafe_load(sc_array_ptr) return unsafe_wrap_sc(T, sc_array_obj) end function unsafe_wrap_sc(::Type{T}, sc_array_obj::sc_array) where {T} elem_count = sc_array_obj.elem_count array = sc_array_obj.array - return unsafe_wrap(Array, Ptr{T}(array), elem_count) end -# Load the ith element (1-indexed) of an sc array of type T -function unsafe_load_sc(::Type{T}, sc_array::Ptr{sc_array}, i = 1) where {T} - sc_array_obj = unsafe_load(sc_array) - return unsafe_load_sc(T, sc_array_obj, i) -end +function unsafe_wrap_sc(::Type{T}, sc_array_pw::PointerWrapper{sc_array}) where {T} + elem_count = sc_array_pw.elem_count[] + array = sc_array_pw.array -function unsafe_load_sc(::Type{T}, sc_array_obj::sc_array, i = 1) where {T} - element_size = sc_array_obj.elem_size - @assert element_size == sizeof(T) + return unsafe_wrap(Array, Ptr{T}(pointer(array)), elem_count) +end - return unsafe_load(Ptr{T}(sc_array_obj.array), i) +# Load the ith element (1-indexed) of an sc array of type T as PointerWrapper +function load_pointerwrapper_sc(::Type{T}, sc_array::PointerWrapper{sc_array}, + i::Integer = 1) where {T} + return PointerWrapper(T, pointer(sc_array.array) + (i - 1) * sizeof(T)) end # Create new `p4est` from a p4est_connectivity # 2D -function new_p4est(connectivity::Ptr{p4est_connectivity_t}, initial_refinement_level) +function new_p4est(connectivity::PointerOrWrapper{p4est_connectivity_t}, + initial_refinement_level) comm = P4est.uses_mpi() ? mpi_comm() : 0 # Use Trixi.jl's MPI communicator if p4est supports MPI p4est_new_ext(comm, connectivity, @@ -65,7 +68,8 @@ function new_p4est(connectivity::Ptr{p4est_connectivity_t}, initial_refinement_l end # 3D -function new_p4est(connectivity::Ptr{p8est_connectivity_t}, initial_refinement_level) +function new_p4est(connectivity::PointerOrWrapper{p8est_connectivity_t}, + initial_refinement_level) comm = P4est.uses_mpi() ? mpi_comm() : 0 # Use Trixi.jl's MPI communicator if p4est supports MPI p8est_new_ext(comm, connectivity, 0, initial_refinement_level, true, 2 * sizeof(Int), C_NULL, C_NULL) @@ -73,13 +77,13 @@ end # Save `p4est` data to file # 2D -function save_p4est!(file, p4est::Ptr{p4est_t}) +function save_p4est!(file, p4est::PointerOrWrapper{p4est_t}) # Don't save user data of the quads p4est_save(file, p4est, false) end # 3D -function save_p4est!(file, p8est::Ptr{p8est_t}) +function save_p4est!(file, p8est::PointerOrWrapper{p8est_t}) # Don't save user data of the quads p8est_save(file, p8est, false) end @@ -107,27 +111,33 @@ read_inp_p4est(meshfile, ::Val{3}) = p8est_connectivity_read_inp(meshfile) # Refine `p4est` if refine_fn_c returns 1 # 2D -function refine_p4est!(p4est::Ptr{p4est_t}, recursive, refine_fn_c, init_fn_c) +function refine_p4est!(p4est::PointerOrWrapper{p4est_t}, recursive, refine_fn_c, + init_fn_c) p4est_refine(p4est, recursive, refine_fn_c, init_fn_c) end # 3D -function refine_p4est!(p8est::Ptr{p8est_t}, recursive, refine_fn_c, init_fn_c) +function refine_p4est!(p8est::PointerOrWrapper{p8est_t}, recursive, refine_fn_c, + init_fn_c) p8est_refine(p8est, recursive, refine_fn_c, init_fn_c) end # Refine `p4est` if coarsen_fn_c returns 1 # 2D -function coarsen_p4est!(p4est::Ptr{p4est_t}, recursive, coarsen_fn_c, init_fn_c) +function coarsen_p4est!(p4est::PointerOrWrapper{p4est_t}, recursive, coarsen_fn_c, + init_fn_c) p4est_coarsen(p4est, recursive, coarsen_fn_c, init_fn_c) end # 3D -function coarsen_p4est!(p8est::Ptr{p8est_t}, recursive, coarsen_fn_c, init_fn_c) +function coarsen_p4est!(p8est::PointerOrWrapper{p8est_t}, recursive, coarsen_fn_c, + init_fn_c) p8est_coarsen(p8est, recursive, coarsen_fn_c, init_fn_c) end # Create new ghost layer from p4est, only connections via faces are relevant # 2D -ghost_new_p4est(p4est::Ptr{p4est_t}) = p4est_ghost_new(p4est, P4est.P4EST_CONNECT_FACE) +function ghost_new_p4est(p4est::PointerOrWrapper{p4est_t}) + p4est_ghost_new(p4est, P4est.P4EST_CONNECT_FACE) +end # 3D # In 3D it is not sufficient to use `P8EST_CONNECT_FACE`. Consider the neighbor elements of a mortar # in 3D. We have to determine which MPI ranks are involved in this mortar. @@ -147,28 +157,37 @@ ghost_new_p4est(p4est::Ptr{p4est_t}) = p4est_ghost_new(p4est, P4est.P4EST_CONNEC # `P8EST_CONNECT_FACE`. But if it is not in the ghost layer, it will not be available in # `iterate_p4est` and thus we cannot determine its MPI rank # (see https://github.com/cburstedde/p4est/blob/439bc9aae849555256ddfe4b03d1f9fe8d18ff0e/src/p8est_iterate.h#L66-L72). -ghost_new_p4est(p8est::Ptr{p8est_t}) = p8est_ghost_new(p8est, P4est.P8EST_CONNECT_FULL) +function ghost_new_p4est(p8est::PointerOrWrapper{p8est_t}) + p8est_ghost_new(p8est, P4est.P8EST_CONNECT_FULL) +end # Check if ghost layer is valid # 2D -function ghost_is_valid_p4est(p4est::Ptr{p4est_t}, ghost_layer::Ptr{p4est_ghost_t}) +function ghost_is_valid_p4est(p4est::PointerOrWrapper{p4est_t}, + ghost_layer::Ptr{p4est_ghost_t}) return p4est_ghost_is_valid(p4est, ghost_layer) end # 3D -function ghost_is_valid_p4est(p4est::Ptr{p8est_t}, ghost_layer::Ptr{p8est_ghost_t}) +function ghost_is_valid_p4est(p4est::PointerOrWrapper{p8est_t}, + ghost_layer::Ptr{p8est_ghost_t}) return p8est_ghost_is_valid(p4est, ghost_layer) end # Destroy ghost layer # 2D -ghost_destroy_p4est(ghost_layer::Ptr{p4est_ghost_t}) = p4est_ghost_destroy(ghost_layer) +function ghost_destroy_p4est(ghost_layer::PointerOrWrapper{p4est_ghost_t}) + p4est_ghost_destroy(ghost_layer) +end # 3D -ghost_destroy_p4est(ghost_layer::Ptr{p8est_ghost_t}) = p8est_ghost_destroy(ghost_layer) +function ghost_destroy_p4est(ghost_layer::PointerOrWrapper{p8est_ghost_t}) + p8est_ghost_destroy(ghost_layer) +end # Let `p4est` iterate over each cell volume and cell face. # Call iter_volume_c for each cell and iter_face_c for each face. # 2D -function iterate_p4est(p4est::Ptr{p4est_t}, user_data; ghost_layer = C_NULL, +function iterate_p4est(p4est::PointerOrWrapper{p4est_t}, user_data; + ghost_layer = C_NULL, iter_volume_c = C_NULL, iter_face_c = C_NULL) if user_data === C_NULL user_data_ptr = user_data @@ -191,7 +210,8 @@ function iterate_p4est(p4est::Ptr{p4est_t}, user_data; ghost_layer = C_NULL, end # 3D -function iterate_p4est(p8est::Ptr{p8est_t}, user_data; ghost_layer = C_NULL, +function iterate_p4est(p8est::PointerOrWrapper{p8est_t}, user_data; + ghost_layer = C_NULL, iter_volume_c = C_NULL, iter_face_c = C_NULL) if user_data === C_NULL user_data_ptr = user_data @@ -216,23 +236,25 @@ end # Load i-th element of the sc_array info.sides of the type p[48]est_iter_face_side_t # 2D version -function unsafe_load_side(info::Ptr{p4est_iter_face_info_t}, i = 1) - return unsafe_load_sc(p4est_iter_face_side_t, unsafe_load(info).sides, i) +function load_pointerwrapper_side(info::PointerWrapper{p4est_iter_face_info_t}, + i::Integer = 1) + return load_pointerwrapper_sc(p4est_iter_face_side_t, info.sides, i) end # 3D version -function unsafe_load_side(info::Ptr{p8est_iter_face_info_t}, i = 1) - return unsafe_load_sc(p8est_iter_face_side_t, unsafe_load(info).sides, i) +function load_pointerwrapper_side(info::PointerWrapper{p8est_iter_face_info_t}, + i::Integer = 1) + return load_pointerwrapper_sc(p8est_iter_face_side_t, info.sides, i) end # Load i-th element of the sc_array p4est.trees of the type p[48]est_tree_t # 2D version -function unsafe_load_tree(p4est::Ptr{p4est_t}, i = 1) - return unsafe_load_sc(p4est_tree_t, unsafe_load(p4est).trees, i) +function load_pointerwrapper_tree(p4est::PointerWrapper{p4est_t}, i::Integer = 1) + return load_pointerwrapper_sc(p4est_tree_t, p4est.trees, i) end # 3D version -function unsafe_load_tree(p8est::Ptr{p8est_t}, i = 1) - return unsafe_load_sc(p8est_tree_t, unsafe_load(p8est).trees, i) +function load_pointerwrapper_tree(p8est::PointerWrapper{p8est_t}, i::Integer = 1) + return load_pointerwrapper_sc(p8est_tree_t, p8est.trees, i) end end # @muladd diff --git a/src/callbacks_step/amr.jl b/src/callbacks_step/amr.jl index d6e19b7988..bef49b4c48 100644 --- a/src/callbacks_step/amr.jl +++ b/src/callbacks_step/amr.jl @@ -348,24 +348,24 @@ end # Copy controller values to quad user data storage, will be called below function copy_to_quad_iter_volume(info, user_data) - info_obj = unsafe_load(info) + info_pw = PointerWrapper(info) # Load tree from global trees array, one-based indexing - tree = unsafe_load_tree(info_obj.p4est, info_obj.treeid + 1) + tree_pw = load_pointerwrapper_tree(info_pw.p4est, info_pw.treeid[] + 1) # Quadrant numbering offset of this quadrant - offset = tree.quadrants_offset + offset = tree_pw.quadrants_offset[] # Global quad ID - quad_id = offset + info_obj.quadid + quad_id = offset + info_pw.quadid[] # Access user_data = lambda - user_data_ptr = Ptr{Int}(user_data) + user_data_pw = PointerWrapper(Int, user_data) # Load controller_value = lambda[quad_id + 1] - controller_value = unsafe_load(user_data_ptr, quad_id + 1) + controller_value = user_data_pw[quad_id + 1] # Access quadrant's user data ([global quad ID, controller_value]) - quad_data_ptr = Ptr{Int}(unsafe_load(info_obj.quad.p.user_data)) + quad_data_pw = PointerWrapper(Int, info_pw.quad.p.user_data[]) # Save controller value to quadrant's user data. - unsafe_store!(quad_data_ptr, controller_value, 2) + quad_data_pw[2] = controller_value return nothing end @@ -599,22 +599,22 @@ function current_element_levels(mesh::TreeMesh, solver, cache) end function extract_levels_iter_volume(info, user_data) - info_obj = unsafe_load(info) + info_pw = PointerWrapper(info) # Load tree from global trees array, one-based indexing - tree = unsafe_load_tree(info_obj.p4est, info_obj.treeid + 1) + tree_pw = load_pointerwrapper_tree(info_pw.p4est, info_pw.treeid[] + 1) # Quadrant numbering offset of this quadrant - offset = tree.quadrants_offset + offset = tree_pw.quadrants_offset[] # Global quad ID - quad_id = offset + info_obj.quadid + quad_id = offset + info_pw.quadid[] # Julia element ID element_id = quad_id + 1 - current_level = unsafe_load(info_obj.quad.level) + current_level = info_pw.quad.level[] # Unpack user_data = current_levels and save current element level - ptr = Ptr{Int}(user_data) - unsafe_store!(ptr, current_level, element_id) + pw = PointerWrapper(Int, user_data) + pw[element_id] = current_level return nothing end diff --git a/src/callbacks_step/amr_dg.jl b/src/callbacks_step/amr_dg.jl index 19bbebd925..1dcfdccdea 100644 --- a/src/callbacks_step/amr_dg.jl +++ b/src/callbacks_step/amr_dg.jl @@ -9,8 +9,7 @@ function rebalance_solver!(u_ode::AbstractVector, mesh::ParallelP4estMesh, equations, dg::DGSEM, cache, old_global_first_quadrant) # mpi ranks are 0-based, this array uses 1-based indices - global_first_quadrant = unsafe_wrap(Array, - unsafe_load(mesh.p4est).global_first_quadrant, + global_first_quadrant = unsafe_wrap(Array, mesh.p4est.global_first_quadrant, mpi_nranks() + 1) if global_first_quadrant[mpi_rank() + 1] == old_global_first_quadrant[mpi_rank() + 1] && diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl index 7fa2e21a24..8cf43a1d15 100644 --- a/src/callbacks_step/analysis.jl +++ b/src/callbacks_step/analysis.jl @@ -508,7 +508,7 @@ function print_amr_information(callbacks, mesh::P4estMesh, solver, cache) elements_per_level = zeros(P4EST_MAXLEVEL + 1) - for tree in unsafe_wrap_sc(p4est_tree_t, unsafe_load(mesh.p4est).trees) + for tree in unsafe_wrap_sc(p4est_tree_t, mesh.p4est.trees) elements_per_level .+= tree.quadrants_per_level end diff --git a/src/equations/acoustic_perturbation_2d.jl b/src/equations/acoustic_perturbation_2d.jl index 786630a14c..f4ce770e1e 100644 --- a/src/equations/acoustic_perturbation_2d.jl +++ b/src/equations/acoustic_perturbation_2d.jl @@ -145,7 +145,7 @@ function initial_condition_convergence_test(x, t, end """ - source_terms_convergence_test(u, x, t, equations::AcousticPerturbationEquations2D) + source_terms_convergence_test(u, x, t, equations::AcousticPerturbationEquations2D) Source terms used for convergence tests in combination with [`initial_condition_convergence_test`](@ref). diff --git a/src/equations/compressible_euler_1d.jl b/src/equations/compressible_euler_1d.jl index f484f26a58..15f7a2cb4c 100644 --- a/src/equations/compressible_euler_1d.jl +++ b/src/equations/compressible_euler_1d.jl @@ -374,7 +374,7 @@ Splitting of the compressible Euler flux of Steger and Warming. Returns a tuple of the fluxes "minus" (associated with waves going into the negative axis direction) and "plus" (associated with waves going into the positive axis direction). If only one of the fluxes is required, use the -function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}`. +function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}()`. !!! warning "Experimental implementation (upwind SBP)" This is an experimental feature and may change in future releases. @@ -462,7 +462,7 @@ it proved the most robust in practice. Returns a tuple of the fluxes "minus" (associated with waves going into the negative axis direction) and "plus" (associated with waves going into the positive axis direction). If only one of the fluxes is required, use the -function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}`. +function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}()`. !!! warning "Experimental implementation (upwind SBP)" This is an experimental feature and may change in future releases. @@ -555,7 +555,7 @@ are to handle flows at the low Mach number limit. Returns a tuple of the fluxes "minus" (associated with waves going into the negative axis direction) and "plus" (associated with waves going into the positive axis direction). If only one of the fluxes is required, use the -function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}`. +function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}()`. !!! warning "Experimental implementation (upwind SBP)" This is an experimental feature and may change in future releases. diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index 66e3c7bff8..05987c510b 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -31,7 +31,7 @@ The compressible Euler equations ``` for an ideal gas with ratio of specific heats `gamma` in two space dimensions. -Here, ``\rho`` is the density, ``v_1``,`v_2` the velocities, ``e`` the specific total energy **rather than** specific internal energy, and +Here, ``\rho`` is the density, ``v_1``, ``v_2`` the velocities, ``e`` the specific total energy **rather than** specific internal energy, and ```math p = (\gamma - 1) \left( \rho e - \frac{1}{2} \rho (v_1^2+v_2^2) \right) ``` @@ -694,7 +694,7 @@ Splitting of the compressible Euler flux of Steger and Warming. Returns a tuple of the fluxes "minus" (associated with waves going into the negative axis direction) and "plus" (associated with waves going into the positive axis direction). If only one of the fluxes is required, use the -function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}`. +function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}()`. !!! warning "Experimental implementation (upwind SBP)" This is an experimental feature and may change in future releases. @@ -826,7 +826,7 @@ it proved the most robust in practice. Returns a tuple of the fluxes "minus" (associated with waves going into the negative axis direction) and "plus" (associated with waves going into the positive axis direction). If only one of the fluxes is required, use the -function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}`. +function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}()`. !!! warning "Experimental implementation (upwind SBP)" This is an experimental feature and may change in future releases. @@ -924,7 +924,7 @@ to Burgers' equation. Returns a tuple of the fluxes "minus" (associated with waves going into the negative axis direction) and "plus" (associated with waves going into the positive axis direction). If only one of the fluxes is required, use the -function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}`. +function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}()`. !!! warning "Experimental implementation (upwind SBP)" This is an experimental feature and may change in future releases. diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl index c16a454b17..2085811f83 100644 --- a/src/equations/compressible_euler_3d.jl +++ b/src/equations/compressible_euler_3d.jl @@ -36,7 +36,7 @@ The compressible Euler equations ``` for an ideal gas with ratio of specific heats `gamma` in three space dimensions. -Here, ``\rho`` is the density, ``v_1``,`v_2`, `v_3` the velocities, ``e`` the specific total energy **rather than** specific internal energy, and +Here, ``\rho`` is the density, ``v_1``, ``v_2``, ``v_3`` the velocities, ``e`` the specific total energy **rather than** specific internal energy, and ```math p = (\gamma - 1) \left( \rho e - \frac{1}{2} \rho (v_1^2+v_2^2+v_3^2) \right) ``` @@ -770,7 +770,7 @@ Splitting of the compressible Euler flux of Steger and Warming. Returns a tuple of the fluxes "minus" (associated with waves going into the negative axis direction) and "plus" (associated with waves going into the positive axis direction). If only one of the fluxes is required, use the -function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}`. +function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}()`. !!! warning "Experimental implementation (upwind SBP)" This is an experimental feature and may change in future releases. diff --git a/src/equations/compressible_euler_multicomponent_1d.jl b/src/equations/compressible_euler_multicomponent_1d.jl index 4a50d60471..23ac222b97 100644 --- a/src/equations/compressible_euler_multicomponent_1d.jl +++ b/src/equations/compressible_euler_multicomponent_1d.jl @@ -44,8 +44,8 @@ specific heat capacity at constant volume of component ``i``. In case of more than one component, the specific heat ratios `gammas` and the gas constants `gas_constants` should be passed as tuples, e.g., `gammas=(1.4, 1.667)`. -The remaining variables like the specific heats at constant volume 'cv' or the specific heats at -constant pressure 'cp' are then calculated considering a calorically perfect gas. +The remaining variables like the specific heats at constant volume `cv` or the specific heats at +constant pressure `cp` are then calculated considering a calorically perfect gas. """ struct CompressibleEulerMulticomponentEquations1D{NVARS, NCOMP, RealT <: Real} <: AbstractCompressibleEulerMulticomponentEquations{1, NVARS, NCOMP} @@ -247,8 +247,8 @@ end Entropy conserving two-point flux by - Ayoub Gouasmi, Karthik Duraisamy (2020) - "Formulation of Entropy-Stable schemes for the multicomponent compressible Euler equations"" - arXiv:1904.00972v3 [math.NA] 4 Feb 2020 + "Formulation of Entropy-Stable schemes for the multicomponent compressible Euler equations" + [arXiv:1904.00972v3](https://arxiv.org/abs/1904.00972) [math.NA] 4 Feb 2020 """ @inline function flux_chandrashekar(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerMulticomponentEquations1D) diff --git a/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl index 5a015777cb..7b437f4a1b 100644 --- a/src/equations/compressible_euler_multicomponent_2d.jl +++ b/src/equations/compressible_euler_multicomponent_2d.jl @@ -48,8 +48,8 @@ specific heat capacity at constant volume of component ``i``. In case of more than one component, the specific heat ratios `gammas` and the gas constants `gas_constants` in [kJ/(kg*K)] should be passed as tuples, e.g., `gammas=(1.4, 1.667)`. -The remaining variables like the specific heats at constant volume 'cv' or the specific heats at -constant pressure 'cp' are then calculated considering a calorically perfect gas. +The remaining variables like the specific heats at constant volume `cv` or the specific heats at +constant pressure `cp` are then calculated considering a calorically perfect gas. """ struct CompressibleEulerMulticomponentEquations2D{NVARS, NCOMP, RealT <: Real} <: AbstractCompressibleEulerMulticomponentEquations{2, NVARS, NCOMP} @@ -275,8 +275,8 @@ end Adaption of the entropy conserving two-point flux by - Ayoub Gouasmi, Karthik Duraisamy (2020) - "Formulation of Entropy-Stable schemes for the multicomponent compressible Euler equations"" - arXiv:1904.00972v3 [math.NA] 4 Feb 2020 + "Formulation of Entropy-Stable schemes for the multicomponent compressible Euler equations" + [arXiv:1904.00972v3](https://arxiv.org/abs/1904.00972) [math.NA] 4 Feb 2020 """ @inline function flux_chandrashekar(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerMulticomponentEquations2D) diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index 33badba15d..9b06e0b5ab 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -73,8 +73,8 @@ where w_2 = \frac{\rho v_1}{p},\, w_3 = \frac{\rho v_2}{p},\, w_4 = -\frac{\rho}{p} ``` -#!!! warning "Experimental code" -# This code is experimental and may be changed or removed in any future release. +!!! warning "Experimental code" + This code is experimental and may be changed or removed in any future release. """ struct CompressibleNavierStokesDiffusion2D{GradientVariables, RealT <: Real, E <: AbstractCompressibleEulerEquations{2}} <: @@ -94,8 +94,8 @@ struct CompressibleNavierStokesDiffusion2D{GradientVariables, RealT <: Real, end """ -#!!! warning "Experimental code" -# This code is experimental and may be changed or removed in any future release. +!!! warning "Experimental code" + This code is experimental and may be changed or removed in any future release. `GradientVariablesPrimitive` and `GradientVariablesEntropy` are gradient variable type parameters for `CompressibleNavierStokesDiffusion2D`. By default, the gradient variables are set to be diff --git a/src/equations/compressible_navier_stokes_3d.jl b/src/equations/compressible_navier_stokes_3d.jl index 8930489295..0b770dff1c 100644 --- a/src/equations/compressible_navier_stokes_3d.jl +++ b/src/equations/compressible_navier_stokes_3d.jl @@ -73,8 +73,8 @@ where w_2 = \frac{\rho v_1}{p},\, w_3 = \frac{\rho v_2}{p},\, w_4 = \frac{\rho v_3}{p},\, w_5 = -\frac{\rho}{p} ``` -#!!! warning "Experimental code" -# This code is experimental and may be changed or removed in any future release. +!!! warning "Experimental code" + This code is experimental and may be changed or removed in any future release. """ struct CompressibleNavierStokesDiffusion3D{GradientVariables, RealT <: Real, E <: AbstractCompressibleEulerEquations{3}} <: diff --git a/src/equations/hyperbolic_diffusion_2d.jl b/src/equations/hyperbolic_diffusion_2d.jl index 25536a060f..511d1b8935 100644 --- a/src/equations/hyperbolic_diffusion_2d.jl +++ b/src/equations/hyperbolic_diffusion_2d.jl @@ -10,7 +10,7 @@ The linear hyperbolic diffusion equations in two space dimensions. A description of this system can be found in Sec. 2.5 of the book "I Do Like CFD, Too: Vol 1". -The book is freely available at http://www.cfdbooks.com/ and further analysis can be found in +The book is freely available at [http://www.cfdbooks.com/](http://www.cfdbooks.com/) and further analysis can be found in the paper by Nishikawa [DOI: 10.1016/j.jcp.2007.07.029](https://doi.org/10.1016/j.jcp.2007.07.029) """ struct HyperbolicDiffusionEquations2D{RealT <: Real} <: diff --git a/src/equations/hyperbolic_diffusion_3d.jl b/src/equations/hyperbolic_diffusion_3d.jl index bf6a00140d..ed807511b6 100644 --- a/src/equations/hyperbolic_diffusion_3d.jl +++ b/src/equations/hyperbolic_diffusion_3d.jl @@ -10,7 +10,7 @@ The linear hyperbolic diffusion equations in three space dimensions. A description of this system can be found in Sec. 2.5 of the book "I Do Like CFD, Too: Vol 1". -The book is freely available at http://www.cfdbooks.com/ and further analysis can be found in +The book is freely available at [http://www.cfdbooks.com/](http://www.cfdbooks.com/) and further analysis can be found in the paper by Nishikawa [DOI: 10.1016/j.jcp.2007.07.029](https://doi.org/10.1016/j.jcp.2007.07.029) """ struct HyperbolicDiffusionEquations3D{RealT <: Real} <: diff --git a/src/equations/ideal_glm_mhd_3d.jl b/src/equations/ideal_glm_mhd_3d.jl index 401fcd2daf..2e149d2849 100644 --- a/src/equations/ideal_glm_mhd_3d.jl +++ b/src/equations/ideal_glm_mhd_3d.jl @@ -41,7 +41,7 @@ end # Set initial conditions at physical location `x` for time `t` """ -initial_condition_constant(x, t, equations::IdealGlmMhdEquations3D) + initial_condition_constant(x, t, equations::IdealGlmMhdEquations3D) A constant initial condition to test free-stream preservation. """ diff --git a/src/equations/inviscid_burgers_1d.jl b/src/equations/inviscid_burgers_1d.jl index 8d4410b6ff..f2387f26ba 100644 --- a/src/equations/inviscid_burgers_1d.jl +++ b/src/equations/inviscid_burgers_1d.jl @@ -132,12 +132,12 @@ end equations::InviscidBurgersEquation1D) Naive local Lax-Friedrichs style flux splitting of the form `f⁺ = 0.5 (f + λ u)` -and `f⁻ = 0.5 (f - λ u)` where λ = abs(u). +and `f⁻ = 0.5 (f - λ u)` where `λ = abs(u)`. Returns a tuple of the fluxes "minus" (associated with waves going into the negative axis direction) and "plus" (associated with waves going into the positive axis direction). If only one of the fluxes is required, use the -function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}`. +function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}()`. !!! warning "Experimental implementation (upwind SBP)" This is an experimental feature and may change in future releases. diff --git a/src/equations/linear_scalar_advection_1d.jl b/src/equations/linear_scalar_advection_1d.jl index 7769cb61fb..6c6b9dd372 100644 --- a/src/equations/linear_scalar_advection_1d.jl +++ b/src/equations/linear_scalar_advection_1d.jl @@ -172,6 +172,44 @@ end return abs.(equation.advection_velocity) end +""" + splitting_lax_friedrichs(u, orientation::Integer, + equations::LinearScalarAdvectionEquation1D) + splitting_lax_friedrichs(u, which::Union{Val{:minus}, Val{:plus}} + orientation::Integer, + equations::LinearScalarAdvectionEquation1D) + +Naive local Lax-Friedrichs style flux splitting of the form `f⁺ = 0.5 (f + λ u)` +and `f⁻ = 0.5 (f - λ u)` where `λ` is the absolute value of the advection +velocity. + +Returns a tuple of the fluxes "minus" (associated with waves going into the +negative axis direction) and "plus" (associated with waves going into the +positive axis direction). If only one of the fluxes is required, use the +function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}()`. + +!!! warning "Experimental implementation (upwind SBP)" + This is an experimental feature and may change in future releases. +""" +@inline function splitting_lax_friedrichs(u, orientation::Integer, + equations::LinearScalarAdvectionEquation1D) + fm = splitting_lax_friedrichs(u, Val{:minus}(), orientation, equations) + fp = splitting_lax_friedrichs(u, Val{:plus}(), orientation, equations) + return fm, fp +end + +@inline function splitting_lax_friedrichs(u, ::Val{:plus}, orientation::Integer, + equations::LinearScalarAdvectionEquation1D) + a = equations.advection_velocity[1] + return a > 0 ? flux(u, orientation, equations) : zero(u) +end + +@inline function splitting_lax_friedrichs(u, ::Val{:minus}, orientation::Integer, + equations::LinearScalarAdvectionEquation1D) + a = equations.advection_velocity[1] + return a < 0 ? flux(u, orientation, equations) : zero(u) +end + # Convert conservative variables to primitive @inline cons2prim(u, equation::LinearScalarAdvectionEquation1D) = u diff --git a/src/equations/shallow_water_two_layer_1d.jl b/src/equations/shallow_water_two_layer_1d.jl index edf7d5e32f..e126eec7c2 100644 --- a/src/equations/shallow_water_two_layer_1d.jl +++ b/src/equations/shallow_water_two_layer_1d.jl @@ -11,28 +11,28 @@ Two-Layer Shallow Water equations (2LSWE) in one space dimension. The equations are given by ```math \begin{alignat*}{4} -&\frac{\partial}{\partial t}h_{upper} -&&+ \frac{\partial}{\partial x}\left(h_{upper} v_{1,upper}\right) +&\frac{\partial}{\partial t}h_{upper} +&&+ \frac{\partial}{\partial x}\left(h_{upper} v_{1,upper}\right) &&= 0 \\ -&\frac{\partial}{\partial t}\left(h_{upper}v_{1,upper}\right) -&&+ \frac{\partial}{\partial x}\left(h_{upper}v_{1,upper}^2 + \dfrac{gh_{upper}^2}{2}\right) +&\frac{\partial}{\partial t}\left(h_{upper}v_{1,upper}\right) +&&+ \frac{\partial}{\partial x}\left(h_{upper}v_{1,upper}^2 + \dfrac{gh_{upper}^2}{2}\right) &&= -gh_{upper}\frac{\partial}{\partial x}\left(b+h_{lower}\right)\\ -&\frac{\partial}{\partial t}h_{lower} -&&+ \frac{\partial}{\partial x}\left(h_{lower}v_{1,lower}\right) +&\frac{\partial}{\partial t}h_{lower} +&&+ \frac{\partial}{\partial x}\left(h_{lower}v_{1,lower}\right) &&= 0 \\ -&\frac{\partial}{\partial t}\left(h_{lower}v_{1,lower}\right) -&&+ \frac{\partial}{\partial x}\left(h_{lower}v_{1,lower}^2 + \dfrac{gh_{lower}^2}{2}\right) +&\frac{\partial}{\partial t}\left(h_{lower}v_{1,lower}\right) +&&+ \frac{\partial}{\partial x}\left(h_{lower}v_{1,lower}^2 + \dfrac{gh_{lower}^2}{2}\right) &&= -gh_{lower}\frac{\partial}{\partial x}\left(b+\dfrac{\rho_{upper}}{\rho_{lower}}h_{upper}\right). \end{alignat*} ``` -The unknown quantities of the 2LSWE are the water heights of the {lower} layer ``h_{lower}`` and the -{upper} layer ``h_{upper}`` with respective velocities ``v_{1,upper}`` and ``v_{1,lower}``. The gravitational constant is -denoted by `g`, the layer densitites by ``\rho_{upper}``and ``\rho_{lower}`` and the (possibly) variable -bottom topography function ``b(x)``. The conservative variable water height ``h_{lower}`` is measured -from the bottom topography ``b`` and ``h_{upper}`` relative to ``h_{lower}``, therefore one also defines the +The unknown quantities of the 2LSWE are the water heights of the {lower} layer ``h_{lower}`` and the +{upper} layer ``h_{upper}`` with respective velocities ``v_{1,upper}`` and ``v_{1,lower}``. The gravitational constant is +denoted by `g`, the layer densitites by ``\rho_{upper}``and ``\rho_{lower}`` and the (possibly) variable +bottom topography function ``b(x)``. The conservative variable water height ``h_{lower}`` is measured +from the bottom topography ``b`` and ``h_{upper}`` relative to ``h_{lower}``, therefore one also defines the total water heights as ``H_{upper} = h_{upper} + h_{upper} + b`` and ``H_{lower} = h_{lower} + b``. -The densities must be chosen such that ``\rho_{upper} < \rho_{lower}``, to make sure that the heavier fluid +The densities must be chosen such that ``\rho_{upper} < \rho_{lower}``, to make sure that the heavier fluid ``\rho_{lower}`` is in the bottom layer and the lighter fluid ``\rho_{upper}`` in the {upper} layer. The additional quantity ``H_0`` is also available to store a reference value for the total water @@ -41,13 +41,13 @@ height that is useful to set initial conditions or test the "lake-at-rest" well- The bottom topography function ``b(x)`` is set inside the initial condition routine for a particular problem setup. -In addition to the unknowns, Trixi currently stores the bottom topography values at the -approximation points despite being fixed in time. This is done for convenience of computing the -bottom topography gradients on the fly during the approximation as well as computing auxiliary +In addition to the unknowns, Trixi currently stores the bottom topography values at the +approximation points despite being fixed in time. This is done for convenience of computing the +bottom topography gradients on the fly during the approximation as well as computing auxiliary quantities like the total water height ``H`` or the entropy variables. This affects the implementation and use of these equations in various ways: * The flux values corresponding to the bottom topography must be zero. -* The bottom topography values must be included when defining initial conditions, boundary +* The bottom topography values must be included when defining initial conditions, boundary conditions or source terms. * [`AnalysisCallback`](@ref) analyzes this variable. * Trixi's visualization tools will visualize the bottom topography by default. @@ -101,7 +101,7 @@ end initial_condition_convergence_test(x, t, equations::ShallowWaterTwoLayerEquations1D) A smooth initial condition used for convergence tests in combination with -[`source_terms_convergence_test`](@ref) (and +[`source_terms_convergence_test`](@ref) (and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains). """ function initial_condition_convergence_test(x, t, @@ -121,9 +121,9 @@ end """ source_terms_convergence_test(u, x, t, equations::ShallowWaterTwoLayerEquations1D) -Source terms used for convergence tests in combination with -[`initial_condition_convergence_test`](@ref) -(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) +Source terms used for convergence tests in combination with +[`initial_condition_convergence_test`](@ref) +(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains). """ @inline function source_terms_convergence_test(u, x, t, @@ -167,8 +167,8 @@ the internal value. For details see Section 9.2.5 of the book: - Eleuterio F. Toro (2001) - Shock-Capturing Methods for Free-Surface Shallow Flows - 1st edition + Shock-Capturing Methods for Free-Surface Shallow Flows + 1st edition ISBN 0471987662 """ @inline function boundary_condition_slip_wall(u_inner, orientation_or_normal, direction, @@ -219,7 +219,7 @@ end Non-symmetric two-point volume flux discretizing the nonconservative (source) term that contains the gradient of the bottom topography [`ShallowWaterTwoLayerEquations2D`](@ref) and an -additional term that couples the momentum of both layers. This is a slightly modified version +additional term that couples the momentum of both layers. This is a slightly modified version to account for the additional source term compared to the standard SWE described in the paper. Further details are available in the paper: @@ -238,7 +238,7 @@ Further details are available in the paper: z = zero(eltype(u_ll)) - # Bottom gradient nonconservative term: (0, g*h_upper*(b+h_lower)_x, + # Bottom gradient nonconservative term: (0, g*h_upper*(b+h_lower)_x, # 0, g*h_lower*(b+r*h_upper)_x, 0) f = SVector(z, equations.gravity * h_upper_ll * (b_rr + h_lower_rr), @@ -254,9 +254,9 @@ end !!! warning "Experimental code" This numerical flux is experimental and may change in any future release. - -Non-symmetric two-point surface flux discretizing the nonconservative (source) term that contains -the gradients of the bottom topography and an additional term that couples the momentum of both + +Non-symmetric two-point surface flux discretizing the nonconservative (source) term that contains +the gradients of the bottom topography and an additional term that couples the momentum of both layers [`ShallowWaterTwoLayerEquations2D`](@ref). Further details are available in the paper: @@ -286,7 +286,7 @@ formulation. z = zero(eltype(u_ll)) - # Bottom gradient nonconservative term: (0, g*h_upper*(b+h_lower)_x, + # Bottom gradient nonconservative term: (0, g*h_upper*(b+h_lower)_x, # 0, g*h_lower*(b+r*h_upper)_x, 0) f = SVector(z, g * h_upper_ll * (b_ll + h_lower_ll) + @@ -303,13 +303,13 @@ end flux_fjordholm_etal(u_ll, u_rr, orientation, equations::ShallowWaterTwoLayerEquations1D) -Total energy conservative (mathematical entropy for shallow water equations). When the bottom -topography is nonzero this should only be used as a surface flux otherwise the scheme will not be +Total energy conservative (mathematical entropy for shallow water equations). When the bottom +topography is nonzero this should only be used as a surface flux otherwise the scheme will not be well-balanced. For well-balancedness in the volume flux use [`flux_wintermeyer_etal`](@ref). Details are available in Eq. (4.1) in the paper: - Ulrik S. Fjordholm, Siddhartha Mishra and Eitan Tadmor (2011) - Well-balanced and energy stable schemes for the shallow water equations with discontinuous + Well-balanced and energy stable schemes for the shallow water equations with discontinuous topography [DOI: 10.1016/j.jcp.2011.03.042](https://doi.org/10.1016/j.jcp.2011.03.042) and the application to two layers is shown in the paper: - Ulrik Skre Fjordholm (2012) @@ -348,11 +348,11 @@ end """ flux_wintermeyer_etal(u_ll, u_rr, orientation, equations::ShallowWaterTwoLayerEquations1D) - + Total energy conservative (mathematical entropy for two-layer shallow water equations) split form. When the bottom topography is nonzero this scheme will be well-balanced when used as a `volume_flux`. The `surface_flux` should still use, e.g., [`flux_fjordholm_etal`](@ref). To obtain the flux for the -two-layer shallow water equations the flux that is described in the paper for the normal shallow +two-layer shallow water equations the flux that is described in the paper for the normal shallow water equations is used within each layer. Further details are available in Theorem 1 of the paper: @@ -391,8 +391,8 @@ end flux_es_fjordholm_etal(u_ll, u_rr, orientation, equations::ShallowWaterTwoLayerEquations1D) -Entropy stable surface flux for the two-layer shallow water equations. Uses the entropy -conservative flux_fjordholm_etal and adds a Lax-Friedrichs type dissipation dependent on the jump +Entropy stable surface flux for the two-layer shallow water equations. Uses the entropy +conservative [`flux_fjordholm_etal`](@ref) and adds a Lax-Friedrichs type dissipation dependent on the jump of entropy variables. Further details are available in the paper: @@ -460,10 +460,10 @@ formulation. end # Calculate approximation for maximum wave speed for local Lax-Friedrichs-type dissipation as the -# maximum velocity magnitude plus the maximum speed of sound. This function uses approximate -# eigenvalues using the speed of the barotropic mode as there is no simple way to calculate them +# maximum velocity magnitude plus the maximum speed of sound. This function uses approximate +# eigenvalues using the speed of the barotropic mode as there is no simple way to calculate them # analytically. -# +# # A good overview of the derivation is given in: # - Jonas Nycander, Andrew McC. Hogg, Leela M. Frankcombe (2008) # Open boundary conditions for nonlinear channel Flows @@ -488,7 +488,7 @@ end return (max(abs(v_m_ll) + c_ll, abs(v_m_rr) + c_rr)) end -# Specialized `DissipationLocalLaxFriedrichs` to avoid spurious dissipation in the bottom +# Specialized `DissipationLocalLaxFriedrichs` to avoid spurious dissipation in the bottom # topography @inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr, orientation_or_normal_direction, @@ -530,7 +530,7 @@ end end # Convert conservative variables to entropy variables -# Note, only the first four are the entropy variables, the fifth entry still just carries the +# Note, only the first four are the entropy variables, the fifth entry still just carries the # bottom topography values for convenience @inline function cons2entropy(u, equations::ShallowWaterTwoLayerEquations1D) h_upper, _, h_lower, _, b = u @@ -567,7 +567,7 @@ end # Calculate total energy for a conservative state `cons` @inline function energy_total(cons, equations::ShallowWaterTwoLayerEquations1D) - h_upper, h_lower, h_v1_upper, h_v2_lower, b = cons + h_upper, h_v1_upper, h_lower, h_v2_lower, b = cons # Set new variables for better readability g = equations.gravity rho_upper = equations.rho_upper diff --git a/src/equations/shallow_water_two_layer_2d.jl b/src/equations/shallow_water_two_layer_2d.jl index b5e52d636e..a54831c711 100644 --- a/src/equations/shallow_water_two_layer_2d.jl +++ b/src/equations/shallow_water_two_layer_2d.jl @@ -695,8 +695,9 @@ end """ flux_es_fjordholm_etal(u_ll, u_rr, orientation_or_normal_direction, equations::ShallowWaterTwoLayerEquations1D) + Entropy stable surface flux for the two-layer shallow water equations. Uses the entropy conservative -flux_fjordholm_etal and adds a Lax-Friedrichs type dissipation dependent on the jump of entropy +[`flux_fjordholm_etal`](@ref) and adds a Lax-Friedrichs type dissipation dependent on the jump of entropy variables. Further details are available in the paper: diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index ddd6cf473e..60db285e04 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -13,9 +13,9 @@ to manage trees and mesh refinement. """ mutable struct P4estMesh{NDIMS, RealT <: Real, IsParallel, P, Ghost, NDIMSP2, NNODES} <: AbstractMesh{NDIMS} - p4est::P # Either Ptr{p4est_t} or Ptr{p8est_t} - is_parallel::IsParallel - ghost::Ghost # Either Ptr{p4est_ghost_t} or Ptr{p8est_ghost_t} + p4est :: P # Either PointerWrapper{p4est_t} or PointerWrapper{p8est_t} + is_parallel :: IsParallel + ghost :: Ghost # Either PointerWrapper{p4est_ghost_t} or PointerWrapper{p8est_ghost_t} # Coordinates at the nodes specified by the tensor product of `nodes` (NDIMS times). # This specifies the geometry interpolation for each tree. tree_node_coordinates::Array{RealT, NDIMSP2} # [dimension, i, j, k, tree] @@ -43,18 +43,21 @@ mutable struct P4estMesh{NDIMS, RealT <: Real, IsParallel, P, Ghost, NDIMSP2, NN is_parallel = False() end + p4est_pw = PointerWrapper(p4est) + ghost = ghost_new_p4est(p4est) + ghost_pw = PointerWrapper(ghost) mesh = new{NDIMS, eltype(tree_node_coordinates), typeof(is_parallel), - typeof(p4est), typeof(ghost), NDIMS + 2, length(nodes)}(p4est, - is_parallel, - ghost, - tree_node_coordinates, - nodes, - boundary_names, - current_filename, - unsaved_changes, - p4est_partition_allow_for_coarsening) + typeof(p4est_pw), typeof(ghost_pw), NDIMS + 2, length(nodes)}(p4est_pw, + is_parallel, + ghost_pw, + tree_node_coordinates, + nodes, + boundary_names, + current_filename, + unsaved_changes, + p4est_partition_allow_for_coarsening) # Destroy `p4est` structs when the mesh is garbage collected finalizer(destroy_mesh, mesh) @@ -70,14 +73,14 @@ const ParallelP4estMesh{NDIMS} = P4estMesh{NDIMS, <:Real, <:True} @inline mpi_parallel(mesh::ParallelP4estMesh) = True() function destroy_mesh(mesh::P4estMesh{2}) - connectivity = unsafe_load(mesh.p4est).connectivity + connectivity = mesh.p4est.connectivity p4est_ghost_destroy(mesh.ghost) p4est_destroy(mesh.p4est) p4est_connectivity_destroy(connectivity) end function destroy_mesh(mesh::P4estMesh{3}) - connectivity = unsafe_load(mesh.p4est).connectivity + connectivity = mesh.p4est.connectivity p8est_ghost_destroy(mesh.ghost) p8est_destroy(mesh.p4est) p8est_connectivity_destroy(connectivity) @@ -87,11 +90,10 @@ end @inline Base.real(::P4estMesh{NDIMS, RealT}) where {NDIMS, RealT} = RealT @inline function ntrees(mesh::P4estMesh) - trees = unsafe_load(mesh.p4est).trees - return unsafe_load(trees).elem_count + return mesh.p4est.trees.elem_count[] end # returns Int32 by default which causes a weird method error when creating the cache -@inline ncells(mesh::P4estMesh) = Int(unsafe_load(mesh.p4est).local_num_quadrants) +@inline ncells(mesh::P4estMesh) = Int(mesh.p4est.local_num_quadrants[]) function Base.show(io::IO, mesh::P4estMesh) print(io, "P4estMesh{", ndims(mesh), ", ", real(mesh), "}") @@ -387,14 +389,14 @@ function p4est_mesh_from_hohqmesh_abaqus(meshfile, initial_refinement_level, n_dimensions, RealT) # Create the mesh connectivity using `p4est` connectivity = read_inp_p4est(meshfile, Val(n_dimensions)) - connectivity_obj = unsafe_load(connectivity) + connectivity_pw = PointerWrapper(connectivity) # These need to be of the type Int for unsafe_wrap below to work - n_trees::Int = connectivity_obj.num_trees - n_vertices::Int = connectivity_obj.num_vertices + n_trees::Int = connectivity_pw.num_trees[] + n_vertices::Int = connectivity_pw.num_vertices[] # Extract a copy of the element vertices to compute the tree node coordinates - vertices = unsafe_wrap(Array, connectivity_obj.vertices, (3, n_vertices)) + vertices = unsafe_wrap(Array, connectivity_pw.vertices, (3, n_vertices)) # Readin all the information from the mesh file into a string array file_lines = readlines(open(meshfile)) @@ -445,14 +447,14 @@ function p4est_mesh_from_standard_abaqus(meshfile, mapping, polydeg, initial_refinement_level, n_dimensions, RealT) # Create the mesh connectivity using `p4est` connectivity = read_inp_p4est(meshfile, Val(n_dimensions)) - connectivity_obj = unsafe_load(connectivity) + connectivity_pw = PointerWrapper(connectivity) # These need to be of the type Int for unsafe_wrap below to work - n_trees::Int = connectivity_obj.num_trees - n_vertices::Int = connectivity_obj.num_vertices + n_trees::Int = connectivity_pw.num_trees[] + n_vertices::Int = connectivity_pw.num_vertices[] - vertices = unsafe_wrap(Array, connectivity_obj.vertices, (3, n_vertices)) - tree_to_vertex = unsafe_wrap(Array, connectivity_obj.tree_to_vertex, + vertices = unsafe_wrap(Array, connectivity_pw.vertices, (3, n_vertices)) + tree_to_vertex = unsafe_wrap(Array, connectivity_pw.tree_to_vertex, (2^n_dimensions, n_trees)) basis = LobattoLegendreBasis(RealT, polydeg) @@ -1511,17 +1513,18 @@ end function update_ghost_layer!(mesh::P4estMesh) ghost_destroy_p4est(mesh.ghost) - mesh.ghost = ghost_new_p4est(mesh.p4est) + mesh.ghost = PointerWrapper(ghost_new_p4est(mesh.p4est)) end function init_fn(p4est, which_tree, quadrant) # Unpack quadrant's user data ([global quad ID, controller_value]) - ptr = Ptr{Int}(unsafe_load(quadrant.p.user_data)) + # Use `unsafe_load` here since `quadrant.p.user_data isa Ptr{Ptr{Nothing}}` + # and we only need the first (only!) entry + pw = PointerWrapper(Int, unsafe_load(quadrant.p.user_data)) # Initialize quad ID as -1 and controller_value as 0 (don't refine or coarsen) - unsafe_store!(ptr, -1, 1) - unsafe_store!(ptr, 0, 2) - + pw[1] = -1 + pw[2] = 0 return nothing end @@ -1539,8 +1542,10 @@ end function refine_fn(p4est, which_tree, quadrant) # Controller value has been copied to the quadrant's user data storage before. # Unpack quadrant's user data ([global quad ID, controller_value]). - ptr = Ptr{Int}(unsafe_load(quadrant.p.user_data)) - controller_value = unsafe_load(ptr, 2) + # Use `unsafe_load` here since `quadrant.p.user_data isa Ptr{Ptr{Nothing}}` + # and we only need the first (only!) entry + pw = PointerWrapper(Int, unsafe_load(quadrant.p.user_data)) + controller_value = pw[2] if controller_value > 0 # return true (refine) @@ -1586,9 +1591,9 @@ function coarsen_fn(p4est, which_tree, quadrants_ptr) # Controller value has been copied to the quadrant's user data storage before. # Load controller value from quadrant's user data ([global quad ID, controller_value]). - function controller_value(i) - unsafe_load(Ptr{Int}(unsafe_load(quadrants[i].p.user_data)), 2) - end + # Use `unsafe_load` here since `quadrant.p.user_data isa Ptr{Ptr{Nothing}}` + # and we only need the first (only!) entry + controller_value(i) = PointerWrapper(Int, unsafe_load(quadrants[i].p.user_data))[2] # `p4est` calls this function for each 2^ndims quads that could be coarsened to a single one. # Only coarsen if all these 2^ndims quads have been marked for coarsening. @@ -1671,20 +1676,19 @@ end # Copy global quad ID to quad's user data storage, will be called below function save_original_id_iter_volume(info, user_data) - info_obj = unsafe_load(info) + info_pw = PointerWrapper(info) # Load tree from global trees array, one-based indexing - tree = unsafe_load_tree(info_obj.p4est, info_obj.treeid + 1) + tree_pw = load_pointerwrapper_tree(info_pw.p4est, info_pw.treeid[] + 1) # Quadrant numbering offset of this quadrant - offset = tree.quadrants_offset + offset = tree_pw.quadrants_offset[] # Global quad ID - quad_id = offset + info_obj.quadid + quad_id = offset + info_pw.quadid[] # Unpack quadrant's user data ([global quad ID, controller_value]) - ptr = Ptr{Int}(unsafe_load(info_obj.quad.p.user_data)) + pw = PointerWrapper(Int, info_pw.quad.p.user_data[]) # Save global quad ID - unsafe_store!(ptr, quad_id, 1) - + pw[1] = quad_id return nothing end @@ -1708,24 +1712,23 @@ end # Extract information about which cells have been changed function collect_changed_iter_volume(info, user_data) - info_obj = unsafe_load(info) + info_pw = PointerWrapper(info) # The original element ID has been saved to user_data before. # Load original quad ID from quad's user data ([global quad ID, controller_value]). - quad_data_ptr = Ptr{Int}(unsafe_load(info_obj.quad.p.user_data)) - original_id = unsafe_load(quad_data_ptr, 1) + quad_data_pw = PointerWrapper(Int, info_pw.quad.p.user_data[]) + original_id = quad_data_pw[1] # original_id of cells that have been newly created is -1 if original_id >= 0 # Unpack user_data = original_cells - user_data_ptr = Ptr{Int}(user_data) + user_data_pw = PointerWrapper(Int, user_data) # If quad has an original_id, it existed before refinement/coarsening, # and therefore wasn't changed. # Mark original_id as "not changed during refinement/coarsening" in original_cells - unsafe_store!(user_data_ptr, 0, original_id + 1) + user_data_pw[original_id + 1] = 0 end - return nothing end @@ -1756,29 +1759,27 @@ end # Extract newly created cells function collect_new_iter_volume(info, user_data) - info_obj = unsafe_load(info) + info_pw = PointerWrapper(info) # The original element ID has been saved to user_data before. # Unpack quadrant's user data ([global quad ID, controller_value]). - quad_data_ptr = Ptr{Int}(unsafe_load(info_obj.quad.p.user_data)) - original_id = unsafe_load(quad_data_ptr, 1) + original_id = PointerWrapper(Int, info_pw.quad.p.user_data[])[1] # original_id of cells that have been newly created is -1 if original_id < 0 # Load tree from global trees array, one-based indexing - tree = unsafe_load_tree(info_obj.p4est, info_obj.treeid + 1) + tree_pw = load_pointerwrapper_tree(info_pw.p4est, info_pw.treeid[] + 1) # Quadrant numbering offset of this quadrant - offset = tree.quadrants_offset + offset = tree_pw.quadrants_offset[] # Global quad ID - quad_id = offset + info_obj.quadid + quad_id = offset + info_pw.quadid[] # Unpack user_data = original_cells - user_data_ptr = Ptr{Int}(user_data) + user_data_pw = PointerWrapper(Int, user_data) # Mark cell as "newly created during refinement/coarsening/balancing" - unsafe_store!(user_data_ptr, 1, quad_id + 1) + user_data_pw[quad_id + 1] = 1 end - return nothing end diff --git a/src/meshes/structured_mesh.jl b/src/meshes/structured_mesh.jl index 5872681933..df067db833 100644 --- a/src/meshes/structured_mesh.jl +++ b/src/meshes/structured_mesh.jl @@ -33,7 +33,7 @@ Create a StructuredMesh of the given size and shape that uses `RealT` as coordin the reference mesh to the physical domain. If no `mapping_as_string` is defined, this function must be defined with the name `mapping` to allow for restarts. - This will be changed in the future, see https://github.com/trixi-framework/Trixi.jl/issues/541. + This will be changed in the future, see [https://github.com/trixi-framework/Trixi.jl/issues/541](https://github.com/trixi-framework/Trixi.jl/issues/541). - `RealT::Type`: the type that should be used for coordinates. - `periodicity`: either a `Bool` deciding if all of the boundaries are periodic or an `NTuple{NDIMS, Bool}` deciding for each dimension if the boundaries in this dimension are periodic. @@ -41,7 +41,7 @@ Create a StructuredMesh of the given size and shape that uses `RealT` as coordin - `mapping_as_string::String`: the code that defines the `mapping`. If `CodeTracking` can't find the function definition, it can be passed directly here. The code string must define the mapping function with the name `mapping`. - This will be changed in the future, see https://github.com/trixi-framework/Trixi.jl/issues/541. + This will be changed in the future, see [https://github.com/trixi-framework/Trixi.jl/issues/541](https://github.com/trixi-framework/Trixi.jl/issues/541). """ function StructuredMesh(cells_per_dimension, mapping; RealT = Float64, periodicity = true, unsaved_changes = true, diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl index 838fa2d581..2536cfe0bf 100644 --- a/src/solvers/dg.jl +++ b/src/solvers/dg.jl @@ -628,6 +628,15 @@ function compute_coefficients!(u, func, t, mesh::AbstractMesh{1}, equations, dg: for i in eachnode(dg) x_node = get_node_coords(cache.elements.node_coordinates, equations, dg, i, element) + # Changing the node positions passed to the initial condition by the minimum + # amount possible with the current type of floating point numbers allows setting + # discontinuous initial data in a simple way. In particular, a check like `if x < x_jump` + # works if the jump location `x_jump` is at the position of an interface. + if i == 1 + x_node = SVector(nextfloat(x_node[1])) + elseif i == nnodes(dg) + x_node = SVector(prevfloat(x_node[1])) + end u_node = func(x_node, t, equations) set_node_vars!(u, u_node, equations, dg, i, element) end diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl index d51c7cabf9..bc76aa1a9d 100644 --- a/src/solvers/dgmulti/dg.jl +++ b/src/solvers/dgmulti/dg.jl @@ -226,6 +226,11 @@ function estimate_dt(mesh::DGMultiMesh, dg::DGMulti) return StartUpDG.estimate_h(rd, mesh.md) / StartUpDG.inverse_trace_constant(rd) end +dt_polydeg_scaling(dg::DGMulti) = inv(dg.basis.N + 1) +function dt_polydeg_scaling(dg::DGMulti{3, <:Wedge, <:TensorProductWedge}) + inv(maximum(dg.basis.N) + 1) +end + # for the stepsize callback function max_dt(u, t, mesh::DGMultiMesh, constant_speed::False, equations, dg::DGMulti{NDIMS}, @@ -247,8 +252,7 @@ function max_dt(u, t, mesh::DGMultiMesh, # `polydeg+1`. This is because `nnodes(dg)` returns the total number of # multi-dimensional nodes for DGMulti solver types, while `nnodes(dg)` returns # the number of 1D nodes for `DGSEM` solvers. - polydeg = rd.N - return 2 * dt_min / (polydeg + 1) + return 2 * dt_min * dt_polydeg_scaling(dg) end function max_dt(u, t, mesh::DGMultiMesh, @@ -270,8 +274,7 @@ function max_dt(u, t, mesh::DGMultiMesh, # `polydeg+1`. This is because `nnodes(dg)` returns the total number of # multi-dimensional nodes for DGMulti solver types, while `nnodes(dg)` returns # the number of 1D nodes for `DGSEM` solvers. - polydeg = rd.N - return 2 * dt_min / (polydeg + 1) + return 2 * dt_min * dt_polydeg_scaling(dg) end # interpolates from solution coefficients to face quadrature points diff --git a/src/solvers/dgmulti/types.jl b/src/solvers/dgmulti/types.jl index c225e334e8..fe6510856b 100644 --- a/src/solvers/dgmulti/types.jl +++ b/src/solvers/dgmulti/types.jl @@ -96,6 +96,21 @@ function DGMulti(; polydeg = nothing, polydeg = polydeg, kwargs...) end +# dispatchable constructor for DGMulti using a TensorProductWedge +function DGMulti(element_type::Wedge, + approximation_type, + volume_integral, + surface_integral; + polydeg::Tuple, + kwargs...) + factor_a = RefElemData(Tri(), approximation_type, polydeg[1]; kwargs...) + factor_b = RefElemData(Line(), approximation_type, polydeg[2]; kwargs...) + + tensor = TensorProductWedge(factor_a, factor_b) + rd = RefElemData(element_type, tensor; kwargs...) + return DG(rd, nothing, surface_integral, volume_integral) +end + # dispatchable constructor for DGMulti to allow for specialization function DGMulti(element_type::AbstractElemShape, approximation_type, @@ -165,9 +180,9 @@ GeometricTermsType(mesh_type::Curved, element_type::AbstractElemShape) = NonAffi # other potential mesh types to add later: Polynomial{polydeg_geo}? """ - DGMultiMesh(dg::DGMulti{NDIMS}, vertex_coordinates, EToV; - is_on_boundary=nothing, - periodicity=ntuple(_->false, NDIMS)) where {NDIMS} + DGMultiMesh(dg::DGMulti{NDIMS}, vertex_coordinates, EToV; + is_on_boundary=nothing, + periodicity=ntuple(_->false, NDIMS)) where {NDIMS} - `dg::DGMulti` contains information associated with to the reference element (e.g., quadrature, basis evaluation, differentiation, etc). diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index 9b87de777a..2b9c6987d2 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -276,18 +276,18 @@ function init_boundaries!(boundaries, mesh::P4estMesh) end # Function barrier for type stability -function init_boundaries_iter_face_inner(info, boundaries, boundary_id, mesh) +function init_boundaries_iter_face_inner(info_pw, boundaries, boundary_id, mesh) # Extract boundary data - side = unsafe_load_side(info) + side_pw = load_pointerwrapper_side(info_pw) # Get local tree, one-based indexing - tree = unsafe_load_tree(mesh.p4est, side.treeid + 1) + tree_pw = load_pointerwrapper_tree(mesh.p4est, side_pw.treeid[] + 1) # Quadrant numbering offset of this quadrant - offset = tree.quadrants_offset + offset = tree_pw.quadrants_offset[] # Verify before accessing is.full, but this should never happen - @assert side.is_hanging == false + @assert side_pw.is_hanging[] == false - local_quad_id = side.is.full.quadid + local_quad_id = side_pw.is.full.quadid[] # Global ID of this quad quad_id = offset + local_quad_id @@ -296,13 +296,13 @@ function init_boundaries_iter_face_inner(info, boundaries, boundary_id, mesh) boundaries.neighbor_ids[boundary_id] = quad_id + 1 # Face at which the boundary lies - face = side.face + face = side_pw.face[] # Save boundaries.node_indices dimension specific in containers_[23]d.jl init_boundary_node_indices!(boundaries, face, boundary_id) # One-based indexing - boundaries.name[boundary_id] = mesh.boundary_names[face + 1, side.treeid + 1] + boundaries.name[boundary_id] = mesh.boundary_names[face + 1, side_pw.treeid[] + 1] return nothing end @@ -479,32 +479,33 @@ end # Function barrier for type stability function init_surfaces_iter_face_inner(info, user_data) @unpack interfaces, mortars, boundaries = user_data - elem_count = unsafe_load(info).sides.elem_count + info_pw = PointerWrapper(info) + elem_count = info_pw.sides.elem_count[] if elem_count == 2 # Two neighboring elements => Interface or mortar # Extract surface data - sides = (unsafe_load_side(info, 1), unsafe_load_side(info, 2)) + sides_pw = (load_pointerwrapper_side(info_pw, 1), + load_pointerwrapper_side(info_pw, 2)) - if sides[1].is_hanging == false && sides[2].is_hanging == false + if sides_pw[1].is_hanging[] == false && sides_pw[2].is_hanging[] == false # No hanging nodes => normal interface if interfaces !== nothing - init_interfaces_iter_face_inner(info, sides, user_data) + init_interfaces_iter_face_inner(info_pw, sides_pw, user_data) end else # Hanging nodes => mortar if mortars !== nothing - init_mortars_iter_face_inner(info, sides, user_data) + init_mortars_iter_face_inner(info_pw, sides_pw, user_data) end end elseif elem_count == 1 # One neighboring elements => boundary if boundaries !== nothing - init_boundaries_iter_face_inner(info, user_data) + init_boundaries_iter_face_inner(info_pw, user_data) end end - return nothing end @@ -519,18 +520,18 @@ function init_surfaces!(interfaces, mortars, boundaries, mesh::P4estMesh) end # Initialization of interfaces after the function barrier -function init_interfaces_iter_face_inner(info, sides, user_data) +function init_interfaces_iter_face_inner(info_pw, sides_pw, user_data) @unpack interfaces, interface_id, mesh = user_data user_data.interface_id += 1 # Get Tuple of local trees, one-based indexing - trees = (unsafe_load_tree(mesh.p4est, sides[1].treeid + 1), - unsafe_load_tree(mesh.p4est, sides[2].treeid + 1)) + trees_pw = (load_pointerwrapper_tree(mesh.p4est, sides_pw[1].treeid[] + 1), + load_pointerwrapper_tree(mesh.p4est, sides_pw[2].treeid[] + 1)) # Quadrant numbering offsets of the quadrants at this interface - offsets = SVector(trees[1].quadrants_offset, - trees[2].quadrants_offset) + offsets = SVector(trees_pw[1].quadrants_offset[], + trees_pw[2].quadrants_offset[]) - local_quad_ids = SVector(sides[1].is.full.quadid, sides[2].is.full.quadid) + local_quad_ids = SVector(sides_pw[1].is.full.quadid[], sides_pw[2].is.full.quadid[]) # Global IDs of the neighboring quads quad_ids = offsets + local_quad_ids @@ -540,31 +541,30 @@ function init_interfaces_iter_face_inner(info, sides, user_data) interfaces.neighbor_ids[2, interface_id] = quad_ids[2] + 1 # Face at which the interface lies - faces = (sides[1].face, sides[2].face) + faces = (sides_pw[1].face[], sides_pw[2].face[]) # Save interfaces.node_indices dimension specific in containers_[23]d.jl - init_interface_node_indices!(interfaces, faces, - unsafe_load(info).orientation, interface_id) + init_interface_node_indices!(interfaces, faces, info_pw.orientation[], interface_id) return nothing end # Initialization of boundaries after the function barrier -function init_boundaries_iter_face_inner(info, user_data) +function init_boundaries_iter_face_inner(info_pw, user_data) @unpack boundaries, boundary_id, mesh = user_data user_data.boundary_id += 1 # Extract boundary data - side = unsafe_load_side(info) + side_pw = load_pointerwrapper_side(info_pw) # Get local tree, one-based indexing - tree = unsafe_load_tree(mesh.p4est, side.treeid + 1) + tree_pw = load_pointerwrapper_tree(mesh.p4est, side_pw.treeid[] + 1) # Quadrant numbering offset of this quadrant - offset = tree.quadrants_offset + offset = tree_pw.quadrants_offset[] # Verify before accessing is.full, but this should never happen - @assert side.is_hanging == false + @assert side_pw.is_hanging[] == false - local_quad_id = side.is.full.quadid + local_quad_id = side_pw.is.full.quadid[] # Global ID of this quad quad_id = offset + local_quad_id @@ -573,52 +573,52 @@ function init_boundaries_iter_face_inner(info, user_data) boundaries.neighbor_ids[boundary_id] = quad_id + 1 # Face at which the boundary lies - face = side.face + face = side_pw.face[] # Save boundaries.node_indices dimension specific in containers_[23]d.jl init_boundary_node_indices!(boundaries, face, boundary_id) # One-based indexing - boundaries.name[boundary_id] = mesh.boundary_names[face + 1, side.treeid + 1] + boundaries.name[boundary_id] = mesh.boundary_names[face + 1, side_pw.treeid[] + 1] return nothing end # Initialization of mortars after the function barrier -function init_mortars_iter_face_inner(info, sides, user_data) +function init_mortars_iter_face_inner(info_pw, sides_pw, user_data) @unpack mortars, mortar_id, mesh = user_data user_data.mortar_id += 1 # Get Tuple of local trees, one-based indexing - trees = (unsafe_load_tree(mesh.p4est, sides[1].treeid + 1), - unsafe_load_tree(mesh.p4est, sides[2].treeid + 1)) + trees_pw = (load_pointerwrapper_tree(mesh.p4est, sides_pw[1].treeid[] + 1), + load_pointerwrapper_tree(mesh.p4est, sides_pw[2].treeid[] + 1)) # Quadrant numbering offsets of the quadrants at this interface - offsets = SVector(trees[1].quadrants_offset, - trees[2].quadrants_offset) + offsets = SVector(trees_pw[1].quadrants_offset[], + trees_pw[2].quadrants_offset[]) - if sides[1].is_hanging == true + if sides_pw[1].is_hanging[] == true # Left is small, right is large - faces = (sides[1].face, sides[2].face) + faces = (sides_pw[1].face[], sides_pw[2].face[]) - local_small_quad_ids = sides[1].is.hanging.quadid + local_small_quad_ids = sides_pw[1].is.hanging.quadid[] # Global IDs of the two small quads small_quad_ids = offsets[1] .+ local_small_quad_ids # Just be sure before accessing is.full - @assert sides[2].is_hanging == false - large_quad_id = offsets[2] + sides[2].is.full.quadid - else # sides[2].is_hanging == true + @assert sides_pw[2].is_hanging[] == false + large_quad_id = offsets[2] + sides_pw[2].is.full.quadid[] + else # sides_pw[2].is_hanging[] == true # Right is small, left is large. # init_mortar_node_indices! below expects side 1 to contain the small elements. - faces = (sides[2].face, sides[1].face) + faces = (sides_pw[2].face[], sides_pw[1].face[]) - local_small_quad_ids = sides[2].is.hanging.quadid + local_small_quad_ids = sides_pw[2].is.hanging.quadid[] # Global IDs of the two small quads small_quad_ids = offsets[2] .+ local_small_quad_ids # Just be sure before accessing is.full - @assert sides[1].is_hanging == false - large_quad_id = offsets[1] + sides[1].is.full.quadid + @assert sides_pw[1].is_hanging[] == false + large_quad_id = offsets[1] + sides_pw[1].is.full.quadid[] end # Write data to mortar container, 1 and 2 are the small elements @@ -627,7 +627,7 @@ function init_mortars_iter_face_inner(info, sides, user_data) # Last entry is the large element mortars.neighbor_ids[end, mortar_id] = large_quad_id + 1 - init_mortar_node_indices!(mortars, faces, unsafe_load(info).orientation, mortar_id) + init_mortar_node_indices!(mortars, faces, info_pw.orientation[], mortar_id) return nothing end @@ -638,34 +638,36 @@ end # - boundaries # and collect the numbers in `user_data` in this order. function count_surfaces_iter_face(info, user_data) - elem_count = unsafe_load(info).sides.elem_count + info_pw = PointerWrapper(info) + elem_count = info_pw.sides.elem_count[] if elem_count == 2 # Two neighboring elements => Interface or mortar # Extract surface data - sides = (unsafe_load_side(info, 1), unsafe_load_side(info, 2)) + sides_pw = (load_pointerwrapper_side(info_pw, 1), + load_pointerwrapper_side(info_pw, 2)) - if sides[1].is_hanging == false && sides[2].is_hanging == false + if sides_pw[1].is_hanging[] == false && sides_pw[2].is_hanging[] == false # No hanging nodes => normal interface # Unpack user_data = [interface_count] and increment interface_count - ptr = Ptr{Int}(user_data) - id = unsafe_load(ptr, 1) - unsafe_store!(ptr, id + 1, 1) + pw = PointerWrapper(Int, user_data) + id = pw[1] + pw[1] = id + 1 else # Hanging nodes => mortar # Unpack user_data = [mortar_count] and increment mortar_count - ptr = Ptr{Int}(user_data) - id = unsafe_load(ptr, 2) - unsafe_store!(ptr, id + 1, 2) + pw = PointerWrapper(Int, user_data) + id = pw[2] + pw[2] = id + 1 end elseif elem_count == 1 # One neighboring elements => boundary # Unpack user_data = [boundary_count] and increment boundary_count - ptr = Ptr{Int}(user_data) - id = unsafe_load(ptr, 3) - unsafe_store!(ptr, id + 1, 3) + pw = PointerWrapper(Int, user_data) + id = pw[3] + pw[3] = id + 1 end return nothing diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index 4f7d903897..11747f1f17 100644 --- a/src/solvers/dgsem_p4est/containers_2d.jl +++ b/src/solvers/dgsem_p4est/containers_2d.jl @@ -52,7 +52,7 @@ function calc_node_coordinates!(node_coordinates, p4est_root_len = 1 << P4EST_MAXLEVEL p4est_quadrant_len(l) = 1 << (P4EST_MAXLEVEL - l) - trees = unsafe_wrap_sc(p4est_tree_t, unsafe_load(mesh.p4est).trees) + trees = unsafe_wrap_sc(p4est_tree_t, mesh.p4est.trees) for tree in eachindex(trees) offset = trees[tree].quadrants_offset diff --git a/src/solvers/dgsem_p4est/containers_3d.jl b/src/solvers/dgsem_p4est/containers_3d.jl index 6cdc2cf961..e9994fe456 100644 --- a/src/solvers/dgsem_p4est/containers_3d.jl +++ b/src/solvers/dgsem_p4est/containers_3d.jl @@ -43,7 +43,7 @@ function calc_node_coordinates!(node_coordinates, p4est_root_len = 1 << P4EST_MAXLEVEL p4est_quadrant_len(l) = 1 << (P4EST_MAXLEVEL - l) - trees = unsafe_wrap_sc(p8est_tree_t, unsafe_load(mesh.p4est).trees) + trees = unsafe_wrap_sc(p8est_tree_t, mesh.p4est.trees) for tree in eachindex(trees) offset = trees[tree].quadrants_offset diff --git a/src/solvers/dgsem_p4est/containers_parallel.jl b/src/solvers/dgsem_p4est/containers_parallel.jl index 42d6ea44c5..e7ee1f8147 100644 --- a/src/solvers/dgsem_p4est/containers_parallel.jl +++ b/src/solvers/dgsem_p4est/containers_parallel.jl @@ -311,21 +311,24 @@ function init_surfaces_iter_face_inner(info, # surfaces at once or any subset of them, some of the unpacked values above may be `nothing` if # they're not supposed to be initialized during this call. That is why we need additional # `!== nothing` checks below before initializing individual faces. - if unsafe_load(info).sides.elem_count == 2 + info_pw = PointerWrapper(info) + if info_pw.sides.elem_count[] == 2 # Two neighboring elements => Interface or mortar # Extract surface data - sides = (unsafe_load_side(info, 1), unsafe_load_side(info, 2)) + sides_pw = (load_pointerwrapper_side(info_pw, 1), + load_pointerwrapper_side(info_pw, 2)) - if sides[1].is_hanging == false && sides[2].is_hanging == false + if sides_pw[1].is_hanging[] == false && sides_pw[2].is_hanging[] == false # No hanging nodes => normal interface or MPI interface - if sides[1].is.full.is_ghost == true || sides[2].is.full.is_ghost == true # remote side => MPI interface + if sides_pw[1].is.full.is_ghost[] == true || + sides_pw[2].is.full.is_ghost[] == true # remote side => MPI interface if mpi_interfaces !== nothing - init_mpi_interfaces_iter_face_inner(info, sides, user_data) + init_mpi_interfaces_iter_face_inner(info_pw, sides_pw, user_data) end else if interfaces !== nothing - init_interfaces_iter_face_inner(info, sides, user_data) + init_interfaces_iter_face_inner(info_pw, sides_pw, user_data) end end else @@ -333,18 +336,18 @@ function init_surfaces_iter_face_inner(info, # First, we check which side is hanging, i.e., on which side we have the refined cells. # Then we check if any of the refined cells or the coarse cell are "ghost" cells, i.e., they # belong to another rank. That way we can determine if this is a regular mortar or MPI mortar - if sides[1].is_hanging == true - @assert sides[2].is_hanging == false - if any(sides[1].is.hanging.is_ghost .== true) || - sides[2].is.full.is_ghost == true + if sides_pw[1].is_hanging[] == true + @assert sides_pw[2].is_hanging[] == false + if any(sides_pw[1].is.hanging.is_ghost[] .== true) || + sides_pw[2].is.full.is_ghost[] == true face_has_ghost_side = true else face_has_ghost_side = false end - else # sides[2].is_hanging == true - @assert sides[1].is_hanging == false - if sides[1].is.full.is_ghost == true || - any(sides[2].is.hanging.is_ghost .== true) + else # sides_pw[2].is_hanging[] == true + @assert sides_pw[1].is_hanging[] == false + if sides_pw[1].is.full.is_ghost[] == true || + any(sides_pw[2].is.hanging.is_ghost[] .== true) face_has_ghost_side = true else face_has_ghost_side = false @@ -352,15 +355,15 @@ function init_surfaces_iter_face_inner(info, end # Initialize mortar or MPI mortar if face_has_ghost_side && mpi_mortars !== nothing - init_mpi_mortars_iter_face_inner(info, sides, user_data) + init_mpi_mortars_iter_face_inner(info_pw, sides_pw, user_data) elseif !face_has_ghost_side && mortars !== nothing - init_mortars_iter_face_inner(info, sides, user_data) + init_mortars_iter_face_inner(info_pw, sides_pw, user_data) end end - elseif unsafe_load(info).sides.elem_count == 1 + elseif info_pw.sides.elem_count[] == 1 # One neighboring elements => boundary if boundaries !== nothing - init_boundaries_iter_face_inner(info, user_data) + init_boundaries_iter_face_inner(info_pw, user_data) end end @@ -381,23 +384,23 @@ function init_surfaces!(interfaces, mortars, boundaries, mpi_interfaces, mpi_mor end # Initialization of MPI interfaces after the function barrier -function init_mpi_interfaces_iter_face_inner(info, sides, user_data) +function init_mpi_interfaces_iter_face_inner(info_pw, sides_pw, user_data) @unpack mpi_interfaces, mpi_interface_id, mesh = user_data user_data.mpi_interface_id += 1 - if sides[1].is.full.is_ghost == true + if sides_pw[1].is.full.is_ghost[] == true local_side = 2 - elseif sides[2].is.full.is_ghost == true + elseif sides_pw[2].is.full.is_ghost[] == true local_side = 1 else error("should not happen") end # Get local tree, one-based indexing - tree = unsafe_load_tree(mesh.p4est, sides[local_side].treeid + 1) + tree_pw = load_pointerwrapper_tree(mesh.p4est, sides_pw[local_side].treeid[] + 1) # Quadrant numbering offset of the local quadrant at this interface - offset = tree.quadrants_offset - tree_quad_id = sides[local_side].is.full.quadid # quadid in the local tree + offset = tree_pw.quadrants_offset[] + tree_quad_id = sides_pw[local_side].is.full.quadid[] # quadid in the local tree # ID of the local neighboring quad, cumulative over local trees local_quad_id = offset + tree_quad_id @@ -406,52 +409,52 @@ function init_mpi_interfaces_iter_face_inner(info, sides, user_data) mpi_interfaces.local_sides[mpi_interface_id] = local_side # Face at which the interface lies - faces = (sides[1].face, sides[2].face) + faces = (sides_pw[1].face[], sides_pw[2].face[]) # Save mpi_interfaces.node_indices dimension specific in containers_[23]d_parallel.jl init_mpi_interface_node_indices!(mpi_interfaces, faces, local_side, - unsafe_load(info).orientation, + info_pw.orientation[], mpi_interface_id) return nothing end # Initialization of MPI mortars after the function barrier -function init_mpi_mortars_iter_face_inner(info, sides, user_data) +function init_mpi_mortars_iter_face_inner(info_pw, sides_pw, user_data) @unpack mpi_mortars, mpi_mortar_id, mesh = user_data user_data.mpi_mortar_id += 1 # Get Tuple of adjacent trees, one-based indexing - trees = (unsafe_load_tree(mesh.p4est, sides[1].treeid + 1), - unsafe_load_tree(mesh.p4est, sides[2].treeid + 1)) + trees_pw = (load_pointerwrapper_tree(mesh.p4est, sides_pw[1].treeid[] + 1), + load_pointerwrapper_tree(mesh.p4est, sides_pw[2].treeid[] + 1)) # Quadrant numbering offsets of the quadrants at this mortar - offsets = SVector(trees[1].quadrants_offset, - trees[2].quadrants_offset) + offsets = SVector(trees_pw[1].quadrants_offset[], + trees_pw[2].quadrants_offset[]) - if sides[1].is_hanging == true + if sides_pw[1].is_hanging[] == true hanging_side = 1 full_side = 2 - else # sides[2].is_hanging == true + else # sides_pw[2].is_hanging[] == true hanging_side = 2 full_side = 1 end # Just be sure before accessing is.full or is.hanging later - @assert sides[full_side].is_hanging == false - @assert sides[hanging_side].is_hanging == true + @assert sides_pw[full_side].is_hanging[] == false + @assert sides_pw[hanging_side].is_hanging[] == true # Find small quads that are locally available - local_small_quad_positions = findall(sides[hanging_side].is.hanging.is_ghost .== + local_small_quad_positions = findall(sides_pw[hanging_side].is.hanging.is_ghost[] .== false) # Get id of local small quadrants within their tree # Indexing CBinding.Caccessor via a Vector does not work here -> use map instead - tree_small_quad_ids = map(p -> sides[hanging_side].is.hanging.quadid[p], + tree_small_quad_ids = map(p -> sides_pw[hanging_side].is.hanging.quadid[][p], local_small_quad_positions) local_small_quad_ids = offsets[hanging_side] .+ tree_small_quad_ids # ids cumulative over local trees # Determine if large quadrant is available and if yes, determine its id - if sides[full_side].is.full.is_ghost == false - local_large_quad_id = offsets[full_side] + sides[full_side].is.full.quadid + if sides_pw[full_side].is.full.is_ghost[] == false + local_large_quad_id = offsets[full_side] + sides_pw[full_side].is.full.quadid[] else local_large_quad_id = -1 # large quad is ghost end @@ -470,9 +473,8 @@ function init_mpi_mortars_iter_face_inner(info, sides, user_data) mpi_mortars.local_neighbor_positions[mpi_mortar_id] = local_neighbor_positions # init_mortar_node_indices! expects side 1 to contain small elements - faces = (sides[hanging_side].face, sides[full_side].face) - init_mortar_node_indices!(mpi_mortars, faces, unsafe_load(info).orientation, - mpi_mortar_id) + faces = (sides_pw[hanging_side].face[], sides_pw[full_side].face[]) + init_mortar_node_indices!(mpi_mortars, faces, info_pw.orientation[], mpi_mortar_id) return nothing end @@ -485,42 +487,45 @@ end # - (MPI) mortars at subdomain boundaries # and collect the numbers in `user_data` in this order. function count_surfaces_iter_face_parallel(info, user_data) - if unsafe_load(info).sides.elem_count == 2 + info_pw = PointerWrapper(info) + if info_pw.sides.elem_count[] == 2 # Two neighboring elements => Interface or mortar # Extract surface data - sides = (unsafe_load_side(info, 1), unsafe_load_side(info, 2)) + sides_pw = (load_pointerwrapper_side(info_pw, 1), + load_pointerwrapper_side(info_pw, 2)) - if sides[1].is_hanging == false && sides[2].is_hanging == false + if sides_pw[1].is_hanging[] == false && sides_pw[2].is_hanging[] == false # No hanging nodes => normal interface or MPI interface - if sides[1].is.full.is_ghost == true || sides[2].is.full.is_ghost == true # remote side => MPI interface + if sides_pw[1].is.full.is_ghost[] == true || + sides_pw[2].is.full.is_ghost[] == true # remote side => MPI interface # Unpack user_data = [mpi_interface_count] and increment mpi_interface_count - ptr = Ptr{Int}(user_data) - id = unsafe_load(ptr, 4) - unsafe_store!(ptr, id + 1, 4) + pw = PointerWrapper(Int, user_data) + id = pw[4] + pw[4] = id + 1 else # Unpack user_data = [interface_count] and increment interface_count - ptr = Ptr{Int}(user_data) - id = unsafe_load(ptr, 1) - unsafe_store!(ptr, id + 1, 1) + pw = PointerWrapper(Int, user_data) + id = pw[1] + pw[1] = id + 1 end else # Hanging nodes => mortar or MPI mortar # First, we check which side is hanging, i.e., on which side we have the refined cells. # Then we check if any of the refined cells or the coarse cell are "ghost" cells, i.e., they # belong to another rank. That way we can determine if this is a regular mortar or MPI mortar - if sides[1].is_hanging == true - @assert sides[2].is_hanging == false - if any(sides[1].is.hanging.is_ghost .== true) || - sides[2].is.full.is_ghost == true + if sides_pw[1].is_hanging[] == true + @assert sides_pw[2].is_hanging[] == false + if any(sides_pw[1].is.hanging.is_ghost[] .== true) || + sides_pw[2].is.full.is_ghost[] == true face_has_ghost_side = true else face_has_ghost_side = false end - else # sides[2].is_hanging == true - @assert sides[1].is_hanging == false - if sides[1].is.full.is_ghost == true || - any(sides[2].is.hanging.is_ghost .== true) + else # sides_pw[2].is_hanging[] == true + @assert sides_pw[1].is_hanging[] == false + if sides_pw[1].is.full.is_ghost[] == true || + any(sides_pw[2].is.hanging.is_ghost[] .== true) face_has_ghost_side = true else face_has_ghost_side = false @@ -528,23 +533,23 @@ function count_surfaces_iter_face_parallel(info, user_data) end if face_has_ghost_side # Unpack user_data = [mpi_mortar_count] and increment mpi_mortar_count - ptr = Ptr{Int}(user_data) - id = unsafe_load(ptr, 5) - unsafe_store!(ptr, id + 1, 5) + pw = PointerWrapper(Int, user_data) + id = pw[5] + pw[5] = id + 1 else # Unpack user_data = [mortar_count] and increment mortar_count - ptr = Ptr{Int}(user_data) - id = unsafe_load(ptr, 2) - unsafe_store!(ptr, id + 1, 2) + pw = PointerWrapper(Int, user_data) + id = pw[2] + pw[2] = id + 1 end end - elseif unsafe_load(info).sides.elem_count == 1 + elseif info_pw.sides.elem_count[] == 1 # One neighboring elements => boundary # Unpack user_data = [boundary_count] and increment boundary_count - ptr = Ptr{Int}(user_data) - id = unsafe_load(ptr, 3) - unsafe_store!(ptr, id + 1, 3) + pw = PointerWrapper(Int, user_data) + id = pw[3] + pw[3] = id + 1 end return nothing diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl index e73a8cda9b..73ac47ed1e 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl @@ -365,6 +365,7 @@ function prolong2interfaces!(cache_parabolic, flux_viscous, return nothing end +# This version is used for divergence flux computations function calc_interface_flux!(surface_flux_values, mesh::P4estMesh{2}, equations_parabolic, dg::DG, cache_parabolic) @@ -405,7 +406,7 @@ function calc_interface_flux!(surface_flux_values, end for node in eachnode(dg) - # We prolong the viscous flux dotted with respect the outward normal on the + # We prolong the viscous flux dotted with respect the outward normal on the # primary element. We assume a BR-1 type of flux. viscous_flux_normal_ll, viscous_flux_normal_rr = get_surface_node_vars(cache_parabolic.interfaces.u, equations_parabolic, @@ -446,6 +447,7 @@ function prolong2boundaries!(cache_parabolic, flux_viscous, # a start value and a step size to get the correct face and orientation. element = boundaries.neighbor_ids[boundary] node_indices = boundaries.node_indices[boundary] + direction = indices2direction(node_indices) i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], index_range) j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], index_range) @@ -454,15 +456,12 @@ function prolong2boundaries!(cache_parabolic, flux_viscous, j_node = j_node_start for i in eachnode(dg) # this is the outward normal direction on the primary element - normal_direction = get_normal_direction(primary_direction, - contravariant_vectors, - i_node, j_node, primary_element) + normal_direction = get_normal_direction(direction, contravariant_vectors, + i_node, j_node, element) for v in eachvariable(equations_parabolic) - flux_viscous = SVector(flux_viscous_x[v, i_primary, j_primary, - primary_element], - flux_viscous_y[v, i_primary, j_primary, - primary_element]) + flux_viscous = SVector(flux_viscous_x[v, i_node, j_node, element], + flux_viscous_y[v, i_node, j_node, element]) boundaries.u[v, i, boundary] = dot(flux_viscous, normal_direction) end @@ -470,6 +469,124 @@ function prolong2boundaries!(cache_parabolic, flux_viscous, j_node += j_node_step end end + return nothing +end + +function calc_boundary_flux_gradients!(cache, t, + boundary_condition::Union{BoundaryConditionPeriodic, + BoundaryConditionDoNothing + }, + mesh::P4estMesh, equations, surface_integral, dg::DG) + @assert isempty(eachboundary(dg, cache)) +end + +# Function barrier for type stability +function calc_boundary_flux_gradients!(cache, t, boundary_conditions, mesh::P4estMesh, + equations, surface_integral, dg::DG) + (; boundary_condition_types, boundary_indices) = boundary_conditions + calc_boundary_flux_by_type!(cache, t, boundary_condition_types, boundary_indices, + Gradient(), mesh, equations, surface_integral, dg) return nothing end + +function calc_boundary_flux_divergence!(cache, t, boundary_conditions, mesh::P4estMesh, + equations, surface_integral, dg::DG) + (; boundary_condition_types, boundary_indices) = boundary_conditions + + calc_boundary_flux_by_type!(cache, t, boundary_condition_types, boundary_indices, + Divergence(), mesh, equations, surface_integral, dg) + return nothing +end + +# Iterate over tuples of boundary condition types and associated indices +# in a type-stable way using "lispy tuple programming". +function calc_boundary_flux_by_type!(cache, t, BCs::NTuple{N, Any}, + BC_indices::NTuple{N, Vector{Int}}, + operator_type, + mesh::P4estMesh, + equations, surface_integral, dg::DG) where {N} + # Extract the boundary condition type and index vector + boundary_condition = first(BCs) + boundary_condition_indices = first(BC_indices) + # Extract the remaining types and indices to be processed later + remaining_boundary_conditions = Base.tail(BCs) + remaining_boundary_condition_indices = Base.tail(BC_indices) + + # process the first boundary condition type + calc_boundary_flux!(cache, t, boundary_condition, boundary_condition_indices, + operator_type, mesh, equations, surface_integral, dg) + + # recursively call this method with the unprocessed boundary types + calc_boundary_flux_by_type!(cache, t, remaining_boundary_conditions, + remaining_boundary_condition_indices, + operator_type, + mesh, equations, surface_integral, dg) + + return nothing +end + +# terminate the type-stable iteration over tuples +function calc_boundary_flux_by_type!(cache, t, BCs::Tuple{}, BC_indices::Tuple{}, + operator_type, mesh::P4estMesh, equations, + surface_integral, dg::DG) + nothing +end + +function calc_boundary_flux!(cache, t, + boundary_condition_parabolic, # works with Dict types + boundary_condition_indices, + operator_type, mesh::P4estMesh{2}, + equations_parabolic::AbstractEquationsParabolic, + surface_integral, dg::DG) + (; boundaries) = cache + (; node_coordinates, surface_flux_values) = cache.elements + (; contravariant_vectors) = cache.elements + index_range = eachnode(dg) + + @threaded for local_index in eachindex(boundary_condition_indices) + # Use the local index to get the global boundary index from the pre-sorted list + boundary_index = boundary_condition_indices[local_index] + + # Get information on the adjacent element, compute the surface fluxes, + # and store them + element = boundaries.neighbor_ids[boundary_index] + node_indices = boundaries.node_indices[boundary_index] + direction_index = indices2direction(node_indices) + + i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], index_range) + j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], index_range) + + i_node = i_node_start + j_node = j_node_start + for node_index in eachnode(dg) + # Extract solution data from boundary container + u_inner = get_node_vars(boundaries.u, equations_parabolic, dg, node_index, + boundary_index) + + # Outward-pointing normal direction (not normalized) + normal_direction = get_normal_direction(direction_index, contravariant_vectors, + i_node, j_node, element) + + # TODO: revisit if we want more general boundary treatments. + # This assumes the gradient numerical flux at the boundary is the gradient variable, + # which is consistent with BR1, LDG. + flux_inner = u_inner + + # Coordinates at boundary node + x = get_node_coords(node_coordinates, equations_parabolic, dg, i_node, j_node, + element) + + flux_ = boundary_condition_parabolic(flux_inner, u_inner, normal_direction, + x, t, operator_type, equations_parabolic) + + # Copy flux to element storage in the correct orientation + for v in eachvariable(equations_parabolic) + surface_flux_values[v, node_index, direction_index, element] = flux_[v] + end + + i_node += i_node_step + j_node += j_node_step + end + end +end diff --git a/src/solvers/dgsem_p4est/dg_parallel.jl b/src/solvers/dgsem_p4est/dg_parallel.jl index ac122d048c..324bc7f3cd 100644 --- a/src/solvers/dgsem_p4est/dg_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_parallel.jl @@ -263,15 +263,13 @@ function init_mpi_cache!(mpi_cache::P4estMPICache, mesh::ParallelP4estMesh, uEltype) # Determine local and total number of elements - n_elements_global = Int(unsafe_load(mesh.p4est).global_num_quadrants) - n_elements_by_rank = vcat(Int.(unsafe_wrap(Array, - unsafe_load(mesh.p4est).global_first_quadrant, + n_elements_global = Int(mesh.p4est.global_num_quadrants[]) + n_elements_by_rank = vcat(Int.(unsafe_wrap(Array, mesh.p4est.global_first_quadrant, mpi_nranks())), n_elements_global) |> diff # diff sufficient due to 0-based quad indices n_elements_by_rank = OffsetArray(n_elements_by_rank, 0:(mpi_nranks() - 1)) # Account for 1-based indexing in Julia - first_element_global_id = Int(unsafe_load(unsafe_load(mesh.p4est).global_first_quadrant, - mpi_rank() + 1)) + 1 + first_element_global_id = Int(mesh.p4est.global_first_quadrant[mpi_rank() + 1]) + 1 @assert n_elements_global==sum(n_elements_by_rank) "error in total number of elements" # TODO reuse existing structures @@ -379,17 +377,19 @@ function init_neighbor_rank_connectivity_iter_face_inner(info, user_data) @unpack interfaces, interface_id, global_interface_ids, neighbor_ranks_interface, mortars, mortar_id, global_mortar_ids, neighbor_ranks_mortar, mesh = user_data + info_pw = PointerWrapper(info) # Get the global interface/mortar ids and neighbor rank if current face belongs to an MPI # interface/mortar - if unsafe_load(info).sides.elem_count == 2 # MPI interfaces/mortars have two neighboring elements + if info_pw.sides.elem_count[] == 2 # MPI interfaces/mortars have two neighboring elements # Extract surface data - sides = (unsafe_load_side(info, 1), unsafe_load_side(info, 2)) + sides_pw = (load_pointerwrapper_side(info_pw, 1), + load_pointerwrapper_side(info_pw, 2)) - if sides[1].is_hanging == false && sides[2].is_hanging == false # No hanging nodes for MPI interfaces - if sides[1].is.full.is_ghost == true + if sides_pw[1].is_hanging[] == false && sides_pw[2].is_hanging[] == false # No hanging nodes for MPI interfaces + if sides_pw[1].is.full.is_ghost[] == true remote_side = 1 local_side = 2 - elseif sides[2].is.full.is_ghost == true + elseif sides_pw[2].is.full.is_ghost[] == true remote_side = 2 local_side = 1 else # both sides are on this rank -> skip since it's a regular interface @@ -397,16 +397,17 @@ function init_neighbor_rank_connectivity_iter_face_inner(info, user_data) end # Sanity check, current face should belong to current MPI interface - local_tree = unsafe_load_tree(mesh.p4est, sides[local_side].treeid + 1) # one-based indexing - local_quad_id = local_tree.quadrants_offset + - sides[local_side].is.full.quadid + local_tree_pw = load_pointerwrapper_tree(mesh.p4est, + sides_pw[local_side].treeid[] + 1) # one-based indexing + local_quad_id = local_tree_pw.quadrants_offset[] + + sides_pw[local_side].is.full.quadid[] @assert interfaces.local_neighbor_ids[interface_id] == local_quad_id + 1 # one-based indexing # Get neighbor ID from ghost layer proc_offsets = unsafe_wrap(Array, - unsafe_load(unsafe_load(info).ghost_layer).proc_offsets, + info_pw.ghost_layer.proc_offsets, mpi_nranks() + 1) - ghost_id = sides[remote_side].is.full.quadid # indexes the ghost layer, 0-based + ghost_id = sides_pw[remote_side].is.full.quadid[] # indexes the ghost layer, 0-based neighbor_rank = findfirst(r -> proc_offsets[r] <= ghost_id < proc_offsets[r + 1], 1:mpi_nranks()) - 1 # MPI ranks are 0-based @@ -415,21 +416,18 @@ function init_neighbor_rank_connectivity_iter_face_inner(info, user_data) # Global interface id is the globally unique quadrant id of the quadrant on the primary # side (1) multiplied by the number of faces per quadrant plus face if local_side == 1 - offset = unsafe_load(unsafe_load(mesh.p4est).global_first_quadrant, - mpi_rank() + 1) # one-based indexing + offset = mesh.p4est.global_first_quadrant[mpi_rank() + 1] # one-based indexing primary_quad_id = offset + local_quad_id else - offset = unsafe_load(unsafe_load(mesh.p4est).global_first_quadrant, - neighbor_rank + 1) # one-based indexing - primary_quad_id = offset + - unsafe_load(sides[1].is.full.quad.p.piggy3.local_num) + offset = mesh.p4est.global_first_quadrant[neighbor_rank + 1] # one-based indexing + primary_quad_id = offset + sides_pw[1].is.full.quad.p.piggy3.local_num[] end - global_interface_id = 2 * ndims(mesh) * primary_quad_id + sides[1].face + global_interface_id = 2 * ndims(mesh) * primary_quad_id + sides_pw[1].face[] global_interface_ids[interface_id] = global_interface_id user_data.interface_id += 1 else # hanging node - if sides[1].is_hanging == true + if sides_pw[1].is_hanging[] == true hanging_side = 1 full_side = 2 else @@ -437,26 +435,26 @@ function init_neighbor_rank_connectivity_iter_face_inner(info, user_data) full_side = 1 end # Verify before accessing is.full / is.hanging - @assert sides[hanging_side].is_hanging == true && - sides[full_side].is_hanging == false + @assert sides_pw[hanging_side].is_hanging[] == true && + sides_pw[full_side].is_hanging[] == false # If all quadrants are locally available, this is a regular mortar -> skip - if sides[full_side].is.full.is_ghost == false && - all(sides[hanging_side].is.hanging.is_ghost .== false) + if sides_pw[full_side].is.full.is_ghost[] == false && + all(sides_pw[hanging_side].is.hanging.is_ghost[] .== false) return nothing end - trees = (unsafe_load_tree(mesh.p4est, sides[1].treeid + 1), - unsafe_load_tree(mesh.p4est, sides[2].treeid + 1)) + trees_pw = (load_pointerwrapper_tree(mesh.p4est, sides_pw[1].treeid[] + 1), + load_pointerwrapper_tree(mesh.p4est, sides_pw[2].treeid[] + 1)) # Find small quads that are remote and determine which rank owns them - remote_small_quad_positions = findall(sides[hanging_side].is.hanging.is_ghost .== + remote_small_quad_positions = findall(sides_pw[hanging_side].is.hanging.is_ghost[] .== true) proc_offsets = unsafe_wrap(Array, - unsafe_load(unsafe_load(info).ghost_layer).proc_offsets, + info_pw.ghost_layer.proc_offsets, mpi_nranks() + 1) # indices of small remote quads inside the ghost layer, 0-based - ghost_ids = map(pos -> sides[hanging_side].is.hanging.quadid[pos], + ghost_ids = map(pos -> sides_pw[hanging_side].is.hanging.quadid[][pos], remote_small_quad_positions) neighbor_ranks = map(ghost_ids) do ghost_id return findfirst(r -> proc_offsets[r] <= ghost_id < proc_offsets[r + 1], @@ -464,28 +462,26 @@ function init_neighbor_rank_connectivity_iter_face_inner(info, user_data) end # Determine global quad id of large element to determine global MPI mortar id # Furthermore, if large element is ghost, add its owner rank to neighbor_ranks - if sides[full_side].is.full.is_ghost == true - ghost_id = sides[full_side].is.full.quadid + if sides_pw[full_side].is.full.is_ghost[] == true + ghost_id = sides_pw[full_side].is.full.quadid[] large_quad_owner_rank = findfirst(r -> proc_offsets[r] <= ghost_id < proc_offsets[r + 1], 1:mpi_nranks()) - 1 # MPI ranks are 0-based push!(neighbor_ranks, large_quad_owner_rank) - offset = unsafe_load(unsafe_load(mesh.p4est).global_first_quadrant, - large_quad_owner_rank + 1) # one-based indexing + offset = mesh.p4est.global_first_quadrant[large_quad_owner_rank + 1] # one-based indexing large_quad_id = offset + - unsafe_load(sides[full_side].is.full.quad.p.piggy3.local_num) + sides_pw[full_side].is.full.quad.p.piggy3.local_num[] else - offset = unsafe_load(unsafe_load(mesh.p4est).global_first_quadrant, - mpi_rank() + 1) # one-based indexing - large_quad_id = offset + trees[full_side].quadrants_offset + - sides[full_side].is.full.quadid + offset = mesh.p4est.global_first_quadrant[mpi_rank() + 1] # one-based indexing + large_quad_id = offset + trees_pw[full_side].quadrants_offset[] + + sides_pw[full_side].is.full.quadid[] end neighbor_ranks_mortar[mortar_id] = neighbor_ranks # Global mortar id is the globally unique quadrant id of the large quadrant multiplied by the # number of faces per quadrant plus face global_mortar_ids[mortar_id] = 2 * ndims(mesh) * large_quad_id + - sides[full_side].face + sides_pw[full_side].face[] user_data.mortar_id += 1 end diff --git a/src/solvers/dgsem_tree/indicators.jl b/src/solvers/dgsem_tree/indicators.jl index 2eb0af8714..b8f8a796f2 100644 --- a/src/solvers/dgsem_tree/indicators.jl +++ b/src/solvers/dgsem_tree/indicators.jl @@ -159,7 +159,7 @@ and `basis` if this indicator should be used for shock capturing. - Löhner (1987) "An adaptive finite element scheme for transient problems in CFD" [doi: 10.1016/0045-7825(87)90098-3](https://doi.org/10.1016/0045-7825(87)90098-3) -- http://flash.uchicago.edu/site/flashcode/user_support/flash4_ug_4p62/node59.html#SECTION05163100000000000000 +- [https://flash.rochester.edu/site/flashcode/user_support/flash4_ug_4p62/node59.html#SECTION05163100000000000000](https://flash.rochester.edu/site/flashcode/user_support/flash4_ug_4p62/node59.html#SECTION05163100000000000000) """ struct IndicatorLöhner{RealT <: Real, Variable, Cache} <: AbstractIndicator f_wave::RealT # TODO: Taal documentation diff --git a/test/test_dgmulti_3d.jl b/test/test_dgmulti_3d.jl index 22c0a0fd3b..68fa1d1330 100644 --- a/test/test_dgmulti_3d.jl +++ b/test/test_dgmulti_3d.jl @@ -135,6 +135,12 @@ isdir(outdir) && rm(outdir, recursive=true) ) end + @trixi_testset "elixir_advection_tensor_wedge.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_tensor_wedge.jl"), + l2 = [2.30487910e-04] , + linf = [6.31795281e-04] ) + end + end # Clean up afterwards: delete Trixi.jl output directory diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index b0ac63d4ce..471b976e99 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -200,6 +200,38 @@ isdir(outdir) && rm(outdir, recursive=true) ) end + @trixi_testset "P4estMesh2D: elixir_advection_diffusion_periodic_curved.jl" begin + @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_advection_diffusion_periodic_curved.jl"), + trees_per_dimension = (1, 1), initial_refinement_level = 2, tspan=(0.0, 0.5), + l2 = [0.012380458938507371], + linf = [0.10860506906472567] + ) + end + + @trixi_testset "P4estMesh2D: elixir_advection_diffusion_nonperiodic_curved.jl" begin + @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_advection_diffusion_nonperiodic_curved.jl"), + trees_per_dimension = (1, 1), initial_refinement_level = 2, tspan=(0.0, 0.5), + l2 = [0.04933902988507035], + linf = [0.2550261714590271] + ) + end + + @trixi_testset "P4estMesh2D: elixir_navierstokes_convergence.jl" begin + @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_navierstokes_convergence.jl"), + initial_refinement_level = 1, tspan=(0.0, 0.2), + l2 = [0.0003811978985836709, 0.0005874314969169538, 0.0009142898787923481, 0.0011613918899727263], + linf = [0.0021633623982135752, 0.009484348274135372, 0.004231572066492217, 0.011661660275365193] + ) + end + + @trixi_testset "P4estMesh2D: elixir_navierstokes_lid_driven_cavity.jl" begin + @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_navierstokes_lid_driven_cavity.jl"), + initial_refinement_level = 2, tspan=(0.0, 0.5), + l2 = [0.00028716166408816073, 0.08101204560401647, 0.02099595625377768, 0.05008149754143295], + linf = [0.014804500261322406, 0.9513271652357098, 0.7223919625994717, 1.4846907331004786] + ) + end + end # Clean up afterwards: delete Trixi.jl output directory diff --git a/test/test_structured_1d.jl b/test/test_structured_1d.jl index a27d3c219e..ec8c7a138d 100644 --- a/test/test_structured_1d.jl +++ b/test/test_structured_1d.jl @@ -27,8 +27,8 @@ isdir(outdir) && rm(outdir, recursive=true) @trixi_testset "elixir_advection_shockcapturing.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_shockcapturing.jl"), - l2 = [0.08004076716881656], - linf = [0.6342577638501385], + l2 = [0.08015029105233593], + linf = [0.610709468736576], atol = 1.0e-5) end diff --git a/test/test_tree_1d.jl b/test/test_tree_1d.jl index e37e1efc3e..7737a93a15 100644 --- a/test/test_tree_1d.jl +++ b/test/test_tree_1d.jl @@ -271,7 +271,7 @@ end sol = solve(ode, Tsit5(), abstol=1.0e-6, reltol=1.0e-6, save_everystep=false, callback=callbacks); - @test analysis_callback(sol).l2 ≈ [0.00029610274971929974, 5.573684084938363e-6] + @test analysis_callback(sol).l2 ≈ [0.00029609575838969394, 5.5681704039507985e-6] end diff --git a/test/test_tree_1d_burgers.jl b/test/test_tree_1d_burgers.jl index 788c7ab419..8c4cfaa406 100644 --- a/test/test_tree_1d_burgers.jl +++ b/test/test_tree_1d_burgers.jl @@ -22,14 +22,14 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_burgers_shock.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_burgers_shock.jl"), - l2 = [0.4407585104869119], - linf = [1.000000000000001]) + l2 = [0.4422505602587537], + linf = [1.0000000000000009]) end @trixi_testset "elixir_burgers_rarefaction.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_burgers_rarefaction.jl"), - l2 = [0.40287062735307044], - linf = [1.0042992585765542]) + l2 = [0.4038224690923722], + linf = [1.0049201454652736]) end end diff --git a/test/test_tree_1d_euler.jl b/test/test_tree_1d_euler.jl index 40f2a38b0e..5fb74b80bc 100644 --- a/test/test_tree_1d_euler.jl +++ b/test/test_tree_1d_euler.jl @@ -22,8 +22,8 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_euler_density_wave.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_density_wave.jl"), - l2 = [0.0011482554820185795, 0.00011482554830363504, 5.741277417754598e-6], - linf = [0.004090978306820037, 0.00040909783134346345, 2.0454891732413216e-5]) + l2 = [0.0011482554820217855, 0.00011482554830323462, 5.741277429325267e-6], + linf = [0.004090978306812376, 0.0004090978313582294, 2.045489210189544e-5]) end @trixi_testset "elixir_euler_density_wave.jl with initial_condition_constant" begin @@ -41,14 +41,14 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_euler_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"), - l2 = [0.11915540925414216, 0.15489191247295198, 0.44543052524765375], - linf = [0.2751485868543495, 0.2712764982000735, 0.9951407418216425]) + l2 = [0.11821957357197649, 0.15330089521538678, 0.4417674632047301], + linf = [0.24280567569982958, 0.29130548795961936, 0.8847009003152442]) end @trixi_testset "elixir_euler_ec.jl with flux_kennedy_gruber" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"), - l2 = [0.07905582221868049, 0.10180958900546237, 0.29596551476711125], - linf = [0.23515297345769826, 0.2958208108392532, 0.8694224308790321], + l2 = [0.07803455838661963, 0.10032577312032283, 0.29228156303827935], + linf = [0.2549869853794955, 0.3376472164661263, 0.9650477546553962], maxiters = 10, surface_flux = flux_kennedy_gruber, volume_flux = flux_kennedy_gruber) @@ -56,8 +56,8 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_euler_ec.jl with flux_shima_etal" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"), - l2 = [0.07909267609417114, 0.1018246500951966, 0.2959649187481973], - linf = [0.23631829743146504, 0.2977756307879202, 0.8642794698697331], + l2 = [0.07800654460172655, 0.10030365573277883, 0.2921481199111959], + linf = [0.25408579350400395, 0.3388657679031271, 0.9776486386921928], maxiters = 10, surface_flux = flux_shima_etal, volume_flux = flux_shima_etal) @@ -65,8 +65,8 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_euler_ec.jl with flux_chandrashekar" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"), - l2 = [0.07905306555214126, 0.10181180378499956, 0.2959171937479504], - linf = [0.24057642004451651, 0.29691454643616433, 0.886425723870524], + l2 = [0.07801923089205756, 0.10039557434912669, 0.2922210399923278], + linf = [0.2576521982607225, 0.3409717926625057, 0.9772961936567048], maxiters = 10, surface_flux = flux_chandrashekar, volume_flux = flux_chandrashekar) @@ -74,8 +74,8 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_euler_ec.jl with flux_hll" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"), - l2 = [0.07959780803600519, 0.10342491934977621, 0.2978851659149904], - linf = [0.19228754121840885, 0.2524152253292552, 0.725604944702432], + l2 = [0.07852272782240548, 0.10209790867523805, 0.293873048809011], + linf = [0.19244768908604093, 0.2515941686151897, 0.7258000837553769], maxiters = 10, surface_flux = flux_hll, volume_flux = flux_ranocha) @@ -83,8 +83,8 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_euler_shockcapturing.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_shockcapturing.jl"), - l2 = [0.11665968950973675, 0.15105507394693413, 0.43503082674771115], - linf = [0.1867400345208743, 0.24621854448555328, 0.703826406555577]) + l2 = [0.11606096465319675, 0.15028768943458806, 0.4328230323046703], + linf = [0.18031710091067965, 0.2351582421501841, 0.6776805692092567]) end @trixi_testset "elixir_euler_sedov_blast_wave.jl" begin @@ -96,8 +96,8 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_euler_sedov_blast_wave_pure_fv.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave_pure_fv.jl"), - l2 = [1.075075094036344, 0.06766902169711514, 0.9221426570128292], - linf = [3.3941512671408542, 0.16862631133303882, 2.6572394126490315], + l2 = [1.0735456065491455, 0.07131078703089379, 0.9205739468590453], + linf = [3.4296365168219216, 0.17635583964559245, 2.6574584326179505], # Let this test run longer to cover some lines in flux_hllc coverage_override = (maxiters=10^5, tspan=(0.0, 0.1))) end @@ -129,22 +129,22 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_euler_blast_wave.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave.jl"), - l2 = [0.21651329948737183, 0.28091709900008616, 0.5580778880050432], - linf = [1.513525457073142, 1.5328754303137992, 2.0467706106669556], + l2 = [0.21934822867340323, 0.28131919126002686, 0.554361702716662], + linf = [1.5180897390290355, 1.3967085956620369, 2.0663825294019595], maxiters = 30) end @trixi_testset "elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl"), - l2 = [2.13605618e-01, 2.79953055e-01, 5.54424459e-01], - linf = [1.55151701e+00, 1.55696782e+00, 2.05525953e+00], + l2 = [0.21814833203212694, 0.2818328665444332, 0.5528379124720818], + linf = [1.5548653877320868, 1.4474018998129738, 2.071919577393772], maxiters = 30) end @trixi_testset "elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl"), - l2 = [2.18148857e-01, 2.83182959e-01, 5.59096194e-01], - linf = [1.62706876e+00, 1.61680275e+00, 2.05876517e+00], + l2 = [0.22054468879127423, 0.2828269190680846, 0.5542369885642424], + linf = [1.5623359741479623, 1.4290121654488288, 2.1040405133123072], maxiters = 30) end end diff --git a/test/test_tree_1d_eulermulti.jl b/test/test_tree_1d_eulermulti.jl index eac54a9372..e880f98e2d 100644 --- a/test/test_tree_1d_eulermulti.jl +++ b/test/test_tree_1d_eulermulti.jl @@ -11,39 +11,47 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_eulermulti_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_ec.jl"), - l2 = [1.54891912e-01, 4.45430525e-01, 1.70222013e-02, 3.40444026e-02, 6.80888053e-02], - linf = [2.71276498e-01, 9.95140742e-01, 3.93069410e-02, 7.86138820e-02, 1.57227764e-01]) + l2 = [0.15330089521538684, 0.4417674632047301, 0.016888510510282385, 0.03377702102056477, + 0.06755404204112954], + linf = [0.29130548795961864, 0.8847009003152357, 0.034686525099975274, 0.06937305019995055, + 0.1387461003999011]) end @trixi_testset "elixir_eulermulti_es.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_es.jl"), - l2 = [1.53387916e-01, 4.41585576e-01, 3.93605635e-02, 7.87211270e-02], - linf = [2.49632117e-01, 7.21088064e-01, 6.38328770e-02, 1.27665754e-01]) + l2 = [0.1522380497572071, 0.43830846465313206, 0.03907262116499431, 0.07814524232998862], + linf = [0.24939193075537294, 0.7139395740052739, 0.06324208768391237, 0.12648417536782475]) end @trixi_testset "elixir_eulermulti_convergence_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_ec.jl"), - l2 = [8.57523604e-05, 1.63878043e-04, 1.94126993e-05, 3.88253986e-05], - linf = [3.05932773e-04, 6.24480393e-04, 7.25312144e-05, 1.45062429e-04]) + l2 = [8.575236038539227e-5, 0.00016387804318585358, 1.9412699303977585e-5, 3.882539860795517e-5], + linf = [0.00030593277277124464, 0.0006244803933350696, 7.253121435135679e-5, 0.00014506242870271358]) end @trixi_testset "elixir_eulermulti_convergence_es.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_es.jl"), - l2 = [1.8983933794407234e-5, 6.207744299844731e-5, 1.5466205761868047e-6, 3.0932411523736094e-6, 6.186482304747219e-6, 1.2372964609494437e-5], - linf = [0.00012014372605895218, 0.0003313207215800418, 6.50836791016296e-6, 1.301673582032592e-5, 2.603347164065184e-5, 5.206694328130368e-5]) + l2 = [1.8983933794407234e-5, 6.207744299844731e-5, 1.5466205761868047e-6, 3.0932411523736094e-6, + 6.186482304747219e-6, 1.2372964609494437e-5], + linf = [0.00012014372605895218, 0.0003313207215800418, 6.50836791016296e-6, 1.301673582032592e-5, + 2.603347164065184e-5, 5.206694328130368e-5]) end @trixi_testset "elixir_eulermulti_convergence_es.jl with flux_chandrashekar" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_es.jl"), - l2 = [1.88845048e-05, 5.49106005e-05, 9.42673716e-07, 1.88534743e-06, 3.77069486e-06, 7.54138973e-06], - linf = [1.16223512e-04, 3.07922197e-04, 3.21774233e-06, 6.43548465e-06, 1.28709693e-05, 2.57419386e-05], + l2 = [1.888450477353845e-5, 5.4910600482795386e-5, 9.426737161533622e-7, 1.8853474323067245e-6, + 3.770694864613449e-6, 7.541389729226898e-6], + linf = [0.00011622351152063004, 0.0003079221967086099, 3.2177423254231563e-6, 6.435484650846313e-6, + 1.2870969301692625e-5, 2.574193860338525e-5], volume_flux = flux_chandrashekar) end @trixi_testset "elixir_eulermulti_two_interacting_blast_waves.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_two_interacting_blast_waves.jl"), - l2 = [1.28886761e+00, 8.27133526e+01, 3.50680272e-03, 1.36987844e-02, 1.91795185e-02], - linf = [2.96413045e+01, 1.32258448e+03, 9.19191937e-02, 3.10929710e-01, 4.41798976e-01], + l2 = [1.288867611915533, 82.71335258388848, 0.00350680272313187, 0.013698784353152794, + 0.019179518517518084], + linf = [29.6413044707026, 1322.5844802186496, 0.09191919374782143, 0.31092970966717925, + 0.4417989757182038], tspan = (0.0, 0.0001)) end diff --git a/test/test_tree_1d_fdsbp.jl b/test/test_tree_1d_fdsbp.jl index a966b3836f..118385c34b 100644 --- a/test/test_tree_1d_fdsbp.jl +++ b/test/test_tree_1d_fdsbp.jl @@ -7,6 +7,24 @@ include("test_trixi.jl") EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_fdsbp") +@testset "Linear scalar advection" begin + @trixi_testset "elixir_advection_upwind.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_upwind.jl"), + l2 = [1.7735637157305526e-6], + linf = [1.0418854521951328e-5], + tspan = (0.0, 0.5)) + + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end +end + @testset "Inviscid Burgers" begin @trixi_testset "elixir_burgers_basic.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_burgers_basic.jl"), diff --git a/test/test_tree_1d_mhd.jl b/test/test_tree_1d_mhd.jl index 938959831c..e3a0cda325 100644 --- a/test/test_tree_1d_mhd.jl +++ b/test/test_tree_1d_mhd.jl @@ -32,14 +32,14 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_mhd_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec.jl"), - l2 = [5.86009540e-02, 8.16048158e-02, 5.46791194e-02, 5.46791194e-02, 1.54509265e-01, 4.13046273e-17, 5.47637521e-02, 5.47637521e-02], - linf = [1.10014999e-01, 1.81982581e-01, 9.13611439e-02, 9.13611439e-02, 4.23831370e-01, 1.11022302e-16, 9.93731761e-02, 9.93731761e-02]) + l2 = [0.05815183849746399, 0.08166807325621023, 0.054659228513541165, 0.054659228513541165, 0.15578125987042743, 4.130462730494e-17, 0.05465258887150046, 0.05465258887150046], + linf = [0.12165312668363826, 0.1901920742264952, 0.10059813883022554, 0.10059813883022554, 0.44079257431070706, 1.1102230246251565e-16, 0.10528911365809579, 0.10528911365809579]) end @trixi_testset "elixir_mhd_briowu_shock_tube.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_briowu_shock_tube.jl"), - l2 = [0.17764301067932906, 0.19693621875378622, 0.3635136528288642, 0.0, 0.3757321708837591, 8.593007507325741e-16, 0.36473438378159656, 0.0], - linf = [0.5601530250396535, 0.43867368105486537, 1.0960903616351099, 0.0, 1.0551794137886303, 4.107825191113079e-15, 1.5374410890043144, 0.0], + l2 = [0.17477712356961989, 0.19489623595086944, 0.3596546157640463, 0.0, 0.3723215736814466, 1.2060075775846403e-15, 0.36276754492568164, 0.0], + linf = [0.5797109945880677, 0.4372991899547103, 1.0906536287185835, 0.0, 1.0526758874956808, 5.995204332975845e-15, 1.5122922036932964, 0.0], coverage_override = (maxiters=6,)) end @@ -51,8 +51,8 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_mhd_ryujones_shock_tube.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ryujones_shock_tube.jl"), - l2 = [2.34809441e-01, 3.92255943e-01, 8.23575546e-02, 1.75599624e-01, 9.61613519e-01, 6.60825891e-17, 2.15346454e-01, 1.07006529e-01], - linf = [6.40732148e-01, 9.44889516e-01, 3.54932707e-01, 8.54060243e-01, 2.07757711e+00, 1.11022302e-16, 4.92584725e-01, 2.49526561e-01], + l2 = [0.23469781891518154, 0.3916675299696121, 0.08245195301016353, 0.1745346945706147, 0.9606363432904367, 6.608258910237605e-17, 0.21542929107153735, 0.10705457908737925], + linf = [0.6447951791685409, 0.9461857095377463, 0.35074627554617605, 0.8515177411529542, 2.0770652030507053, 1.1102230246251565e-16, 0.49670855513788204, 0.24830199967863564], tspan = (0.0, 0.1)) end diff --git a/test/test_tree_1d_mhdmulti.jl b/test/test_tree_1d_mhdmulti.jl index 2985e6d566..5214ed26d3 100644 --- a/test/test_tree_1d_mhdmulti.jl +++ b/test/test_tree_1d_mhdmulti.jl @@ -11,33 +11,33 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_mhdmulti_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_ec.jl"), - l2 = [0.08160481582829862, 0.05467911944326103, 0.05467911944326103, 0.15450926504459692, - 4.130462730494e-17, 0.054763752050210085, 0.054763752050210085, 0.008371564857135208, - 0.016743129714270416, 0.03348625942854083], - linf = [0.18198258075330706, 0.09136114386311774, 0.09136114386311774, 0.423831369951313, - 1.1102230246251565e-16, 0.09937317613143604, 0.09937317613143604, 0.0157164284712992, - 0.0314328569425984, 0.0628657138851968]) + l2 = [0.08166807325620999, 0.054659228513541616, 0.054659228513541616, 0.15578125987042812, + 4.130462730494e-17, 0.054652588871500665, 0.054652588871500665, 0.008307405499637766, + 0.01661481099927553, 0.03322962199855106], + linf = [0.19019207422649645, 0.10059813883022888, 0.10059813883022888, 0.4407925743107146, + 1.1102230246251565e-16, 0.10528911365809623, 0.10528911365809623, 0.01737901809766182, + 0.03475803619532364, 0.06951607239064728]) end @trixi_testset "elixir_mhdmulti_ec.jl with flux_derigs_etal" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_ec.jl"), - l2 = [0.08153372259925547, 0.05464109003345891, 0.05464109003345891, 0.1540576724164453, - 4.130462730494e-17, 0.054734930802131036, 0.054734930802131036, 0.008391254781284321, - 0.016782509562568642, 0.033565019125137284], - linf = [0.17492544007323832, 0.09029632168248182, 0.09029632168248182, 0.40798609353896564, - 1.1102230246251565e-16, 0.09872923637833075, 0.09872923637833075, 0.01609818847160674, - 0.03219637694321348, 0.06439275388642696], + l2 = [0.08151404166186461, 0.054640238302693274, 0.054640238302693274, 0.15536125426328573, + 4.130462730494e-17, 0.054665489963920275, 0.054665489963920275, 0.008308349501359825, + 0.01661669900271965, 0.0332333980054393], + linf = [0.1824424257860952, 0.09734687137001484, 0.09734687137001484, 0.4243089502087325, + 1.1102230246251565e-16, 0.09558639591092555, 0.09558639591092555, 0.017364773041550624, + 0.03472954608310125, 0.0694590921662025], volume_flux = flux_derigs_etal) end @trixi_testset "elixir_mhdmulti_es.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_es.jl"), - l2 = [0.07968782477167513, 0.05398115008116676, 0.05398115008116676, 0.15015281822439228, - 4.130462730494e-17, 0.053629890024921495, 0.053629890024921495, 0.008279068245579706, - 0.016558136491159413, 0.033116272982318826], - linf = [0.14118014632124837, 0.07820697032983395, 0.07820697032983395, 0.3390558674728652, - 1.1102230246251565e-16, 0.06998787893467828, 0.06998787893467828, 0.014943825414763745, - 0.02988765082952749, 0.05977530165905498]) + l2 = [0.07994082660130175, 0.053940174914031976, 0.053940174914031976, 0.15165513559250643, + 4.130462730494e-17, 0.05363207135290325, 0.05363207135290325, 0.008258265884659555, + 0.01651653176931911, 0.03303306353863822], + linf = [0.14101014428198477, 0.07762441749521025, 0.07762441749521025, 0.3381334453289866, + 1.1102230246251565e-16, 0.07003646400675223, 0.07003646400675223, 0.014962483760600165, + 0.02992496752120033, 0.05984993504240066]) end @trixi_testset "elixir_mhdmulti_convergence.jl" begin @@ -52,12 +52,12 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_mhdmulti_briowu_shock_tube.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_briowu_shock_tube.jl"), - l2 = [0.1946577804333822, 0.3591196215528672, 0.0, 0.36875476066849383, - 4.7644020131827105e-16, 0.36668249926193885, 0.0, 0.05775369214541893, - 0.11550738429083786], - linf = [0.4345551123140612, 1.0874941615375844, 0.0, 1.0493729052116585, - 3.219646771412954e-15, 1.5160434573973656, 0.0, 0.18616213071936066, - 0.3723242614387213], + l2 = [0.1877830835572639, 0.3455841730726793, 0.0, 0.35413123388836687, + 8.745556626531982e-16, 0.3629920109231055, 0.0, 0.05329005553971236, + 0.10658011107942472], + linf = [0.4288187627971754, 1.0386547815614993, 0.0, 0.9541678878162702, + 5.773159728050814e-15, 1.4595119339458051, 0.0, 0.18201910908829552, + 0.36403821817659104], coverage_override = (maxiters=6,)) end diff --git a/test/test_tree_1d_shallowwater.jl b/test/test_tree_1d_shallowwater.jl index f8901a3dcb..1c3bac1fab 100644 --- a/test/test_tree_1d_shallowwater.jl +++ b/test/test_tree_1d_shallowwater.jl @@ -10,22 +10,30 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @testset "Shallow Water" begin @trixi_testset "elixir_shallowwater_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec.jl"), - l2 = [0.8122354510732459, 1.01586214815876, 0.43404255061704217], - linf = [1.4883285368551107, 3.8717508164234276, 1.7711213427919539], + l2 = [0.244729018751225, 0.8583565222389505, 0.07330427577586297], + linf = [2.1635021283528504, 3.8717508164234453, 1.7711213427919539], + tspan = (0.0, 0.25)) + end + + @trixi_testset "elixir_shallowwater_ec.jl with initial_condition_weak_blast_wave" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec.jl"), + l2 = [0.39464782107209717, 2.03880864210846, 4.1623084150546725e-10], + linf = [0.778905801278281, 3.2409883402608273, 7.419800190922032e-10], + initial_condition=initial_condition_weak_blast_wave, tspan = (0.0, 0.25)) end @trixi_testset "elixir_shallowwater_well_balanced.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), - l2 = [1.2427984842961743, 1.0332499675061871e-14, 1.2427984842961741], - linf = [1.619041478244762, 1.266865149831811e-14, 1.6190414782447629], + l2 = [0.10416666834254829, 1.4352935256803184e-14, 0.10416666834254838], + linf = [1.9999999999999996, 3.248036646353028e-14, 2.0], tspan = (0.0, 0.25)) end @trixi_testset "elixir_shallowwater_well_balanced.jl with FluxHydrostaticReconstruction" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), - l2 = [1.2427984842961743, 1.2663646513352053e-14, 1.2427984842961741], - linf = [1.619041478244762, 2.4566658711604395e-14, 1.6190414782447629], + l2 = [0.10416666834254835, 1.1891029971551825e-14, 0.10416666834254838], + linf = [2.0000000000000018, 2.4019608337954543e-14, 2.0], surface_flux=(FluxHydrostaticReconstruction(flux_lax_friedrichs, hydrostatic_reconstruction_audusse_etal), flux_nonconservative_audusse_etal), tspan = (0.0, 0.25)) end @@ -59,14 +67,14 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") tspan = (0.0, 0.025)) end - @trixi_testset "elixir_shallowwater_well_balanced_nonperiodic.jl with dirichlet boundary" begin + @trixi_testset "elixir_shallowwater_well_balanced_nonperiodic.jl with Dirichlet boundary" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced_nonperiodic.jl"), l2 = [1.725964362045055e-8, 5.0427180314307505e-16, 1.7259643530442137e-8], linf = [3.844551077492042e-8, 3.469453422316143e-15, 3.844551077492042e-8], tspan = (0.0, 0.25)) end - @trixi_testset "elixir_shallowwater_well_nonperiodic.jl with wall boundary" begin + @trixi_testset "elixir_shallowwater_well_balanced_nonperiodic.jl with wall boundary" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced_nonperiodic.jl"), l2 = [1.7259643614361866e-8, 3.5519018243195145e-16, 1.7259643530442137e-8], linf = [3.844551010878661e-8, 9.846474508971374e-16, 3.844551077492042e-8], @@ -76,8 +84,8 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_shallowwater_shock_capturing.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_shock_capturing.jl"), - l2 = [0.2884024818919076, 0.5252262013521178, 0.2890348477852955], - linf = [0.7565706154863958, 2.076621603471687, 0.8646939843534258], + l2 = [0.07424140641160326, 0.2148642632748155, 0.0372579849000542], + linf = [1.1209754279344226, 1.3230788645853582, 0.8646939843534251], tspan = (0.0, 0.05)) end end diff --git a/test/test_tree_1d_shallowwater_twolayer.jl b/test/test_tree_1d_shallowwater_twolayer.jl index 6c0ad2941c..0d8a83806f 100644 --- a/test/test_tree_1d_shallowwater_twolayer.jl +++ b/test/test_tree_1d_shallowwater_twolayer.jl @@ -38,9 +38,9 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_shallowwater_twolayer_dam_break.jl with flux_lax_friedrichs" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_twolayer_dam_break.jl"), - l2 = [0.35490827242437256, 1.6715402155795918, 0.6960264969949427, - 0.9351481433409805, 0.7938172946965545], - linf = [0.6417127471419837, 1.9742107034120873, 1.135774587483082, 1.236125279347084, 1.1], + l2 = [0.10010269243463918, 0.5668733957648654, 0.08759617327649398, + 0.4538443183566172, 0.013638618139749523], + linf = [0.5854202777756559, 2.1278930820498934, 0.5193686074348809, 1.8071213168086229, 0.5], surface_flux = (flux_lax_friedrichs, flux_nonconservative_fjordholm_etal), tspan = (0.0, 0.25)) end diff --git a/utils/pre-commit b/utils/pre-commit new file mode 100755 index 0000000000..73ad061bae --- /dev/null +++ b/utils/pre-commit @@ -0,0 +1,40 @@ +#!/bin/bash + +# Copy this file into .git/hooks/pre-commit to execute before each commit. +# It checks and corrects the format for each file. +# If incorrect formatting is found you can add the correction via git add -p + +echo "Checking format before committing" + +if git ref-parse --verify HEAD >/dev/null 2>&1 +then + against=HEAD +else + # Initial commit: diff against an empty tree object + against=280fc57fade28e35046c3e884e587ffef05d3867 +fi + +# Redirect output to stderr. +exec 1>&2 + +# Create a list of files to format. +files=() + +for file in `git diff --cached --name-only` +do + # only indent existing files, this is necessary since if we rename or delete + # a file it is added to the committed files and we thus would try to indent a + # nonexisting file. + if [ ! -e $file ] + then + continue + fi + # We only indent .jl files + FILE_ENDING="${file##*.}" + if [ $FILE_ENDING = "jl" ] + then + files+=($file) + fi +done + +julia utils/trixi-format-file.jl "${files[@]}" diff --git a/utils/trixi-format-file.jl b/utils/trixi-format-file.jl new file mode 100755 index 0000000000..c4d8e7c903 --- /dev/null +++ b/utils/trixi-format-file.jl @@ -0,0 +1,42 @@ +#!/usr/bin/env julia + +using Pkg +Pkg.activate(; temp = true, io = devnull) +Pkg.add("JuliaFormatter"; preserve = PRESERVE_ALL, io = devnull) + +using JuliaFormatter: format_file + +function main() + # Show help + if "-h" in ARGS || "--help" in ARGS + println("usage: trixi-format.jl PATH [PATH...]") + println() + println("positional arguments:") + println() + println(" PATH One or more paths (directories or files) to format. Default: '.'") + return nothing + end + + file_list = ARGS + if isempty(ARGS) + exit(0) + end + non_formatted_files = Vector{String}() + for file in file_list + println("Checking file " * file) + if !format_file(file) + push!(non_formatted_files, file) + end + end + if isempty(non_formatted_files) + exit(0) + else + @error "Some files have not been formatted! Formatting has been applied, run 'git add -p' to update changes." + for file in non_formatted_files + println(file) + end + exit(1) + end +end + +main()