From b346c8f08df9b24c1d299ca5388704ccba6ad340 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 31 May 2022 07:29:54 -0500 Subject: [PATCH] Additional work on parabolic terms (#1148) * make "do nothing BC" a struct, move it to basic_types * add varnames for LaplaceDiffusion * add forgotten @threaded * update comment * Revert "add forgotten @threaded" This reverts commit 1564501a64ecba856ef9cc28bda1a34278ff1eb6. * add @threaded (without enumerate) * remove periodic advection diffusion elixir (redundant) * adding test * one more test * Update src/basic_types.jl Co-authored-by: Hendrik Ranocha Co-authored-by: Jesse Chan Co-authored-by: Hendrik Ranocha --- .../elixir_advection_diffusion_periodic.jl | 36 ------------------- src/basic_types.jl | 14 ++++++++ src/equations/laplace_diffusion_2d.jl | 3 ++ src/solvers/dgmulti/dg.jl | 5 ++- test/test_parabolic_2d.jl | 7 ++-- 5 files changed, 24 insertions(+), 41 deletions(-) delete mode 100644 examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl diff --git a/examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl b/examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl deleted file mode 100644 index ceeb01f71f1..00000000000 --- a/examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl +++ /dev/null @@ -1,36 +0,0 @@ - -using Trixi, OrdinaryDiffEq - -dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial(), - surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs), - volume_integral = VolumeIntegralWeakForm()) - -equations = LinearScalarAdvectionEquation2D(1.0, 1.0) -equations_parabolic = LaplaceDiffusion2D(1e-2, equations) - -function initial_condition_sharp_gaussian(x, t, equations::LinearScalarAdvectionEquation2D) - return SVector(exp(-100 * (x[1]^2 + x[2]^2))) -end -initial_condition = initial_condition_sharp_gaussian - -mesh = DGMultiMesh(dg, cells_per_dimension = (16, 16), periodicity=true) -semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), - initial_condition, dg) - -tspan = (0.0, 2.0) -ode = semidiscretize(semi, tspan) - -summary_callback = SummaryCallback() -alive_callback = AliveCallback(alive_interval=10) -analysis_interval = 100 -analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg)) -callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback) - -############################################################################### -# run the simulation - -time_int_tol = 1e-6 -sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol, - save_everystep=false, callback=callbacks) - -summary_callback() # print the timer summary diff --git a/src/basic_types.jl b/src/basic_types.jl index d6ae7f9fceb..1116f3e6b66 100644 --- a/src/basic_types.jl +++ b/src/basic_types.jl @@ -70,5 +70,19 @@ const boundary_condition_periodic = BoundaryConditionPeriodic() Base.show(io::IO, ::BoundaryConditionPeriodic) = print(io, "boundary_condition_periodic") +""" + boundary_condition_do_nothing = BoundaryConditionDoNothing() + +Imposing no boundary condition just evaluates the flux at the inner state. +""" +struct BoundaryConditionDoNothing end + +@inline function (boundary_condition::BoundaryConditionDoNothing)(inner_flux_or_state, other_args...) + return inner_flux_or_state +end + +const boundary_condition_do_nothing = BoundaryConditionDoNothing() + +Base.show(io::IO, ::BoundaryConditionDoNothing) = print(io, "boundary_condition_do_nothing") end # @muladd diff --git a/src/equations/laplace_diffusion_2d.jl b/src/equations/laplace_diffusion_2d.jl index 774d93cf7b3..84374f5e0f2 100644 --- a/src/equations/laplace_diffusion_2d.jl +++ b/src/equations/laplace_diffusion_2d.jl @@ -12,6 +12,9 @@ end LaplaceDiffusion2D(diffusivity, equations) = LaplaceDiffusion2D{typeof(equations), nvariables(equations), typeof(diffusivity)}(diffusivity, equations) +varnames(variable_mapping, equations_parabolic::LaplaceDiffusion2D) = + varnames(variable_mapping, equations_parabolic.equations) + # no orientation specified since the flux is vector-valued function flux(u, grad_u, equations::LaplaceDiffusion2D) dudx, dudy = grad_u diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl index 6f2a835dbd4..c4abe2e91b7 100644 --- a/src/solvers/dgmulti/dg.jl +++ b/src/solvers/dgmulti/dg.jl @@ -31,9 +31,8 @@ mul_by_accum!(A::UniformScaling) = MulByAccumUniformScaling() # solution storage formats. @inline apply_to_each_field(f::MulByUniformScaling, out, x, args...) = copy!(out, x) @inline function apply_to_each_field(f::MulByAccumUniformScaling, out, x, args...) - # TODO: DGMulti speed up using threads - for (i, x_i) in enumerate(x) - out[i] = out[i] + x_i + @threaded for i in eachindex(x) + out[i] = out[i] + x[i] end end diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 8da8468c8c0..8e8df9f43f5 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -30,6 +30,7 @@ isdir(outdir) && rm(outdir, recursive=true) semi = SemidiscretizationHyperbolicParabolic(mesh, equations, equations_parabolic, initial_condition, dg) @test_nowarn_debug show(stdout, semi) @test_nowarn_debug show(stdout, MIME"text/plain"(), semi) + @test_nowarn_debug show(stdout, boundary_condition_do_nothing) @test nvariables(semi)==nvariables(equations) @test Base.ndims(semi)==Base.ndims(mesh) @@ -40,6 +41,9 @@ isdir(outdir) && rm(outdir, recursive=true) Trixi.compute_coefficients!(u0, 0.0, semi) @test u0 ≈ ode.u0 + # test "do nothing" BC just returns first argument + @test boundary_condition_do_nothing(u0, nothing) == u0 + @unpack cache, cache_parabolic, equations_parabolic = semi @unpack u_grad = cache_parabolic for dim in eachindex(u_grad) @@ -47,8 +51,7 @@ isdir(outdir) && rm(outdir, recursive=true) end t = 0.0 - # pass in `boundary_condition_periodic` to fake "do-nothing" - # TODO: DGMulti. Make `boundary_condition_do_nothing` a callable struct like BoundaryConditionPeriodic + # pass in `boundary_condition_periodic` to skip boundary flux/integral evaluation Trixi.calc_gradient!(u_grad, ode.u0, t, mesh, equations_parabolic, boundary_condition_periodic, dg, cache, cache_parabolic) @unpack x, y = mesh.md