diff --git a/NEWS.md b/NEWS.md index 5dd911391a6..265979c3508 100644 --- a/NEWS.md +++ b/NEWS.md @@ -7,6 +7,7 @@ for human readability. ## Changes when updating to v0.6 from v0.5.x #### Added +- AMR for hyperbolic-parabolic equations on 2D `P4estMesh` #### Changed diff --git a/Project.toml b/Project.toml index 17c02a2adac..4b7b69af93b 100644 --- a/Project.toml +++ b/Project.toml @@ -66,7 +66,7 @@ Makie = "0.19" MuladdMacro = "0.2.2" Octavian = "0.3.5" OffsetArrays = "1.3" -P4est = "0.4" +P4est = "0.4.9" Polyester = "0.7.5" PrecompileTools = "1.1" Printf = "1" @@ -84,7 +84,7 @@ StaticArrays = "1" StrideArrays = "0.1.18" StructArrays = "0.6" SummationByPartsOperators = "0.5.41" -T8code = "0.4.1" +T8code = "0.4.3" TimerOutputs = "0.5" Triangulate = "2.0" TriplotBase = "0.1" diff --git a/docs/src/parallelization.md b/docs/src/parallelization.md index fa6fc1a5d32..3cce7c381b2 100644 --- a/docs/src/parallelization.md +++ b/docs/src/parallelization.md @@ -53,30 +53,43 @@ a system-provided MPI installation with Trixi.jl can be found in the following s ### [Using a system-provided MPI installation](@id parallel_system_MPI) -When using Trixi.jl with a system-provided MPI backend the underlying -[`p4est`](https://github.com/cburstedde/p4est) and [`t8code`](https://github.com/DLR-AMR/t8code) -libraries need to be compiled with the same MPI installation. Therefore, you also need to -use system-provided `p4est` and `t8code` installations (for notes on how to install `p4est` -and `t8code` see e.g. [here](https://github.com/cburstedde/p4est/blob/master/README) and -[here](https://github.com/DLR-AMR/t8code/wiki/Installation), use the configure option -`--enable-mpi`). Note that `t8code` already comes with a `p4est` installation, so it suffices -to install `t8code`. In addition, [P4est.jl](https://github.com/trixi-framework/P4est.jl) and -[T8code.jl](https://github.com/DLR-AMR/T8code.jl) need to be configured to use the custom -installations. Follow the steps described -[here](https://github.com/DLR-AMR/T8code.jl/blob/main/README.md#installation) and -[here](https://github.com/trixi-framework/P4est.jl/blob/main/README.md#installation) for the -configuration. The paths that point to `libp4est.so` (and potentially to `libsc.so`) need to be -the same for P4est.jl and T8code.jl. This could e.g. be `libp4est.so` that usually can be found -in `lib/` or `local/lib/` in the installation directory of `t8code`. -In total, in your active Julia project you should have a LocalPreferences.toml file with sections -`[MPIPreferences]`, `[T8code]` and `[P4est]` as well as an entry `MPIPreferences` in your -Project.toml to use a custom MPI installation. A `LocalPreferences.toml` file +When using Trixi.jl with a system-provided MPI backend, the underlying +[`p4est`](https://github.com/cburstedde/p4est), [`t8code`](https://github.com/DLR-AMR/t8code) +and [`HDF5`](https://github.com/HDFGroup/hdf5) libraries need to be compiled with the same MPI +installation. If you want to use `p4est` (via the `P4estMesh`) or `t8code` (via the `T8codeMesh`) +from Trixi.jl, you also need to use system-provided `p4est` or `t8code` installations +(for notes on how to install `p4est` and `t8code` see, e.g., [here](https://github.com/cburstedde/p4est/blob/master/README) +and [here](https://github.com/DLR-AMR/t8code/wiki/Installation), use the configure option +`--enable-mpi`). Otherwise, there will be warnings that no preference is set for P4est.jl and +T8code.jl that can be ignored if you do not use these libraries from Trixi.jl. Note that +`t8code` already comes with a `p4est` installation, so it suffices to install `t8code`. +In order to use system-provided `p4est` and `t8code` installations, [P4est.jl](https://github.com/trixi-framework/P4est.jl) +and [T8code.jl](https://github.com/DLR-AMR/T8code.jl) need to be configured to use the custom +installations. Follow the steps described [here](https://github.com/DLR-AMR/T8code.jl/blob/main/README.md#installation) and +[here](https://github.com/trixi-framework/P4est.jl/blob/main/README.md#installation). +for the configuration. The paths that point to `libp4est.so` (and potentially to `libsc.so`) need to be +the same for P4est.jl and T8code.jl. This could, e.g., be `libp4est.so` that usually can be found +in `lib/` or `local/lib/` in the installation directory of `t8code`. Note that the `T8codeMesh`, however, +does not support MPI yet. +The preferences for [HDF5.jl](https://github.com/JuliaIO/HDF5.jl) always need to be set, even if you +do not want to use `HDF5` from Trixi.jl, see also https://github.com/JuliaIO/HDF5.jl/issues/1079. +To set the preferences for HDF5.jl, follow the instructions described +[here](https://trixi-framework.github.io/Trixi.jl/stable/parallelization/#Using-parallel-input-and-output). + +In total, in your active Julia project you should have a `LocalPreferences.toml` file with sections +`[MPIPreferences]`, `[T8code]` (only needed if `T8codeMesh` is used), `[P4est]` (only needed if +`P4estMesh` is used), and `[HDF5]` as well as an entry `MPIPreferences` in your +`Project.toml` to use a custom MPI installation. A `LocalPreferences.toml` file created as described above might look something like the following: ```toml [HDF5] libhdf5 = "/usr/lib/x86_64-linux-gnu/hdf5/openmpi/libhdf5.so" libhdf5_hl = "/usr/lib/x86_64-linux-gnu/hdf5/openmpi/libhdf5_hl.so" +[HDF5_jll] +libhdf5_hl_path = "/usr/lib/x86_64-linux-gnu/hdf5/openmpi/libhdf5_hl.so" +libhdf5_path = "/usr/lib/x86_64-linux-gnu/hdf5/openmpi/libhdf5.so" + [MPIPreferences] __clear__ = ["preloads_env_switch"] _format = "1.0" @@ -97,6 +110,22 @@ libsc = "/home/mschlott/hackathon/libtrixi/t8code/install/lib/libsc.so" libt8 = "/home/mschlott/hackathon/libtrixi/t8code/install/lib/libt8.so" ``` +This file is created with the following sequence of commands: +```julia +julia> using MPIPreferences +julia> MPIPreferences.use_system_binary() +``` +Restart the Julia REPL +```julia +julia> using P4est +julia> P4est.set_library_p4est!("/home/mschlott/hackathon/libtrixi/t8code/install/lib/libp4est.so") +julia> P4est.set_library_sc!("/home/mschlott/hackathon/libtrixi/t8code/install/lib/libsc.so") +julia> using T8code +julia> T8code.set_libraries_path!("/home/mschlott/hackathon/libtrixi/t8code/install/lib/") +julia> using HDF5 +julia> HDF5.API.set_libraries!("/usr/lib/x86_64-linux-gnu/hdf5/openmpi/libhdf5.so", "/usr/lib/x86_64-linux-gnu/hdf5/openmpi/libhdf5_hl.so") +``` +After the preferences are set, restart the Julia REPL again. ### [Usage](@id parallel_usage) @@ -218,7 +247,7 @@ julia> HDF5.API.set_libraries!("/path/to/your/libhdf5.so", "/path/to/your/libhdf ``` For more information see also the [documentation of HDF5.jl](https://juliaio.github.io/HDF5.jl/stable/mpi/). In total, you should -have a file called LocalPreferences.toml in the project directory that contains a section +have a file called `LocalPreferences.toml` in the project directory that contains a section `[MPIPreferences]`, a section `[HDF5]` with entries `libhdf5` and `libhdf5_hl`, a section `[P4est]` with the entry `libp4est` as well as a section `[T8code]` with the entries `libt8`, `libp4est` and `libsc`. diff --git a/examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl b/examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl new file mode 100644 index 00000000000..f87c0e056ca --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl @@ -0,0 +1,98 @@ +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) + +# 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) # minimum coordinates (min(x), min(y)) +coordinates_max = (0.0, 0.5) # maximum coordinates (max(x), max(y)) + +trees_per_dimension = (4, 4) +mesh = P4estMesh(trees_per_dimension, + polydeg = 3, initial_refinement_level = 2, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, + periodicity = false) + +# 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)) + +# 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, 0.5) +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 = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +# The AliveCallback prints short status information in regular intervals +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first), + base_level = 1, + med_level = 2, med_threshold = 0.9, + max_level = 3, max_threshold = 1.0) + +amr_callback = AMRCallback(semi, amr_controller, + interval = 50) + +# 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, amr_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, dt = 1e-7, 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_advection_diffusion_periodic_amr.jl b/examples/p4est_2d_dgsem/elixir_advection_diffusion_periodic_amr.jl new file mode 100644 index 00000000000..fdecb05f7bd --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_advection_diffusion_periodic_amr.jl @@ -0,0 +1,83 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection-diffusion equation + +advection_velocity = (1.5, 1.0) +equations = LinearScalarAdvectionEquation2D(advection_velocity) +diffusivity() = 5.0e-2 +equations_parabolic = LaplaceDiffusion2D(diffusivity(), equations) + +# 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)) + +trees_per_dimension = (4, 4) +mesh = P4estMesh(trees_per_dimension, + polydeg = 3, initial_refinement_level = 1, + coordinates_min = coordinates_min, coordinates_max = coordinates_max) + +# Define initial condition +function initial_condition_diffusive_convergence_test(x, t, + equation::LinearScalarAdvectionEquation2D) + # Store translated coordinate for easy use of exact solution + x_trans = x - equation.advection_velocity * t + + nu = diffusivity() + c = 1.0 + A = 0.5 + L = 2 + f = 1 / L + omega = 2 * pi * f + scalar = c + A * sin(omega * sum(x_trans)) * exp(-2 * nu * omega^2 * t) + return SVector(scalar) +end +initial_condition = initial_condition_diffusive_convergence_test + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolicParabolic(mesh, + (equations, equations_parabolic), + initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span `tspan` +tspan = (0.0, 0.5) +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) + +amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first), + base_level = 1, + med_level = 2, med_threshold = 1.25, + max_level = 3, max_threshold = 1.45) + +amr_callback = AMRCallback(semi, amr_controller, + interval = 20) + +# 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, amr_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_navierstokes_lid_driven_cavity.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl index e6566edb18a..bc28ae6ffb3 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl @@ -41,7 +41,6 @@ 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, diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity_amr.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity_amr.jl new file mode 100644 index 00000000000..898366969a4 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity_amr.jl @@ -0,0 +1,94 @@ +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 = (6, 6) +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) + +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 = 2000) +analysis_interval = 2000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +amr_indicator = IndicatorLöhner(semi, variable = Trixi.density) + +amr_controller = ControllerThreeLevel(semi, amr_indicator, + base_level = 0, + med_level = 1, med_threshold = 0.005, + max_level = 2, max_threshold = 0.01) + +amr_callback = AMRCallback(semi, amr_controller, + interval = 50, + adapt_initial_condition = true, + adapt_initial_condition_only_refine = true) + +callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback, amr_callback) +# 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/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl b/examples/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl index 7bd1ec0c647..70d76fc9075 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl @@ -40,7 +40,6 @@ 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 = boundary_condition_slip_wall boundary_conditions_parabolic = (; x_neg = boundary_condition_cavity, diff --git a/src/auxiliary/p4est.jl b/src/auxiliary/p4est.jl index 968af339cbd..0b826254129 100644 --- a/src/auxiliary/p4est.jl +++ b/src/auxiliary/p4est.jl @@ -13,14 +13,20 @@ This function will check if `p4est` is already initialized and if yes, do nothing, thus it is safe to call it multiple times. """ function init_p4est() - p4est_package_id = P4est.package_id() - if p4est_package_id >= 0 - return nothing + # Only initialize p4est if P4est.jl can be used + if P4est.preferences_set_correctly() + p4est_package_id = P4est.package_id() + if p4est_package_id >= 0 + return nothing + end + + # Initialize `p4est` with log level ERROR to prevent a lot of output in AMR simulations + p4est_init(C_NULL, SC_LP_ERROR) + else + @warn "Preferences for P4est.jl are not set correctly. Until fixed, using `P4estMesh` will result in a crash. " * + "See also https://trixi-framework.github.io/Trixi.jl/stable/parallelization/#parallel_system_MPI" end - # Initialize `p4est` with log level ERROR to prevent a lot of output in AMR simulations - p4est_init(C_NULL, SC_LP_ERROR) - return nothing end diff --git a/src/auxiliary/t8code.jl b/src/auxiliary/t8code.jl index 37cb782bb93..bd781b21c1e 100644 --- a/src/auxiliary/t8code.jl +++ b/src/auxiliary/t8code.jl @@ -7,34 +7,40 @@ is already initialized and if yes, do nothing, thus it is safe to call it multiple times. """ function init_t8code() - t8code_package_id = t8_get_package_id() - if t8code_package_id >= 0 - return nothing - end + # Only initialize t8code if T8code.jl can be used + if T8code.preferences_set_correctly() + t8code_package_id = t8_get_package_id() + if t8code_package_id >= 0 + return nothing + end - # Initialize the sc library, has to happen before we initialize t8code. - let catch_signals = 0, print_backtrace = 0, log_handler = C_NULL - T8code.Libt8.sc_init(mpi_comm(), catch_signals, print_backtrace, log_handler, - T8code.Libt8.SC_LP_ERROR) - end + # Initialize the sc library, has to happen before we initialize t8code. + let catch_signals = 0, print_backtrace = 0, log_handler = C_NULL + T8code.Libt8.sc_init(mpi_comm(), catch_signals, print_backtrace, log_handler, + T8code.Libt8.SC_LP_ERROR) + end - if T8code.Libt8.p4est_is_initialized() == 0 - # Initialize `p4est` with log level ERROR to prevent a lot of output in AMR simulations - T8code.Libt8.p4est_init(C_NULL, T8code.Libt8.SC_LP_ERROR) - end + if T8code.Libt8.p4est_is_initialized() == 0 + # Initialize `p4est` with log level ERROR to prevent a lot of output in AMR simulations + T8code.Libt8.p4est_init(C_NULL, T8code.Libt8.SC_LP_ERROR) + end - # Initialize t8code with log level ERROR to prevent a lot of output in AMR simulations. - t8_init(T8code.Libt8.SC_LP_ERROR) - - if haskey(ENV, "TRIXI_T8CODE_SC_FINALIZE") - # Normally, `sc_finalize` should always be called during shutdown of an - # application. It checks whether there is still un-freed memory by t8code - # and/or T8code.jl and throws an exception if this is the case. For - # production runs this is not mandatory, but is helpful during - # development. Hence, this option is only activated when environment - # variable TRIXI_T8CODE_SC_FINALIZE exists. - @warn "T8code.jl: sc_finalize will be called during shutdown of Trixi.jl." - MPI.add_finalize_hook!(T8code.Libt8.sc_finalize) + # Initialize t8code with log level ERROR to prevent a lot of output in AMR simulations. + t8_init(T8code.Libt8.SC_LP_ERROR) + + if haskey(ENV, "TRIXI_T8CODE_SC_FINALIZE") + # Normally, `sc_finalize` should always be called during shutdown of an + # application. It checks whether there is still un-freed memory by t8code + # and/or T8code.jl and throws an exception if this is the case. For + # production runs this is not mandatory, but is helpful during + # development. Hence, this option is only activated when environment + # variable TRIXI_T8CODE_SC_FINALIZE exists. + @warn "T8code.jl: sc_finalize will be called during shutdown of Trixi.jl." + MPI.add_finalize_hook!(T8code.Libt8.sc_finalize) + end + else + @warn "Preferences for T8code.jl are not set correctly. Until fixed, using `T8codeMesh` will result in a crash. " * + "See also https://trixi-framework.github.io/Trixi.jl/stable/parallelization/#parallel_system_MPI" end return nothing diff --git a/src/callbacks_step/amr.jl b/src/callbacks_step/amr.jl index ba840ff9675..5854c8617c3 100644 --- a/src/callbacks_step/amr.jl +++ b/src/callbacks_step/amr.jl @@ -528,6 +528,103 @@ function copy_to_quad_iter_volume(info, user_data) return nothing end +# specialized callback which includes the `cache_parabolic` argument +function (amr_callback::AMRCallback)(u_ode::AbstractVector, mesh::P4estMesh, + equations, dg::DG, cache, cache_parabolic, + semi, + t, iter; + only_refine = false, only_coarsen = false, + passive_args = ()) + @unpack controller, adaptor = amr_callback + + u = wrap_array(u_ode, mesh, equations, dg, cache) + lambda = @trixi_timeit timer() "indicator" controller(u, mesh, equations, dg, cache, + t = t, iter = iter) + + @boundscheck begin + @assert axes(lambda)==(Base.OneTo(ncells(mesh)),) ("Indicator array (axes = $(axes(lambda))) and mesh cells (axes = $(Base.OneTo(ncells(mesh)))) have different axes") + end + + # Copy controller value of each quad to the quad's user data storage + iter_volume_c = cfunction(copy_to_quad_iter_volume, Val(ndims(mesh))) + + # The pointer to lambda will be interpreted as Ptr{Int} below + @assert lambda isa Vector{Int} + iterate_p4est(mesh.p4est, lambda; iter_volume_c = iter_volume_c) + + @trixi_timeit timer() "refine" if !only_coarsen + # Refine mesh + refined_original_cells = @trixi_timeit timer() "mesh" refine!(mesh) + + # Refine solver + @trixi_timeit timer() "solver" refine!(u_ode, adaptor, mesh, equations, dg, + cache, cache_parabolic, + refined_original_cells) + for (p_u_ode, p_mesh, p_equations, p_dg, p_cache) in passive_args + @trixi_timeit timer() "passive solver" refine!(p_u_ode, adaptor, p_mesh, + p_equations, + p_dg, p_cache, + refined_original_cells) + end + else + # If there is nothing to refine, create empty array for later use + refined_original_cells = Int[] + end + + @trixi_timeit timer() "coarsen" if !only_refine + # Coarsen mesh + coarsened_original_cells = @trixi_timeit timer() "mesh" coarsen!(mesh) + + # coarsen solver + @trixi_timeit timer() "solver" coarsen!(u_ode, adaptor, mesh, equations, dg, + cache, cache_parabolic, + coarsened_original_cells) + for (p_u_ode, p_mesh, p_equations, p_dg, p_cache) in passive_args + @trixi_timeit timer() "passive solver" coarsen!(p_u_ode, adaptor, p_mesh, + p_equations, + p_dg, p_cache, + coarsened_original_cells) + end + else + # If there is nothing to coarsen, create empty array for later use + coarsened_original_cells = Int[] + end + + # Store whether there were any cells coarsened or refined and perform load balancing + has_changed = !isempty(refined_original_cells) || !isempty(coarsened_original_cells) + # Check if mesh changed on other processes + if mpi_isparallel() + has_changed = MPI.Allreduce!(Ref(has_changed), |, mpi_comm())[] + end + + if has_changed # TODO: Taal decide, where shall we set this? + # don't set it to has_changed since there can be changes from earlier calls + mesh.unsaved_changes = true + + if mpi_isparallel() && amr_callback.dynamic_load_balancing + @trixi_timeit timer() "dynamic load balancing" begin + global_first_quadrant = unsafe_wrap(Array, + unsafe_load(mesh.p4est).global_first_quadrant, + mpi_nranks() + 1) + old_global_first_quadrant = copy(global_first_quadrant) + partition!(mesh) + rebalance_solver!(u_ode, mesh, equations, dg, cache, + old_global_first_quadrant) + end + end + + reinitialize_boundaries!(semi.boundary_conditions, cache) + # if the semidiscretization also stores parabolic boundary conditions, + # reinitialize them after each refinement step as well. + if hasproperty(semi, :boundary_conditions_parabolic) + reinitialize_boundaries!(semi.boundary_conditions_parabolic, cache) + end + end + + # Return true if there were any cells coarsened or refined, otherwise false + return has_changed +end + # 2D function cfunction(::typeof(copy_to_quad_iter_volume), ::Val{2}) @cfunction(copy_to_quad_iter_volume, Cvoid, diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index 6395a9f348f..969f9c564f3 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -136,8 +136,8 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, P4estM return nothing end -# AMR for hyperbolic-parabolic equations currently only supported on TreeMeshes -function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, TreeMesh{3}}, +function refine!(u_ode::AbstractVector, adaptor, + mesh::Union{TreeMesh{2}, P4estMesh{2}, TreeMesh{3}}, equations, dg::DGSEM, cache, cache_parabolic, elements_to_refine) # Call `refine!` for the hyperbolic part, which does the heavy lifting of @@ -298,8 +298,8 @@ function coarsen!(u_ode::AbstractVector, adaptor, return nothing end -# AMR for hyperbolic-parabolic equations currently only supported on TreeMeshes -function coarsen!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, TreeMesh{3}}, +function coarsen!(u_ode::AbstractVector, adaptor, + mesh::Union{TreeMesh{2}, P4estMesh{2}, TreeMesh{3}}, equations, dg::DGSEM, cache, cache_parabolic, elements_to_remove) # Call `coarsen!` for the hyperbolic part, which does the heavy lifting of diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index 0176f5c6346..5fe68e06710 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -429,13 +429,17 @@ function reinitialize_containers!(mesh::P4estMesh, equations, dg::DGSEM, cache) @unpack boundaries = cache resize!(boundaries, required.boundaries) - # resize mortars container - @unpack mortars = cache - resize!(mortars, required.mortars) - - # re-initialize containers together to reduce - # the number of iterations over the mesh in `p4est` - init_surfaces!(interfaces, mortars, boundaries, mesh) + # re-initialize mortars container + if hasproperty(cache, :mortars) # cache_parabolic does not carry mortars + @unpack mortars = cache + resize!(mortars, required.mortars) + + # re-initialize containers together to reduce + # the number of iterations over the mesh in `p4est` + init_surfaces!(interfaces, mortars, boundaries, mesh) + else + init_surfaces!(interfaces, nothing, boundaries, mesh) + end end # A helper struct used in initialization methods below diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl index cf07645b949..9dd10df16ae 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl @@ -94,8 +94,19 @@ function rhs_parabolic!(du, u, t, mesh::P4estMesh{2}, dg.surface_integral, dg) end - # TODO: parabolic; extend to mortars - @assert nmortars(dg, cache) == 0 + # Prolong solution to mortars (specialized for AbstractEquationsParabolic) + # !!! NOTE: we reuse the hyperbolic cache here since it contains "mortars" and "u_threaded". See https://github.com/trixi-framework/Trixi.jl/issues/1674 for a discussion + @trixi_timeit timer() "prolong2mortars" begin + prolong2mortars_divergence!(cache, flux_viscous, mesh, equations_parabolic, + dg.mortar, dg.surface_integral, dg) + end + + # Calculate mortar fluxes (specialized for AbstractEquationsParabolic) + @trixi_timeit timer() "mortar flux" begin + calc_mortar_flux_divergence!(cache_parabolic.elements.surface_flux_values, + mesh, equations_parabolic, dg.mortar, + dg.surface_integral, dg, cache) + end # Calculate surface integrals @trixi_timeit timer() "surface integral" begin @@ -176,14 +187,15 @@ function calc_gradient!(gradients, u_transformed, t, end end - # Prolong solution to interfaces + # Prolong solution to interfaces. + # This reuses `prolong2interfaces` for the purely hyperbolic case. @trixi_timeit timer() "prolong2interfaces" begin prolong2interfaces!(cache_parabolic, u_transformed, mesh, equations_parabolic, dg.surface_integral, dg) end - # Calculate interface fluxes for the gradient. This reuses P4est `calc_interface_flux!` along with a - # specialization for AbstractEquationsParabolic. + # Calculate interface fluxes for the gradient. + # This reuses `calc_interface_flux!` for the purely hyperbolic case. @trixi_timeit timer() "interface flux" begin calc_interface_flux!(cache_parabolic.elements.surface_flux_values, mesh, False(), # False() = no nonconservative terms @@ -202,8 +214,21 @@ function calc_gradient!(gradients, u_transformed, t, mesh, equations_parabolic, dg.surface_integral, dg) end - # TODO: parabolic; mortars - @assert nmortars(dg, cache) == 0 + # Prolong solution to mortars. This resues the hyperbolic version of `prolong2mortars` + @trixi_timeit timer() "prolong2mortars" begin + prolong2mortars!(cache, u_transformed, mesh, equations_parabolic, + dg.mortar, dg.surface_integral, dg) + end + + # Calculate mortar fluxes. This reuses the hyperbolic version of `calc_mortar_flux`, + # along with a specialization on `calc_mortar_flux!(fstar, ...)` and `mortar_fluxes_to_elements!` for + # AbstractEquationsParabolic. + @trixi_timeit timer() "mortar flux" begin + calc_mortar_flux!(cache_parabolic.elements.surface_flux_values, + mesh, False(), # False() = no nonconservative terms + equations_parabolic, + dg.mortar, dg.surface_integral, dg, cache) + end # Calculate surface integrals @trixi_timeit timer() "surface integral" begin @@ -303,6 +328,64 @@ function calc_gradient!(gradients, u_transformed, t, return nothing end +# This version is called during `calc_gradients!` and must be specialized because the +# flux for the gradient is {u}, which doesn't depend on the outward normal. Thus, +# you don't need to scale by 2 (e.g., the scaling factor in the normals (and in the +# contravariant vectors) along large/small elements across a non-conforming +# interface in 2D) and flip the sign when storing the mortar fluxes back +# into `surface_flux_values`. +@inline function mortar_fluxes_to_elements!(surface_flux_values, + mesh::Union{P4estMesh{2}, T8codeMesh{2}}, + equations::AbstractEquationsParabolic, + mortar_l2::LobattoLegendreMortarL2, + dg::DGSEM, cache, mortar, fstar, u_buffer) + @unpack neighbor_ids, node_indices = cache.mortars + # Copy solution small to small + small_indices = node_indices[1, mortar] + small_direction = indices2direction(small_indices) + + for position in 1:2 + element = neighbor_ids[position, mortar] + for i in eachnode(dg) + for v in eachvariable(equations) + surface_flux_values[v, i, small_direction, element] = fstar[position][v, + i] + end + end + end + + # Project small fluxes to large element. + multiply_dimensionwise!(u_buffer, + mortar_l2.reverse_upper, fstar[2], + mortar_l2.reverse_lower, fstar[1]) + + # Copy interpolated flux values from buffer to large element face in the + # correct orientation. + # Note that the index of the small sides will always run forward but + # the index of the large side might need to run backwards for flipped sides. + large_element = neighbor_ids[3, mortar] + large_indices = node_indices[2, mortar] + large_direction = indices2direction(large_indices) + + if :i_backward in large_indices + for i in eachnode(dg) + for v in eachvariable(equations) + surface_flux_values[v, end + 1 - i, large_direction, large_element] = u_buffer[v, + i] + end + end + else + for i in eachnode(dg) + for v in eachvariable(equations) + surface_flux_values[v, i, large_direction, large_element] = u_buffer[v, + i] + end + end + end + + return nothing +end + # This version is used for parabolic gradient computations @inline function calc_interface_flux!(surface_flux_values, mesh::P4estMesh{2}, nonconservative_terms::False, @@ -321,7 +404,7 @@ end flux_ = 0.5 * (u_ll + u_rr) # we assume that the gradient computations utilize a central flux - # Note that we don't flip the sign on the secondondary flux. This is because for parabolic terms, + # Note that we don't flip the sign on the secondary flux. This is because for parabolic terms, # the normals are not embedded in `flux_` for the parabolic gradient computations. for v in eachvariable(equations) surface_flux_values[v, primary_node_index, primary_direction_index, primary_element_index] = flux_[v] @@ -520,6 +603,163 @@ function calc_interface_flux!(surface_flux_values, return nothing end +function prolong2mortars_divergence!(cache, flux_viscous::Vector{Array{uEltype, 4}}, + mesh::Union{P4estMesh{2}, T8codeMesh{2}}, equations, + mortar_l2::LobattoLegendreMortarL2, + surface_integral, dg::DGSEM) where {uEltype <: Real} + @unpack neighbor_ids, node_indices = cache.mortars + @unpack contravariant_vectors = cache.elements + index_range = eachnode(dg) + + flux_viscous_x, flux_viscous_y = flux_viscous + + @threaded for mortar in eachmortar(dg, cache) + # Copy solution data from the small elements using "delayed indexing" with + # a start value and a step size to get the correct face and orientation. + small_indices = node_indices[1, mortar] + direction_index = indices2direction(small_indices) + + i_small_start, i_small_step = index_to_start_step_2d(small_indices[1], + index_range) + j_small_start, j_small_step = index_to_start_step_2d(small_indices[2], + index_range) + + for position in 1:2 + i_small = i_small_start + j_small = j_small_start + element = neighbor_ids[position, mortar] + for i in eachnode(dg) + normal_direction = get_normal_direction(direction_index, + contravariant_vectors, + i_small, j_small, element) + + for v in eachvariable(equations) + flux_viscous = SVector(flux_viscous_x[v, i_small, j_small, element], + flux_viscous_y[v, i_small, j_small, element]) + + cache.mortars.u[1, v, position, i, mortar] = dot(flux_viscous, + normal_direction) + end + i_small += i_small_step + j_small += j_small_step + end + end + + # Buffer to copy solution values of the large element in the correct orientation + # before interpolating + u_buffer = cache.u_threaded[Threads.threadid()] + + # Copy solution of large element face to buffer in the + # correct orientation + large_indices = node_indices[2, mortar] + direction_index = indices2direction(large_indices) + + i_large_start, i_large_step = index_to_start_step_2d(large_indices[1], + index_range) + j_large_start, j_large_step = index_to_start_step_2d(large_indices[2], + index_range) + + i_large = i_large_start + j_large = j_large_start + element = neighbor_ids[3, mortar] + for i in eachnode(dg) + normal_direction = get_normal_direction(direction_index, contravariant_vectors, + i_large, j_large, element) + + for v in eachvariable(equations) + flux_viscous = SVector(flux_viscous_x[v, i_large, j_large, element], + flux_viscous_y[v, i_large, j_large, element]) + + # We prolong the viscous flux dotted with respect the outward normal + # on the small element. We scale by -1/2 here because the normal + # direction on the large element is negative 2x that of the small + # element (these normal directions are "scaled" by the surface Jacobian) + u_buffer[v, i] = -0.5 * dot(flux_viscous, normal_direction) + end + i_large += i_large_step + j_large += j_large_step + end + + # Interpolate large element face data from buffer to small face locations + multiply_dimensionwise!(view(cache.mortars.u, 2, :, 1, :, mortar), + mortar_l2.forward_lower, + u_buffer) + multiply_dimensionwise!(view(cache.mortars.u, 2, :, 2, :, mortar), + mortar_l2.forward_upper, + u_buffer) + end + + return nothing +end + +# We specialize `calc_mortar_flux!` for the divergence part of +# the parabolic terms. +function calc_mortar_flux_divergence!(surface_flux_values, + mesh::Union{P4estMesh{2}, T8codeMesh{2}}, + equations::AbstractEquationsParabolic, + mortar_l2::LobattoLegendreMortarL2, + surface_integral, dg::DG, cache) + @unpack neighbor_ids, node_indices = cache.mortars + @unpack contravariant_vectors = cache.elements + @unpack fstar_upper_threaded, fstar_lower_threaded = cache + index_range = eachnode(dg) + + @threaded for mortar in eachmortar(dg, cache) + # Choose thread-specific pre-allocated container + fstar = (fstar_lower_threaded[Threads.threadid()], + fstar_upper_threaded[Threads.threadid()]) + + for position in 1:2 + for node in eachnode(dg) + for v in eachvariable(equations) + viscous_flux_normal_ll = cache.mortars.u[1, v, position, node, mortar] + viscous_flux_normal_rr = cache.mortars.u[2, v, position, node, mortar] + + # TODO: parabolic; only BR1 at the moment + fstar[position][v, node] = 0.5 * (viscous_flux_normal_ll + + viscous_flux_normal_rr) + end + end + end + + # Buffer to interpolate flux values of the large element to before + # copying in the correct orientation + u_buffer = cache.u_threaded[Threads.threadid()] + + # this reuses the hyperbolic version of `mortar_fluxes_to_elements!` + mortar_fluxes_to_elements!(surface_flux_values, + mesh, equations, mortar_l2, dg, cache, + mortar, fstar, u_buffer) + end + + return nothing +end + +# We structure `calc_interface_flux!` similarly to "calc_mortar_flux!" for +# hyperbolic equations with no nonconservative terms. +# The reasoning is that parabolic fluxes are treated like conservative +# terms (e.g., we compute a viscous conservative "flux") and thus no +# non-conservative terms are present. +@inline function calc_mortar_flux!(fstar, + mesh::Union{P4estMesh{2}, T8codeMesh{2}}, + nonconservative_terms::False, + equations::AbstractEquationsParabolic, + surface_integral, dg::DG, cache, + mortar_index, position_index, normal_direction, + node_index) + @unpack u = cache.mortars + @unpack surface_flux = surface_integral + + u_ll, u_rr = get_surface_node_vars(u, equations, dg, position_index, node_index, + mortar_index) + + # TODO: parabolic; only BR1 at the moment + flux_ = 0.5 * (u_ll + u_rr) + + # Copy flux to buffer + set_node_vars!(fstar[position_index], flux_, equations, dg, node_index) +end + # TODO: parabolic, finish implementing `calc_boundary_flux_gradients!` and `calc_boundary_flux_divergence!` function prolong2boundaries!(cache_parabolic, flux_viscous, mesh::P4estMesh{2}, diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 152ca52a6ca..6632cd0bb27 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -541,6 +541,38 @@ end end end +@trixi_testset "P4estMesh2D: elixir_advection_diffusion_periodic_amr.jl" begin + @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", + "elixir_advection_diffusion_periodic_amr.jl"), + tspan=(0.0, 0.01), + l2=[0.014715887539773128], + linf=[0.2285802791900049]) + # 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 + +@trixi_testset "P4estMesh2D: elixir_advection_diffusion_nonperiodic_amr.jl" begin + @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", + "elixir_advection_diffusion_nonperiodic_amr.jl"), + tspan=(0.0, 0.01), + l2=[0.00793438523666649], + linf=[0.11030633127144573]) + # 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 + @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"), @@ -635,6 +667,28 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + +@trixi_testset "P4estMesh2D: elixir_navierstokes_lid_driven_cavity_amr.jl" begin + @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", + "elixir_navierstokes_lid_driven_cavity_amr.jl"), + tspan=(0.0, 1.0), + l2=[ + 0.0005323841980601085, 0.07892044543547208, + 0.02909671646389337, 0.11717468256112017, + ], + linf=[ + 0.006045292737899444, 0.9233292581786228, + 0.7982129977236198, 1.6864546235292153, + ]) + # 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 # Clean up afterwards: delete Trixi.jl output directory