From ba4702a5ad38fba502d10cb8f018e364a2947ca6 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 25 Aug 2022 13:25:59 -0500 Subject: [PATCH 1/9] fixing a BC in parabolic tutorial --- docs/literate/src/files/parabolic_terms.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/literate/src/files/parabolic_terms.jl b/docs/literate/src/files/parabolic_terms.jl index ff972fc75fe..aeceb7b7e6f 100644 --- a/docs/literate/src/files/parabolic_terms.jl +++ b/docs/literate/src/files/parabolic_terms.jl @@ -1,4 +1,4 @@ -#src # Adding parabolic terms (advection-diffusion). +#src # Parabolic terms (advection-diffusion). # Experimental support for parabolic diffusion terms is available in Trixi.jl. # This demo illustrates parabolic terms for the advection-diffusion equation. @@ -36,7 +36,7 @@ boundary_condition_zero_dirichlet = BoundaryConditionDirichlet((x, t, equations) boundary_conditions_hyperbolic = (; x_neg = BoundaryConditionDirichlet((x, t, equations) -> SVector(1 + 0.5 * x[2])), y_neg = boundary_condition_zero_dirichlet, - y_pos = boundary_condition_zero_dirichlet, + y_pos = boundary_condition_do_nothing, x_pos = boundary_condition_do_nothing) boundary_conditions_parabolic = (; x_neg = BoundaryConditionDirichlet((x, t, equations) -> SVector(1 + 0.5 * x[2])), From 362de9c420391d279b8f90b77285ea2433fb8e4c Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 25 Aug 2022 13:26:46 -0500 Subject: [PATCH 2/9] "adding new parabolic terms" tutorial --- .../src/files/adding_new_parabolic_terms.jl | 160 ++++++++++++++++++ docs/make.jl | 1 + 2 files changed, 161 insertions(+) create mode 100644 docs/literate/src/files/adding_new_parabolic_terms.jl diff --git a/docs/literate/src/files/adding_new_parabolic_terms.jl b/docs/literate/src/files/adding_new_parabolic_terms.jl new file mode 100644 index 00000000000..761b63cb89a --- /dev/null +++ b/docs/literate/src/files/adding_new_parabolic_terms.jl @@ -0,0 +1,160 @@ +#src # Adding new parabolic terms. + +# This demo illustrates the steps involved in adding new parabolic terms for the scalar +# advection equation. In particular, we will add an anisotropic diffusion. We begin by +# defining the hyperbolic (advection) part of the advection-diffusion equation. + +using OrdinaryDiffEq +using Trixi + + +advection_velocity = (1.0, 1.0) +equations_hyperbolic = LinearScalarAdvectionEquation2D(advection_velocity); + +# Next, we define a 2D parabolic diffusion term type. This is similar to [`LaplaceDiffusion2D`](@ref) +# except that the `diffusivity` field refers to a spatially constant diffusivity matrix now. Note that +# `ConstantAnisotropicDiffusion2D` has a field for `equations_hyperbolic`. It is useful to have +# information about the hyperbolic system available to the parabolic part so that we can reuse +# functions defined for hyperbolic equations (such as `varnames`). + +struct ConstantAnisotropicDiffusion2D{E, T} <: Trixi.AbstractEquationsParabolic{2, 1} + diffusivity::T + equations_hyperbolic::E +end + +varnames(variable_mapping, equations_parabolic::ConstantAnisotropicDiffusion2D) = + varnames(variable_mapping, equations_parabolic.equations_hyperbolic) + +# Next, we define the viscous flux function. We assume that the mixed hyperbolic-parabolic system +# is of the form +# ```math +# \partial_t u(t,x) + \partial_x (f_x(u) - g_x(u, \nabla u)) +# + \partial_y (f_y(u) - g_y(u, \nabla u)) = 0 +# ``` +# where `f_x(u)`, `f_y(u)` are the hyperbolic fluxes and `g_x(u, \nabla u)`, `g_y(u, \nabla u)` denote +# the viscous fluxes. For anisotropic diffusion, the viscous 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`. + +function Trixi.flux(u, gradients, orientation::Integer, equations_parabolic::ConstantAnisotropicDiffusion2D) + @unpack diffusivity = equations_parabolic + dudx, dudy = gradients + if orientation == 1 + return SVector(diffusivity[1, 1] * dudx + diffusivity[1, 2] * dudy) + else # if orientation == 2 + return SVector(diffusivity[2, 1] * dudx + diffusivity[2, 2] * dudy) + end +end + +# ## Defining boundary conditions + +# Trixi.jl's implementation of parabolic terms discretizes both the gradient and divergence +# using weak formulation. In other words, we discretize the system +# ```math +# \begin{aligned} +# \bm{q} &= \nabla u \\ +# \bm{sigma} &= \begin{pmatrix} g_x(u, \bm{q}) \\ g_y(u, \bm{q}) \end{pmatrix} \\ +# \text{viscous 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 +# ``\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. +# +# As an example, let us introduce a Dirichlet boundary condition with constant boundary data. + +struct BoundaryConditionConstantDirichlet{T <: Real} + boundary_value::T +end + +# This boundary condition contains only the field `boundary_value`, which we assume to be some +# real-valued constant which we will impose as the Dirichlet data on the boundary. +# +# Boundary conditions have generally been defined as "callable structs" (also known as "functors). +# For each boundary condition, we need to specify the appropriate boundary data to return for both +# the `Gradient` and `Divergence`. Since the gradient is operating on the solution `u`, the boundary +# data should be the value of `u`, and we can directly impose Dirichlet data. + +@inline function (boundary_condition::BoundaryConditionConstantDirichlet)(flux_inner, u_inner, normal::AbstractVector, + x, t, operator_type::Trixi.Gradient, + equations_parabolic::ConstantAnisotropicDiffusion2D) + return boundary_condition.boundary_value +end + +# While the gradient acts on the solution `u`, the divergence acts on the viscous 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. +# +# Thus, for the `Divergence` boundary data under a Dirichlet boundary condition, we simply return +# `flux_inner`, which is boundary data for ``\bm{\sigma}`` computed using the "inner" or interior solution. +# This way, we supply boundary data for the divergence operation without imposing any additional conditions. + +@inline function (boundary_condition::BoundaryConditionConstantDirichlet)(flux_inner, u_inner, normal::AbstractVector, + x, t, operator_type::Trixi.Divergence, + equations_parabolic::ConstantAnisotropicDiffusion2D) + return flux_inner +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) +# for an example of this). If this is done, then the boundary condition for the `Gradient` operator +# should be modified accordingly as well. +# +# ## Putting things together +# +# Finally, we can instantiate our new parabolic equation type, define boundary conditions, +# and run a simulation. The specific anisotropic diffusion matrix we use produces more +# dissipation in the direction (1, -1) than an isotropic diffusion. +# +# For boundary conditions, we impose that ``u=1`` on the left wall, ``u=2`` on the bottom +# wall, and ``u = 0`` on the outflow walls. The initial condition is taken to be ``u = 0``. +# Note that we use `BoundaryConditionConstantDirichlet` only for the parabolic boundary +# conditions, since we have not defined its behavior for the hyperbolic part. + +diffusivity = 5e-2 * SMatrix{2, 2}([2 -1; -1 2]) +equations_parabolic = ConstantAnisotropicDiffusion2D(diffusivity, equations_hyperbolic); + +u_left = 1.0 +u_bottom = 2.0 + +boundary_conditions_hyperbolic = (; x_neg = BoundaryConditionDirichlet((x, t, equations) -> SVector(u_left)), + y_neg = BoundaryConditionDirichlet((x, t, equations) -> SVector(u_bottom)), + y_pos = boundary_condition_do_nothing, + x_pos = boundary_condition_do_nothing) + +boundary_conditions_parabolic = (; x_neg = BoundaryConditionConstantDirichlet(u_left), + y_neg = BoundaryConditionConstantDirichlet(u_bottom), + y_pos = BoundaryConditionConstantDirichlet(0.0), + x_pos = BoundaryConditionConstantDirichlet(0.0)); + +solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs) +coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y)) +coordinates_max = ( 1.0, 1.0) # maximum coordinates (max(x), max(y)) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=4, + periodicity=false, n_cells_max=30_000) # set maximum capacity of tree data structure + +initial_condition = (x, t, equations) -> SVector(0.0) + +semi = SemidiscretizationHyperbolicParabolic(mesh, + (equations_hyperbolic, equations_parabolic), + initial_condition, solver; + boundary_conditions=(boundary_conditions_hyperbolic, + boundary_conditions_parabolic)) + +tspan = (0.0, 2.0) +ode = semidiscretize(semi, tspan) +callbacks = CallbackSet(SummaryCallback()) +time_int_tol = 1.0e-6 +sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol, + save_everystep=false, callback=callbacks); + +using Plots +plot(sol) + diff --git a/docs/make.jl b/docs/make.jl index 1f2791d7054..d8d1298fba6 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -40,6 +40,7 @@ files = [ "Adding a new scalar conservation law" => "adding_new_scalar_equations.jl", "Adding a non-conservative equation" => "adding_nonconservative_equation.jl", "Parabolic terms" => "parabolic_terms.jl", + "Adding new parabolic terms" => "adding_new_parabolic_terms.jl", # Topic: meshes "Adaptive mesh refinement" => "adaptive_mesh_refinement.jl", "Structured mesh with curvilinear mapping" => "structured_mesh_mapping.jl", From 8807e2ae1613ec73db0c26893e92b297d2725388 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 25 Aug 2022 13:41:59 -0500 Subject: [PATCH 3/9] import SMatrix --- docs/literate/src/files/adding_new_parabolic_terms.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/literate/src/files/adding_new_parabolic_terms.jl b/docs/literate/src/files/adding_new_parabolic_terms.jl index 761b63cb89a..9ce075db40c 100644 --- a/docs/literate/src/files/adding_new_parabolic_terms.jl +++ b/docs/literate/src/files/adding_new_parabolic_terms.jl @@ -117,6 +117,7 @@ end # Note that we use `BoundaryConditionConstantDirichlet` only for the parabolic boundary # conditions, since we have not defined its behavior for the hyperbolic part. +using StaticArrays: SMatrix diffusivity = 5e-2 * SMatrix{2, 2}([2 -1; -1 2]) equations_parabolic = ConstantAnisotropicDiffusion2D(diffusivity, equations_hyperbolic); From e3b9b45e099d6d0c42bbc8e751a0f7ad3fa7935b Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 25 Aug 2022 15:06:22 -0500 Subject: [PATCH 4/9] fix using Trixi --- docs/literate/src/files/adding_new_parabolic_terms.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/literate/src/files/adding_new_parabolic_terms.jl b/docs/literate/src/files/adding_new_parabolic_terms.jl index 9ce075db40c..a6068edbf11 100644 --- a/docs/literate/src/files/adding_new_parabolic_terms.jl +++ b/docs/literate/src/files/adding_new_parabolic_terms.jl @@ -117,7 +117,7 @@ end # Note that we use `BoundaryConditionConstantDirichlet` only for the parabolic boundary # conditions, since we have not defined its behavior for the hyperbolic part. -using StaticArrays: SMatrix +using Trixi: SMatrix diffusivity = 5e-2 * SMatrix{2, 2}([2 -1; -1 2]) equations_parabolic = ConstantAnisotropicDiffusion2D(diffusivity, equations_hyperbolic); From 10d84fdd387c575459836fd03ed9378b8ff1209c Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 25 Aug 2022 15:36:27 -0500 Subject: [PATCH 5/9] add extra header --- docs/literate/src/files/adding_new_parabolic_terms.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/literate/src/files/adding_new_parabolic_terms.jl b/docs/literate/src/files/adding_new_parabolic_terms.jl index a6068edbf11..8f55ffc37ab 100644 --- a/docs/literate/src/files/adding_new_parabolic_terms.jl +++ b/docs/literate/src/files/adding_new_parabolic_terms.jl @@ -11,6 +11,8 @@ using Trixi advection_velocity = (1.0, 1.0) equations_hyperbolic = LinearScalarAdvectionEquation2D(advection_velocity); +# ## Define a new parabolic equation type +# # Next, we define a 2D parabolic diffusion term type. This is similar to [`LaplaceDiffusion2D`](@ref) # except that the `diffusivity` field refers to a spatially constant diffusivity matrix now. Note that # `ConstantAnisotropicDiffusion2D` has a field for `equations_hyperbolic`. It is useful to have From 82f5e6b713f229d65b83b4ec36af1c09c608b167 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 25 Aug 2022 20:46:35 -0500 Subject: [PATCH 6/9] fix some math --- docs/literate/src/files/adding_new_parabolic_terms.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/literate/src/files/adding_new_parabolic_terms.jl b/docs/literate/src/files/adding_new_parabolic_terms.jl index 8f55ffc37ab..def56d6bdcf 100644 --- a/docs/literate/src/files/adding_new_parabolic_terms.jl +++ b/docs/literate/src/files/adding_new_parabolic_terms.jl @@ -33,7 +33,7 @@ varnames(variable_mapping, equations_parabolic::ConstantAnisotropicDiffusion2D) # \partial_t u(t,x) + \partial_x (f_x(u) - g_x(u, \nabla u)) # + \partial_y (f_y(u) - g_y(u, \nabla u)) = 0 # ``` -# where `f_x(u)`, `f_y(u)` are the hyperbolic fluxes and `g_x(u, \nabla u)`, `g_y(u, \nabla u)` denote +# where ``f_x(u)``, ``f_y(u)`` are the hyperbolic fluxes and ``g_x(u, \nabla u)``, ``g_y(u, \nabla u)`` denote # the viscous fluxes. For anisotropic diffusion, the viscous fluxes are the first and second components # of the matrix-vector product involving `diffusivity` and the gradient vector. # @@ -56,8 +56,8 @@ end # ```math # \begin{aligned} # \bm{q} &= \nabla u \\ -# \bm{sigma} &= \begin{pmatrix} g_x(u, \bm{q}) \\ g_y(u, \bm{q}) \end{pmatrix} \\ -# \text{viscous contribution } = \nabla \cdot \bm{\sigma} +# \bm{\sigma} &= \begin{pmatrix} g_x(u, \bm{q}) \\ g_y(u, \bm{q}) \end{pmatrix} \\ +# \text{viscous contribution } &= \nabla \cdot \bm{\sigma} # \end{aligned} # ``` # From 3f034db80559a4ed96efed4f278c79c9d8950679 Mon Sep 17 00:00:00 2001 From: Jesse Chan <1156048+jlchan@users.noreply.github.com> Date: Fri, 26 Aug 2022 08:33:45 -0500 Subject: [PATCH 7/9] Apply suggestions from code review Co-authored-by: Andrew Winters Co-authored-by: Michael Schlottke-Lakemper --- docs/literate/src/files/adding_new_parabolic_terms.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/docs/literate/src/files/adding_new_parabolic_terms.jl b/docs/literate/src/files/adding_new_parabolic_terms.jl index def56d6bdcf..af9cd996faa 100644 --- a/docs/literate/src/files/adding_new_parabolic_terms.jl +++ b/docs/literate/src/files/adding_new_parabolic_terms.jl @@ -33,7 +33,7 @@ varnames(variable_mapping, equations_parabolic::ConstantAnisotropicDiffusion2D) # \partial_t u(t,x) + \partial_x (f_x(u) - g_x(u, \nabla u)) # + \partial_y (f_y(u) - g_y(u, \nabla u)) = 0 # ``` -# where ``f_x(u)``, ``f_y(u)`` are the hyperbolic fluxes and ``g_x(u, \nabla u)``, ``g_y(u, \nabla u)`` denote +# 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 # of the matrix-vector product involving `diffusivity` and the gradient vector. # @@ -56,7 +56,7 @@ end # ```math # \begin{aligned} # \bm{q} &= \nabla u \\ -# \bm{\sigma} &= \begin{pmatrix} g_x(u, \bm{q}) \\ g_y(u, \bm{q}) \end{pmatrix} \\ +# \bm{\sigma} &= \begin{pmatrix} g_1(u, \bm{q}) \\ g_2(u, \bm{q}) \end{pmatrix} \\ # \text{viscous contribution } &= \nabla \cdot \bm{\sigma} # \end{aligned} # ``` @@ -75,7 +75,7 @@ end # This boundary condition contains only the field `boundary_value`, which we assume to be some # real-valued constant which we will impose as the Dirichlet data on the boundary. # -# Boundary conditions have generally been defined as "callable structs" (also known as "functors). +# Boundary conditions have generally been defined as "callable structs" (also known as "functors"). # For each boundary condition, we need to specify the appropriate boundary data to return for both # the `Gradient` and `Divergence`. Since the gradient is operating on the solution `u`, the boundary # data should be the value of `u`, and we can directly impose Dirichlet data. @@ -112,7 +112,7 @@ end # # Finally, we can instantiate our new parabolic equation type, define boundary conditions, # and run a simulation. The specific anisotropic diffusion matrix we use produces more -# dissipation in the direction (1, -1) than an isotropic diffusion. +# dissipation in the direction ``(1, -1)`` as an isotropic diffusion. # # For boundary conditions, we impose that ``u=1`` on the left wall, ``u=2`` on the bottom # wall, and ``u = 0`` on the outflow walls. The initial condition is taken to be ``u = 0``. @@ -120,7 +120,7 @@ end # conditions, since we have not defined its behavior for the hyperbolic part. using Trixi: SMatrix -diffusivity = 5e-2 * SMatrix{2, 2}([2 -1; -1 2]) +diffusivity = 5.0e-2 * SMatrix{2, 2}([2 -1; -1 2]) equations_parabolic = ConstantAnisotropicDiffusion2D(diffusivity, equations_hyperbolic); u_left = 1.0 From f03e6eaa0a1d3e402c9b2753d08679dc32b867bf Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 26 Aug 2022 08:35:31 -0500 Subject: [PATCH 8/9] removing u_left, u_bottom --- docs/literate/src/files/adding_new_parabolic_terms.jl | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/docs/literate/src/files/adding_new_parabolic_terms.jl b/docs/literate/src/files/adding_new_parabolic_terms.jl index def56d6bdcf..be7ae8b2f2d 100644 --- a/docs/literate/src/files/adding_new_parabolic_terms.jl +++ b/docs/literate/src/files/adding_new_parabolic_terms.jl @@ -123,16 +123,13 @@ using Trixi: SMatrix diffusivity = 5e-2 * SMatrix{2, 2}([2 -1; -1 2]) equations_parabolic = ConstantAnisotropicDiffusion2D(diffusivity, equations_hyperbolic); -u_left = 1.0 -u_bottom = 2.0 - -boundary_conditions_hyperbolic = (; x_neg = BoundaryConditionDirichlet((x, t, equations) -> SVector(u_left)), - y_neg = BoundaryConditionDirichlet((x, t, equations) -> SVector(u_bottom)), +boundary_conditions_hyperbolic = (; x_neg = BoundaryConditionDirichlet((x, t, equations) -> SVector(1.0)), + y_neg = BoundaryConditionDirichlet((x, t, equations) -> SVector(2.0)), y_pos = boundary_condition_do_nothing, x_pos = boundary_condition_do_nothing) -boundary_conditions_parabolic = (; x_neg = BoundaryConditionConstantDirichlet(u_left), - y_neg = BoundaryConditionConstantDirichlet(u_bottom), +boundary_conditions_parabolic = (; x_neg = BoundaryConditionConstantDirichlet(1.0), + y_neg = BoundaryConditionConstantDirichlet(2.0), y_pos = BoundaryConditionConstantDirichlet(0.0), x_pos = BoundaryConditionConstantDirichlet(0.0)); From 2089716ea88edde56c400b6c9ee8fc743216ea75 Mon Sep 17 00:00:00 2001 From: Jesse Chan <1156048+jlchan@users.noreply.github.com> Date: Fri, 26 Aug 2022 09:14:38 -0500 Subject: [PATCH 9/9] Update docs/literate/src/files/adding_new_parabolic_terms.jl Co-authored-by: Andrew Winters --- docs/literate/src/files/adding_new_parabolic_terms.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/literate/src/files/adding_new_parabolic_terms.jl b/docs/literate/src/files/adding_new_parabolic_terms.jl index 12359853893..a6cecb351be 100644 --- a/docs/literate/src/files/adding_new_parabolic_terms.jl +++ b/docs/literate/src/files/adding_new_parabolic_terms.jl @@ -30,8 +30,8 @@ varnames(variable_mapping, equations_parabolic::ConstantAnisotropicDiffusion2D) # Next, we define the viscous flux function. We assume that the mixed hyperbolic-parabolic system # is of the form # ```math -# \partial_t u(t,x) + \partial_x (f_x(u) - g_x(u, \nabla u)) -# + \partial_y (f_y(u) - g_y(u, \nabla u)) = 0 +# \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