From 26afcbea2bd3c588724b316b152f92880242f240 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 17 Apr 2023 14:56:18 -0500 Subject: [PATCH 01/36] factor volume integral out of calc_gradient! and calc_divergence! --- src/solvers/dgmulti/dg_parabolic.jl | 45 ++++++++++++++++++----------- 1 file changed, 28 insertions(+), 17 deletions(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index 50cfd8ab17d..740c7d38691 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -71,15 +71,10 @@ function calc_gradient_surface_integral(gradients, u, scalar_flux_face_values, end end -function calc_gradient!(gradients, u::StructArray, t, mesh::DGMultiMesh, - equations::AbstractEquationsParabolic, - boundary_conditions, dg::DGMulti, cache, cache_parabolic) +function calc_gradient_volume_integral!(gradients, u, mesh::DGMultiMesh, equations, + dg::DGMulti, cache, cache_parabolic) - @unpack weak_differentiation_matrices = cache_parabolic - - for dim in eachindex(gradients) - reset_du!(gradients[dim], dg) - end + (; weak_differentiation_matrices) = cache_parabolic # compute volume contributions to gradients @threaded for e in eachelement(mesh, dg) @@ -89,6 +84,17 @@ function calc_gradient!(gradients, u::StructArray, t, mesh::DGMultiMesh, view(gradients[i], :, e), view(u, :, e)) end end +end + +function calc_gradient!(gradients, u::StructArray, t, mesh::DGMultiMesh, + equations::AbstractEquationsParabolic, + boundary_conditions, dg::DGMulti, cache, cache_parabolic) + + for dim in eachindex(gradients) + reset_du!(gradients[dim], dg) + end + + calc_gradient_volume_integral!(gradients, u, mesh, equations, dg, cache, cache_parabolic) @unpack u_face_values = cache_parabolic prolong2interfaces!(u_face_values, u, mesh, equations, dg.surface_integral, dg, cache) @@ -189,7 +195,8 @@ function calc_viscous_fluxes!(flux_viscous, u, gradients, mesh::DGMultiMesh, end # interpolate u and gradient to quadrature points, store in `local_flux_viscous` - apply_to_each_field(mul_by!(dg.basis.Vq), local_u_values, view(u, :, e)) # TODO: DGMulti. Specialize for nodal collocation methods (SBP, GaussSBP) + # TODO: DGMulti. Specialize for nodal collocation methods (SBP, GaussSBP)? + apply_to_each_field(mul_by!(dg.basis.Vq), local_u_values, view(u, :, e)) for dim in eachdim(mesh) apply_to_each_field(mul_by!(dg.basis.Vq), local_flux_viscous[dim], view(gradients[dim], :, e)) end @@ -233,14 +240,9 @@ function calc_viscous_penalty!(scalar_flux_face_values, u_face_values, t, bounda return nothing end - -function calc_divergence!(du, u::StructArray, t, flux_viscous, mesh::DGMultiMesh, - equations::AbstractEquationsParabolic, - boundary_conditions, dg::DGMulti, parabolic_scheme, cache, cache_parabolic) - - @unpack weak_differentiation_matrices = cache_parabolic - - reset_du!(du, dg) +function calc_divergence_volume_integral!(du, u, flux_viscous, mesh::DGMultiMesh, equations, + dg::DGMulti, cache, cache_parabolic) + (; weak_differentiation_matrices) = cache_parabolic # compute volume contributions to divergence @threaded for e in eachelement(mesh, dg) @@ -250,6 +252,15 @@ function calc_divergence!(du, u::StructArray, t, flux_viscous, mesh::DGMultiMesh view(du, :, e), view(flux_viscous[i], :, e)) end end +end + +function calc_divergence!(du, u::StructArray, t, flux_viscous, mesh::DGMultiMesh, + equations::AbstractEquationsParabolic, + boundary_conditions, dg::DGMulti, parabolic_scheme, cache, cache_parabolic) + + reset_du!(du, dg) + + calc_divergence_volume_integral!(du, u, flux_viscous, mesh, equations, dg, cache, cache_parabolic) # interpolates from solution coefficients to face quadrature points flux_viscous_face_values = cache_parabolic.gradients_face_values # reuse storage From 2f5732feb57082b3cc08d8ad1f54006ba4291fcb Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 17 Apr 2023 16:29:00 -0500 Subject: [PATCH 02/36] add "!" for consistency --- src/solvers/dgmulti/dg_parabolic.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index 740c7d38691..f2c441b445d 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -55,9 +55,9 @@ function prolong2interfaces!(u_face_values, u, mesh::DGMultiMesh, equations::Abs apply_to_each_field(mul_by!(dg.basis.Vf), u_face_values, u) end -function calc_gradient_surface_integral(gradients, u, scalar_flux_face_values, - mesh, equations::AbstractEquationsParabolic, - dg::DGMulti, cache, cache_parabolic) +function calc_gradient_surface_integral!(gradients, u, scalar_flux_face_values, + mesh, equations::AbstractEquationsParabolic, + dg::DGMulti, cache, cache_parabolic) @unpack local_flux_face_values_threaded = cache_parabolic @threaded for e in eachelement(mesh, dg) local_flux_values = local_flux_face_values_threaded[Threads.threadid()] @@ -113,8 +113,8 @@ function calc_gradient!(gradients, u::StructArray, t, mesh::DGMultiMesh, mesh, equations, dg, cache, cache_parabolic) # compute surface contributions - calc_gradient_surface_integral(gradients, u, scalar_flux_face_values, - mesh, equations, dg, cache, cache_parabolic) + calc_gradient_surface_integral!(gradients, u, scalar_flux_face_values, + mesh, equations, dg, cache, cache_parabolic) for dim in eachdim(mesh) invert_jacobian!(gradients[dim], mesh, equations, dg, cache; scaling=1.0) From 20df8de66075042f52d46ba0ac72ae697d4d1fc9 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 17 Apr 2023 20:44:11 -0500 Subject: [PATCH 03/36] minor cleanup --- src/solvers/dgmulti/dg_parabolic.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index f2c441b445d..3c6ee2b1995 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -6,11 +6,11 @@ function create_cache_parabolic(mesh::DGMultiMesh, # TODO: parabolic; utilize the parabolic variables in `equations_parabolic` to reduce memory usage in the parabolic cache nvars = nvariables(equations_hyperbolic) - @unpack M, Drst = dg.basis + (; M, Drst) = dg.basis weak_differentiation_matrices = map(A -> -M \ (A' * M), Drst) # u_transformed stores "transformed" variables for computing the gradient - @unpack md = mesh + (; md) = mesh u_transformed = allocate_nested_array(uEltype, nvars, size(md.x), dg) gradients = ntuple(_ -> similar(u_transformed), ndims(mesh)) flux_viscous = similar.(gradients) @@ -71,7 +71,8 @@ function calc_gradient_surface_integral!(gradients, u, scalar_flux_face_values, end end -function calc_gradient_volume_integral!(gradients, u, mesh::DGMultiMesh, equations, +function calc_gradient_volume_integral!(gradients, u, mesh::DGMultiMesh, + equations::AbstractEquationsParabolic, dg::DGMulti, cache, cache_parabolic) (; weak_differentiation_matrices) = cache_parabolic @@ -240,7 +241,7 @@ function calc_viscous_penalty!(scalar_flux_face_values, u_face_values, t, bounda return nothing end -function calc_divergence_volume_integral!(du, u, flux_viscous, mesh::DGMultiMesh, equations, +function calc_divergence_volume_integral!(du, u, flux_viscous, mesh::DGMultiMesh, dg::DGMulti, cache, cache_parabolic) (; weak_differentiation_matrices) = cache_parabolic From e689a724a143350d638ee92e0fbb9850775fa9ee Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 17 Apr 2023 20:52:03 -0500 Subject: [PATCH 04/36] remove inv_h from penalty(...) - it's not necessary for coercivity (BR1 is coercive with any penalty parameter > 0; a factor of 1/h is only necessary for IPDG) - simplifies the implementation for curved meshes - it's not in CI tests yet (should be added soon) --- src/equations/laplace_diffusion_2d.jl | 4 ++-- src/solvers/dgmulti/dg_parabolic.jl | 18 +++++------------- 2 files changed, 7 insertions(+), 15 deletions(-) diff --git a/src/equations/laplace_diffusion_2d.jl b/src/equations/laplace_diffusion_2d.jl index 2f1afe25a6d..1dc843ca6c7 100644 --- a/src/equations/laplace_diffusion_2d.jl +++ b/src/equations/laplace_diffusion_2d.jl @@ -28,8 +28,8 @@ end # TODO: parabolic; should this remain in the equations file, be moved to solvers, or live in the elixir? # 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) - return dg.penalty_parameter * (u_outer - u_inner) * equations_parabolic.diffusivity * inv_h +function penalty(u_outer, u_inner, equations_parabolic::LaplaceDiffusion2D, dg::ViscousFormulationLocalDG) + return dg.penalty_parameter * (u_outer - u_inner) * equations_parabolic.diffusivity end # Dirichlet-type boundary condition for use with a parabolic solver in weak form diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index 3c6ee2b1995..f35f5063798 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -2,6 +2,7 @@ function create_cache_parabolic(mesh::DGMultiMesh, equations_hyperbolic::AbstractEquations, equations_parabolic::AbstractEquationsParabolic, dg::DGMulti, parabolic_scheme, RealT, uEltype) + # default to taking derivatives of all hyperbolic terms # TODO: parabolic; utilize the parabolic variables in `equations_parabolic` to reduce memory usage in the parabolic cache nvars = nvariables(equations_hyperbolic) @@ -23,17 +24,8 @@ function create_cache_parabolic(mesh::DGMultiMesh, local_flux_viscous_threaded = [ntuple(_ -> similar(u_transformed, dg.basis.Nq), ndims(mesh)) for _ in 1:Threads.nthreads()] local_flux_face_values_threaded = [similar(scalar_flux_face_values[:, 1]) for _ in 1:Threads.nthreads()] - # precompute 1 / h for penalty terms - inv_h = similar(mesh.md.Jf) - J = dg.basis.Vf * mesh.md.J # interp to face nodes - for e in eachelement(mesh, dg) - for i in each_face_node(mesh, dg) - inv_h[i, e] = mesh.md.Jf[i, e] / J[i, e] - end - end - return (; u_transformed, gradients, flux_viscous, - weak_differentiation_matrices, inv_h, + weak_differentiation_matrices, u_face_values, gradients_face_values, scalar_flux_face_values, local_u_values_threaded, local_flux_viscous_threaded, local_flux_face_values_threaded) end @@ -230,18 +222,18 @@ function calc_viscous_penalty!(scalar_flux_face_values, u_face_values, t, bounda mesh, equations::AbstractEquationsParabolic, dg::DGMulti, parabolic_scheme, cache, cache_parabolic) # compute fluxes at interfaces - @unpack scalar_flux_face_values, inv_h = cache_parabolic + @unpack scalar_flux_face_values = cache_parabolic @unpack mapM, mapP = mesh.md @threaded for face_node_index in each_face_node_global(mesh, dg) idM, idP = mapM[face_node_index], mapP[face_node_index] uM, uP = u_face_values[idM], u_face_values[idP] - inv_h_face = inv_h[face_node_index] - scalar_flux_face_values[idM] = scalar_flux_face_values[idM] + penalty(uP, uM, inv_h_face, equations, parabolic_scheme) + scalar_flux_face_values[idM] = scalar_flux_face_values[idM] + penalty(uP, uM, equations, parabolic_scheme) end return nothing end function calc_divergence_volume_integral!(du, u, flux_viscous, mesh::DGMultiMesh, + equations::AbstractEquationsParabolic, dg::DGMulti, cache, cache_parabolic) (; weak_differentiation_matrices) = cache_parabolic From 8b914def333d9893e744676b71c3eacc0ea9affb Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 24 Apr 2023 16:11:15 -0500 Subject: [PATCH 05/36] generalized invert_jacobian --- src/solvers/dgmulti/dg.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl index dd6d9f43363..f9e30f8f871 100644 --- a/src/solvers/dgmulti/dg.jl +++ b/src/solvers/dgmulti/dg.jl @@ -489,8 +489,9 @@ end # inverts Jacobian and scales by -1.0 function invert_jacobian!(du, mesh::DGMultiMesh, equations, dg::DGMulti, cache; scaling=-1) @threaded for e in eachelement(mesh, dg, cache) + invJ = cache.invJ[1, e] for i in axes(du, 1) - du[i, e] *= scaling * cache.invJ[i, e] + du[i, e] *= scaling * invJ end end end From 9d90992399e7bf0e6906f9c40b193ea05771a1b4 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 24 Apr 2023 16:12:41 -0500 Subject: [PATCH 06/36] remove "prolong2interfaces" the call isn't that complex, and it's easier for me to debug --- src/solvers/dgmulti/dg_parabolic.jl | 14 ++++---------- 1 file changed, 4 insertions(+), 10 deletions(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index f35f5063798..5d3a12d02b9 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -40,12 +40,7 @@ function transform_variables!(u_transformed, u, mesh, equations_parabolic::Abstr end end -# interpolates from solution coefficients to face quadrature points -# We pass the `surface_integral` argument solely for dispatch -function prolong2interfaces!(u_face_values, u, mesh::DGMultiMesh, equations::AbstractEquationsParabolic, - surface_integral, dg::DGMulti, cache) - apply_to_each_field(mul_by!(dg.basis.Vf), u_face_values, u) -end +# TODO: reuse entropy projection computations for DGMultiFluxDiff{<:Polynomial} (including `GaussSBP` solvers) function calc_gradient_surface_integral!(gradients, u, scalar_flux_face_values, mesh, equations::AbstractEquationsParabolic, @@ -89,8 +84,8 @@ function calc_gradient!(gradients, u::StructArray, t, mesh::DGMultiMesh, calc_gradient_volume_integral!(gradients, u, mesh, equations, dg, cache, cache_parabolic) - @unpack u_face_values = cache_parabolic - prolong2interfaces!(u_face_values, u, mesh, equations, dg.surface_integral, dg, cache) + (; u_face_values) = cache_parabolic + apply_to_each_field(mul_by!(dg.basis.Vf), u_face_values, u) # compute fluxes at interfaces @unpack scalar_flux_face_values = cache_parabolic @@ -258,8 +253,7 @@ function calc_divergence!(du, u::StructArray, t, flux_viscous, mesh::DGMultiMesh # interpolates from solution coefficients to face quadrature points flux_viscous_face_values = cache_parabolic.gradients_face_values # reuse storage for dim in eachdim(mesh) - prolong2interfaces!(flux_viscous_face_values[dim], flux_viscous[dim], mesh, equations, - dg.surface_integral, dg, cache) + apply_to_each_field(mul_by!(dg.basis.Vf), flux_viscous_face_values[dim], flux_viscous[dim]) end # compute fluxes at interfaces From 4b032ecc274531d49783f15a4497a4a2898b5b28 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 24 Apr 2023 16:13:16 -0500 Subject: [PATCH 07/36] switch to "strong-weak" formulation --- src/solvers/dgmulti/dg_parabolic.jl | 68 +++++++++++++++++------------ 1 file changed, 41 insertions(+), 27 deletions(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index 5d3a12d02b9..eac164a62b0 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -1,18 +1,22 @@ +# version for standard (e.g., non-entropy stable or flux differencing) schemes function create_cache_parabolic(mesh::DGMultiMesh, equations_hyperbolic::AbstractEquations, equations_parabolic::AbstractEquationsParabolic, dg::DGMulti, parabolic_scheme, RealT, uEltype) - # default to taking derivatives of all hyperbolic terms + # default to taking derivatives of all hyperbolic variables # TODO: parabolic; utilize the parabolic variables in `equations_parabolic` to reduce memory usage in the parabolic cache nvars = nvariables(equations_hyperbolic) - (; M, Drst) = dg.basis - weak_differentiation_matrices = map(A -> -M \ (A' * M), Drst) + (; M, Vq, Pq, Drst) = dg.basis + weak_differentiation_matrices = map(A -> (M \ (-A' * M)), Drst) + strong_differentiation_matrices = Drst + lift_matrix = dg.basis.LIFT # u_transformed stores "transformed" variables for computing the gradient (; md) = mesh u_transformed = allocate_nested_array(uEltype, nvars, size(md.x), dg) + #gradients = ntuple(_ -> similar(u_transformed, (dg.basis.Nq, mesh.md.num_elements)), ndims(mesh)) gradients = ntuple(_ -> similar(u_transformed), ndims(mesh)) flux_viscous = similar.(gradients) @@ -25,7 +29,7 @@ function create_cache_parabolic(mesh::DGMultiMesh, local_flux_face_values_threaded = [similar(scalar_flux_face_values[:, 1]) for _ in 1:Threads.nthreads()] return (; u_transformed, gradients, flux_viscous, - weak_differentiation_matrices, + weak_differentiation_matrices, strong_differentiation_matrices, lift_matrix, u_face_values, gradients_face_values, scalar_flux_face_values, local_u_values_threaded, local_flux_viscous_threaded, local_flux_face_values_threaded) end @@ -45,7 +49,7 @@ end function calc_gradient_surface_integral!(gradients, u, scalar_flux_face_values, mesh, equations::AbstractEquationsParabolic, dg::DGMulti, cache, cache_parabolic) - @unpack local_flux_face_values_threaded = cache_parabolic + @unpack lift_matrix, local_flux_face_values_threaded = cache_parabolic @threaded for e in eachelement(mesh, dg) local_flux_values = local_flux_face_values_threaded[Threads.threadid()] for dim in eachdim(mesh) @@ -53,7 +57,7 @@ function calc_gradient_surface_integral!(gradients, u, scalar_flux_face_values, # compute flux * (nx, ny, nz) local_flux_values[i] = scalar_flux_face_values[i, e] * mesh.md.nxyzJ[dim][i, e] end - apply_to_each_field(mul_by_accum!(dg.basis.LIFT), view(gradients[dim], :, e), local_flux_values) + apply_to_each_field(mul_by_accum!(lift_matrix), view(gradients[dim], :, e), local_flux_values) end end end @@ -62,13 +66,13 @@ function calc_gradient_volume_integral!(gradients, u, mesh::DGMultiMesh, equations::AbstractEquationsParabolic, dg::DGMulti, cache, cache_parabolic) - (; weak_differentiation_matrices) = cache_parabolic + (; strong_differentiation_matrices) = cache_parabolic # compute volume contributions to gradients @threaded for e in eachelement(mesh, dg) for i in eachdim(mesh), j in eachdim(mesh) dxidxhatj = mesh.md.rstxyzJ[i, j][1, e] # TODO: DGMulti. Assumes mesh is affine here. - apply_to_each_field(mul_by_accum!(weak_differentiation_matrices[j], dxidxhatj), + apply_to_each_field(mul_by_accum!(strong_differentiation_matrices[j], dxidxhatj), view(gradients[i], :, e), view(u, :, e)) end end @@ -88,13 +92,13 @@ function calc_gradient!(gradients, u::StructArray, t, mesh::DGMultiMesh, apply_to_each_field(mul_by!(dg.basis.Vf), u_face_values, u) # compute fluxes at interfaces - @unpack scalar_flux_face_values = cache_parabolic - @unpack mapM, mapP, Jf = mesh.md + (; scalar_flux_face_values) = cache_parabolic + (; mapM, mapP, Jf) = mesh.md @threaded for face_node_index in each_face_node_global(mesh, dg) idM, idP = mapM[face_node_index], mapP[face_node_index] uM = u_face_values[idM] uP = u_face_values[idP] - scalar_flux_face_values[idM] = 0.5 * (uP + uM) # TODO: use strong/weak formulation for curved meshes? + scalar_flux_face_values[idM] = 0.5 * (uP - uM) # TODO: use strong/weak formulation for curved meshes? end calc_boundary_flux!(scalar_flux_face_values, u_face_values, t, Gradient(), boundary_conditions, @@ -105,7 +109,12 @@ function calc_gradient!(gradients, u::StructArray, t, mesh::DGMultiMesh, mesh, equations, dg, cache, cache_parabolic) for dim in eachdim(mesh) - invert_jacobian!(gradients[dim], mesh, equations, dg, cache; scaling=1.0) + @threaded for e in eachelement(mesh, dg) + invJ = cache.invJ[1, e] + for i in axes(gradients[dim], 1) + gradients[dim][i, e] = gradients[dim][i, e] * invJ + end + end end end @@ -141,7 +150,7 @@ function calc_single_boundary_flux!(flux_face_values, u_face_values, t, md = mesh.md num_pts_per_face = rd.Nfq ÷ rd.Nfaces - @unpack xyzf, nxyzJ, Jf = md + @unpack xyzf, nxyz = md for f in mesh.boundary_faces[boundary_key] for i in Base.OneTo(num_pts_per_face) @@ -149,7 +158,7 @@ function calc_single_boundary_flux!(flux_face_values, u_face_values, t, e = ((f-1) ÷ rd.Nfaces) + 1 fid = i + ((f-1) % rd.Nfaces) * num_pts_per_face - face_normal = SVector{NDIMS}(getindex.(nxyzJ, fid, e)) / Jf[fid,e] + face_normal = SVector{NDIMS}(getindex.(nxyz, fid, e)) face_coordinates = SVector{NDIMS}(getindex.(xyzf, fid, e)) # for both the gradient and the divergence, the boundary flux is scalar valued. @@ -157,6 +166,12 @@ function calc_single_boundary_flux!(flux_face_values, u_face_values, t, flux_face_values[fid,e] = boundary_condition(flux_face_values[fid,e], u_face_values[fid,e], face_normal, face_coordinates, t, operator_type, equations) + + # use "strong form" for the Gradient ("weak form" for Divergence) + if operator_type isa Gradient + flux_face_values[fid, e] = flux_face_values[fid, e] - u_face_values[fid, e] + end + end end return nothing @@ -174,18 +189,15 @@ function calc_viscous_fluxes!(flux_viscous, u, gradients, mesh::DGMultiMesh, @threaded for e in eachelement(mesh, dg) - # reset local storage for each element - local_flux_viscous = local_flux_viscous_threaded[Threads.threadid()] + # reset local storage for each element, interpolate u to quadrature points + # TODO: DGMulti. Specialize for nodal collocation methods (SBP, GaussSBP)? local_u_values = local_u_values_threaded[Threads.threadid()] fill!(local_u_values, zero(eltype(local_u_values))) - for dim in eachdim(mesh) - fill!(local_flux_viscous[dim], zero(eltype(local_flux_viscous[dim]))) - end - - # interpolate u and gradient to quadrature points, store in `local_flux_viscous` - # TODO: DGMulti. Specialize for nodal collocation methods (SBP, GaussSBP)? apply_to_each_field(mul_by!(dg.basis.Vq), local_u_values, view(u, :, e)) + + local_flux_viscous = local_flux_viscous_threaded[Threads.threadid()] for dim in eachdim(mesh) + fill!(local_flux_viscous[dim], zero(eltype(local_flux_viscous[dim]))) apply_to_each_field(mul_by!(dg.basis.Vq), local_flux_viscous[dim], view(gradients[dim], :, e)) end @@ -257,18 +269,19 @@ function calc_divergence!(du, u::StructArray, t, flux_viscous, mesh::DGMultiMesh end # compute fluxes at interfaces - @unpack scalar_flux_face_values = cache_parabolic - @unpack mapM, mapP, nxyzJ = mesh.md + (; scalar_flux_face_values) = cache_parabolic + (; mapM, mapP, nxyzJ) = mesh.md + @threaded for face_node_index in each_face_node_global(mesh, dg, cache, cache_parabolic) idM, idP = mapM[face_node_index], mapP[face_node_index] # compute f(u, ∇u) ⋅ n flux_face_value = zero(eltype(scalar_flux_face_values)) for dim in eachdim(mesh) - uM = flux_viscous_face_values[dim][idM] - uP = flux_viscous_face_values[dim][idP] + fM = flux_viscous_face_values[dim][idM] + fP = flux_viscous_face_values[dim][idP] # TODO: use strong/weak formulation to ensure stability on curved meshes? - flux_face_value = flux_face_value + 0.5 * (uP + uM) * nxyzJ[dim][face_node_index] + flux_face_value = flux_face_value + 0.5 * (fP + fM) * nxyzJ[dim][face_node_index] end scalar_flux_face_values[idM] = flux_face_value end @@ -276,6 +289,7 @@ function calc_divergence!(du, u::StructArray, t, flux_viscous, mesh::DGMultiMesh calc_boundary_flux!(scalar_flux_face_values, cache_parabolic.u_face_values, t, 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) From 5a2384e8f9dcdde4d4642bc3c53d42f71e6b5dc7 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 24 Apr 2023 16:39:56 -0500 Subject: [PATCH 08/36] precompute divergence interp matrix --- src/solvers/dgmulti/dg_parabolic.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index eac164a62b0..0886ee28c11 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -12,6 +12,7 @@ function create_cache_parabolic(mesh::DGMultiMesh, weak_differentiation_matrices = map(A -> (M \ (-A' * M)), Drst) strong_differentiation_matrices = Drst lift_matrix = dg.basis.LIFT + projection_face_interpolation_matrix = dg.basis.Vf * dg.basis.Pq # u_transformed stores "transformed" variables for computing the gradient (; md) = mesh @@ -263,9 +264,11 @@ function calc_divergence!(du, u::StructArray, t, flux_viscous, mesh::DGMultiMesh calc_divergence_volume_integral!(du, u, flux_viscous, 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 for dim in eachdim(mesh) - apply_to_each_field(mul_by!(dg.basis.Vf), flux_viscous_face_values[dim], flux_viscous[dim]) + # TODO: optimize this!! + apply_to_each_field(mul_by!(projection_face_interpolation_matrix), flux_viscous_face_values[dim], flux_viscous[dim]) end # compute fluxes at interfaces From a6931aea257a0f9d5733f0bb5455e0126e1671da Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 24 Apr 2023 16:40:29 -0500 Subject: [PATCH 09/36] specialize differentiation matrices for Gradient vs Divergence --- src/solvers/dgmulti/dg_parabolic.jl | 46 +++++++++++++---------------- 1 file changed, 20 insertions(+), 26 deletions(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index 0886ee28c11..6d89455caef 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -2,35 +2,40 @@ function create_cache_parabolic(mesh::DGMultiMesh, equations_hyperbolic::AbstractEquations, equations_parabolic::AbstractEquationsParabolic, - dg::DGMulti, parabolic_scheme, RealT, uEltype) + dg::DGMulti{NDIMS}, parabolic_scheme, RealT, uEltype) where {NDIMS} # default to taking derivatives of all hyperbolic variables # TODO: parabolic; utilize the parabolic variables in `equations_parabolic` to reduce memory usage in the parabolic cache nvars = nvariables(equations_hyperbolic) (; M, Vq, Pq, Drst) = dg.basis - weak_differentiation_matrices = map(A -> (M \ (-A' * M)), Drst) - strong_differentiation_matrices = Drst - lift_matrix = dg.basis.LIFT + + # gradient operators: map from nodes to quadrature + strong_differentiation_matrices = map(A -> Vq * A, Drst) + gradient_lift_matrix = Vq * dg.basis.LIFT + + # divergence operators: map from quadrature to nodes + weak_differentiation_matrices = map(A -> (M \ (-A' * M * Pq)), Drst) + divergence_lift_matrix = dg.basis.LIFT projection_face_interpolation_matrix = dg.basis.Vf * dg.basis.Pq # u_transformed stores "transformed" variables for computing the gradient (; md) = mesh u_transformed = allocate_nested_array(uEltype, nvars, size(md.x), dg) - #gradients = ntuple(_ -> similar(u_transformed, (dg.basis.Nq, mesh.md.num_elements)), ndims(mesh)) - gradients = ntuple(_ -> similar(u_transformed), ndims(mesh)) + gradients = SVector{NDIMS}(ntuple(_ -> similar(u_transformed, (dg.basis.Nq, mesh.md.num_elements)), NDIMS)) flux_viscous = similar.(gradients) u_face_values = allocate_nested_array(uEltype, nvars, size(md.xf), dg) scalar_flux_face_values = similar(u_face_values) - gradients_face_values = ntuple(_ -> similar(u_face_values), ndims(mesh)) + gradients_face_values = ntuple(_ -> similar(u_face_values), NDIMS) local_u_values_threaded = [similar(u_transformed, dg.basis.Nq) for _ in 1:Threads.nthreads()] - local_flux_viscous_threaded = [ntuple(_ -> similar(u_transformed, dg.basis.Nq), ndims(mesh)) for _ in 1:Threads.nthreads()] + local_flux_viscous_threaded = [ntuple(_ -> similar(u_transformed, dg.basis.Nq), NDIMS) for _ in 1:Threads.nthreads()] local_flux_face_values_threaded = [similar(scalar_flux_face_values[:, 1]) for _ in 1:Threads.nthreads()] return (; u_transformed, gradients, flux_viscous, - weak_differentiation_matrices, strong_differentiation_matrices, lift_matrix, + weak_differentiation_matrices, strong_differentiation_matrices, + gradient_lift_matrix, projection_face_interpolation_matrix, divergence_lift_matrix, u_face_values, gradients_face_values, scalar_flux_face_values, local_u_values_threaded, local_flux_viscous_threaded, local_flux_face_values_threaded) end @@ -50,7 +55,7 @@ end function calc_gradient_surface_integral!(gradients, u, scalar_flux_face_values, mesh, equations::AbstractEquationsParabolic, dg::DGMulti, cache, cache_parabolic) - @unpack lift_matrix, local_flux_face_values_threaded = cache_parabolic + (; gradient_lift_matrix, local_flux_face_values_threaded) = cache_parabolic @threaded for e in eachelement(mesh, dg) local_flux_values = local_flux_face_values_threaded[Threads.threadid()] for dim in eachdim(mesh) @@ -58,7 +63,7 @@ function calc_gradient_surface_integral!(gradients, u, scalar_flux_face_values, # compute flux * (nx, ny, nz) local_flux_values[i] = scalar_flux_face_values[i, e] * mesh.md.nxyzJ[dim][i, e] end - apply_to_each_field(mul_by_accum!(lift_matrix), view(gradients[dim], :, e), local_flux_values) + apply_to_each_field(mul_by_accum!(gradient_lift_matrix), view(gradients[dim], :, e), local_flux_values) end end end @@ -186,7 +191,7 @@ function calc_viscous_fluxes!(flux_viscous, u, gradients, mesh::DGMultiMesh, reset_du!(flux_viscous[dim], dg) end - @unpack local_flux_viscous_threaded, local_u_values_threaded = cache_parabolic + (; local_u_values_threaded) = cache_parabolic @threaded for e in eachelement(mesh, dg) @@ -196,26 +201,16 @@ 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)) - local_flux_viscous = local_flux_viscous_threaded[Threads.threadid()] - for dim in eachdim(mesh) - fill!(local_flux_viscous[dim], zero(eltype(local_flux_viscous[dim]))) - apply_to_each_field(mul_by!(dg.basis.Vq), local_flux_viscous[dim], view(gradients[dim], :, e)) - end - # compute viscous flux at quad points for i in eachindex(local_u_values) u_i = local_u_values[i] - gradients_i = getindex.(local_flux_viscous, i) + gradients_i = getindex.(gradients, i, e) for dim in eachdim(mesh) flux_viscous_i = flux(u_i, gradients_i, dim, equations) - setindex!(local_flux_viscous[dim], flux_viscous_i, i) + setindex!(flux_viscous[dim], flux_viscous_i, i, e) end end - # project back to the DG approximation space - for dim in eachdim(mesh) - apply_to_each_field(mul_by!(dg.basis.Pq), view(flux_viscous[dim], :, e), local_flux_viscous[dim]) - end end end @@ -292,13 +287,12 @@ function calc_divergence!(du, u::StructArray, t, flux_viscous, mesh::DGMultiMesh calc_boundary_flux!(scalar_flux_face_values, cache_parabolic.u_face_values, t, 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) # surface contributions - apply_to_each_field(mul_by_accum!(dg.basis.LIFT), du, scalar_flux_face_values) + apply_to_each_field(mul_by_accum!(cache_parabolic.divergence_lift_matrix), du, scalar_flux_face_values) # 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 From 6d3cd0dc76606d8ca09e55f868e0c7c2862251b0 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 24 Apr 2023 16:47:45 -0500 Subject: [PATCH 10/36] remove some @unpacks --- src/solvers/dgmulti/dg_parabolic.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index 6d89455caef..3b46897954f 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -156,7 +156,7 @@ function calc_single_boundary_flux!(flux_face_values, u_face_values, t, md = mesh.md num_pts_per_face = rd.Nfq ÷ rd.Nfaces - @unpack xyzf, nxyz = md + (; xyzf, nxyz) = md for f in mesh.boundary_faces[boundary_key] for i in Base.OneTo(num_pts_per_face) @@ -225,8 +225,8 @@ function calc_viscous_penalty!(scalar_flux_face_values, u_face_values, t, bounda mesh, equations::AbstractEquationsParabolic, dg::DGMulti, parabolic_scheme, cache, cache_parabolic) # compute fluxes at interfaces - @unpack scalar_flux_face_values = cache_parabolic - @unpack mapM, mapP = mesh.md + (; scalar_flux_face_values) = cache_parabolic + (; mapM, mapP) = mesh.md @threaded for face_node_index in each_face_node_global(mesh, dg) idM, idP = mapM[face_node_index], mapP[face_node_index] uM, uP = u_face_values[idM], u_face_values[idP] @@ -313,7 +313,7 @@ function rhs_parabolic!(du, u, t, mesh::DGMultiMesh, equations_parabolic::Abstra reset_du!(du, dg) - @unpack u_transformed, gradients, flux_viscous = cache_parabolic + (; u_transformed, gradients, flux_viscous) = cache_parabolic transform_variables!(u_transformed, u, mesh, equations_parabolic, dg, parabolic_scheme, cache, cache_parabolic) From 2f134732c5a1a1d28ea66d2a69f34e0792574435 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 24 Apr 2023 17:41:45 -0500 Subject: [PATCH 11/36] update parabolic tests --- test/test_parabolic_2d.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index ac4adbfd69b..27b4c140f6d 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -53,9 +53,9 @@ isdir(outdir) && rm(outdir, recursive=true) # pass in `boundary_condition_periodic` to skip boundary flux/integral evaluation Trixi.calc_gradient!(gradients, ode.u0, t, mesh, equations_parabolic, boundary_condition_periodic, dg, cache, cache_parabolic) - @unpack x, y = mesh.md - @test getindex.(gradients[1], 1) ≈ 2 * x .* y - @test getindex.(gradients[2], 1) ≈ x.^2 + @unpack xq, yq = mesh.md + @test getindex.(gradients[1], 1) ≈ 2 * xq .* yq + @test getindex.(gradients[2], 1) ≈ xq.^2 u_flux = similar.(gradients) Trixi.calc_viscous_fluxes!(u_flux, ode.u0, gradients, mesh, equations_parabolic, From 09554f19c209da705e01000e13c223fa30d91a23 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 24 Apr 2023 21:56:36 -0500 Subject: [PATCH 12/36] fix test again --- test/test_parabolic_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 27b4c140f6d..0b73208c2eb 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -66,7 +66,7 @@ isdir(outdir) && rm(outdir, recursive=true) du = similar(ode.u0) Trixi.calc_divergence!(du, ode.u0, t, u_flux, mesh, equations_parabolic, boundary_condition_periodic, dg, semi.solver_parabolic, cache, cache_parabolic) - @test getindex.(du, 1) ≈ 2 * y + @test getindex.(du, 1) ≈ 2 * yq end @trixi_testset "DGMulti: elixir_advection_diffusion.jl" begin From 99706c71acfe878f438e0d02e918c2d0dd9b8f16 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 25 Apr 2023 06:45:27 -0500 Subject: [PATCH 13/36] fix test again --- test/test_parabolic_2d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 0b73208c2eb..44b6d107d74 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -53,7 +53,7 @@ isdir(outdir) && rm(outdir, recursive=true) # pass in `boundary_condition_periodic` to skip boundary flux/integral evaluation Trixi.calc_gradient!(gradients, ode.u0, t, mesh, equations_parabolic, boundary_condition_periodic, dg, cache, cache_parabolic) - @unpack xq, yq = mesh.md + @unpack x, y, xq, yq = mesh.md @test getindex.(gradients[1], 1) ≈ 2 * xq .* yq @test getindex.(gradients[2], 1) ≈ xq.^2 @@ -66,7 +66,7 @@ isdir(outdir) && rm(outdir, recursive=true) du = similar(ode.u0) Trixi.calc_divergence!(du, ode.u0, t, u_flux, mesh, equations_parabolic, boundary_condition_periodic, dg, semi.solver_parabolic, cache, cache_parabolic) - @test getindex.(du, 1) ≈ 2 * yq + @test getindex.(du, 1) ≈ 2 * y end @trixi_testset "DGMulti: elixir_advection_diffusion.jl" begin From e14f3ce5965ff9f90507d329d71d84d6848f37dd Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 25 Apr 2023 11:54:54 -0500 Subject: [PATCH 14/36] define md earlier --- src/solvers/dgmulti/dg_parabolic.jl | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index 3b46897954f..e48b1101429 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -19,8 +19,13 @@ function create_cache_parabolic(mesh::DGMultiMesh, divergence_lift_matrix = dg.basis.LIFT projection_face_interpolation_matrix = dg.basis.Vf * dg.basis.Pq - # u_transformed stores "transformed" variables for computing the gradient + # evaluate geometric terms at quadrature points in case the mesh is curved (; md) = mesh + J = dg.basis.Vq * md.J + invJ = inv.(J) + dxidxhatj = map(x -> dg.basis.Vq * x, md.rstxyzJ) + + # u_transformed stores "transformed" variables for computing the gradient u_transformed = allocate_nested_array(uEltype, nvars, size(md.x), dg) gradients = SVector{NDIMS}(ntuple(_ -> similar(u_transformed, (dg.basis.Nq, mesh.md.num_elements)), NDIMS)) flux_viscous = similar.(gradients) From e1f1d01a081f0dd640d25ec4cfc00f205ab141ea Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 25 Apr 2023 11:55:23 -0500 Subject: [PATCH 15/36] use SVector for consistency/efficiency --- src/solvers/dgmulti/dg_parabolic.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index e48b1101429..e3bf065d8dc 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -35,12 +35,13 @@ function create_cache_parabolic(mesh::DGMultiMesh, gradients_face_values = ntuple(_ -> similar(u_face_values), NDIMS) local_u_values_threaded = [similar(u_transformed, dg.basis.Nq) for _ in 1:Threads.nthreads()] - local_flux_viscous_threaded = [ntuple(_ -> similar(u_transformed, dg.basis.Nq), NDIMS) for _ in 1:Threads.nthreads()] + local_flux_viscous_threaded = [SVector{NDIMS}(ntuple(_ -> similar(u_transformed, dg.basis.Nq), NDIMS)) for _ in 1:Threads.nthreads()] local_flux_face_values_threaded = [similar(scalar_flux_face_values[:, 1]) for _ in 1:Threads.nthreads()] return (; u_transformed, gradients, flux_viscous, 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_flux_face_values_threaded) end From 0c2e7ac61fa291701bb505553e7980a7044f47d3 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 25 Apr 2023 11:55:48 -0500 Subject: [PATCH 16/36] add curved gradient volume integral --- src/solvers/dgmulti/dg_parabolic.jl | 24 +++++++++++++++++++++++- 1 file changed, 23 insertions(+), 1 deletion(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index e3bf065d8dc..4be5ca9e167 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -57,7 +57,6 @@ function transform_variables!(u_transformed, u, mesh, equations_parabolic::Abstr end # TODO: reuse entropy projection computations for DGMultiFluxDiff{<:Polynomial} (including `GaussSBP` solvers) - function calc_gradient_surface_integral!(gradients, u, scalar_flux_face_values, mesh, equations::AbstractEquationsParabolic, dg::DGMulti, cache, cache_parabolic) @@ -90,6 +89,29 @@ function calc_gradient_volume_integral!(gradients, u, mesh::DGMultiMesh, end end +function calc_gradient_volume_integral!(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 + + # 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()] + for i in eachdim(mesh) + apply_to_each_field(mul_by!(strong_differentiation_matrices[i]), + local_reference_gradients[i], view(u, :, e)) + end + + # rotate to physical frame on each element + for i in eachdim(mesh), j in eachdim(mesh) + @. gradients[i][:, e] = gradients[i][:, e] + dxidxhatj[i, j][:, e] * local_reference_gradients[j] + end + end +end + function calc_gradient!(gradients, u::StructArray, t, mesh::DGMultiMesh, equations::AbstractEquationsParabolic, boundary_conditions, dg::DGMulti, cache, cache_parabolic) From dc78c5a5308c4f57e735cbeae19c86aaabb0ee75 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 25 Apr 2023 11:56:03 -0500 Subject: [PATCH 17/36] add curved "invert_jacobian_gradient!" --- src/solvers/dgmulti/dg_parabolic.jl | 26 +++++++++++++++++++++++--- 1 file changed, 23 insertions(+), 3 deletions(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index 4be5ca9e167..c4ee4846d74 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -142,9 +142,16 @@ function calc_gradient!(gradients, u::StructArray, t, mesh::DGMultiMesh, calc_gradient_surface_integral!(gradients, u, scalar_flux_face_values, mesh, equations, dg, cache, cache_parabolic) - for dim in eachdim(mesh) - @threaded for e in eachelement(mesh, dg) - invJ = cache.invJ[1, e] + invert_jacobian_gradient!(gradients, mesh, equations, dg, cache, cache_parabolic) + +end + +# affine mesh - constant Jacobian version +function invert_jacobian_gradient!(gradients, mesh::DGMultiMesh, equations, dg::DGMulti, + cache, cache_parabolic) + @threaded for e in eachelement(mesh, dg) + invJ = cache_parabolic.invJ[1, e] + for dim in eachdim(mesh) for i in axes(gradients[dim], 1) gradients[dim][i, e] = gradients[dim][i, e] * invJ end @@ -152,6 +159,19 @@ function calc_gradient!(gradients, u::StructArray, t, mesh::DGMultiMesh, end end +# non-affine mesh - variable Jacobian version +function invert_jacobian_gradient!(gradients, mesh::DGMultiMesh{NDIMS, <:NonAffine}, equations, + dg::DGMulti, cache, cache_parabolic) where {NDIMS} + (; invJ) = cache_parabolic + @threaded for e in eachelement(mesh, dg) + for dim in eachdim(mesh) + for i in axes(gradients[dim], 1) + gradients[dim][i, e] = gradients[dim][i, e] * invJ[i, e] + end + end + end +end + # do nothing for periodic domains function calc_boundary_flux!(flux, u, t, operator_type, ::BoundaryConditionPeriodic, mesh, equations::AbstractEquationsParabolic, dg::DGMulti, From 387b6e45e78ab4388a0a0ec8693ff77e968c56fd Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 25 Apr 2023 11:56:17 -0500 Subject: [PATCH 18/36] add curved divergence volume integral --- src/solvers/dgmulti/dg_parabolic.jl | 25 ++++++++++++++++++++++++- 1 file changed, 24 insertions(+), 1 deletion(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index c4ee4846d74..81e9b26fff3 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -293,7 +293,30 @@ function calc_divergence_volume_integral!(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_viscous[i], :, e)) + end + end +end + +function calc_divergence_volume_integral!(du, u, flux_viscous, mesh::DGMultiMesh{NDIMS, <:NonAffine}, + equations::AbstractEquationsParabolic, + dg::DGMulti, cache, cache_parabolic) where {NDIMS} + (; weak_differentiation_matrices, dxidxhatj, local_flux_viscous_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] + for i in eachdim(mesh) + # rotate flux to reference coordinates + fill!(local_viscous_flux, zero(eltype(local_viscous_flux))) + for j in eachdim(mesh) + @. local_viscous_flux = local_viscous_flux + dxidxhatj[j, i][:, e] * flux_viscous[j][:, e] + end + + # differentiate with respect to reference coordinates + apply_to_each_field(mul_by_accum!(weak_differentiation_matrices[i]), + view(du, :, e), local_viscous_flux) end end end From dd216dae58a815ec81f3dae65e6c9349468ffe83 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 25 Apr 2023 11:56:23 -0500 Subject: [PATCH 19/36] add curved convergence elixir --- .../elixir_navierstokes_convergence_curved.jl | 212 ++++++++++++++++++ 1 file changed, 212 insertions(+) create mode 100644 examples/dgmulti_2d/elixir_navierstokes_convergence_curved.jl diff --git a/examples/dgmulti_2d/elixir_navierstokes_convergence_curved.jl b/examples/dgmulti_2d/elixir_navierstokes_convergence_curved.jl new file mode 100644 index 00000000000..854c2af773d --- /dev/null +++ b/examples/dgmulti_2d/elixir_navierstokes_convergence_curved.jl @@ -0,0 +1,212 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the ideal compressible Navier-Stokes equations + +prandtl_number() = 0.72 +mu() = 0.01 + +equations = CompressibleEulerEquations2D(1.4) +# Note: If you change the Navier-Stokes parameters here, also change them in the initial condition +# I really do not like this structure but it should work for now +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu=mu(), Prandtl=prandtl_number(), + gradient_variables=GradientVariablesPrimitive()) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial(), + surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs), + volume_integral = VolumeIntegralWeakForm()) + +top_bottom(x, tol=50*eps()) = abs(abs(x[2]) - 1) < tol +is_on_boundary = Dict(:top_bottom => top_bottom) + +function mapping(xi, eta) + x = xi + 0.1 * sin(pi * xi) * sin(pi * eta) + y = eta + 0.1 * sin(pi * xi) * sin(pi * eta) + return SVector(x, y) +end +cells_per_dimension = (16, 16) +mesh = DGMultiMesh(dg, cells_per_dimension, mapping; periodicity=(true, false), is_on_boundary) + +# Note: the initial condition cannot be specialized to `CompressibleNavierStokesDiffusion2D` +# since it is called by both the parabolic solver (which passes in `CompressibleNavierStokesDiffusion2D`) +# and by the initial condition (which passes in `CompressibleEulerEquations2D`). +# This convergence test setup was originally derived by Andrew Winters (@andrewwinters5000) +function initial_condition_navier_stokes_convergence_test(x, t, equations) + # Amplitude and shift + A = 0.5 + c = 2.0 + + # convenience values for trig. functions + pi_x = pi * x[1] + pi_y = pi * x[2] + pi_t = pi * t + + rho = c + A * sin(pi_x) * cos(pi_y) * cos(pi_t) + v1 = sin(pi_x) * log(x[2] + 2.0) * (1.0 - exp(-A * (x[2] - 1.0)) ) * cos(pi_t) + v2 = v1 + p = rho^2 + + return prim2cons(SVector(rho, v1, v2, p), equations) +end + +@inline function source_terms_navier_stokes_convergence_test(u, x, t, equations) + y = x[2] + + # TODO: parabolic + # we currently need to hardcode these parameters until we fix the "combined equation" issue + # see also https://github.com/trixi-framework/Trixi.jl/pull/1160 + inv_gamma_minus_one = inv(equations.gamma - 1) + Pr = prandtl_number() + mu_ = mu() + + # Same settings as in `initial_condition` + # Amplitude and shift + A = 0.5 + c = 2.0 + + # convenience values for trig. functions + pi_x = pi * x[1] + pi_y = pi * x[2] + pi_t = pi * t + + # compute the manufactured solution and all necessary derivatives + rho = c + A * sin(pi_x) * cos(pi_y) * cos(pi_t) + rho_t = -pi * A * sin(pi_x) * cos(pi_y) * sin(pi_t) + rho_x = pi * A * cos(pi_x) * cos(pi_y) * cos(pi_t) + rho_y = -pi * A * sin(pi_x) * sin(pi_y) * cos(pi_t) + rho_xx = -pi * pi * A * sin(pi_x) * cos(pi_y) * cos(pi_t) + rho_yy = -pi * pi * A * sin(pi_x) * cos(pi_y) * cos(pi_t) + + v1 = sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t) + v1_t = -pi * sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * sin(pi_t) + v1_x = pi * cos(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t) + v1_y = sin(pi_x) * (A * log(y + 2.0) * exp(-A * (y - 1.0)) + (1.0 - exp(-A * (y - 1.0))) / (y + 2.0)) * cos(pi_t) + v1_xx = -pi * pi * sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t) + v1_xy = pi * cos(pi_x) * (A * log(y + 2.0) * exp(-A * (y - 1.0)) + (1.0 - exp(-A * (y - 1.0))) / (y + 2.0)) * cos(pi_t) + v1_yy = (sin(pi_x) * ( 2.0 * A * exp(-A * (y - 1.0)) / (y + 2.0) + - A * A * log(y + 2.0) * exp(-A * (y - 1.0)) + - (1.0 - exp(-A * (y - 1.0))) / ((y + 2.0) * (y + 2.0))) * cos(pi_t)) + v2 = v1 + v2_t = v1_t + v2_x = v1_x + v2_y = v1_y + v2_xx = v1_xx + v2_xy = v1_xy + v2_yy = v1_yy + + p = rho * rho + p_t = 2.0 * rho * rho_t + p_x = 2.0 * rho * rho_x + p_y = 2.0 * rho * rho_y + p_xx = 2.0 * rho * rho_xx + 2.0 * rho_x * rho_x + p_yy = 2.0 * rho * rho_yy + 2.0 * rho_y * rho_y + + # Note this simplifies slightly because the ansatz assumes that v1 = v2 + E = p * inv_gamma_minus_one + 0.5 * rho * (v1^2 + v2^2) + E_t = p_t * inv_gamma_minus_one + rho_t * v1^2 + 2.0 * rho * v1 * v1_t + E_x = p_x * inv_gamma_minus_one + rho_x * v1^2 + 2.0 * rho * v1 * v1_x + E_y = p_y * inv_gamma_minus_one + rho_y * v1^2 + 2.0 * rho * v1 * v1_y + + # Some convenience constants + T_const = equations.gamma * inv_gamma_minus_one / Pr + inv_rho_cubed = 1.0 / (rho^3) + + # compute the source terms + # density equation + du1 = rho_t + rho_x * v1 + rho * v1_x + rho_y * v2 + rho * v2_y + + # x-momentum equation + du2 = ( rho_t * v1 + rho * v1_t + p_x + rho_x * v1^2 + + 2.0 * rho * v1 * v1_x + + rho_y * v1 * v2 + + rho * v1_y * v2 + + rho * v1 * v2_y + # stress tensor from x-direction + - 4.0 / 3.0 * v1_xx * mu_ + + 2.0 / 3.0 * v2_xy * mu_ + - v1_yy * mu_ + - v2_xy * mu_ ) + # y-momentum equation + du3 = ( rho_t * v2 + rho * v2_t + p_y + rho_x * v1 * v2 + + rho * v1_x * v2 + + rho * v1 * v2_x + + rho_y * v2^2 + + 2.0 * rho * v2 * v2_y + # stress tensor from y-direction + - v1_xy * mu_ + - v2_xx * mu_ + - 4.0 / 3.0 * v2_yy * mu_ + + 2.0 / 3.0 * v1_xy * mu_ ) + # total energy equation + du4 = ( E_t + v1_x * (E + p) + v1 * (E_x + p_x) + + v2_y * (E + p) + v2 * (E_y + p_y) + # stress tensor and temperature gradient terms from x-direction + - 4.0 / 3.0 * v1_xx * v1 * mu_ + + 2.0 / 3.0 * v2_xy * v1 * mu_ + - 4.0 / 3.0 * v1_x * v1_x * mu_ + + 2.0 / 3.0 * v2_y * v1_x * mu_ + - v1_xy * v2 * mu_ + - v2_xx * v2 * mu_ + - v1_y * v2_x * mu_ + - v2_x * v2_x * mu_ + - T_const * inv_rho_cubed * ( p_xx * rho * rho + - 2.0 * p_x * rho * rho_x + + 2.0 * p * rho_x * rho_x + - p * rho * rho_xx ) * mu_ + # stress tensor and temperature gradient terms from y-direction + - v1_yy * v1 * mu_ + - v2_xy * v1 * mu_ + - v1_y * v1_y * mu_ + - v2_x * v1_y * mu_ + - 4.0 / 3.0 * v2_yy * v2 * mu_ + + 2.0 / 3.0 * v1_xy * v2 * mu_ + - 4.0 / 3.0 * v2_y * v2_y * mu_ + + 2.0 / 3.0 * v1_x * v2_y * mu_ + - T_const * inv_rho_cubed * ( p_yy * rho * rho + - 2.0 * p_y * rho * rho_y + + 2.0 * p * rho_y * rho_y + - p * rho * rho_yy ) * mu_ ) + + return SVector(du1, du2, du3, du4) +end + +initial_condition = initial_condition_navier_stokes_convergence_test + +# BC types +velocity_bc_top_bottom = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x, t, equations)[2:3]) +heat_bc_top_bottom = Adiabatic((x, t, equations) -> 0.0) +boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom, heat_bc_top_bottom) + +# define inviscid boundary conditions +boundary_conditions = (; :top_bottom => boundary_condition_slip_wall) + +# define viscous boundary conditions +boundary_conditions_parabolic = (; :top_bottom => boundary_condition_top_bottom) + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, dg; + boundary_conditions=(boundary_conditions, boundary_conditions_parabolic), + source_terms=source_terms_navier_stokes_convergence_test) + + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span `tspan` +tspan = (0.0, 0.5) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() +alive_callback = AliveCallback(alive_interval=10) +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg)) +callbacks = CallbackSet(summary_callback, alive_callback) + +############################################################################### +# run the simulation + +time_int_tol = 1e-8 +sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol, + ode_default_options()..., callback=callbacks) +summary_callback() # print the timer summary From 2db366a81092e3166597a786a98e74558ad522ef Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 25 Apr 2023 11:57:35 -0500 Subject: [PATCH 20/36] add curved test --- test/test_parabolic_2d.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 44b6d107d74..588f43e4543 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -101,6 +101,14 @@ isdir(outdir) && rm(outdir, recursive=true) ) end + @trixi_testset "DGMulti: elixir_navierstokes_convergence_curved.jl" begin + @test_trixi_include(joinpath(examples_dir(), "dgmulti_2d", "elixir_navierstokes_convergence_curved.jl"), + cells_per_dimension = (4, 4), tspan=(0.0, 0.1), + l2 = [0.004255101916146187, 0.011118488923215765, 0.011281831283462686, 0.03573656447388509], + linf = [0.015071710669706473, 0.04103132025858458, 0.03990424085750277, 0.1309401718598764], + ) + end + @trixi_testset "DGMulti: elixir_navierstokes_lid_driven_cavity.jl" begin @test_trixi_include(joinpath(examples_dir(), "dgmulti_2d", "elixir_navierstokes_lid_driven_cavity.jl"), cells_per_dimension = (4, 4), tspan=(0.0, 0.5), From bfbbb7357f9ea2d9b88994cce0673c0d4aeabc71 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 25 Apr 2023 12:00:33 -0500 Subject: [PATCH 21/36] add 3D test --- .../elixir_navierstokes_convergence_curved.jl | 261 ++++++++++++++++++ test/test_parabolic_3d.jl | 8 + 2 files changed, 269 insertions(+) create mode 100644 examples/dgmulti_3d/elixir_navierstokes_convergence_curved.jl diff --git a/examples/dgmulti_3d/elixir_navierstokes_convergence_curved.jl b/examples/dgmulti_3d/elixir_navierstokes_convergence_curved.jl new file mode 100644 index 00000000000..b8971a645e2 --- /dev/null +++ b/examples/dgmulti_3d/elixir_navierstokes_convergence_curved.jl @@ -0,0 +1,261 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the ideal compressible Navier-Stokes equations + +prandtl_number() = 0.72 +mu() = 0.01 + +equations = CompressibleEulerEquations3D(1.4) +equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu=mu(), Prandtl=prandtl_number(), + gradient_variables=GradientVariablesPrimitive()) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +dg = DGMulti(polydeg = 3, element_type = Hex(), approximation_type = Polynomial(), + surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs), + volume_integral = VolumeIntegralWeakForm()) + +top_bottom(x, tol=50*eps()) = abs(abs(x[2]) - 1) < tol +is_on_boundary = Dict(:top_bottom => top_bottom) + +function mapping(xi, eta, zeta) + x = xi + 0.1 * sin(pi * xi) * sin(pi * eta) + y = eta + 0.1 * sin(pi * xi) * sin(pi * eta) + z = zeta + 0.1 * sin(pi * xi) * sin(pi * eta) + return SVector(x, y, z) +end +cells_per_dimension = (8, 8, 8) +mesh = DGMultiMesh(dg, cells_per_dimension, mapping; periodicity=(true, false, true), is_on_boundary) + +# Note: the initial condition cannot be specialized to `CompressibleNavierStokesDiffusion3D` +# since it is called by both the parabolic solver (which passes in `CompressibleNavierStokesDiffusion3D`) +# and by the initial condition (which passes in `CompressibleEulerEquations3D`). +# This convergence test setup was originally derived by Andrew Winters (@andrewwinters5000) +function initial_condition_navier_stokes_convergence_test(x, t, equations) + # Constants. OBS! Must match those in `source_terms_navier_stokes_convergence_test` + c = 2.0 + A1 = 0.5 + A2 = 1.0 + A3 = 0.5 + + # Convenience values for trig. functions + pi_x = pi * x[1] + pi_y = pi * x[2] + pi_z = pi * x[3] + pi_t = pi * t + + rho = c + A1 * sin(pi_x) * cos(pi_y) * sin(pi_z) * cos(pi_t) + v1 = A2 * sin(pi_x) * log(x[2] + 2.0) * (1.0 - exp(-A3 * (x[2] - 1.0))) * sin(pi_z) * cos(pi_t) + v2 = v1 + v3 = v1 + p = rho^2 + + return prim2cons(SVector(rho, v1, v2, v3, p), equations) +end + +@inline function source_terms_navier_stokes_convergence_test(u, x, t, equations) + # TODO: parabolic + # we currently need to hardcode these parameters until we fix the "combined equation" issue + # see also https://github.com/trixi-framework/Trixi.jl/pull/1160 + inv_gamma_minus_one = inv(equations.gamma - 1) + Pr = prandtl_number() + mu_ = mu() + + # Constants. OBS! Must match those in `initial_condition_navier_stokes_convergence_test` + c = 2.0 + A1 = 0.5 + A2 = 1.0 + A3 = 0.5 + + # Convenience values for trig. functions + pi_x = pi * x[1] + pi_y = pi * x[2] + pi_z = pi * x[3] + pi_t = pi * t + + # Define auxiliary functions for the strange function of the y variable + # to make expressions easier to read + g = log(x[2] + 2.0) * (1.0 - exp(-A3 * (x[2] - 1.0))) + g_y = ( A3 * log(x[2] + 2.0) * exp(-A3 * (x[2] - 1.0)) + + (1.0 - exp(-A3 * (x[2] - 1.0))) / (x[2] + 2.0) ) + g_yy = ( 2.0 * A3 * exp(-A3 * (x[2] - 1.0)) / (x[2] + 2.0) + - (1.0 - exp(-A3 * (x[2] - 1.0))) / ((x[2] + 2.0)^2) + - A3^2 * log(x[2] + 2.0) * exp(-A3 * (x[2] - 1.0)) ) + + # Density and its derivatives + rho = c + A1 * sin(pi_x) * cos(pi_y) * sin(pi_z) * cos(pi_t) + rho_t = -pi * A1 * sin(pi_x) * cos(pi_y) * sin(pi_z) * sin(pi_t) + rho_x = pi * A1 * cos(pi_x) * cos(pi_y) * sin(pi_z) * cos(pi_t) + rho_y = -pi * A1 * sin(pi_x) * sin(pi_y) * sin(pi_z) * cos(pi_t) + rho_z = pi * A1 * sin(pi_x) * cos(pi_y) * cos(pi_z) * cos(pi_t) + rho_xx = -pi^2 * (rho - c) + rho_yy = -pi^2 * (rho - c) + rho_zz = -pi^2 * (rho - c) + + # Velocities and their derivatives + # v1 terms + v1 = A2 * sin(pi_x) * g * sin(pi_z) * cos(pi_t) + v1_t = -pi * A2 * sin(pi_x) * g * sin(pi_z) * sin(pi_t) + v1_x = pi * A2 * cos(pi_x) * g * sin(pi_z) * cos(pi_t) + v1_y = A2 * sin(pi_x) * g_y * sin(pi_z) * cos(pi_t) + v1_z = pi * A2 * sin(pi_x) * g * cos(pi_z) * cos(pi_t) + v1_xx = -pi^2 * v1 + v1_yy = A2 * sin(pi_x) * g_yy * sin(pi_z) * cos(pi_t) + v1_zz = -pi^2 * v1 + v1_xy = pi * A2 * cos(pi_x) * g_y * sin(pi_z) * cos(pi_t) + v1_xz = pi^2 * A2 * cos(pi_x) * g * cos(pi_z) * cos(pi_t) + v1_yz = pi * A2 * sin(pi_x) * g_y * cos(pi_z) * cos(pi_t) + # v2 terms (simplifies from ansatz) + v2 = v1 + v2_t = v1_t + v2_x = v1_x + v2_y = v1_y + v2_z = v1_z + v2_xx = v1_xx + v2_yy = v1_yy + v2_zz = v1_zz + v2_xy = v1_xy + v2_yz = v1_yz + # v3 terms (simplifies from ansatz) + v3 = v1 + v3_t = v1_t + v3_x = v1_x + v3_y = v1_y + v3_z = v1_z + v3_xx = v1_xx + v3_yy = v1_yy + v3_zz = v1_zz + v3_xz = v1_xz + v3_yz = v1_yz + + # Pressure and its derivatives + p = rho^2 + p_t = 2.0 * rho * rho_t + p_x = 2.0 * rho * rho_x + p_y = 2.0 * rho * rho_y + p_z = 2.0 * rho * rho_z + + # Total energy and its derivatives; simiplifies from ansatz that v2 = v1 and v3 = v1 + E = p * inv_gamma_minus_one + 1.5 * rho * v1^2 + E_t = p_t * inv_gamma_minus_one + 1.5 * rho_t * v1^2 + 3.0 * rho * v1 * v1_t + E_x = p_x * inv_gamma_minus_one + 1.5 * rho_x * v1^2 + 3.0 * rho * v1 * v1_x + E_y = p_y * inv_gamma_minus_one + 1.5 * rho_y * v1^2 + 3.0 * rho * v1 * v1_y + E_z = p_z * inv_gamma_minus_one + 1.5 * rho_z * v1^2 + 3.0 * rho * v1 * v1_z + + # Divergence of Fick's law ∇⋅∇q = kappa ∇⋅∇T; simplifies because p = rho², so T = p/rho = rho + kappa = equations.gamma * inv_gamma_minus_one / Pr + q_xx = kappa * rho_xx # kappa T_xx + q_yy = kappa * rho_yy # kappa T_yy + q_zz = kappa * rho_zz # kappa T_zz + + # Stress tensor and its derivatives (exploit symmetry) + tau11 = 4.0 / 3.0 * v1_x - 2.0 / 3.0 * (v2_y + v3_z) + tau12 = v1_y + v2_x + tau13 = v1_z + v3_x + tau22 = 4.0 / 3.0 * v2_y - 2.0 / 3.0 * (v1_x + v3_z) + tau23 = v2_z + v3_y + tau33 = 4.0 / 3.0 * v3_z - 2.0 / 3.0 * (v1_x + v2_y) + + tau11_x = 4.0 / 3.0 * v1_xx - 2.0 / 3.0 * (v2_xy + v3_xz) + tau12_x = v1_xy + v2_xx + tau13_x = v1_xz + v3_xx + + tau12_y = v1_yy + v2_xy + tau22_y = 4.0 / 3.0 * v2_yy - 2.0 / 3.0 * (v1_xy + v3_yz) + tau23_y = v2_yz + v3_yy + + tau13_z = v1_zz + v3_xz + tau23_z = v2_zz + v3_yz + tau33_z = 4.0 / 3.0 * v3_zz - 2.0 / 3.0 * (v1_xz + v2_yz) + + # Compute the source terms + # Density equation + du1 = ( rho_t + rho_x * v1 + rho * v1_x + + rho_y * v2 + rho * v2_y + + rho_z * v3 + rho * v3_z ) + # x-momentum equation + du2 = ( rho_t * v1 + rho * v1_t + p_x + rho_x * v1^2 + + 2.0 * rho * v1 * v1_x + + rho_y * v1 * v2 + + rho * v1_y * v2 + + rho * v1 * v2_y + + rho_z * v1 * v3 + + rho * v1_z * v3 + + rho * v1 * v3_z + - mu_ * (tau11_x + tau12_y + tau13_z) ) + # y-momentum equation + du3 = ( rho_t * v2 + rho * v2_t + p_y + rho_x * v1 * v2 + + rho * v1_x * v2 + + rho * v1 * v2_x + + rho_y * v2^2 + + 2.0 * rho * v2 * v2_y + + rho_z * v2 * v3 + + rho * v2_z * v3 + + rho * v2 * v3_z + - mu_ * (tau12_x + tau22_y + tau23_z) ) + # z-momentum equation + du4 = ( rho_t * v3 + rho * v3_t + p_z + rho_x * v1 * v3 + + rho * v1_x * v3 + + rho * v1 * v3_x + + rho_y * v2 * v3 + + rho * v2_y * v3 + + rho * v2 * v3_y + + rho_z * v3^2 + + 2.0 * rho * v3 * v3_z + - mu_ * (tau13_x + tau23_y + tau33_z) ) + # Total energy equation + du5 = ( E_t + v1_x * (E + p) + v1 * (E_x + p_x) + + v2_y * (E + p) + v2 * (E_y + p_y) + + v3_z * (E + p) + v3 * (E_z + p_z) + # stress tensor and temperature gradient from x-direction + - mu_ * ( q_xx + v1_x * tau11 + v2_x * tau12 + v3_x * tau13 + + v1 * tau11_x + v2 * tau12_x + v3 * tau13_x) + # stress tensor and temperature gradient terms from y-direction + - mu_ * ( q_yy + v1_y * tau12 + v2_y * tau22 + v3_y * tau23 + + v1 * tau12_y + v2 * tau22_y + v3 * tau23_y) + # stress tensor and temperature gradient terms from z-direction + - mu_ * ( q_zz + v1_z * tau13 + v2_z * tau23 + v3_z * tau33 + + v1 * tau13_z + v2 * tau23_z + v3 * tau33_z) ) + + return SVector(du1, du2, du3, du4, du5) +end + +initial_condition = initial_condition_navier_stokes_convergence_test + +# BC types +velocity_bc_top_bottom = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x, t, equations)[2:4]) +heat_bc_top_bottom = Adiabatic((x, t, equations) -> 0.0) +boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom, heat_bc_top_bottom) + +# define inviscid boundary conditions +boundary_conditions = (; :top_bottom => boundary_condition_slip_wall) + +# define viscous boundary conditions +boundary_conditions_parabolic = (; :top_bottom => boundary_condition_top_bottom) + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, dg; + boundary_conditions=(boundary_conditions, boundary_conditions_parabolic), + source_terms=source_terms_navier_stokes_convergence_test) + + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span `tspan` +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() +alive_callback = AliveCallback(alive_interval=10) +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg)) +callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback) + +############################################################################### +# run the simulation + +time_int_tol = 1e-8 +sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol, + ode_default_options()..., callback=callbacks) +summary_callback() # print the timer summary diff --git a/test/test_parabolic_3d.jl b/test/test_parabolic_3d.jl index 8799e5ccdae..1ae5eed44ae 100644 --- a/test/test_parabolic_3d.jl +++ b/test/test_parabolic_3d.jl @@ -19,6 +19,14 @@ isdir(outdir) && rm(outdir, recursive=true) ) end + @trixi_testset "DGMulti: elixir_navierstokes_convergence_curved.jl" begin + @test_trixi_include(joinpath(examples_dir(), "dgmulti_3d", "elixir_navierstokes_convergence_curved.jl"), + cells_per_dimension = (4, 4, 4), tspan=(0.0, 0.1), + l2 = [0.0014027227251207474, 0.0021322235533273513, 0.0027873741447455194, 0.0024587473070627423, 0.00997836818019202], + linf = [0.006341750402837576, 0.010306014252246865, 0.01520740250924979, 0.010968264045485565, 0.047454389831591115] + ) + end + @trixi_testset "DGMulti: elixir_navierstokes_taylor_green_vortex.jl" begin @test_trixi_include(joinpath(examples_dir(), "dgmulti_3d", "elixir_navierstokes_taylor_green_vortex.jl"), cells_per_dimension = (4, 4, 4), tspan=(0.0, 0.25), From 0ce82023eb260d041a2fe9f2cca7de42ecc499ac Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Wed, 26 Apr 2023 10:38:39 -0500 Subject: [PATCH 22/36] add invJ comment for affine elements --- src/solvers/dgmulti/dg_parabolic.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index 81e9b26fff3..25c9698a37f 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -150,7 +150,11 @@ end function invert_jacobian_gradient!(gradients, mesh::DGMultiMesh, equations, dg::DGMulti, cache, cache_parabolic) @threaded for e in eachelement(mesh, dg) + + # Here, we exploit the fact that J is constant on for affine elements, + # so we only have to access invJ once per element. invJ = cache_parabolic.invJ[1, e] + for dim in eachdim(mesh) for i in axes(gradients[dim], 1) gradients[dim][i, e] = gradients[dim][i, e] * invJ From faefe3fb48f6f923e4be11a73011455620287fd0 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Wed, 26 Apr 2023 10:38:47 -0500 Subject: [PATCH 23/36] NDIMS -> ndims(mesh) --- src/solvers/dgmulti/dg_parabolic.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index 25c9698a37f..073954549c2 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -2,7 +2,7 @@ function create_cache_parabolic(mesh::DGMultiMesh, equations_hyperbolic::AbstractEquations, equations_parabolic::AbstractEquationsParabolic, - dg::DGMulti{NDIMS}, parabolic_scheme, RealT, uEltype) where {NDIMS} + dg::DGMulti, parabolic_scheme, RealT, uEltype) # default to taking derivatives of all hyperbolic variables # TODO: parabolic; utilize the parabolic variables in `equations_parabolic` to reduce memory usage in the parabolic cache @@ -27,15 +27,15 @@ function create_cache_parabolic(mesh::DGMultiMesh, # u_transformed stores "transformed" variables for computing the gradient u_transformed = allocate_nested_array(uEltype, nvars, size(md.x), dg) - gradients = SVector{NDIMS}(ntuple(_ -> similar(u_transformed, (dg.basis.Nq, mesh.md.num_elements)), NDIMS)) + gradients = SVector{ndims(mesh)}(ntuple(_ -> similar(u_transformed, (dg.basis.Nq, mesh.md.num_elements)), ndims(mesh))) flux_viscous = similar.(gradients) u_face_values = allocate_nested_array(uEltype, nvars, size(md.xf), dg) scalar_flux_face_values = similar(u_face_values) - gradients_face_values = ntuple(_ -> similar(u_face_values), NDIMS) + gradients_face_values = ntuple(_ -> similar(u_face_values), ndims(mesh)) local_u_values_threaded = [similar(u_transformed, dg.basis.Nq) for _ in 1:Threads.nthreads()] - local_flux_viscous_threaded = [SVector{NDIMS}(ntuple(_ -> similar(u_transformed, dg.basis.Nq), NDIMS)) for _ in 1:Threads.nthreads()] + local_flux_viscous_threaded = [SVector{ndims(mesh)}(ntuple(_ -> similar(u_transformed, dg.basis.Nq), ndims(mesh))) for _ in 1:Threads.nthreads()] local_flux_face_values_threaded = [similar(scalar_flux_face_values[:, 1]) for _ in 1:Threads.nthreads()] return (; u_transformed, gradients, flux_viscous, From 14a45a50165ec814021e65215cade15836fbb3d7 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Wed, 26 Apr 2023 10:43:27 -0500 Subject: [PATCH 24/36] comments --- src/solvers/dgmulti/dg_parabolic.jl | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index 073954549c2..d562065cfd9 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -82,7 +82,10 @@ function calc_gradient_volume_integral!(gradients, u, mesh::DGMultiMesh, # compute volume contributions to gradients @threaded for e in eachelement(mesh, dg) for i in eachdim(mesh), j in eachdim(mesh) - dxidxhatj = mesh.md.rstxyzJ[i, j][1, e] # TODO: DGMulti. Assumes mesh is affine here. + + # We assume each element is affine (e.g., constant geometric terms) here. + dxidxhatj = mesh.md.rstxyzJ[i, j][1, e] + apply_to_each_field(mul_by_accum!(strong_differentiation_matrices[j], dxidxhatj), view(gradients[i], :, e), view(u, :, e)) end @@ -132,7 +135,9 @@ function calc_gradient!(gradients, u::StructArray, t, mesh::DGMultiMesh, idM, idP = mapM[face_node_index], mapP[face_node_index] uM = u_face_values[idM] uP = u_face_values[idP] - scalar_flux_face_values[idM] = 0.5 * (uP - uM) # TODO: use strong/weak formulation for curved meshes? + # Here, we use the "strong" formulation to compute the gradient. This guarantees that the parabolic + # formulation is symmetric and stable on curved meshes with variable geometric terms. + scalar_flux_face_values[idM] = 0.5 * (uP - uM) end calc_boundary_flux!(scalar_flux_face_values, u_face_values, t, Gradient(), boundary_conditions, From 0463178cbdecbf805b3738109457d58570dc5b4e Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Wed, 26 Apr 2023 10:43:33 -0500 Subject: [PATCH 25/36] remove outdated TODO --- src/solvers/dgmulti/dg_parabolic.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index d562065cfd9..8b8f1ff2c86 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -205,7 +205,6 @@ end calc_boundary_flux!(flux, u, t, operator_type, boundary_conditions::NamedTuple{(),Tuple{}}, mesh, equations, dg::DGMulti, cache, cache_parabolic) = nothing -# TODO: DGMulti. Decide if we want to use the input `u_face_values` (currently unused) function calc_single_boundary_flux!(flux_face_values, u_face_values, t, operator_type, boundary_condition, boundary_key, mesh, equations, dg::DGMulti{NDIMS}, cache, cache_parabolic) where {NDIMS} From f31adb92db3b1eda2950b5be1cba27a4725b3824 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Wed, 26 Apr 2023 10:47:12 -0500 Subject: [PATCH 26/36] more comments --- src/solvers/dgmulti/dg_parabolic.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index 8b8f1ff2c86..d7e6509c1e6 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -229,7 +229,10 @@ function calc_single_boundary_flux!(flux_face_values, u_face_values, t, face_normal, face_coordinates, t, operator_type, equations) - # use "strong form" for the Gradient ("weak form" for Divergence) + # Here, we use the "strong form" for the Gradient (and the "weak form" for Divergence). + # `flux_face_values` should contain the boundary values for `u`, and we + # subtract off `u_face_values[fid, e]` because we are using the strong formulation to + # compute the gradient. if operator_type isa Gradient flux_face_values[fid, e] = flux_face_values[fid, e] - u_face_values[fid, e] end From f78abbe005b50e498a4d06b54c77de13193db4f4 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Wed, 26 Apr 2023 10:48:10 -0500 Subject: [PATCH 27/36] remove outdated TODOs --- src/solvers/dgmulti/dg_parabolic.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index d7e6509c1e6..97547df18a1 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -344,7 +344,6 @@ function calc_divergence!(du, u::StructArray, t, flux_viscous, mesh::DGMultiMesh (; projection_face_interpolation_matrix) = cache_parabolic flux_viscous_face_values = cache_parabolic.gradients_face_values # reuse storage for dim in eachdim(mesh) - # TODO: optimize this!! apply_to_each_field(mul_by!(projection_face_interpolation_matrix), flux_viscous_face_values[dim], flux_viscous[dim]) end @@ -360,7 +359,7 @@ function calc_divergence!(du, u::StructArray, t, flux_viscous, mesh::DGMultiMesh for dim in eachdim(mesh) fM = flux_viscous_face_values[dim][idM] fP = flux_viscous_face_values[dim][idP] - # TODO: use strong/weak formulation to ensure stability on curved meshes? + # Here, we use the "weak" formulation to compute the divergence (to ensure stability on curved meshes). flux_face_value = flux_face_value + 0.5 * (fP + fM) * nxyzJ[dim][face_node_index] end scalar_flux_face_values[idM] = flux_face_value From 22f469135f1df56d2304d9fce69ef6eb3a5ff0d4 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Wed, 26 Apr 2023 11:05:07 -0500 Subject: [PATCH 28/36] remove broadcasts to reduce allocations --- src/solvers/dgmulti/dg_parabolic.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index 97547df18a1..7e78037768e 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -110,7 +110,9 @@ function calc_gradient_volume_integral!(gradients, u, mesh::DGMultiMesh{NDIMS, < # rotate to physical frame on each element for i in eachdim(mesh), j in eachdim(mesh) - @. gradients[i][:, e] = gradients[i][:, e] + dxidxhatj[i, j][:, e] * local_reference_gradients[j] + for node in eachindex(local_reference_gradients[j]) + gradients[i][node, e] = gradients[i][node, e] + dxidxhatj[i, j][node, e] * local_reference_gradients[j][node] + end end end end @@ -322,7 +324,9 @@ function calc_divergence_volume_integral!(du, u, flux_viscous, mesh::DGMultiMesh # rotate flux to reference coordinates fill!(local_viscous_flux, zero(eltype(local_viscous_flux))) for j in eachdim(mesh) - @. local_viscous_flux = local_viscous_flux + dxidxhatj[j, i][:, e] * flux_viscous[j][:, e] + 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] + end end # differentiate with respect to reference coordinates From 5681534ca03466e6b098ea6dd831b3b2c05bdc39 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Wed, 26 Apr 2023 11:05:23 -0500 Subject: [PATCH 29/36] add timers to rhs_parabolic --- src/solvers/dgmulti/dg_parabolic.jl | 27 +++++++++++++++++---------- 1 file changed, 17 insertions(+), 10 deletions(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index 7e78037768e..d06712909ca 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -398,19 +398,26 @@ function rhs_parabolic!(du, u, t, mesh::DGMultiMesh, equations_parabolic::Abstra reset_du!(du, dg) - (; u_transformed, gradients, flux_viscous) = cache_parabolic - transform_variables!(u_transformed, u, mesh, equations_parabolic, - dg, parabolic_scheme, cache, cache_parabolic) - - calc_gradient!(gradients, u_transformed, t, mesh, equations_parabolic, - boundary_conditions, dg, cache, cache_parabolic) + @trixi_timeit timer() "transform variables" begin + (; u_transformed, gradients, flux_viscous) = cache_parabolic + transform_variables!(u_transformed, u, mesh, equations_parabolic, + dg, parabolic_scheme, cache, cache_parabolic) + end - calc_viscous_fluxes!(flux_viscous, u_transformed, gradients, - mesh, equations_parabolic, dg, cache, cache_parabolic) + @trixi_timeit timer() "calc gradient" begin + calc_gradient!(gradients, u_transformed, t, mesh, equations_parabolic, + boundary_conditions, dg, cache, cache_parabolic) + end - calc_divergence!(du, u_transformed, t, flux_viscous, mesh, equations_parabolic, - boundary_conditions, dg, parabolic_scheme, cache, cache_parabolic) + @trixi_timeit timer() "calc viscous fluxes" begin + calc_viscous_fluxes!(flux_viscous, 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, + boundary_conditions, dg, parabolic_scheme, cache, cache_parabolic) + end return nothing end From 398be9a9623502456e49254dc46e94329849e11a Mon Sep 17 00:00:00 2001 From: Jesse Chan <1156048+jlchan@users.noreply.github.com> Date: Thu, 27 Apr 2023 08:56:15 -0500 Subject: [PATCH 30/36] Update src/solvers/dgmulti/dg_parabolic.jl Co-authored-by: Hendrik Ranocha --- src/solvers/dgmulti/dg_parabolic.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index d06712909ca..7b4a876b487 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -158,7 +158,7 @@ function invert_jacobian_gradient!(gradients, mesh::DGMultiMesh, equations, dg:: cache, cache_parabolic) @threaded for e in eachelement(mesh, dg) - # Here, we exploit the fact that J is constant on for affine elements, + # Here, we exploit the fact that J is constant on affine elements, # so we only have to access invJ once per element. invJ = cache_parabolic.invJ[1, e] From 68522b2c95fb3486d2e8fe8034639ada095dffc9 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 27 Apr 2023 21:18:29 -0500 Subject: [PATCH 31/36] remove unused variable --- src/solvers/dgmulti/dg_parabolic.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index 7b4a876b487..ef5153f421f 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -132,7 +132,7 @@ function calc_gradient!(gradients, u::StructArray, t, mesh::DGMultiMesh, # compute fluxes at interfaces (; scalar_flux_face_values) = cache_parabolic - (; mapM, mapP, Jf) = mesh.md + (; mapM, mapP) = mesh.md @threaded for face_node_index in each_face_node_global(mesh, dg) idM, idP = mapM[face_node_index], mapP[face_node_index] uM = u_face_values[idM] From adc0fdec01d75c45c82f3984d6cbf728e846cf60 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 27 Apr 2023 22:05:20 -0500 Subject: [PATCH 32/36] replace extra curved elixirs with "trixi_include" calls --- .../elixir_navierstokes_convergence_curved.jl | 212 -------------- .../elixir_navierstokes_convergence_curved.jl | 261 ------------------ test/test_parabolic_2d.jl | 11 +- test/test_parabolic_3d.jl | 12 +- 4 files changed, 17 insertions(+), 479 deletions(-) delete mode 100644 examples/dgmulti_2d/elixir_navierstokes_convergence_curved.jl delete mode 100644 examples/dgmulti_3d/elixir_navierstokes_convergence_curved.jl diff --git a/examples/dgmulti_2d/elixir_navierstokes_convergence_curved.jl b/examples/dgmulti_2d/elixir_navierstokes_convergence_curved.jl deleted file mode 100644 index 854c2af773d..00000000000 --- a/examples/dgmulti_2d/elixir_navierstokes_convergence_curved.jl +++ /dev/null @@ -1,212 +0,0 @@ -using OrdinaryDiffEq -using Trixi - -############################################################################### -# semidiscretization of the ideal compressible Navier-Stokes equations - -prandtl_number() = 0.72 -mu() = 0.01 - -equations = CompressibleEulerEquations2D(1.4) -# Note: If you change the Navier-Stokes parameters here, also change them in the initial condition -# I really do not like this structure but it should work for now -equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu=mu(), Prandtl=prandtl_number(), - gradient_variables=GradientVariablesPrimitive()) - -# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux -dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial(), - surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs), - volume_integral = VolumeIntegralWeakForm()) - -top_bottom(x, tol=50*eps()) = abs(abs(x[2]) - 1) < tol -is_on_boundary = Dict(:top_bottom => top_bottom) - -function mapping(xi, eta) - x = xi + 0.1 * sin(pi * xi) * sin(pi * eta) - y = eta + 0.1 * sin(pi * xi) * sin(pi * eta) - return SVector(x, y) -end -cells_per_dimension = (16, 16) -mesh = DGMultiMesh(dg, cells_per_dimension, mapping; periodicity=(true, false), is_on_boundary) - -# Note: the initial condition cannot be specialized to `CompressibleNavierStokesDiffusion2D` -# since it is called by both the parabolic solver (which passes in `CompressibleNavierStokesDiffusion2D`) -# and by the initial condition (which passes in `CompressibleEulerEquations2D`). -# This convergence test setup was originally derived by Andrew Winters (@andrewwinters5000) -function initial_condition_navier_stokes_convergence_test(x, t, equations) - # Amplitude and shift - A = 0.5 - c = 2.0 - - # convenience values for trig. functions - pi_x = pi * x[1] - pi_y = pi * x[2] - pi_t = pi * t - - rho = c + A * sin(pi_x) * cos(pi_y) * cos(pi_t) - v1 = sin(pi_x) * log(x[2] + 2.0) * (1.0 - exp(-A * (x[2] - 1.0)) ) * cos(pi_t) - v2 = v1 - p = rho^2 - - return prim2cons(SVector(rho, v1, v2, p), equations) -end - -@inline function source_terms_navier_stokes_convergence_test(u, x, t, equations) - y = x[2] - - # TODO: parabolic - # we currently need to hardcode these parameters until we fix the "combined equation" issue - # see also https://github.com/trixi-framework/Trixi.jl/pull/1160 - inv_gamma_minus_one = inv(equations.gamma - 1) - Pr = prandtl_number() - mu_ = mu() - - # Same settings as in `initial_condition` - # Amplitude and shift - A = 0.5 - c = 2.0 - - # convenience values for trig. functions - pi_x = pi * x[1] - pi_y = pi * x[2] - pi_t = pi * t - - # compute the manufactured solution and all necessary derivatives - rho = c + A * sin(pi_x) * cos(pi_y) * cos(pi_t) - rho_t = -pi * A * sin(pi_x) * cos(pi_y) * sin(pi_t) - rho_x = pi * A * cos(pi_x) * cos(pi_y) * cos(pi_t) - rho_y = -pi * A * sin(pi_x) * sin(pi_y) * cos(pi_t) - rho_xx = -pi * pi * A * sin(pi_x) * cos(pi_y) * cos(pi_t) - rho_yy = -pi * pi * A * sin(pi_x) * cos(pi_y) * cos(pi_t) - - v1 = sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t) - v1_t = -pi * sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * sin(pi_t) - v1_x = pi * cos(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t) - v1_y = sin(pi_x) * (A * log(y + 2.0) * exp(-A * (y - 1.0)) + (1.0 - exp(-A * (y - 1.0))) / (y + 2.0)) * cos(pi_t) - v1_xx = -pi * pi * sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t) - v1_xy = pi * cos(pi_x) * (A * log(y + 2.0) * exp(-A * (y - 1.0)) + (1.0 - exp(-A * (y - 1.0))) / (y + 2.0)) * cos(pi_t) - v1_yy = (sin(pi_x) * ( 2.0 * A * exp(-A * (y - 1.0)) / (y + 2.0) - - A * A * log(y + 2.0) * exp(-A * (y - 1.0)) - - (1.0 - exp(-A * (y - 1.0))) / ((y + 2.0) * (y + 2.0))) * cos(pi_t)) - v2 = v1 - v2_t = v1_t - v2_x = v1_x - v2_y = v1_y - v2_xx = v1_xx - v2_xy = v1_xy - v2_yy = v1_yy - - p = rho * rho - p_t = 2.0 * rho * rho_t - p_x = 2.0 * rho * rho_x - p_y = 2.0 * rho * rho_y - p_xx = 2.0 * rho * rho_xx + 2.0 * rho_x * rho_x - p_yy = 2.0 * rho * rho_yy + 2.0 * rho_y * rho_y - - # Note this simplifies slightly because the ansatz assumes that v1 = v2 - E = p * inv_gamma_minus_one + 0.5 * rho * (v1^2 + v2^2) - E_t = p_t * inv_gamma_minus_one + rho_t * v1^2 + 2.0 * rho * v1 * v1_t - E_x = p_x * inv_gamma_minus_one + rho_x * v1^2 + 2.0 * rho * v1 * v1_x - E_y = p_y * inv_gamma_minus_one + rho_y * v1^2 + 2.0 * rho * v1 * v1_y - - # Some convenience constants - T_const = equations.gamma * inv_gamma_minus_one / Pr - inv_rho_cubed = 1.0 / (rho^3) - - # compute the source terms - # density equation - du1 = rho_t + rho_x * v1 + rho * v1_x + rho_y * v2 + rho * v2_y - - # x-momentum equation - du2 = ( rho_t * v1 + rho * v1_t + p_x + rho_x * v1^2 - + 2.0 * rho * v1 * v1_x - + rho_y * v1 * v2 - + rho * v1_y * v2 - + rho * v1 * v2_y - # stress tensor from x-direction - - 4.0 / 3.0 * v1_xx * mu_ - + 2.0 / 3.0 * v2_xy * mu_ - - v1_yy * mu_ - - v2_xy * mu_ ) - # y-momentum equation - du3 = ( rho_t * v2 + rho * v2_t + p_y + rho_x * v1 * v2 - + rho * v1_x * v2 - + rho * v1 * v2_x - + rho_y * v2^2 - + 2.0 * rho * v2 * v2_y - # stress tensor from y-direction - - v1_xy * mu_ - - v2_xx * mu_ - - 4.0 / 3.0 * v2_yy * mu_ - + 2.0 / 3.0 * v1_xy * mu_ ) - # total energy equation - du4 = ( E_t + v1_x * (E + p) + v1 * (E_x + p_x) - + v2_y * (E + p) + v2 * (E_y + p_y) - # stress tensor and temperature gradient terms from x-direction - - 4.0 / 3.0 * v1_xx * v1 * mu_ - + 2.0 / 3.0 * v2_xy * v1 * mu_ - - 4.0 / 3.0 * v1_x * v1_x * mu_ - + 2.0 / 3.0 * v2_y * v1_x * mu_ - - v1_xy * v2 * mu_ - - v2_xx * v2 * mu_ - - v1_y * v2_x * mu_ - - v2_x * v2_x * mu_ - - T_const * inv_rho_cubed * ( p_xx * rho * rho - - 2.0 * p_x * rho * rho_x - + 2.0 * p * rho_x * rho_x - - p * rho * rho_xx ) * mu_ - # stress tensor and temperature gradient terms from y-direction - - v1_yy * v1 * mu_ - - v2_xy * v1 * mu_ - - v1_y * v1_y * mu_ - - v2_x * v1_y * mu_ - - 4.0 / 3.0 * v2_yy * v2 * mu_ - + 2.0 / 3.0 * v1_xy * v2 * mu_ - - 4.0 / 3.0 * v2_y * v2_y * mu_ - + 2.0 / 3.0 * v1_x * v2_y * mu_ - - T_const * inv_rho_cubed * ( p_yy * rho * rho - - 2.0 * p_y * rho * rho_y - + 2.0 * p * rho_y * rho_y - - p * rho * rho_yy ) * mu_ ) - - return SVector(du1, du2, du3, du4) -end - -initial_condition = initial_condition_navier_stokes_convergence_test - -# BC types -velocity_bc_top_bottom = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x, t, equations)[2:3]) -heat_bc_top_bottom = Adiabatic((x, t, equations) -> 0.0) -boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom, heat_bc_top_bottom) - -# define inviscid boundary conditions -boundary_conditions = (; :top_bottom => boundary_condition_slip_wall) - -# define viscous boundary conditions -boundary_conditions_parabolic = (; :top_bottom => boundary_condition_top_bottom) - -semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, dg; - boundary_conditions=(boundary_conditions, boundary_conditions_parabolic), - source_terms=source_terms_navier_stokes_convergence_test) - - -############################################################################### -# ODE solvers, callbacks etc. - -# Create ODE problem with time span `tspan` -tspan = (0.0, 0.5) -ode = semidiscretize(semi, tspan) - -summary_callback = SummaryCallback() -alive_callback = AliveCallback(alive_interval=10) -analysis_interval = 100 -analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg)) -callbacks = CallbackSet(summary_callback, alive_callback) - -############################################################################### -# run the simulation - -time_int_tol = 1e-8 -sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol, - ode_default_options()..., callback=callbacks) -summary_callback() # print the timer summary diff --git a/examples/dgmulti_3d/elixir_navierstokes_convergence_curved.jl b/examples/dgmulti_3d/elixir_navierstokes_convergence_curved.jl deleted file mode 100644 index b8971a645e2..00000000000 --- a/examples/dgmulti_3d/elixir_navierstokes_convergence_curved.jl +++ /dev/null @@ -1,261 +0,0 @@ -using OrdinaryDiffEq -using Trixi - -############################################################################### -# semidiscretization of the ideal compressible Navier-Stokes equations - -prandtl_number() = 0.72 -mu() = 0.01 - -equations = CompressibleEulerEquations3D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu=mu(), Prandtl=prandtl_number(), - gradient_variables=GradientVariablesPrimitive()) - -# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux -dg = DGMulti(polydeg = 3, element_type = Hex(), approximation_type = Polynomial(), - surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs), - volume_integral = VolumeIntegralWeakForm()) - -top_bottom(x, tol=50*eps()) = abs(abs(x[2]) - 1) < tol -is_on_boundary = Dict(:top_bottom => top_bottom) - -function mapping(xi, eta, zeta) - x = xi + 0.1 * sin(pi * xi) * sin(pi * eta) - y = eta + 0.1 * sin(pi * xi) * sin(pi * eta) - z = zeta + 0.1 * sin(pi * xi) * sin(pi * eta) - return SVector(x, y, z) -end -cells_per_dimension = (8, 8, 8) -mesh = DGMultiMesh(dg, cells_per_dimension, mapping; periodicity=(true, false, true), is_on_boundary) - -# Note: the initial condition cannot be specialized to `CompressibleNavierStokesDiffusion3D` -# since it is called by both the parabolic solver (which passes in `CompressibleNavierStokesDiffusion3D`) -# and by the initial condition (which passes in `CompressibleEulerEquations3D`). -# This convergence test setup was originally derived by Andrew Winters (@andrewwinters5000) -function initial_condition_navier_stokes_convergence_test(x, t, equations) - # Constants. OBS! Must match those in `source_terms_navier_stokes_convergence_test` - c = 2.0 - A1 = 0.5 - A2 = 1.0 - A3 = 0.5 - - # Convenience values for trig. functions - pi_x = pi * x[1] - pi_y = pi * x[2] - pi_z = pi * x[3] - pi_t = pi * t - - rho = c + A1 * sin(pi_x) * cos(pi_y) * sin(pi_z) * cos(pi_t) - v1 = A2 * sin(pi_x) * log(x[2] + 2.0) * (1.0 - exp(-A3 * (x[2] - 1.0))) * sin(pi_z) * cos(pi_t) - v2 = v1 - v3 = v1 - p = rho^2 - - return prim2cons(SVector(rho, v1, v2, v3, p), equations) -end - -@inline function source_terms_navier_stokes_convergence_test(u, x, t, equations) - # TODO: parabolic - # we currently need to hardcode these parameters until we fix the "combined equation" issue - # see also https://github.com/trixi-framework/Trixi.jl/pull/1160 - inv_gamma_minus_one = inv(equations.gamma - 1) - Pr = prandtl_number() - mu_ = mu() - - # Constants. OBS! Must match those in `initial_condition_navier_stokes_convergence_test` - c = 2.0 - A1 = 0.5 - A2 = 1.0 - A3 = 0.5 - - # Convenience values for trig. functions - pi_x = pi * x[1] - pi_y = pi * x[2] - pi_z = pi * x[3] - pi_t = pi * t - - # Define auxiliary functions for the strange function of the y variable - # to make expressions easier to read - g = log(x[2] + 2.0) * (1.0 - exp(-A3 * (x[2] - 1.0))) - g_y = ( A3 * log(x[2] + 2.0) * exp(-A3 * (x[2] - 1.0)) - + (1.0 - exp(-A3 * (x[2] - 1.0))) / (x[2] + 2.0) ) - g_yy = ( 2.0 * A3 * exp(-A3 * (x[2] - 1.0)) / (x[2] + 2.0) - - (1.0 - exp(-A3 * (x[2] - 1.0))) / ((x[2] + 2.0)^2) - - A3^2 * log(x[2] + 2.0) * exp(-A3 * (x[2] - 1.0)) ) - - # Density and its derivatives - rho = c + A1 * sin(pi_x) * cos(pi_y) * sin(pi_z) * cos(pi_t) - rho_t = -pi * A1 * sin(pi_x) * cos(pi_y) * sin(pi_z) * sin(pi_t) - rho_x = pi * A1 * cos(pi_x) * cos(pi_y) * sin(pi_z) * cos(pi_t) - rho_y = -pi * A1 * sin(pi_x) * sin(pi_y) * sin(pi_z) * cos(pi_t) - rho_z = pi * A1 * sin(pi_x) * cos(pi_y) * cos(pi_z) * cos(pi_t) - rho_xx = -pi^2 * (rho - c) - rho_yy = -pi^2 * (rho - c) - rho_zz = -pi^2 * (rho - c) - - # Velocities and their derivatives - # v1 terms - v1 = A2 * sin(pi_x) * g * sin(pi_z) * cos(pi_t) - v1_t = -pi * A2 * sin(pi_x) * g * sin(pi_z) * sin(pi_t) - v1_x = pi * A2 * cos(pi_x) * g * sin(pi_z) * cos(pi_t) - v1_y = A2 * sin(pi_x) * g_y * sin(pi_z) * cos(pi_t) - v1_z = pi * A2 * sin(pi_x) * g * cos(pi_z) * cos(pi_t) - v1_xx = -pi^2 * v1 - v1_yy = A2 * sin(pi_x) * g_yy * sin(pi_z) * cos(pi_t) - v1_zz = -pi^2 * v1 - v1_xy = pi * A2 * cos(pi_x) * g_y * sin(pi_z) * cos(pi_t) - v1_xz = pi^2 * A2 * cos(pi_x) * g * cos(pi_z) * cos(pi_t) - v1_yz = pi * A2 * sin(pi_x) * g_y * cos(pi_z) * cos(pi_t) - # v2 terms (simplifies from ansatz) - v2 = v1 - v2_t = v1_t - v2_x = v1_x - v2_y = v1_y - v2_z = v1_z - v2_xx = v1_xx - v2_yy = v1_yy - v2_zz = v1_zz - v2_xy = v1_xy - v2_yz = v1_yz - # v3 terms (simplifies from ansatz) - v3 = v1 - v3_t = v1_t - v3_x = v1_x - v3_y = v1_y - v3_z = v1_z - v3_xx = v1_xx - v3_yy = v1_yy - v3_zz = v1_zz - v3_xz = v1_xz - v3_yz = v1_yz - - # Pressure and its derivatives - p = rho^2 - p_t = 2.0 * rho * rho_t - p_x = 2.0 * rho * rho_x - p_y = 2.0 * rho * rho_y - p_z = 2.0 * rho * rho_z - - # Total energy and its derivatives; simiplifies from ansatz that v2 = v1 and v3 = v1 - E = p * inv_gamma_minus_one + 1.5 * rho * v1^2 - E_t = p_t * inv_gamma_minus_one + 1.5 * rho_t * v1^2 + 3.0 * rho * v1 * v1_t - E_x = p_x * inv_gamma_minus_one + 1.5 * rho_x * v1^2 + 3.0 * rho * v1 * v1_x - E_y = p_y * inv_gamma_minus_one + 1.5 * rho_y * v1^2 + 3.0 * rho * v1 * v1_y - E_z = p_z * inv_gamma_minus_one + 1.5 * rho_z * v1^2 + 3.0 * rho * v1 * v1_z - - # Divergence of Fick's law ∇⋅∇q = kappa ∇⋅∇T; simplifies because p = rho², so T = p/rho = rho - kappa = equations.gamma * inv_gamma_minus_one / Pr - q_xx = kappa * rho_xx # kappa T_xx - q_yy = kappa * rho_yy # kappa T_yy - q_zz = kappa * rho_zz # kappa T_zz - - # Stress tensor and its derivatives (exploit symmetry) - tau11 = 4.0 / 3.0 * v1_x - 2.0 / 3.0 * (v2_y + v3_z) - tau12 = v1_y + v2_x - tau13 = v1_z + v3_x - tau22 = 4.0 / 3.0 * v2_y - 2.0 / 3.0 * (v1_x + v3_z) - tau23 = v2_z + v3_y - tau33 = 4.0 / 3.0 * v3_z - 2.0 / 3.0 * (v1_x + v2_y) - - tau11_x = 4.0 / 3.0 * v1_xx - 2.0 / 3.0 * (v2_xy + v3_xz) - tau12_x = v1_xy + v2_xx - tau13_x = v1_xz + v3_xx - - tau12_y = v1_yy + v2_xy - tau22_y = 4.0 / 3.0 * v2_yy - 2.0 / 3.0 * (v1_xy + v3_yz) - tau23_y = v2_yz + v3_yy - - tau13_z = v1_zz + v3_xz - tau23_z = v2_zz + v3_yz - tau33_z = 4.0 / 3.0 * v3_zz - 2.0 / 3.0 * (v1_xz + v2_yz) - - # Compute the source terms - # Density equation - du1 = ( rho_t + rho_x * v1 + rho * v1_x - + rho_y * v2 + rho * v2_y - + rho_z * v3 + rho * v3_z ) - # x-momentum equation - du2 = ( rho_t * v1 + rho * v1_t + p_x + rho_x * v1^2 - + 2.0 * rho * v1 * v1_x - + rho_y * v1 * v2 - + rho * v1_y * v2 - + rho * v1 * v2_y - + rho_z * v1 * v3 - + rho * v1_z * v3 - + rho * v1 * v3_z - - mu_ * (tau11_x + tau12_y + tau13_z) ) - # y-momentum equation - du3 = ( rho_t * v2 + rho * v2_t + p_y + rho_x * v1 * v2 - + rho * v1_x * v2 - + rho * v1 * v2_x - + rho_y * v2^2 - + 2.0 * rho * v2 * v2_y - + rho_z * v2 * v3 - + rho * v2_z * v3 - + rho * v2 * v3_z - - mu_ * (tau12_x + tau22_y + tau23_z) ) - # z-momentum equation - du4 = ( rho_t * v3 + rho * v3_t + p_z + rho_x * v1 * v3 - + rho * v1_x * v3 - + rho * v1 * v3_x - + rho_y * v2 * v3 - + rho * v2_y * v3 - + rho * v2 * v3_y - + rho_z * v3^2 - + 2.0 * rho * v3 * v3_z - - mu_ * (tau13_x + tau23_y + tau33_z) ) - # Total energy equation - du5 = ( E_t + v1_x * (E + p) + v1 * (E_x + p_x) - + v2_y * (E + p) + v2 * (E_y + p_y) - + v3_z * (E + p) + v3 * (E_z + p_z) - # stress tensor and temperature gradient from x-direction - - mu_ * ( q_xx + v1_x * tau11 + v2_x * tau12 + v3_x * tau13 - + v1 * tau11_x + v2 * tau12_x + v3 * tau13_x) - # stress tensor and temperature gradient terms from y-direction - - mu_ * ( q_yy + v1_y * tau12 + v2_y * tau22 + v3_y * tau23 - + v1 * tau12_y + v2 * tau22_y + v3 * tau23_y) - # stress tensor and temperature gradient terms from z-direction - - mu_ * ( q_zz + v1_z * tau13 + v2_z * tau23 + v3_z * tau33 - + v1 * tau13_z + v2 * tau23_z + v3 * tau33_z) ) - - return SVector(du1, du2, du3, du4, du5) -end - -initial_condition = initial_condition_navier_stokes_convergence_test - -# BC types -velocity_bc_top_bottom = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x, t, equations)[2:4]) -heat_bc_top_bottom = Adiabatic((x, t, equations) -> 0.0) -boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom, heat_bc_top_bottom) - -# define inviscid boundary conditions -boundary_conditions = (; :top_bottom => boundary_condition_slip_wall) - -# define viscous boundary conditions -boundary_conditions_parabolic = (; :top_bottom => boundary_condition_top_bottom) - -semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, dg; - boundary_conditions=(boundary_conditions, boundary_conditions_parabolic), - source_terms=source_terms_navier_stokes_convergence_test) - - -############################################################################### -# ODE solvers, callbacks etc. - -# Create ODE problem with time span `tspan` -tspan = (0.0, 1.0) -ode = semidiscretize(semi, tspan) - -summary_callback = SummaryCallback() -alive_callback = AliveCallback(alive_interval=10) -analysis_interval = 100 -analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg)) -callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback) - -############################################################################### -# run the simulation - -time_int_tol = 1e-8 -sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol, - ode_default_options()..., callback=callbacks) -summary_callback() # print the timer summary diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 588f43e4543..12d579008d1 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -101,9 +101,14 @@ isdir(outdir) && rm(outdir, recursive=true) ) end - @trixi_testset "DGMulti: elixir_navierstokes_convergence_curved.jl" begin - @test_trixi_include(joinpath(examples_dir(), "dgmulti_2d", "elixir_navierstokes_convergence_curved.jl"), - cells_per_dimension = (4, 4), tspan=(0.0, 0.1), + @trixi_testset "DGMulti: elixir_navierstokes_convergence.jl (curved)" begin + function mapping(xi, eta) + x = xi + 0.1 * sin(pi * xi) * sin(pi * eta) + y = eta + 0.1 * sin(pi * xi) * sin(pi * eta) + return SVector(x, y) + end + @test_trixi_include(joinpath(examples_dir(), "dgmulti_2d", "elixir_navierstokes_convergence.jl"), + cells_per_dimension = (4, 4), tspan=(0.0, 0.1), mesh = DGMultiMesh(dg, cells_per_dimension, mapping; periodicity=(true, false), is_on_boundary), l2 = [0.004255101916146187, 0.011118488923215765, 0.011281831283462686, 0.03573656447388509], linf = [0.015071710669706473, 0.04103132025858458, 0.03990424085750277, 0.1309401718598764], ) diff --git a/test/test_parabolic_3d.jl b/test/test_parabolic_3d.jl index 1ae5eed44ae..210e571e8d4 100644 --- a/test/test_parabolic_3d.jl +++ b/test/test_parabolic_3d.jl @@ -19,9 +19,15 @@ isdir(outdir) && rm(outdir, recursive=true) ) end - @trixi_testset "DGMulti: elixir_navierstokes_convergence_curved.jl" begin - @test_trixi_include(joinpath(examples_dir(), "dgmulti_3d", "elixir_navierstokes_convergence_curved.jl"), - cells_per_dimension = (4, 4, 4), tspan=(0.0, 0.1), + @trixi_testset "DGMulti: elixir_navierstokes_convergence.jl (curved)" begin + function mapping(xi, eta, zeta) + x = xi + 0.1 * sin(pi * xi) * sin(pi * eta) + y = eta + 0.1 * sin(pi * xi) * sin(pi * eta) + z = zeta + 0.1 * sin(pi * xi) * sin(pi * eta) + return SVector(x, y, z) + end + @test_trixi_include(joinpath(examples_dir(), "dgmulti_3d", "elixir_navierstokes_convergence.jl"), + cells_per_dimension = (4, 4, 4), tspan=(0.0, 0.1), mesh = DGMultiMesh(dg, cells_per_dimension, mapping; periodicity=(true, false, true), is_on_boundary), l2 = [0.0014027227251207474, 0.0021322235533273513, 0.0027873741447455194, 0.0024587473070627423, 0.00997836818019202], linf = [0.006341750402837576, 0.010306014252246865, 0.01520740250924979, 0.010968264045485565, 0.047454389831591115] ) From d1a42f2569d3cbf40a1d6d49a0e4c0de0836338f Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 27 Apr 2023 22:08:19 -0500 Subject: [PATCH 33/36] revert removal of "inv_h" in `penalty` this makes it so this PR is no longer a breaking change --- src/equations/laplace_diffusion_2d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/equations/laplace_diffusion_2d.jl b/src/equations/laplace_diffusion_2d.jl index 1dc843ca6c7..2f1afe25a6d 100644 --- a/src/equations/laplace_diffusion_2d.jl +++ b/src/equations/laplace_diffusion_2d.jl @@ -28,8 +28,8 @@ end # TODO: parabolic; should this remain in the equations file, be moved to solvers, or live in the elixir? # 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, equations_parabolic::LaplaceDiffusion2D, dg::ViscousFormulationLocalDG) - return dg.penalty_parameter * (u_outer - u_inner) * equations_parabolic.diffusivity +function penalty(u_outer, u_inner, inv_h, equations_parabolic::LaplaceDiffusion2D, dg::ViscousFormulationLocalDG) + return dg.penalty_parameter * (u_outer - u_inner) * equations_parabolic.diffusivity * inv_h end # Dirichlet-type boundary condition for use with a parabolic solver in weak form From 02c7fe9a714973671fe6474a979ad64acfdad705 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 27 Apr 2023 22:05:20 -0500 Subject: [PATCH 34/36] Revert "replace extra curved elixirs with "trixi_include" calls" This reverts commit adc0fdec01d75c45c82f3984d6cbf728e846cf60. --- .../elixir_navierstokes_convergence_curved.jl | 212 ++++++++++++++ .../elixir_navierstokes_convergence_curved.jl | 261 ++++++++++++++++++ test/test_parabolic_2d.jl | 11 +- test/test_parabolic_3d.jl | 12 +- 4 files changed, 479 insertions(+), 17 deletions(-) create mode 100644 examples/dgmulti_2d/elixir_navierstokes_convergence_curved.jl create mode 100644 examples/dgmulti_3d/elixir_navierstokes_convergence_curved.jl diff --git a/examples/dgmulti_2d/elixir_navierstokes_convergence_curved.jl b/examples/dgmulti_2d/elixir_navierstokes_convergence_curved.jl new file mode 100644 index 00000000000..854c2af773d --- /dev/null +++ b/examples/dgmulti_2d/elixir_navierstokes_convergence_curved.jl @@ -0,0 +1,212 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the ideal compressible Navier-Stokes equations + +prandtl_number() = 0.72 +mu() = 0.01 + +equations = CompressibleEulerEquations2D(1.4) +# Note: If you change the Navier-Stokes parameters here, also change them in the initial condition +# I really do not like this structure but it should work for now +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu=mu(), Prandtl=prandtl_number(), + gradient_variables=GradientVariablesPrimitive()) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial(), + surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs), + volume_integral = VolumeIntegralWeakForm()) + +top_bottom(x, tol=50*eps()) = abs(abs(x[2]) - 1) < tol +is_on_boundary = Dict(:top_bottom => top_bottom) + +function mapping(xi, eta) + x = xi + 0.1 * sin(pi * xi) * sin(pi * eta) + y = eta + 0.1 * sin(pi * xi) * sin(pi * eta) + return SVector(x, y) +end +cells_per_dimension = (16, 16) +mesh = DGMultiMesh(dg, cells_per_dimension, mapping; periodicity=(true, false), is_on_boundary) + +# Note: the initial condition cannot be specialized to `CompressibleNavierStokesDiffusion2D` +# since it is called by both the parabolic solver (which passes in `CompressibleNavierStokesDiffusion2D`) +# and by the initial condition (which passes in `CompressibleEulerEquations2D`). +# This convergence test setup was originally derived by Andrew Winters (@andrewwinters5000) +function initial_condition_navier_stokes_convergence_test(x, t, equations) + # Amplitude and shift + A = 0.5 + c = 2.0 + + # convenience values for trig. functions + pi_x = pi * x[1] + pi_y = pi * x[2] + pi_t = pi * t + + rho = c + A * sin(pi_x) * cos(pi_y) * cos(pi_t) + v1 = sin(pi_x) * log(x[2] + 2.0) * (1.0 - exp(-A * (x[2] - 1.0)) ) * cos(pi_t) + v2 = v1 + p = rho^2 + + return prim2cons(SVector(rho, v1, v2, p), equations) +end + +@inline function source_terms_navier_stokes_convergence_test(u, x, t, equations) + y = x[2] + + # TODO: parabolic + # we currently need to hardcode these parameters until we fix the "combined equation" issue + # see also https://github.com/trixi-framework/Trixi.jl/pull/1160 + inv_gamma_minus_one = inv(equations.gamma - 1) + Pr = prandtl_number() + mu_ = mu() + + # Same settings as in `initial_condition` + # Amplitude and shift + A = 0.5 + c = 2.0 + + # convenience values for trig. functions + pi_x = pi * x[1] + pi_y = pi * x[2] + pi_t = pi * t + + # compute the manufactured solution and all necessary derivatives + rho = c + A * sin(pi_x) * cos(pi_y) * cos(pi_t) + rho_t = -pi * A * sin(pi_x) * cos(pi_y) * sin(pi_t) + rho_x = pi * A * cos(pi_x) * cos(pi_y) * cos(pi_t) + rho_y = -pi * A * sin(pi_x) * sin(pi_y) * cos(pi_t) + rho_xx = -pi * pi * A * sin(pi_x) * cos(pi_y) * cos(pi_t) + rho_yy = -pi * pi * A * sin(pi_x) * cos(pi_y) * cos(pi_t) + + v1 = sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t) + v1_t = -pi * sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * sin(pi_t) + v1_x = pi * cos(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t) + v1_y = sin(pi_x) * (A * log(y + 2.0) * exp(-A * (y - 1.0)) + (1.0 - exp(-A * (y - 1.0))) / (y + 2.0)) * cos(pi_t) + v1_xx = -pi * pi * sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t) + v1_xy = pi * cos(pi_x) * (A * log(y + 2.0) * exp(-A * (y - 1.0)) + (1.0 - exp(-A * (y - 1.0))) / (y + 2.0)) * cos(pi_t) + v1_yy = (sin(pi_x) * ( 2.0 * A * exp(-A * (y - 1.0)) / (y + 2.0) + - A * A * log(y + 2.0) * exp(-A * (y - 1.0)) + - (1.0 - exp(-A * (y - 1.0))) / ((y + 2.0) * (y + 2.0))) * cos(pi_t)) + v2 = v1 + v2_t = v1_t + v2_x = v1_x + v2_y = v1_y + v2_xx = v1_xx + v2_xy = v1_xy + v2_yy = v1_yy + + p = rho * rho + p_t = 2.0 * rho * rho_t + p_x = 2.0 * rho * rho_x + p_y = 2.0 * rho * rho_y + p_xx = 2.0 * rho * rho_xx + 2.0 * rho_x * rho_x + p_yy = 2.0 * rho * rho_yy + 2.0 * rho_y * rho_y + + # Note this simplifies slightly because the ansatz assumes that v1 = v2 + E = p * inv_gamma_minus_one + 0.5 * rho * (v1^2 + v2^2) + E_t = p_t * inv_gamma_minus_one + rho_t * v1^2 + 2.0 * rho * v1 * v1_t + E_x = p_x * inv_gamma_minus_one + rho_x * v1^2 + 2.0 * rho * v1 * v1_x + E_y = p_y * inv_gamma_minus_one + rho_y * v1^2 + 2.0 * rho * v1 * v1_y + + # Some convenience constants + T_const = equations.gamma * inv_gamma_minus_one / Pr + inv_rho_cubed = 1.0 / (rho^3) + + # compute the source terms + # density equation + du1 = rho_t + rho_x * v1 + rho * v1_x + rho_y * v2 + rho * v2_y + + # x-momentum equation + du2 = ( rho_t * v1 + rho * v1_t + p_x + rho_x * v1^2 + + 2.0 * rho * v1 * v1_x + + rho_y * v1 * v2 + + rho * v1_y * v2 + + rho * v1 * v2_y + # stress tensor from x-direction + - 4.0 / 3.0 * v1_xx * mu_ + + 2.0 / 3.0 * v2_xy * mu_ + - v1_yy * mu_ + - v2_xy * mu_ ) + # y-momentum equation + du3 = ( rho_t * v2 + rho * v2_t + p_y + rho_x * v1 * v2 + + rho * v1_x * v2 + + rho * v1 * v2_x + + rho_y * v2^2 + + 2.0 * rho * v2 * v2_y + # stress tensor from y-direction + - v1_xy * mu_ + - v2_xx * mu_ + - 4.0 / 3.0 * v2_yy * mu_ + + 2.0 / 3.0 * v1_xy * mu_ ) + # total energy equation + du4 = ( E_t + v1_x * (E + p) + v1 * (E_x + p_x) + + v2_y * (E + p) + v2 * (E_y + p_y) + # stress tensor and temperature gradient terms from x-direction + - 4.0 / 3.0 * v1_xx * v1 * mu_ + + 2.0 / 3.0 * v2_xy * v1 * mu_ + - 4.0 / 3.0 * v1_x * v1_x * mu_ + + 2.0 / 3.0 * v2_y * v1_x * mu_ + - v1_xy * v2 * mu_ + - v2_xx * v2 * mu_ + - v1_y * v2_x * mu_ + - v2_x * v2_x * mu_ + - T_const * inv_rho_cubed * ( p_xx * rho * rho + - 2.0 * p_x * rho * rho_x + + 2.0 * p * rho_x * rho_x + - p * rho * rho_xx ) * mu_ + # stress tensor and temperature gradient terms from y-direction + - v1_yy * v1 * mu_ + - v2_xy * v1 * mu_ + - v1_y * v1_y * mu_ + - v2_x * v1_y * mu_ + - 4.0 / 3.0 * v2_yy * v2 * mu_ + + 2.0 / 3.0 * v1_xy * v2 * mu_ + - 4.0 / 3.0 * v2_y * v2_y * mu_ + + 2.0 / 3.0 * v1_x * v2_y * mu_ + - T_const * inv_rho_cubed * ( p_yy * rho * rho + - 2.0 * p_y * rho * rho_y + + 2.0 * p * rho_y * rho_y + - p * rho * rho_yy ) * mu_ ) + + return SVector(du1, du2, du3, du4) +end + +initial_condition = initial_condition_navier_stokes_convergence_test + +# BC types +velocity_bc_top_bottom = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x, t, equations)[2:3]) +heat_bc_top_bottom = Adiabatic((x, t, equations) -> 0.0) +boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom, heat_bc_top_bottom) + +# define inviscid boundary conditions +boundary_conditions = (; :top_bottom => boundary_condition_slip_wall) + +# define viscous boundary conditions +boundary_conditions_parabolic = (; :top_bottom => boundary_condition_top_bottom) + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, dg; + boundary_conditions=(boundary_conditions, boundary_conditions_parabolic), + source_terms=source_terms_navier_stokes_convergence_test) + + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span `tspan` +tspan = (0.0, 0.5) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() +alive_callback = AliveCallback(alive_interval=10) +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg)) +callbacks = CallbackSet(summary_callback, alive_callback) + +############################################################################### +# run the simulation + +time_int_tol = 1e-8 +sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol, + ode_default_options()..., callback=callbacks) +summary_callback() # print the timer summary diff --git a/examples/dgmulti_3d/elixir_navierstokes_convergence_curved.jl b/examples/dgmulti_3d/elixir_navierstokes_convergence_curved.jl new file mode 100644 index 00000000000..b8971a645e2 --- /dev/null +++ b/examples/dgmulti_3d/elixir_navierstokes_convergence_curved.jl @@ -0,0 +1,261 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the ideal compressible Navier-Stokes equations + +prandtl_number() = 0.72 +mu() = 0.01 + +equations = CompressibleEulerEquations3D(1.4) +equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu=mu(), Prandtl=prandtl_number(), + gradient_variables=GradientVariablesPrimitive()) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +dg = DGMulti(polydeg = 3, element_type = Hex(), approximation_type = Polynomial(), + surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs), + volume_integral = VolumeIntegralWeakForm()) + +top_bottom(x, tol=50*eps()) = abs(abs(x[2]) - 1) < tol +is_on_boundary = Dict(:top_bottom => top_bottom) + +function mapping(xi, eta, zeta) + x = xi + 0.1 * sin(pi * xi) * sin(pi * eta) + y = eta + 0.1 * sin(pi * xi) * sin(pi * eta) + z = zeta + 0.1 * sin(pi * xi) * sin(pi * eta) + return SVector(x, y, z) +end +cells_per_dimension = (8, 8, 8) +mesh = DGMultiMesh(dg, cells_per_dimension, mapping; periodicity=(true, false, true), is_on_boundary) + +# Note: the initial condition cannot be specialized to `CompressibleNavierStokesDiffusion3D` +# since it is called by both the parabolic solver (which passes in `CompressibleNavierStokesDiffusion3D`) +# and by the initial condition (which passes in `CompressibleEulerEquations3D`). +# This convergence test setup was originally derived by Andrew Winters (@andrewwinters5000) +function initial_condition_navier_stokes_convergence_test(x, t, equations) + # Constants. OBS! Must match those in `source_terms_navier_stokes_convergence_test` + c = 2.0 + A1 = 0.5 + A2 = 1.0 + A3 = 0.5 + + # Convenience values for trig. functions + pi_x = pi * x[1] + pi_y = pi * x[2] + pi_z = pi * x[3] + pi_t = pi * t + + rho = c + A1 * sin(pi_x) * cos(pi_y) * sin(pi_z) * cos(pi_t) + v1 = A2 * sin(pi_x) * log(x[2] + 2.0) * (1.0 - exp(-A3 * (x[2] - 1.0))) * sin(pi_z) * cos(pi_t) + v2 = v1 + v3 = v1 + p = rho^2 + + return prim2cons(SVector(rho, v1, v2, v3, p), equations) +end + +@inline function source_terms_navier_stokes_convergence_test(u, x, t, equations) + # TODO: parabolic + # we currently need to hardcode these parameters until we fix the "combined equation" issue + # see also https://github.com/trixi-framework/Trixi.jl/pull/1160 + inv_gamma_minus_one = inv(equations.gamma - 1) + Pr = prandtl_number() + mu_ = mu() + + # Constants. OBS! Must match those in `initial_condition_navier_stokes_convergence_test` + c = 2.0 + A1 = 0.5 + A2 = 1.0 + A3 = 0.5 + + # Convenience values for trig. functions + pi_x = pi * x[1] + pi_y = pi * x[2] + pi_z = pi * x[3] + pi_t = pi * t + + # Define auxiliary functions for the strange function of the y variable + # to make expressions easier to read + g = log(x[2] + 2.0) * (1.0 - exp(-A3 * (x[2] - 1.0))) + g_y = ( A3 * log(x[2] + 2.0) * exp(-A3 * (x[2] - 1.0)) + + (1.0 - exp(-A3 * (x[2] - 1.0))) / (x[2] + 2.0) ) + g_yy = ( 2.0 * A3 * exp(-A3 * (x[2] - 1.0)) / (x[2] + 2.0) + - (1.0 - exp(-A3 * (x[2] - 1.0))) / ((x[2] + 2.0)^2) + - A3^2 * log(x[2] + 2.0) * exp(-A3 * (x[2] - 1.0)) ) + + # Density and its derivatives + rho = c + A1 * sin(pi_x) * cos(pi_y) * sin(pi_z) * cos(pi_t) + rho_t = -pi * A1 * sin(pi_x) * cos(pi_y) * sin(pi_z) * sin(pi_t) + rho_x = pi * A1 * cos(pi_x) * cos(pi_y) * sin(pi_z) * cos(pi_t) + rho_y = -pi * A1 * sin(pi_x) * sin(pi_y) * sin(pi_z) * cos(pi_t) + rho_z = pi * A1 * sin(pi_x) * cos(pi_y) * cos(pi_z) * cos(pi_t) + rho_xx = -pi^2 * (rho - c) + rho_yy = -pi^2 * (rho - c) + rho_zz = -pi^2 * (rho - c) + + # Velocities and their derivatives + # v1 terms + v1 = A2 * sin(pi_x) * g * sin(pi_z) * cos(pi_t) + v1_t = -pi * A2 * sin(pi_x) * g * sin(pi_z) * sin(pi_t) + v1_x = pi * A2 * cos(pi_x) * g * sin(pi_z) * cos(pi_t) + v1_y = A2 * sin(pi_x) * g_y * sin(pi_z) * cos(pi_t) + v1_z = pi * A2 * sin(pi_x) * g * cos(pi_z) * cos(pi_t) + v1_xx = -pi^2 * v1 + v1_yy = A2 * sin(pi_x) * g_yy * sin(pi_z) * cos(pi_t) + v1_zz = -pi^2 * v1 + v1_xy = pi * A2 * cos(pi_x) * g_y * sin(pi_z) * cos(pi_t) + v1_xz = pi^2 * A2 * cos(pi_x) * g * cos(pi_z) * cos(pi_t) + v1_yz = pi * A2 * sin(pi_x) * g_y * cos(pi_z) * cos(pi_t) + # v2 terms (simplifies from ansatz) + v2 = v1 + v2_t = v1_t + v2_x = v1_x + v2_y = v1_y + v2_z = v1_z + v2_xx = v1_xx + v2_yy = v1_yy + v2_zz = v1_zz + v2_xy = v1_xy + v2_yz = v1_yz + # v3 terms (simplifies from ansatz) + v3 = v1 + v3_t = v1_t + v3_x = v1_x + v3_y = v1_y + v3_z = v1_z + v3_xx = v1_xx + v3_yy = v1_yy + v3_zz = v1_zz + v3_xz = v1_xz + v3_yz = v1_yz + + # Pressure and its derivatives + p = rho^2 + p_t = 2.0 * rho * rho_t + p_x = 2.0 * rho * rho_x + p_y = 2.0 * rho * rho_y + p_z = 2.0 * rho * rho_z + + # Total energy and its derivatives; simiplifies from ansatz that v2 = v1 and v3 = v1 + E = p * inv_gamma_minus_one + 1.5 * rho * v1^2 + E_t = p_t * inv_gamma_minus_one + 1.5 * rho_t * v1^2 + 3.0 * rho * v1 * v1_t + E_x = p_x * inv_gamma_minus_one + 1.5 * rho_x * v1^2 + 3.0 * rho * v1 * v1_x + E_y = p_y * inv_gamma_minus_one + 1.5 * rho_y * v1^2 + 3.0 * rho * v1 * v1_y + E_z = p_z * inv_gamma_minus_one + 1.5 * rho_z * v1^2 + 3.0 * rho * v1 * v1_z + + # Divergence of Fick's law ∇⋅∇q = kappa ∇⋅∇T; simplifies because p = rho², so T = p/rho = rho + kappa = equations.gamma * inv_gamma_minus_one / Pr + q_xx = kappa * rho_xx # kappa T_xx + q_yy = kappa * rho_yy # kappa T_yy + q_zz = kappa * rho_zz # kappa T_zz + + # Stress tensor and its derivatives (exploit symmetry) + tau11 = 4.0 / 3.0 * v1_x - 2.0 / 3.0 * (v2_y + v3_z) + tau12 = v1_y + v2_x + tau13 = v1_z + v3_x + tau22 = 4.0 / 3.0 * v2_y - 2.0 / 3.0 * (v1_x + v3_z) + tau23 = v2_z + v3_y + tau33 = 4.0 / 3.0 * v3_z - 2.0 / 3.0 * (v1_x + v2_y) + + tau11_x = 4.0 / 3.0 * v1_xx - 2.0 / 3.0 * (v2_xy + v3_xz) + tau12_x = v1_xy + v2_xx + tau13_x = v1_xz + v3_xx + + tau12_y = v1_yy + v2_xy + tau22_y = 4.0 / 3.0 * v2_yy - 2.0 / 3.0 * (v1_xy + v3_yz) + tau23_y = v2_yz + v3_yy + + tau13_z = v1_zz + v3_xz + tau23_z = v2_zz + v3_yz + tau33_z = 4.0 / 3.0 * v3_zz - 2.0 / 3.0 * (v1_xz + v2_yz) + + # Compute the source terms + # Density equation + du1 = ( rho_t + rho_x * v1 + rho * v1_x + + rho_y * v2 + rho * v2_y + + rho_z * v3 + rho * v3_z ) + # x-momentum equation + du2 = ( rho_t * v1 + rho * v1_t + p_x + rho_x * v1^2 + + 2.0 * rho * v1 * v1_x + + rho_y * v1 * v2 + + rho * v1_y * v2 + + rho * v1 * v2_y + + rho_z * v1 * v3 + + rho * v1_z * v3 + + rho * v1 * v3_z + - mu_ * (tau11_x + tau12_y + tau13_z) ) + # y-momentum equation + du3 = ( rho_t * v2 + rho * v2_t + p_y + rho_x * v1 * v2 + + rho * v1_x * v2 + + rho * v1 * v2_x + + rho_y * v2^2 + + 2.0 * rho * v2 * v2_y + + rho_z * v2 * v3 + + rho * v2_z * v3 + + rho * v2 * v3_z + - mu_ * (tau12_x + tau22_y + tau23_z) ) + # z-momentum equation + du4 = ( rho_t * v3 + rho * v3_t + p_z + rho_x * v1 * v3 + + rho * v1_x * v3 + + rho * v1 * v3_x + + rho_y * v2 * v3 + + rho * v2_y * v3 + + rho * v2 * v3_y + + rho_z * v3^2 + + 2.0 * rho * v3 * v3_z + - mu_ * (tau13_x + tau23_y + tau33_z) ) + # Total energy equation + du5 = ( E_t + v1_x * (E + p) + v1 * (E_x + p_x) + + v2_y * (E + p) + v2 * (E_y + p_y) + + v3_z * (E + p) + v3 * (E_z + p_z) + # stress tensor and temperature gradient from x-direction + - mu_ * ( q_xx + v1_x * tau11 + v2_x * tau12 + v3_x * tau13 + + v1 * tau11_x + v2 * tau12_x + v3 * tau13_x) + # stress tensor and temperature gradient terms from y-direction + - mu_ * ( q_yy + v1_y * tau12 + v2_y * tau22 + v3_y * tau23 + + v1 * tau12_y + v2 * tau22_y + v3 * tau23_y) + # stress tensor and temperature gradient terms from z-direction + - mu_ * ( q_zz + v1_z * tau13 + v2_z * tau23 + v3_z * tau33 + + v1 * tau13_z + v2 * tau23_z + v3 * tau33_z) ) + + return SVector(du1, du2, du3, du4, du5) +end + +initial_condition = initial_condition_navier_stokes_convergence_test + +# BC types +velocity_bc_top_bottom = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x, t, equations)[2:4]) +heat_bc_top_bottom = Adiabatic((x, t, equations) -> 0.0) +boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom, heat_bc_top_bottom) + +# define inviscid boundary conditions +boundary_conditions = (; :top_bottom => boundary_condition_slip_wall) + +# define viscous boundary conditions +boundary_conditions_parabolic = (; :top_bottom => boundary_condition_top_bottom) + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, dg; + boundary_conditions=(boundary_conditions, boundary_conditions_parabolic), + source_terms=source_terms_navier_stokes_convergence_test) + + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span `tspan` +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() +alive_callback = AliveCallback(alive_interval=10) +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg)) +callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback) + +############################################################################### +# run the simulation + +time_int_tol = 1e-8 +sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol, + ode_default_options()..., callback=callbacks) +summary_callback() # print the timer summary diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 12d579008d1..588f43e4543 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -101,14 +101,9 @@ isdir(outdir) && rm(outdir, recursive=true) ) end - @trixi_testset "DGMulti: elixir_navierstokes_convergence.jl (curved)" begin - function mapping(xi, eta) - x = xi + 0.1 * sin(pi * xi) * sin(pi * eta) - y = eta + 0.1 * sin(pi * xi) * sin(pi * eta) - return SVector(x, y) - end - @test_trixi_include(joinpath(examples_dir(), "dgmulti_2d", "elixir_navierstokes_convergence.jl"), - cells_per_dimension = (4, 4), tspan=(0.0, 0.1), mesh = DGMultiMesh(dg, cells_per_dimension, mapping; periodicity=(true, false), is_on_boundary), + @trixi_testset "DGMulti: elixir_navierstokes_convergence_curved.jl" begin + @test_trixi_include(joinpath(examples_dir(), "dgmulti_2d", "elixir_navierstokes_convergence_curved.jl"), + cells_per_dimension = (4, 4), tspan=(0.0, 0.1), l2 = [0.004255101916146187, 0.011118488923215765, 0.011281831283462686, 0.03573656447388509], linf = [0.015071710669706473, 0.04103132025858458, 0.03990424085750277, 0.1309401718598764], ) diff --git a/test/test_parabolic_3d.jl b/test/test_parabolic_3d.jl index 210e571e8d4..1ae5eed44ae 100644 --- a/test/test_parabolic_3d.jl +++ b/test/test_parabolic_3d.jl @@ -19,15 +19,9 @@ isdir(outdir) && rm(outdir, recursive=true) ) end - @trixi_testset "DGMulti: elixir_navierstokes_convergence.jl (curved)" begin - function mapping(xi, eta, zeta) - x = xi + 0.1 * sin(pi * xi) * sin(pi * eta) - y = eta + 0.1 * sin(pi * xi) * sin(pi * eta) - z = zeta + 0.1 * sin(pi * xi) * sin(pi * eta) - return SVector(x, y, z) - end - @test_trixi_include(joinpath(examples_dir(), "dgmulti_3d", "elixir_navierstokes_convergence.jl"), - cells_per_dimension = (4, 4, 4), tspan=(0.0, 0.1), mesh = DGMultiMesh(dg, cells_per_dimension, mapping; periodicity=(true, false, true), is_on_boundary), + @trixi_testset "DGMulti: elixir_navierstokes_convergence_curved.jl" begin + @test_trixi_include(joinpath(examples_dir(), "dgmulti_3d", "elixir_navierstokes_convergence_curved.jl"), + cells_per_dimension = (4, 4, 4), tspan=(0.0, 0.1), l2 = [0.0014027227251207474, 0.0021322235533273513, 0.0027873741447455194, 0.0024587473070627423, 0.00997836818019202], linf = [0.006341750402837576, 0.010306014252246865, 0.01520740250924979, 0.010968264045485565, 0.047454389831591115] ) From 5ce9eef5d7fed7d51260bc349bfbe10c7faea037 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 27 Apr 2023 22:51:38 -0500 Subject: [PATCH 35/36] add comments --- examples/dgmulti_2d/elixir_navierstokes_convergence_curved.jl | 2 ++ examples/dgmulti_3d/elixir_navierstokes_convergence_curved.jl | 2 ++ 2 files changed, 4 insertions(+) diff --git a/examples/dgmulti_2d/elixir_navierstokes_convergence_curved.jl b/examples/dgmulti_2d/elixir_navierstokes_convergence_curved.jl index 854c2af773d..86b5ae64348 100644 --- a/examples/dgmulti_2d/elixir_navierstokes_convergence_curved.jl +++ b/examples/dgmulti_2d/elixir_navierstokes_convergence_curved.jl @@ -29,6 +29,8 @@ end cells_per_dimension = (16, 16) mesh = DGMultiMesh(dg, cells_per_dimension, mapping; periodicity=(true, false), is_on_boundary) +# This initial condition is taken from `examples/dgmulti_2d/elixir_navierstokes_convergence.jl` + # Note: the initial condition cannot be specialized to `CompressibleNavierStokesDiffusion2D` # since it is called by both the parabolic solver (which passes in `CompressibleNavierStokesDiffusion2D`) # and by the initial condition (which passes in `CompressibleEulerEquations2D`). diff --git a/examples/dgmulti_3d/elixir_navierstokes_convergence_curved.jl b/examples/dgmulti_3d/elixir_navierstokes_convergence_curved.jl index b8971a645e2..c14d6620803 100644 --- a/examples/dgmulti_3d/elixir_navierstokes_convergence_curved.jl +++ b/examples/dgmulti_3d/elixir_navierstokes_convergence_curved.jl @@ -28,6 +28,8 @@ end cells_per_dimension = (8, 8, 8) mesh = DGMultiMesh(dg, cells_per_dimension, mapping; periodicity=(true, false, true), is_on_boundary) +# This initial condition is taken from `examples/dgmulti_3d/elixir_navierstokes_convergence.jl` + # Note: the initial condition cannot be specialized to `CompressibleNavierStokesDiffusion3D` # since it is called by both the parabolic solver (which passes in `CompressibleNavierStokesDiffusion3D`) # and by the initial condition (which passes in `CompressibleEulerEquations3D`). From a8c0933bd0077c45372c3937c40fce772e348886 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 28 Apr 2023 06:41:48 -0500 Subject: [PATCH 36/36] remove inv_h from penalty formula --- src/equations/laplace_diffusion_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equations/laplace_diffusion_2d.jl b/src/equations/laplace_diffusion_2d.jl index 2f1afe25a6d..3963c616af2 100644 --- a/src/equations/laplace_diffusion_2d.jl +++ b/src/equations/laplace_diffusion_2d.jl @@ -29,7 +29,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) - return dg.penalty_parameter * (u_outer - u_inner) * equations_parabolic.diffusivity * inv_h + return dg.penalty_parameter * (u_outer - u_inner) * equations_parabolic.diffusivity end # Dirichlet-type boundary condition for use with a parabolic solver in weak form