diff --git a/NEWS.md b/NEWS.md index 24c7e28af6f..58b3783f533 100644 --- a/NEWS.md +++ b/NEWS.md @@ -7,6 +7,23 @@ for human readability. ## Changes when updating to v0.16 from v0.15.x +#### Changed + +- The implementation of the local DG (`ViscousFormulationLocalDG`) `solver_parabolic` has been changed for the `P4estMesh`. +In particular, instead of computing the `ldg_switch` as the dot product of the normal direction with ones, +i.e., summing up the normal components, the `ldg_switch` is now selected as +the sign of the maximum (in absolute value sense) normal direction component, +which corresponds to the dominant direction of the interface normal. +This might change results slightly for some meshes where the sum of the normal might be close to zero, +thus introducing some spurious switch assignments ([#2871]). +- The word "viscous" is now used only where it refers specifically to fluid viscosity. +The word "parabolic" is used in more general contexts. +In particular, viscosity is no longer used as a proxy for any parabolic/diffusive process such as heat conduction. +For example, `ViscousFormulationLocalDG` is now `ParabolicFormulationLocalDG` and +`ViscousFormulationBassiRebay1` is now `ParabolicFormulationBassiRebay1`. +For consistency, `cfl_advective` and `cfl_diffusive` have also been renamed `cfl_hyperbolic` and `cfl_parabolic` ([#2868]). +Moreover, some internal functions have been renamed accordingly, including the results shown by the timer outputs after running a simulation. + #### Added - Introducing GPU support: Based on work by Jan Kraus and Lars Christmann, Trixi.jl can @@ -22,15 +39,6 @@ for human readability. GPU kernels are currently CI-tested on NVIDIA GPUs in a buildkite workflow using `TRIXI_TEST=CUDA` ([#2590]). -#### Changed - -- The implementation of the local DG (`ViscousFormulationLocalDG`) `solver_parabolic` has been changed for the `P4estMesh`. -In particular, instead of computing the `ldg_switch` as the dot product of the normal direction with ones, -i.e., summing up the normal components, the `ldg_switch` is now selected as -the sign of the maximum (in absolute value sense) normal direction component, -which corresponds to the dominant direction of the interface normal. -This might change results slightly for some meshes where the sum of the normal might be close to zero, -thus introducing some spurious switch assignments ([#2871]). ## Changes in the v0.15 lifecycle diff --git a/docs/literate/src/files/adding_new_parabolic_terms.jl b/docs/literate/src/files/adding_new_parabolic_terms.jl index 9ec30998eb8..1b252c89785 100644 --- a/docs/literate/src/files/adding_new_parabolic_terms.jl +++ b/docs/literate/src/files/adding_new_parabolic_terms.jl @@ -35,14 +35,14 @@ function varnames(variable_mapping, equations_parabolic::ConstantAnisotropicDiff return varnames(variable_mapping, equations_parabolic.equations_hyperbolic) end -# Next, we define the viscous flux function. We assume that the mixed hyperbolic-parabolic system +# Next, we define the parabolic flux function. We assume that the mixed hyperbolic-parabolic system # is of the form # ```math # \partial_t u(t,x) + \partial_x (f_1(u) - g_1(u, \nabla u)) # + \partial_y (f_2(u) - g_2(u, \nabla u)) = 0 # ``` # where ``f_1(u)``, ``f_2(u)`` are the hyperbolic fluxes and ``g_1(u, \nabla u)``, ``g_2(u, \nabla u)`` denote -# the viscous fluxes. For anisotropic diffusion, the viscous fluxes are the first and second components +# the parabolic fluxes. For anisotropic diffusion, the parabolic fluxes are the first and second components # of the matrix-vector product involving `diffusivity` and the gradient vector. # # Here, we specialize the flux to our new parabolic equation type `ConstantAnisotropicDiffusion2D`. @@ -66,12 +66,12 @@ end # \begin{aligned} # \bm{q} &= \nabla u \\ # \bm{\sigma} &= \begin{pmatrix} g_1(u, \bm{q}) \\ g_2(u, \bm{q}) \end{pmatrix} \\ -# \text{viscous contribution } &= \nabla \cdot \bm{\sigma} +# \text{parabolic contribution} &= \nabla \cdot \bm{\sigma} # \end{aligned} # ``` # # Boundary data must be specified for all spatial derivatives, e.g., for both the gradient -# equation ``\bm{q} = \nabla u`` and the divergence of the viscous flux +# equation ``\bm{q} = \nabla u`` and the divergence of the parabolic flux # ``\nabla \cdot \bm{\sigma}``. We account for this by introducing internal `Gradient` # and `Divergence` types which are used to dispatch on each type of boundary condition. # @@ -98,7 +98,7 @@ end return boundary_condition.boundary_value end -# While the gradient acts on the solution `u`, the divergence acts on the viscous flux ``\bm{\sigma}``. +# While the gradient acts on the solution `u`, the divergence acts on the parabolic flux ``\bm{\sigma}``. # Thus, we have to supply boundary data for the `Divergence` operator that corresponds to ``\bm{\sigma}``. # However, we've already imposed boundary data on `u` for a Dirichlet boundary condition, and imposing # boundary data for ``\bm{\sigma}`` might overconstrain our problem. @@ -119,7 +119,7 @@ end # ### A note on the choice of gradient variables # # It is often simpler to transform the solution variables (and solution gradients) to another set of -# variables prior to computing the viscous fluxes (see [`CompressibleNavierStokesDiffusion2D`](@ref) +# variables prior to computing the parabolic fluxes (see [`CompressibleNavierStokesDiffusion2D`](@ref) # for an example of this). If this is done, then the boundary condition for the `Gradient` operator # should be modified accordingly as well. # @@ -182,7 +182,7 @@ plot(sol) # To be able to do so, we need to define [`max_diffusivity`](@ref) and # [`have_constant_diffusivity`](@ref) for the new parabolic terms. # In Trixi.jl, currently only the standard Laplace Diffusion and Compressible Navier-Stokes-Fourier -# viscous terms are implemented. +# parabolic terms are implemented. # Since these equations have **isotropic** diffusivity, i.e., direction-independent coefficients, # [`max_diffusivity`](@ref) is expected to return a scalar value. # @@ -204,11 +204,11 @@ end return lambda_max() end -# We supply now the advective(hyperbolic) and diffusive(parabolic) CFL numbers -cfl_advective = 2.0 # Not restrictive for this example -cfl_diffusive = 0.21 # Restricts the timestep -stepsize_callback = StepsizeCallback(cfl = cfl_advective, - cfl_diffusive = cfl_diffusive) +# We now supply the hyperbolic and parabolic CFL numbers +cfl_hyperbolic = 2.0 # Not restrictive for this example +cfl_parabolic = 0.21 # Restricts the timestep +stepsize_callback = StepsizeCallback(cfl = cfl_hyperbolic, + cfl_parabolic = cfl_parabolic) # Add the stepsize callback to the existing callbacks callbacks = CallbackSet(SummaryCallback(), stepsize_callback); diff --git a/docs/literate/src/files/parabolic_source_terms.jl b/docs/literate/src/files/parabolic_source_terms.jl index 458510ce6b0..f98f340da7f 100644 --- a/docs/literate/src/files/parabolic_source_terms.jl +++ b/docs/literate/src/files/parabolic_source_terms.jl @@ -62,7 +62,7 @@ end # to OrdinaryDiffEq.jl. # # Note that for this problem, since viscosity `nu` is relatively large, we utilize -# `ViscousFormulationLocalDG` instead of the default `ViscousFormulationBassiRebay1` +# `ParabolicFormulationLocalDG` instead of the default `ParabolicFormulationBassiRebay1` # parabolic solver, since the Bassi-Rebay 1 formulation is not accurate when the # diffusivity is large relative to the mesh size. @@ -76,7 +76,7 @@ boundary_conditions_parabolic = boundary_condition_periodic semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; - solver_parabolic = ViscousFormulationLocalDG(), + solver_parabolic = ParabolicFormulationLocalDG(), source_terms = source_terms, source_terms_parabolic = source_terms_parabolic, boundary_conditions = (boundary_conditions, @@ -89,10 +89,10 @@ ode = semidiscretize(semi, tspan) # stable time-step is $O(h^2)$ due to the dominant parabolic term. We enforce this more stringent # parabolic CFL condition using a diffusion-aware `StepsizeCallback`. -cfl_advective = 0.5 -cfl_diffusive = 0.05 -stepsize_callback = StepsizeCallback(cfl = cfl_advective, - cfl_diffusive = cfl_diffusive) +cfl_hyperbolic = 0.5 +cfl_parabolic = 0.05 +stepsize_callback = StepsizeCallback(cfl = cfl_hyperbolic, + cfl_parabolic = cfl_parabolic) callbacks = CallbackSet(SummaryCallback(), stepsize_callback) sol = solve(ode, RDPK3SpFSAL35(); adaptive = false, dt = stepsize_callback(ode), ode_default_options()..., callback = callbacks) diff --git a/examples/dgmulti_1d/elixir_advection_diffusion_gradient_source_terms.jl b/examples/dgmulti_1d/elixir_advection_diffusion_gradient_source_terms.jl index 29ee42de3f3..f2e7b8a542e 100644 --- a/examples/dgmulti_1d/elixir_advection_diffusion_gradient_source_terms.jl +++ b/examples/dgmulti_1d/elixir_advection_diffusion_gradient_source_terms.jl @@ -58,10 +58,10 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval, uEltype alive_callback = AliveCallback(analysis_interval = 100) -cfl_advective = 0.5 # Not restrictive for this example -cfl_diffusive = 0.025 # Restricts the timestep -stepsize_callback = StepsizeCallback(cfl = cfl_advective, - cfl_diffusive = cfl_diffusive) +cfl_hyperbolic = 0.5 # Not restrictive for this example +cfl_parabolic = 0.025 # Restricts the timestep +stepsize_callback = StepsizeCallback(cfl = cfl_hyperbolic, + cfl_parabolic = cfl_parabolic) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, stepsize_callback) diff --git a/examples/dgmulti_2d/elixir_advection_diffusion.jl b/examples/dgmulti_2d/elixir_advection_diffusion.jl index 8eed206f315..df9935f1d54 100644 --- a/examples/dgmulti_2d/elixir_advection_diffusion.jl +++ b/examples/dgmulti_2d/elixir_advection_diffusion.jl @@ -34,7 +34,7 @@ boundary_conditions = (; left = boundary_condition_left, top = boundary_condition_do_nothing, right = boundary_condition_do_nothing) -# define viscous boundary conditions +# define parabolic boundary conditions boundary_conditions_parabolic = (; left = boundary_condition_left, bottom = boundary_condition_zero, top = boundary_condition_zero, diff --git a/examples/dgmulti_2d/elixir_advection_diffusion_nonperiodic.jl b/examples/dgmulti_2d/elixir_advection_diffusion_nonperiodic.jl index 59ffbf181ea..4f339968256 100644 --- a/examples/dgmulti_2d/elixir_advection_diffusion_nonperiodic.jl +++ b/examples/dgmulti_2d/elixir_advection_diffusion_nonperiodic.jl @@ -53,7 +53,7 @@ boundary_conditions = (; left = boundary_condition, bottom = boundary_condition, right = boundary_condition_do_nothing) -# define viscous boundary conditions +# define parabolic boundary conditions boundary_conditions_parabolic = (; entire_boundary = boundary_condition) semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), diff --git a/examples/dgmulti_2d/elixir_navierstokes_convergence.jl b/examples/dgmulti_2d/elixir_navierstokes_convergence.jl index 63e4285d6b9..c7f0c5ca498 100644 --- a/examples/dgmulti_2d/elixir_navierstokes_convergence.jl +++ b/examples/dgmulti_2d/elixir_navierstokes_convergence.jl @@ -198,7 +198,7 @@ boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_to # define inviscid boundary conditions boundary_conditions = (; top_bottom = boundary_condition_slip_wall) -# define viscous boundary conditions +# define parabolic boundary conditions boundary_conditions_parabolic = (; top_bottom = boundary_condition_top_bottom) semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), diff --git a/examples/dgmulti_2d/elixir_navierstokes_convergence_curved.jl b/examples/dgmulti_2d/elixir_navierstokes_convergence_curved.jl index a4ca0c60791..b0b7601a9a9 100644 --- a/examples/dgmulti_2d/elixir_navierstokes_convergence_curved.jl +++ b/examples/dgmulti_2d/elixir_navierstokes_convergence_curved.jl @@ -206,7 +206,7 @@ boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_to # define inviscid boundary conditions boundary_conditions = (; top_bottom = boundary_condition_slip_wall) -# define viscous boundary conditions +# define parabolic boundary conditions boundary_conditions_parabolic = (; top_bottom = boundary_condition_top_bottom) semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), diff --git a/examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl b/examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl index 38c8234b839..167bc54e57c 100644 --- a/examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl +++ b/examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl @@ -51,7 +51,7 @@ boundary_condition_cavity = BoundaryConditionNavierStokesWall(velocity_bc_cavity boundary_conditions = (; top = boundary_condition_slip_wall, rest_of_boundary = boundary_condition_slip_wall) -# define viscous boundary conditions +# define parabolic boundary conditions boundary_conditions_parabolic = (; top = boundary_condition_lid, rest_of_boundary = boundary_condition_cavity) diff --git a/examples/dgmulti_3d/elixir_navierstokes_convergence.jl b/examples/dgmulti_3d/elixir_navierstokes_convergence.jl index 9adb48efa0a..3c174196457 100644 --- a/examples/dgmulti_3d/elixir_navierstokes_convergence.jl +++ b/examples/dgmulti_3d/elixir_navierstokes_convergence.jl @@ -241,7 +241,7 @@ boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_to # define inviscid boundary conditions boundary_conditions = (; top_bottom = boundary_condition_slip_wall) -# define viscous boundary conditions +# define parabolic boundary conditions boundary_conditions_parabolic = (; top_bottom = boundary_condition_top_bottom) semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), diff --git a/examples/dgmulti_3d/elixir_navierstokes_convergence_curved.jl b/examples/dgmulti_3d/elixir_navierstokes_convergence_curved.jl index 96469770bd1..bdc997f8953 100644 --- a/examples/dgmulti_3d/elixir_navierstokes_convergence_curved.jl +++ b/examples/dgmulti_3d/elixir_navierstokes_convergence_curved.jl @@ -249,7 +249,7 @@ boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_to # define inviscid boundary conditions boundary_conditions = (; top_bottom = boundary_condition_slip_wall) -# define viscous boundary conditions +# define parabolic boundary conditions boundary_conditions_parabolic = (; top_bottom = boundary_condition_top_bottom) semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), diff --git a/examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl b/examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl index 5490425f558..eec4ea0d79b 100644 --- a/examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl +++ b/examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl @@ -50,7 +50,7 @@ boundary_conditions_parabolic = BoundaryConditionDirichlet(initial_condition) semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; - solver_parabolic = ViscousFormulationBassiRebay1(), + solver_parabolic = ParabolicFormulationBassiRebay1(), boundary_conditions = (boundary_conditions, boundary_conditions_parabolic)) diff --git a/examples/p4est_2d_dgsem/elixir_advection_diffusion_rotated.jl b/examples/p4est_2d_dgsem/elixir_advection_diffusion_rotated.jl index 5bf67a09000..110e87937e6 100644 --- a/examples/p4est_2d_dgsem/elixir_advection_diffusion_rotated.jl +++ b/examples/p4est_2d_dgsem/elixir_advection_diffusion_rotated.jl @@ -32,7 +32,7 @@ mesh = P4estMesh(trees_per_dimension, semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; - solver_parabolic = ViscousFormulationLocalDG(), + solver_parabolic = ParabolicFormulationLocalDG(), boundary_conditions = (boundary_condition_periodic, boundary_condition_periodic)) diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl index 11ae50eec90..5341299e572 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl @@ -199,7 +199,7 @@ boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_to boundary_conditions = (; y_neg = boundary_condition_slip_wall, y_pos = boundary_condition_slip_wall) -# define viscous boundary conditions +# define parabolic boundary conditions boundary_conditions_parabolic = (; y_neg = boundary_condition_top_bottom, y_pos = boundary_condition_top_bottom) diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic.jl index a5025116087..4e11d052d96 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic.jl @@ -203,7 +203,7 @@ boundary_conditions = (; x_neg = boundary_condition_left_right, y_neg = boundary_condition_slip_wall, y_pos = boundary_condition_slip_wall) -# define viscous boundary conditions +# define parabolic boundary conditions boundary_conditions_parabolic = (; x_neg = boundary_condition_left_right, x_pos = boundary_condition_left_right, y_neg = boundary_condition_top_bottom, diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_freestream_ldg.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_freestream_ldg.jl index a5280d1a3e8..3efd44f1afb 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_freestream_ldg.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_freestream_ldg.jl @@ -19,7 +19,7 @@ initial_condition = initial_condition_const polydeg = 3 solver = DGSEM(polydeg = polydeg, surface_flux = flux_lax_friedrichs, volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha)) -solver_parabolic = ViscousFormulationLocalDG() +solver_parabolic = ParabolicFormulationLocalDG() mu() = 0.5 prandtl_number() = 0.72 diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_shearlayer_nonconforming.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_shearlayer_nonconforming.jl index c4595b85180..70fe918a476 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_shearlayer_nonconforming.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_shearlayer_nonconforming.jl @@ -68,7 +68,7 @@ Trixi.refine_p4est!(mesh.p4est, true, refine_fn_c, C_NULL) semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; - solver_parabolic = ViscousFormulationBassiRebay1(), + solver_parabolic = ParabolicFormulationBassiRebay1(), boundary_conditions = (boundary_condition_periodic, boundary_condition_periodic)) diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_vortex_street.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_vortex_street.jl index c206202af03..932926770e9 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_vortex_street.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_vortex_street.jl @@ -119,8 +119,8 @@ function Trixi.get_node_variable(::Val{:vorticity}, u, mesh, equations, dg, cach n_nodes, n_nodes, # equivalent: `ntuple(_ -> n_nodes, ndims(mesh))...,` n_elements) - @unpack viscous_container = cache_parabolic - @unpack gradients = viscous_container + @unpack parabolic_container = cache_parabolic + @unpack gradients = parabolic_container gradients_x, gradients_y = gradients # We can accelerate the computation by thread-parallelizing the loop over elements diff --git a/examples/p4est_3d_dgsem/elixir_advection_diffusion_amr_curved.jl b/examples/p4est_3d_dgsem/elixir_advection_diffusion_amr_curved.jl index 4e39dac1e69..aacca586c23 100644 --- a/examples/p4est_3d_dgsem/elixir_advection_diffusion_amr_curved.jl +++ b/examples/p4est_3d_dgsem/elixir_advection_diffusion_amr_curved.jl @@ -64,7 +64,7 @@ mesh = P4estMesh{3}(mesh_file, polydeg = 2, semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; - solver_parabolic = ViscousFormulationBassiRebay1(), + solver_parabolic = ParabolicFormulationBassiRebay1(), boundary_conditions = (boundary_conditions, boundary_conditions)) diff --git a/examples/p4est_3d_dgsem/elixir_advection_diffusion_nonconforming.jl b/examples/p4est_3d_dgsem/elixir_advection_diffusion_nonconforming.jl index 4498efe0936..68e6c6ffb6e 100644 --- a/examples/p4est_3d_dgsem/elixir_advection_diffusion_nonconforming.jl +++ b/examples/p4est_3d_dgsem/elixir_advection_diffusion_nonconforming.jl @@ -70,7 +70,7 @@ boundary_conditions_parabolic = BoundaryConditionDirichlet(initial_condition) semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; - solver_parabolic = ViscousFormulationBassiRebay1(), + solver_parabolic = ParabolicFormulationBassiRebay1(), boundary_conditions = (boundary_conditions, boundary_conditions_parabolic)) diff --git a/examples/p4est_3d_dgsem/elixir_advection_diffusion_nonperiodic.jl b/examples/p4est_3d_dgsem/elixir_advection_diffusion_nonperiodic.jl index d098dc0579a..a440f570767 100644 --- a/examples/p4est_3d_dgsem/elixir_advection_diffusion_nonperiodic.jl +++ b/examples/p4est_3d_dgsem/elixir_advection_diffusion_nonperiodic.jl @@ -58,7 +58,7 @@ boundary_conditions = BoundaryConditionDirichlet(initial_condition) semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; - solver_parabolic = ViscousFormulationBassiRebay1(), + solver_parabolic = ParabolicFormulationBassiRebay1(), boundary_conditions = (boundary_conditions, boundary_conditions)) @@ -76,7 +76,7 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval) alive_callback = AliveCallback(analysis_interval = analysis_interval) stepsize_callback = StepsizeCallback(cfl = 1.6, - cfl_diffusive = 0.25) + cfl_parabolic = 0.25) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, stepsize_callback) diff --git a/examples/p4est_3d_dgsem/elixir_navierstokes_blast_wave_amr.jl b/examples/p4est_3d_dgsem/elixir_navierstokes_blast_wave_amr.jl index 497c0a741e5..2c3b563df00 100644 --- a/examples/p4est_3d_dgsem/elixir_navierstokes_blast_wave_amr.jl +++ b/examples/p4est_3d_dgsem/elixir_navierstokes_blast_wave_amr.jl @@ -67,7 +67,7 @@ mesh = P4estMesh(trees_per_dimension, polydeg = 3, semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; - solver_parabolic = ViscousFormulationBassiRebay1(), + solver_parabolic = ParabolicFormulationBassiRebay1(), boundary_conditions = (boundary_condition_periodic, boundary_condition_periodic)) diff --git a/examples/p4est_3d_dgsem/elixir_navierstokes_convergence.jl b/examples/p4est_3d_dgsem/elixir_navierstokes_convergence.jl index 27ccde81b66..f2e7eb08028 100644 --- a/examples/p4est_3d_dgsem/elixir_navierstokes_convergence.jl +++ b/examples/p4est_3d_dgsem/elixir_navierstokes_convergence.jl @@ -243,7 +243,7 @@ boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_to boundary_conditions = (; y_neg = boundary_condition_slip_wall, y_pos = boundary_condition_slip_wall) -# define viscous boundary conditions +# define parabolic boundary conditions boundary_conditions_parabolic = (; y_neg = boundary_condition_top_bottom, y_pos = boundary_condition_top_bottom) diff --git a/examples/p4est_3d_dgsem/elixir_navierstokes_freestream_boundaries.jl b/examples/p4est_3d_dgsem/elixir_navierstokes_freestream_boundaries.jl index 51d19f5302e..4ca0cb88c41 100644 --- a/examples/p4est_3d_dgsem/elixir_navierstokes_freestream_boundaries.jl +++ b/examples/p4est_3d_dgsem/elixir_navierstokes_freestream_boundaries.jl @@ -19,7 +19,7 @@ initial_condition = initial_condition_const polydeg = 3 solver = DGSEM(polydeg = polydeg, surface_flux = flux_lax_friedrichs) -solver_parabolic = ViscousFormulationBassiRebay1() +solver_parabolic = ParabolicFormulationBassiRebay1() mu() = 0.5 prandtl_number() = 0.72 diff --git a/examples/tree_1d_dgsem/elixir_advection_diffusion_cfl.jl b/examples/tree_1d_dgsem/elixir_advection_diffusion_cfl.jl index d06efc92334..d9969ccade6 100644 --- a/examples/tree_1d_dgsem/elixir_advection_diffusion_cfl.jl +++ b/examples/tree_1d_dgsem/elixir_advection_diffusion_cfl.jl @@ -62,10 +62,10 @@ analysis_callback = AnalysisCallback(semi, interval = 100) alive_callback = AliveCallback(analysis_interval = 100) # Stepsize callback which selects the timestep according to the most restrictive CFL condition. -# For coarser grids, linear stability is governed by the advective CFL condition, +# For coarser grids, linear stability is governed by the hyperbolic CFL condition, # while for high refinements the flow becomes diffusion-dominated. stepsize_callback = StepsizeCallback(cfl = 1.6, - cfl_diffusive = 0.3) + cfl_parabolic = 0.3) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, stepsize_callback) diff --git a/examples/tree_1d_dgsem/elixir_advection_diffusion_dirichlet_amr.jl b/examples/tree_1d_dgsem/elixir_advection_diffusion_dirichlet_amr.jl index a2caf612b24..13d8d03f784 100644 --- a/examples/tree_1d_dgsem/elixir_advection_diffusion_dirichlet_amr.jl +++ b/examples/tree_1d_dgsem/elixir_advection_diffusion_dirichlet_amr.jl @@ -42,7 +42,7 @@ boundary_conditions_parabolic = BoundaryConditionDirichlet(initial_condition) semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; - solver_parabolic = ViscousFormulationBassiRebay1(), + solver_parabolic = ParabolicFormulationBassiRebay1(), boundary_conditions = (boundary_conditions, boundary_conditions_parabolic)) diff --git a/examples/tree_1d_dgsem/elixir_advection_diffusion_gradient_source_terms.jl b/examples/tree_1d_dgsem/elixir_advection_diffusion_gradient_source_terms.jl index 7908945fdae..83459b535bc 100644 --- a/examples/tree_1d_dgsem/elixir_advection_diffusion_gradient_source_terms.jl +++ b/examples/tree_1d_dgsem/elixir_advection_diffusion_gradient_source_terms.jl @@ -40,7 +40,7 @@ boundary_conditions_parabolic = boundary_condition_periodic semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; - solver_parabolic = ViscousFormulationLocalDG(), + solver_parabolic = ParabolicFormulationLocalDG(), source_terms = source_terms, source_terms_parabolic = source_terms_parabolic, boundary_conditions = (boundary_conditions, @@ -58,10 +58,10 @@ analysis_callback = AnalysisCallback(semi, interval = 100) alive_callback = AliveCallback(analysis_interval = 100) -cfl_advective = 0.5 -cfl_diffusive = 0.05 -stepsize_callback = StepsizeCallback(cfl = cfl_advective, - cfl_diffusive = cfl_diffusive) +cfl_hyperbolic = 0.5 +cfl_parabolic = 0.05 +stepsize_callback = StepsizeCallback(cfl = cfl_hyperbolic, + cfl_parabolic = cfl_parabolic) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, stepsize_callback) diff --git a/examples/tree_1d_dgsem/elixir_advection_diffusion_ldg.jl b/examples/tree_1d_dgsem/elixir_advection_diffusion_ldg.jl index 4a5c9776e62..8cc4fa8aad9 100644 --- a/examples/tree_1d_dgsem/elixir_advection_diffusion_ldg.jl +++ b/examples/tree_1d_dgsem/elixir_advection_diffusion_ldg.jl @@ -54,7 +54,7 @@ boundary_conditions_parabolic = boundary_condition_periodic semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; - solver_parabolic = ViscousFormulationLocalDG(), + solver_parabolic = ParabolicFormulationLocalDG(), boundary_conditions = (boundary_conditions, boundary_conditions_parabolic)) diff --git a/examples/tree_1d_dgsem/elixir_advection_diffusion_neumann_amr.jl b/examples/tree_1d_dgsem/elixir_advection_diffusion_neumann_amr.jl index 0a8531679b2..e34a594c4fd 100644 --- a/examples/tree_1d_dgsem/elixir_advection_diffusion_neumann_amr.jl +++ b/examples/tree_1d_dgsem/elixir_advection_diffusion_neumann_amr.jl @@ -26,7 +26,7 @@ boundary_condition_neumann_zero = BoundaryConditionNeumann((x, t, equations_para boundary_conditions = (; x_neg = boundary_condition_left, x_pos = boundary_condition_do_nothing) -# define viscous boundary conditions +# define parabolic boundary conditions boundary_conditions_parabolic = (; x_neg = boundary_condition_left, x_pos = boundary_condition_neumann_zero) diff --git a/examples/tree_1d_dgsem/elixir_diffusion_ldg.jl b/examples/tree_1d_dgsem/elixir_diffusion_ldg.jl index 143a9cc9840..8f79de2d3c4 100644 --- a/examples/tree_1d_dgsem/elixir_diffusion_ldg.jl +++ b/examples/tree_1d_dgsem/elixir_diffusion_ldg.jl @@ -41,7 +41,7 @@ boundary_conditions = boundary_condition_periodic boundary_conditions_parabolic = boundary_condition_periodic # A semidiscretization collects data structures and functions for the spatial discretization -solver_parabolic = ViscousFormulationLocalDG() +solver_parabolic = ParabolicFormulationLocalDG() semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; solver_parabolic, diff --git a/examples/tree_1d_dgsem/elixir_diffusion_ldg_newton_krylov.jl b/examples/tree_1d_dgsem/elixir_diffusion_ldg_newton_krylov.jl index 8fc08add88a..f52b62738cd 100644 --- a/examples/tree_1d_dgsem/elixir_diffusion_ldg_newton_krylov.jl +++ b/examples/tree_1d_dgsem/elixir_diffusion_ldg_newton_krylov.jl @@ -33,7 +33,7 @@ function initial_condition_pure_diffusion_1d_convergence_test(x, t, end initial_condition = initial_condition_pure_diffusion_1d_convergence_test -solver_parabolic = ViscousFormulationLocalDG() +solver_parabolic = ParabolicFormulationLocalDG() semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; solver_parabolic, diff --git a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic_cfl.jl b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic_cfl.jl index ea61ce8c547..7f0006b9c1b 100644 --- a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic_cfl.jl +++ b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic_cfl.jl @@ -125,10 +125,10 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval) alive_callback = AliveCallback(analysis_interval = analysis_interval) # Stepsize callback which selects the timestep according to the most restrictive CFL condition. -# For coarser grids, linear stability is governed by the advective/convective CFL condition, +# For coarser grids, linear stability is governed by the hyperbolic CFL condition, # while for high refinements (e.g. initial_refinement_level = 8) the flow becomes diffusion-dominated. stepsize_callback = StepsizeCallback(cfl = 2.7, - cfl_diffusive = 0.2) + cfl_parabolic = 0.2) callbacks = CallbackSet(summary_callback, analysis_callback, diff --git a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls.jl b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls.jl index c0239f5c08d..f31dd16be61 100644 --- a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls.jl +++ b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls.jl @@ -152,7 +152,7 @@ boundary_condition_right = BoundaryConditionNavierStokesWall(velocity_bc_left_ri boundary_conditions = (; x_neg = boundary_condition_slip_wall, x_pos = boundary_condition_slip_wall) -# define viscous boundary conditions +# define parabolic boundary conditions boundary_conditions_parabolic = (; x_neg = boundary_condition_left, x_pos = boundary_condition_right) diff --git a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls_amr.jl b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls_amr.jl index 998f5bf107d..0a387965160 100644 --- a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls_amr.jl +++ b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls_amr.jl @@ -152,7 +152,7 @@ boundary_condition_right = BoundaryConditionNavierStokesWall(velocity_bc_left_ri boundary_conditions = (; x_neg = boundary_condition_slip_wall, x_pos = boundary_condition_slip_wall) -# define viscous boundary conditions +# define parabolic boundary conditions boundary_conditions_parabolic = (; x_neg = boundary_condition_left, x_pos = boundary_condition_right) diff --git a/examples/tree_1d_dgsem/elixir_navierstokes_viscous_shock.jl b/examples/tree_1d_dgsem/elixir_navierstokes_viscous_shock.jl index a7a91921604..5707c6b3efd 100644 --- a/examples/tree_1d_dgsem/elixir_navierstokes_viscous_shock.jl +++ b/examples/tree_1d_dgsem/elixir_navierstokes_viscous_shock.jl @@ -141,7 +141,7 @@ boundary_conditions_parabolic = (; x_neg = boundary_condition_parabolic, # Since this is a diffusion-dominated problem, using the LDG scheme should achieve optimal rates of convergence. # In contrast, BR-1 may achieve suboptimal rates of convergence in diffusion-dominated regimes. # The LDG scheme can be used by specifying the keyword -# solver_parabolic = ViscousFormulationLocalDG() +# solver_parabolic = ParabolicFormulationLocalDG() # in the semidiscretization call below. semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; diff --git a/examples/tree_1d_dgsem/elixir_navierstokes_viscous_shock_imex.jl b/examples/tree_1d_dgsem/elixir_navierstokes_viscous_shock_imex.jl index c3166cab961..428ba79535d 100644 --- a/examples/tree_1d_dgsem/elixir_navierstokes_viscous_shock_imex.jl +++ b/examples/tree_1d_dgsem/elixir_navierstokes_viscous_shock_imex.jl @@ -133,7 +133,7 @@ boundary_conditions_parabolic = (; x_neg = boundary_condition_parabolic, semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; - solver_parabolic = ViscousFormulationLocalDG(), + solver_parabolic = ParabolicFormulationLocalDG(), boundary_conditions = (boundary_conditions, boundary_conditions_parabolic)) diff --git a/examples/tree_1d_dgsem/elixir_viscous_burgers_n_wave.jl b/examples/tree_1d_dgsem/elixir_viscous_burgers_n_wave.jl index 091ab163eae..6555c3fcba2 100644 --- a/examples/tree_1d_dgsem/elixir_viscous_burgers_n_wave.jl +++ b/examples/tree_1d_dgsem/elixir_viscous_burgers_n_wave.jl @@ -53,8 +53,8 @@ analysis_callback = AnalysisCallback(semi, interval = 100) alive_callback = AliveCallback(analysis_interval = 100) -# Timestep is limited by standard/advective/convective CFL -stepsize_callback = StepsizeCallback(cfl = 0.6, cfl_diffusive = 0.1) +# Timestep is limited by the hyperbolic CFL +stepsize_callback = StepsizeCallback(cfl = 0.6, cfl_parabolic = 0.1) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, diff --git a/examples/tree_1d_dgsem/elixir_viscous_burgers_shock.jl b/examples/tree_1d_dgsem/elixir_viscous_burgers_shock.jl index 83bc8498b76..e04f14d640d 100644 --- a/examples/tree_1d_dgsem/elixir_viscous_burgers_shock.jl +++ b/examples/tree_1d_dgsem/elixir_viscous_burgers_shock.jl @@ -51,8 +51,8 @@ analysis_callback = AnalysisCallback(semi, interval = 100) alive_callback = AliveCallback(analysis_interval = 100) -# Timestep is limited by diffusive CFL -stepsize_callback = StepsizeCallback(cfl = 0.8, cfl_diffusive = 0.15) +# Timestep is limited by parabolic CFL +stepsize_callback = StepsizeCallback(cfl = 0.8, cfl_parabolic = 0.15) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, diff --git a/examples/tree_2d_dgsem/elixir_advection_diffusion.jl b/examples/tree_2d_dgsem/elixir_advection_diffusion.jl index ed4ed2432e5..f1d209c7ea6 100644 --- a/examples/tree_2d_dgsem/elixir_advection_diffusion.jl +++ b/examples/tree_2d_dgsem/elixir_advection_diffusion.jl @@ -47,7 +47,7 @@ boundary_conditions_parabolic = boundary_condition_periodic semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; - solver_parabolic = ViscousFormulationBassiRebay1(), + solver_parabolic = ParabolicFormulationBassiRebay1(), boundary_conditions = (boundary_conditions, boundary_conditions_parabolic)) diff --git a/examples/tree_2d_dgsem/elixir_advection_diffusion_amr.jl b/examples/tree_2d_dgsem/elixir_advection_diffusion_amr.jl index b4277434129..d91fd2226bf 100644 --- a/examples/tree_2d_dgsem/elixir_advection_diffusion_amr.jl +++ b/examples/tree_2d_dgsem/elixir_advection_diffusion_amr.jl @@ -43,7 +43,7 @@ boundary_conditions_parabolic = boundary_condition_periodic semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; - solver_parabolic = ViscousFormulationLocalDG(), + solver_parabolic = ParabolicFormulationLocalDG(), boundary_conditions = (boundary_conditions, boundary_conditions_parabolic)) diff --git a/examples/tree_2d_dgsem/elixir_advection_diffusion_gradient_source_terms.jl b/examples/tree_2d_dgsem/elixir_advection_diffusion_gradient_source_terms.jl index 740594ee954..a68f0c1eeda 100644 --- a/examples/tree_2d_dgsem/elixir_advection_diffusion_gradient_source_terms.jl +++ b/examples/tree_2d_dgsem/elixir_advection_diffusion_gradient_source_terms.jl @@ -49,7 +49,7 @@ boundary_conditions_parabolic = boundary_condition_periodic semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; - solver_parabolic = ViscousFormulationLocalDG(), + solver_parabolic = ParabolicFormulationLocalDG(), source_terms = source_terms, source_terms_parabolic = source_terms_parabolic, boundary_conditions = (boundary_conditions, @@ -67,10 +67,10 @@ analysis_callback = AnalysisCallback(semi, interval = 100) alive_callback = AliveCallback(analysis_interval = 100) -cfl_advective = 0.5 -cfl_diffusive = 0.01 -stepsize_callback = StepsizeCallback(cfl = cfl_advective, - cfl_diffusive = cfl_diffusive) +cfl_hyperbolic = 0.5 +cfl_parabolic = 0.01 +stepsize_callback = StepsizeCallback(cfl = cfl_hyperbolic, + cfl_parabolic = cfl_parabolic) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, stepsize_callback) diff --git a/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl b/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl index 93abfd49a3a..2234d8fdf76 100644 --- a/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl +++ b/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl @@ -51,7 +51,7 @@ boundary_conditions_parabolic = BoundaryConditionDirichlet(initial_condition) semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; - solver_parabolic = ViscousFormulationBassiRebay1(), + solver_parabolic = ParabolicFormulationBassiRebay1(), boundary_conditions = (boundary_conditions, boundary_conditions_parabolic)) diff --git a/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl b/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl index 95f9c46f2b8..e7667f9808d 100644 --- a/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl +++ b/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl @@ -49,7 +49,7 @@ boundary_conditions_parabolic = BoundaryConditionDirichlet(initial_condition) semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; - solver_parabolic = ViscousFormulationBassiRebay1(), + solver_parabolic = ParabolicFormulationBassiRebay1(), boundary_conditions = (boundary_conditions, boundary_conditions_parabolic)) diff --git a/examples/tree_2d_dgsem/elixir_diffusion_steady_state_linear_map.jl b/examples/tree_2d_dgsem/elixir_diffusion_steady_state_linear_map.jl index c92140a6613..51530f3286f 100644 --- a/examples/tree_2d_dgsem/elixir_diffusion_steady_state_linear_map.jl +++ b/examples/tree_2d_dgsem/elixir_diffusion_steady_state_linear_map.jl @@ -44,11 +44,11 @@ boundary_conditions = (; x_neg = bc_homogeneous_dirichlet, y_pos = bc_sin_dirichlet, x_pos = bc_homogeneous_dirichlet) -# `solver_parabolic = ViscousFormulationLocalDG()` strictly required for elliptic/diffusion-dominated problem +# `solver_parabolic = ParabolicFormulationLocalDG()` strictly required for elliptic/diffusion-dominated problem semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; - solver_parabolic = ViscousFormulationLocalDG(), + solver_parabolic = ParabolicFormulationLocalDG(), boundary_conditions = (boundary_conditions, boundary_conditions)) diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl b/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl index 16ea4742315..38aed3028c6 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl @@ -205,7 +205,7 @@ boundary_conditions = (; x_neg = boundary_condition_periodic, y_neg = boundary_condition_slip_wall, y_pos = boundary_condition_slip_wall) -# define viscous boundary conditions +# define parabolic boundary conditions boundary_conditions_parabolic = (; x_neg = boundary_condition_periodic, x_pos = boundary_condition_periodic, y_neg = boundary_condition_top_bottom, diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_nonconforming.jl b/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_nonconforming.jl index 257d3293448..ef623715a20 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_nonconforming.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_nonconforming.jl @@ -55,7 +55,7 @@ mesh = TreeMesh(coordinates_min, coordinates_max, semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; - solver_parabolic = ViscousFormulationBassiRebay1(), + solver_parabolic = ParabolicFormulationBassiRebay1(), boundary_conditions = (boundary_condition_periodic, boundary_condition_periodic)) diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_viscous_shock.jl b/examples/tree_2d_dgsem/elixir_navierstokes_viscous_shock.jl index 5f34292a840..4f45baa4a88 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_viscous_shock.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_viscous_shock.jl @@ -153,7 +153,7 @@ boundary_conditions_parabolic = (x_neg = boundary_condition_parabolic, semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; - solver_parabolic = ViscousFormulationBassiRebay1(), + solver_parabolic = ParabolicFormulationBassiRebay1(), boundary_conditions = (boundary_conditions, boundary_conditions_parabolic)) @@ -171,12 +171,12 @@ alive_callback = AliveCallback(alive_interval = 10) analysis_interval = 100 analysis_callback = AnalysisCallback(semi, interval = analysis_interval) -# Admissible stepsize is governed by the diffusive CFL condition. -# Unless the advective cfl number `cfl` is not reduced to e.g. `0.1` +# Admissible stepsize is governed by the parabolic CFL condition. +# Unless the hyperbolic CFL number `cfl` is reduced to e.g. `0.1` # (which is overly restrictive for this problem), -# the diffusive CFL restricts the timestep for this problem. +# the parabolic CFL restricts the timestep for this problem. stepsize_callback = StepsizeCallback(cfl = 0.2, - cfl_diffusive = 0.2) + cfl_parabolic = 0.2) callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback, stepsize_callback) diff --git a/examples/tree_3d_dgsem/elixir_advection_diffusion_amr.jl b/examples/tree_3d_dgsem/elixir_advection_diffusion_amr.jl index 5d886b01cf4..e79f495bef5 100644 --- a/examples/tree_3d_dgsem/elixir_advection_diffusion_amr.jl +++ b/examples/tree_3d_dgsem/elixir_advection_diffusion_amr.jl @@ -43,7 +43,7 @@ boundary_conditions_parabolic = boundary_condition_periodic semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; - solver_parabolic = ViscousFormulationBassiRebay1(), + solver_parabolic = ParabolicFormulationBassiRebay1(), boundary_conditions = (boundary_conditions, boundary_conditions_parabolic)) diff --git a/examples/tree_3d_dgsem/elixir_advection_diffusion_gradient_source_terms.jl b/examples/tree_3d_dgsem/elixir_advection_diffusion_gradient_source_terms.jl index 1a107511344..d8379ba0972 100644 --- a/examples/tree_3d_dgsem/elixir_advection_diffusion_gradient_source_terms.jl +++ b/examples/tree_3d_dgsem/elixir_advection_diffusion_gradient_source_terms.jl @@ -49,7 +49,7 @@ boundary_conditions_parabolic = boundary_condition_periodic semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; - solver_parabolic = ViscousFormulationLocalDG(), + solver_parabolic = ParabolicFormulationLocalDG(), source_terms = source_terms, source_terms_parabolic = source_terms_parabolic, boundary_conditions = (boundary_conditions, @@ -67,10 +67,10 @@ analysis_callback = AnalysisCallback(semi, interval = 100) alive_callback = AliveCallback(analysis_interval = 100) -cfl_advective = 0.5 # Not restrictive for this example -cfl_diffusive = 0.01 # Restricts the timestep -stepsize_callback = StepsizeCallback(cfl = cfl_advective, - cfl_diffusive = cfl_diffusive) +cfl_hyperbolic = 0.5 # Not restrictive for this example +cfl_parabolic = 0.01 # Restricts the timestep +stepsize_callback = StepsizeCallback(cfl = cfl_hyperbolic, + cfl_parabolic = cfl_parabolic) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, stepsize_callback) diff --git a/examples/tree_3d_dgsem/elixir_advection_diffusion_nonconforming.jl b/examples/tree_3d_dgsem/elixir_advection_diffusion_nonconforming.jl index 9ddd6f875e3..226d583b105 100644 --- a/examples/tree_3d_dgsem/elixir_advection_diffusion_nonconforming.jl +++ b/examples/tree_3d_dgsem/elixir_advection_diffusion_nonconforming.jl @@ -55,7 +55,7 @@ boundary_conditions_parabolic = BoundaryConditionDirichlet(initial_condition) semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; - solver_parabolic = ViscousFormulationBassiRebay1(), + solver_parabolic = ParabolicFormulationBassiRebay1(), boundary_conditions = (boundary_conditions, boundary_conditions_parabolic)) diff --git a/examples/tree_3d_dgsem/elixir_advection_diffusion_nonperiodic.jl b/examples/tree_3d_dgsem/elixir_advection_diffusion_nonperiodic.jl index a8a6d27edb0..03d2ffff6c7 100644 --- a/examples/tree_3d_dgsem/elixir_advection_diffusion_nonperiodic.jl +++ b/examples/tree_3d_dgsem/elixir_advection_diffusion_nonperiodic.jl @@ -52,7 +52,7 @@ boundary_conditions_parabolic = BoundaryConditionDirichlet(initial_condition) semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; - solver_parabolic = ViscousFormulationBassiRebay1(), + solver_parabolic = ParabolicFormulationBassiRebay1(), boundary_conditions = (boundary_conditions, boundary_conditions_parabolic)) diff --git a/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl b/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl index 1f1c734ba00..83f60647e61 100644 --- a/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl +++ b/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl @@ -246,7 +246,7 @@ boundary_conditions = (; x_neg = boundary_condition_periodic, z_neg = boundary_condition_periodic, z_pos = boundary_condition_periodic) -# define viscous boundary conditions +# define parabolic boundary conditions boundary_conditions_parabolic = (; x_neg = boundary_condition_periodic, x_pos = boundary_condition_periodic, y_neg = boundary_condition_top_bottom, diff --git a/examples/tree_3d_dgsem/elixir_navierstokes_viscous_shock.jl b/examples/tree_3d_dgsem/elixir_navierstokes_viscous_shock.jl index 8df31a0561a..858029fdd96 100644 --- a/examples/tree_3d_dgsem/elixir_navierstokes_viscous_shock.jl +++ b/examples/tree_3d_dgsem/elixir_navierstokes_viscous_shock.jl @@ -174,10 +174,10 @@ alive_callback = AliveCallback(alive_interval = 10) analysis_interval = 100 analysis_callback = AnalysisCallback(semi, interval = analysis_interval) -# For this setup, both advective and diffusive time step restrictions are relevant, i.e., +# For this setup, both hyperbolic and parabolic timestep restrictions are relevant, i.e., # may not be increased beyond the given values. stepsize_callback = StepsizeCallback(cfl = 0.4, - cfl_diffusive = 0.2) + cfl_parabolic = 0.2) callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback, stepsize_callback) diff --git a/src/Trixi.jl b/src/Trixi.jl index 26b9d17fec7..42340659520 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -337,7 +337,7 @@ export convergence_test, export DGMulti, DGMultiBasis, estimate_dt, DGMultiMesh, GaussSBP -export ViscousFormulationBassiRebay1, ViscousFormulationLocalDG +export ParabolicFormulationBassiRebay1, ParabolicFormulationLocalDG # Visualization-related exports export PlotData1D, PlotData2D, ScalarPlotData2D, getmesh, adapt_to_mesh_level!, diff --git a/src/callbacks_step/amr_dg1d.jl b/src/callbacks_step/amr_dg1d.jl index 8b36edfee91..346b8d6da69 100644 --- a/src/callbacks_step/amr_dg1d.jl +++ b/src/callbacks_step/amr_dg1d.jl @@ -70,8 +70,8 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1}, refine!(u_ode, adaptor, mesh, equations, dg, cache, elements_to_refine) # Resize parabolic helper variables - @unpack viscous_container = cache_parabolic - resize!(viscous_container, equations, dg, cache) + @unpack parabolic_container = cache_parabolic + resize!(parabolic_container, equations, dg, cache) return nothing end @@ -195,8 +195,8 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1}, coarsen!(u_ode, adaptor, mesh, equations, dg, cache, elements_to_remove) # Resize parabolic helper variables - @unpack viscous_container = cache_parabolic - resize!(viscous_container, equations, dg, cache) + @unpack parabolic_container = cache_parabolic + resize!(parabolic_container, equations, dg, cache) return nothing end diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index 2141c2713bc..fd4cceb9efc 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -190,8 +190,8 @@ function refine!(u_ode::AbstractVector, adaptor, refine!(u_ode, adaptor, mesh, equations, dg, cache, elements_to_refine) # Resize parabolic helper variables - @unpack viscous_container = cache_parabolic - resize!(viscous_container, equations, dg, cache) + @unpack parabolic_container = cache_parabolic + resize!(parabolic_container, equations, dg, cache) return nothing end @@ -385,8 +385,8 @@ function coarsen!(u_ode::AbstractVector, adaptor, coarsen!(u_ode, adaptor, mesh, equations, dg, cache, elements_to_remove) # Resize parabolic helper variables - @unpack viscous_container = cache_parabolic - resize!(viscous_container, equations, dg, cache) + @unpack parabolic_container = cache_parabolic + resize!(parabolic_container, equations, dg, cache) return nothing end diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl index 8659119fa36..ab4736582ec 100644 --- a/src/callbacks_step/analysis.jl +++ b/src/callbacks_step/analysis.jl @@ -702,11 +702,11 @@ end # Special analyze for `SemidiscretizationHyperbolicParabolic` such that # precomputed gradients are available. Required for `enstrophy` (see above) and viscous forces. # Note that this needs to be included after `analysis_surface_integral_2d.jl` to -# have `VariableViscous` available. +# have `VariableParabolic` available. function analyze(quantity::AnalysisSurfaceIntegral{Variable}, du, u, t, semi::SemidiscretizationHyperbolicParabolic) where {Variable <: - VariableViscous} + VariableParabolic} mesh, equations, solver, cache = mesh_equations_solver_cache(semi) equations_parabolic = semi.equations_parabolic cache_parabolic = semi.cache_parabolic diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index ccf8843b8aa..368b7f5e4cc 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -397,7 +397,7 @@ function integrate(func::Func, u, mesh::Union{TreeMesh{2}, P4estMesh{2}}, equations, equations_parabolic, dg::DGSEM, cache, cache_parabolic; normalize = true) where {Func} - gradients_x, gradients_y = cache_parabolic.viscous_container.gradients + gradients_x, gradients_y = cache_parabolic.parabolic_container.gradients integrate_via_indices(u, mesh, equations, dg, cache; normalize = normalize) do u, i, j, element, equations, dg u_local = get_node_vars(u, equations, dg, i, j, element) diff --git a/src/callbacks_step/analysis_dg3d.jl b/src/callbacks_step/analysis_dg3d.jl index 09b166646f5..7708a32e6ba 100644 --- a/src/callbacks_step/analysis_dg3d.jl +++ b/src/callbacks_step/analysis_dg3d.jl @@ -445,7 +445,7 @@ function integrate(func::Func, u, mesh::Union{TreeMesh{3}, P4estMesh{3}}, equations, equations_parabolic, dg::DGSEM, cache, cache_parabolic; normalize = true) where {Func} - gradients_x, gradients_y, gradients_z = cache_parabolic.viscous_container.gradients + gradients_x, gradients_y, gradients_z = cache_parabolic.parabolic_container.gradients integrate_via_indices(u, mesh, equations, dg, cache; normalize = normalize) do u, i, j, k, element, equations, dg u_local = get_node_vars(u, equations, dg, i, j, k, element) diff --git a/src/callbacks_step/analysis_surface_integral.jl b/src/callbacks_step/analysis_surface_integral.jl index 1366c8c73a9..33d4262429a 100644 --- a/src/callbacks_step/analysis_surface_integral.jl +++ b/src/callbacks_step/analysis_surface_integral.jl @@ -54,7 +54,7 @@ end # Abstract base type used for dispatch of `analyze` for quantities # requiring gradients of the velocity field. -abstract type VariableViscous end +abstract type VariableParabolic end struct LiftCoefficientPressure{RealT <: Real, NDIMS} force_state::ForceState{RealT, NDIMS} @@ -64,11 +64,11 @@ struct DragCoefficientPressure{RealT <: Real, NDIMS} force_state::ForceState{RealT, NDIMS} end -struct LiftCoefficientShearStress{RealT <: Real, NDIMS} <: VariableViscous +struct LiftCoefficientShearStress{RealT <: Real, NDIMS} <: VariableParabolic force_state::ForceState{RealT, NDIMS} end -struct DragCoefficientShearStress{RealT <: Real, NDIMS} <: VariableViscous +struct DragCoefficientShearStress{RealT <: Real, NDIMS} <: VariableParabolic force_state::ForceState{RealT, NDIMS} end diff --git a/src/callbacks_step/analysis_surface_integral_2d.jl b/src/callbacks_step/analysis_surface_integral_2d.jl index fef4b9872d1..9c1899ef750 100644 --- a/src/callbacks_step/analysis_surface_integral_2d.jl +++ b/src/callbacks_step/analysis_surface_integral_2d.jl @@ -261,7 +261,7 @@ function analyze(surface_variable::AnalysisSurfaceIntegral{Variable}, du, u, t, mesh::P4estMesh{2}, equations, equations_parabolic, dg::DGSEM, cache, semi, - cache_parabolic) where {Variable <: VariableViscous} + cache_parabolic) where {Variable <: VariableParabolic} @unpack boundaries = cache @unpack node_coordinates, contravariant_vectors = cache.elements @unpack weights = dg.basis @@ -271,8 +271,8 @@ function analyze(surface_variable::AnalysisSurfaceIntegral{Variable}, du, u, t, boundary_indices = get_boundary_indices(boundary_symbols, boundary_symbol_indices) # Additions for parabolic - @unpack viscous_container = cache_parabolic - @unpack gradients = viscous_container + @unpack parabolic_container = cache_parabolic + @unpack gradients = parabolic_container gradients_x, gradients_y = gradients diff --git a/src/callbacks_step/stepsize.jl b/src/callbacks_step/stepsize.jl index e81cf027745..aca3098af22 100644 --- a/src/callbacks_step/stepsize.jl +++ b/src/callbacks_step/stepsize.jl @@ -6,30 +6,30 @@ #! format: noindent """ - StepsizeCallback(; cfl=1.0, cfl_diffusive = 0.0, + StepsizeCallback(; cfl=1.0, cfl_parabolic = 0.0, interval = 1) -Set the time step size according to a CFL condition with CFL number `cfl` +Set the time step size according to a CFL condition with hyperbolic CFL number `cfl` if the time integration method isn't adaptive itself. -The keyword argument `cfl` must be either a `Real` number, corresponding to a constant +The hyperbolic CFL number `cfl` must be either a `Real` number, corresponding to a constant CFL number, or a function of time `t` returning a `Real` number. The latter approach allows for variable CFL numbers that can be used to realize, e.g., a ramp-up of the time step. -One can additionally supply a diffusive CFL number `cfl_diffusive` to -limit the admissible timestep also respecting diffusive restrictions. +One can additionally supply a parabolic CFL number `cfl_parabolic` to +limit the admissible timestep also respecting parabolic restrictions. This is only applicable for semidiscretizations of type [`SemidiscretizationHyperbolicParabolic`](@ref). -To enable checking for diffusive timestep restrictions, provide a value greater than zero for `cfl_diffusive`. -By default, `cfl_diffusive` is set to zero which means that only the advective/convective CFL number is considered. -The keyword argument `cfl_diffusive` must be either a `Real` number, corresponding to a constant -diffusive CFL number, or a function of time `t` returning a `Real` number. +To enable checking for parabolic timestep restrictions, provide a value greater than zero for `cfl_parabolic`. +By default, `cfl_parabolic` is set to zero which means that only the hyperbolic CFL number `cfl` is considered. +The keyword argument `cfl_parabolic` must be either a `Real` number, corresponding to a constant +parabolic CFL number, or a function of time `t` returning a `Real` number. By default, the timestep will be adjusted at every step. For different values of `interval`, the timestep will be adjusted every `interval` steps. """ -struct StepsizeCallback{CflAdvectiveType, CflDiffusiveType} - cfl_advective::CflAdvectiveType - cfl_diffusive::CflDiffusiveType +struct StepsizeCallback{CflHyperbolicType, CflParabolicType} + cfl_hyperbolic::CflHyperbolicType + cfl_parabolic::CflParabolicType interval::Int end @@ -37,10 +37,10 @@ function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:StepsizeCallback}) @nospecialize cb # reduce precompilation time stepsize_callback = cb.affect! - @unpack cfl_advective, cfl_diffusive, interval = stepsize_callback + @unpack cfl_hyperbolic, cfl_parabolic, interval = stepsize_callback print(io, "StepsizeCallback(", - "cfl_advective=", cfl_advective, ", ", - "cfl_diffusive=", cfl_diffusive, ", ", + "cfl_hyperbolic=", cfl_hyperbolic, ", ", + "cfl_parabolic=", cfl_parabolic, ", ", "interval=", interval, ")") return nothing end @@ -55,22 +55,22 @@ function Base.show(io::IO, ::MIME"text/plain", stepsize_callback = cb.affect! setup = [ - "CFL Advective" => stepsize_callback.cfl_advective, - "CFL Diffusive" => stepsize_callback.cfl_diffusive, + "CFL Hyperbolic" => stepsize_callback.cfl_hyperbolic, + "CFL Parabolic" => stepsize_callback.cfl_parabolic, "Interval" => stepsize_callback.interval ] summary_box(io, "StepsizeCallback", setup) end end -function StepsizeCallback(; cfl = 1.0, cfl_diffusive = 0.0, +function StepsizeCallback(; cfl = 1.0, cfl_parabolic = 0.0, interval = 1) # Convert plain real numbers to functions for unified treatment - cfl_conv = isa(cfl, Real) ? Returns(cfl) : cfl - cfl_diff = isa(cfl_diffusive, Real) ? Returns(cfl_diffusive) : cfl_diffusive - stepsize_callback = StepsizeCallback{typeof(cfl_conv), typeof(cfl_diff)}(cfl_conv, - cfl_diff, - interval) + cfl_hyp = isa(cfl, Real) ? Returns(cfl) : cfl + cfl_para = isa(cfl_parabolic, Real) ? Returns(cfl_parabolic) : cfl_parabolic + stepsize_callback = StepsizeCallback{typeof(cfl_hyp), typeof(cfl_para)}(cfl_hyp, + cfl_para, + interval) return DiscreteCallback(stepsize_callback, stepsize_callback, # the first one is the condition, the second the affect! save_positions = (false, false), @@ -78,9 +78,9 @@ function StepsizeCallback(; cfl = 1.0, cfl_diffusive = 0.0, end # Compatibility constructor used in `EulerAcousticsCouplingCallback` -function StepsizeCallback(cfl_advective) - RealT = typeof(cfl_advective) - return StepsizeCallback{RealT, RealT}(cfl_advective, zero(RealT), 1) +function StepsizeCallback(cfl_hyperbolic) + RealT = typeof(cfl_hyperbolic) + return StepsizeCallback{RealT, RealT}(cfl_hyperbolic, zero(RealT), 1) end function initialize!(cb::DiscreteCallback{Condition, Affect!}, u, t, @@ -106,11 +106,11 @@ end t = integrator.t u_ode = integrator.u semi = integrator.p - @unpack cfl_advective, cfl_diffusive = stepsize_callback + @unpack cfl_hyperbolic, cfl_parabolic = stepsize_callback # Dispatch based on semidiscretization - dt = @trixi_timeit timer() "calculate dt" calculate_dt(u_ode, t, cfl_advective, - cfl_diffusive, semi) + dt = @trixi_timeit timer() "calculate dt" calculate_dt(u_ode, t, cfl_hyperbolic, + cfl_parabolic, semi) set_proposed_dt!(integrator, dt) integrator.opts.dtmax = dt @@ -130,57 +130,57 @@ function (cb::DiscreteCallback{Condition, Affect!})(ode::ODEProblem) where {Cond StepsizeCallback } stepsize_callback = cb.affect! - @unpack cfl_advective, cfl_diffusive = stepsize_callback + @unpack cfl_hyperbolic, cfl_parabolic = stepsize_callback u_ode = ode.u0 t = first(ode.tspan) semi = ode.p - return calculate_dt(u_ode, t, cfl_advective, cfl_diffusive, semi) + return calculate_dt(u_ode, t, cfl_hyperbolic, cfl_parabolic, semi) end # General case for an abstract single (i.e., non-coupled) semidiscretization -function calculate_dt(u_ode, t, cfl_advective, cfl_diffusive, +function calculate_dt(u_ode, t, cfl_hyperbolic, cfl_parabolic, semi::AbstractSemidiscretization) mesh, equations, solver, cache = mesh_equations_solver_cache(semi) u = wrap_array(u_ode, mesh, equations, solver, cache) - return cfl_advective(t) * max_dt(u, t, mesh, + return cfl_hyperbolic(t) * max_dt(u, t, mesh, have_constant_speed(equations), equations, solver, cache) end # For Euler-Acoustic simulations with `EulerAcousticsCouplingCallback` -function calculate_dt(u_ode, t, cfl_advective::Real, cfl_diffusive::Real, +function calculate_dt(u_ode, t, cfl_hyperbolic::Real, cfl_parabolic::Real, semi::AbstractSemidiscretization) mesh, equations, solver, cache = mesh_equations_solver_cache(semi) u = wrap_array(u_ode, mesh, equations, solver, cache) - return cfl_advective * max_dt(u, t, mesh, + return cfl_hyperbolic * max_dt(u, t, mesh, have_constant_speed(equations), equations, solver, cache) end # Case for a hyperbolic-parabolic semidiscretization -function calculate_dt(u_ode, t, cfl_advective, cfl_diffusive, +function calculate_dt(u_ode, t, cfl_hyperbolic, cfl_parabolic, semi::SemidiscretizationHyperbolicParabolic) mesh, equations, solver, cache = mesh_equations_solver_cache(semi) equations_parabolic = semi.equations_parabolic u = wrap_array(u_ode, mesh, equations, solver, cache) - dt_advective = cfl_advective(t) * max_dt(u, t, mesh, - have_constant_speed(equations), equations, - solver, cache) + dt_hyperbolic = cfl_hyperbolic(t) * max_dt(u, t, mesh, + have_constant_speed(equations), equations, + solver, cache) - cfl_diff = cfl_diffusive(t) - if cfl_diff > 0 # Check if diffusive CFL should be considered - dt_diffusive = cfl_diff * max_dt(u, t, mesh, + cfl_para = cfl_parabolic(t) + if cfl_para > 0 # Check if parabolic CFL should be considered + dt_parabolic = cfl_para * max_dt(u, t, mesh, have_constant_diffusivity(equations_parabolic), equations, equations_parabolic, solver, cache) - return min(dt_advective, dt_diffusive) + return min(dt_hyperbolic, dt_parabolic) else - return dt_advective + return dt_hyperbolic end end diff --git a/src/equations/compressible_navier_stokes.jl b/src/equations/compressible_navier_stokes.jl index 3cfeaaf5018..e15c11f3a02 100644 --- a/src/equations/compressible_navier_stokes.jl +++ b/src/equations/compressible_navier_stokes.jl @@ -115,7 +115,7 @@ dynamic_viscosity(u, mu::T, equations) where {T} = mu(u, equations) # Returns - `False()` -Used in diffusive CFL condition computation (see [`StepsizeCallback`](@ref)) to indicate that the +Used in parabolic CFL condition computation (see [`StepsizeCallback`](@ref)) to indicate that the diffusivity is not constant in space and that [`max_diffusivity`](@ref) needs to be computed at every node in every element. diff --git a/src/equations/compressible_navier_stokes_1d.jl b/src/equations/compressible_navier_stokes_1d.jl index f37a751d3a5..3802c23afd8 100644 --- a/src/equations/compressible_navier_stokes_1d.jl +++ b/src/equations/compressible_navier_stokes_1d.jl @@ -96,7 +96,7 @@ struct CompressibleNavierStokesDiffusion1D{GradientVariables, RealT <: Real, Mu, mu::Mu # viscosity Pr::RealT # Prandtl number kappa::RealT # thermal diffusivity for Fick's law - max_1_kappa::RealT # max(1, kappa) used for diffusive CFL => `max_diffusivity` + max_1_kappa::RealT # max(1, kappa) used for parabolic CFL => `max_diffusivity` equations_hyperbolic::E # CompressibleEulerEquations1D gradient_variables::GradientVariables # GradientVariablesPrimitive or GradientVariablesEntropy @@ -176,7 +176,7 @@ function flux(u, gradients, orientation::Integer, # by dispatching on the type of `equations.mu`. mu = dynamic_viscosity(u, equations) - # viscous flux components in the x-direction + # parabolic flux components in the x-direction f1 = 0 f2 = tau_11 * mu f3 = (v1 * tau_11 + q1) * mu @@ -242,7 +242,7 @@ function entropy2cons(w, equations::CompressibleNavierStokesDiffusion1D) end # the `flux` function takes in transformed variables `u` which depend on the type of the gradient variables. -# For CNS, it is simplest to formulate the viscous terms in primitive variables, so we transform the transformed +# For CNS, it is simplest to formulate the parabolic terms in primitive variables, so we transform the transformed # variables into primitive variables. @inline function convert_transformed_to_primitive(u_transformed, equations::CompressibleNavierStokesDiffusion1D{GradientVariablesPrimitive}) @@ -260,7 +260,7 @@ end # reverse engineers the gradients to be terms of the primitive variables (v1, T). # Helpful because then the diffusive fluxes have the same form as on paper. # Note, the first component of `gradient_entropy_vars` contains gradient(rho) which is unused. -# TODO: parabolic; entropy stable viscous terms +# TODO: parabolic; entropy stable parabolic terms @inline function convert_derivative_to_primitive(u, gradient, ::CompressibleNavierStokesDiffusion1D{GradientVariablesPrimitive}) return gradient diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index 0deb1599742..2e7625838da 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -96,7 +96,7 @@ struct CompressibleNavierStokesDiffusion2D{GradientVariables, RealT <: Real, Mu, mu::Mu # viscosity Pr::RealT # Prandtl number kappa::RealT # thermal diffusivity for Fick's law - max_4over3_kappa::RealT # max(4/3, kappa) used for diffusive CFL => `max_diffusivity` + max_4over3_kappa::RealT # max(4/3, kappa) used for parabolic CFL => `max_diffusivity` equations_hyperbolic::E # CompressibleEulerEquations2D gradient_variables::GradientVariables # GradientVariablesPrimitive or GradientVariablesEntropy @@ -185,7 +185,7 @@ function flux(u, gradients, orientation::Integer, mu = dynamic_viscosity(u, equations) if orientation == 1 - # viscous flux components in the x-direction + # parabolic flux components in the x-direction f1 = 0 f2 = tau_11 * mu f3 = tau_12 * mu @@ -193,7 +193,7 @@ function flux(u, gradients, orientation::Integer, return SVector(f1, f2, f3, f4) else # if orientation == 2 - # viscous flux components in the y-direction + # parabolic flux components in the y-direction # Note, symmetry is exploited for tau_12 = tau_21 g1 = 0 g2 = tau_12 * mu # tau_21 * mu @@ -263,7 +263,7 @@ function entropy2cons(w, equations::CompressibleNavierStokesDiffusion2D) end # the `flux` function takes in transformed variables `u` which depend on the type of the gradient variables. -# For CNS, it is simplest to formulate the viscous terms in primitive variables, so we transform the transformed +# For CNS, it is simplest to formulate the parabolic terms in primitive variables, so we transform the transformed # variables into primitive variables. @inline function convert_transformed_to_primitive(u_transformed, equations::CompressibleNavierStokesDiffusion2D{GradientVariablesPrimitive}) @@ -281,7 +281,7 @@ end # reverse engineers the gradients to be terms of the primitive variables (v1, v2, T). # Helpful because then the diffusive fluxes have the same form as on paper. # Note, the first component of `gradient_entropy_vars` contains gradient(rho) which is unused. -# TODO: parabolic; entropy stable viscous terms +# TODO: parabolic; entropy stable parabolic terms @inline function convert_derivative_to_primitive(u, gradient, ::CompressibleNavierStokesDiffusion2D{GradientVariablesPrimitive}) return gradient @@ -579,7 +579,7 @@ end x, t, operator_type::Divergence, equations::CompressibleNavierStokesDiffusion2D{GradientVariablesPrimitive}) - # for Dirichlet boundary conditions, we do not impose any conditions on the viscous fluxes + # for Dirichlet boundary conditions, we do not impose any conditions on the parabolic fluxes return flux_inner end end # @muladd diff --git a/src/equations/compressible_navier_stokes_3d.jl b/src/equations/compressible_navier_stokes_3d.jl index b24f8467033..874408e1dc8 100644 --- a/src/equations/compressible_navier_stokes_3d.jl +++ b/src/equations/compressible_navier_stokes_3d.jl @@ -96,7 +96,7 @@ struct CompressibleNavierStokesDiffusion3D{GradientVariables, RealT <: Real, Mu, mu::Mu # viscosity Pr::RealT # Prandtl number kappa::RealT # thermal diffusivity for Fick's law - max_4over3_kappa::RealT # max(4/3, kappa) used for diffusive CFL => `max_diffusivity` + max_4over3_kappa::RealT # max(4/3, kappa) used for parabolic CFL => `max_diffusivity` equations_hyperbolic::E # CompressibleEulerEquations3D gradient_variables::GradientVariables # GradientVariablesPrimitive or GradientVariablesEntropy @@ -198,7 +198,7 @@ function flux(u, gradients, orientation::Integer, mu = dynamic_viscosity(u, equations) if orientation == 1 - # viscous flux components in the x-direction + # parabolic flux components in the x-direction f1 = 0 f2 = tau_11 * mu f3 = tau_12 * mu @@ -207,7 +207,7 @@ function flux(u, gradients, orientation::Integer, return SVector(f1, f2, f3, f4, f5) elseif orientation == 2 - # viscous flux components in the y-direction + # parabolic flux components in the y-direction # Note, symmetry is exploited for tau_12 = tau_21 g1 = 0 g2 = tau_12 * mu # tau_21 * mu @@ -217,7 +217,7 @@ function flux(u, gradients, orientation::Integer, return SVector(g1, g2, g3, g4, g5) else # if orientation == 3 - # viscous flux components in the z-direction + # parabolic flux components in the z-direction # Note, symmetry is exploited for tau_13 = tau_31, tau_23 = tau_32 h1 = 0 h2 = tau_13 * mu # tau_31 * mu @@ -289,7 +289,7 @@ function entropy2cons(w, equations::CompressibleNavierStokesDiffusion3D) end # the `flux` function takes in transformed variables `u` which depend on the type of the gradient variables. -# For CNS, it is simplest to formulate the viscous terms in primitive variables, so we transform the transformed +# For CNS, it is simplest to formulate the parabolic terms in primitive variables, so we transform the transformed # variables into primitive variables. @inline function convert_transformed_to_primitive(u_transformed, equations::CompressibleNavierStokesDiffusion3D{GradientVariablesPrimitive}) @@ -307,7 +307,7 @@ end # reverse engineers the gradients to be terms of the primitive variables (v1, v2, v3, T). # Helpful because then the diffusive fluxes have the same form as on paper. # Note, the first component of `gradient_entropy_vars` contains gradient(rho) which is unused. -# TODO: parabolic; entropy stable viscous terms +# TODO: parabolic; entropy stable parabolic terms @inline function convert_derivative_to_primitive(u, gradient, ::CompressibleNavierStokesDiffusion3D{GradientVariablesPrimitive}) return gradient @@ -619,7 +619,7 @@ end x, t, operator_type::Divergence, equations::CompressibleNavierStokesDiffusion3D{GradientVariablesPrimitive}) - # for Dirichlet boundary conditions, we do not impose any conditions on the viscous fluxes + # for Dirichlet boundary conditions, we do not impose any conditions on the parabolic fluxes return flux_inner end end # @muladd diff --git a/src/equations/equations_parabolic.jl b/src/equations/equations_parabolic.jl index bfd33cfae78..716aba0ef2a 100644 --- a/src/equations/equations_parabolic.jl +++ b/src/equations/equations_parabolic.jl @@ -17,7 +17,7 @@ abstract type AbstractLaplaceDiffusion{NDIMS, NVARS} <: # Returns - `True()` -Used in diffusive CFL condition computation (see [`StepsizeCallback`](@ref)) to indicate that the +Used in parabolic CFL condition computation (see [`StepsizeCallback`](@ref)) to indicate that the diffusivity is constant in space and that [`max_diffusivity`](@ref) needs **not** to be re-computed at every node in every element. @@ -32,7 +32,7 @@ if the diffusion term is linear in the variables/constant. # Returns - `equations_parabolic.diffusivity` -Returns isotropic diffusion coefficient for use in diffusive CFL condition computation, +Returns isotropic diffusion coefficient for use in parabolic CFL condition computation, see [`StepsizeCallback`](@ref). """ @inline function max_diffusivity(equations_parabolic::AbstractLaplaceDiffusion) diff --git a/src/equations/laplace_diffusion_2d.jl b/src/equations/laplace_diffusion_2d.jl index 3741116f1bb..c243f3de364 100644 --- a/src/equations/laplace_diffusion_2d.jl +++ b/src/equations/laplace_diffusion_2d.jl @@ -32,7 +32,7 @@ end # The penalization depends on the solver, but also depends explicitly on physical parameters, # and would probably need to be specialized for every different equation. function penalty(u_outer, u_inner, inv_h, equations_parabolic::LaplaceDiffusion2D, - dg::ViscousFormulationLocalDG) + dg::ParabolicFormulationLocalDG) return dg.penalty_parameter * (u_outer - u_inner) * equations_parabolic.diffusivity end diff --git a/src/equations/laplace_diffusion_3d.jl b/src/equations/laplace_diffusion_3d.jl index c2e3c49afee..ec3957676a2 100644 --- a/src/equations/laplace_diffusion_3d.jl +++ b/src/equations/laplace_diffusion_3d.jl @@ -35,7 +35,7 @@ end # The penalization depends on the solver, but also depends explicitly on physical parameters, # and would probably need to be specialized for every different equation. function penalty(u_outer, u_inner, inv_h, equations_parabolic::LaplaceDiffusion3D, - dg::ViscousFormulationLocalDG) + dg::ParabolicFormulationLocalDG) return dg.penalty_parameter * (u_outer - u_inner) * equations_parabolic.diffusivity end diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 01315331bdd..dc86a8bac84 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -357,11 +357,12 @@ end ### StepsizeCallback ################################################################################ # In case of coupled system, use minimum timestep over all systems -function calculate_dt(u_ode, t, cfl_advective, cfl_diffusive, +function calculate_dt(u_ode, t, cfl_hyperbolic, cfl_parabolic, semi::SemidiscretizationCoupled) dt = minimum(eachsystem(semi)) do i u_ode_slice = get_system_u_ode(u_ode, i, semi) - return calculate_dt(u_ode_slice, t, cfl_advective, cfl_diffusive, semi.semis[i]) + return calculate_dt(u_ode_slice, t, cfl_hyperbolic, cfl_parabolic, + semi.semis[i]) end return dt diff --git a/src/semidiscretization/semidiscretization_coupled_p4est.jl b/src/semidiscretization/semidiscretization_coupled_p4est.jl index 18fc21c8f50..9de383d5b1c 100644 --- a/src/semidiscretization/semidiscretization_coupled_p4est.jl +++ b/src/semidiscretization/semidiscretization_coupled_p4est.jl @@ -323,11 +323,11 @@ end # In case of coupled system, use minimum timestep over all systems # Case for constant `cfl_number`. -function calculate_dt(u_ode, t, cfl_advective, cfl_diffusive, +function calculate_dt(u_ode, t, cfl_hyperbolic, cfl_parabolic, semi::SemidiscretizationCoupledP4est) dt = minimum(eachsystem(semi)) do i u_ode_slice = get_system_u_ode(u_ode, i, semi) - calculate_dt(u_ode_slice, t, cfl_advective, cfl_diffusive, semi.semis[i]) + calculate_dt(u_ode_slice, t, cfl_hyperbolic, cfl_parabolic, semi.semis[i]) end return dt diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index c730439017c..bb6487d46f4 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -79,7 +79,7 @@ function SemidiscretizationHyperbolicParabolic(mesh, equations::Tuple, @assert ndims(mesh) == ndims(equations_parabolic) if !(nvariables(equations) == nvariables(equations_parabolic)) - throw(ArgumentError("Current implementation of viscous terms requires the same number of conservative and gradient variables.")) + throw(ArgumentError("Current implementation of parabolic terms requires the same number of conservative and gradient variables.")) end boundary_conditions, boundary_conditions_parabolic = boundary_conditions diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index cb500094ca4..507f91c7890 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -30,7 +30,7 @@ function create_cache_parabolic(mesh::DGMultiMesh, (dg.basis.Nq, mesh.md.num_elements)), ndims(mesh))) - flux_viscous = similar.(gradients) + flux_parabolic = similar.(gradients) u_face_values = allocate_nested_array(uEltype, nvars, size(md.xf), dg) scalar_flux_face_values = similar(u_face_values) @@ -38,20 +38,20 @@ function create_cache_parabolic(mesh::DGMultiMesh, local_u_values_threaded = [similar(u_transformed, dg.basis.Nq) for _ in 1:Threads.maxthreadid()] - local_flux_viscous_threaded = [SVector{ndims(mesh)}(ntuple(_ -> similar(u_transformed, - dg.basis.Nq), - ndims(mesh))) - for _ in 1:Threads.maxthreadid()] + local_flux_parabolic_threaded = [SVector{ndims(mesh)}(ntuple(_ -> similar(u_transformed, + dg.basis.Nq), + ndims(mesh))) + for _ in 1:Threads.maxthreadid()] local_flux_face_values_threaded = [similar(scalar_flux_face_values[:, 1]) for _ in 1:Threads.maxthreadid()] - return (; u_transformed, gradients, flux_viscous, + return (; u_transformed, gradients, flux_parabolic, weak_differentiation_matrices, strong_differentiation_matrices, gradient_lift_matrix, projection_face_interpolation_matrix, divergence_lift_matrix, dxidxhatj, J, invJ, # geometric terms u_face_values, gradients_face_values, scalar_flux_face_values, - local_u_values_threaded, local_flux_viscous_threaded, + local_u_values_threaded, local_flux_parabolic_threaded, local_flux_face_values_threaded) end @@ -115,13 +115,13 @@ end function calc_volume_integral_gradient!(gradients, u, mesh::DGMultiMesh{NDIMS, <:NonAffine}, equations::AbstractEquationsParabolic, dg::DGMulti, cache, cache_parabolic) where {NDIMS} - (; strong_differentiation_matrices, dxidxhatj, local_flux_viscous_threaded) = cache_parabolic + (; strong_differentiation_matrices, dxidxhatj, local_flux_parabolic_threaded) = cache_parabolic # compute volume contributions to gradients @threaded for e in eachelement(mesh, dg) # compute gradients with respect to reference coordinates - local_reference_gradients = local_flux_viscous_threaded[Threads.threadid()] + local_reference_gradients = local_flux_parabolic_threaded[Threads.threadid()] for i in eachdim(mesh) apply_to_each_field(mul_by!(strong_differentiation_matrices[i]), local_reference_gradients[i], view(u, :, e)) @@ -143,7 +143,7 @@ end function calc_interface_flux_gradient!(scalar_flux_face_values, mesh::DGMultiMesh, equations, dg::DGMulti, - parabolic_scheme::ViscousFormulationBassiRebay1, + parabolic_scheme::ParabolicFormulationBassiRebay1, cache, cache_parabolic) (; u_face_values) = cache_parabolic (; mapM, mapP) = mesh.md @@ -298,11 +298,11 @@ function calc_single_boundary_flux!(flux_face_values, u_face_values, t, return nothing end -function calc_viscous_fluxes!(flux_viscous, u, gradients, mesh::DGMultiMesh, - equations::AbstractEquationsParabolic, - dg::DGMulti, cache, cache_parabolic) +function calc_parabolic_fluxes!(flux_parabolic, u, gradients, mesh::DGMultiMesh, + equations::AbstractEquationsParabolic, + dg::DGMulti, cache, cache_parabolic) for dim in eachdim(mesh) - set_zero!(flux_viscous[dim], dg) + set_zero!(flux_parabolic[dim], dg) end (; local_u_values_threaded) = cache_parabolic @@ -315,13 +315,13 @@ function calc_viscous_fluxes!(flux_viscous, u, gradients, mesh::DGMultiMesh, fill!(local_u_values, zero(eltype(local_u_values))) apply_to_each_field(mul_by!(dg.basis.Vq), local_u_values, view(u, :, e)) - # compute viscous flux at quad points + # compute parabolic flux at quad points for i in eachindex(local_u_values) u_i = local_u_values[i] gradients_i = getindex.(gradients, i, e) for dim in eachdim(mesh) - flux_viscous_i = flux(u_i, gradients_i, dim, equations) - setindex!(flux_viscous[dim], flux_viscous_i, i, e) + flux_parabolic_i = flux(u_i, gradients_i, dim, equations) + setindex!(flux_parabolic[dim], flux_parabolic_i, i, e) end end end @@ -330,18 +330,19 @@ function calc_viscous_fluxes!(flux_viscous, u, gradients, mesh::DGMultiMesh, end # no penalization for a BR1 parabolic solver -function calc_viscous_penalty!(scalar_flux_face_values, u_face_values, t, - boundary_conditions, - mesh, equations::AbstractEquationsParabolic, - dg::DGMulti, parabolic_scheme::ViscousFormulationBassiRebay1, - cache, cache_parabolic) +function calc_parabolic_penalty!(scalar_flux_face_values, u_face_values, t, + boundary_conditions, + mesh, equations::AbstractEquationsParabolic, + dg::DGMulti, + parabolic_scheme::ParabolicFormulationBassiRebay1, + cache, cache_parabolic) return nothing end -function calc_viscous_penalty!(scalar_flux_face_values, u_face_values, t, - boundary_conditions, mesh, - equations::AbstractEquationsParabolic, - dg::DGMulti, parabolic_scheme, cache, cache_parabolic) +function calc_parabolic_penalty!(scalar_flux_face_values, u_face_values, t, + boundary_conditions, mesh, + equations::AbstractEquationsParabolic, + dg::DGMulti, parabolic_scheme, cache, cache_parabolic) # compute fluxes at interfaces (; scalar_flux_face_values) = cache_parabolic (; mapM, mapP) = mesh.md @@ -354,7 +355,7 @@ function calc_viscous_penalty!(scalar_flux_face_values, u_face_values, t, return nothing end -function calc_volume_integral_divergence!(du, u, flux_viscous, mesh::DGMultiMesh, +function calc_volume_integral_divergence!(du, u, flux_parabolic, mesh::DGMultiMesh, equations::AbstractEquationsParabolic, dg::DGMulti, cache, cache_parabolic) (; weak_differentiation_matrices) = cache_parabolic @@ -364,36 +365,36 @@ function calc_volume_integral_divergence!(du, u, flux_viscous, mesh::DGMultiMesh for i in eachdim(mesh), j in eachdim(mesh) dxidxhatj = mesh.md.rstxyzJ[i, j][1, e] # assumes mesh is affine apply_to_each_field(mul_by_accum!(weak_differentiation_matrices[j], dxidxhatj), - view(du, :, e), view(flux_viscous[i], :, e)) + view(du, :, e), view(flux_parabolic[i], :, e)) end end return nothing end -function calc_volume_integral_divergence!(du, u, flux_viscous, +function calc_volume_integral_divergence!(du, u, flux_parabolic, mesh::DGMultiMesh{NDIMS, <:NonAffine}, equations::AbstractEquationsParabolic, dg::DGMulti, cache, cache_parabolic) where {NDIMS} - (; weak_differentiation_matrices, dxidxhatj, local_flux_viscous_threaded) = cache_parabolic + (; weak_differentiation_matrices, dxidxhatj, local_flux_parabolic_threaded) = cache_parabolic # compute volume contributions to divergence @threaded for e in eachelement(mesh, dg) - local_viscous_flux = local_flux_viscous_threaded[Threads.threadid()][1] + local_parabolic_flux = local_flux_parabolic_threaded[Threads.threadid()][1] for i in eachdim(mesh) # rotate flux to reference coordinates - fill!(local_viscous_flux, zero(eltype(local_viscous_flux))) + fill!(local_parabolic_flux, zero(eltype(local_parabolic_flux))) for j in eachdim(mesh) - for node in eachindex(local_viscous_flux) - local_viscous_flux[node] = local_viscous_flux[node] + - dxidxhatj[j, i][node, e] * - flux_viscous[j][node, e] + for node in eachindex(local_parabolic_flux) + local_parabolic_flux[node] = local_parabolic_flux[node] + + dxidxhatj[j, i][node, e] * + flux_parabolic[j][node, e] end end # differentiate with respect to reference coordinates apply_to_each_field(mul_by_accum!(weak_differentiation_matrices[i]), - view(du, :, e), local_viscous_flux) + view(du, :, e), local_parabolic_flux) end end @@ -402,9 +403,9 @@ end function calc_interface_flux_divergence!(scalar_flux_face_values, mesh, equations, dg, - parabolic_scheme::ViscousFormulationBassiRebay1, + parabolic_scheme::ParabolicFormulationBassiRebay1, cache, cache_parabolic) - flux_viscous_face_values = cache_parabolic.gradients_face_values # reuse storage + flux_parabolic_face_values = cache_parabolic.gradients_face_values # reuse storage (; mapM, mapP, nxyzJ) = mesh.md @threaded for face_node_index in each_face_node_global(mesh, dg, cache, cache_parabolic) @@ -413,8 +414,8 @@ function calc_interface_flux_divergence!(scalar_flux_face_values, # compute f(u, ∇u) ⋅ n flux_face_value = zero(eltype(scalar_flux_face_values)) for dim in eachdim(mesh) - fM = flux_viscous_face_values[dim][idM] - fP = flux_viscous_face_values[dim][idP] + fM = flux_parabolic_face_values[dim][idM] + fP = flux_parabolic_face_values[dim][idP] # Here, we use the "weak" formulation to compute the divergence (to ensure stability on curved meshes). flux_face_value = flux_face_value + 0.5f0 * (fP + fM) * nxyzJ[dim][face_node_index] @@ -425,21 +426,21 @@ function calc_interface_flux_divergence!(scalar_flux_face_values, return nothing end -function calc_divergence!(du, u::StructArray, t, flux_viscous, mesh::DGMultiMesh, +function calc_divergence!(du, u::StructArray, t, flux_parabolic, mesh::DGMultiMesh, equations::AbstractEquationsParabolic, boundary_conditions, dg::DGMulti, parabolic_scheme, cache, cache_parabolic) set_zero!(du, dg) - calc_volume_integral_divergence!(du, u, flux_viscous, mesh, equations, dg, cache, + calc_volume_integral_divergence!(du, u, flux_parabolic, mesh, equations, dg, cache, cache_parabolic) # interpolates from solution coefficients to face quadrature points (; projection_face_interpolation_matrix) = cache_parabolic - flux_viscous_face_values = cache_parabolic.gradients_face_values # reuse storage + flux_parabolic_face_values = cache_parabolic.gradients_face_values # reuse storage for dim in eachdim(mesh) apply_to_each_field(mul_by!(projection_face_interpolation_matrix), - flux_viscous_face_values[dim], flux_viscous[dim]) + flux_parabolic_face_values[dim], flux_parabolic[dim]) end # compute fluxes at interfaces @@ -452,9 +453,9 @@ function calc_divergence!(du, u::StructArray, t, flux_viscous, mesh::DGMultiMesh Divergence(), boundary_conditions, mesh, equations, dg, cache, cache_parabolic) - calc_viscous_penalty!(scalar_flux_face_values, cache_parabolic.u_face_values, t, - boundary_conditions, mesh, equations, dg, parabolic_scheme, - cache, cache_parabolic) + calc_parabolic_penalty!(scalar_flux_face_values, cache_parabolic.u_face_values, t, + boundary_conditions, mesh, equations, dg, parabolic_scheme, + cache, cache_parabolic) # surface contributions apply_to_each_field(mul_by_accum!(cache_parabolic.divergence_lift_matrix), du, @@ -476,7 +477,7 @@ function rhs_parabolic!(du, u, t, mesh::DGMultiMesh, set_zero!(du, dg) @trixi_timeit timer() "transform variables" begin - (; u_transformed, gradients, flux_viscous) = cache_parabolic + (; u_transformed, gradients, flux_parabolic) = cache_parabolic transform_variables!(u_transformed, u, mesh, equations_parabolic, dg, cache) end @@ -486,13 +487,13 @@ function rhs_parabolic!(du, u, t, mesh::DGMultiMesh, boundary_conditions, dg, parabolic_scheme, cache, cache_parabolic) end - @trixi_timeit timer() "calc viscous fluxes" begin - calc_viscous_fluxes!(flux_viscous, u_transformed, gradients, - mesh, equations_parabolic, dg, cache, cache_parabolic) + @trixi_timeit timer() "calc parabolic fluxes" begin + calc_parabolic_fluxes!(flux_parabolic, u_transformed, gradients, + mesh, equations_parabolic, dg, cache, cache_parabolic) end @trixi_timeit timer() "calc divergence" begin - calc_divergence!(du, u_transformed, t, flux_viscous, mesh, equations_parabolic, + calc_divergence!(du, u_transformed, t, flux_parabolic, mesh, equations_parabolic, boundary_conditions, dg, parabolic_scheme, cache, cache_parabolic) end @@ -500,7 +501,7 @@ function rhs_parabolic!(du, u, t, mesh::DGMultiMesh, # Note: we do not flip the sign of the geometric Jacobian here. # This is because the parabolic fluxes are assumed to be of the form # `du/dt + df/dx = dg/dx + source(x,t)`, - # where f(u) is the inviscid flux and g(u) is the viscous flux. + # where f(u) is the inviscid flux and g(u) is the parabolic flux. invert_jacobian!(du, mesh, equations_parabolic, dg, cache; scaling = 1) end diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl index d8565b6d22a..21153d2deca 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl @@ -10,7 +10,7 @@ Reusing `rhs_parabolic!` for `P4estMesh`es is not easily possible as for `P4estMesh`es we call ``` - prolong2mortars_divergence!(cache, flux_viscous, mesh, equations_parabolic, + prolong2mortars_divergence!(cache, flux_parabolic, mesh, equations_parabolic, dg.mortar, dg) calc_mortar_flux_divergence!(cache_parabolic.elements.surface_flux_values, @@ -18,7 +18,7 @@ for `P4estMesh`es we call ``` instead of ``` - prolong2mortars!(cache, flux_viscous, mesh, equations_parabolic, + prolong2mortars!(cache, flux_parabolic, mesh, equations_parabolic, dg.mortar, dg) calc_mortar_flux!(cache_parabolic.elements.surface_flux_values, @@ -30,10 +30,10 @@ function rhs_parabolic!(du, u, t, mesh::Union{P4estMesh{2}, P4estMesh{3}}, equations_parabolic::AbstractEquationsParabolic, boundary_conditions_parabolic, source_terms_parabolic, dg::DG, parabolic_scheme, cache, cache_parabolic) - @unpack viscous_container = cache_parabolic - @unpack u_transformed, gradients, flux_viscous = viscous_container + @unpack parabolic_container = cache_parabolic + @unpack u_transformed, gradients, flux_parabolic = parabolic_container - # Convert conservative variables to a form more suitable for viscous flux calculations + # Convert conservative variables to a form more suitable for parabolic flux calculations @trixi_timeit timer() "transform variables" begin transform_variables!(u_transformed, u, mesh, equations_parabolic, dg, cache) @@ -46,20 +46,20 @@ function rhs_parabolic!(du, u, t, mesh::Union{P4estMesh{2}, P4estMesh{3}}, dg, parabolic_scheme, cache) end - # Compute and store the viscous fluxes - @trixi_timeit timer() "calculate viscous fluxes" begin - calc_viscous_fluxes!(flux_viscous, gradients, u_transformed, mesh, - equations_parabolic, dg, cache) + # Compute and store the parabolic fluxes + @trixi_timeit timer() "calculate parabolic fluxes" begin + calc_parabolic_fluxes!(flux_parabolic, gradients, u_transformed, mesh, + equations_parabolic, dg, cache) end # The remainder of this function is essentially a regular rhs! for parabolic - # equations (i.e., it computes the divergence of the viscous fluxes) + # equations (i.e., it computes the divergence of the parabolic fluxes) # - # OBS! In `calc_viscous_fluxes!`, the viscous flux values at the volume nodes of each element have - # been computed and stored in `fluxes_viscous`. In the following, we *reuse* (abuse) the + # OBS! In `calc_parabolic_fluxes!`, the parabolic flux values at the volume nodes of each element have + # been computed and stored in `flux_parabolic`. In the following, we *reuse* (abuse) the # `interfaces` and `boundaries` containers in `cache` to interpolate and store the # *fluxes* at the element surfaces, as opposed to interpolating and storing the *solution* (as it - # is done in the hyperbolic operator). That is, `interfaces.u`/`boundaries.u` store *viscous flux values* + # is done in the hyperbolic operator). That is, `interfaces.u`/`boundaries.u` store *parabolic flux values* # and *not the solution*. The advantage is that a) we do not need to allocate more storage, b) we # do not need to recreate the existing data structure only with a different name, and c) we do not # need to interpolate solutions *and* gradients to the surfaces. @@ -68,32 +68,32 @@ function rhs_parabolic!(du, u, t, mesh::Union{P4estMesh{2}, P4estMesh{3}}, @trixi_timeit timer() "reset ∂u/∂t" set_zero!(du, dg, cache) # Calculate volume integral - # This calls the specialized version for the viscous fluxes from + # This calls the specialized version for the parabolic fluxes from # `dg_2d_parabolic.jl` or `dg_3d_parabolic.jl`. @trixi_timeit timer() "volume integral" begin - calc_volume_integral!(du, flux_viscous, mesh, equations_parabolic, dg, cache) + calc_volume_integral!(du, flux_parabolic, mesh, equations_parabolic, dg, cache) end # Prolong solution to interfaces. - # This calls the specialized version for the viscous fluxes from + # This calls the specialized version for the parabolic fluxes from # `dg_2d_parabolic.jl` or `dg_3d_parabolic.jl`. @trixi_timeit timer() "prolong2interfaces" begin - prolong2interfaces!(cache, flux_viscous, mesh, equations_parabolic, dg) + prolong2interfaces!(cache, flux_parabolic, mesh, equations_parabolic, dg) end # Calculate interface fluxes - # This calls the specialized version for the viscous fluxes from + # This calls the specialized version for the parabolic fluxes from # `dg_2d_parabolic.jl` or `dg_3d_parabolic.jl`. @trixi_timeit timer() "interface flux" begin calc_interface_flux!(cache.elements.surface_flux_values, mesh, equations_parabolic, dg, parabolic_scheme, cache) end - # Prolong viscous fluxes to boundaries. - # This calls the specialized version for the viscous fluxes from + # Prolong parabolic fluxes to boundaries. + # This calls the specialized version for the parabolic fluxes from # `dg_2d_parabolic.jl` or `dg_3d_parabolic.jl`. @trixi_timeit timer() "prolong2boundaries" begin - prolong2boundaries!(cache, flux_viscous, mesh, equations_parabolic, dg) + prolong2boundaries!(cache, flux_parabolic, mesh, equations_parabolic, dg) end # Calculate boundary fluxes. @@ -108,7 +108,7 @@ function rhs_parabolic!(du, u, t, mesh::Union{P4estMesh{2}, P4estMesh{3}}, # Prolong solution to mortars. # This calls the specialized version for parabolic equations. @trixi_timeit timer() "prolong2mortars" begin - prolong2mortars_divergence!(cache, flux_viscous, mesh, equations_parabolic, + prolong2mortars_divergence!(cache, flux_parabolic, mesh, equations_parabolic, dg.mortar, dg) end @@ -337,7 +337,7 @@ function calc_interface_flux_gradient!(surface_flux_values, return nothing end -# This is the version used when calculating the gradient of the viscous fluxes (called from above) +# This is the version used when calculating the gradient of the parabolic fluxes (called from above) @inline function calc_interface_flux_gradient!(surface_flux_values, mesh::P4estMesh{2}, equations_parabolic, dg::DG, parabolic_scheme, cache, @@ -366,22 +366,22 @@ end return nothing end -# This is the version used when calculating the divergence of the viscous fluxes. +# This is the version used when calculating the divergence of the parabolic fluxes. # Identical to weak-form volume integral/kernel for the purely hyperbolic case, -# except that the fluxes are here already precomputed in `calc_viscous_fluxes!` -function calc_volume_integral!(du, flux_viscous, mesh::P4estMesh{2}, +# except that the fluxes are here already precomputed in `calc_parabolic_fluxes!` +function calc_volume_integral!(du, flux_parabolic, mesh::P4estMesh{2}, equations_parabolic::AbstractEquationsParabolic, dg::DGSEM, cache) (; derivative_hat) = dg.basis (; contravariant_vectors) = cache.elements - flux_viscous_x, flux_viscous_y = flux_viscous + flux_parabolic_x, flux_parabolic_y = flux_parabolic @threaded for element in eachelement(dg, cache) # Calculate volume terms in one element for j in eachnode(dg), i in eachnode(dg) - flux1 = get_node_vars(flux_viscous_x, equations_parabolic, dg, + flux1 = get_node_vars(flux_parabolic_x, equations_parabolic, dg, i, j, element) - flux2 = get_node_vars(flux_viscous_y, equations_parabolic, dg, + flux2 = get_node_vars(flux_parabolic_y, equations_parabolic, dg, i, j, element) # Compute the contravariant flux by taking the scalar product of the @@ -411,18 +411,18 @@ function calc_volume_integral!(du, flux_viscous, mesh::P4estMesh{2}, return nothing end -# This is the version used when calculating the divergence of the viscous fluxes. -# Specialization `flux_viscous::Tuple` needed to +# This is the version used when calculating the divergence of the parabolic fluxes. +# Specialization `flux_parabolic::Tuple` needed to # avoid amibiguity with the hyperbolic version of `prolong2interfaces!` in dg_2d.jl # which is for the variables itself, i.e., `u::Array{uEltype, 4}`. -function prolong2interfaces!(cache, flux_viscous::Tuple, +function prolong2interfaces!(cache, flux_parabolic::Tuple, mesh::Union{P4estMesh{2}, P4estMeshView{2}}, equations_parabolic::AbstractEquationsParabolic, dg::DG) (; interfaces) = cache (; contravariant_vectors) = cache.elements index_range = eachnode(dg) - flux_viscous_x, flux_viscous_y = flux_viscous + flux_parabolic_x, flux_parabolic_y = flux_parabolic @threaded for interface in eachinterface(dg, cache) # Copy solution data from the primary element using "delayed indexing" with @@ -450,12 +450,12 @@ function prolong2interfaces!(cache, flux_viscous::Tuple, for v in eachvariable(equations_parabolic) # OBS! `interfaces.u` stores the interpolated *fluxes* and *not the solution*! - flux_viscous = SVector(flux_viscous_x[v, i_primary, j_primary, - primary_element], - flux_viscous_y[v, i_primary, j_primary, - primary_element]) + flux_parabolic = SVector(flux_parabolic_x[v, i_primary, j_primary, + primary_element], + flux_parabolic_y[v, i_primary, j_primary, + primary_element]) - interfaces.u[1, v, i, interface] = dot(flux_viscous, normal_direction) + interfaces.u[1, v, i, interface] = dot(flux_parabolic, normal_direction) end i_primary += i_primary_step j_primary += j_primary_step @@ -485,13 +485,14 @@ function prolong2interfaces!(cache, flux_viscous::Tuple, for v in eachvariable(equations_parabolic) # OBS! `interfaces.u` stores the interpolated *fluxes* and *not the solution*! - flux_viscous = SVector(flux_viscous_x[v, i_secondary, j_secondary, - secondary_element], - flux_viscous_y[v, i_secondary, j_secondary, - secondary_element]) + flux_parabolic = SVector(flux_parabolic_x[v, i_secondary, j_secondary, + secondary_element], + flux_parabolic_y[v, i_secondary, j_secondary, + secondary_element]) # store the normal flux with respect to the primary normal direction, # which is the negative of the secondary normal direction - interfaces.u[2, v, i, interface] = -dot(flux_viscous, normal_direction) + interfaces.u[2, v, i, interface] = -dot(flux_parabolic, + normal_direction) end i_secondary += i_secondary_step j_secondary += j_secondary_step @@ -547,15 +548,15 @@ function calc_interface_flux!(surface_flux_values, mesh::P4estMesh{2}, i_primary, j_primary, primary_element) - # We prolong the viscous flux dotted with respect the outward normal on the + # We prolong the parabolic flux dotted with respect the outward normal on the # primary element. - viscous_flux_normal_ll, viscous_flux_normal_rr = get_surface_node_vars(cache.interfaces.u, - equations_parabolic, - dg, - i, - interface) + parabolic_flux_normal_ll, parabolic_flux_normal_rr = get_surface_node_vars(cache.interfaces.u, + equations_parabolic, + dg, + i, + interface) - flux_ = flux_parabolic(viscous_flux_normal_ll, viscous_flux_normal_rr, + flux_ = flux_parabolic(parabolic_flux_normal_ll, parabolic_flux_normal_rr, normal_direction, Divergence(), equations_parabolic, parabolic_scheme) @@ -577,7 +578,7 @@ function calc_interface_flux!(surface_flux_values, mesh::P4estMesh{2}, return nothing end -function prolong2mortars_divergence!(cache, flux_viscous, +function prolong2mortars_divergence!(cache, flux_parabolic, mesh::P4estMesh{2}, equations_parabolic, mortar_l2::LobattoLegendreMortarL2, dg::DGSEM) @@ -585,7 +586,7 @@ function prolong2mortars_divergence!(cache, flux_viscous, @unpack contravariant_vectors = cache.elements index_range = eachnode(dg) - flux_viscous_x, flux_viscous_y = flux_viscous + flux_parabolic_x, flux_parabolic_y = flux_parabolic @threaded for mortar in eachmortar(dg, cache) # Copy solution data from the small elements using "delayed indexing" with @@ -608,10 +609,12 @@ function prolong2mortars_divergence!(cache, flux_viscous, i_small, j_small, element) for v in eachvariable(equations_parabolic) - flux_viscous = SVector(flux_viscous_x[v, i_small, j_small, element], - flux_viscous_y[v, i_small, j_small, element]) + flux_parabolic = SVector(flux_parabolic_x[v, i_small, j_small, + element], + flux_parabolic_y[v, i_small, j_small, + element]) - cache.mortars.u[1, v, position, i, mortar] = dot(flux_viscous, + cache.mortars.u[1, v, position, i, mortar] = dot(flux_parabolic, normal_direction) end i_small += i_small_step @@ -642,14 +645,14 @@ function prolong2mortars_divergence!(cache, flux_viscous, i_large, j_large, element) for v in eachvariable(equations_parabolic) - flux_viscous = SVector(flux_viscous_x[v, i_large, j_large, element], - flux_viscous_y[v, i_large, j_large, element]) + flux_parabolic = SVector(flux_parabolic_x[v, i_large, j_large, element], + flux_parabolic_y[v, i_large, j_large, element]) - # We prolong the viscous flux dotted with respect the outward normal + # We prolong the parabolic flux dotted with respect the outward normal # on the small element. We scale by -1/2 here because the normal # direction on the large element is negative 2x that of the small # element (these normal directions are "scaled" by the surface Jacobian) - u_buffer[v, i] = -0.5f0 * dot(flux_viscous, normal_direction) + u_buffer[v, i] = -0.5f0 * dot(flux_parabolic, normal_direction) end i_large += i_large_step j_large += j_large_step @@ -699,11 +702,13 @@ function calc_mortar_flux_divergence!(surface_flux_values, mesh::P4estMesh{2}, i_small, j_small, element) for v in eachvariable(equations_parabolic) - viscous_flux_normal_ll = cache.mortars.u[1, v, position, i, mortar] - viscous_flux_normal_rr = cache.mortars.u[2, v, position, i, mortar] + parabolic_flux_normal_ll = cache.mortars.u[1, v, position, i, + mortar] + parabolic_flux_normal_rr = cache.mortars.u[2, v, position, i, + mortar] - flux_ = flux_parabolic(viscous_flux_normal_ll, - viscous_flux_normal_rr, + flux_ = flux_parabolic(parabolic_flux_normal_ll, + parabolic_flux_normal_rr, normal_direction, Divergence(), equations_parabolic, parabolic_scheme) @@ -792,7 +797,7 @@ end # We structure `calc_mortar_flux_gradient!` similarly to "calc_mortar_flux!" for # hyperbolic equations with no nonconservative terms. # The reasoning is that parabolic fluxes are treated like conservative -# terms (e.g., we compute a viscous conservative "flux") and thus no +# terms (e.g., we compute a parabolic conservative "flux") and thus no # non-conservative terms are present. @inline function calc_mortar_flux_gradient!(fstar_primary, fstar_secondary, mesh::P4estMesh{2}, equations_parabolic, @@ -818,10 +823,10 @@ end return nothing end -# Specialization `flux_viscous::Tuple` needed to +# Specialization `flux_parabolic::Tuple` needed to # avoid amibiguity with the hyperbolic version of `prolong2boundaries!` in dg_2d.jl # which is for the variables itself, i.e., `u::Array{uEltype, 4}`. -function prolong2boundaries!(cache, flux_viscous::Tuple, +function prolong2boundaries!(cache, flux_parabolic::Tuple, mesh::P4estMesh{2}, equations_parabolic::AbstractEquationsParabolic, dg::DG) @@ -829,7 +834,7 @@ function prolong2boundaries!(cache, flux_viscous::Tuple, (; contravariant_vectors) = cache.elements index_range = eachnode(dg) - flux_viscous_x, flux_viscous_y = flux_viscous + flux_parabolic_x, flux_parabolic_y = flux_parabolic @threaded for boundary in eachboundary(dg, cache) # Copy solution data from the element using "delayed indexing" with @@ -849,10 +854,10 @@ function prolong2boundaries!(cache, flux_viscous::Tuple, i_node, j_node, element) for v in eachvariable(equations_parabolic) - flux_viscous = SVector(flux_viscous_x[v, i_node, j_node, element], - flux_viscous_y[v, i_node, j_node, element]) + flux_parabolic = SVector(flux_parabolic_x[v, i_node, j_node, element], + flux_parabolic_y[v, i_node, j_node, element]) - boundaries.u[v, i, boundary] = dot(flux_viscous, normal_direction) + boundaries.u[v, i, boundary] = dot(flux_parabolic, normal_direction) end i_node += i_node_step j_node += j_node_step @@ -1190,7 +1195,7 @@ end # Needed to *not* flip the sign of the inverse Jacobian. # This is because the parabolic fluxes are assumed to be of the form # `du/dt + df/dx = dg/dx + source(x,t)`, -# where f(u) is the inviscid flux and g(u) is the viscous flux. +# where f(u) is the inviscid flux and g(u) is the parabolic flux. function apply_jacobian_parabolic!(du::AbstractArray, mesh::P4estMesh{2}, equations_parabolic::AbstractEquationsParabolic, dg::DG, cache) diff --git a/src/solvers/dgsem_p4est/dg_3d_parabolic.jl b/src/solvers/dgsem_p4est/dg_3d_parabolic.jl index 635fb9b72bc..4dd650e5c4f 100644 --- a/src/solvers/dgsem_p4est/dg_3d_parabolic.jl +++ b/src/solvers/dgsem_p4est/dg_3d_parabolic.jl @@ -195,24 +195,24 @@ end return nothing end -# This is the version used when calculating the divergence of the viscous fluxes. +# This is the version used when calculating the divergence of the parabolic fluxes. # Identical to weak-form volume integral/kernel for the purely hyperbolic case, -# except that the fluxes are here already precomputed in `calc_viscous_fluxes!` -function calc_volume_integral!(du, flux_viscous, mesh::P4estMesh{3}, +# except that the fluxes are here already precomputed in `calc_parabolic_fluxes!` +function calc_volume_integral!(du, flux_parabolic, mesh::P4estMesh{3}, equations_parabolic::AbstractEquationsParabolic, dg::DGSEM, cache) (; derivative_hat) = dg.basis (; contravariant_vectors) = cache.elements - flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous + flux_parabolic_x, flux_parabolic_y, flux_parabolic_z = flux_parabolic @threaded for element in eachelement(dg, cache) # Calculate volume terms in one element for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) - flux1 = get_node_vars(flux_viscous_x, equations_parabolic, dg, + flux1 = get_node_vars(flux_parabolic_x, equations_parabolic, dg, i, j, k, element) - flux2 = get_node_vars(flux_viscous_y, equations_parabolic, dg, + flux2 = get_node_vars(flux_parabolic_y, equations_parabolic, dg, i, j, k, element) - flux3 = get_node_vars(flux_viscous_z, equations_parabolic, dg, + flux3 = get_node_vars(flux_parabolic_z, equations_parabolic, dg, i, j, k, element) # Compute the contravariant flux by taking the scalar product of the @@ -253,18 +253,18 @@ function calc_volume_integral!(du, flux_viscous, mesh::P4estMesh{3}, return nothing end -# This is the version used when calculating the divergence of the viscous fluxes. -# Specialization `flux_viscous::Tuple` needed to +# This is the version used when calculating the divergence of the parabolic fluxes. +# Specialization `flux_parabolic::Tuple` needed to # avoid amibiguity with the hyperbolic version of `prolong2interfaces!` in dg_3d.jl # which is for the variables itself, i.e., `u::Array{uEltype, 5}`. -function prolong2interfaces!(cache, flux_viscous::Tuple, +function prolong2interfaces!(cache, flux_parabolic::Tuple, mesh::P4estMesh{3}, equations_parabolic::AbstractEquationsParabolic, dg::DG) (; interfaces) = cache (; contravariant_vectors) = cache.elements index_range = eachnode(dg) - flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous + flux_parabolic_x, flux_parabolic_y, flux_parabolic_z = flux_parabolic @threaded for interface in eachinterface(dg, cache) # Copy solution data from the primary element using "delayed indexing" with @@ -298,14 +298,23 @@ function prolong2interfaces!(cache, flux_viscous::Tuple, for v in eachvariable(equations_parabolic) # OBS! `interfaces.u` stores the interpolated *fluxes* and *not the solution*! - flux_viscous = SVector(flux_viscous_x[v, i_primary, j_primary, - k_primary, primary_element], - flux_viscous_y[v, i_primary, j_primary, - k_primary, primary_element], - flux_viscous_z[v, i_primary, j_primary, - k_primary, primary_element]) - - interfaces.u[1, v, i, j, interface] = dot(flux_viscous, + flux_parabolic = SVector(flux_parabolic_x[v, + i_primary, + j_primary, + k_primary, + primary_element], + flux_parabolic_y[v, + i_primary, + j_primary, + k_primary, + primary_element], + flux_parabolic_z[v, + i_primary, + j_primary, + k_primary, + primary_element]) + + interfaces.u[1, v, i, j, interface] = dot(flux_parabolic, normal_direction) end i_primary += i_primary_step_i @@ -346,18 +355,24 @@ function prolong2interfaces!(cache, flux_viscous::Tuple, for v in eachvariable(equations_parabolic) # OBS! `interfaces.u` stores the interpolated *fluxes* and *not the solution*! - flux_viscous = SVector(flux_viscous_x[v, i_secondary, j_secondary, - k_secondary, - secondary_element], - flux_viscous_y[v, i_secondary, j_secondary, - k_secondary, - secondary_element], - flux_viscous_z[v, i_secondary, j_secondary, - k_secondary, - secondary_element]) + flux_parabolic = SVector(flux_parabolic_x[v, + i_secondary, + j_secondary, + k_secondary, + secondary_element], + flux_parabolic_y[v, + i_secondary, + j_secondary, + k_secondary, + secondary_element], + flux_parabolic_z[v, + i_secondary, + j_secondary, + k_secondary, + secondary_element]) # store the normal flux with respect to the primary normal direction, # which is the negative of the secondary normal direction - interfaces.u[2, v, i, j, interface] = -dot(flux_viscous, + interfaces.u[2, v, i, j, interface] = -dot(flux_parabolic, normal_direction) end i_secondary += i_secondary_step_i @@ -423,16 +438,17 @@ function calc_interface_flux!(surface_flux_values, mesh::P4estMesh{3}, contravariant_vectors, i_primary, j_primary, k_primary, primary_element) - # We prolong the viscous flux dotted with respect the outward normal on the + # We prolong the parabolic flux dotted with respect the outward normal on the # primary element. - viscous_flux_normal_ll, viscous_flux_normal_rr = get_surface_node_vars(cache.interfaces.u, - equations_parabolic, - dg, - i, - j, - interface) - - flux_ = flux_parabolic(viscous_flux_normal_ll, viscous_flux_normal_rr, + parabolic_flux_normal_ll, parabolic_flux_normal_rr = get_surface_node_vars(cache.interfaces.u, + equations_parabolic, + dg, + i, + j, + interface) + + flux_ = flux_parabolic(parabolic_flux_normal_ll, + parabolic_flux_normal_rr, normal_direction, Divergence(), equations_parabolic, parabolic_scheme) @@ -464,7 +480,7 @@ function calc_interface_flux!(surface_flux_values, mesh::P4estMesh{3}, return nothing end -function prolong2mortars_divergence!(cache, flux_viscous, +function prolong2mortars_divergence!(cache, flux_parabolic, mesh::P4estMesh{3}, equations_parabolic, mortar_l2::LobattoLegendreMortarL2, dg::DGSEM) @@ -473,7 +489,7 @@ function prolong2mortars_divergence!(cache, flux_viscous, @unpack contravariant_vectors = cache.elements index_range = eachnode(dg) - flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous + flux_parabolic_x, flux_parabolic_y, flux_parabolic_z = flux_parabolic @threaded for mortar in eachmortar(dg, cache) # Copy solution data from the small elements using "delayed indexing" with @@ -501,14 +517,14 @@ function prolong2mortars_divergence!(cache, flux_viscous, element) for v in eachvariable(equations_parabolic) - flux_viscous = SVector(flux_viscous_x[v, i_small, j_small, - k_small, element], - flux_viscous_y[v, i_small, j_small, - k_small, element], - flux_viscous_z[v, i_small, j_small, - k_small, element]) - - cache.mortars.u[1, v, position, i, j, mortar] = dot(flux_viscous, + flux_parabolic = SVector(flux_parabolic_x[v, i_small, j_small, + k_small, element], + flux_parabolic_y[v, i_small, j_small, + k_small, element], + flux_parabolic_z[v, i_small, j_small, + k_small, element]) + + cache.mortars.u[1, v, position, i, j, mortar] = dot(flux_parabolic, normal_direction) end i_small += i_small_step_i @@ -551,18 +567,27 @@ function prolong2mortars_divergence!(cache, flux_viscous, element) for v in eachvariable(equations_parabolic) - flux_viscous = SVector(flux_viscous_x[v, i_large, j_large, k_large, - element], - flux_viscous_y[v, i_large, j_large, k_large, - element], - flux_viscous_z[v, i_large, j_large, k_large, - element]) - - # We prolong the viscous flux dotted with respect the outward normal + flux_parabolic = SVector(flux_parabolic_x[v, + i_large, + j_large, + k_large, + element], + flux_parabolic_y[v, + i_large, + j_large, + k_large, + element], + flux_parabolic_z[v, + i_large, + j_large, + k_large, + element]) + + # We prolong the parabolic flux dotted with respect the outward normal # on the small element. We scale by -1/2 here because the normal # direction on the large element is negative 2x that of the small # element (these normal directions are "scaled" by the surface Jacobian) - u_buffer[v, i, j] = -0.5f0 * dot(flux_viscous, normal_direction) + u_buffer[v, i, j] = -0.5f0 * dot(flux_parabolic, normal_direction) end i_large += i_large_step_i j_large += j_large_step_i @@ -634,13 +659,13 @@ function calc_mortar_flux_divergence!(surface_flux_values, element) for v in eachvariable(equations_parabolic) - viscous_flux_normal_ll = cache.mortars.u[1, v, position, i, j, - mortar] - viscous_flux_normal_rr = cache.mortars.u[2, v, position, i, j, - mortar] + parabolic_flux_normal_ll = cache.mortars.u[1, v, position, i, j, + mortar] + parabolic_flux_normal_rr = cache.mortars.u[2, v, position, i, j, + mortar] - flux_ = flux_parabolic(viscous_flux_normal_ll, - viscous_flux_normal_rr, + flux_ = flux_parabolic(parabolic_flux_normal_ll, + parabolic_flux_normal_rr, normal_direction, Divergence(), equations_parabolic, parabolic_scheme) @@ -748,7 +773,7 @@ end # We structure `calc_mortar_flux_gradient!` similarly to "calc_mortar_flux!" for # hyperbolic equations with no nonconservative terms. # The reasoning is that parabolic fluxes are treated like conservative -# terms (e.g., we compute a viscous conservative "flux") and thus no +# terms (e.g., we compute a parabolic conservative "flux") and thus no # non-conservative terms are present. @inline function calc_mortar_flux_gradient!(fstar_primary, fstar_secondary, mesh::P4estMesh{3}, @@ -852,10 +877,10 @@ function calc_volume_integral_gradient!(gradients, u_transformed, return nothing end -# Specialization `flux_viscous::Tuple` needed to +# Specialization `flux_parabolic::Tuple` needed to # avoid amibiguity with the hyperbolic version of `prolong2boundaries!` in dg_3d.jl # which is for the variables itself, i.e., `u::Array{uEltype, 5}`. -function prolong2boundaries!(cache, flux_viscous::Tuple, +function prolong2boundaries!(cache, flux_parabolic::Tuple, mesh::P4estMesh{3}, equations_parabolic::AbstractEquationsParabolic, dg::DG) @@ -863,7 +888,7 @@ function prolong2boundaries!(cache, flux_viscous::Tuple, (; contravariant_vectors) = cache.elements index_range = eachnode(dg) - flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous + flux_parabolic_x, flux_parabolic_y, flux_parabolic_z = flux_parabolic @threaded for boundary in eachboundary(dg, cache) # Copy solution data from the element using "delayed indexing" with @@ -891,14 +916,14 @@ function prolong2boundaries!(cache, flux_viscous::Tuple, i_node, j_node, k_node, element) for v in eachvariable(equations_parabolic) - flux_viscous = SVector(flux_viscous_x[v, i_node, j_node, k_node, - element], - flux_viscous_y[v, i_node, j_node, k_node, - element], - flux_viscous_z[v, i_node, j_node, k_node, - element]) - - boundaries.u[v, i, j, boundary] = dot(flux_viscous, + flux_parabolic = SVector(flux_parabolic_x[v, i_node, j_node, k_node, + element], + flux_parabolic_y[v, i_node, j_node, k_node, + element], + flux_parabolic_z[v, i_node, j_node, k_node, + element]) + + boundaries.u[v, i, j, boundary] = dot(flux_parabolic, normal_direction) end i_node += i_node_step_i @@ -1092,7 +1117,7 @@ end # Needed to *not* flip the sign of the inverse Jacobian. # This is because the parabolic fluxes are assumed to be of the form # `du/dt + df/dx = dg/dx + source(x,t)`, -# where f(u) is the inviscid flux and g(u) is the viscous flux. +# where f(u) is the inviscid flux and g(u) is the parabolic flux. function apply_jacobian_parabolic!(du::AbstractArray, mesh::P4estMesh{3}, equations_parabolic::AbstractEquationsParabolic, dg::DG, cache) diff --git a/src/solvers/dgsem_tree/container_parabolic_1d.jl b/src/solvers/dgsem_tree/container_parabolic_1d.jl new file mode 100644 index 00000000000..af9a6841636 --- /dev/null +++ b/src/solvers/dgsem_tree/container_parabolic_1d.jl @@ -0,0 +1,58 @@ +mutable struct ParabolicContainer1D{uEltype <: Real} + u_transformed::Array{uEltype, 3} + gradients::Array{uEltype, 3} + flux_parabolic::Array{uEltype, 3} + + # internal `resize!`able storage + _u_transformed::Vector{uEltype} + _gradients::Vector{uEltype} + _flux_parabolic::Vector{uEltype} + + function ParabolicContainer1D{uEltype}(n_vars::Integer, n_nodes::Integer, + n_elements::Integer) where {uEltype <: Real} + return new(Array{uEltype, 3}(undef, n_vars, n_nodes, n_elements), + Array{uEltype, 3}(undef, n_vars, n_nodes, n_elements), + Array{uEltype, 3}(undef, n_vars, n_nodes, n_elements), + Vector{uEltype}(undef, n_vars * n_nodes * n_elements), + Vector{uEltype}(undef, n_vars * n_nodes * n_elements), + Vector{uEltype}(undef, n_vars * n_nodes * n_elements)) + end +end + +function init_parabolic_container_1d(n_vars::Integer, n_nodes::Integer, + n_elements::Integer, + ::Type{uEltype}) where {uEltype <: Real} + return ParabolicContainer1D{uEltype}(n_vars, n_nodes, n_elements) +end + +# Only one-dimensional `Array`s are `resize!`able in Julia. +# Hence, we use `Vector`s as internal storage and `resize!` +# them whenever needed. Then, we reuse the same memory by +# `unsafe_wrap`ping multi-dimensional `Array`s around the +# internal storage. +function Base.resize!(parabolic_container::ParabolicContainer1D, equations, dg, cache) + capacity = nvariables(equations) * nnodes(dg) * nelements(dg, cache) + resize!(parabolic_container._u_transformed, capacity) + resize!(parabolic_container._gradients, capacity) + resize!(parabolic_container._flux_parabolic, capacity) + + parabolic_container.u_transformed = unsafe_wrap(Array, + pointer(parabolic_container._u_transformed), + (nvariables(equations), + nnodes(dg), + nelements(dg, cache))) + + parabolic_container.gradients = unsafe_wrap(Array, + pointer(parabolic_container._gradients), + (nvariables(equations), + nnodes(dg), + nelements(dg, cache))) + + parabolic_container.flux_parabolic = unsafe_wrap(Array, + pointer(parabolic_container._flux_parabolic), + (nvariables(equations), + nnodes(dg), + nelements(dg, cache))) + + return nothing +end diff --git a/src/solvers/dgsem_tree/container_parabolic_2d.jl b/src/solvers/dgsem_tree/container_parabolic_2d.jl new file mode 100644 index 00000000000..7eee9c1302d --- /dev/null +++ b/src/solvers/dgsem_tree/container_parabolic_2d.jl @@ -0,0 +1,84 @@ +mutable struct ParabolicContainer2D{uEltype <: Real} + u_transformed::Array{uEltype, 4} + gradients::NTuple{2, Array{uEltype, 4}} + flux_parabolic::NTuple{2, Array{uEltype, 4}} + + # internal `resize!`able storage + _u_transformed::Vector{uEltype} + # Use Tuple for outer, fixed-size datastructure + _gradients::Tuple{Vector{uEltype}, Vector{uEltype}} + _flux_parabolic::Tuple{Vector{uEltype}, Vector{uEltype}} + + function ParabolicContainer2D{uEltype}(n_vars::Integer, n_nodes::Integer, + n_elements::Integer) where {uEltype <: Real} + return new(Array{uEltype, 4}(undef, n_vars, n_nodes, n_nodes, n_elements), # `u_transformed` + # `gradients` + (Array{uEltype, 4}(undef, n_vars, n_nodes, n_nodes, n_elements), + Array{uEltype, 4}(undef, n_vars, n_nodes, n_nodes, n_elements)), + # `flux_parabolic` + (Array{uEltype, 4}(undef, n_vars, n_nodes, n_nodes, n_elements), + Array{uEltype, 4}(undef, n_vars, n_nodes, n_nodes, n_elements)), + # `_u_transformed` + Vector{uEltype}(undef, n_vars * n_nodes^2 * n_elements), + # `_gradients` + (Vector{uEltype}(undef, n_vars * n_nodes^2 * n_elements), + Vector{uEltype}(undef, n_vars * n_nodes^2 * n_elements)), + # `_flux_parabolic` + (Vector{uEltype}(undef, n_vars * n_nodes^2 * n_elements), + Vector{uEltype}(undef, n_vars * n_nodes^2 * n_elements))) + end +end + +function init_parabolic_container_2d(n_vars::Integer, n_nodes::Integer, + n_elements::Integer, + ::Type{uEltype}) where {uEltype <: Real} + return ParabolicContainer2D{uEltype}(n_vars, n_nodes, n_elements) +end + +# Only one-dimensional `Array`s are `resize!`able in Julia. +# Hence, we use `Vector`s as internal storage and `resize!` +# them whenever needed. Then, we reuse the same memory by +# `unsafe_wrap`ping multi-dimensional `Array`s around the +# internal storage. +function Base.resize!(parabolic_container::ParabolicContainer2D, equations, dg, cache) + capacity = nvariables(equations) * nnodes(dg)^2 * nelements(dg, cache) + resize!(parabolic_container._u_transformed, capacity) + for dim in 1:2 + resize!(parabolic_container._gradients[dim], capacity) + resize!(parabolic_container._flux_parabolic[dim], capacity) + end + + parabolic_container.u_transformed = unsafe_wrap(Array, + pointer(parabolic_container._u_transformed), + (nvariables(equations), + nnodes(dg), nnodes(dg), + nelements(dg, cache))) + + gradients_1 = unsafe_wrap(Array, + pointer(parabolic_container._gradients[1]), + (nvariables(equations), + nnodes(dg), nnodes(dg), + nelements(dg, cache))) + gradients_2 = unsafe_wrap(Array, + pointer(parabolic_container._gradients[2]), + (nvariables(equations), + nnodes(dg), nnodes(dg), + nelements(dg, cache))) + + parabolic_container.gradients = (gradients_1, gradients_2) + + flux_parabolic_1 = unsafe_wrap(Array, + pointer(parabolic_container._flux_parabolic[1]), + (nvariables(equations), + nnodes(dg), nnodes(dg), + nelements(dg, cache))) + flux_parabolic_2 = unsafe_wrap(Array, + pointer(parabolic_container._flux_parabolic[2]), + (nvariables(equations), + nnodes(dg), nnodes(dg), + nelements(dg, cache))) + + parabolic_container.flux_parabolic = (flux_parabolic_1, flux_parabolic_2) + + return nothing +end diff --git a/src/solvers/dgsem_tree/container_parabolic_3d.jl b/src/solvers/dgsem_tree/container_parabolic_3d.jl new file mode 100644 index 00000000000..07a1403b438 --- /dev/null +++ b/src/solvers/dgsem_tree/container_parabolic_3d.jl @@ -0,0 +1,99 @@ +mutable struct ParabolicContainer3D{uEltype <: Real} + u_transformed::Array{uEltype, 5} + gradients::NTuple{3, Array{uEltype, 5}} + flux_parabolic::NTuple{3, Array{uEltype, 5}} + + # internal `resize!`able storage + _u_transformed::Vector{uEltype} + # Use Tuple for outer, fixed-size datastructure + _gradients::Tuple{Vector{uEltype}, Vector{uEltype}, Vector{uEltype}} + _flux_parabolic::Tuple{Vector{uEltype}, Vector{uEltype}, Vector{uEltype}} + + function ParabolicContainer3D{uEltype}(n_vars::Integer, n_nodes::Integer, + n_elements::Integer) where {uEltype <: Real} + return new(Array{uEltype, 5}(undef, n_vars, n_nodes, n_nodes, n_nodes, n_elements), # `u_transformed` + # `gradients` + (Array{uEltype, 5}(undef, n_vars, n_nodes, n_nodes, n_nodes, n_elements), + Array{uEltype, 5}(undef, n_vars, n_nodes, n_nodes, n_nodes, n_elements), + Array{uEltype, 5}(undef, n_vars, n_nodes, n_nodes, n_nodes, n_elements)), + # `flux_parabolic` + (Array{uEltype, 5}(undef, n_vars, n_nodes, n_nodes, n_nodes, n_elements), + Array{uEltype, 5}(undef, n_vars, n_nodes, n_nodes, n_nodes, n_elements), + Array{uEltype, 5}(undef, n_vars, n_nodes, n_nodes, n_nodes, n_elements)), + # `u_transformed` + Vector{uEltype}(undef, n_vars * n_nodes^3 * n_elements), + # `_gradients` + (Vector{uEltype}(undef, n_vars * n_nodes^3 * n_elements), + Vector{uEltype}(undef, n_vars * n_nodes^3 * n_elements), + Vector{uEltype}(undef, n_vars * n_nodes^3 * n_elements)), + # `_flux_parabolic` + (Vector{uEltype}(undef, n_vars * n_nodes^3 * n_elements), + Vector{uEltype}(undef, n_vars * n_nodes^3 * n_elements), + Vector{uEltype}(undef, n_vars * n_nodes^3 * n_elements))) + end +end + +function init_parabolic_container_3d(n_vars::Integer, n_nodes::Integer, + n_elements::Integer, + ::Type{uEltype}) where {uEltype <: Real} + return ParabolicContainer3D{uEltype}(n_vars, n_nodes, n_elements) +end + +# Only one-dimensional `Array`s are `resize!`able in Julia. +# Hence, we use `Vector`s as internal storage and `resize!` +# them whenever needed. Then, we reuse the same memory by +# `unsafe_wrap`ping multi-dimensional `Array`s around the +# internal storage. +function Base.resize!(parabolic_container::ParabolicContainer3D, equations, dg, cache) + capacity = nvariables(equations) * nnodes(dg)^3 * nelements(dg, cache) + resize!(parabolic_container._u_transformed, capacity) + for dim in 1:3 + resize!(parabolic_container._gradients[dim], capacity) + resize!(parabolic_container._flux_parabolic[dim], capacity) + end + + parabolic_container.u_transformed = unsafe_wrap(Array, + pointer(parabolic_container._u_transformed), + (nvariables(equations), + nnodes(dg), nnodes(dg), nnodes(dg), + nelements(dg, cache))) + + gradients_1 = unsafe_wrap(Array, + pointer(parabolic_container._gradients[1]), + (nvariables(equations), + nnodes(dg), nnodes(dg), nnodes(dg), + nelements(dg, cache))) + gradients_2 = unsafe_wrap(Array, + pointer(parabolic_container._gradients[2]), + (nvariables(equations), + nnodes(dg), nnodes(dg), nnodes(dg), + nelements(dg, cache))) + gradients_3 = unsafe_wrap(Array, + pointer(parabolic_container._gradients[3]), + (nvariables(equations), + nnodes(dg), nnodes(dg), nnodes(dg), + nelements(dg, cache))) + + parabolic_container.gradients = (gradients_1, gradients_2, gradients_3) + + flux_parabolic_1 = unsafe_wrap(Array, + pointer(parabolic_container._flux_parabolic[1]), + (nvariables(equations), + nnodes(dg), nnodes(dg), nnodes(dg), + nelements(dg, cache))) + flux_parabolic_2 = unsafe_wrap(Array, + pointer(parabolic_container._flux_parabolic[2]), + (nvariables(equations), + nnodes(dg), nnodes(dg), nnodes(dg), + nelements(dg, cache))) + flux_parabolic_3 = unsafe_wrap(Array, + pointer(parabolic_container._flux_parabolic[3]), + (nvariables(equations), + nnodes(dg), nnodes(dg), nnodes(dg), + nelements(dg, cache))) + + parabolic_container.flux_parabolic = (flux_parabolic_1, flux_parabolic_2, + flux_parabolic_3) + + return nothing +end diff --git a/src/solvers/dgsem_tree/container_viscous_1d.jl b/src/solvers/dgsem_tree/container_viscous_1d.jl deleted file mode 100644 index 661fbfb237f..00000000000 --- a/src/solvers/dgsem_tree/container_viscous_1d.jl +++ /dev/null @@ -1,58 +0,0 @@ -mutable struct ViscousContainer1D{uEltype <: Real} - u_transformed::Array{uEltype, 3} - gradients::Array{uEltype, 3} - flux_viscous::Array{uEltype, 3} - - # internal `resize!`able storage - _u_transformed::Vector{uEltype} - _gradients::Vector{uEltype} - _flux_viscous::Vector{uEltype} - - function ViscousContainer1D{uEltype}(n_vars::Integer, n_nodes::Integer, - n_elements::Integer) where {uEltype <: Real} - return new(Array{uEltype, 3}(undef, n_vars, n_nodes, n_elements), - Array{uEltype, 3}(undef, n_vars, n_nodes, n_elements), - Array{uEltype, 3}(undef, n_vars, n_nodes, n_elements), - Vector{uEltype}(undef, n_vars * n_nodes * n_elements), - Vector{uEltype}(undef, n_vars * n_nodes * n_elements), - Vector{uEltype}(undef, n_vars * n_nodes * n_elements)) - end -end - -function init_viscous_container_1d(n_vars::Integer, n_nodes::Integer, - n_elements::Integer, - ::Type{uEltype}) where {uEltype <: Real} - return ViscousContainer1D{uEltype}(n_vars, n_nodes, n_elements) -end - -# Only one-dimensional `Array`s are `resize!`able in Julia. -# Hence, we use `Vector`s as internal storage and `resize!` -# them whenever needed. Then, we reuse the same memory by -# `unsafe_wrap`ping multi-dimensional `Array`s around the -# internal storage. -function Base.resize!(viscous_container::ViscousContainer1D, equations, dg, cache) - capacity = nvariables(equations) * nnodes(dg) * nelements(dg, cache) - resize!(viscous_container._u_transformed, capacity) - resize!(viscous_container._gradients, capacity) - resize!(viscous_container._flux_viscous, capacity) - - viscous_container.u_transformed = unsafe_wrap(Array, - pointer(viscous_container._u_transformed), - (nvariables(equations), - nnodes(dg), - nelements(dg, cache))) - - viscous_container.gradients = unsafe_wrap(Array, - pointer(viscous_container._gradients), - (nvariables(equations), - nnodes(dg), - nelements(dg, cache))) - - viscous_container.flux_viscous = unsafe_wrap(Array, - pointer(viscous_container._flux_viscous), - (nvariables(equations), - nnodes(dg), - nelements(dg, cache))) - - return nothing -end diff --git a/src/solvers/dgsem_tree/container_viscous_2d.jl b/src/solvers/dgsem_tree/container_viscous_2d.jl deleted file mode 100644 index a4e69643d3f..00000000000 --- a/src/solvers/dgsem_tree/container_viscous_2d.jl +++ /dev/null @@ -1,84 +0,0 @@ -mutable struct ViscousContainer2D{uEltype <: Real} - u_transformed::Array{uEltype, 4} - gradients::NTuple{2, Array{uEltype, 4}} - flux_viscous::NTuple{2, Array{uEltype, 4}} - - # internal `resize!`able storage - _u_transformed::Vector{uEltype} - # Use Tuple for outer, fixed-size datastructure - _gradients::Tuple{Vector{uEltype}, Vector{uEltype}} - _flux_viscous::Tuple{Vector{uEltype}, Vector{uEltype}} - - function ViscousContainer2D{uEltype}(n_vars::Integer, n_nodes::Integer, - n_elements::Integer) where {uEltype <: Real} - return new(Array{uEltype, 4}(undef, n_vars, n_nodes, n_nodes, n_elements), # `u_transformed` - # `gradients` - (Array{uEltype, 4}(undef, n_vars, n_nodes, n_nodes, n_elements), - Array{uEltype, 4}(undef, n_vars, n_nodes, n_nodes, n_elements)), - # `flux_viscous` - (Array{uEltype, 4}(undef, n_vars, n_nodes, n_nodes, n_elements), - Array{uEltype, 4}(undef, n_vars, n_nodes, n_nodes, n_elements)), - # `_u_transformed` - Vector{uEltype}(undef, n_vars * n_nodes^2 * n_elements), - # `_gradients` - (Vector{uEltype}(undef, n_vars * n_nodes^2 * n_elements), - Vector{uEltype}(undef, n_vars * n_nodes^2 * n_elements)), - # `_flux_viscous` - (Vector{uEltype}(undef, n_vars * n_nodes^2 * n_elements), - Vector{uEltype}(undef, n_vars * n_nodes^2 * n_elements))) - end -end - -function init_viscous_container_2d(n_vars::Integer, n_nodes::Integer, - n_elements::Integer, - ::Type{uEltype}) where {uEltype <: Real} - return ViscousContainer2D{uEltype}(n_vars, n_nodes, n_elements) -end - -# Only one-dimensional `Array`s are `resize!`able in Julia. -# Hence, we use `Vector`s as internal storage and `resize!` -# them whenever needed. Then, we reuse the same memory by -# `unsafe_wrap`ping multi-dimensional `Array`s around the -# internal storage. -function Base.resize!(viscous_container::ViscousContainer2D, equations, dg, cache) - capacity = nvariables(equations) * nnodes(dg)^2 * nelements(dg, cache) - resize!(viscous_container._u_transformed, capacity) - for dim in 1:2 - resize!(viscous_container._gradients[dim], capacity) - resize!(viscous_container._flux_viscous[dim], capacity) - end - - viscous_container.u_transformed = unsafe_wrap(Array, - pointer(viscous_container._u_transformed), - (nvariables(equations), - nnodes(dg), nnodes(dg), - nelements(dg, cache))) - - gradients_1 = unsafe_wrap(Array, - pointer(viscous_container._gradients[1]), - (nvariables(equations), - nnodes(dg), nnodes(dg), - nelements(dg, cache))) - gradients_2 = unsafe_wrap(Array, - pointer(viscous_container._gradients[2]), - (nvariables(equations), - nnodes(dg), nnodes(dg), - nelements(dg, cache))) - - viscous_container.gradients = (gradients_1, gradients_2) - - flux_viscous_1 = unsafe_wrap(Array, - pointer(viscous_container._flux_viscous[1]), - (nvariables(equations), - nnodes(dg), nnodes(dg), - nelements(dg, cache))) - flux_viscous_2 = unsafe_wrap(Array, - pointer(viscous_container._flux_viscous[2]), - (nvariables(equations), - nnodes(dg), nnodes(dg), - nelements(dg, cache))) - - viscous_container.flux_viscous = (flux_viscous_1, flux_viscous_2) - - return nothing -end diff --git a/src/solvers/dgsem_tree/container_viscous_3d.jl b/src/solvers/dgsem_tree/container_viscous_3d.jl deleted file mode 100644 index a55fc5147d9..00000000000 --- a/src/solvers/dgsem_tree/container_viscous_3d.jl +++ /dev/null @@ -1,98 +0,0 @@ -mutable struct ViscousContainer3D{uEltype <: Real} - u_transformed::Array{uEltype, 5} - gradients::NTuple{3, Array{uEltype, 5}} - flux_viscous::NTuple{3, Array{uEltype, 5}} - - # internal `resize!`able storage - _u_transformed::Vector{uEltype} - # Use Tuple for outer, fixed-size datastructure - _gradients::Tuple{Vector{uEltype}, Vector{uEltype}, Vector{uEltype}} - _flux_viscous::Tuple{Vector{uEltype}, Vector{uEltype}, Vector{uEltype}} - - function ViscousContainer3D{uEltype}(n_vars::Integer, n_nodes::Integer, - n_elements::Integer) where {uEltype <: Real} - return new(Array{uEltype, 5}(undef, n_vars, n_nodes, n_nodes, n_nodes, n_elements), # `u_transformed` - # `gradients` - (Array{uEltype, 5}(undef, n_vars, n_nodes, n_nodes, n_nodes, n_elements), - Array{uEltype, 5}(undef, n_vars, n_nodes, n_nodes, n_nodes, n_elements), - Array{uEltype, 5}(undef, n_vars, n_nodes, n_nodes, n_nodes, n_elements)), - # `flux_viscous` - (Array{uEltype, 5}(undef, n_vars, n_nodes, n_nodes, n_nodes, n_elements), - Array{uEltype, 5}(undef, n_vars, n_nodes, n_nodes, n_nodes, n_elements), - Array{uEltype, 5}(undef, n_vars, n_nodes, n_nodes, n_nodes, n_elements)), - # `u_transformed` - Vector{uEltype}(undef, n_vars * n_nodes^3 * n_elements), - # `_gradients` - (Vector{uEltype}(undef, n_vars * n_nodes^3 * n_elements), - Vector{uEltype}(undef, n_vars * n_nodes^3 * n_elements), - Vector{uEltype}(undef, n_vars * n_nodes^3 * n_elements)), - # `_flux_viscous` - (Vector{uEltype}(undef, n_vars * n_nodes^3 * n_elements), - Vector{uEltype}(undef, n_vars * n_nodes^3 * n_elements), - Vector{uEltype}(undef, n_vars * n_nodes^3 * n_elements))) - end -end - -function init_viscous_container_3d(n_vars::Integer, n_nodes::Integer, - n_elements::Integer, - ::Type{uEltype}) where {uEltype <: Real} - return ViscousContainer3D{uEltype}(n_vars, n_nodes, n_elements) -end - -# Only one-dimensional `Array`s are `resize!`able in Julia. -# Hence, we use `Vector`s as internal storage and `resize!` -# them whenever needed. Then, we reuse the same memory by -# `unsafe_wrap`ping multi-dimensional `Array`s around the -# internal storage. -function Base.resize!(viscous_container::ViscousContainer3D, equations, dg, cache) - capacity = nvariables(equations) * nnodes(dg)^3 * nelements(dg, cache) - resize!(viscous_container._u_transformed, capacity) - for dim in 1:3 - resize!(viscous_container._gradients[dim], capacity) - resize!(viscous_container._flux_viscous[dim], capacity) - end - - viscous_container.u_transformed = unsafe_wrap(Array, - pointer(viscous_container._u_transformed), - (nvariables(equations), - nnodes(dg), nnodes(dg), nnodes(dg), - nelements(dg, cache))) - - gradients_1 = unsafe_wrap(Array, - pointer(viscous_container._gradients[1]), - (nvariables(equations), - nnodes(dg), nnodes(dg), nnodes(dg), - nelements(dg, cache))) - gradients_2 = unsafe_wrap(Array, - pointer(viscous_container._gradients[2]), - (nvariables(equations), - nnodes(dg), nnodes(dg), nnodes(dg), - nelements(dg, cache))) - gradients_3 = unsafe_wrap(Array, - pointer(viscous_container._gradients[3]), - (nvariables(equations), - nnodes(dg), nnodes(dg), nnodes(dg), - nelements(dg, cache))) - - viscous_container.gradients = (gradients_1, gradients_2, gradients_3) - - flux_viscous_1 = unsafe_wrap(Array, - pointer(viscous_container._flux_viscous[1]), - (nvariables(equations), - nnodes(dg), nnodes(dg), nnodes(dg), - nelements(dg, cache))) - flux_viscous_2 = unsafe_wrap(Array, - pointer(viscous_container._flux_viscous[2]), - (nvariables(equations), - nnodes(dg), nnodes(dg), nnodes(dg), - nelements(dg, cache))) - flux_viscous_3 = unsafe_wrap(Array, - pointer(viscous_container._flux_viscous[3]), - (nvariables(equations), - nnodes(dg), nnodes(dg), nnodes(dg), - nelements(dg, cache))) - - viscous_container.flux_viscous = (flux_viscous_1, flux_viscous_2, flux_viscous_3) - - return nothing -end diff --git a/src/solvers/dgsem_tree/containers_parabolic.jl b/src/solvers/dgsem_tree/containers_parabolic.jl new file mode 100644 index 00000000000..f90eac9dcb0 --- /dev/null +++ b/src/solvers/dgsem_tree/containers_parabolic.jl @@ -0,0 +1,4 @@ +# Dimension-specific implementations +include("container_parabolic_1d.jl") +include("container_parabolic_2d.jl") +include("container_parabolic_3d.jl") diff --git a/src/solvers/dgsem_tree/containers_viscous.jl b/src/solvers/dgsem_tree/containers_viscous.jl deleted file mode 100644 index 444f2cb7303..00000000000 --- a/src/solvers/dgsem_tree/containers_viscous.jl +++ /dev/null @@ -1,4 +0,0 @@ -# Dimension-specific implementations -include("container_viscous_1d.jl") -include("container_viscous_2d.jl") -include("container_viscous_3d.jl") diff --git a/src/solvers/dgsem_tree/dg.jl b/src/solvers/dgsem_tree/dg.jl index c389af89880..5887367bc2b 100644 --- a/src/solvers/dgsem_tree/dg.jl +++ b/src/solvers/dgsem_tree/dg.jl @@ -35,7 +35,7 @@ include("containers.jl") include("dg_parallel.jl") # Helper structs for parabolic AMR -include("containers_viscous.jl") +include("containers_parabolic.jl") # Some functions for a second-order Finite-Volume (MUSCL) alike # scheme on DG-subcells. diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl index 79b649e7955..65383168835 100644 --- a/src/solvers/dgsem_tree/dg_1d.jl +++ b/src/solvers/dgsem_tree/dg_1d.jl @@ -385,8 +385,8 @@ end end # Used for both the purely hyperbolic conserved variables `u` -# and the viscous flux in x-direction in the 1D parabolic case. -function prolong2interfaces!(cache, u_or_flux_viscous, +# and the parabolic flux in x-direction in the 1D parabolic case. +function prolong2interfaces!(cache, u_or_flux_parabolic, mesh::TreeMesh{1}, equations, dg::DG) @unpack interfaces = cache @unpack neighbor_ids = interfaces @@ -398,16 +398,16 @@ function prolong2interfaces!(cache, u_or_flux_viscous, # interface in x-direction for v in eachvariable(equations) - interfaces_u[1, v, interface] = u_or_flux_viscous[v, nnodes(dg), - left_element] - interfaces_u[2, v, interface] = u_or_flux_viscous[v, 1, right_element] + interfaces_u[1, v, interface] = u_or_flux_parabolic[v, nnodes(dg), + left_element] + interfaces_u[2, v, interface] = u_or_flux_parabolic[v, 1, right_element] end end return nothing end -function prolong2interfaces!(cache, u_or_flux_viscous, +function prolong2interfaces!(cache, u_or_flux_parabolic, mesh::TreeMesh{1}, equations, dg::DGSEM{<:GaussLegendreBasis}) @unpack interfaces = cache @@ -430,11 +430,11 @@ function prolong2interfaces!(cache, u_or_flux_viscous, # (see comment at the top of the file) # Need `boundary_interpolation` at right (+1) node for left element interface_u_1 = (interface_u_1 + - u_or_flux_viscous[v, ii, left_element] * + u_or_flux_parabolic[v, ii, left_element] * boundary_interpolation[ii, 2]) # Need `boundary_interpolation` at left (-1) node for right element interface_u_2 = (interface_u_2 + - u_or_flux_viscous[v, ii, right_element] * + u_or_flux_parabolic[v, ii, right_element] * boundary_interpolation[ii, 1]) end interfaces_u[1, v, interface] = interface_u_1 @@ -519,8 +519,8 @@ function calc_interface_flux!(surface_flux_values, end # Used for both the purely hyperbolic conserved variables `u` -# and the viscous flux in x-direction in the 1D parabolic case. -function prolong2boundaries!(cache, u_or_flux_viscous, +# and the parabolic flux in x-direction in the 1D parabolic case. +function prolong2boundaries!(cache, u_or_flux_parabolic, mesh::TreeMesh{1}, equations, dg::DG) @unpack boundaries = cache @unpack neighbor_sides = boundaries @@ -532,11 +532,12 @@ function prolong2boundaries!(cache, u_or_flux_viscous, if neighbor_sides[boundary] == 1 # element in -x direction of boundary for v in eachvariable(equations) - boundaries.u[1, v, boundary] = u_or_flux_viscous[v, nnodes(dg), element] + boundaries.u[1, v, boundary] = u_or_flux_parabolic[v, nnodes(dg), + element] end else # Element in +x direction of boundary for v in eachvariable(equations) - boundaries.u[2, v, boundary] = u_or_flux_viscous[v, 1, element] + boundaries.u[2, v, boundary] = u_or_flux_parabolic[v, 1, element] end end end @@ -544,7 +545,7 @@ function prolong2boundaries!(cache, u_or_flux_viscous, return nothing end -function prolong2boundaries!(cache, u_or_flux_viscous, +function prolong2boundaries!(cache, u_or_flux_parabolic, mesh::TreeMesh{1}, equations, dg::DGSEM{<:GaussLegendreBasis}) @unpack boundaries = cache @@ -565,7 +566,7 @@ function prolong2boundaries!(cache, u_or_flux_viscous, # Not += to allow `@muladd` to turn these into FMAs # (see comment at the top of the file) boundary_u_1 = (boundary_u_1 + - u_or_flux_viscous[v, ii, element] * + u_or_flux_parabolic[v, ii, element] * boundary_interpolation[ii, 2]) end boundaries.u[1, v, boundary] = boundary_u_1 @@ -575,7 +576,7 @@ function prolong2boundaries!(cache, u_or_flux_viscous, boundary_u_2 = zero(eltype(boundaries.u)) for ii in eachnode(dg) boundary_u_2 = (boundary_u_2 + - u_or_flux_viscous[v, ii, element] * + u_or_flux_parabolic[v, ii, element] * boundary_interpolation[ii, 1]) end boundaries.u[2, v, boundary] = boundary_u_2 diff --git a/src/solvers/dgsem_tree/dg_1d_parabolic.jl b/src/solvers/dgsem_tree/dg_1d_parabolic.jl index d359469f823..f234b23eb70 100644 --- a/src/solvers/dgsem_tree/dg_1d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_1d_parabolic.jl @@ -11,11 +11,11 @@ function create_cache_parabolic(mesh::TreeMesh{1}, equations_hyperbolic::AbstractEquations, dg::DG, n_elements, uEltype) - viscous_container = init_viscous_container_1d(nvariables(equations_hyperbolic), - nnodes(dg), n_elements, - uEltype) + parabolic_container = init_parabolic_container_1d(nvariables(equations_hyperbolic), + nnodes(dg), n_elements, + uEltype) - cache_parabolic = (; viscous_container) + cache_parabolic = (; parabolic_container) return cache_parabolic end @@ -32,10 +32,10 @@ function rhs_parabolic!(du, u, t, mesh::TreeMesh{1}, equations_parabolic::AbstractEquationsParabolic, boundary_conditions_parabolic, source_terms_parabolic, dg::DG, parabolic_scheme, cache, cache_parabolic) - @unpack viscous_container = cache_parabolic - @unpack u_transformed, gradients, flux_viscous = viscous_container + @unpack parabolic_container = cache_parabolic + @unpack u_transformed, gradients, flux_parabolic = parabolic_container - # Convert conservative variables to a form more suitable for viscous flux calculations + # Convert conservative variables to a form more suitable for parabolic flux calculations @trixi_timeit timer() "transform variables" begin transform_variables!(u_transformed, u, mesh, equations_parabolic, dg, cache) @@ -48,20 +48,20 @@ function rhs_parabolic!(du, u, t, mesh::TreeMesh{1}, parabolic_scheme, cache) end - # Compute and store the viscous fluxes - @trixi_timeit timer() "calculate viscous fluxes" begin - calc_viscous_fluxes!(flux_viscous, gradients, u_transformed, mesh, - equations_parabolic, dg, cache) + # Compute and store the parabolic fluxes + @trixi_timeit timer() "calculate parabolic fluxes" begin + calc_parabolic_fluxes!(flux_parabolic, gradients, u_transformed, mesh, + equations_parabolic, dg, cache) end # The remainder of this function is essentially a regular rhs! for - # parabolic equations (i.e., it computes the divergence of the viscous fluxes) + # parabolic equations (i.e., it computes the divergence of the parabolic fluxes) # - # OBS! In `calc_viscous_fluxes!`, the viscous flux values at the volume nodes of each element have - # been computed and stored in `fluxes_viscous`. In the following, we *reuse* (abuse) the + # OBS! In `calc_parabolic_fluxes!`, the parabolic flux values at the volume nodes of each element have + # been computed and stored in `flux_parabolic`. In the following, we *reuse* (abuse) the # `interfaces` and `boundaries` containers in `cache` to interpolate and store the # *fluxes* at the element surfaces, as opposed to interpolating and storing the *solution* (as it - # is done in the hyperbolic operator). That is, `interfaces.u`/`boundaries.u` store *viscous flux values* + # is done in the hyperbolic operator). That is, `interfaces.u`/`boundaries.u` store *parabolic flux values* # and *not the solution*. The advantage is that a) we do not need to allocate more storage, b) we # do not need to recreate the existing data structure only with a different name, and c) we do not # need to interpolate solutions *and* gradients to the surfaces. @@ -70,19 +70,19 @@ function rhs_parabolic!(du, u, t, mesh::TreeMesh{1}, @trixi_timeit timer() "reset ∂u/∂t" set_zero!(du, dg, cache) # Calculate volume integral - # This calls the specialized version for the viscous flux. + # This calls the specialized version for the parabolic flux. @trixi_timeit timer() "volume integral" begin - calc_volume_integral!(du, flux_viscous, mesh, equations_parabolic, dg, cache) + calc_volume_integral!(du, flux_parabolic, mesh, equations_parabolic, dg, cache) end # Prolong solution to interfaces # This reuses `prolong2interfaces!` for the purely hyperbolic case. @trixi_timeit timer() "prolong2interfaces" begin - prolong2interfaces!(cache, flux_viscous, mesh, equations_parabolic, dg) + prolong2interfaces!(cache, flux_parabolic, mesh, equations_parabolic, dg) end # Calculate interface fluxes. - # This calls the specialized version for the viscous flux. + # This calls the specialized version for the parabolic flux. @trixi_timeit timer() "interface flux" begin calc_interface_flux!(cache.elements.surface_flux_values, mesh, equations_parabolic, dg, @@ -92,7 +92,7 @@ function rhs_parabolic!(du, u, t, mesh::TreeMesh{1}, # Prolong solution to boundaries. # This reuses `prolong2boundaries!` for the purely hyperbolic case. @trixi_timeit timer() "prolong2boundaries" begin - prolong2boundaries!(cache, flux_viscous, mesh, equations_parabolic, dg) + prolong2boundaries!(cache, flux_parabolic, mesh, equations_parabolic, dg) end # Calculate boundary fluxes. @@ -146,10 +146,10 @@ function transform_variables!(u_transformed, u, mesh::TreeMesh{1}, return nothing end -# This is the version used when calculating the divergence of the viscous fluxes. +# This is the version used when calculating the divergence of the parabolic fluxes. # Identical to weak-form volume integral/kernel for the purely hyperbolic case, -# except that the fluxes are here already precomputed in `calc_viscous_fluxes!` -function calc_volume_integral!(du, flux_viscous, mesh::TreeMesh{1}, +# except that the fluxes are here already precomputed in `calc_parabolic_fluxes!` +function calc_volume_integral!(du, flux_parabolic, mesh::TreeMesh{1}, equations_parabolic::AbstractEquationsParabolic, dg::DGSEM, cache) @unpack derivative_hat = dg.basis @@ -157,7 +157,7 @@ function calc_volume_integral!(du, flux_viscous, mesh::TreeMesh{1}, @threaded for element in eachelement(dg, cache) # Calculate volume terms in one element for i in eachnode(dg) - flux_1_node = get_node_vars(flux_viscous, equations_parabolic, dg, i, + flux_1_node = get_node_vars(flux_parabolic, equations_parabolic, dg, i, element) for ii in eachnode(dg) @@ -170,7 +170,7 @@ function calc_volume_integral!(du, flux_viscous, mesh::TreeMesh{1}, return nothing end -# This is the version used when calculating the divergence of the viscous fluxes +# This is the version used when calculating the divergence of the parabolic fluxes function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{1}, equations_parabolic, dg::DG, parabolic_scheme, cache) @@ -205,9 +205,10 @@ function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{1}, return nothing end -function calc_viscous_fluxes!(flux_viscous, gradients, u_transformed, mesh::TreeMesh{1}, - equations_parabolic::AbstractEquationsParabolic, - dg::DG, cache) +function calc_parabolic_fluxes!(flux_parabolic, gradients, u_transformed, + mesh::TreeMesh{1}, + equations_parabolic::AbstractEquationsParabolic, + dg::DG, cache) @threaded for element in eachelement(dg, cache) for i in eachnode(dg) # Get solution and gradients @@ -215,10 +216,10 @@ function calc_viscous_fluxes!(flux_viscous, gradients, u_transformed, mesh::Tree gradients_1_node = get_node_vars(gradients, equations_parabolic, dg, i, element) - # Calculate viscous flux and store each component for later use - flux_viscous_node = flux(u_node, (gradients_1_node,), 1, - equations_parabolic) - set_node_vars!(flux_viscous, flux_viscous_node, equations_parabolic, dg, + # Calculate parabolic flux and store each component for later use + flux_parabolic_node = flux(u_node, (gradients_1_node,), 1, + equations_parabolic) + set_node_vars!(flux_parabolic, flux_parabolic_node, equations_parabolic, dg, i, element) end end @@ -348,14 +349,14 @@ function calc_boundary_flux_by_direction_divergence!(surface_flux_values::Abstra @unpack surface_flux = surface_integral # Note: cache.boundaries.u contains the unsigned normal component (using "orientation", not "direction") - # of the viscous flux, as computed in `prolong2boundaries!` + # of the parabolic flux, as computed in `prolong2boundaries!` @unpack u, neighbor_ids, neighbor_sides, node_coordinates, orientations = cache.boundaries @threaded for boundary in first_boundary:last_boundary # Get neighboring element neighbor = neighbor_ids[boundary] - # Get viscous boundary fluxes + # Get parabolic boundary fluxes flux_ll, flux_rr = get_surface_node_vars(u, equations_parabolic, dg, boundary) if neighbor_sides[boundary] == 1 # Element is on the left, boundary on the right flux_inner = flux_ll @@ -566,7 +567,7 @@ end # Needed to *not* flip the sign of the inverse Jacobian. # This is because the parabolic fluxes are assumed to be of the form # `du/dt + df/dx = dg/dx + source(x,t)`, -# where f(u) is the inviscid flux and g(u) is the viscous flux. +# where f(u) is the inviscid flux and g(u) is the parabolic flux. function apply_jacobian_parabolic!(du::AbstractArray, mesh::TreeMesh{1}, equations_parabolic::AbstractEquationsParabolic, dg::DG, cache) diff --git a/src/solvers/dgsem_tree/dg_2d_parabolic.jl b/src/solvers/dgsem_tree/dg_2d_parabolic.jl index 53ecc51228a..bd2d4e5173b 100644 --- a/src/solvers/dgsem_tree/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_2d_parabolic.jl @@ -11,11 +11,11 @@ function create_cache_parabolic(mesh::Union{TreeMesh{2}, P4estMesh{2}}, equations_hyperbolic::AbstractEquations, dg::DG, n_elements, uEltype) - viscous_container = init_viscous_container_2d(nvariables(equations_hyperbolic), - nnodes(dg), n_elements, - uEltype) + parabolic_container = init_parabolic_container_2d(nvariables(equations_hyperbolic), + nnodes(dg), n_elements, + uEltype) - cache_parabolic = (; viscous_container) + cache_parabolic = (; parabolic_container) return cache_parabolic end @@ -32,10 +32,10 @@ function rhs_parabolic!(du, u, t, mesh::Union{TreeMesh{2}, TreeMesh{3}}, equations_parabolic::AbstractEquationsParabolic, boundary_conditions_parabolic, source_terms_parabolic, dg::DG, parabolic_scheme, cache, cache_parabolic) - @unpack viscous_container = cache_parabolic - @unpack u_transformed, gradients, flux_viscous = viscous_container + @unpack parabolic_container = cache_parabolic + @unpack u_transformed, gradients, flux_parabolic = parabolic_container - # Convert conservative variables to a form more suitable for viscous flux calculations + # Convert conservative variables to a form more suitable for parabolic flux calculations @trixi_timeit timer() "transform variables" begin transform_variables!(u_transformed, u, mesh, equations_parabolic, dg, cache) @@ -48,20 +48,20 @@ function rhs_parabolic!(du, u, t, mesh::Union{TreeMesh{2}, TreeMesh{3}}, dg, parabolic_scheme, cache) end - # Compute and store the viscous fluxes - @trixi_timeit timer() "calculate viscous fluxes" begin - calc_viscous_fluxes!(flux_viscous, gradients, u_transformed, mesh, - equations_parabolic, dg, cache) + # Compute and store the parabolic fluxes + @trixi_timeit timer() "calculate parabolic fluxes" begin + calc_parabolic_fluxes!(flux_parabolic, gradients, u_transformed, mesh, + equations_parabolic, dg, cache) end # The remainder of this function is essentially a regular rhs! for parabolic - # equations (i.e., it computes the divergence of the viscous fluxes) + # equations (i.e., it computes the divergence of the parabolic fluxes) # - # OBS! In `calc_viscous_fluxes!`, the viscous flux values at the volume nodes of each element have - # been computed and stored in `fluxes_viscous`. In the following, we *reuse* (abuse) the + # OBS! In `calc_parabolic_fluxes!`, the parabolic flux values at the volume nodes of each element have + # been computed and stored in `flux_parabolic`. In the following, we *reuse* (abuse) the # `interfaces` and `boundaries` containers in `cache` to interpolate and store the # *fluxes* at the element surfaces, as opposed to interpolating and storing the *solution* (as it - # is done in the hyperbolic operator). That is, `interfaces.u`/`boundaries.u` store *viscous flux values* + # is done in the hyperbolic operator). That is, `interfaces.u`/`boundaries.u` store *parabolic flux values* # and *not the solution*. The advantage is that a) we do not need to allocate more storage, b) we # do not need to recreate the existing data structure only with a different name, and c) we do not # need to interpolate solutions *and* gradients to the surfaces. @@ -70,21 +70,21 @@ function rhs_parabolic!(du, u, t, mesh::Union{TreeMesh{2}, TreeMesh{3}}, @trixi_timeit timer() "reset ∂u/∂t" set_zero!(du, dg, cache) # Calculate volume integral. - # This calls the specialized version for the viscous fluxes from + # This calls the specialized version for the parabolic fluxes from # `dg_2d_parabolic.jl` or `dg_3d_parabolic.jl`. @trixi_timeit timer() "volume integral" begin - calc_volume_integral!(du, flux_viscous, mesh, equations_parabolic, dg, cache) + calc_volume_integral!(du, flux_parabolic, mesh, equations_parabolic, dg, cache) end # Prolong solution to interfaces. - # This calls the specialized version for the viscous fluxes from + # This calls the specialized version for the parabolic fluxes from # `dg_2d_parabolic.jl` or `dg_3d_parabolic.jl`. @trixi_timeit timer() "prolong2interfaces" begin - prolong2interfaces!(cache, flux_viscous, mesh, equations_parabolic, dg) + prolong2interfaces!(cache, flux_parabolic, mesh, equations_parabolic, dg) end # Calculate interface fluxes - # This calls the specialized version for the viscous fluxes from + # This calls the specialized version for the parabolic fluxes from # `dg_2d_parabolic.jl` or `dg_3d_parabolic.jl`. @trixi_timeit timer() "interface flux" begin calc_interface_flux!(cache.elements.surface_flux_values, @@ -92,11 +92,11 @@ function rhs_parabolic!(du, u, t, mesh::Union{TreeMesh{2}, TreeMesh{3}}, parabolic_scheme, cache) end - # Prolong viscous fluxes to boundaries. - # This calls the specialized version for the viscous fluxes from + # Prolong parabolic fluxes to boundaries. + # This calls the specialized version for the parabolic fluxes from # `dg_2d_parabolic.jl` or `dg_3d_parabolic.jl`. @trixi_timeit timer() "prolong2boundaries" begin - prolong2boundaries!(cache, flux_viscous, mesh, equations_parabolic, dg) + prolong2boundaries!(cache, flux_parabolic, mesh, equations_parabolic, dg) end # Calculate boundary fluxes. @@ -108,11 +108,11 @@ function rhs_parabolic!(du, u, t, mesh::Union{TreeMesh{2}, TreeMesh{3}}, dg.surface_integral, dg) end - # Prolong viscous fluxes to mortars. - # This calls the specialized version for the viscous fluxes from + # Prolong parabolic fluxes to mortars. + # This calls the specialized version for the parabolic fluxes from # `dg_2d_parabolic.jl` or `dg_3d_parabolic.jl`. @trixi_timeit timer() "prolong2mortars" begin - prolong2mortars!(cache, flux_viscous, mesh, equations_parabolic, + prolong2mortars!(cache, flux_parabolic, mesh, equations_parabolic, dg.mortar, dg) end @@ -166,21 +166,21 @@ function transform_variables!(u_transformed, u, mesh::Union{TreeMesh{2}, P4estMe return nothing end -# This is the version used when calculating the divergence of the viscous fluxes. +# This is the version used when calculating the divergence of the parabolic fluxes. # Identical to weak-form volume integral/kernel for the purely hyperbolic case, -# except that the fluxes are here already precomputed in `calc_viscous_fluxes!` -function calc_volume_integral!(du, flux_viscous, mesh::TreeMesh{2}, +# except that the fluxes are here already precomputed in `calc_parabolic_fluxes!` +function calc_volume_integral!(du, flux_parabolic, mesh::TreeMesh{2}, equations_parabolic::AbstractEquationsParabolic, dg::DGSEM, cache) @unpack derivative_hat = dg.basis - flux_viscous_x, flux_viscous_y = flux_viscous + flux_parabolic_x, flux_parabolic_y = flux_parabolic @threaded for element in eachelement(dg, cache) # Calculate volume terms in one element for j in eachnode(dg), i in eachnode(dg) - flux_1_node = get_node_vars(flux_viscous_x, equations_parabolic, dg, + flux_1_node = get_node_vars(flux_parabolic_x, equations_parabolic, dg, i, j, element) - flux_2_node = get_node_vars(flux_viscous_y, equations_parabolic, dg, + flux_2_node = get_node_vars(flux_parabolic_y, equations_parabolic, dg, i, j, element) for ii in eachnode(dg) @@ -198,11 +198,11 @@ function calc_volume_integral!(du, flux_viscous, mesh::TreeMesh{2}, return nothing end -# This is the version used when calculating the divergence of the viscous fluxes. -# Specialization `flux_viscous::Tuple` needed to +# This is the version used when calculating the divergence of the parabolic fluxes. +# Specialization `flux_parabolic::Tuple` needed to # avoid amibiguity with the hyperbolic version of `prolong2interfaces!` in dg_2d.jl # which is for the variables itself, i.e., `u::Array{uEltype, 4}`. -function prolong2interfaces!(cache, flux_viscous::Tuple, +function prolong2interfaces!(cache, flux_parabolic::Tuple, mesh::TreeMesh{2}, equations_parabolic::AbstractEquationsParabolic, dg::DG) @@ -212,7 +212,7 @@ function prolong2interfaces!(cache, flux_viscous::Tuple, # OBS! `interfaces_u` stores the interpolated *fluxes* and *not the solution*! interfaces_u = interfaces.u - flux_viscous_x, flux_viscous_y = flux_viscous + flux_parabolic_x, flux_parabolic_y = flux_parabolic @threaded for interface in eachinterface(dg, cache) left_element = neighbor_ids[1, interface] @@ -221,18 +221,18 @@ function prolong2interfaces!(cache, flux_viscous::Tuple, if orientations[interface] == 1 # interface in x-direction for j in eachnode(dg), v in eachvariable(equations_parabolic) - interfaces_u[1, v, j, interface] = flux_viscous_x[v, nnodes(dg), j, - left_element] - interfaces_u[2, v, j, interface] = flux_viscous_x[v, 1, j, - right_element] + interfaces_u[1, v, j, interface] = flux_parabolic_x[v, nnodes(dg), j, + left_element] + interfaces_u[2, v, j, interface] = flux_parabolic_x[v, 1, j, + right_element] end else # if orientations[interface] == 2 # interface in y-direction for i in eachnode(dg), v in eachvariable(equations_parabolic) - interfaces_u[1, v, i, interface] = flux_viscous_y[v, i, nnodes(dg), - left_element] - interfaces_u[2, v, i, interface] = flux_viscous_y[v, i, 1, - right_element] + interfaces_u[1, v, i, interface] = flux_parabolic_y[v, i, nnodes(dg), + left_element] + interfaces_u[2, v, i, interface] = flux_parabolic_y[v, i, 1, + right_element] end end end @@ -240,11 +240,11 @@ function prolong2interfaces!(cache, flux_viscous::Tuple, return nothing end -# This is the version used when calculating the divergence of the viscous fluxes. -# Specialization `flux_viscous::Tuple` needed to +# This is the version used when calculating the divergence of the parabolic fluxes. +# Specialization `flux_parabolic::Tuple` needed to # avoid amibiguity with the hyperbolic version of `prolong2interfaces!` in dg_2d.jl # which is for the variables itself, i.e., `u::Array{uEltype, 4}`. -function prolong2interfaces!(cache, flux_viscous::Tuple, +function prolong2interfaces!(cache, flux_parabolic::Tuple, mesh::TreeMesh{2}, equations_parabolic::AbstractEquationsParabolic, dg::DGSEM{<:GaussLegendreBasis}) @@ -255,7 +255,7 @@ function prolong2interfaces!(cache, flux_viscous::Tuple, # OBS! `interfaces_u` stores the interpolated *fluxes* and *not the solution*! interfaces_u = interfaces.u - flux_viscous_x, flux_viscous_y = flux_viscous + flux_parabolic_x, flux_parabolic_y = flux_parabolic @threaded for interface in eachinterface(dg, cache) left_element = neighbor_ids[1, interface] @@ -274,11 +274,11 @@ function prolong2interfaces!(cache, flux_viscous::Tuple, # (see comment at the top of the file) # Need `boundary_interpolation` at right (+1) node for left element interface_u_1 = (interface_u_1 + - flux_viscous_x[v, ii, j, left_element] * + flux_parabolic_x[v, ii, j, left_element] * boundary_interpolation[ii, 2]) # Need `boundary_interpolation` at left (-1) node for right element interface_u_2 = (interface_u_2 + - flux_viscous_x[v, ii, j, right_element] * + flux_parabolic_x[v, ii, j, right_element] * boundary_interpolation[ii, 1]) end interfaces_u[1, v, j, interface] = interface_u_1 @@ -294,11 +294,11 @@ function prolong2interfaces!(cache, flux_viscous::Tuple, for jj in eachnode(dg) # Need `boundary_interpolation` at right (+1) node for left element interface_u_1 = (interface_u_1 + - flux_viscous_y[v, i, jj, left_element] * + flux_parabolic_y[v, i, jj, left_element] * boundary_interpolation[jj, 2]) # Need `boundary_interpolation` at left (-1) node for right element interface_u_2 = (interface_u_2 + - flux_viscous_y[v, i, jj, right_element] * + flux_parabolic_y[v, i, jj, right_element] * boundary_interpolation[jj, 1]) end interfaces_u[1, v, i, interface] = interface_u_1 @@ -311,7 +311,7 @@ function prolong2interfaces!(cache, flux_viscous::Tuple, return nothing end -# This is the version used when calculating the divergence of the viscous fluxes +# This is the version used when calculating the divergence of the parabolic fluxes function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{2}, equations_parabolic, dg::DG, parabolic_scheme, cache) @@ -348,11 +348,11 @@ function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{2}, return nothing end -# This is the version used when calculating the divergence of the viscous fluxes. -# Specialization `flux_viscous::Tuple` needed to +# This is the version used when calculating the divergence of the parabolic fluxes. +# Specialization `flux_parabolic::Tuple` needed to # avoid amibiguity with the hyperbolic version of `prolong2boundaries!` in dg_2d.jl # which is for the variables itself, i.e., `u::Array{uEltype, 4}`. -function prolong2boundaries!(cache, flux_viscous::Tuple, +function prolong2boundaries!(cache, flux_parabolic::Tuple, mesh::TreeMesh{2}, equations_parabolic::AbstractEquationsParabolic, dg::DG) @@ -361,7 +361,7 @@ function prolong2boundaries!(cache, flux_viscous::Tuple, # OBS! `boundaries_u` stores the "interpolated" *fluxes* and *not the solution*! boundaries_u = boundaries.u - flux_viscous_x, flux_viscous_y = flux_viscous + flux_parabolic_x, flux_parabolic_y = flux_parabolic @threaded for boundary in eachboundary(dg, cache) element = neighbor_ids[boundary] @@ -371,13 +371,13 @@ function prolong2boundaries!(cache, flux_viscous::Tuple, if neighbor_sides[boundary] == 1 # element in -x direction of boundary for l in eachnode(dg), v in eachvariable(equations_parabolic) - boundaries_u[1, v, l, boundary] = flux_viscous_x[v, nnodes(dg), l, - element] + boundaries_u[1, v, l, boundary] = flux_parabolic_x[v, nnodes(dg), l, + element] end else # Element in +x direction of boundary for l in eachnode(dg), v in eachvariable(equations_parabolic) - boundaries_u[2, v, l, boundary] = flux_viscous_x[v, 1, l, - element] + boundaries_u[2, v, l, boundary] = flux_parabolic_x[v, 1, l, + element] end end else # if orientations[boundary] == 2 @@ -385,14 +385,14 @@ function prolong2boundaries!(cache, flux_viscous::Tuple, if neighbor_sides[boundary] == 1 # element in -y direction of boundary for l in eachnode(dg), v in eachvariable(equations_parabolic) - boundaries_u[1, v, l, boundary] = flux_viscous_y[v, l, nnodes(dg), - element] + boundaries_u[1, v, l, boundary] = flux_parabolic_y[v, l, nnodes(dg), + element] end else # element in +y direction of boundary for l in eachnode(dg), v in eachvariable(equations_parabolic) - boundaries_u[2, v, l, boundary] = flux_viscous_y[v, l, 1, - element] + boundaries_u[2, v, l, boundary] = flux_parabolic_y[v, l, 1, + element] end end end @@ -401,11 +401,11 @@ function prolong2boundaries!(cache, flux_viscous::Tuple, return nothing end -# This is the version used when calculating the divergence of the viscous fluxes. -# Specialization `flux_viscous::Tuple` needed to +# This is the version used when calculating the divergence of the parabolic fluxes. +# Specialization `flux_parabolic::Tuple` needed to # avoid amibiguity with the hyperbolic version of `prolong2boundaries!` in dg_2d.jl # which is for the variables itself, i.e., `u::Array{uEltype, 4}`. -function prolong2boundaries!(cache, flux_viscous::Tuple, +function prolong2boundaries!(cache, flux_parabolic::Tuple, mesh::TreeMesh{2}, equations_parabolic::AbstractEquationsParabolic, dg::DGSEM{<:GaussLegendreBasis}) @@ -415,7 +415,7 @@ function prolong2boundaries!(cache, flux_viscous::Tuple, # OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*! boundaries_u = boundaries.u - flux_viscous_x, flux_viscous_y = flux_viscous + flux_parabolic_x, flux_parabolic_y = flux_parabolic @threaded for boundary in eachboundary(dg, cache) element = neighbor_ids[boundary] @@ -433,7 +433,7 @@ function prolong2boundaries!(cache, flux_viscous::Tuple, # Not += to allow `@muladd` to turn these into FMAs # (see comment at the top of the file) boundary_u = (boundary_u + - flux_viscous_x[v, ii, l, element] * + flux_parabolic_x[v, ii, l, element] * boundary_interpolation[ii, 2]) end boundaries_u[1, v, l, boundary] = boundary_u @@ -445,7 +445,7 @@ function prolong2boundaries!(cache, flux_viscous::Tuple, boundary_u = zero(eltype(boundaries_u)) for ii in eachnode(dg) boundary_u = (boundary_u + - flux_viscous_x[v, ii, l, element] * + flux_parabolic_x[v, ii, l, element] * boundary_interpolation[ii, 1]) end boundaries_u[2, v, l, boundary] = boundary_u @@ -461,7 +461,7 @@ function prolong2boundaries!(cache, flux_viscous::Tuple, boundary_u = zero(eltype(boundaries_u)) for jj in eachnode(dg) boundary_u = (boundary_u + - flux_viscous_y[v, l, jj, element] * + flux_parabolic_y[v, l, jj, element] * boundary_interpolation[jj, 2]) end boundaries_u[1, v, l, boundary] = boundary_u @@ -473,7 +473,7 @@ function prolong2boundaries!(cache, flux_viscous::Tuple, boundary_u = zero(eltype(boundaries_u)) for jj in eachnode(dg) boundary_u = (boundary_u + - flux_viscous_y[v, l, jj, element] * + flux_parabolic_y[v, l, jj, element] * boundary_interpolation[jj, 1]) end boundaries_u[2, v, l, boundary] = boundary_u @@ -486,13 +486,13 @@ function prolong2boundaries!(cache, flux_viscous::Tuple, return nothing end -function calc_viscous_fluxes!(flux_viscous, - gradients, u_transformed, - mesh::Union{TreeMesh{2}, P4estMesh{2}}, - equations_parabolic::AbstractEquationsParabolic, - dg::DG, cache) +function calc_parabolic_fluxes!(flux_parabolic, + gradients, u_transformed, + mesh::Union{TreeMesh{2}, P4estMesh{2}}, + equations_parabolic::AbstractEquationsParabolic, + dg::DG, cache) gradients_x, gradients_y = gradients - flux_viscous_x, flux_viscous_y = flux_viscous # output arrays + flux_parabolic_x, flux_parabolic_y = flux_parabolic # output arrays @threaded for element in eachelement(dg, cache) for j in eachnode(dg), i in eachnode(dg) @@ -504,14 +504,16 @@ function calc_viscous_fluxes!(flux_viscous, gradients_2_node = get_node_vars(gradients_y, equations_parabolic, dg, i, j, element) - # Calculate viscous flux and store each component for later use - flux_viscous_node_x = flux(u_node, (gradients_1_node, gradients_2_node), 1, - equations_parabolic) - flux_viscous_node_y = flux(u_node, (gradients_1_node, gradients_2_node), 2, - equations_parabolic) - set_node_vars!(flux_viscous_x, flux_viscous_node_x, equations_parabolic, dg, + # Calculate parabolic flux and store each component for later use + flux_parabolic_node_x = flux(u_node, (gradients_1_node, gradients_2_node), + 1, equations_parabolic) + flux_parabolic_node_y = flux(u_node, (gradients_1_node, gradients_2_node), + 2, equations_parabolic) + set_node_vars!(flux_parabolic_x, flux_parabolic_node_x, + equations_parabolic, dg, i, j, element) - set_node_vars!(flux_viscous_y, flux_viscous_node_y, equations_parabolic, dg, + set_node_vars!(flux_parabolic_y, flux_parabolic_node_y, + equations_parabolic, dg, i, j, element) end end @@ -674,7 +676,7 @@ function calc_boundary_flux_by_direction_divergence!(surface_flux_values::Abstra @unpack surface_flux = surface_integral # Note: cache.boundaries.u contains the unsigned normal component (using "orientation", not "direction") - # of the viscous flux, as computed in `prolong2boundaries!` + # of the parabolic flux, as computed in `prolong2boundaries!` @unpack u, neighbor_ids, neighbor_sides, node_coordinates, orientations = cache.boundaries @threaded for boundary in first_boundary:last_boundary @@ -682,7 +684,7 @@ function calc_boundary_flux_by_direction_divergence!(surface_flux_values::Abstra neighbor = neighbor_ids[boundary] for i in eachnode(dg) - # Get viscous boundary fluxes + # Get parabolic boundary fluxes flux_ll, flux_rr = get_surface_node_vars(u, equations_parabolic, dg, i, boundary) if neighbor_sides[boundary] == 1 # Element is on the left, boundary on the right @@ -713,15 +715,15 @@ function calc_boundary_flux_by_direction_divergence!(surface_flux_values::Abstra return nothing end -# Specialization `flux_viscous::Tuple` needed to +# Specialization `flux_parabolic::Tuple` needed to # avoid amibiguity with the hyperbolic version of `prolong2mortars!` in dg_2d.jl # which is for the variables itself, i.e., `u::Array{uEltype, 4}`. -function prolong2mortars!(cache, flux_viscous::Tuple, +function prolong2mortars!(cache, flux_parabolic::Tuple, mesh::TreeMesh{2}, equations_parabolic::AbstractEquationsParabolic, mortar_l2::LobattoLegendreMortarL2, dg::DGSEM) - flux_viscous_x, flux_viscous_y = flux_viscous + flux_parabolic_x, flux_parabolic_y = flux_parabolic @threaded for mortar in eachmortar(dg, cache) large_element = cache.mortars.neighbor_ids[3, mortar] upper_element = cache.mortars.neighbor_ids[2, mortar] @@ -733,28 +735,28 @@ function prolong2mortars!(cache, flux_viscous::Tuple, # L2 mortars in x-direction for l in eachnode(dg) for v in eachvariable(equations_parabolic) - cache.mortars.u_upper[2, v, l, mortar] = flux_viscous_x[v, - 1, - l, - upper_element] - cache.mortars.u_lower[2, v, l, mortar] = flux_viscous_x[v, - 1, - l, - lower_element] + cache.mortars.u_upper[2, v, l, mortar] = flux_parabolic_x[v, + 1, + l, + upper_element] + cache.mortars.u_lower[2, v, l, mortar] = flux_parabolic_x[v, + 1, + l, + lower_element] end end else # L2 mortars in y-direction for l in eachnode(dg) for v in eachvariable(equations_parabolic) - cache.mortars.u_upper[2, v, l, mortar] = flux_viscous_y[v, - l, - 1, - upper_element] - cache.mortars.u_lower[2, v, l, mortar] = flux_viscous_y[v, - l, - 1, - lower_element] + cache.mortars.u_upper[2, v, l, mortar] = flux_parabolic_y[v, + l, + 1, + upper_element] + cache.mortars.u_lower[2, v, l, mortar] = flux_parabolic_y[v, + l, + 1, + lower_element] end end end @@ -763,28 +765,28 @@ function prolong2mortars!(cache, flux_viscous::Tuple, # L2 mortars in x-direction for l in eachnode(dg) for v in eachvariable(equations_parabolic) - cache.mortars.u_upper[1, v, l, mortar] = flux_viscous_x[v, - nnodes(dg), - l, - upper_element] - cache.mortars.u_lower[1, v, l, mortar] = flux_viscous_x[v, - nnodes(dg), - l, - lower_element] + cache.mortars.u_upper[1, v, l, mortar] = flux_parabolic_x[v, + nnodes(dg), + l, + upper_element] + cache.mortars.u_lower[1, v, l, mortar] = flux_parabolic_x[v, + nnodes(dg), + l, + lower_element] end end else # L2 mortars in y-direction for l in eachnode(dg) for v in eachvariable(equations_parabolic) - cache.mortars.u_upper[1, v, l, mortar] = flux_viscous_y[v, - l, - nnodes(dg), - upper_element] - cache.mortars.u_lower[1, v, l, mortar] = flux_viscous_y[v, - l, - nnodes(dg), - lower_element] + cache.mortars.u_upper[1, v, l, mortar] = flux_parabolic_y[v, + l, + nnodes(dg), + upper_element] + cache.mortars.u_lower[1, v, l, mortar] = flux_parabolic_y[v, + l, + nnodes(dg), + lower_element] end end end @@ -795,12 +797,12 @@ function prolong2mortars!(cache, flux_viscous::Tuple, leftright = 1 if cache.mortars.orientations[mortar] == 1 # L2 mortars in x-direction - u_large = view(flux_viscous_x, :, nnodes(dg), :, large_element) + u_large = view(flux_parabolic_x, :, nnodes(dg), :, large_element) element_solutions_to_mortars!(cache.mortars, mortar_l2, leftright, mortar, u_large) else # L2 mortars in y-direction - u_large = view(flux_viscous_y, :, :, nnodes(dg), large_element) + u_large = view(flux_parabolic_y, :, :, nnodes(dg), large_element) element_solutions_to_mortars!(cache.mortars, mortar_l2, leftright, mortar, u_large) end @@ -808,12 +810,12 @@ function prolong2mortars!(cache, flux_viscous::Tuple, leftright = 2 if cache.mortars.orientations[mortar] == 1 # L2 mortars in x-direction - u_large = view(flux_viscous_x, :, 1, :, large_element) + u_large = view(flux_parabolic_x, :, 1, :, large_element) element_solutions_to_mortars!(cache.mortars, mortar_l2, leftright, mortar, u_large) else # L2 mortars in y-direction - u_large = view(flux_viscous_y, :, :, 1, large_element) + u_large = view(flux_parabolic_y, :, :, 1, large_element) element_solutions_to_mortars!(cache.mortars, mortar_l2, leftright, mortar, u_large) end @@ -1241,7 +1243,7 @@ end # Needed to *not* flip the sign of the inverse Jacobian. # This is because the parabolic fluxes are assumed to be of the form # `du/dt + df/dx = dg/dx + source(x,t)`, -# where f(u) is the inviscid flux and g(u) is the viscous flux. +# where f(u) is the inviscid flux and g(u) is the parabolic flux. function apply_jacobian_parabolic!(du::AbstractArray, mesh::TreeMesh{2}, equations_parabolic::AbstractEquationsParabolic, dg::DG, cache) diff --git a/src/solvers/dgsem_tree/dg_3d_parabolic.jl b/src/solvers/dgsem_tree/dg_3d_parabolic.jl index 79b9aaddcec..986778085ec 100644 --- a/src/solvers/dgsem_tree/dg_3d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_3d_parabolic.jl @@ -11,11 +11,11 @@ function create_cache_parabolic(mesh::Union{TreeMesh{3}, P4estMesh{3}}, equations_hyperbolic::AbstractEquations, dg::DG, n_elements, uEltype) - viscous_container = init_viscous_container_3d(nvariables(equations_hyperbolic), - nnodes(dg), n_elements, - uEltype) + parabolic_container = init_parabolic_container_3d(nvariables(equations_hyperbolic), + nnodes(dg), n_elements, + uEltype) - cache_parabolic = (; viscous_container) + cache_parabolic = (; parabolic_container) return cache_parabolic end @@ -51,23 +51,23 @@ function reset_gradients!(gradients::NTuple{3}, dg::DG, cache) return nothing end -# This is the version used when calculating the divergence of the viscous fluxes. +# This is the version used when calculating the divergence of the parabolic fluxes. # Identical to weak-form volume integral/kernel for the purely hyperbolic case, -# except that the fluxes are here already precomputed in `calc_viscous_fluxes!` -function calc_volume_integral!(du, flux_viscous, mesh::TreeMesh{3}, +# except that the fluxes are here already precomputed in `calc_parabolic_fluxes!` +function calc_volume_integral!(du, flux_parabolic, mesh::TreeMesh{3}, equations_parabolic::AbstractEquationsParabolic, dg::DGSEM, cache) @unpack derivative_hat = dg.basis - flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous + flux_parabolic_x, flux_parabolic_y, flux_parabolic_z = flux_parabolic @threaded for element in eachelement(dg, cache) # Calculate volume terms in one element for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) - flux_1_node = get_node_vars(flux_viscous_x, equations_parabolic, dg, + flux_1_node = get_node_vars(flux_parabolic_x, equations_parabolic, dg, i, j, k, element) - flux_2_node = get_node_vars(flux_viscous_y, equations_parabolic, dg, + flux_2_node = get_node_vars(flux_parabolic_y, equations_parabolic, dg, i, j, k, element) - flux_3_node = get_node_vars(flux_viscous_z, equations_parabolic, dg, + flux_3_node = get_node_vars(flux_parabolic_z, equations_parabolic, dg, i, j, k, element) for ii in eachnode(dg) @@ -90,11 +90,11 @@ function calc_volume_integral!(du, flux_viscous, mesh::TreeMesh{3}, return nothing end -# This is the version used when calculating the divergence of the viscous fluxes. -# Specialization `flux_viscous::Tuple` needed to +# This is the version used when calculating the divergence of the parabolic fluxes. +# Specialization `flux_parabolic::Tuple` needed to # avoid amibiguity with the hyperbolic version of `prolong2interfaces!` in dg_3d.jl # which is for the variables itself, i.e., `u::Array{uEltype, 5}`. -function prolong2interfaces!(cache, flux_viscous::Tuple, +function prolong2interfaces!(cache, flux_parabolic::Tuple, mesh::TreeMesh{3}, equations_parabolic::AbstractEquationsParabolic, dg::DG) @@ -104,7 +104,7 @@ function prolong2interfaces!(cache, flux_viscous::Tuple, # OBS! `interfaces_u` stores the interpolated *fluxes* and *not the solution*! interfaces_u = interfaces.u - flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous + flux_parabolic_x, flux_parabolic_y, flux_parabolic_z = flux_parabolic @threaded for interface in eachinterface(dg, cache) left_element = neighbor_ids[1, interface] @@ -115,36 +115,36 @@ function prolong2interfaces!(cache, flux_viscous::Tuple, for k in eachnode(dg), j in eachnode(dg), v in eachvariable(equations_parabolic) - interfaces_u[1, v, j, k, interface] = flux_viscous_x[v, - nnodes(dg), j, k, - left_element] - interfaces_u[2, v, j, k, interface] = flux_viscous_x[v, - 1, j, k, - right_element] + interfaces_u[1, v, j, k, interface] = flux_parabolic_x[v, + nnodes(dg), j, k, + left_element] + interfaces_u[2, v, j, k, interface] = flux_parabolic_x[v, + 1, j, k, + right_element] end elseif orientations[interface] == 2 # interface in y-direction for k in eachnode(dg), i in eachnode(dg), v in eachvariable(equations_parabolic) - interfaces_u[1, v, i, k, interface] = flux_viscous_y[v, - i, nnodes(dg), k, - left_element] - interfaces_u[2, v, i, k, interface] = flux_viscous_y[v, - i, 1, k, - right_element] + interfaces_u[1, v, i, k, interface] = flux_parabolic_y[v, + i, nnodes(dg), k, + left_element] + interfaces_u[2, v, i, k, interface] = flux_parabolic_y[v, + i, 1, k, + right_element] end else # if orientations[interface] == 3 # interface in z-direction for j in eachnode(dg), i in eachnode(dg), v in eachvariable(equations_parabolic) - interfaces_u[1, v, i, j, interface] = flux_viscous_z[v, - i, j, nnodes(dg), - left_element] - interfaces_u[2, v, i, j, interface] = flux_viscous_z[v, - i, j, 1, - right_element] + interfaces_u[1, v, i, j, interface] = flux_parabolic_z[v, + i, j, nnodes(dg), + left_element] + interfaces_u[2, v, i, j, interface] = flux_parabolic_z[v, + i, j, 1, + right_element] end end end @@ -152,7 +152,7 @@ function prolong2interfaces!(cache, flux_viscous::Tuple, return nothing end -# This is the version used when calculating the divergence of the viscous fluxes +# This is the version used when calculating the divergence of the parabolic fluxes function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{3}, equations_parabolic, dg::DG, parabolic_scheme, cache) @@ -190,11 +190,11 @@ function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{3}, return nothing end -# This is the version used when calculating the divergence of the viscous fluxes. -# Specialization `flux_viscous::Tuple` needed to +# This is the version used when calculating the divergence of the parabolic fluxes. +# Specialization `flux_parabolic::Tuple` needed to # avoid amibiguity with the hyperbolic version of `prolong2boundaries!` in dg_3d.jl # which is for the variables itself, i.e., `u::Array{uEltype, 5}`. -function prolong2boundaries!(cache, flux_viscous::Tuple, +function prolong2boundaries!(cache, flux_parabolic::Tuple, mesh::TreeMesh{3}, equations_parabolic::AbstractEquationsParabolic, dg::DG) @@ -203,7 +203,7 @@ function prolong2boundaries!(cache, flux_viscous::Tuple, # OBS! `boundaries_u` stores the "interpolated" *fluxes* and *not the solution*! boundaries_u = boundaries.u - flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous + flux_parabolic_x, flux_parabolic_y, flux_parabolic_z = flux_parabolic @threaded for boundary in eachboundary(dg, cache) element = neighbor_ids[boundary] @@ -215,21 +215,21 @@ function prolong2boundaries!(cache, flux_viscous::Tuple, for k in eachnode(dg), j in eachnode(dg), v in eachvariable(equations_parabolic) - boundaries_u[1, v, j, k, boundary] = flux_viscous_x[v, - nnodes(dg), - j, - k, - element] + boundaries_u[1, v, j, k, boundary] = flux_parabolic_x[v, + nnodes(dg), + j, + k, + element] end else # Element in +x direction of boundary for k in eachnode(dg), j in eachnode(dg), v in eachvariable(equations_parabolic) - boundaries_u[2, v, j, k, boundary] = flux_viscous_x[v, - 1, - j, - k, - element] + boundaries_u[2, v, j, k, boundary] = flux_parabolic_x[v, + 1, + j, + k, + element] end end elseif orientations[boundary] == 2 @@ -239,22 +239,22 @@ function prolong2boundaries!(cache, flux_viscous::Tuple, for k in eachnode(dg), i in eachnode(dg), v in eachvariable(equations_parabolic) - boundaries_u[1, v, i, k, boundary] = flux_viscous_y[v, - i, - nnodes(dg), - k, - element] + boundaries_u[1, v, i, k, boundary] = flux_parabolic_y[v, + i, + nnodes(dg), + k, + element] end else # element in +y direction of boundary for k in eachnode(dg), i in eachnode(dg), v in eachvariable(equations_parabolic) - boundaries_u[2, v, i, k, boundary] = flux_viscous_y[v, - i, - 1, - k, - element] + boundaries_u[2, v, i, k, boundary] = flux_parabolic_y[v, + i, + 1, + k, + element] end end else # if orientations[boundary] == 3 @@ -264,22 +264,22 @@ function prolong2boundaries!(cache, flux_viscous::Tuple, for j in eachnode(dg), i in eachnode(dg), v in eachvariable(equations_parabolic) - boundaries_u[1, v, i, j, boundary] = flux_viscous_z[v, - i, - j, - nnodes(dg), - element] + boundaries_u[1, v, i, j, boundary] = flux_parabolic_z[v, + i, + j, + nnodes(dg), + element] end else # element in +z direction of boundary for j in eachnode(dg), i in eachnode(dg), v in eachvariable(equations_parabolic) - boundaries_u[2, v, i, j, boundary] = flux_viscous_z[v, - i, - j, - 1, - element] + boundaries_u[2, v, i, j, boundary] = flux_parabolic_z[v, + i, + j, + 1, + element] end end end @@ -288,13 +288,13 @@ function prolong2boundaries!(cache, flux_viscous::Tuple, return nothing end -function calc_viscous_fluxes!(flux_viscous, - gradients, u_transformed, - mesh::Union{TreeMesh{3}, P4estMesh{3}}, - equations_parabolic::AbstractEquationsParabolic, - dg::DG, cache) +function calc_parabolic_fluxes!(flux_parabolic, + gradients, u_transformed, + mesh::Union{TreeMesh{3}, P4estMesh{3}}, + equations_parabolic::AbstractEquationsParabolic, + dg::DG, cache) gradients_x, gradients_y, gradients_z = gradients - flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous # output arrays + flux_parabolic_x, flux_parabolic_y, flux_parabolic_z = flux_parabolic # output arrays @threaded for element in eachelement(dg, cache) for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) @@ -308,21 +308,24 @@ function calc_viscous_fluxes!(flux_viscous, gradients_3_node = get_node_vars(gradients_z, equations_parabolic, dg, i, j, k, element) - # Calculate viscous flux and store each component for later use - flux_viscous_node_x = flux(u_node, - (gradients_1_node, gradients_2_node, - gradients_3_node), 1, equations_parabolic) - flux_viscous_node_y = flux(u_node, - (gradients_1_node, gradients_2_node, - gradients_3_node), 2, equations_parabolic) - flux_viscous_node_z = flux(u_node, - (gradients_1_node, gradients_2_node, - gradients_3_node), 3, equations_parabolic) - set_node_vars!(flux_viscous_x, flux_viscous_node_x, equations_parabolic, dg, + # Calculate parabolic flux and store each component for later use + flux_parabolic_node_x = flux(u_node, + (gradients_1_node, gradients_2_node, + gradients_3_node), 1, equations_parabolic) + flux_parabolic_node_y = flux(u_node, + (gradients_1_node, gradients_2_node, + gradients_3_node), 2, equations_parabolic) + flux_parabolic_node_z = flux(u_node, + (gradients_1_node, gradients_2_node, + gradients_3_node), 3, equations_parabolic) + set_node_vars!(flux_parabolic_x, flux_parabolic_node_x, + equations_parabolic, dg, i, j, k, element) - set_node_vars!(flux_viscous_y, flux_viscous_node_y, equations_parabolic, dg, + set_node_vars!(flux_parabolic_y, flux_parabolic_node_y, + equations_parabolic, dg, i, j, k, element) - set_node_vars!(flux_viscous_z, flux_viscous_node_z, equations_parabolic, dg, + set_node_vars!(flux_parabolic_z, flux_parabolic_node_z, + equations_parabolic, dg, i, j, k, element) end end @@ -510,7 +513,7 @@ function calc_boundary_flux_by_direction_divergence!(surface_flux_values::Abstra @unpack surface_flux = surface_integral # Note: cache.boundaries.u contains the unsigned normal component (using "orientation", not "direction") - # of the viscous flux, as computed in `prolong2boundaries!` + # of the parabolic flux, as computed in `prolong2boundaries!` @unpack u, neighbor_ids, neighbor_sides, node_coordinates, orientations = cache.boundaries @threaded for boundary in first_boundary:last_boundary @@ -518,7 +521,7 @@ function calc_boundary_flux_by_direction_divergence!(surface_flux_values::Abstra neighbor = neighbor_ids[boundary] for j in eachnode(dg), i in eachnode(dg) - # Get viscous boundary fluxes + # Get parabolic boundary fluxes flux_ll, flux_rr = get_surface_node_vars(u, equations_parabolic, dg, i, j, boundary) if neighbor_sides[boundary] == 1 # Element is on the left, boundary on the right @@ -550,17 +553,17 @@ function calc_boundary_flux_by_direction_divergence!(surface_flux_values::Abstra return nothing end -# Specialization `flux_viscous::Tuple` needed to +# Specialization `flux_parabolic::Tuple` needed to # avoid amibiguity with the hyperbolic version of `prolong2mortars!` in dg_3d.jl # which is for the variables itself, i.e., `u::Array{uEltype, 5}`. -function prolong2mortars!(cache, flux_viscous::Tuple, +function prolong2mortars!(cache, flux_parabolic::Tuple, mesh::TreeMesh{3}, equations_parabolic::AbstractEquationsParabolic, mortar_l2::LobattoLegendreMortarL2, dg::DGSEM) # temporary buffer for projections @unpack fstar_tmp1_threaded = cache - flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous + flux_parabolic_x, flux_parabolic_y, flux_parabolic_z = flux_parabolic @threaded for mortar in eachmortar(dg, cache) fstar_tmp1 = fstar_tmp1_threaded[Threads.threadid()] @@ -576,78 +579,78 @@ function prolong2mortars!(cache, flux_viscous::Tuple, # L2 mortars in x-direction for k in eachnode(dg), j in eachnode(dg) for v in eachvariable(equations_parabolic) - cache.mortars.u_upper_left[2, v, j, k, mortar] = flux_viscous_x[v, - 1, - j, - k, - upper_left_element] - cache.mortars.u_upper_right[2, v, j, k, mortar] = flux_viscous_x[v, - 1, - j, - k, - upper_right_element] - cache.mortars.u_lower_left[2, v, j, k, mortar] = flux_viscous_x[v, - 1, - j, - k, - lower_left_element] - cache.mortars.u_lower_right[2, v, j, k, mortar] = flux_viscous_x[v, - 1, - j, - k, - lower_right_element] + cache.mortars.u_upper_left[2, v, j, k, mortar] = flux_parabolic_x[v, + 1, + j, + k, + upper_left_element] + cache.mortars.u_upper_right[2, v, j, k, mortar] = flux_parabolic_x[v, + 1, + j, + k, + upper_right_element] + cache.mortars.u_lower_left[2, v, j, k, mortar] = flux_parabolic_x[v, + 1, + j, + k, + lower_left_element] + cache.mortars.u_lower_right[2, v, j, k, mortar] = flux_parabolic_x[v, + 1, + j, + k, + lower_right_element] end end elseif cache.mortars.orientations[mortar] == 2 # L2 mortars in y-direction for k in eachnode(dg), i in eachnode(dg) for v in eachvariable(equations_parabolic) - cache.mortars.u_upper_left[2, v, i, k, mortar] = flux_viscous_y[v, - i, - 1, - k, - upper_left_element] - cache.mortars.u_upper_right[2, v, i, k, mortar] = flux_viscous_y[v, - i, - 1, - k, - upper_right_element] - cache.mortars.u_lower_left[2, v, i, k, mortar] = flux_viscous_y[v, - i, - 1, - k, - lower_left_element] - cache.mortars.u_lower_right[2, v, i, k, mortar] = flux_viscous_y[v, - i, - 1, - k, - lower_right_element] + cache.mortars.u_upper_left[2, v, i, k, mortar] = flux_parabolic_y[v, + i, + 1, + k, + upper_left_element] + cache.mortars.u_upper_right[2, v, i, k, mortar] = flux_parabolic_y[v, + i, + 1, + k, + upper_right_element] + cache.mortars.u_lower_left[2, v, i, k, mortar] = flux_parabolic_y[v, + i, + 1, + k, + lower_left_element] + cache.mortars.u_lower_right[2, v, i, k, mortar] = flux_parabolic_y[v, + i, + 1, + k, + lower_right_element] end end else # orientations[mortar] == 3 # L2 mortars in z-direction for j in eachnode(dg), i in eachnode(dg) for v in eachvariable(equations_parabolic) - cache.mortars.u_upper_left[2, v, i, j, mortar] = flux_viscous_z[v, - i, - j, - 1, - upper_left_element] - cache.mortars.u_upper_right[2, v, i, j, mortar] = flux_viscous_z[v, - i, - j, - 1, - upper_right_element] - cache.mortars.u_lower_left[2, v, i, j, mortar] = flux_viscous_z[v, - i, - j, - 1, - lower_left_element] - cache.mortars.u_lower_right[2, v, i, j, mortar] = flux_viscous_z[v, - i, - j, - 1, - lower_right_element] + cache.mortars.u_upper_left[2, v, i, j, mortar] = flux_parabolic_z[v, + i, + j, + 1, + upper_left_element] + cache.mortars.u_upper_right[2, v, i, j, mortar] = flux_parabolic_z[v, + i, + j, + 1, + upper_right_element] + cache.mortars.u_lower_left[2, v, i, j, mortar] = flux_parabolic_z[v, + i, + j, + 1, + lower_left_element] + cache.mortars.u_lower_right[2, v, i, j, mortar] = flux_parabolic_z[v, + i, + j, + 1, + lower_right_element] end end end @@ -656,78 +659,78 @@ function prolong2mortars!(cache, flux_viscous::Tuple, # L2 mortars in x-direction for k in eachnode(dg), j in eachnode(dg) for v in eachvariable(equations_parabolic) - cache.mortars.u_upper_left[1, v, j, k, mortar] = flux_viscous_x[v, - nnodes(dg), - j, - k, - upper_left_element] - cache.mortars.u_upper_right[1, v, j, k, mortar] = flux_viscous_x[v, - nnodes(dg), - j, - k, - upper_right_element] - cache.mortars.u_lower_left[1, v, j, k, mortar] = flux_viscous_x[v, - nnodes(dg), - j, - k, - lower_left_element] - cache.mortars.u_lower_right[1, v, j, k, mortar] = flux_viscous_x[v, - nnodes(dg), - j, - k, - lower_right_element] + cache.mortars.u_upper_left[1, v, j, k, mortar] = flux_parabolic_x[v, + nnodes(dg), + j, + k, + upper_left_element] + cache.mortars.u_upper_right[1, v, j, k, mortar] = flux_parabolic_x[v, + nnodes(dg), + j, + k, + upper_right_element] + cache.mortars.u_lower_left[1, v, j, k, mortar] = flux_parabolic_x[v, + nnodes(dg), + j, + k, + lower_left_element] + cache.mortars.u_lower_right[1, v, j, k, mortar] = flux_parabolic_x[v, + nnodes(dg), + j, + k, + lower_right_element] end end elseif cache.mortars.orientations[mortar] == 2 # L2 mortars in y-direction for k in eachnode(dg), i in eachnode(dg) for v in eachvariable(equations_parabolic) - cache.mortars.u_upper_left[1, v, i, k, mortar] = flux_viscous_y[v, - i, - nnodes(dg), - k, - upper_left_element] - cache.mortars.u_upper_right[1, v, i, k, mortar] = flux_viscous_y[v, - i, - nnodes(dg), - k, - upper_right_element] - cache.mortars.u_lower_left[1, v, i, k, mortar] = flux_viscous_y[v, - i, - nnodes(dg), - k, - lower_left_element] - cache.mortars.u_lower_right[1, v, i, k, mortar] = flux_viscous_y[v, - i, - nnodes(dg), - k, - lower_right_element] + cache.mortars.u_upper_left[1, v, i, k, mortar] = flux_parabolic_y[v, + i, + nnodes(dg), + k, + upper_left_element] + cache.mortars.u_upper_right[1, v, i, k, mortar] = flux_parabolic_y[v, + i, + nnodes(dg), + k, + upper_right_element] + cache.mortars.u_lower_left[1, v, i, k, mortar] = flux_parabolic_y[v, + i, + nnodes(dg), + k, + lower_left_element] + cache.mortars.u_lower_right[1, v, i, k, mortar] = flux_parabolic_y[v, + i, + nnodes(dg), + k, + lower_right_element] end end else # if cache.mortars.orientations[mortar] == 3 # L2 mortars in z-direction for j in eachnode(dg), i in eachnode(dg) for v in eachvariable(equations_parabolic) - cache.mortars.u_upper_left[1, v, i, j, mortar] = flux_viscous_z[v, - i, - j, - nnodes(dg), - upper_left_element] - cache.mortars.u_upper_right[1, v, i, j, mortar] = flux_viscous_z[v, - i, - j, - nnodes(dg), - upper_right_element] - cache.mortars.u_lower_left[1, v, i, j, mortar] = flux_viscous_z[v, - i, - j, - nnodes(dg), - lower_left_element] - cache.mortars.u_lower_right[1, v, i, j, mortar] = flux_viscous_z[v, - i, - j, - nnodes(dg), - lower_right_element] + cache.mortars.u_upper_left[1, v, i, j, mortar] = flux_parabolic_z[v, + i, + j, + nnodes(dg), + upper_left_element] + cache.mortars.u_upper_right[1, v, i, j, mortar] = flux_parabolic_z[v, + i, + j, + nnodes(dg), + upper_right_element] + cache.mortars.u_lower_left[1, v, i, j, mortar] = flux_parabolic_z[v, + i, + j, + nnodes(dg), + lower_left_element] + cache.mortars.u_lower_right[1, v, i, j, mortar] = flux_parabolic_z[v, + i, + j, + nnodes(dg), + lower_right_element] end end end @@ -738,17 +741,17 @@ function prolong2mortars!(cache, flux_viscous::Tuple, leftright = 1 if cache.mortars.orientations[mortar] == 1 # L2 mortars in x-direction - u_large = view(flux_viscous_x, :, nnodes(dg), :, :, large_element) + u_large = view(flux_parabolic_x, :, nnodes(dg), :, :, large_element) element_solutions_to_mortars!(cache.mortars, mortar_l2, leftright, mortar, u_large, fstar_tmp1) elseif cache.mortars.orientations[mortar] == 2 # L2 mortars in y-direction - u_large = view(flux_viscous_y, :, :, nnodes(dg), :, large_element) + u_large = view(flux_parabolic_y, :, :, nnodes(dg), :, large_element) element_solutions_to_mortars!(cache.mortars, mortar_l2, leftright, mortar, u_large, fstar_tmp1) else # cache.mortars.orientations[mortar] == 3 # L2 mortars in z-direction - u_large = view(flux_viscous_z, :, :, :, nnodes(dg), large_element) + u_large = view(flux_parabolic_z, :, :, :, nnodes(dg), large_element) element_solutions_to_mortars!(cache.mortars, mortar_l2, leftright, mortar, u_large, fstar_tmp1) end @@ -756,17 +759,17 @@ function prolong2mortars!(cache, flux_viscous::Tuple, leftright = 2 if cache.mortars.orientations[mortar] == 1 # L2 mortars in x-direction - u_large = view(flux_viscous_x, :, 1, :, :, large_element) + u_large = view(flux_parabolic_x, :, 1, :, :, large_element) element_solutions_to_mortars!(cache.mortars, mortar_l2, leftright, mortar, u_large, fstar_tmp1) elseif cache.mortars.orientations[mortar] == 2 # L2 mortars in y-direction - u_large = view(flux_viscous_y, :, :, 1, :, large_element) + u_large = view(flux_parabolic_y, :, :, 1, :, large_element) element_solutions_to_mortars!(cache.mortars, mortar_l2, leftright, mortar, u_large, fstar_tmp1) else # cache.mortars.orientations[mortar] == 3 # L2 mortars in z-direction - u_large = view(flux_viscous_z, :, :, :, 1, large_element) + u_large = view(flux_parabolic_z, :, :, :, 1, large_element) element_solutions_to_mortars!(cache.mortars, mortar_l2, leftright, mortar, u_large, fstar_tmp1) end @@ -1128,7 +1131,7 @@ end # Needed to *not* flip the sign of the inverse Jacobian. # This is because the parabolic fluxes are assumed to be of the form # `du/dt + df/dx = dg/dx + source(x,t)`, -# where f(u) is the inviscid flux and g(u) is the viscous flux. +# where f(u) is the inviscid flux and g(u) is the parabolic flux. function apply_jacobian_parabolic!(du::AbstractArray, mesh::TreeMesh{3}, equations_parabolic::AbstractEquationsParabolic, dg::DG, cache) diff --git a/src/solvers/solvers_parabolic.jl b/src/solvers/solvers_parabolic.jl index 9f502547d36..d183c1f04f3 100644 --- a/src/solvers/solvers_parabolic.jl +++ b/src/solvers/solvers_parabolic.jl @@ -1,5 +1,5 @@ """ - ViscousFormulationBassiRebay1() + ParabolicFormulationBassiRebay1() The classical BR1 flux from @@ -15,59 +15,59 @@ A more detailed study of the BR1 scheme for the DGSEM can be found in The BR1 scheme works well for convection-dominated problems, but may cause instabilities or reduced convergence for diffusion-dominated problems. -In the latter case, the [`ViscousFormulationLocalDG`](@ref) scheme is recommended. +In the latter case, the [`ParabolicFormulationLocalDG`](@ref) scheme is recommended. """ -struct ViscousFormulationBassiRebay1 end +struct ParabolicFormulationBassiRebay1 end """ flux_parabolic(u_ll, u_rr, gradient_or_divergence, equations_parabolic, - parabolic_scheme::ViscousFormulationBassiRebay1) + parabolic_scheme::ParabolicFormulationBassiRebay1) flux_parabolic(u_ll, u_rr, normal_direction::AbstractVector, gradient_or_divergence, equations_parabolic, - parabolic_scheme::ViscousFormulationBassiRebay1) + parabolic_scheme::ParabolicFormulationBassiRebay1) This computes the classical BR1 flux. Since the interface flux for both the DG gradient and DG divergence under BR1 are identical, this function does not need to be specialized for `Gradient` and `Divergence`. `normal_direction` is not used in the BR1 flux, -but is included as an argument for consistency with the [`ViscousFormulationLocalDG`](@ref) flux, +but is included as an argument for consistency with the [`ParabolicFormulationLocalDG`](@ref) flux, which does use the `normal_direction` to compute the LDG "switch" on the generally non-Cartesian [`P4estMesh`](@ref). """ function flux_parabolic(u_ll, u_rr, # Version for `TreeMesh` gradient_or_divergence, equations_parabolic, - parabolic_scheme::ViscousFormulationBassiRebay1) + parabolic_scheme::ParabolicFormulationBassiRebay1) return 0.5f0 * (u_ll + u_rr) end # Version for `P4estMesh` function flux_parabolic(u_ll, u_rr, normal_direction::AbstractVector, gradient_or_divergence, equations_parabolic, - parabolic_scheme::ViscousFormulationBassiRebay1) + parabolic_scheme::ParabolicFormulationBassiRebay1) return 0.5f0 * (u_ll + u_rr) end """ - ViscousFormulationLocalDG(penalty_parameter) + ParabolicFormulationLocalDG(penalty_parameter) The local DG (LDG) flux from "The Local Discontinuous Galerkin Method for Time-Dependent Convection-Diffusion Systems" by Cockburn and Shu (1998). The parabolic "upwinding" vector is currently implemented for `TreeMesh`; for all other mesh types, -the LDG solver is equivalent to [`ViscousFormulationBassiRebay1`](@ref) with an LDG-type penalization. +the LDG solver is equivalent to [`ParabolicFormulationBassiRebay1`](@ref) with an LDG-type penalization. - Cockburn and Shu (1998). The Local Discontinuous Galerkin Method for Time-Dependent Convection-Diffusion Systems [DOI: 10.1137/S0036142997316712](https://doi.org/10.1137/S0036142997316712) """ -struct ViscousFormulationLocalDG{P} +struct ParabolicFormulationLocalDG{P} penalty_parameter::P end """ - ViscousFormulationLocalDG() + ParabolicFormulationLocalDG() The minimum dissipation local DG (LDG) flux from "An Analysis of the Minimal Dissipation Local Discontinuous Galerkin Method for Convection–Diffusion Problems" by Cockburn and Dong (2007). @@ -79,16 +79,16 @@ Cockburn and Dong proved that this scheme is still stable despite the zero penal Galerkin Method for Convection–Diffusion Problems. [DOI: 10.1007/s10915-007-9130-3](https://doi.org/10.1007/s10915-007-9130-3) """ -ViscousFormulationLocalDG() = ViscousFormulationLocalDG(nothing) +ParabolicFormulationLocalDG() = ParabolicFormulationLocalDG(nothing) @doc raw""" flux_parabolic(u_ll, u_rr, ::Gradient, equations_parabolic, - parabolic_scheme::ViscousFormulationLocalDG) + parabolic_scheme::ParabolicFormulationLocalDG) flux_parabolic(u_ll, u_rr, normal_direction, ::Gradient, equations_parabolic, - parabolic_scheme::ViscousFormulationLocalDG) + parabolic_scheme::ParabolicFormulationLocalDG) These fluxes computes the gradient and divergence interface fluxes for the local DG method. The local DG method uses an "upwind/downwind" flux for the @@ -110,7 +110,7 @@ i = \text{argmax} \{ \begin{pmatrix} \vert n_1 \vert \\ \vert n_2 \vert \\ \dots """ function flux_parabolic(u_ll, u_rr, # Version for `TreeMesh` ::Gradient, equations_parabolic, - parabolic_scheme::ViscousFormulationLocalDG) + parabolic_scheme::ParabolicFormulationLocalDG) # The LDG flux is {{f}} + beta * [[f]], where beta is the LDG "switch", # which we set to -1 on the left and +1 on the right in 1D. The sign of the # jump term should be opposite that of the sign used in the divergence flux. @@ -121,7 +121,7 @@ end # Version for `P4estMesh` function flux_parabolic(u_ll, u_rr, normal_direction, ::Gradient, equations_parabolic, - parabolic_scheme::ViscousFormulationLocalDG) + parabolic_scheme::ParabolicFormulationLocalDG) # Use "Upwind in dominant direction" for LDG switch abs_max_dir = argmax(abs.(normal_direction)) ldg_switch = sign(normal_direction[abs_max_dir]) @@ -131,11 +131,11 @@ end @doc raw""" flux_parabolic(u_ll, u_rr, ::Divergence, equations_parabolic, - parabolic_scheme::ViscousFormulationLocalDG) + parabolic_scheme::ParabolicFormulationLocalDG) flux_parabolic(u_ll, u_rr, normal_direction, ::Divergence, equations_parabolic, - parabolic_scheme::ViscousFormulationLocalDG) + parabolic_scheme::ParabolicFormulationLocalDG) These fluxes computes the gradient and divergence interface fluxes for the local DG method. The local DG method uses an "upwind/downwind" flux for the @@ -158,17 +158,17 @@ i = \text{argmax} \{ \begin{pmatrix} \vert n_1 \vert \\ \vert n_2 \vert \\ \dots """ function flux_parabolic(u_ll, u_rr, # Version for `TreeMesh` ::Divergence, equations_parabolic, - parabolic_scheme::ViscousFormulationLocalDG) + parabolic_scheme::ParabolicFormulationLocalDG) return u_rr # Use the downwind value for the divergence interface flux end # Version for `P4estMesh` function flux_parabolic(u_ll, u_rr, normal_direction, ::Divergence, equations_parabolic, - parabolic_scheme::ViscousFormulationLocalDG) + parabolic_scheme::ParabolicFormulationLocalDG) # Use "Downwind in dominant direction" for LDG switch abs_max_dir = argmax(abs.(normal_direction)) ldg_switch = -sign(normal_direction[abs_max_dir]) return 0.5f0 * (u_ll + u_rr - ldg_switch * (u_rr - u_ll)) end -default_parabolic_solver() = ViscousFormulationBassiRebay1() +default_parabolic_solver() = ParabolicFormulationBassiRebay1() diff --git a/test/test_parabolic_1d.jl b/test/test_parabolic_1d.jl index 839d7e4ccdd..9fef07fed7a 100644 --- a/test/test_parabolic_1d.jl +++ b/test/test_parabolic_1d.jl @@ -487,7 +487,7 @@ end callbacks=CallbackSet(summary_callback, alive_callback, analysis_callback, StepsizeCallback(cfl = 0.5, - cfl_diffusive = 0.1)), + cfl_parabolic = 0.1)), adaptive=false, l2=[ 3.804624387087144e-5, diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index cafb1812161..a0ae0d0864a 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -68,9 +68,9 @@ isdir(outdir) && rm(outdir, recursive = true) @test getindex.(gradients[2], 1) ≈ xq .^ 2 u_flux = similar.(gradients) - Trixi.calc_viscous_fluxes!(u_flux, u0, gradients, mesh, - equations_parabolic, - dg, cache, cache_parabolic) + Trixi.calc_parabolic_fluxes!(u_flux, u0, gradients, mesh, + equations_parabolic, + dg, cache, cache_parabolic) @test u_flux[1] ≈ gradients[1] @test u_flux[2] ≈ gradients[2] @@ -223,7 +223,7 @@ end @trixi_testset "TreeMesh2D: elixir_advection_diffusion.jl (LDG)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_2d_dgsem", "elixir_advection_diffusion.jl"), - solver_parabolic=ViscousFormulationLocalDG(), + solver_parabolic=ParabolicFormulationLocalDG(), initial_refinement_level=2, tspan=(0.0, 0.4), polydeg=5, l2=[6.193056910594806e-6], linf=[4.918855889635143e-5]) # Ensure that we do not have excessive memory allocations @@ -247,7 +247,7 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_2d_dgsem", "elixir_advection_diffusion_gradient_source_terms.jl"), initial_refinement_level=2, tspan=(0.0, 0.4), - solver_parabolic=ViscousFormulationBassiRebay1(), nu=1e-3, + solver_parabolic=ParabolicFormulationBassiRebay1(), nu=1e-3, stepsize_callback=TrivialCallback(), dt=1e-1, l2=[0.0017395186758592685], linf=[0.007481527467476025]) # Ensure that we do not have excessive memory allocations @@ -264,7 +264,7 @@ end coordinates_max = coordinates_max, periodicity = true), tspan=(0.0, 0.4), - solver_parabolic=ViscousFormulationBassiRebay1(), nu=1e-3, + solver_parabolic=ParabolicFormulationBassiRebay1(), nu=1e-3, stepsize_callback=TrivialCallback(), dt=1e-1, l2=[0.0017395186758592685], linf=[0.007481527467476025]) # Ensure that we do not have excessive memory allocations @@ -358,7 +358,7 @@ end @trixi_testset "TreeMesh2D: elixir_advection_diffusion_nonperiodic_amr.jl (LDG)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_2d_dgsem", "elixir_advection_diffusion_nonperiodic_amr.jl"), - solver_parabolic=ViscousFormulationLocalDG(), + solver_parabolic=ParabolicFormulationLocalDG(), tspan=(0.0, 0.01), l2=[0.000684755734524055], linf=[0.01141444199847298]) @@ -385,7 +385,7 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_2d_dgsem", "elixir_advection_diffusion_nonperiodic.jl"), initial_refinement_level=2, tspan=(0.0, 0.1), - solver_parabolic=ViscousFormulationLocalDG(), + solver_parabolic=ParabolicFormulationLocalDG(), l2=[0.007009146246373517], linf=[0.09535203925012649]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -637,7 +637,7 @@ end @trixi_testset "TreeMesh2D: elixir_navierstokes_shearlayer_nonconforming.jl (LDG)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_2d_dgsem", "elixir_navierstokes_shearlayer_nonconforming.jl"), - solver_parabolic=ViscousFormulationLocalDG(), + solver_parabolic=ParabolicFormulationLocalDG(), l2=[ 0.005352370793583371, 0.5969444914287823, @@ -704,8 +704,8 @@ end "elixir_navierstokes_viscous_shock.jl"), solver=DGSEM(polydeg = 3, surface_flux = flux_hlle, basis_type = GaussLegendreBasis), - solver_parabolic=ViscousFormulationLocalDG(), - cfl_diffusive=0.04, + solver_parabolic=ParabolicFormulationLocalDG(), + cfl_parabolic=0.04, l2=[ 6.599006355897759e-6, 4.514805201434994e-6, @@ -796,14 +796,14 @@ end @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1000) end -@trixi_testset "P4estMesh2D: elixir_advection_diffusion_nonperiodic_amr.jl (Diffusive CFL)" begin +@trixi_testset "P4estMesh2D: elixir_advection_diffusion_nonperiodic_amr.jl (Parabolic CFL)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_2d_dgsem", "elixir_advection_diffusion_nonperiodic_amr.jl"), initial_refinement_level=2, callbacks=CallbackSet(summary_callback, analysis_callback, alive_callback, StepsizeCallback(cfl = 1.6, - cfl_diffusive = 0.2)), + cfl_parabolic = 0.2)), ode_alg=CarpenterKennedy2N54(williamson_condition = false), dt=1.0, # will be overwritten l2=[0.00010850375815619432], @@ -817,7 +817,7 @@ end @trixi_testset "P4estMesh2D: elixir_advection_diffusion_nonperiodic_amr.jl (LDG)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_2d_dgsem", "elixir_advection_diffusion_nonperiodic_amr.jl"), - solver_parabolic=ViscousFormulationLocalDG(), + solver_parabolic=ParabolicFormulationLocalDG(), tspan=(0.0, 0.01), l2=[0.0006847533999311489], linf=[0.01141430509080712]) @@ -948,7 +948,7 @@ end @trixi_testset "P4estMesh2D: elixir_navierstokes_shearlayer_nonconforming.jl (LDG)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_2d_dgsem", "elixir_navierstokes_shearlayer_nonconforming.jl"), - solver_parabolic=ViscousFormulationLocalDG(), + solver_parabolic=ParabolicFormulationLocalDG(), l2=[ 0.0053523707935916025, 0.5969444914278867, @@ -1155,7 +1155,7 @@ end callbacks=CallbackSet(summary_callback, analysis_callback, alive_callback, StepsizeCallback(cfl = 2.3, - cfl_diffusive = 1.0)), + cfl_parabolic = 1.0)), adaptive=false, # respect CFL ode_alg=CKLLSRK95_4S(), l2=[ diff --git a/test/test_parabolic_3d.jl b/test/test_parabolic_3d.jl index 4e86de3cfa7..8b192e015ce 100644 --- a/test/test_parabolic_3d.jl +++ b/test/test_parabolic_3d.jl @@ -408,7 +408,7 @@ end @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1000) end -@trixi_testset "P4estMesh3D: elixir_navierstokes_taylor_green_vortex.jl (Diffusive CFL)" begin +@trixi_testset "P4estMesh3D: elixir_navierstokes_taylor_green_vortex.jl (Parabolic CFL)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_3d_dgsem", "elixir_navierstokes_taylor_green_vortex.jl"), tspan=(0.0, 0.1), @@ -416,7 +416,7 @@ end callbacks=CallbackSet(summary_callback, analysis_callback, alive_callback, StepsizeCallback(cfl = 2.3, - cfl_diffusive = 0.4)), + cfl_parabolic = 0.4)), adaptive=false, # respect CFL ode_alg=CKLLSRK95_4S(), l2=[ @@ -457,7 +457,7 @@ end @trixi_testset "TreeMesh3D: elixir_advection_diffusion_amr.jl (LDG)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_3d_dgsem", "elixir_advection_diffusion_amr.jl"), - solver_parabolic=ViscousFormulationLocalDG(), + solver_parabolic=ParabolicFormulationLocalDG(), initial_refinement_level=2, base_level=2, med_level=3, @@ -495,7 +495,7 @@ end @trixi_testset "TreeMesh3D: elixir_advection_diffusion_nonperiodic.jl (LDG)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_3d_dgsem", "elixir_advection_diffusion_nonperiodic.jl"), - solver_parabolic=ViscousFormulationLocalDG(), + solver_parabolic=ParabolicFormulationLocalDG(), l2=[0.0009432415534931421], linf=[0.016955330290404563]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -525,8 +525,8 @@ end @trixi_testset "P4estMesh3D: elixir_advection_diffusion_nonperiodic.jl (LDG)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_3d_dgsem", "elixir_advection_diffusion_nonperiodic.jl"), - solver_parabolic=ViscousFormulationLocalDG(), - cfl_diffusive=0.07, + solver_parabolic=ParabolicFormulationLocalDG(), + cfl_parabolic=0.07, l2=[0.004185076476662267], linf=[0.05166349548111486]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -547,7 +547,7 @@ end @trixi_testset "P4estMesh3D: elixir_advection_diffusion_amr_curved.jl (LDG)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_3d_dgsem", "elixir_advection_diffusion_amr_curved.jl"), - solver_parabolic=ViscousFormulationLocalDG(), + solver_parabolic=ParabolicFormulationLocalDG(), l2=[0.0006853004145232737], linf=[0.02352694543085776]) # 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 2427593fe24..06d9c4b1192 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -1878,7 +1878,7 @@ isdir(outdir) && rm(outdir, recursive = true) RealT end - parabolic_solver = ViscousFormulationLocalDG(RealT(0.1)) + parabolic_solver = ParabolicFormulationLocalDG(RealT(0.1)) @test eltype(@inferred Trixi.penalty(u_outer, u_inner, inv_h, equations_parabolic, parabolic_solver)) == RealT @@ -1901,7 +1901,7 @@ isdir(outdir) && rm(outdir, recursive = true) RealT end - parabolic_solver = ViscousFormulationLocalDG(RealT(0.1)) + parabolic_solver = ParabolicFormulationLocalDG(RealT(0.1)) @test eltype(@inferred Trixi.penalty(u_outer, u_inner, inv_h, equations_parabolic, parabolic_solver)) == RealT