From b9ace6d7ff34eec684eaaba337f66fdfa571a99a Mon Sep 17 00:00:00 2001 From: Huiyu Xie Date: Fri, 16 Aug 2024 04:23:48 -1000 Subject: [PATCH 1/6] Add numerical support of other real types (`polytropic_euler`) (#2015) * start * fix flux tests * keep conflicts * complete equations * complete unit test --- src/equations/polytropic_euler_2d.jl | 52 +++++++++---------- test/test_type.jl | 74 ++++++++++++++++++++++++++++ 2 files changed, 101 insertions(+), 25 deletions(-) diff --git a/src/equations/polytropic_euler_2d.jl b/src/equations/polytropic_euler_2d.jl index e900fd64073..571d4723f6f 100644 --- a/src/equations/polytropic_euler_2d.jl +++ b/src/equations/polytropic_euler_2d.jl @@ -62,7 +62,7 @@ in combination with [`source_terms_convergence_test`](@ref). function initial_condition_convergence_test(x, t, equations::PolytropicEulerEquations2D) # manufactured initial condition from Winters (2019) [0.1007/s10543-019-00789-w] # domain must be set to [0, 1] x [0, 1] - h = 8 + cos(2 * pi * x[1]) * sin(2 * pi * x[2]) * cos(2 * pi * t) + h = 8 + cospi(2 * x[1]) * sinpi(2 * x[2]) * cospi(2 * t) return SVector(h, h / 2, 3 * h / 2) end @@ -78,19 +78,20 @@ Source terms used for convergence tests in combination with rho, v1, v2 = cons2prim(u, equations) # Residual from Winters (2019) [0.1007/s10543-019-00789-w] eq. (5.2). - h = 8 + cos(2 * pi * x[1]) * sin(2 * pi * x[2]) * cos(2 * pi * t) - h_t = -2 * pi * cos(2 * pi * x[1]) * sin(2 * pi * x[2]) * sin(2 * pi * t) - h_x = -2 * pi * sin(2 * pi * x[1]) * sin(2 * pi * x[2]) * cos(2 * pi * t) - h_y = 2 * pi * cos(2 * pi * x[1]) * cos(2 * pi * x[2]) * cos(2 * pi * t) + RealT = eltype(u) + h = 8 + cospi(2 * x[1]) * sinpi(2 * x[2]) * cospi(2 * t) + h_t = -2 * convert(RealT, pi) * cospi(2 * x[1]) * sinpi(2 * x[2]) * sinpi(2 * t) + h_x = -2 * convert(RealT, pi) * sinpi(2 * x[1]) * sinpi(2 * x[2]) * cospi(2 * t) + h_y = 2 * convert(RealT, pi) * cospi(2 * x[1]) * cospi(2 * x[2]) * cospi(2 * t) rho_x = h_x rho_y = h_y b = equations.kappa * equations.gamma * h^(equations.gamma - 1) - r_1 = h_t + h_x / 2 + 3 / 2 * h_y - r_2 = h_t / 2 + h_x / 4 + b * rho_x + 3 / 4 * h_y - r_3 = 3 / 2 * h_t + 3 / 4 * h_x + 9 / 4 * h_y + b * rho_y + r_1 = h_t + h_x / 2 + 3 * h_y / 2 + r_2 = h_t / 2 + h_x / 4 + b * rho_x + 3 * h_y / 4 + r_3 = 3 * h_t / 2 + 3 * h_x / 4 + 9 * h_y / 4 + b * rho_y return SVector(r_1, r_2, r_3) end @@ -113,9 +114,10 @@ function initial_condition_weak_blast_wave(x, t, equations::PolytropicEulerEquat phi = atan(y_norm, x_norm) # Calculate primitive variables - rho = r > 0.5 ? 1.0 : 1.1691 - v1 = r > 0.5 ? 0.0 : 0.1882 * cos(phi) - v2 = r > 0.5 ? 0.0 : 0.1882 * sin(phi) + RealT = eltype(x) + rho = r > 0.5f0 ? one(RealT) : convert(RealT, 1.1691) + v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos(phi) + v2 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * sin(phi) return prim2cons(SVector(rho, v1, v2), equations) end @@ -181,17 +183,17 @@ For details see Section 3.2 of the following reference v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] # Compute the necessary mean values - if equations.gamma == 1.0 # isothermal gas + if equations.gamma == 1 # isothermal gas rho_mean = ln_mean(rho_ll, rho_rr) else # equations.gamma > 1 # polytropic gas rho_mean = stolarsky_mean(rho_ll, rho_rr, equations.gamma) end - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - p_avg = 0.5 * (p_ll + p_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + p_avg = 0.5f0 * (p_ll + p_rr) # Calculate fluxes depending on normal_direction - f1 = rho_mean * 0.5 * (v_dot_n_ll + v_dot_n_rr) + f1 = rho_mean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr) f2 = f1 * v1_avg + p_avg * normal_direction[1] f3 = f1 * v2_avg + p_avg * normal_direction[2] @@ -207,21 +209,21 @@ end p_rr = equations.kappa * rho_rr^equations.gamma # Compute the necessary mean values - if equations.gamma == 1.0 # isothermal gas + if equations.gamma == 1 # isothermal gas rho_mean = ln_mean(rho_ll, rho_rr) else # equations.gamma > 1 # polytropic gas rho_mean = stolarsky_mean(rho_ll, rho_rr, equations.gamma) end - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - p_avg = 0.5 * (p_ll + p_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + p_avg = 0.5f0 * (p_ll + p_rr) if orientation == 1 # x-direction - f1 = rho_mean * 0.5 * (v1_ll + v1_rr) + f1 = rho_mean * 0.5f0 * (v1_ll + v1_rr) f2 = f1 * v1_avg + p_avg f3 = f1 * v2_avg else # y-direction - f1 = rho_mean * 0.5 * (v2_ll + v2_rr) + f1 = rho_mean * 0.5f0 * (v2_ll + v2_rr) f2 = f1 * v1_avg f3 = f1 * v2_avg + p_avg end @@ -360,14 +362,14 @@ end v_square = v1^2 + v2^2 p = pressure(u, equations) # Form of the internal energy depends on gas type - if equations.gamma == 1.0 # isothermal gas + if equations.gamma == 1 # isothermal gas internal_energy = equations.kappa * log(rho) else # equations.gamma > 1 # polytropic gas internal_energy = equations.kappa * rho^(equations.gamma - 1) / - (equations.gamma - 1.0) + (equations.gamma - 1) end - w1 = internal_energy + p / rho - 0.5 * v_square + w1 = internal_energy + p / rho - 0.5f0 * v_square w2 = v1 w3 = v2 diff --git a/test/test_type.jl b/test/test_type.jl index 9b8c3366cfb..d22afa65c0a 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -1926,6 +1926,80 @@ isdir(outdir) && rm(outdir, recursive = true) end end + @timed_testset "Polytropic Euler 2D" begin + for RealT in (Float32, Float64) + equations1 = @inferred PolytropicEulerEquations2D(RealT(1), + RealT(1)) # equations.gamma == 1 + equations2 = @inferred PolytropicEulerEquations2D(RealT(1.4), RealT(0.5)) # equations.gamma > 1 + + for equations in (equations1, equations2) + x = SVector(zero(RealT), zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = SVector(one(RealT), one(RealT), one(RealT)) + orientations = [1, 2] + directions = [1, 2, 3, 4] + normal_direction = SVector(one(RealT), zero(RealT)) + + @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) == + RealT + @test eltype(@inferred source_terms_convergence_test(u, x, t, equations)) == + RealT + @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) == + RealT + + @test eltype(@inferred flux(u, normal_direction, equations)) == RealT + if RealT == Float32 + # check `ln_mean` and `stolarsky_mean` (test broken) + @test_broken eltype(@inferred flux_winters_etal(u_ll, u_rr, + normal_direction, + equations)) == + RealT + else + @test eltype(@inferred flux_winters_etal(u_ll, u_rr, normal_direction, + equations)) == + RealT + end + @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, normal_direction, + equations)) == + RealT + @test eltype(@inferred min_max_speed_davis(u_ll, u_rr, normal_direction, + equations)) == + RealT + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, + equations)) == + RealT + + for orientation in orientations + @test eltype(@inferred flux(u, orientation, equations)) == RealT + if RealT == Float32 + # check `ln_mean` and `stolarsky_mean` (test broken) + @test_broken eltype(@inferred flux_winters_etal(u_ll, u_rr, + orientation, + equations)) == + RealT + else + @test eltype(@inferred flux_winters_etal(u_ll, u_rr, orientation, + equations)) == + RealT + end + @test eltype(@inferred min_max_speed_davis(u_ll, u_rr, orientation, + equations)) == + RealT + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, + equations)) == + RealT + end + + @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred prim2cons(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test typeof(@inferred density(u, equations)) == RealT + @test typeof(@inferred pressure(u, equations)) == RealT + end + end + end + @timed_testset "Shallow Water 1D" begin for RealT in (Float32, Float64) equations = @inferred ShallowWaterEquations1D(gravity_constant = RealT(9.81)) From b0d3a44dc55b6b61148947cc9fc2d8832fb3a33f Mon Sep 17 00:00:00 2001 From: Patrick Ersing <114223904+patrickersing@users.noreply.github.com> Date: Mon, 19 Aug 2024 10:47:46 +0200 Subject: [PATCH 2/6] Update NC-fluxes for the SWE (#2038) * rewrite nonconservative fluxes * update test values for unstructured2d * update test values for unstructured 2d --------- Co-authored-by: Andrew Winters --- .../tree_2d_dgsem/elixir_shallowwater_wall.jl | 4 +- src/Trixi.jl | 1 - src/equations/shallow_water_1d.jl | 69 ++------ src/equations/shallow_water_2d.jl | 151 +++--------------- test/test_tree_1d_shallowwater.jl | 32 ++-- test/test_tree_2d_shallowwater.jl | 26 ++- test/test_type.jl | 10 -- test/test_unstructured_2d.jl | 40 ++--- 8 files changed, 88 insertions(+), 245 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl b/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl index b8dbad50680..f8f601d4120 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl @@ -29,8 +29,8 @@ boundary_condition = boundary_condition_slip_wall ############################################################################### # Get the DG approximation space -volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) -surface_flux = (flux_lax_friedrichs, flux_nonconservative_ersing_etal) +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +surface_flux = (flux_lax_friedrichs, flux_nonconservative_wintermeyer_etal) solver = DGSEM(polydeg = 3, surface_flux = surface_flux, volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) diff --git a/src/Trixi.jl b/src/Trixi.jl index 23a8cfe1d0a..90fa590c50a 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -179,7 +179,6 @@ export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle, flux_kennedy_gruber, flux_shima_etal, flux_ec, flux_fjordholm_etal, flux_nonconservative_fjordholm_etal, flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal, - flux_nonconservative_ersing_etal, flux_chan_etal, flux_nonconservative_chan_etal, flux_winters_etal, hydrostatic_reconstruction_audusse_etal, flux_nonconservative_audusse_etal, FluxPlusDissipation, DissipationGlobalLaxFriedrichs, DissipationLocalLaxFriedrichs, diff --git a/src/equations/shallow_water_1d.jl b/src/equations/shallow_water_1d.jl index 998deb04de2..3c218eee9d9 100644 --- a/src/equations/shallow_water_1d.jl +++ b/src/equations/shallow_water_1d.jl @@ -201,20 +201,27 @@ end Non-symmetric two-point volume flux discretizing the nonconservative (source) term that contains the gradient of the bottom topography [`ShallowWaterEquations1D`](@ref). -Further details are available in the paper: +Gives entropy conservation and well-balancedness on both the volume and surface when combined with +[`flux_wintermeyer_etal`](@ref). + +Further details are available in the papers: - Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017) An entropy stable nodal discontinuous Galerkin method for the two dimensional shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry [DOI: 10.1016/j.jcp.2017.03.036](https://doi.org/10.1016/j.jcp.2017.03.036) +- Patrick Ersing, Andrew R. Winters (2023) + An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on + curvilinear meshes + [DOI: 10.48550/arXiv.2306.12699](https://doi.org/10.48550/arXiv.2306.12699) """ @inline function flux_nonconservative_wintermeyer_etal(u_ll, u_rr, orientation::Integer, equations::ShallowWaterEquations1D) # Pull the necessary left and right state information h_ll = waterheight(u_ll, equations) - b_rr = u_rr[3] + b_jump = u_rr[3] - u_ll[3] # Bottom gradient nonconservative term: (0, g h b_x, 0) - f = SVector(0, equations.gravity * h_ll * b_rr, 0) + f = SVector(0, equations.gravity * h_ll * b_jump, 0) return f end @@ -226,11 +233,8 @@ end Non-symmetric two-point surface flux discretizing the nonconservative (source) term of that contains the gradient of the bottom topography [`ShallowWaterEquations1D`](@ref). -This contains additional terms compared to [`flux_nonconservative_wintermeyer_etal`](@ref) -that account for possible discontinuities in the bottom topography function. -Thus, this flux should be used in general at interfaces. For flux differencing volume terms, -[`flux_nonconservative_wintermeyer_etal`](@ref) is analytically equivalent but slightly -cheaper. +This flux can be used together with [`flux_fjordholm_etal`](@ref) at interfaces to ensure entropy +conservation and well-balancedness. Further details for the original finite volume formulation are available in - Ulrik S. Fjordholm, Siddhartha Mishr and Eitan Tadmor (2011) @@ -256,7 +260,6 @@ and for curvilinear 2D case in the paper: # cross-averaging across a discontinuous bottom topography # (ii) True surface part that uses `h_average` and `b_jump` to handle discontinuous bathymetry f = SVector(0, - equations.gravity * h_ll * b_ll + equations.gravity * h_average * b_jump, 0) @@ -285,7 +288,7 @@ Further details on the hydrostatic reconstruction and its motivation can be foun orientation::Integer, equations::ShallowWaterEquations1D) # Pull the water height and bottom topography on the left - h_ll, _, b_ll = u_ll + h_ll, _, _ = u_ll # Create the hydrostatic reconstruction for the left solution state u_ll_star, _ = hydrostatic_reconstruction_audusse_etal(u_ll, u_rr, equations) @@ -293,52 +296,11 @@ Further details on the hydrostatic reconstruction and its motivation can be foun # Copy the reconstructed water height for easier to read code h_ll_star = u_ll_star[1] - # Includes two parts: - # (i) Diagonal (consistent) term from the volume flux that uses `b_ll` to avoid - # cross-averaging across a discontinuous bottom topography - # (ii) True surface part that uses `h_ll` and `h_ll_star` to handle discontinuous bathymetry return SVector(0, - equations.gravity * h_ll * b_ll + equations.gravity * (h_ll^2 - h_ll_star^2), 0) end -""" - flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer, - equations::ShallowWaterEquations1D) - -!!! warning "Experimental code" - This numerical flux is experimental and may change in any future release. - -Non-symmetric path-conservative two-point volume flux discretizing the nonconservative (source) term -that contains the gradient of the bottom topography [`ShallowWaterEquations1D`](@ref). - -This is a modified version of [`flux_nonconservative_wintermeyer_etal`](@ref) that gives entropy -conservation and well-balancedness in both the volume and surface when combined with -[`flux_wintermeyer_etal`](@ref). - -For further details see: -- Patrick Ersing, Andrew R. Winters (2023) - An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on - curvilinear meshes - [DOI: 10.48550/arXiv.2306.12699](https://doi.org/10.48550/arXiv.2306.12699) -""" -@inline function flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer, - equations::ShallowWaterEquations1D) - # Pull the necessary left and right state information - h_ll = waterheight(u_ll, equations) - b_rr = u_rr[3] - b_ll = u_ll[3] - - # Calculate jump - b_jump = b_rr - b_ll - - # Bottom gradient nonconservative term: (0, g h b_x, 0) - f = SVector(0, equations.gravity * h_ll * b_jump, 0) - - return f -end - """ flux_fjordholm_etal(u_ll, u_rr, orientation, equations::ShallowWaterEquations1D) @@ -378,7 +340,8 @@ end Total energy conservative (mathematical entropy for 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). +For the `surface_flux` either [`flux_wintermeyer_etal`](@ref) or [`flux_fjordholm_etal`](@ref) can +be used to ensure well-balancedness and entropy conservation. Further details are available in Theorem 1 of the paper: - Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017) @@ -571,7 +534,7 @@ end # Note, only the first two are the entropy variables, the third entry still # just carries the bottom topography values for convenience @inline function cons2entropy(u, equations::ShallowWaterEquations1D) - h, h_v, b = u + h, _, b = u v = velocity(u, equations) diff --git a/src/equations/shallow_water_2d.jl b/src/equations/shallow_water_2d.jl index db8f00fc15b..4ecaf3b6e14 100644 --- a/src/equations/shallow_water_2d.jl +++ b/src/equations/shallow_water_2d.jl @@ -274,28 +274,30 @@ end Non-symmetric two-point volume flux discretizing the nonconservative (source) term that contains the gradient of the bottom topography [`ShallowWaterEquations2D`](@ref). -On curvilinear meshes, this nonconservative flux depends on both the -contravariant vector (normal direction) at the current node and the averaged -one. This is different from numerical fluxes used to discretize conservative -terms. +For the `surface_flux` either [`flux_wintermeyer_etal`](@ref) or [`flux_fjordholm_etal`](@ref) can +be used to ensure well-balancedness and entropy conservation. -Further details are available in the paper: +Further details are available in the papers: - Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017) An entropy stable nodal discontinuous Galerkin method for the two dimensional shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry [DOI: 10.1016/j.jcp.2017.03.036](https://doi.org/10.1016/j.jcp.2017.03.036) +- Patrick Ersing, Andrew R. Winters (2023) + An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on + curvilinear meshes + [DOI: 10.48550/arXiv.2306.12699](https://doi.org/10.48550/arXiv.2306.12699) """ @inline function flux_nonconservative_wintermeyer_etal(u_ll, u_rr, orientation::Integer, equations::ShallowWaterEquations2D) # Pull the necessary left and right state information h_ll = waterheight(u_ll, equations) - b_rr = u_rr[4] + b_jump = u_rr[4] - u_ll[4] # Bottom gradient nonconservative term: (0, g h b_x, g h b_y, 0) if orientation == 1 - f = SVector(0, equations.gravity * h_ll * b_rr, 0, 0) + f = SVector(0, equations.gravity * h_ll * b_jump, 0, 0) else # orientation == 2 - f = SVector(0, 0, equations.gravity * h_ll * b_rr, 0) + f = SVector(0, 0, equations.gravity * h_ll * b_jump, 0) end return f end @@ -306,12 +308,12 @@ end equations::ShallowWaterEquations2D) # Pull the necessary left and right state information h_ll = waterheight(u_ll, equations) - b_rr = u_rr[4] - # Note this routine only uses the `normal_direction_average` and the average of the - # bottom topography to get a quadratic split form DG gradient on curved elements + b_jump = u_rr[4] - u_ll[4] + + # Bottom gradient nonconservative term: (0, g h b_x, g h b_y, 0) return SVector(0, - normal_direction_average[1] * equations.gravity * h_ll * b_rr, - normal_direction_average[2] * equations.gravity * h_ll * b_rr, + normal_direction_average[1] * equations.gravity * h_ll * b_jump, + normal_direction_average[2] * equations.gravity * h_ll * b_jump, 0) end @@ -326,16 +328,8 @@ end Non-symmetric two-point surface flux discretizing the nonconservative (source) term of that contains the gradient of the bottom topography [`ShallowWaterEquations2D`](@ref). -On curvilinear meshes, this nonconservative flux depends on both the -contravariant vector (normal direction) at the current node and the averaged -one. This is different from numerical fluxes used to discretize conservative -terms. - -This contains additional terms compared to [`flux_nonconservative_wintermeyer_etal`](@ref) -that account for possible discontinuities in the bottom topography function. -Thus, this flux should be used in general at interfaces. For flux differencing volume terms, -[`flux_nonconservative_wintermeyer_etal`](@ref) is analytically equivalent but slightly -cheaper. +This flux can be used together with [`flux_fjordholm_etal`](@ref) at interfaces to ensure entropy +conservation and well-balancedness. Further details for the original finite volume formulation are available in - Ulrik S. Fjordholm, Siddhartha Mishr and Eitan Tadmor (2011) @@ -356,18 +350,13 @@ and for curvilinear 2D case in the paper: h_average = 0.5f0 * (h_ll + h_rr) b_jump = b_rr - b_ll - # Includes two parts: - # (i) Diagonal (consistent) term from the volume flux that uses `b_ll` to avoid - # cross-averaging across a discontinuous bottom topography - # (ii) True surface part that uses `h_average` and `b_jump` to handle discontinuous bathymetry + # Bottom gradient nonconservative term: (0, g h b_x, g h b_y, 0) if orientation == 1 f = SVector(0, - equations.gravity * h_ll * b_ll + equations.gravity * h_average * b_jump, 0, 0) else # orientation == 2 f = SVector(0, 0, - equations.gravity * h_ll * b_ll + equations.gravity * h_average * b_jump, 0) end @@ -383,20 +372,12 @@ end h_ll, _, _, b_ll = u_ll h_rr, _, _, b_rr = u_rr - # Comes in two parts: - # (i) Diagonal (consistent) term from the volume flux that uses `normal_direction_average` - # but we use `b_ll` to avoid cross-averaging across a discontinuous bottom topography - - f2 = normal_direction_average[1] * equations.gravity * h_ll * b_ll - f3 = normal_direction_average[2] * equations.gravity * h_ll * b_ll - - # (ii) True surface part that uses `normal_direction_ll`, `h_average` and `b_jump` - # to handle discontinuous bathymetry h_average = 0.5f0 * (h_ll + h_rr) b_jump = b_rr - b_ll - f2 += normal_direction_ll[1] * equations.gravity * h_average * b_jump - f3 += normal_direction_ll[2] * equations.gravity * h_average * b_jump + # Bottom gradient nonconservative term: (0, g h b_x, g h b_y, 0) + f2 = normal_direction_average[1] * equations.gravity * h_average * b_jump + f3 = normal_direction_average[2] * equations.gravity * h_average * b_jump # First and last equations do not have a nonconservative flux f1 = f4 = 0 @@ -472,18 +453,12 @@ Further details for the hydrostatic reconstruction and its motivation can be fou # Copy the reconstructed water height for easier to read code h_ll_star = u_ll_star[1] - # Includes two parts: - # (i) Diagonal (consistent) term from the volume flux that uses `b_ll` to avoid - # cross-averaging across a discontinuous bottom topography - # (ii) True surface part that uses `h_ll` and `h_ll_star` to handle discontinuous bathymetry if orientation == 1 f = SVector(0, - equations.gravity * h_ll * b_ll + equations.gravity * (h_ll^2 - h_ll_star^2), 0, 0) else # orientation == 2 f = SVector(0, 0, - equations.gravity * h_ll * b_ll + equations.gravity * (h_ll^2 - h_ll_star^2), 0) end @@ -504,18 +479,8 @@ end # Copy the reconstructed water height for easier to read code h_ll_star = u_ll_star[1] - # Comes in two parts: - # (i) Diagonal (consistent) term from the volume flux that uses `normal_direction_average` - # but we use `b_ll` to avoid cross-averaging across a discontinuous bottom topography - - f2 = normal_direction_average[1] * equations.gravity * h_ll * b_ll - f3 = normal_direction_average[2] * equations.gravity * h_ll * b_ll - - # (ii) True surface part that uses `normal_direction_ll`, `h_ll` and `h_ll_star` - # to handle discontinuous bathymetry - - f2 += normal_direction_ll[1] * equations.gravity * (h_ll^2 - h_ll_star^2) - f3 += normal_direction_ll[2] * equations.gravity * (h_ll^2 - h_ll_star^2) + f2 = normal_direction_average[1] * equations.gravity * (h_ll^2 - h_ll_star^2) + f3 = normal_direction_average[2] * equations.gravity * (h_ll^2 - h_ll_star^2) # First and last equations do not have a nonconservative flux f1 = f4 = 0 @@ -523,73 +488,6 @@ end return SVector(f1, f2, f3, f4) end -""" - flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer, - equations::ShallowWaterEquations2D) - flux_nonconservative_ersing_etal(u_ll, u_rr, - normal_direction_ll::AbstractVector, - normal_direction_average::AbstractVector, - equations::ShallowWaterEquations2D) - -!!! warning "Experimental code" - This numerical flux is experimental and may change in any future release. - -Non-symmetric path-conservative two-point volume flux discretizing the nonconservative (source) term -that contains the gradient of the bottom topography [`ShallowWaterEquations2D`](@ref). - -On curvilinear meshes, this nonconservative flux depends on both the -contravariant vector (normal direction) at the current node and the averaged -one. This is different from numerical fluxes used to discretize conservative -terms. - -This is a modified version of [`flux_nonconservative_wintermeyer_etal`](@ref) that gives entropy -conservation and well-balancedness in both the volume and surface when combined with -[`flux_wintermeyer_etal`](@ref). - -For further details see: -- Patrick Ersing, Andrew R. Winters (2023) - An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on - curvilinear meshes - [DOI: 10.48550/arXiv.2306.12699](https://doi.org/10.48550/arXiv.2306.12699) -""" -@inline function flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer, - equations::ShallowWaterEquations2D) - # Pull the necessary left and right state information - h_ll = waterheight(u_ll, equations) - b_rr = u_rr[4] - b_ll = u_ll[4] - - # Calculate jump - b_jump = b_rr - b_ll - - # Bottom gradient nonconservative term: (0, g h b_x, g h b_y, 0) - if orientation == 1 - f = SVector(0, equations.gravity * h_ll * b_jump, 0, 0) - else # orientation == 2 - f = SVector(0, 0, equations.gravity * h_ll * b_jump, 0) - end - return f -end - -@inline function flux_nonconservative_ersing_etal(u_ll, u_rr, - normal_direction_ll::AbstractVector, - normal_direction_average::AbstractVector, - equations::ShallowWaterEquations2D) - # Pull the necessary left and right state information - h_ll = waterheight(u_ll, equations) - b_rr = u_rr[4] - b_ll = u_ll[4] - - # Calculate jump - b_jump = b_rr - b_ll - # Note this routine only uses the `normal_direction_average` and the average of the - # bottom topography to get a quadratic split form DG gradient on curved elements - return SVector(0, - normal_direction_average[1] * equations.gravity * h_ll * b_jump, - normal_direction_average[2] * equations.gravity * h_ll * b_jump, - 0) -end - """ flux_fjordholm_etal(u_ll, u_rr, orientation_or_normal_direction, equations::ShallowWaterEquations2D) @@ -664,7 +562,8 @@ end Total energy conservative (mathematical entropy for 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). +For the `surface_flux` either [`flux_wintermeyer_etal`](@ref) or [`flux_fjordholm_etal`](@ref) can +be used to ensure well-balancedness and entropy conservation. Further details are available in Theorem 1 of the paper: - Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017) diff --git a/test/test_tree_1d_shallowwater.jl b/test/test_tree_1d_shallowwater.jl index 8fe3291a938..42a91e578e4 100644 --- a/test/test_tree_1d_shallowwater.jl +++ b/test/test_tree_1d_shallowwater.jl @@ -98,7 +98,7 @@ end end end -@trixi_testset "elixir_shallowwater_well_balanced.jl with flux_nonconservative_ersing_etal" begin +@trixi_testset "elixir_shallowwater_well_balanced.jl with flux_nonconservative_wintermeyer_etal" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), l2=[ 0.10416666834254838, @@ -107,9 +107,7 @@ end ], linf=[2.0000000000000004, 3.0610625110157164e-14, 2.0], surface_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), - volume_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), + flux_nonconservative_wintermeyer_etal), tspan=(0.0, 0.25)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -169,7 +167,7 @@ end end end -@trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_ersing_etal" begin +@trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_wintermeyer_etal" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), l2=[ 0.005774284062933275, @@ -182,9 +180,7 @@ end 9.098379777450205e-5, ], surface_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), - volume_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), + flux_nonconservative_wintermeyer_etal), tspan=(0.0, 0.025)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -200,13 +196,13 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms_dirichlet.jl"), l2=[ - 0.0022851099219788917, - 0.01560453773635554, - 4.43649172558535e-5, + 0.0022667320585353927, + 0.01571629729279524, + 4.4364917255842716e-5, ], linf=[ - 0.008934615705174398, - 0.059403169140869405, + 0.008945234652224965, + 0.059403165802872415, 9.098379777405796e-5, ], tspan=(0.0, 0.025)) @@ -224,13 +220,13 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms_dirichlet.jl"), l2=[ - 0.0022956052733432287, - 0.015540053559855601, - 4.43649172558535e-5, + 0.0022774071143995952, + 0.01566214422689219, + 4.4364917255842716e-5, ], linf=[ - 0.008460440313118323, - 0.05720939349382359, + 0.008451721489057373, + 0.05720939055279217, 9.098379777405796e-5, ], surface_flux=(FluxHydrostaticReconstruction(FluxHLL(min_max_speed_naive), diff --git a/test/test_tree_2d_shallowwater.jl b/test/test_tree_2d_shallowwater.jl index 9a3ba36c7d4..bcad663008c 100644 --- a/test/test_tree_2d_shallowwater.jl +++ b/test/test_tree_2d_shallowwater.jl @@ -114,7 +114,7 @@ end end end -@trixi_testset "elixir_shallowwater_well_balanced.jl with flux_nonconservative_ersing_etal" begin +@trixi_testset "elixir_shallowwater_well_balanced.jl with flux_nonconservative_wintermeyer_etal" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), l2=[ 0.9130579602987146, @@ -129,9 +129,7 @@ end 2.1130620376156584, ], surface_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), - volume_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), + flux_nonconservative_wintermeyer_etal), tspan=(0.0, 0.25)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -172,15 +170,15 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms_dirichlet.jl"), l2=[ - 0.0018746929418489125, - 0.017332321628469628, - 0.01634953679145536, - 6.274146767717023e-5, + 0.0018596727473552813, + 0.017306217777629147, + 0.016367646997420396, + 6.274146767723934e-5, ], linf=[ - 0.016262353691956388, - 0.08726160620859424, - 0.09043621801418844, + 0.016548007102923368, + 0.08726160568822783, + 0.09043621622245013, 0.0001819675955490041, ], tspan=(0.0, 0.025)) @@ -248,7 +246,7 @@ end end end -@trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_ersing_etal" begin +@trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_wintermeyer_etal" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), l2=[ 0.002471853426064005, @@ -263,9 +261,7 @@ end 0.0001819675955490041, ], surface_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), - volume_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), + flux_nonconservative_wintermeyer_etal), tspan=(0.0, 0.25)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) diff --git a/test/test_type.jl b/test/test_type.jl index d22afa65c0a..c15ac5f78c6 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -2045,8 +2045,6 @@ isdir(outdir) && rm(outdir, recursive = true) @test eltype(@inferred flux_nonconservative_audusse_etal(u_ll, u_rr, orientation, equations)) == RealT - @test eltype(@inferred flux_nonconservative_ersing_etal(u_ll, u_rr, orientation, - equations)) == RealT @test eltype(@inferred flux_fjordholm_etal(u_ll, u_rr, orientation, equations)) == RealT @test eltype(@inferred flux_wintermeyer_etal(u_ll, u_rr, orientation, @@ -2133,10 +2131,6 @@ isdir(outdir) && rm(outdir, recursive = true) normal_direction_ll, normal_direction_average, equations)) == RealT - @test eltype(@inferred flux_nonconservative_ersing_etal(u_ll, u_rr, - normal_direction_ll, - normal_direction_average, - equations)) == RealT @test eltype(@inferred flux_fjordholm_etal(u_ll, u_rr, normal_direction, equations)) == RealT @test eltype(@inferred flux_wintermeyer_etal(u_ll, u_rr, normal_direction, @@ -2181,10 +2175,6 @@ isdir(outdir) && rm(outdir, recursive = true) orientation, equations)) == RealT - @test eltype(eltype(@inferred flux_nonconservative_ersing_etal(u_ll, u_rr, - orientation, - equations))) == - RealT @test eltype(@inferred flux_fjordholm_etal(u_ll, u_rr, orientation, equations)) == RealT @test eltype(@inferred flux_wintermeyer_etal(u_ll, u_rr, orientation, diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index 5c228d1e04c..563f99792d0 100644 --- a/test/test_unstructured_2d.jl +++ b/test/test_unstructured_2d.jl @@ -351,16 +351,16 @@ end @trixi_testset "elixir_shallowwater_well_balanced.jl with FluxHydrostaticReconstruction" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), l2=[ - 1.2164292510839085, - 1.2643106818778908e-12, - 1.269230436589819e-12, - 1.2164292510839079, + 1.2164292510839063, + 1.2676379081600215e-12, + 1.255855785593831e-12, + 1.2164292510839074, ], linf=[ - 1.513851228231562, - 1.6670644673575802e-11, - 1.8426585188623954e-11, - 1.513851228231574, + 1.5138512282315604, + 1.658245722058109e-11, + 1.8665562182185795e-11, + 1.5138512282315737, ], surface_flux=(FluxHydrostaticReconstruction(flux_lax_friedrichs, hydrostatic_reconstruction_audusse_etal), @@ -376,7 +376,7 @@ end end end -@trixi_testset "elixir_shallowwater_well_balanced.jl with flux_nonconservative_ersing_etal" begin +@trixi_testset "elixir_shallowwater_well_balanced.jl with flux_nonconservative_wintermeyer_etal" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), l2=[ 1.2164292510839083, @@ -391,9 +391,7 @@ end 1.513851228231574, ], surface_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), - volume_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), + flux_nonconservative_wintermeyer_etal), tspan=(0.0, 0.25)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -458,7 +456,7 @@ end end end -@trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_ersing_etal" begin +@trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_wintermeyer_etal" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), l2=[ 0.001118046975499805, @@ -473,9 +471,7 @@ end 2.6407324614341476e-5, ], surface_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), - volume_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), + flux_nonconservative_wintermeyer_etal), tspan=(0.0, 0.025)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -517,12 +513,16 @@ end @trixi_testset "elixir_shallowwater_dirichlet.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_dirichlet.jl"), l2=[ - 1.1577518608938916e-5, 4.859252379740366e-13, - 4.639600837197925e-13, 1.1577518608952174e-5, + 1.1577518608950964e-5, + 4.761947272222427e-13, + 4.546045873135486e-13, + 1.157751860893347e-5, ], linf=[ - 8.3940638787805e-5, 1.1446362498574484e-10, - 1.1124515748367981e-10, 8.39406387962427e-5, + 8.394063879002545e-5, + 1.1211566736150389e-10, + 1.0890426250906834e-10, + 8.394063879602065e-5, ], tspan=(0.0, 2.0)) # Ensure that we do not have excessive memory allocations From ef1100d4f1152fa48d674b875e529761eac8966e Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Mon, 19 Aug 2024 10:48:40 +0200 Subject: [PATCH 3/6] set version to v0.8.7 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index d5627651d6d..7800b904c89 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.8.7-DEV" +version = "0.8.7" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" From 26ab5fa1230fb2ef761bb6bb6553f99d9c7e99d3 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Mon, 19 Aug 2024 10:48:58 +0200 Subject: [PATCH 4/6] set development version to v0.8.8-DEV --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 7800b904c89..923cf5248a8 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.8.7" +version = "0.8.8-DEV" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" From ba82d36cc7f7aa9d78a55e92cb94cd5215923d06 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Mon, 19 Aug 2024 13:40:08 +0200 Subject: [PATCH 5/6] Conservative AMR (#2028) * small typo fixes * Interpolate and project Ju instead of u to retain discrete conservation * add specialized routines for TreeMesh and P4estMesh. Hopefully some tests pass now * remove unused copy * remove false comment * simplify the method on the children * update 2d p4est tests * update 3d p4est tests * update 2d t8code tests * update 3d t8code tests * update mpi test values * fix formatting * remove unnecessary comment * add new 2d tests including conservation * add new 3d tests including conservation * update t8code 2d mpi test * remove unnecessary Union * add comments regarding how element ids are incremented during refine or coarsen * typo fixes * another typo fix * fix comment * Apply suggestions from code review Co-authored-by: Daniel Doehring * Apply suggestions from code review Co-authored-by: Hendrik Ranocha * combine GC.preserve blocks * fix formatting in elixir * simplify implementation and add if blocks for P4estMesh specialities * unify structure across P4estMesh and T8codeMesh * fix some of the bad formatting * Apply suggestions from code review Co-authored-by: Hendrik Ranocha * add news * Update NEWS.md Co-authored-by: Andrew Winters --------- Co-authored-by: Daniel Doehring Co-authored-by: Hendrik Ranocha --- NEWS.md | 8 + .../elixir_euler_weak_blast_wave_amr.jl | 117 +++++++++++++ .../elixir_euler_weak_blast_wave_amr.jl | 114 ++++++++++++ .../elixir_advection_amr_unstructured_flag.jl | 3 +- .../elixir_euler_weak_blast_wave_amr.jl | 114 ++++++++++++ .../elixir_euler_weak_blast_wave_amr.jl | 116 +++++++++++++ src/callbacks_step/amr_dg2d.jl | 152 ++++++++++++++-- src/callbacks_step/amr_dg3d.jl | 163 +++++++++++++++--- src/meshes/structured_mesh.jl | 4 +- src/meshes/surface_interpolant.jl | 2 +- test/test_mpi_p4est_2d.jl | 4 +- test/test_mpi_p4est_3d.jl | 4 +- test/test_mpi_t8code_2d.jl | 5 +- test/test_mpi_t8code_3d.jl | 4 +- test/test_p4est_2d.jl | 127 +++++++++----- test/test_p4est_3d.jl | 61 +++++-- test/test_t8code_2d.jl | 52 +++++- test/test_t8code_3d.jl | 59 +++++-- 18 files changed, 980 insertions(+), 129 deletions(-) create mode 100644 examples/p4est_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl create mode 100644 examples/p4est_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl create mode 100644 examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl create mode 100644 examples/t8code_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl diff --git a/NEWS.md b/NEWS.md index bca7d03e954..9371cafa07f 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,6 +4,14 @@ Trixi.jl follows the interpretation of [semantic versioning (semver)](https://ju used in the Julia ecosystem. Notable changes will be documented in this file for human readability. +## Changes in the v0.8 lifecycle + +#### Changed + +- The AMR routines for `P4estMesh` and `T8codeMesh` were changed to work on the product + of the Jacobian and the conserved variables instead of the conserved variables only + to make AMR fully conservative ([#2028]). This may change AMR results slightly. + ## Changes when updating to v0.8 from v0.7.x #### Added diff --git a/examples/p4est_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl b/examples/p4est_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl new file mode 100644 index 00000000000..68680673712 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl @@ -0,0 +1,117 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations2D(1.4) + +function initial_condition_weak_blast_wave(x, t, equations::CompressibleEulerEquations2D) + # Set up polar coordinates + inicenter = SVector(0.0, 0.0) + x_norm = x[1] - inicenter[1] + y_norm = x[2] - inicenter[2] + r = sqrt(x_norm^2 + y_norm^2) + + r0 = 0.2 + E = 1 + p0_inner = 3 + p0_outer = 1 + + # Calculate primitive variables + rho = 1.1 + v1 = 0.0 + v2 = 0.0 + p = r > r0 ? p0_outer : p0_inner + + return prim2cons(SVector(rho, v1, v2, p), equations) +end + +initial_condition = initial_condition_weak_blast_wave + +# Get the DG approximation space + +# Activate the shock capturing + flux differencing +surface_flux = flux_lax_friedrichs +volume_flux = flux_ranocha +polydeg = 4 +basis = LobattoLegendreBasis(polydeg) +indicator_sc = IndicatorHennemannGassner(equations, basis, + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + variable = density_pressure) +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) + +############################################################################### + +# Affine type mapping to take the [-1,1]^2 domain +# and warp it as described in https://arxiv.org/abs/2012.12040 +# Warping with the coefficient 0.2 is even more extreme. +function mapping_twist(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.0 * pi * y) + return SVector(x, y) +end + +# The mesh below can be made periodic +# Create P4estMesh with 8 x 8 trees +trees_per_dimension = (8, 8) +mesh = P4estMesh(trees_per_dimension, polydeg = 4, + mapping = mapping_twist, + initial_refinement_level = 0, + periodicity = true) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.2) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 400 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + save_analysis = true, + extra_analysis_errors = (:conservation_error,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(dt = 0.2, + save_initial_solution = true, + save_final_solution = true) + +amr_indicator = IndicatorLöhner(semi, variable = Trixi.density) +amr_controller = ControllerThreeLevel(semi, amr_indicator, + base_level = 0, + med_level = 1, med_threshold = 0.05, + max_level = 2, max_threshold = 0.1) +amr_callback = AMRCallback(semi, amr_controller, + interval = 5, + adapt_initial_condition = true, + adapt_initial_condition_only_refine = true) + +stepsize_callback = StepsizeCallback(cfl = 0.5) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + save_solution, + amr_callback, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks);#, maxiters=4); +summary_callback() # print the timer summary diff --git a/examples/p4est_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl b/examples/p4est_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl new file mode 100644 index 00000000000..b5b56220004 --- /dev/null +++ b/examples/p4est_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl @@ -0,0 +1,114 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations3D(1.4) + +function initial_condition_weak_blast_wave(x, t, + equations::CompressibleEulerEquations3D) + # Set up polar coordinates + inicenter = SVector(0.0, 0.0, 0.0) + x_norm = x[1] - inicenter[1] + y_norm = x[2] - inicenter[2] + z_norm = x[3] - inicenter[3] + r = sqrt(x_norm^2 + y_norm^2 + z_norm^2) + + r0 = 0.2 + E = 1.0 + p0_inner = 3 + p0_outer = 1 + + # Calculate primitive variables + rho = 1.1 + v1 = 0.0 + v2 = 0.0 + v3 = 0.0 + p = r > r0 ? p0_outer : p0_inner + + return prim2cons(SVector(rho, v1, v2, v3, p), equations) +end + +initial_condition = initial_condition_weak_blast_wave + +surface_flux = flux_lax_friedrichs +volume_flux = flux_ranocha +polydeg = 4 +basis = LobattoLegendreBasis(polydeg) +indicator_sc = IndicatorHennemannGassner(equations, basis, + alpha_max = 1.0, + alpha_min = 0.001, + alpha_smooth = true, + variable = density_pressure) +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) + +# Setup a periodic mesh with 4 x 4 x 4 trees and 8 x 8 x 8 elements +trees_per_dimension = (4, 4, 4) + +# Affine type mapping to take the [-1,1]^3 domain +# and warp it as described in https://arxiv.org/abs/2012.12040 +function mapping_twist(xi, eta, zeta) + y = eta + 1 / 6 * (cos(1.5 * pi * xi) * cos(0.5 * pi * eta) * cos(0.5 * pi * zeta)) + + x = xi + 1 / 6 * (cos(0.5 * pi * xi) * cos(2 * pi * y) * cos(0.5 * pi * zeta)) + + z = zeta + 1 / 6 * (cos(0.5 * pi * x) * cos(pi * y) * cos(0.5 * pi * zeta)) + + return SVector(x, y, z) +end + +mesh = P4estMesh(trees_per_dimension, + polydeg = 2, + mapping = mapping_twist, + initial_refinement_level = 1, + periodicity = true) + +# Create the semidiscretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_errors = (:conservation_error,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +amr_indicator = IndicatorLöhner(semi, variable = Trixi.density) +amr_controller = ControllerThreeLevel(semi, amr_indicator, + base_level = 1, + med_level = 2, med_threshold = 0.05, + max_level = 3, max_threshold = 0.15) +amr_callback = AMRCallback(semi, amr_controller, + interval = 1, + adapt_initial_condition = false, + adapt_initial_condition_only_refine = false) + +stepsize_callback = StepsizeCallback(cfl = 0.5) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + amr_callback, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl b/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl index 9138586cccf..d8893854811 100644 --- a/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl +++ b/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl @@ -61,7 +61,8 @@ amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first) amr_callback = AMRCallback(semi, amr_controller, interval = 5, adapt_initial_condition = true, - adapt_initial_condition_only_refine = true) + adapt_initial_condition_only_refine = true, + dynamic_load_balancing = true) stepsize_callback = StepsizeCallback(cfl = 0.7) diff --git a/examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl b/examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl new file mode 100644 index 00000000000..a3366caa317 --- /dev/null +++ b/examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl @@ -0,0 +1,114 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations2D(1.4) + +function initial_condition_weak_blast_wave(x, t, equations::CompressibleEulerEquations2D) + # Set up polar coordinates + inicenter = SVector(0.0, 0.0) + x_norm = x[1] - inicenter[1] + y_norm = x[2] - inicenter[2] + r = sqrt(x_norm^2 + y_norm^2) + + r0 = 0.2 + E = 1 + p0_inner = 3 + p0_outer = 1 + + # Calculate primitive variables + rho = 1.1 + v1 = 0.0 + v2 = 0.0 + p = r > r0 ? p0_outer : p0_inner + + return prim2cons(SVector(rho, v1, v2, p), equations) +end + +initial_condition = initial_condition_weak_blast_wave + +# Get the DG approximation space + +# Activate the shock capturing + flux differencing +surface_flux = flux_lax_friedrichs +volume_flux = flux_ranocha +polydeg = 4 +basis = LobattoLegendreBasis(polydeg) +indicator_sc = IndicatorHennemannGassner(equations, basis, + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + variable = density_pressure) +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) + +# Affine type mapping to take the [-1,1]^2 domain +# and warp it as described in https://arxiv.org/abs/2012.12040 +# Warping with the coefficient 0.2 is even more extreme. +function mapping_twist(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.0 * pi * y) + return SVector(x, y) +end + +# The mesh below can be made periodic +# Create T8codeMesh with 8 x 8 trees +trees_per_dimension = (8, 8) +mesh = T8codeMesh(trees_per_dimension, polydeg = 4, + mapping = mapping_twist, + initial_refinement_level = 0, + periodicity = true) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.2) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 400 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + save_analysis = true, + extra_analysis_errors = (:conservation_error,)) + +alive_callback = AliveCallback(analysis_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.05, + max_level = 2, max_threshold = 0.1) +amr_callback = AMRCallback(semi, amr_controller, + interval = 5, + adapt_initial_condition = true, + adapt_initial_condition_only_refine = true) + +stepsize_callback = StepsizeCallback(cfl = 0.5) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + amr_callback, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks);#, maxiters=4); +summary_callback() # print the timer summary + +# Finalize `T8codeMesh` to make sure MPI related objects in t8code are +# released before `MPI` finalizes. +!isinteractive() && finalize(mesh) diff --git a/examples/t8code_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl b/examples/t8code_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl new file mode 100644 index 00000000000..106b4821144 --- /dev/null +++ b/examples/t8code_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl @@ -0,0 +1,116 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations3D(1.4) + +function initial_condition_weak_blast_wave(x, t, + equations::CompressibleEulerEquations3D) + # Set up polar coordinates + inicenter = SVector(0.0, 0.0, 0.0) + x_norm = x[1] - inicenter[1] + y_norm = x[2] - inicenter[2] + z_norm = x[3] - inicenter[3] + r = sqrt(x_norm^2 + y_norm^2 + z_norm^2) + + r0 = 0.2 + E = 1.0 + p0_inner = 3 + p0_outer = 1 + + # Calculate primitive variables + rho = 1.1 + v1 = 0.0 + v2 = 0.0 + v3 = 0.0 + p = r > r0 ? p0_outer : p0_inner + + return prim2cons(SVector(rho, v1, v2, v3, p), equations) +end + +initial_condition = initial_condition_weak_blast_wave + +surface_flux = flux_lax_friedrichs +volume_flux = flux_ranocha +polydeg = 4 +basis = LobattoLegendreBasis(polydeg) +indicator_sc = IndicatorHennemannGassner(equations, basis, + alpha_max = 1.0, + alpha_min = 0.001, + alpha_smooth = true, + variable = density_pressure) +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) + +# Setup a periodic mesh with 4 x 4 x 4 trees and 8 x 8 x 8 elements +trees_per_dimension = (4, 4, 4) + +function mapping_twist(xi, eta, zeta) + y = eta + 1 / 6 * (cos(1.5 * pi * xi) * cos(0.5 * pi * eta) * cos(0.5 * pi * zeta)) + + x = xi + 1 / 6 * (cos(0.5 * pi * xi) * cos(2 * pi * y) * cos(0.5 * pi * zeta)) + + z = zeta + 1 / 6 * (cos(0.5 * pi * x) * cos(pi * y) * cos(0.5 * pi * zeta)) + + return SVector(x, y, z) +end + +mesh = T8codeMesh(trees_per_dimension, + polydeg = 2, + mapping = mapping_twist, + initial_refinement_level = 1, + periodicity = true) + +# Create the semidiscretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_errors = (:conservation_error,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +amr_indicator = IndicatorLöhner(semi, variable = Trixi.density) +amr_controller = ControllerThreeLevel(semi, amr_indicator, + base_level = 1, + med_level = 2, med_threshold = 0.05, + max_level = 3, max_threshold = 0.15) +amr_callback = AMRCallback(semi, amr_controller, + interval = 1, + adapt_initial_condition = false, + adapt_initial_condition_only_refine = false) + +stepsize_callback = StepsizeCallback(cfl = 0.5) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + amr_callback, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary + +# Finalize `T8codeMesh` to make sure MPI related objects in t8code are +# released before `MPI` finalizes. +!isinteractive() && finalize(mesh) diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index 94524b23a3a..062020b6d42 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -76,7 +76,10 @@ function rebalance_solver!(u_ode::AbstractVector, mesh::TreeMesh{2}, equations, end # GC.@preserve old_u_ode end -# Refine elements in the DG solver based on a list of cell_ids that should be refined +# Refine elements in the DG solver based on a list of cell_ids that should be refined. +# On `P4estMesh`, if an element refines the solution scaled by the Jacobian `J*u` is interpolated +# from the parent element into the four children elements. The solution on each child +# element is then recovered by dividing by the new element Jacobians. function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, P4estMesh{2}}, equations, dg::DGSEM, cache, elements_to_refine) # Return early if there is nothing to do @@ -96,9 +99,23 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, P4estM # Retain current solution data old_n_elements = nelements(dg, cache) old_u_ode = copy(u_ode) - GC.@preserve old_u_ode begin # OBS! If we don't GC.@preserve old_u_ode, it might be GC'ed + old_inverse_jacobian = copy(cache.elements.inverse_jacobian) + # OBS! If we don't GC.@preserve old_u_ode and old_inverse_jacobian, they might be GC'ed + GC.@preserve old_u_ode old_inverse_jacobian begin old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) + if mesh isa P4estMesh + # Loop over all elements in old container and scale the old solution by the Jacobian + # prior to projection + for old_element_id in 1:old_n_elements + for v in eachvariable(equations) + old_u[v, .., old_element_id] .= (old_u[v, .., old_element_id] ./ + old_inverse_jacobian[.., + old_element_id]) + end + end + end + reinitialize_containers!(mesh, equations, dg, cache) resize!(u_ode, @@ -112,10 +129,34 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, P4estM # Refine element and store solution directly in new data structure refine_element!(u, element_id, old_u, old_element_id, adaptor, equations, dg) + + if mesh isa P4estMesh + # Before `element_id` is incremented, divide by the new Jacobians on each + # child element and save the result + for m in 0:3 # loop over the children + for v in eachvariable(equations) + u[v, .., element_id + m] .*= (0.25f0 .* + cache.elements.inverse_jacobian[.., + element_id + m]) + end + end + end + + # Increment `element_id` on the refined mesh with the number + # of children, i.e., 4 in 2D element_id += 2^ndims(mesh) else - # Copy old element data to new element container - @views u[:, .., element_id] .= old_u[:, .., old_element_id] + if mesh isa P4estMesh + # Copy old element data to new element container and remove Jacobian scaling + for v in eachvariable(equations) + u[v, .., element_id] .= (old_u[v, .., old_element_id] .* + old_inverse_jacobian[.., + old_element_id]) + end + else # isa TreeMesh + @views u[:, .., element_id] .= old_u[:, .., old_element_id] + end + # No refinement occurred, so increment `element_id` on the new mesh by one element_id += 1 end end @@ -125,7 +166,7 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, P4estM @assert element_id == nelements(dg, cache) + 1||element_id == nelements(dg, cache) + 2^ndims(mesh) "element_id = $element_id, nelements(dg, cache) = $(nelements(dg, cache))" - end # GC.@preserve old_u_ode + end # GC.@preserve old_u_ode old_inverse_jacobian # Sanity check if mesh isa TreeMesh && isperiodic(mesh.tree) && nmortars(cache.mortars) == 0 && @@ -228,7 +269,10 @@ function refine_element!(u::AbstractArray{<:Any, 4}, element_id, return nothing end -# Coarsen elements in the DG solver based on a list of cell_ids that should be removed +# Coarsen elements in the DG solver based on a list of cell_ids that should be removed. +# On `P4estMesh`, if an element coarsens the solution scaled by the Jacobian `J*u` is projected +# from the four children elements back onto the parent element. The solution on the parent +# element is then recovered by dividing by the new element Jacobian. function coarsen!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, P4estMesh{2}}, equations, dg::DGSEM, cache, elements_to_remove) @@ -246,12 +290,26 @@ function coarsen!(u_ode::AbstractVector, adaptor, to_be_removed = falses(nelements(dg, cache)) to_be_removed[elements_to_remove] .= true - # Retain current solution data + # Retain current solution data and Jacobians old_n_elements = nelements(dg, cache) old_u_ode = copy(u_ode) - GC.@preserve old_u_ode begin # OBS! If we don't GC.@preserve old_u_ode, it might be GC'ed + old_inverse_jacobian = copy(cache.elements.inverse_jacobian) + # OBS! If we don't GC.@preserve old_u_ode and old_inverse_jacobian, they might be GC'ed + GC.@preserve old_u_ode old_inverse_jacobian begin old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) + if mesh isa P4estMesh + # Loop over all elements in old container and scale the old solution by the Jacobian + # prior to projection + for old_element_id in 1:old_n_elements + for v in eachvariable(equations) + old_u[v, .., old_element_id] .= (old_u[v, .., old_element_id] ./ + old_inverse_jacobian[.., + old_element_id]) + end + end + end + reinitialize_containers!(mesh, equations, dg, cache) resize!(u_ode, @@ -277,17 +335,39 @@ function coarsen!(u_ode::AbstractVector, adaptor, # Coarsen elements and store solution directly in new data structure coarsen_elements!(u, element_id, old_u, old_element_id, adaptor, equations, dg) + + if mesh isa P4estMesh + # Before `element_id` is incremented, divide by the new Jacobian and save + # the result in the parent element + for v in eachvariable(equations) + u[v, .., element_id] .*= (4 .* + cache.elements.inverse_jacobian[.., + element_id]) + end + end + + # Increment `element_id` on the coarsened mesh by one and `skip` = 3 in 2D + # because 4 children elements become 1 parent element element_id += 1 skip = 2^ndims(mesh) - 1 else - # Copy old element data to new element container - @views u[:, .., element_id] .= old_u[:, .., old_element_id] + if mesh isa P4estMesh + # Copy old element data to new element container and remove Jacobian scaling + for v in eachvariable(equations) + u[v, .., element_id] .= (old_u[v, .., old_element_id] .* + old_inverse_jacobian[.., + old_element_id]) + end + else # isa TreeMesh + @views u[:, .., element_id] .= old_u[:, .., old_element_id] + end + # No coarsening occurred, so increment `element_id` on the new mesh by one element_id += 1 end end # If everything is correct, we should have processed all elements. @assert element_id==nelements(dg, cache) + 1 "element_id = $element_id, nelements(dg, cache) = $(nelements(dg, cache))" - end # GC.@preserve old_u_ode + end # GC.@preserve old_u_ode old_inverse_jacobian # Sanity check if mesh isa TreeMesh && isperiodic(mesh.tree) && nmortars(cache.mortars) == 0 && @@ -404,16 +484,27 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{2}, equations, # Note: This is true for `quads`. T8_CHILDREN = 4 - # Retain current solution data. + # Retain current solution and inverse Jacobian data. old_u_ode = copy(u_ode) + old_inverse_jacobian = copy(cache.elements.inverse_jacobian) + # OBS! If we don't GC.@preserve old_u_ode and old_inverse_jacobian, they might be GC'ed GC.@preserve old_u_ode begin old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) + # Loop over all elements in old container and scale the old solution by the Jacobian + # prior to interpolation or projection + for old_element_id in 1:old_nelems + for v in eachvariable(equations) + old_u[v, .., old_element_id] .= (old_u[v, .., old_element_id] ./ + old_inverse_jacobian[.., + old_element_id]) + end + end + reinitialize_containers!(mesh, equations, dg, cache) - resize!(u_ode, - nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache)) + resize!(u_ode, nvariables(equations) * ndofs(mesh, dg, cache)) u = wrap_array(u_ode, mesh, equations, dg, cache) while old_index <= old_nelems && new_index <= new_nelems @@ -422,6 +513,18 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{2}, equations, # Refine element and store solution directly in new data structure. refine_element!(u, new_index, old_u, old_index, adaptor, equations, dg) + # Before indices are incremented divide by the new Jacobians on each + # child element and save the result + for m in 0:3 # loop over the children + for v in eachvariable(equations) + u[v, .., new_index + m] .*= (0.25f0 .* + cache.elements.inverse_jacobian[.., + new_index + m]) + end + end + + # Increment `old_index` on the original mesh and the `new_index` + # on the refined mesh with the number of children, i.e., T8_CHILDREN = 4 old_index += 1 new_index += T8_CHILDREN @@ -436,19 +539,34 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{2}, equations, coarsen_elements!(u, new_index, old_u, old_index, adaptor, equations, dg) + # Before the indices are incremented divide by the new Jacobian and save + # the result again in the parent element + for v in eachvariable(equations) + u[v, .., new_index] .*= (4 .* cache.elements.inverse_jacobian[.., + new_index]) + end + + # Increment `old_index` on the original mesh with the number of children + # (T8_CHILDREN = 4 in 2D) and the `new_index` by one for the single + # coarsened element old_index += T8_CHILDREN new_index += 1 else # No changes. - # Copy old element data to new element container. - @views u[:, .., new_index] .= old_u[:, .., old_index] + # Copy old element data to new element container and remove Jacobian scaling + for v in eachvariable(equations) + u[v, .., new_index] .= (old_u[v, .., old_index] .* + old_inverse_jacobian[.., old_index]) + end + # No refinement / coarsening occurred, so increment element index + # on each mesh by one old_index += 1 new_index += 1 end end # while - end # GC.@preserve old_u_ode + end # GC.@preserve old_u_ode old_inverse_jacobian return nothing end diff --git a/src/callbacks_step/amr_dg3d.jl b/src/callbacks_step/amr_dg3d.jl index 3f67951bafe..594c30dcca5 100644 --- a/src/callbacks_step/amr_dg3d.jl +++ b/src/callbacks_step/amr_dg3d.jl @@ -5,9 +5,11 @@ @muladd begin #! format: noindent -# Refine elements in the DG solver based on a list of cell_ids that should be refined -function refine!(u_ode::AbstractVector, adaptor, - mesh::Union{TreeMesh{3}, P4estMesh{3}}, +# Refine elements in the DG solver based on a list of cell_ids that should be refined. +# On `P4estMesh`, if an element refines the solution scaled by the Jacobian `J*u` is interpolated +# from the parent element into the eight children elements. The solution on each child +# element is then recovered by dividing by the new element Jacobians. +function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{3}, P4estMesh{3}}, equations, dg::DGSEM, cache, elements_to_refine) # Return early if there is nothing to do if isempty(elements_to_refine) @@ -26,9 +28,23 @@ function refine!(u_ode::AbstractVector, adaptor, # Retain current solution data old_n_elements = nelements(dg, cache) old_u_ode = copy(u_ode) - GC.@preserve old_u_ode begin # OBS! If we don't GC.@preserve old_u_ode, it might be GC'ed + old_inverse_jacobian = copy(cache.elements.inverse_jacobian) + # OBS! If we don't GC.@preserve old_u_ode and old_inverse_jacobian, they might be GC'ed + GC.@preserve old_u_ode old_inverse_jacobian begin old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) + if mesh isa P4estMesh + # Loop over all elements in old container and scale the old solution by the Jacobian + # prior to projection + for old_element_id in 1:old_n_elements + for v in eachvariable(equations) + old_u[v, .., old_element_id] .= (old_u[v, .., old_element_id] ./ + old_inverse_jacobian[.., + old_element_id]) + end + end + end + reinitialize_containers!(mesh, equations, dg, cache) resize!(u_ode, @@ -44,12 +60,36 @@ function refine!(u_ode::AbstractVector, adaptor, for old_element_id in 1:old_n_elements if needs_refinement[old_element_id] # Refine element and store solution directly in new data structure - refine_element!(u, element_id, old_u, old_element_id, - adaptor, equations, dg, u_tmp1, u_tmp2) + refine_element!(u, element_id, old_u, old_element_id, adaptor, + equations, dg, u_tmp1, u_tmp2) + + if mesh isa P4estMesh + # Before `element_id` is incremented, divide by the new Jacobians on each + # child element and save the result + for m in 0:7 # loop over the children + for v in eachvariable(equations) + u[v, .., element_id + m] .*= (0.125f0 .* + cache.elements.inverse_jacobian[.., + element_id + m]) + end + end + end + + # Increment `element_id` on the refined mesh with the number + # of children, i.e., 8 in 3D element_id += 2^ndims(mesh) else - # Copy old element data to new element container - @views u[:, .., element_id] .= old_u[:, .., old_element_id] + if mesh isa P4estMesh + # Copy old element data to new element container and remove Jacobian scaling + for v in eachvariable(equations) + u[v, .., element_id] .= (old_u[v, .., old_element_id] .* + old_inverse_jacobian[.., + old_element_id]) + end + else # isa TreeMesh + @views u[:, .., element_id] .= old_u[:, .., old_element_id] + end + # No refinement occurred, so increment `element_id` on the new mesh by one element_id += 1 end end @@ -59,7 +99,7 @@ function refine!(u_ode::AbstractVector, adaptor, @assert element_id == nelements(dg, cache) + 1||element_id == nelements(dg, cache) + 2^ndims(mesh) "element_id = $element_id, nelements(dg, cache) = $(nelements(dg, cache))" - end # GC.@preserve old_u_ode + end # GC.@preserve old_u_ode old_inverse_jacobian # Sanity check if mesh isa TreeMesh && isperiodic(mesh.tree) && nmortars(cache.mortars) == 0 @@ -145,7 +185,10 @@ function refine_element!(u::AbstractArray{<:Any, 5}, element_id, return nothing end -# Coarsen elements in the DG solver based on a list of cell_ids that should be removed +# Coarsen elements in the DG solver based on a list of cell_ids that should be removed. +# On `P4estMesh`, if an element coarsens the solution scaled by the Jacobian `J*u` is projected +# from the eight children elements back onto the parent element. The solution on the parent +# element is then recovered by dividing by the new element Jacobian. function coarsen!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{3}, P4estMesh{3}}, equations, dg::DGSEM, cache, elements_to_remove) @@ -163,12 +206,26 @@ function coarsen!(u_ode::AbstractVector, adaptor, to_be_removed = falses(nelements(dg, cache)) to_be_removed[elements_to_remove] .= true - # Retain current solution data + # Retain current solution data and Jacobians old_n_elements = nelements(dg, cache) old_u_ode = copy(u_ode) - GC.@preserve old_u_ode begin # OBS! If we don't GC.@preserve old_u_ode, it might be GC'ed + old_inverse_jacobian = copy(cache.elements.inverse_jacobian) + # OBS! If we don't GC.@preserve old_u_ode and old_inverse_jacobian, they might be GC'ed + GC.@preserve old_u_ode old_inverse_jacobian begin old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) + if mesh isa P4estMesh + # Loop over all elements in old container and scale the old solution by the Jacobian + # prior to projection + for old_element_id in 1:old_n_elements + for v in eachvariable(equations) + old_u[v, .., old_element_id] .= (old_u[v, .., old_element_id] ./ + old_inverse_jacobian[.., + old_element_id]) + end + end + end + reinitialize_containers!(mesh, equations, dg, cache) resize!(u_ode, @@ -196,19 +253,41 @@ function coarsen!(u_ode::AbstractVector, adaptor, @assert all(to_be_removed[old_element_id:(old_element_id + 2^ndims(mesh) - 1)]) "bad cell/element order" # Coarsen elements and store solution directly in new data structure - coarsen_elements!(u, element_id, old_u, old_element_id, - adaptor, equations, dg, u_tmp1, u_tmp2) + coarsen_elements!(u, element_id, old_u, old_element_id, adaptor, + equations, dg, u_tmp1, u_tmp2) + + if mesh isa P4estMesh + # Before `element_id` is incremented, divide by the new Jacobian and save + # the result in the parent element + for v in eachvariable(equations) + u[v, .., element_id] .*= (8 .* + cache.elements.inverse_jacobian[.., + element_id]) + end + end + + # Increment `element_id` on the coarsened mesh by one and `skip` = 7 in 3D + # because 8 children elements become 1 parent element element_id += 1 skip = 2^ndims(mesh) - 1 else - # Copy old element data to new element container - @views u[:, .., element_id] .= old_u[:, .., old_element_id] + if mesh isa P4estMesh + # Copy old element data to new element container and remove Jacobian scaling + for v in eachvariable(equations) + u[v, .., element_id] .= (old_u[v, .., old_element_id] .* + old_inverse_jacobian[.., + old_element_id]) + end + else # isa TreeMesh + @views u[:, .., element_id] .= old_u[:, .., old_element_id] + end + # No coarsening occurred, so increment `element_id` on the new mesh by one element_id += 1 end end # If everything is correct, we should have processed all elements. @assert element_id==nelements(dg, cache) + 1 "element_id = $element_id, nelements(dg, cache) = $(nelements(dg, cache))" - end # GC.@preserve old_u_ode + end # GC.@preserve old_u_ode old_inverse_jacobian # Sanity check if mesh isa TreeMesh && isperiodic(mesh.tree) && nmortars(cache.mortars) == 0 @@ -335,16 +414,27 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{3}, equations, # Note: This is only true for `hexs`. T8_CHILDREN = 8 - # Retain current solution data. + # Retain current solution and inverse Jacobian data. old_u_ode = copy(u_ode) + old_inverse_jacobian = copy(cache.elements.inverse_jacobian) + # OBS! If we don't GC.@preserve old_u_ode and old_inverse_jacobian, they might be GC'ed GC.@preserve old_u_ode begin old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) + # Loop over all elements in old container and scale the old solution by the Jacobian + # prior to interpolation or projection + for old_element_id in 1:old_nelems + for v in eachvariable(equations) + old_u[v, .., old_element_id] .= (old_u[v, .., old_element_id] ./ + old_inverse_jacobian[.., + old_element_id]) + end + end + reinitialize_containers!(mesh, equations, dg, cache) - resize!(u_ode, - nvariables(equations) * ndofs(mesh, dg, cache)) + resize!(u_ode, nvariables(equations) * ndofs(mesh, dg, cache)) u = wrap_array(u_ode, mesh, equations, dg, cache) u_tmp1 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), @@ -359,6 +449,18 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{3}, equations, refine_element!(u, new_index, old_u, old_index, adaptor, equations, dg, u_tmp1, u_tmp2) + # Before indices are incremented divide by the new Jacobians on each + # child element and save the result + for m in 0:7 # loop over the children + for v in eachvariable(equations) + u[v, .., new_index + m] .*= (0.125f0 .* + cache.elements.inverse_jacobian[.., + new_index + m]) + end + end + + # Increment `old_index` on the original mesh and the `new_index` + # on the refined mesh with the number of children, i.e., T8_CHILDREN = 8 old_index += 1 new_index += T8_CHILDREN @@ -373,19 +475,34 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{3}, equations, coarsen_elements!(u, new_index, old_u, old_index, adaptor, equations, dg, u_tmp1, u_tmp2) + # Before the indices are incremented divide by the new Jacobian and save + # the result again in the parent element + for v in eachvariable(equations) + u[v, .., new_index] .*= (8 .* cache.elements.inverse_jacobian[.., + new_index]) + end + + # Increment `old_index` on the original mesh with the number of children + # (T8_CHILDREN = 8 in 3D) and the `new_index` by one for the single + # coarsened element old_index += T8_CHILDREN new_index += 1 else # No changes. - # Copy old element data to new element container. - @views u[:, .., new_index] .= old_u[:, .., old_index] + # Copy old element data to new element container and remove Jacobian scaling + for v in eachvariable(equations) + u[v, .., new_index] .= (old_u[v, .., old_index] .* + old_inverse_jacobian[.., old_index]) + end + # No refinement / coarsening occurred, so increment element index + # on each mesh by one old_index += 1 new_index += 1 end end # while - end # GC.@preserve old_u_ode + end # GC.@preserve old_u_ode old_inverse_jacobian return nothing end diff --git a/src/meshes/structured_mesh.jl b/src/meshes/structured_mesh.jl index 553aabbbc20..9e5c55197d8 100644 --- a/src/meshes/structured_mesh.jl +++ b/src/meshes/structured_mesh.jl @@ -255,13 +255,15 @@ function correction_term_3d(x, y, z, faces) linear_interpolate(x, faces[1](1, z), faces[2](1, z)) + linear_interpolate(z, faces[5](x, 1), faces[6](x, 1))) - # Correction for x-terms + # Correction for z-terms c_z = linear_interpolate(z, linear_interpolate(x, faces[1](y, -1), faces[2](y, -1)) + linear_interpolate(y, faces[3](x, -1), faces[4](x, -1)), linear_interpolate(x, faces[1](y, 1), faces[2](y, 1)) + linear_interpolate(y, faces[3](x, 1), faces[4](x, 1))) + # Each of the 12 edges are counted twice above + # so we divide the correction term by two return 0.5 * (c_x + c_y + c_z) end diff --git a/src/meshes/surface_interpolant.jl b/src/meshes/surface_interpolant.jl index 22d14e38c5c..6a7c5c7834b 100644 --- a/src/meshes/surface_interpolant.jl +++ b/src/meshes/surface_interpolant.jl @@ -52,7 +52,7 @@ function derivative_at(s, boundary_curve::CurvedSurface) y_coordinate_at_s_on_boundary_curve_prime end -# Chebysehv-Gauss-Lobatto nodes and weights for use with curved boundaries +# Chebyshev-Gauss-Lobatto nodes and weights for use with curved boundaries function chebyshev_gauss_lobatto_nodes_weights(n_nodes::Integer) # Initialize output diff --git a/test/test_mpi_p4est_2d.jl b/test/test_mpi_p4est_2d.jl index 6d66bc68a26..29de4efc6d0 100644 --- a/test/test_mpi_p4est_2d.jl +++ b/test/test_mpi_p4est_2d.jl @@ -97,8 +97,8 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "p4est_2d_dgsem") @trixi_testset "elixir_advection_amr_unstructured_flag.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr_unstructured_flag.jl"), - l2=[0.0012766060609964525], - linf=[0.01750280631586159], + l2=[0.0012808538770535593], + linf=[0.01752690016659812], coverage_override=(maxiters = 6,)) # Ensure that we do not have excessive memory allocations diff --git a/test/test_mpi_p4est_3d.jl b/test/test_mpi_p4est_3d.jl index cca9093ec51..4f9465b85dc 100644 --- a/test/test_mpi_p4est_3d.jl +++ b/test/test_mpi_p4est_3d.jl @@ -69,8 +69,8 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "p4est_3d_dgsem") @trixi_testset "elixir_advection_amr_unstructured_curved.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr_unstructured_curved.jl"), - l2=[1.6236411810065552e-5], - linf=[0.0010554006923731395], + l2=[1.6163120948209677e-5], + linf=[0.0010572201890564834], tspan=(0.0, 1.0), coverage_override=(maxiters = 6, initial_refinement_level = 0, diff --git a/test/test_mpi_t8code_2d.jl b/test/test_mpi_t8code_2d.jl index 7c7fc03898c..75e65c8c380 100644 --- a/test/test_mpi_t8code_2d.jl +++ b/test/test_mpi_t8code_2d.jl @@ -97,8 +97,9 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "t8code_2d_dgsem") @trixi_testset "elixir_advection_amr_unstructured_flag.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr_unstructured_flag.jl"), - l2=[0.001980652042312077], - linf=[0.0328882442132265], + l2=[0.002019623611753929], + linf=[0.03542375961299987], + dynamic_load_balancing=false, coverage_override=(maxiters = 6,)) # Ensure that we do not have excessive memory allocations diff --git a/test/test_mpi_t8code_3d.jl b/test/test_mpi_t8code_3d.jl index a15690a7629..2e837f79ad8 100644 --- a/test/test_mpi_t8code_3d.jl +++ b/test/test_mpi_t8code_3d.jl @@ -69,8 +69,8 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "t8code_3d_dgsem") @trixi_testset "elixir_advection_amr_unstructured_curved.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr_unstructured_curved.jl"), - l2=[2.0556575425846923e-5], - linf=[0.00105682693484822], + l2=[2.0535121347526814e-5], + linf=[0.0010586603797777504], tspan=(0.0, 1.0), coverage_override=(maxiters = 6, initial_refinement_level = 0, diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index ef5ac2c9de3..f56dc9cd769 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -78,8 +78,8 @@ end @trixi_testset "elixir_advection_amr_unstructured_flag.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr_unstructured_flag.jl"), - l2=[0.0012766060609964525], - linf=[0.01750280631586159], + l2=[0.0012808538770535593], + linf=[0.01752690016659812], coverage_override=(maxiters = 6,)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -301,16 +301,16 @@ end @trixi_testset "elixir_euler_blast_wave_amr.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_amr.jl"), l2=[ - 6.32183914e-01, - 3.86914231e-01, - 3.86869171e-01, - 1.06575688e+00, + 0.6321850210104147, + 0.38691446170269167, + 0.3868695626809587, + 1.0657553825683956, ], linf=[ - 2.76020890e+00, - 2.32659890e+00, - 2.32580837e+00, - 2.15778188e+00, + 2.7602280007469666, + 2.3265993814913672, + 2.3258078438689673, + 2.1577683028925416, ], tspan=(0.0, 0.3), coverage_override=(maxiters = 6,)) @@ -327,16 +327,16 @@ end @trixi_testset "elixir_euler_wall_bc_amr.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_wall_bc_amr.jl"), l2=[ - 0.020291447969983396, - 0.017479614254319948, - 0.011387644425613437, - 0.0514420126021293, + 0.02026685991647352, + 0.017467584076280237, + 0.011378371604813321, + 0.05138942558296091, ], linf=[ - 0.3582779022370579, - 0.32073537890751663, - 0.221818049107692, - 0.9209559420400415, + 0.35924402060711524, + 0.32068389566068806, + 0.2361141752119986, + 0.9289840057748628, ], tspan=(0.0, 0.15)) # Ensure that we do not have excessive memory allocations @@ -352,16 +352,16 @@ end @trixi_testset "elixir_euler_forward_step_amr.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_forward_step_amr.jl"), l2=[ - 0.004194875320833303, - 0.003785140699353966, - 0.0013696609105790351, - 0.03265268616046424, + 0.004191480950848891, + 0.003781298410569231, + 0.0013470418422981045, + 0.03262817609394949, ], linf=[ - 2.0585399781442852, - 2.213428805506876, - 3.862362410419163, - 17.75187237459251, + 2.0581500751947113, + 2.2051301367971288, + 3.8502467979250254, + 17.750333649853616, ], tspan=(0.0, 0.0001), rtol=1.0e-7, @@ -409,16 +409,16 @@ end @trixi_testset "elixir_euler_supersonic_cylinder.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_supersonic_cylinder.jl"), l2=[ - 0.026798021911954406, - 0.05118546368109259, - 0.03206703583774831, - 0.19680026567208672, + 0.02676082999794676, + 0.05110830068968181, + 0.03205164257040607, + 0.1965981012724311, ], linf=[ - 3.653905721692421, - 4.285035711361009, - 6.8544353186357645, - 31.748244912257533, + 3.6830683476364476, + 4.284442685012427, + 6.857777546171545, + 31.749285097390576, ], tspan=(0.0, 0.001), skip_coverage=true) @@ -529,16 +529,17 @@ end @trixi_testset "elixir_mhd_rotor.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_rotor.jl"), - l2=[0.4552084651735862, 0.8918048264575757, 0.832471223081887, + l2=[0.4552094211937344, 0.8918052934760807, 0.8324715234680768, 0.0, - 0.9801264164951583, 0.10475690769435382, 0.1555132409149897, + 0.9801268321975978, 0.10475722739111007, + 0.15551326369033164, 0.0, - 2.0597079362763556e-5], - linf=[10.194181233788775, 18.25472397868819, 10.031307436191334, + 2.0602990858239632e-5], + linf=[10.19421969147307, 18.254409357804683, 10.031954811332596, 0.0, - 19.647239392277378, 1.3938810140985936, 1.8724965294853084, + 19.646870938371492, 1.3938679692894465, 1.8725058401937984, 0.0, - 0.0016290067532561904], + 0.0016201762010257296], tspan=(0.0, 0.02)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -617,12 +618,12 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_NACA0012airfoil_mach085.jl"), l2=[ - 5.634402680811982e-7, 6.748066107517321e-6, - 1.091879472416885e-5, 0.0006686372064029146, + 5.56114097044427e-7, 6.62284247153255e-6, + 1.0823259724601275e-5, 0.000659804574787503, ], linf=[ - 0.0021456247890772823, 0.03957142889488085, - 0.03832024233032798, 2.6628739573358495, + 0.002157589754528455, 0.039163189253511164, + 0.038386804399707625, 2.6685831417913914, ], amr_interval=1, base_level=0, med_level=1, max_level=1, @@ -652,8 +653,8 @@ end lift = Trixi.analyze(lift_coefficient, du, u, tspan[2], mesh, equations, solver, semi.cache, semi) - @test isapprox(lift, 0.029076443678087403, atol = 1e-13) - @test isapprox(drag, 0.13564720009197903, atol = 1e-13) + @test isapprox(lift, 0.029094009322876882, atol = 1e-13) + @test isapprox(drag, 0.13579200776643238, atol = 1e-13) end @trixi_testset "elixir_euler_blast_wave_pure_fv.jl" begin @@ -684,6 +685,40 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + +@trixi_testset "elixir_euler_weak_blast_wave_amr.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weak_blast_wave_amr.jl"), + l2=[ + 0.11134260363848127, + 0.11752357091804219, + 0.11829112104640764, + 0.7557891142955036, + ], + linf=[ + 0.5728647031475109, + 0.8353132977670252, + 0.8266797080712205, + 3.9792506230548317, + ], + tspan=(0.0, 0.1), + coverage_override=(maxiters = 6,)) + # 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 + # Check for conservation + state_integrals = Trixi.integrate(sol.u[2], semi) + initial_state_integrals = analysis_callback.affect!.initial_state_integrals + + @test isapprox(state_integrals[1], initial_state_integrals[1], atol = 1e-13) + @test isapprox(state_integrals[2], initial_state_integrals[2], atol = 1e-13) + @test isapprox(state_integrals[3], initial_state_integrals[3], atol = 1e-13) + @test isapprox(state_integrals[4], initial_state_integrals[4], atol = 1e-13) +end end # Clean up afterwards: delete Trixi.jl output directory diff --git a/test/test_p4est_3d.jl b/test/test_p4est_3d.jl index 7483cde2752..7b69d1c0cf2 100644 --- a/test/test_p4est_3d.jl +++ b/test/test_p4est_3d.jl @@ -78,8 +78,8 @@ end @trixi_testset "elixir_advection_amr_unstructured_curved.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr_unstructured_curved.jl"), - l2=[1.6236411810065552e-5], - linf=[0.0010554006923731395], + l2=[1.6163120948209677e-5], + linf=[0.0010572201890564834], tspan=(0.0, 1.0), coverage_override=(maxiters = 6, initial_refinement_level = 0, base_level = 0, med_level = 1, max_level = 2)) @@ -542,16 +542,16 @@ end @trixi_testset "elixir_mhd_shockcapturing_amr.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_shockcapturing_amr.jl"), - l2=[0.006298541670176575, 0.0064368506652601265, - 0.007108729762852636, 0.006530420607206385, - 0.02061185869237284, 0.005562033787605515, - 0.007571716276627825, 0.005571862660453231, - 3.909755063709152e-6], - linf=[0.20904054009050665, 0.18622917151105936, - 0.2347957890323218, 0.19432508025509926, - 0.6858860133405615, 0.15172116633332622, - 0.22432820727833747, 0.16805989780225183, - 0.000535219040687628], + l2=[0.006297229188299052, 0.0064363477630573936, + 0.007109134822960387, 0.0065295379843073945, + 0.02061487028361094, 0.005561406556868266, + 0.007570747563219415, 0.005571060186624124, + 3.910359570546058e-6], + linf=[0.20904050617411984, 0.18630026905465372, + 0.23476537952044518, 0.19430178061639747, + 0.6858488631108304, 0.15169972134884624, + 0.22431157069631724, 0.16823638724229162, + 0.0005352202836463904], tspan=(0.0, 0.04), coverage_override=(maxiters = 6, initial_refinement_level = 1, base_level = 1, max_level = 2)) @@ -586,6 +586,43 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + +@trixi_testset "elixir_euler_weak_blast_wave_amr.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weak_blast_wave_amr.jl"), + l2=[ + 0.011345993108796831, + 0.018525073963833696, + 0.019102348105917946, + 0.01920515438943838, + 0.15060493968460148, + ], + linf=[ + 0.2994949779783401, + 0.5530175050084679, + 0.5335803757792128, + 0.5647252867336123, + 3.6462732329242566, + ], + tspan=(0.0, 0.025), + coverage_override=(maxiters = 6,)) + # 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 + # Check for conservation + state_integrals = Trixi.integrate(sol.u[2], semi) + initial_state_integrals = analysis_callback.affect!.initial_state_integrals + + @test isapprox(state_integrals[1], initial_state_integrals[1], atol = 1e-13) + @test isapprox(state_integrals[2], initial_state_integrals[2], atol = 1e-13) + @test isapprox(state_integrals[3], initial_state_integrals[3], atol = 1e-13) + @test isapprox(state_integrals[4], initial_state_integrals[4], atol = 1e-13) + @test isapprox(state_integrals[5], initial_state_integrals[5], atol = 1e-13) +end end # Clean up afterwards: delete Trixi.jl output directory diff --git a/test/test_t8code_2d.jl b/test/test_t8code_2d.jl index b63d2a105ac..c1fcc355218 100644 --- a/test/test_t8code_2d.jl +++ b/test/test_t8code_2d.jl @@ -121,8 +121,8 @@ end @trixi_testset "elixir_advection_amr_unstructured_flag.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr_unstructured_flag.jl"), - l2=[0.001993165013217687], - linf=[0.032891018571625796], + l2=[0.002019623611753929], + linf=[0.03542375961299987], coverage_override=(maxiters = 6,)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -305,16 +305,16 @@ end # This test is identical to the one in `test_p4est_2d.jl` besides minor # deviations in the expected error norms. @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_rotor.jl"), - l2=[0.44211360369891683, 0.8805178316216257, 0.8262710688468049, + l2=[0.44207324634847545, 0.8804644301177857, 0.8262542320669426, 0.0, - 0.9616090460973586, 0.10386643568745411, - 0.15403457366543802, 0.0, - 2.8399715649715473e-5], - linf=[10.04369305341599, 17.995640564998403, 9.576041548174265, + 0.9615023124189027, 0.10386709616755131, 0.1540308191628843, 0.0, - 19.429658884314534, 1.3821395681242314, 1.818559351543182, + 2.8350276854372125e-5], + linf=[10.04548675437385, 17.998696852394836, 9.575802136190026, 0.0, - 0.002261930217575465], + 19.431290746184473, 1.3821685018474321, 1.8186235976551453, + 0.0, + 0.002309422702635547], tspan=(0.0, 0.02)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -325,6 +325,40 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + +@trixi_testset "elixir_euler_weak_blast_wave_amr.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weak_blast_wave_amr.jl"), + l2=[ + 0.10823279736983638, + 0.1158152939803735, + 0.11633970342992006, + 0.751152651902375, + ], + linf=[ + 0.5581611332828653, + 0.8354026029724041, + 0.834485181423738, + 3.923553028014343, + ], + tspan=(0.0, 0.1), + coverage_override=(maxiters = 6,)) + # 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 + # Check for conservation + state_integrals = Trixi.integrate(sol.u[2], semi) + initial_state_integrals = analysis_callback.affect!.initial_state_integrals + + @test isapprox(state_integrals[1], initial_state_integrals[1], atol = 1e-13) + @test isapprox(state_integrals[2], initial_state_integrals[2], atol = 1e-13) + @test isapprox(state_integrals[3], initial_state_integrals[3], atol = 1e-13) + @test isapprox(state_integrals[4], initial_state_integrals[4], atol = 1e-13) +end end # Clean up afterwards: delete Trixi.jl output directory diff --git a/test/test_t8code_3d.jl b/test/test_t8code_3d.jl index 5c452444be0..940d2c43372 100644 --- a/test/test_t8code_3d.jl +++ b/test/test_t8code_3d.jl @@ -61,7 +61,7 @@ mkdir(outdir) @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_nonconforming.jl"), l2=[0.00253595715323843], linf=[0.016486952252155795]) - # Ensure that we do not have excessive memory allocations + # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let t = sol.t[end] @@ -80,7 +80,7 @@ mkdir(outdir) linf=[0.0007889950196294793], coverage_override=(maxiters = 6, initial_refinement_level = 1, base_level = 1, med_level = 2, max_level = 3)) - # Ensure that we do not have excessive memory allocations + # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let t = sol.t[end] @@ -95,12 +95,12 @@ mkdir(outdir) @trixi_testset "elixir_advection_amr_unstructured_curved.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr_unstructured_curved.jl"), - l2=[2.0556575425846923e-5], - linf=[0.00105682693484822], + l2=[2.0535121347526814e-5], + linf=[0.0010586603797777504], tspan=(0.0, 1.0), coverage_override=(maxiters = 6, initial_refinement_level = 0, base_level = 0, med_level = 1, max_level = 2)) - # Ensure that we do not have excessive memory allocations + # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let t = sol.t[end] @@ -129,7 +129,7 @@ mkdir(outdir) 0.008526972236273522, ], tspan=(0.0, 0.01)) - # Ensure that we do not have excessive memory allocations + # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let t = sol.t[end] @@ -158,7 +158,7 @@ mkdir(outdir) 0.01562861968368434, ], tspan=(0.0, 1.0)) - # Ensure that we do not have excessive memory allocations + # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let t = sol.t[end] @@ -186,7 +186,7 @@ mkdir(outdir) 9.412914891981927e-12, ], tspan=(0.0, 0.03)) - # Ensure that we do not have excessive memory allocations + # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let t = sol.t[end] @@ -214,7 +214,7 @@ mkdir(outdir) 9.592326932761353e-13, ], tspan=(0.0, 0.1), atol=5.0e-13,) - # Ensure that we do not have excessive memory allocations + # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let t = sol.t[end] @@ -243,7 +243,7 @@ mkdir(outdir) ], tspan=(0.0, 0.2), coverage_override=(polydeg = 3,)) # Prevent long compile time in CI - # Ensure that we do not have excessive memory allocations + # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let t = sol.t[end] @@ -273,7 +273,7 @@ mkdir(outdir) ], tspan=(0.0, 0.3), coverage_override=(polydeg = 3,)) # Prevent long compile time in CI - # Ensure that we do not have excessive memory allocations + # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let t = sol.t[end] @@ -316,6 +316,43 @@ mkdir(outdir) @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + + @trixi_testset "elixir_euler_weak_blast_wave_amr.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weak_blast_wave_amr.jl"), + l2=[ + 0.010014531529951328, + 0.0176268986746271, + 0.01817514447099777, + 0.018271085903740675, + 0.15193033077438198, + ], + linf=[ + 0.2898958869606375, + 0.529717119064458, + 0.5567193302705906, + 0.570663236219957, + 3.5496520808512027, + ], + tspan=(0.0, 0.025), + coverage_override=(maxiters = 6,)) + # 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 + # Check for conservation + state_integrals = Trixi.integrate(sol.u[2], semi) + initial_state_integrals = analysis_callback.affect!.initial_state_integrals + + @test isapprox(state_integrals[1], initial_state_integrals[1], atol = 1e-13) + @test isapprox(state_integrals[2], initial_state_integrals[2], atol = 1e-13) + @test isapprox(state_integrals[3], initial_state_integrals[3], atol = 1e-13) + @test isapprox(state_integrals[4], initial_state_integrals[4], atol = 1e-13) + @test isapprox(state_integrals[5], initial_state_integrals[5], atol = 1e-13) + end end # Clean up afterwards: delete Trixi.jl output directory From 1c1d5ea0f7fc763dd3954121fe665c6001e21df1 Mon Sep 17 00:00:00 2001 From: Huiyu Xie Date: Mon, 19 Aug 2024 20:18:01 -1000 Subject: [PATCH 6/6] Add numerical support of other real types (`laplace`) (#2025) * start * reformat * Update test/test_type.jl Co-authored-by: Hendrik Ranocha * Update test/test_type.jl Co-authored-by: Hendrik Ranocha * Update test/test_type.jl Co-authored-by: Hendrik Ranocha --------- Co-authored-by: Hendrik Ranocha --- test/test_type.jl | 63 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 63 insertions(+) diff --git a/test/test_type.jl b/test/test_type.jl index c15ac5f78c6..d13b626b060 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -1459,6 +1459,69 @@ isdir(outdir) && rm(outdir, recursive = true) end end + @timed_testset "Laplace Diffusion 1D" begin + for RealT in (Float32, Float64) + equations = @inferred LinearScalarAdvectionEquation1D(RealT(1)) + equations_parabolic = @inferred LaplaceDiffusion1D(RealT(0.1), equations) + + x = SVector(zero(RealT)) + t = zero(RealT) + u = gradients = SVector(one(RealT)) + orientation = 1 + + @test eltype(@inferred flux(u, gradients, orientation, equations_parabolic)) == + RealT + + # TODO: BC tests for BoundaryConditionDirichlet + # TODO: BC tests for BoundaryConditionNeumann + end + end + + @timed_testset "Laplace Diffusion 2D" begin + for RealT in (Float32, Float64) + equations = LinearScalarAdvectionEquation2D(RealT(1), RealT(1)) + equations_parabolic = LaplaceDiffusion2D(RealT(0.1), equations) + + x = SVector(zero(RealT), zero(RealT)) + t = zero(RealT) + u = u_inner = u_outer = inv_h = gradients = SVector(one(RealT), one(RealT)) + orientations = [1, 2] + + for orientation in orientations + @test eltype(@inferred flux(u, gradients, orientation, equations_parabolic)) == + RealT + end + + parabolic_solver = ViscousFormulationLocalDG(RealT(0.1)) + @test eltype(@inferred Trixi.penalty(u_outer, u_inner, inv_h, + equations_parabolic, parabolic_solver)) == + RealT + end + end + + @timed_testset "Laplace Diffusion 3D" begin + for RealT in (Float32, Float64) + equations = LinearScalarAdvectionEquation3D(RealT(1), RealT(1), RealT(1)) + equations_parabolic = LaplaceDiffusion3D(RealT(0.1), equations) + + x = SVector(zero(RealT), zero(RealT), zero(RealT)) + t = zero(RealT) + u = u_inner = u_outer = inv_h = gradients = SVector(one(RealT), one(RealT), + one(RealT)) + orientations = [1, 2, 3] + + for orientation in orientations + @test eltype(@inferred flux(u, gradients, orientation, equations_parabolic)) == + RealT + end + + parabolic_solver = ViscousFormulationLocalDG(RealT(0.1)) + @test eltype(@inferred Trixi.penalty(u_outer, u_inner, inv_h, + equations_parabolic, parabolic_solver)) == + RealT + end + end + @timed_testset "Lattice Boltzmann 2D" begin for RealT in (Float32, Float64) equations = @inferred LatticeBoltzmannEquations2D(Ma = RealT(0.1), Re = 1000)