From f956ac1fba335a1e01b2a2978b66c438c43cc1e9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Thu, 13 Jul 2023 09:14:40 +0200 Subject: [PATCH 01/19] Started to implement mimetic metrics --- ...r_advection_free_stream_mimetic_metrics.jl | 75 +++++++++++++++++++ src/solvers/dgsem_structured/containers_3d.jl | 25 +++++++ 2 files changed, 100 insertions(+) create mode 100644 examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl diff --git a/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl b/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl new file mode 100644 index 00000000000..399edb74349 --- /dev/null +++ b/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl @@ -0,0 +1,75 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection equation + +advection_velocity = (0.2, -0.7, 0.5) +equations = LinearScalarAdvectionEquation3D(advection_velocity) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +solver = DGSEM(20, flux_lax_friedrichs) + +# Mapping as described in https://arxiv.org/abs/2012.12040 +function mapping(xi, eta, zeta) + # Transform input variables between -1 and 1 onto [0,3] + + y = eta + 0.1 * (cos(pi * xi) * + cos(pi * eta) * + cos(pi * zeta)) + + x = xi + 0.1 * (cos(pi * xi) * + cos(pi * eta) * + cos(pi * zeta)) + + z = zeta + 0.1 * (cos(pi * xi) * + cos(pi * eta) * + cos(pi * zeta)) + + return SVector(x, y, z) +end + +cells_per_dimension = (1, 1, 1) + +# Create curved mesh with 8 x 8 x 8 elements +mesh = StructuredMesh(cells_per_dimension, mapping) + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_constant, solver) + + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to 1.0 +ode = semidiscretize(semi, (0.0, 1.0)); + +# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup +# and resets the timers +summary_callback = SummaryCallback() + +# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results +analysis_callback = AnalysisCallback(semi, interval=100) + +# The SaveSolutionCallback allows to save the solution to a file in regular intervals +save_solution = SaveSolutionCallback(interval=100, + solution_variables=cons2prim) + +# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step +stepsize_callback = StepsizeCallback(cfl=2.0) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback) + + +############################################################################### +# run the simulation + +# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks +sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), + dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep=false, callback=callbacks); + +# Print the timer summary +summary_callback() diff --git a/src/solvers/dgsem_structured/containers_3d.jl b/src/solvers/dgsem_structured/containers_3d.jl index e843e869bf5..0fc85958602 100644 --- a/src/solvers/dgsem_structured/containers_3d.jl +++ b/src/solvers/dgsem_structured/containers_3d.jl @@ -242,6 +242,31 @@ function calc_contravariant_vectors!(contravariant_vectors::AbstractArray{<:Any, return contravariant_vectors end +""" +New function to compute contravariant vectors +""" +function calc_contravariant_vectors!(contravariant_vectors::AbstractArray{<:Any, 6}, + element, + jacobian_matrix, node_coordinates, + basis::LobattoLegendreBasis) + @unpack derivative_matrix = basis + + # Edge basis functions + V = zeros(Float64, polydeg(basis) + 1, polydeg(basis)) + for j in 1:polydeg(basis) + for i in 1:polydeg(basis)+1 + for k in 1:j + V[i, j] -= derivative_matrix[i, k] + end + end + end +end + +theta(xi, eta, zeta) = 0.1 * cos(pi * xi) * cos(pi * eta) * cos(pi * zeta) +theta_int1(xi, eta, zeta) = (0.1 / pi) * sin(pi * xi) * cos(pi * eta) * cos(pi * zeta) +theta_int2(xi, eta, zeta) = (0.1 / pi) * cos(pi * xi) * sin(pi * eta) * cos(pi * zeta) +theta_int3(xi, eta, zeta) = (0.1 / pi) * cos(pi * xi) * cos(pi * eta) * sin(pi * zeta) + # Calculate inverse Jacobian (determinant of Jacobian matrix of the mapping) in each node function calc_inverse_jacobian!(inverse_jacobian::AbstractArray{<:Any, 4}, element, jacobian_matrix, basis) From 6273c2a1be74212e71b9563f0501dc0666d8a302 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Thu, 13 Jul 2023 10:32:55 +0200 Subject: [PATCH 02/19] Added most of the mimetic metrics computations (only analytic projection is needed) --- src/solvers/dgsem_structured/containers_3d.jl | 108 +++++++++++++++++- 1 file changed, 105 insertions(+), 3 deletions(-) diff --git a/src/solvers/dgsem_structured/containers_3d.jl b/src/solvers/dgsem_structured/containers_3d.jl index 0fc85958602..106e1cb7db2 100644 --- a/src/solvers/dgsem_structured/containers_3d.jl +++ b/src/solvers/dgsem_structured/containers_3d.jl @@ -245,14 +245,14 @@ end """ New function to compute contravariant vectors """ -function calc_contravariant_vectors!(contravariant_vectors::AbstractArray{<:Any, 6}, +function calc_contravariant_vectors_mimetic!(contravariant_vectors::AbstractArray{<:Any, 6}, element, jacobian_matrix, node_coordinates, basis::LobattoLegendreBasis) @unpack derivative_matrix = basis - # Edge basis functions - V = zeros(Float64, polydeg(basis) + 1, polydeg(basis)) + # Define histopolation (edge) basis functions: V[i,j] = hⱼ(ξᵢ) ... TODO: initialize beforehand... + V = zero(MMatrix{polydeg(basis) + 1, polydeg(basis), eltype(derivative_matrix)}) for j in 1:polydeg(basis) for i in 1:polydeg(basis)+1 for k in 1:j @@ -260,6 +260,108 @@ function calc_contravariant_vectors!(contravariant_vectors::AbstractArray{<:Any, end end end + + # Project the mapping potential \vec{g} \in H_{curl} to \vec{G} \in V_1 + Gbar = zeros(3, 3, nnodes(basis), nnodes(basis), nnodes(basis), eltype(derivative_matrix)) # Attention: here I'm allocating N+1 nodes in each direction. We only need N in some directions!! + # ( here we compute Gbar[:, :, i, j, k] ) + + # Evaluate the mapping potential at the Lagrange points + G = zeros(3, 3, nnodes(basis), nnodes(basis), nnodes(basis), eltype(derivative_matrix)) + for k in eachnode(basis) + for j in eachnode(basis) + for i in eachnode(basis) + for ii in 1:polydeg(basis) + G[:, 1, i, j, k] = Gbar[:, 1, ii, j, k] * V[i, ii] + G[:, 2, i, j, k] = Gbar[:, 2, i, ii, k] * V[j, ii] + G[:, 3, i, j, k] = Gbar[:, 3, i, j, ii] * V[k, ii] + end + end + end + end + + # Compute the contravariant vectors as the curl of the mapping potential (at the discrete level!) + # Jaⁱₙ = ( ∇ × gₙ )ᵢ where ∇ = (∂/∂ξ, ∂/∂η, ∂/∂ζ)ᵀ + for n in 1:3 + # Calculate Ja¹ₙ = (g³ₙ)_η - (g²ₙ)_ζ + # For each of these, the first and second summand are computed in separate loops + # for performance reasons. + + # First summand (g³ₙ)_η + @turbo for k in eachnode(basis), j in eachnode(basis), i in eachnode(basis) + result = zero(eltype(contravariant_vectors)) + + for ii in eachnode(basis) + # Multiply derivative_matrix to j-dimension to differentiate wrt η + result += derivative_matrix[j, ii] * G[n, 1, i, ii, k] + end + + contravariant_vectors[n, 1, i, j, k, element] = result + end + + # Second summand -(g²ₙ)_ζ + @turbo for k in eachnode(basis), j in eachnode(basis), i in eachnode(basis) + result = zero(eltype(contravariant_vectors)) + + for ii in eachnode(basis) + # Multiply derivative_matrix to k-dimension to differentiate wrt ζ + result += derivative_matrix[k, ii] * G[n, 2, i, j, ii] + end + + contravariant_vectors[n, 1, i, j, k, element] -= result + end + + # Calculate Ja²ₙ =(g¹ₙ)_ζ - (g³ₙ)_ξ + + # First summand (g¹ₙ)_ζ + @turbo for k in eachnode(basis), j in eachnode(basis), i in eachnode(basis) + result = zero(eltype(contravariant_vectors)) + + for ii in eachnode(basis) + # Multiply derivative_matrix to k-dimension to differentiate wrt ζ + result += derivative_matrix[k, ii] * G[n, 1, i, j, ii] + end + + contravariant_vectors[n, 2, i, j, k, element] = result + end + + # Second summand -(g³ₙ)_ξ + @turbo for k in eachnode(basis), j in eachnode(basis), i in eachnode(basis) + result = zero(eltype(contravariant_vectors)) + + for ii in eachnode(basis) + # Multiply derivative_matrix to i-dimension to differentiate wrt ξ + result += derivative_matrix[i, ii] * G[n, 3, ii, j, k] + end + + contravariant_vectors[n, 2, i, j, k, element] -= result + end + + # Calculate Ja³ₙ = (g²ₙ)_ξ - (g¹ₙ)_η + + # First summand (g²ₙ)_ξ + @turbo for k in eachnode(basis), j in eachnode(basis), i in eachnode(basis) + result = zero(eltype(contravariant_vectors)) + + for ii in eachnode(basis) + # Multiply derivative_matrix to i-dimension to differentiate wrt ξ + result += derivative_matrix[i, ii] * G[n, 2, ii, j, k] + end + + contravariant_vectors[n, 3, i, j, k, element] = result + end + + # Second summand -(g¹ₙ)_η + @turbo for k in eachnode(basis), j in eachnode(basis), i in eachnode(basis) + result = zero(eltype(contravariant_vectors)) + + for ii in eachnode(basis) + # Multiply derivative_matrix to j-dimension to differentiate wrt η + result += derivative_matrix[j, ii] * G[n, 1, i, ii, k] + end + + contravariant_vectors[n, 3, i, j, k, element] -= result + end + end end theta(xi, eta, zeta) = 0.1 * cos(pi * xi) * cos(pi * eta) * cos(pi * zeta) From 31e226af3afce378beaf606b4841a96e5dbe5ccc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Thu, 13 Jul 2023 11:18:02 +0200 Subject: [PATCH 03/19] Fixed some bugs --- .../elixir_advection_free_stream_mimetic_metrics.jl | 2 +- src/solvers/dgsem_structured/containers_3d.jl | 12 ++++++++++-- 2 files changed, 11 insertions(+), 3 deletions(-) diff --git a/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl b/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl index 399edb74349..0c83721918b 100644 --- a/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl +++ b/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl @@ -13,7 +13,7 @@ solver = DGSEM(20, flux_lax_friedrichs) # Mapping as described in https://arxiv.org/abs/2012.12040 function mapping(xi, eta, zeta) - # Transform input variables between -1 and 1 onto [0,3] + # Transform input variables between -1 and 1 onto our crazy domain y = eta + 0.1 * (cos(pi * xi) * cos(pi * eta) * diff --git a/src/solvers/dgsem_structured/containers_3d.jl b/src/solvers/dgsem_structured/containers_3d.jl index 106e1cb7db2..6358a1003bd 100644 --- a/src/solvers/dgsem_structured/containers_3d.jl +++ b/src/solvers/dgsem_structured/containers_3d.jl @@ -263,7 +263,13 @@ function calc_contravariant_vectors_mimetic!(contravariant_vectors::AbstractArra # Project the mapping potential \vec{g} \in H_{curl} to \vec{G} \in V_1 Gbar = zeros(3, 3, nnodes(basis), nnodes(basis), nnodes(basis), eltype(derivative_matrix)) # Attention: here I'm allocating N+1 nodes in each direction. We only need N in some directions!! - # ( here we compute Gbar[:, :, i, j, k] ) + for k in eachnode(basis) + for j in eachnode(basis) + for i in eachnode(basis) + #G[1, 1, ] + end + end + end # Evaluate the mapping potential at the Lagrange points G = zeros(3, 3, nnodes(basis), nnodes(basis), nnodes(basis), eltype(derivative_matrix)) @@ -292,7 +298,7 @@ function calc_contravariant_vectors_mimetic!(contravariant_vectors::AbstractArra for ii in eachnode(basis) # Multiply derivative_matrix to j-dimension to differentiate wrt η - result += derivative_matrix[j, ii] * G[n, 1, i, ii, k] + result += derivative_matrix[j, ii] * G[n, 3, i, ii, k] end contravariant_vectors[n, 1, i, j, k, element] = result @@ -362,6 +368,8 @@ function calc_contravariant_vectors_mimetic!(contravariant_vectors::AbstractArra contravariant_vectors[n, 3, i, j, k, element] -= result end end + + return contravariant_vectors end theta(xi, eta, zeta) = 0.1 * cos(pi * xi) * cos(pi * eta) * cos(pi * zeta) From 873762a6a36512f55ebb769b9614dd7586c894ba Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Thu, 13 Jul 2023 11:56:29 +0200 Subject: [PATCH 04/19] Added analytical projection of metric potential (for a particular example) --- src/solvers/dgsem_structured/containers_3d.jl | 42 ++++++++++++++++++- 1 file changed, 40 insertions(+), 2 deletions(-) diff --git a/src/solvers/dgsem_structured/containers_3d.jl b/src/solvers/dgsem_structured/containers_3d.jl index 6358a1003bd..b5a5afce16a 100644 --- a/src/solvers/dgsem_structured/containers_3d.jl +++ b/src/solvers/dgsem_structured/containers_3d.jl @@ -249,7 +249,7 @@ function calc_contravariant_vectors_mimetic!(contravariant_vectors::AbstractArra element, jacobian_matrix, node_coordinates, basis::LobattoLegendreBasis) - @unpack derivative_matrix = basis + @unpack derivative_matrix, nodes = basis # Define histopolation (edge) basis functions: V[i,j] = hⱼ(ξᵢ) ... TODO: initialize beforehand... V = zero(MMatrix{polydeg(basis) + 1, polydeg(basis), eltype(derivative_matrix)}) @@ -264,9 +264,47 @@ function calc_contravariant_vectors_mimetic!(contravariant_vectors::AbstractArra # Project the mapping potential \vec{g} \in H_{curl} to \vec{G} \in V_1 Gbar = zeros(3, 3, nnodes(basis), nnodes(basis), nnodes(basis), eltype(derivative_matrix)) # Attention: here I'm allocating N+1 nodes in each direction. We only need N in some directions!! for k in eachnode(basis) + for j in eachnode(basis) + for i in 1:polydeg(basis) + Gbar[1, 1, i, j, k] = (nodes[k] * (theta(nodes[i+1], nodes[j], nodes[k]) - theta(nodes[i], nodes[j], nodes[k]) ) + + 0.5 * (theta(nodes[i+1], nodes[j], nodes[k])^2 - theta(nodes[i], nodes[j], nodes[k])^2 ) ) + Gbar[2, 1, i, j, k] = (nodes[i+1] * theta(nodes[i+1], nodes[j], nodes[k]) - nodes[i] * theta(nodes[i], nodes[j], nodes[k]) + + 0.5 * (theta(nodes[i+1], nodes[j], nodes[k])^2 - theta(nodes[i], nodes[j], nodes[k])^2 ) + - (theta_int1(nodes[i+1], nodes[j], nodes[k]) - theta_int1(nodes[i], nodes[j], nodes[k]))) + Gbar[3, 1, i, j, k] = ( nodes[j] * (nodes[i+1] - nodes[i]) + + nodes[j] * ( theta(nodes[i+1], nodes[j], nodes[k]) - theta(nodes[i], nodes[j], nodes[k]) ) + + 0.5 * (theta(nodes[i+1], nodes[j], nodes[k])^2 - theta(nodes[i], nodes[j], nodes[k])^2 ) + + (theta_int1(nodes[i+1], nodes[j], nodes[k]) - theta_int1(nodes[i], nodes[j], nodes[k]))) + end + end + end + for k in eachnode(basis) + for j in 1:polydeg(basis) + for i in eachnode(basis) + Gbar[1, 2, i, j, k] = ( nodes[k] * (nodes[j+1] - nodes[j]) + + nodes[k] * ( theta(nodes[i], nodes[j+1], nodes[k]) - theta(nodes[i], nodes[j], nodes[k]) ) + + 0.5 * (theta(nodes[i], nodes[j+1], nodes[k])^2 - theta(nodes[i], nodes[j], nodes[k])^2 ) + + (theta_int2(nodes[i], nodes[j+1], nodes[k]) - theta_int2(nodes[i], nodes[j], nodes[k]))) + Gbar[2, 2, i, j, k] = (nodes[i] * (theta(nodes[i], nodes[j+1], nodes[k]) - theta(nodes[i], nodes[j], nodes[k]) ) + + 0.5 * (theta(nodes[i], nodes[j+1], nodes[k])^2 - theta(nodes[i], nodes[j], nodes[k])^2 ) ) + Gbar[3, 2, i, j, k] = (nodes[j+1] * theta(nodes[i], nodes[j+1], nodes[k]) - nodes[j] * theta(nodes[i], nodes[j], nodes[k]) + + 0.5 * (theta(nodes[i], nodes[j+1], nodes[k])^2 - theta(nodes[i], nodes[j], nodes[k])^2 ) + - (theta_int2(nodes[i], nodes[j+1], nodes[k]) - theta_int2(nodes[i], nodes[j], nodes[k]))) + end + end + end + for k in 1:polydeg(basis) for j in eachnode(basis) for i in eachnode(basis) - #G[1, 1, ] + Gbar[1, 3, i, j, k] = (nodes[k+1] * theta(nodes[i], nodes[j], nodes[k+1]) - nodes[k] * theta(nodes[i], nodes[j], nodes[k]) + + 0.5 * (theta(nodes[i], nodes[j], nodes[k+1])^2 - theta(nodes[i], nodes[j], nodes[k])^2 ) + - (theta_int3(nodes[i], nodes[j], nodes[k+1]) - theta_int3(nodes[i], nodes[j], nodes[k]))) + Gbar[2, 3, i, j, k] = ( nodes[i] * (nodes[k+1] - nodes[k]) + + nodes[i] * ( theta(nodes[i], nodes[j], nodes[k+1]) - theta(nodes[i], nodes[j], nodes[k]) ) + + 0.5 * (theta(nodes[i], nodes[j], nodes[k+1])^2 - theta(nodes[i], nodes[j], nodes[k])^2 ) + + (theta_int3(nodes[i], nodes[j], nodes[k+1]) - theta_int3(nodes[i], nodes[j], nodes[k]))) + Gbar[3, 3, i, j, k] = (nodes[j] * (theta(nodes[i], nodes[j], nodes[k+1]) - theta(nodes[i], nodes[j], nodes[k]) ) + + 0.5 * (theta(nodes[i], nodes[j], nodes[k+1])^2 - theta(nodes[i], nodes[j], nodes[k])^2 ) ) end end end From 34e155041e466bbc6b7c7ad9ade8e33df374fbbd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Thu, 13 Jul 2023 12:03:12 +0200 Subject: [PATCH 05/19] Fixed some bugs --- .../elixir_advection_free_stream_mimetic_metrics.jl | 6 +++--- src/Trixi.jl | 2 +- src/solvers/dgsem_structured/containers_3d.jl | 10 ++++++---- 3 files changed, 10 insertions(+), 8 deletions(-) diff --git a/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl b/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl index 0c83721918b..424fe15b7a8 100644 --- a/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl +++ b/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl @@ -9,7 +9,7 @@ advection_velocity = (0.2, -0.7, 0.5) equations = LinearScalarAdvectionEquation3D(advection_velocity) # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux -solver = DGSEM(20, flux_lax_friedrichs) +solver = DGSEM(10, flux_lax_friedrichs) # Mapping as described in https://arxiv.org/abs/2012.12040 function mapping(xi, eta, zeta) @@ -30,7 +30,7 @@ function mapping(xi, eta, zeta) return SVector(x, y, z) end -cells_per_dimension = (1, 1, 1) +cells_per_dimension = (1,1,1) # Create curved mesh with 8 x 8 x 8 elements mesh = StructuredMesh(cells_per_dimension, mapping) @@ -57,7 +57,7 @@ save_solution = SaveSolutionCallback(interval=100, solution_variables=cons2prim) # The StepsizeCallback handles the re-calculation of the maximum Δt after each time step -stepsize_callback = StepsizeCallback(cfl=2.0) +stepsize_callback = StepsizeCallback(cfl=0.1) # Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback) diff --git a/src/Trixi.jl b/src/Trixi.jl index 66878f4b459..f140d5758f6 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -59,7 +59,7 @@ using RecipesBase: RecipesBase using Requires: @require using Static: Static, One, True, False @reexport using StaticArrays: SVector -using StaticArrays: StaticArrays, MVector, MArray, SMatrix, @SMatrix +using StaticArrays: StaticArrays, MVector, MArray, SMatrix, @SMatrix, MMatrix using StrideArrays: PtrArray, StrideArray, StaticInt @reexport using StructArrays: StructArrays, StructArray using TimerOutputs: TimerOutputs, @notimeit, TimerOutput, print_timer, reset_timer! diff --git a/src/solvers/dgsem_structured/containers_3d.jl b/src/solvers/dgsem_structured/containers_3d.jl index b5a5afce16a..da9b74d0e23 100644 --- a/src/solvers/dgsem_structured/containers_3d.jl +++ b/src/solvers/dgsem_structured/containers_3d.jl @@ -120,7 +120,7 @@ end # Calculate contravariant vectors, multiplied by the Jacobian determinant J of the transformation mapping, # using the invariant curl form. # These are called Ja^i in Kopriva's blue book. -function calc_contravariant_vectors!(contravariant_vectors::AbstractArray{<:Any, 6}, +function calc_contravariant_vectors_standard!(contravariant_vectors::AbstractArray{<:Any, 6}, element, jacobian_matrix, node_coordinates, basis::LobattoLegendreBasis) @@ -245,12 +245,14 @@ end """ New function to compute contravariant vectors """ -function calc_contravariant_vectors_mimetic!(contravariant_vectors::AbstractArray{<:Any, 6}, +function calc_contravariant_vectors!(contravariant_vectors::AbstractArray{<:Any, 6}, element, jacobian_matrix, node_coordinates, basis::LobattoLegendreBasis) @unpack derivative_matrix, nodes = basis + + println("Hello world ;)") # Define histopolation (edge) basis functions: V[i,j] = hⱼ(ξᵢ) ... TODO: initialize beforehand... V = zero(MMatrix{polydeg(basis) + 1, polydeg(basis), eltype(derivative_matrix)}) for j in 1:polydeg(basis) @@ -262,7 +264,7 @@ function calc_contravariant_vectors_mimetic!(contravariant_vectors::AbstractArra end # Project the mapping potential \vec{g} \in H_{curl} to \vec{G} \in V_1 - Gbar = zeros(3, 3, nnodes(basis), nnodes(basis), nnodes(basis), eltype(derivative_matrix)) # Attention: here I'm allocating N+1 nodes in each direction. We only need N in some directions!! + Gbar = zeros(eltype(derivative_matrix), 3, 3, nnodes(basis), nnodes(basis), nnodes(basis)) # Attention: here I'm allocating N+1 nodes in each direction. We only need N in some directions!! for k in eachnode(basis) for j in eachnode(basis) for i in 1:polydeg(basis) @@ -310,7 +312,7 @@ function calc_contravariant_vectors_mimetic!(contravariant_vectors::AbstractArra end # Evaluate the mapping potential at the Lagrange points - G = zeros(3, 3, nnodes(basis), nnodes(basis), nnodes(basis), eltype(derivative_matrix)) + G = zeros(eltype(derivative_matrix), 3, 3, nnodes(basis), nnodes(basis), nnodes(basis)) for k in eachnode(basis) for j in eachnode(basis) for i in eachnode(basis) From 2955482ff11daf84bc4d800192dbc8a5ff0b2e47 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Thu, 13 Jul 2023 14:35:30 +0200 Subject: [PATCH 06/19] Fixed bugs. FIRST WORKING VERSION. --- src/solvers/dgsem_structured/containers_3d.jl | 26 +++++++++---------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/solvers/dgsem_structured/containers_3d.jl b/src/solvers/dgsem_structured/containers_3d.jl index da9b74d0e23..df63dd9d352 100644 --- a/src/solvers/dgsem_structured/containers_3d.jl +++ b/src/solvers/dgsem_structured/containers_3d.jl @@ -317,18 +317,18 @@ function calc_contravariant_vectors!(contravariant_vectors::AbstractArray{<:Any, for j in eachnode(basis) for i in eachnode(basis) for ii in 1:polydeg(basis) - G[:, 1, i, j, k] = Gbar[:, 1, ii, j, k] * V[i, ii] - G[:, 2, i, j, k] = Gbar[:, 2, i, ii, k] * V[j, ii] - G[:, 3, i, j, k] = Gbar[:, 3, i, j, ii] * V[k, ii] + G[:, 1, i, j, k] += Gbar[:, 1, ii, j, k] * V[i, ii] + G[:, 2, i, j, k] += Gbar[:, 2, i, ii, k] * V[j, ii] + G[:, 3, i, j, k] += Gbar[:, 3, i, j, ii] * V[k, ii] end end end end # Compute the contravariant vectors as the curl of the mapping potential (at the discrete level!) - # Jaⁱₙ = ( ∇ × gₙ )ᵢ where ∇ = (∂/∂ξ, ∂/∂η, ∂/∂ζ)ᵀ + # Jaⁱₙ = -( ∇ × gₙ )ᵢ where ∇ = (∂/∂ξ, ∂/∂η, ∂/∂ζ)ᵀ for n in 1:3 - # Calculate Ja¹ₙ = (g³ₙ)_η - (g²ₙ)_ζ + # Calculate Ja¹ₙ = -(g³ₙ)_η + (g²ₙ)_ζ # For each of these, the first and second summand are computed in separate loops # for performance reasons. @@ -341,7 +341,7 @@ function calc_contravariant_vectors!(contravariant_vectors::AbstractArray{<:Any, result += derivative_matrix[j, ii] * G[n, 3, i, ii, k] end - contravariant_vectors[n, 1, i, j, k, element] = result + contravariant_vectors[n, 1, i, j, k, element] = -result end # Second summand -(g²ₙ)_ζ @@ -353,10 +353,10 @@ function calc_contravariant_vectors!(contravariant_vectors::AbstractArray{<:Any, result += derivative_matrix[k, ii] * G[n, 2, i, j, ii] end - contravariant_vectors[n, 1, i, j, k, element] -= result + contravariant_vectors[n, 1, i, j, k, element] += result end - # Calculate Ja²ₙ =(g¹ₙ)_ζ - (g³ₙ)_ξ + # Calculate Ja²ₙ = -(g¹ₙ)_ζ + (g³ₙ)_ξ # First summand (g¹ₙ)_ζ @turbo for k in eachnode(basis), j in eachnode(basis), i in eachnode(basis) @@ -367,7 +367,7 @@ function calc_contravariant_vectors!(contravariant_vectors::AbstractArray{<:Any, result += derivative_matrix[k, ii] * G[n, 1, i, j, ii] end - contravariant_vectors[n, 2, i, j, k, element] = result + contravariant_vectors[n, 2, i, j, k, element] = -result end # Second summand -(g³ₙ)_ξ @@ -379,10 +379,10 @@ function calc_contravariant_vectors!(contravariant_vectors::AbstractArray{<:Any, result += derivative_matrix[i, ii] * G[n, 3, ii, j, k] end - contravariant_vectors[n, 2, i, j, k, element] -= result + contravariant_vectors[n, 2, i, j, k, element] += result end - # Calculate Ja³ₙ = (g²ₙ)_ξ - (g¹ₙ)_η + # Calculate Ja³ₙ = -(g²ₙ)_ξ + (g¹ₙ)_η # First summand (g²ₙ)_ξ @turbo for k in eachnode(basis), j in eachnode(basis), i in eachnode(basis) @@ -393,7 +393,7 @@ function calc_contravariant_vectors!(contravariant_vectors::AbstractArray{<:Any, result += derivative_matrix[i, ii] * G[n, 2, ii, j, k] end - contravariant_vectors[n, 3, i, j, k, element] = result + contravariant_vectors[n, 3, i, j, k, element] = -result end # Second summand -(g¹ₙ)_η @@ -405,7 +405,7 @@ function calc_contravariant_vectors!(contravariant_vectors::AbstractArray{<:Any, result += derivative_matrix[j, ii] * G[n, 1, i, ii, k] end - contravariant_vectors[n, 3, i, j, k, element] -= result + contravariant_vectors[n, 3, i, j, k, element] += result end end From 134ff6914156e6ab8039ac6bfe26ec9d29757573 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Thu, 13 Jul 2023 17:13:12 +0200 Subject: [PATCH 07/19] Added scripts to test the metric terms --- ...r_advection_free_stream_mimetic_metrics.jl | 81 ++++++++++++++++--- src/meshes/structured_mesh.jl | 4 +- src/solvers/dgsem_structured/containers_3d.jl | 13 +-- 3 files changed, 83 insertions(+), 15 deletions(-) diff --git a/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl b/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl index 424fe15b7a8..e28a853aa1b 100644 --- a/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl +++ b/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl @@ -1,6 +1,8 @@ - +using LinearAlgebra using OrdinaryDiffEq using Trixi +using StaticArrays +using Plots ############################################################################### # semidiscretization of the linear advection equation @@ -8,9 +10,6 @@ using Trixi advection_velocity = (0.2, -0.7, 0.5) equations = LinearScalarAdvectionEquation3D(advection_velocity) -# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux -solver = DGSEM(10, flux_lax_friedrichs) - # Mapping as described in https://arxiv.org/abs/2012.12040 function mapping(xi, eta, zeta) # Transform input variables between -1 and 1 onto our crazy domain @@ -27,18 +26,82 @@ function mapping(xi, eta, zeta) cos(pi * eta) * cos(pi * zeta)) + #= x = xi + y = eta + z = zeta =# return SVector(x, y, z) end +theta_der1(xi, eta, zeta) = -(0.1 * pi) * sin(pi * xi) * cos(pi * eta) * cos(pi * zeta) +theta_der2(xi, eta, zeta) = -(0.1 * pi) * cos(pi * xi) * sin(pi * eta) * cos(pi * zeta) +theta_der3(xi, eta, zeta) = -(0.1 * pi) * cos(pi * xi) * cos(pi * eta) * sin(pi * zeta) + +function exact_contravariant_vectors!(Ja, xi, eta, zeta) + theta_xi = theta_der1(xi, eta, zeta) + theta_eta = theta_der2(xi, eta, zeta) + theta_zeta = theta_der3(xi, eta, zeta) + Ja[1,1] = 1 + theta_eta + theta_zeta + Ja[1,2] = -theta_xi + Ja[1,3] = -theta_xi + Ja[2,1] = -theta_eta + Ja[2,2] = 1 + theta_xi + theta_zeta + Ja[2,3] = -theta_eta + Ja[3,1] = -theta_zeta + Ja[3,2] = -theta_zeta + Ja[3,3] = 1 + theta_xi + theta_eta +end + +function compute_error(solver, semi) + @unpack nodes, weights = solver.basis + exact_Ja = zero(MMatrix{3, 3, Float64}) + error = zero(Float64) + error_L2 = zero(Float64) + for k in eachnode(solver.basis) + for j in eachnode(solver.basis) + for i in eachnode(solver.basis) + exact_contravariant_vectors!(exact_Ja, nodes[i], nodes[j], nodes[k]) + error = max(error, maximum(abs.(semi.cache.elements.contravariant_vectors[:,:,i,j,k,1] - exact_Ja))) + error_L2 += norm(semi.cache.elements.contravariant_vectors[:,1,i,j,k,1] - exact_Ja[:,1]) * weights[i] * weights[j] * weights[k] + end + end + end + return error, error_L2 / 8 +end + cells_per_dimension = (1,1,1) -# Create curved mesh with 8 x 8 x 8 elements -mesh = StructuredMesh(cells_per_dimension, mapping) +errors_inf = zeros(20,2) +errors_L2 = zeros(20,2) +for polydeg in 1:20 + println("Computing polydeg = ", polydeg) + # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux + solver = DGSEM(polydeg, flux_lax_friedrichs) -# A semidiscretization collects data structures and functions for the spatial discretization -semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_constant, solver) + # Create curved mesh with 8 x 8 x 8 elements + mesh = StructuredMesh(cells_per_dimension, mapping; mimetic = false) + # A semidiscretization collects data structures and functions for the spatial discretization + semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_constant, solver) + + error_inf, error_L2 = compute_error(solver, semi) + errors_inf[polydeg,1] = error_inf + errors_L2[polydeg,1] = error_L2 + + # Create curved mesh with 8 x 8 x 8 elements + mesh = StructuredMesh(cells_per_dimension, mapping; mimetic = true) + + # A semidiscretization collects data structures and functions for the spatial discretization + semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_constant, solver) + + error_inf, error_L2 = compute_error(solver, semi) + errors_inf[polydeg,2] = error_inf + errors_L2[polydeg,2] = error_L2 + +end +plot(errors_L2[:,1], yaxis=:log, title = "standard") +plot!(errors_L2[:,2], yaxis=:log, title = "mimetic") +#= ############################################################################### # ODE solvers, callbacks etc. @@ -72,4 +135,4 @@ sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), save_everystep=false, callback=callbacks); # Print the timer summary -summary_callback() +summary_callback() =# \ No newline at end of file diff --git a/src/meshes/structured_mesh.jl b/src/meshes/structured_mesh.jl index df067db833d..55e696b4b21 100644 --- a/src/meshes/structured_mesh.jl +++ b/src/meshes/structured_mesh.jl @@ -20,6 +20,7 @@ mutable struct StructuredMesh{NDIMS, RealT <: Real} <: AbstractMesh{NDIMS} periodicity::NTuple{NDIMS, Bool} current_filename::String unsaved_changes::Bool + mimetic::Bool end """ @@ -44,6 +45,7 @@ Create a StructuredMesh of the given size and shape that uses `RealT` as coordin This will be changed in the future, see [https://github.com/trixi-framework/Trixi.jl/issues/541](https://github.com/trixi-framework/Trixi.jl/issues/541). """ function StructuredMesh(cells_per_dimension, mapping; RealT = Float64, + mimetic = false, periodicity = true, unsaved_changes = true, mapping_as_string = mapping2string(mapping, length(cells_per_dimension))) @@ -63,7 +65,7 @@ function StructuredMesh(cells_per_dimension, mapping; RealT = Float64, return StructuredMesh{NDIMS, RealT}(Tuple(cells_per_dimension), mapping, mapping_as_string, periodicity, "", - unsaved_changes) + unsaved_changes, mimetic) end """ diff --git a/src/solvers/dgsem_structured/containers_3d.jl b/src/solvers/dgsem_structured/containers_3d.jl index df63dd9d352..b1b8c3e225d 100644 --- a/src/solvers/dgsem_structured/containers_3d.jl +++ b/src/solvers/dgsem_structured/containers_3d.jl @@ -21,8 +21,13 @@ function init_elements!(elements, mesh::StructuredMesh{3}, basis::LobattoLegendr calc_jacobian_matrix!(jacobian_matrix, element, node_coordinates, basis) - calc_contravariant_vectors!(contravariant_vectors, element, jacobian_matrix, - node_coordinates, basis) + if mesh.mimetic + calc_contravariant_vectors_mimetic!(contravariant_vectors, element, jacobian_matrix, + node_coordinates, basis) + else + calc_contravariant_vectors_standard!(contravariant_vectors, element, jacobian_matrix, + node_coordinates, basis) + end calc_inverse_jacobian!(inverse_jacobian, element, jacobian_matrix, basis) end @@ -245,14 +250,12 @@ end """ New function to compute contravariant vectors """ -function calc_contravariant_vectors!(contravariant_vectors::AbstractArray{<:Any, 6}, +function calc_contravariant_vectors_mimetic!(contravariant_vectors::AbstractArray{<:Any, 6}, element, jacobian_matrix, node_coordinates, basis::LobattoLegendreBasis) @unpack derivative_matrix, nodes = basis - - println("Hello world ;)") # Define histopolation (edge) basis functions: V[i,j] = hⱼ(ξᵢ) ... TODO: initialize beforehand... V = zero(MMatrix{polydeg(basis) + 1, polydeg(basis), eltype(derivative_matrix)}) for j in 1:polydeg(basis) From ba66691e6b5b4139668d0afe6ba131cbbf382338 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Fri, 14 Jul 2023 13:23:46 +0200 Subject: [PATCH 08/19] Added computation of exact jacobian matrix (potentential) before projection for standard metric terms... And test with periodic and Dirichlent BCs --- ...r_advection_free_stream_mimetic_metrics.jl | 93 +++++++-- ...n_free_stream_mimetic_metrics_dirichlet.jl | 195 ++++++++++++++++++ src/Trixi.jl | 2 + src/meshes/structured_mesh.jl | 4 +- src/solvers/dgsem_structured/containers_3d.jl | 55 ++++- 5 files changed, 328 insertions(+), 21 deletions(-) create mode 100644 examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics_dirichlet.jl diff --git a/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl b/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl index e28a853aa1b..528eb4c6f90 100644 --- a/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl +++ b/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl @@ -32,10 +32,6 @@ function mapping(xi, eta, zeta) return SVector(x, y, z) end -theta_der1(xi, eta, zeta) = -(0.1 * pi) * sin(pi * xi) * cos(pi * eta) * cos(pi * zeta) -theta_der2(xi, eta, zeta) = -(0.1 * pi) * cos(pi * xi) * sin(pi * eta) * cos(pi * zeta) -theta_der3(xi, eta, zeta) = -(0.1 * pi) * cos(pi * xi) * cos(pi * eta) * sin(pi * zeta) - function exact_contravariant_vectors!(Ja, xi, eta, zeta) theta_xi = theta_der1(xi, eta, zeta) theta_eta = theta_der2(xi, eta, zeta) @@ -70,37 +66,98 @@ end cells_per_dimension = (1,1,1) -errors_inf = zeros(20,2) -errors_L2 = zeros(20,2) -for polydeg in 1:20 +max_polydeg = 25 + +errors_normals_inf = zeros(max_polydeg,2) +errors_normals_L2 = zeros(max_polydeg,2) +errors_sol_inf = zeros(max_polydeg,2) +errors_sol_L2 = zeros(max_polydeg,2) +exact_jacobian = true +final_time = 3.00636132e-03 +initial_condition = initial_condition_constant +for polydeg in 1:max_polydeg println("Computing polydeg = ", polydeg) # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux solver = DGSEM(polydeg, flux_lax_friedrichs) # Create curved mesh with 8 x 8 x 8 elements - mesh = StructuredMesh(cells_per_dimension, mapping; mimetic = false) + mesh = StructuredMesh(cells_per_dimension, mapping; mimetic = false, exact_jacobian = exact_jacobian) - # A semidiscretization collects data structures and functions for the spatial discretization - semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_constant, solver) + # A semidiscre tization collects data structures and functions for the spatial discretization + semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) error_inf, error_L2 = compute_error(solver, semi) - errors_inf[polydeg,1] = error_inf - errors_L2[polydeg,1] = error_L2 + errors_normals_inf[polydeg,1] = error_inf + errors_normals_L2[polydeg,1] = error_L2 + + # Create ODE problem with time span from 0.0 to 1.0 + ode = semidiscretize(semi, (0.0, final_time)); + + # The AnalysisCallback allows to analyse the solution in regular intervals and prints the results + analysis_callback = AnalysisCallback(semi, interval=100) + + # The StepsizeCallback handles the re-calculation of the maximum Δt after each time step + stepsize_callback = StepsizeCallback(cfl=0.1) + + # Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver + callbacks = CallbackSet(analysis_callback, stepsize_callback) + + + ############################################################################### + # run the simulation + + # OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks + sol = solve(ode, Euler(), #CarpenterKennedy2N54(williamson_condition=false), + dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep=false, callback=callbacks); + + errors = analysis_callback(sol) + + errors_sol_L2[polydeg, 1] = errors.l2[1] + errors_sol_inf[polydeg, 1] = errors.linf[1] # Create curved mesh with 8 x 8 x 8 elements - mesh = StructuredMesh(cells_per_dimension, mapping; mimetic = true) + mesh = StructuredMesh(cells_per_dimension, mapping; mimetic = true, exact_jacobian = exact_jacobian) # A semidiscretization collects data structures and functions for the spatial discretization - semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_constant, solver) + semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) error_inf, error_L2 = compute_error(solver, semi) - errors_inf[polydeg,2] = error_inf - errors_L2[polydeg,2] = error_L2 + errors_normals_inf[polydeg,2] = error_inf + errors_normals_L2[polydeg,2] = error_L2 + + # Create ODE problem with time span from 0.0 to 1.0 + ode = semidiscretize(semi, (0.0, final_time)); + + # The AnalysisCallback allows to analyse the solution in regular intervals and prints the results + analysis_callback = AnalysisCallback(semi, interval=100) + + # The StepsizeCallback handles the re-calculation of the maximum Δt after each time step + stepsize_callback = StepsizeCallback(cfl=0.1) + + # Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver + callbacks = CallbackSet(analysis_callback, stepsize_callback) + + ############################################################################### + # run the simulation + + # OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks + sol = solve(ode, Euler(), #, CarpenterKennedy2N54(williamson_condition=false), + dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep=false, callback=callbacks); + + errors = analysis_callback(sol) + + errors_sol_L2[polydeg, 2] = errors.l2[1] + errors_sol_inf[polydeg, 2] = errors.linf[1] end -plot(errors_L2[:,1], yaxis=:log, title = "standard") -plot!(errors_L2[:,2], yaxis=:log, title = "mimetic") +#= plot(errors_normals_inf[:,1], yaxis=:log, label = "standard", linewidth=2, thickness_scaling = 1) +plot!(errors_normals_inf[:,2], yaxis=:log, label = "mimetic", linewidth=2, thickness_scaling = 1) =# + +plot(3:max_polydeg,errors_sol_inf[3:end,1], xaxis=:log, yaxis=:log, label = "standard", linewidth=2, thickness_scaling = 1) +plot!(3:max_polydeg,errors_sol_inf[3:end,2], xaxis=:log, yaxis=:log, label = "mimetic", linewidth=2, thickness_scaling = 1) #= ############################################################################### # ODE solvers, callbacks etc. diff --git a/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics_dirichlet.jl b/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics_dirichlet.jl new file mode 100644 index 00000000000..16e924e2db6 --- /dev/null +++ b/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics_dirichlet.jl @@ -0,0 +1,195 @@ +using LinearAlgebra +using OrdinaryDiffEq +using Trixi +using StaticArrays +using Plots + +############################################################################### +# semidiscretization of the linear advection equation + +advection_velocity = (0.2, -0.7, 0.5) +equations = LinearScalarAdvectionEquation3D(advection_velocity) + +# Mapping as described in https://arxiv.org/abs/2012.12040 +function mapping(xi, eta, zeta) + # Transform input variables between -1 and 1 onto our crazy domain + + y = eta + 0.1 * (cos(pi * xi) * + cos(pi * eta) * + cos(pi * zeta)) + + x = xi + 0.1 * (cos(pi * xi) * + cos(pi * eta) * + cos(pi * zeta)) + + z = zeta + 0.1 * (cos(pi * xi) * + cos(pi * eta) * + cos(pi * zeta)) + + #= x = xi + y = eta + z = zeta =# + return SVector(x, y, z) +end + +function exact_contravariant_vectors!(Ja, xi, eta, zeta) + theta_xi = theta_der1(xi, eta, zeta) + theta_eta = theta_der2(xi, eta, zeta) + theta_zeta = theta_der3(xi, eta, zeta) + Ja[1,1] = 1 + theta_eta + theta_zeta + Ja[1,2] = -theta_xi + Ja[1,3] = -theta_xi + Ja[2,1] = -theta_eta + Ja[2,2] = 1 + theta_xi + theta_zeta + Ja[2,3] = -theta_eta + Ja[3,1] = -theta_zeta + Ja[3,2] = -theta_zeta + Ja[3,3] = 1 + theta_xi + theta_eta +end + +function compute_error(solver, semi) + @unpack nodes, weights = solver.basis + exact_Ja = zero(MMatrix{3, 3, Float64}) + error = zero(Float64) + error_L2 = zero(Float64) + for k in eachnode(solver.basis) + for j in eachnode(solver.basis) + for i in eachnode(solver.basis) + exact_contravariant_vectors!(exact_Ja, nodes[i], nodes[j], nodes[k]) + error = max(error, maximum(abs.(semi.cache.elements.contravariant_vectors[:,:,i,j,k,1] - exact_Ja))) + error_L2 += norm(semi.cache.elements.contravariant_vectors[:,1,i,j,k,1] - exact_Ja[:,1]) * weights[i] * weights[j] * weights[k] + end + end + end + return error, error_L2 / 8 +end + +cells_per_dimension = (1,1,1) + +max_polydeg = 25 + +errors_normals_inf = zeros(max_polydeg,2) +errors_normals_L2 = zeros(max_polydeg,2) +errors_sol_inf = zeros(max_polydeg,2) +errors_sol_L2 = zeros(max_polydeg,2) +exact_jacobian = true +final_time = 3.00636132e-03 +initial_condition = initial_condition_constant +for polydeg in 1:max_polydeg + println("Computing polydeg = ", polydeg) + # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux + solver = DGSEM(polydeg, flux_lax_friedrichs) + + # Create curved mesh with 8 x 8 x 8 elements + mesh = StructuredMesh(cells_per_dimension, mapping; mimetic = false, exact_jacobian = exact_jacobian) + + # A semidiscre tization collects data structures and functions for the spatial discretization + semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions = BoundaryConditionDirichlet(initial_condition)) + + error_inf, error_L2 = compute_error(solver, semi) + errors_normals_inf[polydeg,1] = error_inf + errors_normals_L2[polydeg,1] = error_L2 + + # Create ODE problem with time span from 0.0 to 1.0 + ode = semidiscretize(semi, (0.0, final_time)); + + # The AnalysisCallback allows to analyse the solution in regular intervals and prints the results + analysis_callback = AnalysisCallback(semi, interval=100) + + # The StepsizeCallback handles the re-calculation of the maximum Δt after each time step + stepsize_callback = StepsizeCallback(cfl=0.1) + + # Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver + callbacks = CallbackSet(analysis_callback, stepsize_callback) + + + ############################################################################### + # run the simulation + + # OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks + sol = solve(ode, Euler(), #CarpenterKennedy2N54(williamson_condition=false), + dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep=false, callback=callbacks); + + errors = analysis_callback(sol) + + errors_sol_L2[polydeg, 1] = errors.l2[1] + errors_sol_inf[polydeg, 1] = errors.linf[1] + + # Create curved mesh with 8 x 8 x 8 elements + mesh = StructuredMesh(cells_per_dimension, mapping; mimetic = true, exact_jacobian = exact_jacobian) + + # A semidiscretization collects data structures and functions for the spatial discretization + semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions = BoundaryConditionDirichlet(initial_condition)) + + error_inf, error_L2 = compute_error(solver, semi) + errors_normals_inf[polydeg,2] = error_inf + errors_normals_L2[polydeg,2] = error_L2 + + # Create ODE problem with time span from 0.0 to 1.0 + ode = semidiscretize(semi, (0.0, final_time)); + + # The AnalysisCallback allows to analyse the solution in regular intervals and prints the results + analysis_callback = AnalysisCallback(semi, interval=100) + + # The StepsizeCallback handles the re-calculation of the maximum Δt after each time step + stepsize_callback = StepsizeCallback(cfl=0.1) + + # Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver + callbacks = CallbackSet(analysis_callback, stepsize_callback) + + + ############################################################################### + # run the simulation + + # OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks + sol = solve(ode, Euler(), #, CarpenterKennedy2N54(williamson_condition=false), + dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep=false, callback=callbacks); + + errors = analysis_callback(sol) + + errors_sol_L2[polydeg, 2] = errors.l2[1] + errors_sol_inf[polydeg, 2] = errors.linf[1] +end +#= plot(errors_normals_inf[:,1], yaxis=:log, label = "standard", linewidth=2, thickness_scaling = 1) +plot!(errors_normals_inf[:,2], yaxis=:log, label = "mimetic", linewidth=2, thickness_scaling = 1) =# + + +plot(3:max_polydeg,errors_sol_inf[3:end,1], xaxis=:log, yaxis=:log, label = "standard", linewidth=2, thickness_scaling = 1) +plot!(3:max_polydeg,errors_sol_inf[3:end,2], xaxis=:log, yaxis=:log, label = "mimetic", linewidth=2, thickness_scaling = 1) +#= +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to 1.0 +ode = semidiscretize(semi, (0.0, 1.0)); + +# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup +# and resets the timers +summary_callback = SummaryCallback() + +# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results +analysis_callback = AnalysisCallback(semi, interval=100) + +# The SaveSolutionCallback allows to save the solution to a file in regular intervals +save_solution = SaveSolutionCallback(interval=100, + solution_variables=cons2prim) + +# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step +stepsize_callback = StepsizeCallback(cfl=0.1) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback) + + +############################################################################### +# run the simulation + +# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks +sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), + dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep=false, callback=callbacks); + +# Print the timer summary +summary_callback() =# \ No newline at end of file diff --git a/src/Trixi.jl b/src/Trixi.jl index f140d5758f6..804730f3f46 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -220,6 +220,8 @@ export DG, SurfaceIntegralUpwind, MortarL2 +export theta_der1, theta_der2, theta_der3 + export nelements, nnodes, nvariables, eachelement, eachnode, eachvariable diff --git a/src/meshes/structured_mesh.jl b/src/meshes/structured_mesh.jl index 55e696b4b21..1824d8230bb 100644 --- a/src/meshes/structured_mesh.jl +++ b/src/meshes/structured_mesh.jl @@ -21,6 +21,7 @@ mutable struct StructuredMesh{NDIMS, RealT <: Real} <: AbstractMesh{NDIMS} current_filename::String unsaved_changes::Bool mimetic::Bool + exact_jacobian::Bool end """ @@ -46,6 +47,7 @@ Create a StructuredMesh of the given size and shape that uses `RealT` as coordin """ function StructuredMesh(cells_per_dimension, mapping; RealT = Float64, mimetic = false, + exact_jacobian = false, periodicity = true, unsaved_changes = true, mapping_as_string = mapping2string(mapping, length(cells_per_dimension))) @@ -65,7 +67,7 @@ function StructuredMesh(cells_per_dimension, mapping; RealT = Float64, return StructuredMesh{NDIMS, RealT}(Tuple(cells_per_dimension), mapping, mapping_as_string, periodicity, "", - unsaved_changes, mimetic) + unsaved_changes, mimetic, exact_jacobian) end """ diff --git a/src/solvers/dgsem_structured/containers_3d.jl b/src/solvers/dgsem_structured/containers_3d.jl index b1b8c3e225d..34353e94451 100644 --- a/src/solvers/dgsem_structured/containers_3d.jl +++ b/src/solvers/dgsem_structured/containers_3d.jl @@ -19,7 +19,11 @@ function init_elements!(elements, mesh::StructuredMesh{3}, basis::LobattoLegendr calc_node_coordinates!(node_coordinates, element, cell_x, cell_y, cell_z, mesh.mapping, mesh, basis) - calc_jacobian_matrix!(jacobian_matrix, element, node_coordinates, basis) + if mesh.exact_jacobian + calc_jacobian_matrix_exact!(jacobian_matrix, element, node_coordinates, basis) + else + calc_jacobian_matrix_standard!(jacobian_matrix, element, node_coordinates, basis) + end if mesh.mimetic calc_contravariant_vectors_mimetic!(contravariant_vectors, element, jacobian_matrix, @@ -66,8 +70,55 @@ function calc_node_coordinates!(node_coordinates, element, end end +theta_der1(xi, eta, zeta) = -(0.1 * pi) * sin(pi * xi) * cos(pi * eta) * cos(pi * zeta) +theta_der2(xi, eta, zeta) = -(0.1 * pi) * cos(pi * xi) * sin(pi * eta) * cos(pi * zeta) +theta_der3(xi, eta, zeta) = -(0.1 * pi) * cos(pi * xi) * cos(pi * eta) * sin(pi * zeta) + +# Calculate Jacobian matrix of the mapping from the reference element to the element in the physical domain +function calc_jacobian_matrix_exact!(jacobian_matrix::AbstractArray{<:Any, 6}, element, + node_coordinates, basis) + # for dim in 1:3, j in eachnode(basis), i in eachnode(basis) + # # ∂/∂ξ + # jacobian_matrix[dim, 1, :, i, j, element] = basis.derivative_matrix * node_coordinates[dim, :, i, j, element] + # # ∂/∂η + # jacobian_matrix[dim, 2, i, :, j, element] = basis.derivative_matrix * node_coordinates[dim, i, :, j, element] + # # ∂/∂ζ + # jacobian_matrix[dim, 3, i, j, :, element] = basis.derivative_matrix * node_coordinates[dim, i, j, :, element] + # end + @unpack nodes = basis + @turbo for k in eachnode(basis), j in eachnode(basis), i in eachnode(basis) + + jacobian_matrix[1, 1, i, j, k, element] = 1 + theta_der1(nodes[i], nodes[j], nodes[k]) + jacobian_matrix[2, 1, i, j, k, element] = theta_der1(nodes[i], nodes[j], nodes[k]) + jacobian_matrix[3, 1, i, j, k, element] = theta_der1(nodes[i], nodes[j], nodes[k]) + + jacobian_matrix[1, 2, i, j, k, element] = theta_der2(nodes[i], nodes[j], nodes[k]) + jacobian_matrix[2, 2, i, j, k, element] = 1 + theta_der2(nodes[i], nodes[j], nodes[k]) + jacobian_matrix[3, 2, i, j, k, element] = theta_der2(nodes[i], nodes[j], nodes[k]) + + jacobian_matrix[1, 3, i, j, k, element] = theta_der3(nodes[i], nodes[j], nodes[k]) + jacobian_matrix[2, 3, i, j, k, element] = theta_der3(nodes[i], nodes[j], nodes[k]) + jacobian_matrix[3, 3, i, j, k, element] = 1 + theta_der3(nodes[i], nodes[j], nodes[k]) + + #= jacobian_matrix[1, 1, i, j, k, element] = 1 + theta_der1(nodes[i], nodes[j], nodes[k]) + jacobian_matrix[1, 2, i, j, k, element] = theta_der1(nodes[i], nodes[j], nodes[k]) + jacobian_matrix[1, 3, i, j, k, element] = theta_der1(nodes[i], nodes[j], nodes[k]) + + jacobian_matrix[2, 1, i, j, k, element] = theta_der2(nodes[i], nodes[j], nodes[k]) + jacobian_matrix[2, 2, i, j, k, element] = 1 + theta_der2(nodes[i], nodes[j], nodes[k]) + jacobian_matrix[2, 3, i, j, k, element] = theta_der2(nodes[i], nodes[j], nodes[k]) + + jacobian_matrix[3, 1, i, j, k, element] = theta_der3(nodes[i], nodes[j], nodes[k]) + jacobian_matrix[3, 2, i, j, k, element] = theta_der3(nodes[i], nodes[j], nodes[k]) + jacobian_matrix[3, 3, i, j, k, element] = 1 + theta_der3(nodes[i], nodes[j], nodes[k]) =# + + #jacobian_matrix[:, :, i, j, k, element] = transpose(jacobian_matrix[:, :, i, j, k, element]) + end + +end + # Calculate Jacobian matrix of the mapping from the reference element to the element in the physical domain -function calc_jacobian_matrix!(jacobian_matrix::AbstractArray{<:Any, 6}, element, +function calc_jacobian_matrix_standard!(jacobian_matrix::AbstractArray{<:Any, 6}, element, node_coordinates, basis) # The code below is equivalent to the following matrix multiplications but much faster. # From 9e0235bce29410ad0f71d23fea56fe527920c2b1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Fri, 14 Jul 2023 15:25:49 +0200 Subject: [PATCH 09/19] Implemented multi-element mappings --- src/solvers/dgsem_structured/containers_3d.jl | 129 +++++++++++------- 1 file changed, 77 insertions(+), 52 deletions(-) diff --git a/src/solvers/dgsem_structured/containers_3d.jl b/src/solvers/dgsem_structured/containers_3d.jl index 34353e94451..9cf29171bef 100644 --- a/src/solvers/dgsem_structured/containers_3d.jl +++ b/src/solvers/dgsem_structured/containers_3d.jl @@ -12,6 +12,8 @@ function init_elements!(elements, mesh::StructuredMesh{3}, basis::LobattoLegendr linear_indices = LinearIndices(size(mesh)) + node_coordinates_comp = zeros(3, nnodes(basis)) + # Calculate node coordinates, Jacobian matrix, and inverse Jacobian determinant for cell_z in 1:size(mesh, 3), cell_y in 1:size(mesh, 2), cell_x in 1:size(mesh, 1) element = linear_indices[cell_x, cell_y, cell_z] @@ -19,15 +21,17 @@ function init_elements!(elements, mesh::StructuredMesh{3}, basis::LobattoLegendr calc_node_coordinates!(node_coordinates, element, cell_x, cell_y, cell_z, mesh.mapping, mesh, basis) + calc_node_coordinates_computational!(node_coordinates_comp, cell_x, cell_y, cell_z, mesh, basis) + if mesh.exact_jacobian - calc_jacobian_matrix_exact!(jacobian_matrix, element, node_coordinates, basis) + calc_jacobian_matrix_exact!(jacobian_matrix, element, node_coordinates, basis, node_coordinates_comp) else calc_jacobian_matrix_standard!(jacobian_matrix, element, node_coordinates, basis) end if mesh.mimetic calc_contravariant_vectors_mimetic!(contravariant_vectors, element, jacobian_matrix, - node_coordinates, basis) + node_coordinates, basis, node_coordinates_comp) else calc_contravariant_vectors_standard!(contravariant_vectors, element, jacobian_matrix, node_coordinates, basis) @@ -70,13 +74,34 @@ function calc_node_coordinates!(node_coordinates, element, end end +function calc_node_coordinates_computational!(node_coordinates_comp, cell_x, cell_y, cell_z, mesh, basis) + @unpack nodes = basis + + # Get cell length in reference mesh + dx = 2 / size(mesh, 1) + dy = 2 / size(mesh, 2) + dz = 2 / size(mesh, 3) + + # Calculate node coordinates of reference mesh + cell_x_offset = -1 + (cell_x - 1) * dx + dx / 2 + cell_y_offset = -1 + (cell_y - 1) * dy + dy / 2 + cell_z_offset = -1 + (cell_z - 1) * dz + dz / 2 + + for i in eachnode(basis) + # node_coordinates are the mapped reference node_coordinates + node_coordinates_comp[1, i] = cell_x_offset + dx / 2 * nodes[i] + node_coordinates_comp[2, i] = cell_y_offset + dy / 2 * nodes[i] + node_coordinates_comp[3, i] = cell_z_offset + dz / 2 * nodes[i] + end +end + theta_der1(xi, eta, zeta) = -(0.1 * pi) * sin(pi * xi) * cos(pi * eta) * cos(pi * zeta) theta_der2(xi, eta, zeta) = -(0.1 * pi) * cos(pi * xi) * sin(pi * eta) * cos(pi * zeta) theta_der3(xi, eta, zeta) = -(0.1 * pi) * cos(pi * xi) * cos(pi * eta) * sin(pi * zeta) # Calculate Jacobian matrix of the mapping from the reference element to the element in the physical domain function calc_jacobian_matrix_exact!(jacobian_matrix::AbstractArray{<:Any, 6}, element, - node_coordinates, basis) + node_coordinates, basis, nodes) # for dim in 1:3, j in eachnode(basis), i in eachnode(basis) # # ∂/∂ξ # jacobian_matrix[dim, 1, :, i, j, element] = basis.derivative_matrix * node_coordinates[dim, :, i, j, element] @@ -85,35 +110,34 @@ function calc_jacobian_matrix_exact!(jacobian_matrix::AbstractArray{<:Any, 6}, e # # ∂/∂ζ # jacobian_matrix[dim, 3, i, j, :, element] = basis.derivative_matrix * node_coordinates[dim, i, j, :, element] # end - @unpack nodes = basis @turbo for k in eachnode(basis), j in eachnode(basis), i in eachnode(basis) + jacobian_matrix[1, 1, i, j, k, element] = 1 + theta_der1(nodes[1,i], nodes[2,j], nodes[3,k]) + jacobian_matrix[2, 1, i, j, k, element] = theta_der1(nodes[1,i], nodes[2,j], nodes[3,k]) + jacobian_matrix[3, 1, i, j, k, element] = theta_der1(nodes[1,i], nodes[2,j], nodes[3,k]) - jacobian_matrix[1, 1, i, j, k, element] = 1 + theta_der1(nodes[i], nodes[j], nodes[k]) - jacobian_matrix[2, 1, i, j, k, element] = theta_der1(nodes[i], nodes[j], nodes[k]) - jacobian_matrix[3, 1, i, j, k, element] = theta_der1(nodes[i], nodes[j], nodes[k]) + jacobian_matrix[1, 2, i, j, k, element] = theta_der2(nodes[1,i], nodes[2,j], nodes[3,k]) + jacobian_matrix[2, 2, i, j, k, element] = 1 + theta_der2(nodes[1,i], nodes[2,j], nodes[3,k]) + jacobian_matrix[3, 2, i, j, k, element] = theta_der2(nodes[1,i], nodes[2,j], nodes[3,k]) - jacobian_matrix[1, 2, i, j, k, element] = theta_der2(nodes[i], nodes[j], nodes[k]) - jacobian_matrix[2, 2, i, j, k, element] = 1 + theta_der2(nodes[i], nodes[j], nodes[k]) - jacobian_matrix[3, 2, i, j, k, element] = theta_der2(nodes[i], nodes[j], nodes[k]) + jacobian_matrix[1, 3, i, j, k, element] = theta_der3(nodes[1,i], nodes[2,j], nodes[3,k]) + jacobian_matrix[2, 3, i, j, k, element] = theta_der3(nodes[1,i], nodes[2,j], nodes[3,k]) + jacobian_matrix[3, 3, i, j, k, element] = 1 + theta_der3(nodes[1,i], nodes[2,j], nodes[3,k]) - jacobian_matrix[1, 3, i, j, k, element] = theta_der3(nodes[i], nodes[j], nodes[k]) - jacobian_matrix[2, 3, i, j, k, element] = theta_der3(nodes[i], nodes[j], nodes[k]) - jacobian_matrix[3, 3, i, j, k, element] = 1 + theta_der3(nodes[i], nodes[j], nodes[k]) + #= jacobian_matrix[1, 1, i, j, k, element] = 1 + theta_der1(nodes[1,i], nodes[2,j], nodes[3,k]) + jacobian_matrix[1, 2, i, j, k, element] = theta_der1(nodes[1,i], nodes[2,j], nodes[3,k]) + jacobian_matrix[1, 3, i, j, k, element] = theta_der1(nodes[1,i], nodes[2,j], nodes[3,k]) - #= jacobian_matrix[1, 1, i, j, k, element] = 1 + theta_der1(nodes[i], nodes[j], nodes[k]) - jacobian_matrix[1, 2, i, j, k, element] = theta_der1(nodes[i], nodes[j], nodes[k]) - jacobian_matrix[1, 3, i, j, k, element] = theta_der1(nodes[i], nodes[j], nodes[k]) + jacobian_matrix[2, 1, i, j, k, element] = theta_der2(nodes[1,i], nodes[2,j], nodes[3,k]) + jacobian_matrix[2, 2, i, j, k, element] = 1 + theta_der2(nodes[1,i], nodes[2,j], nodes[3,k]) + jacobian_matrix[2, 3, i, j, k, element] = theta_der2(nodes[1,i], nodes[2,j], nodes[3,k]) - jacobian_matrix[2, 1, i, j, k, element] = theta_der2(nodes[i], nodes[j], nodes[k]) - jacobian_matrix[2, 2, i, j, k, element] = 1 + theta_der2(nodes[i], nodes[j], nodes[k]) - jacobian_matrix[2, 3, i, j, k, element] = theta_der2(nodes[i], nodes[j], nodes[k]) - - jacobian_matrix[3, 1, i, j, k, element] = theta_der3(nodes[i], nodes[j], nodes[k]) - jacobian_matrix[3, 2, i, j, k, element] = theta_der3(nodes[i], nodes[j], nodes[k]) - jacobian_matrix[3, 3, i, j, k, element] = 1 + theta_der3(nodes[i], nodes[j], nodes[k]) =# + jacobian_matrix[3, 1, i, j, k, element] = theta_der3(nodes[1,i], nodes[2,j], nodes[3,k]) + jacobian_matrix[3, 2, i, j, k, element] = theta_der3(nodes[1,i], nodes[2,j], nodes[3,k]) + jacobian_matrix[3, 3, i, j, k, element] = 1 + theta_der3(nodes[1,i], nodes[2,j], nodes[3,k]) =# #jacobian_matrix[:, :, i, j, k, element] = transpose(jacobian_matrix[:, :, i, j, k, element]) end + end @@ -304,8 +328,9 @@ New function to compute contravariant vectors function calc_contravariant_vectors_mimetic!(contravariant_vectors::AbstractArray{<:Any, 6}, element, jacobian_matrix, node_coordinates, - basis::LobattoLegendreBasis) - @unpack derivative_matrix, nodes = basis + basis::LobattoLegendreBasis, + nodes) + @unpack derivative_matrix = basis # Define histopolation (edge) basis functions: V[i,j] = hⱼ(ξᵢ) ... TODO: initialize beforehand... V = zero(MMatrix{polydeg(basis) + 1, polydeg(basis), eltype(derivative_matrix)}) @@ -322,45 +347,45 @@ function calc_contravariant_vectors_mimetic!(contravariant_vectors::AbstractArra for k in eachnode(basis) for j in eachnode(basis) for i in 1:polydeg(basis) - Gbar[1, 1, i, j, k] = (nodes[k] * (theta(nodes[i+1], nodes[j], nodes[k]) - theta(nodes[i], nodes[j], nodes[k]) ) - + 0.5 * (theta(nodes[i+1], nodes[j], nodes[k])^2 - theta(nodes[i], nodes[j], nodes[k])^2 ) ) - Gbar[2, 1, i, j, k] = (nodes[i+1] * theta(nodes[i+1], nodes[j], nodes[k]) - nodes[i] * theta(nodes[i], nodes[j], nodes[k]) - + 0.5 * (theta(nodes[i+1], nodes[j], nodes[k])^2 - theta(nodes[i], nodes[j], nodes[k])^2 ) - - (theta_int1(nodes[i+1], nodes[j], nodes[k]) - theta_int1(nodes[i], nodes[j], nodes[k]))) - Gbar[3, 1, i, j, k] = ( nodes[j] * (nodes[i+1] - nodes[i]) - + nodes[j] * ( theta(nodes[i+1], nodes[j], nodes[k]) - theta(nodes[i], nodes[j], nodes[k]) ) - + 0.5 * (theta(nodes[i+1], nodes[j], nodes[k])^2 - theta(nodes[i], nodes[j], nodes[k])^2 ) - + (theta_int1(nodes[i+1], nodes[j], nodes[k]) - theta_int1(nodes[i], nodes[j], nodes[k]))) + Gbar[1, 1, i, j, k] = (nodes[3, k] * (theta(nodes[1,i+1], nodes[2,j], nodes[3,k]) - theta(nodes[1,i], nodes[2,j], nodes[3,k]) ) + + 0.5 * (theta(nodes[1,i+1], nodes[2,j], nodes[3,k])^2 - theta(nodes[1,i], nodes[2,j], nodes[3,k])^2 ) ) + Gbar[2, 1, i, j, k] = (nodes[1, i+1] * theta(nodes[1,i+1], nodes[2,j], nodes[3,k]) - nodes[1, i] * theta(nodes[1,i], nodes[2,j], nodes[3,k]) + + 0.5 * (theta(nodes[1,i+1], nodes[2,j], nodes[3,k])^2 - theta(nodes[1,i], nodes[2,j], nodes[3,k])^2 ) + - (theta_int1(nodes[1,i+1], nodes[2,j], nodes[3,k]) - theta_int1(nodes[1,i], nodes[2,j], nodes[3,k]))) + Gbar[3, 1, i, j, k] = ( nodes[2, j] * (nodes[1, i+1] - nodes[1, i]) + + nodes[2, j] * ( theta(nodes[1,i+1], nodes[2,j], nodes[3,k]) - theta(nodes[1,i], nodes[2,j], nodes[3,k]) ) + + 0.5 * (theta(nodes[1,i+1], nodes[2,j], nodes[3,k])^2 - theta(nodes[1,i], nodes[2,j], nodes[3,k])^2 ) + + (theta_int1(nodes[1,i+1], nodes[2,j], nodes[3,k]) - theta_int1(nodes[1,i], nodes[2,j], nodes[3,k]))) end end end for k in eachnode(basis) for j in 1:polydeg(basis) for i in eachnode(basis) - Gbar[1, 2, i, j, k] = ( nodes[k] * (nodes[j+1] - nodes[j]) - + nodes[k] * ( theta(nodes[i], nodes[j+1], nodes[k]) - theta(nodes[i], nodes[j], nodes[k]) ) - + 0.5 * (theta(nodes[i], nodes[j+1], nodes[k])^2 - theta(nodes[i], nodes[j], nodes[k])^2 ) - + (theta_int2(nodes[i], nodes[j+1], nodes[k]) - theta_int2(nodes[i], nodes[j], nodes[k]))) - Gbar[2, 2, i, j, k] = (nodes[i] * (theta(nodes[i], nodes[j+1], nodes[k]) - theta(nodes[i], nodes[j], nodes[k]) ) - + 0.5 * (theta(nodes[i], nodes[j+1], nodes[k])^2 - theta(nodes[i], nodes[j], nodes[k])^2 ) ) - Gbar[3, 2, i, j, k] = (nodes[j+1] * theta(nodes[i], nodes[j+1], nodes[k]) - nodes[j] * theta(nodes[i], nodes[j], nodes[k]) - + 0.5 * (theta(nodes[i], nodes[j+1], nodes[k])^2 - theta(nodes[i], nodes[j], nodes[k])^2 ) - - (theta_int2(nodes[i], nodes[j+1], nodes[k]) - theta_int2(nodes[i], nodes[j], nodes[k]))) + Gbar[1, 2, i, j, k] = ( nodes[3, k] * (nodes[2, j+1] - nodes[2, j]) + + nodes[3, k] * ( theta(nodes[1,i], nodes[2,j+1], nodes[3,k]) - theta(nodes[1,i], nodes[2,j], nodes[3,k]) ) + + 0.5 * (theta(nodes[1,i], nodes[2,j+1], nodes[3,k])^2 - theta(nodes[1,i], nodes[2,j], nodes[3,k])^2 ) + + (theta_int2(nodes[1,i], nodes[2,j+1], nodes[3,k]) - theta_int2(nodes[1,i], nodes[2,j], nodes[3,k]))) + Gbar[2, 2, i, j, k] = (nodes[1, i] * (theta(nodes[1,i], nodes[2,j+1], nodes[3,k]) - theta(nodes[1,i], nodes[2,j], nodes[3,k]) ) + + 0.5 * (theta(nodes[1,i], nodes[2,j+1], nodes[3,k])^2 - theta(nodes[1,i], nodes[2,j], nodes[3,k])^2 ) ) + Gbar[3, 2, i, j, k] = (nodes[2, j+1] * theta(nodes[1,i], nodes[2,j+1], nodes[3,k]) - nodes[2, j] * theta(nodes[1,i], nodes[2,j], nodes[3,k]) + + 0.5 * (theta(nodes[1,i], nodes[2,j+1], nodes[3,k])^2 - theta(nodes[1,i], nodes[2,j], nodes[3,k])^2 ) + - (theta_int2(nodes[1,i], nodes[2,j+1], nodes[3,k]) - theta_int2(nodes[1,i], nodes[2,j], nodes[3,k]))) end end end for k in 1:polydeg(basis) for j in eachnode(basis) for i in eachnode(basis) - Gbar[1, 3, i, j, k] = (nodes[k+1] * theta(nodes[i], nodes[j], nodes[k+1]) - nodes[k] * theta(nodes[i], nodes[j], nodes[k]) - + 0.5 * (theta(nodes[i], nodes[j], nodes[k+1])^2 - theta(nodes[i], nodes[j], nodes[k])^2 ) - - (theta_int3(nodes[i], nodes[j], nodes[k+1]) - theta_int3(nodes[i], nodes[j], nodes[k]))) - Gbar[2, 3, i, j, k] = ( nodes[i] * (nodes[k+1] - nodes[k]) - + nodes[i] * ( theta(nodes[i], nodes[j], nodes[k+1]) - theta(nodes[i], nodes[j], nodes[k]) ) - + 0.5 * (theta(nodes[i], nodes[j], nodes[k+1])^2 - theta(nodes[i], nodes[j], nodes[k])^2 ) - + (theta_int3(nodes[i], nodes[j], nodes[k+1]) - theta_int3(nodes[i], nodes[j], nodes[k]))) - Gbar[3, 3, i, j, k] = (nodes[j] * (theta(nodes[i], nodes[j], nodes[k+1]) - theta(nodes[i], nodes[j], nodes[k]) ) - + 0.5 * (theta(nodes[i], nodes[j], nodes[k+1])^2 - theta(nodes[i], nodes[j], nodes[k])^2 ) ) + Gbar[1, 3, i, j, k] = (nodes[3, k+1] * theta(nodes[1,i], nodes[2,j], nodes[3,k+1]) - nodes[3, k] * theta(nodes[1,i], nodes[2,j], nodes[3,k]) + + 0.5 * (theta(nodes[1,i], nodes[2,j], nodes[3,k+1])^2 - theta(nodes[1,i], nodes[2,j], nodes[3,k])^2 ) + - (theta_int3(nodes[1,i], nodes[2,j], nodes[3,k+1]) - theta_int3(nodes[1,i], nodes[2,j], nodes[3,k]))) + Gbar[2, 3, i, j, k] = ( nodes[1, i] * (nodes[3, k+1] - nodes[3, k]) + + nodes[1, i] * ( theta(nodes[1,i], nodes[2,j], nodes[3,k+1]) - theta(nodes[1,i], nodes[2,j], nodes[3,k]) ) + + 0.5 * (theta(nodes[1,i], nodes[2,j], nodes[3,k+1])^2 - theta(nodes[1,i], nodes[2,j], nodes[3,k])^2 ) + + (theta_int3(nodes[1,i], nodes[2,j], nodes[3,k+1]) - theta_int3(nodes[1,i], nodes[2,j], nodes[3,k]))) + Gbar[3, 3, i, j, k] = (nodes[2, j] * (theta(nodes[1,i], nodes[2,j], nodes[3,k+1]) - theta(nodes[1,i], nodes[2,j], nodes[3,k]) ) + + 0.5 * (theta(nodes[1,i], nodes[2,j], nodes[3,k+1])^2 - theta(nodes[1,i], nodes[2,j], nodes[3,k])^2 ) ) end end end From d9ccd49b1fa584c6dca21600320f885fb4401cbe Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Fri, 14 Jul 2023 16:14:05 +0200 Subject: [PATCH 10/19] Added computation of the metric terms with standard curl (as opposed to defaut curl invariant) form --- src/solvers/dgsem_structured/containers_3d.jl | 113 ++++++++++++++++++ 1 file changed, 113 insertions(+) diff --git a/src/solvers/dgsem_structured/containers_3d.jl b/src/solvers/dgsem_structured/containers_3d.jl index 9cf29171bef..f7d7b4477e6 100644 --- a/src/solvers/dgsem_structured/containers_3d.jl +++ b/src/solvers/dgsem_structured/containers_3d.jl @@ -322,6 +322,119 @@ function calc_contravariant_vectors_standard!(contravariant_vectors::AbstractArr return contravariant_vectors end +# Calculate contravariant vectors, multiplied by the Jacobian determinant J of the transformation mapping, +# using the invariant curl form. +# These are called Ja^i in Kopriva's blue book. +function calc_contravariant_vectors_standard_curl!(contravariant_vectors::AbstractArray{<:Any, 6}, + element, + jacobian_matrix, node_coordinates, + basis::LobattoLegendreBasis) + @unpack derivative_matrix = basis + + # The general form is + # Jaⁱₙ = 0.5 * ( ∇ × (Xₘ ∇ Xₗ - Xₗ ∇ Xₘ) )ᵢ where (n, m, l) cyclic and ∇ = (∂/∂ξ, ∂/∂η, ∂/∂ζ)ᵀ + + for n in 1:3 + # (n, m, l) cyclic + m = (n % 3) + 1 + l = ((n + 1) % 3) + 1 + + # Calculate Ja¹ₙ = 0.5 * [ (Xₘ Xₗ_ζ - Xₗ Xₘ_ζ)_η - (Xₘ Xₗ_η - Xₗ Xₘ_η)_ζ ] + # For each of these, the first and second summand are computed in separate loops + # for performance reasons. + + # First summand 0.5 * (Xₘ Xₗ_ζ - Xₗ Xₘ_ζ)_η + @turbo for k in eachnode(basis), j in eachnode(basis), i in eachnode(basis) + result = zero(eltype(contravariant_vectors)) + + for ii in eachnode(basis) + # Multiply derivative_matrix to j-dimension to differentiate wrt η + result += derivative_matrix[j, ii] * + ( - node_coordinates[l, i, ii, k, element] * + jacobian_matrix[m, 3, i, ii, k, element]) + end + + contravariant_vectors[n, 1, i, j, k, element] = result + end + + # Second summand -0.5 * (Xₘ Xₗ_η - Xₗ Xₘ_η)_ζ + @turbo for k in eachnode(basis), j in eachnode(basis), i in eachnode(basis) + result = zero(eltype(contravariant_vectors)) + + for ii in eachnode(basis) + # Multiply derivative_matrix to k-dimension to differentiate wrt ζ + result += derivative_matrix[k, ii] * + ( -node_coordinates[l, i, j, ii, element] * + jacobian_matrix[m, 2, i, j, ii, element]) + end + + contravariant_vectors[n, 1, i, j, k, element] -= result + end + + # Calculate Ja²ₙ = 0.5 * [ (Xₘ Xₗ_ξ - Xₗ Xₘ_ξ)_ζ - (Xₘ Xₗ_ζ - Xₗ Xₘ_ζ)_ξ ] + + # First summand 0.5 * (Xₘ Xₗ_ξ - Xₗ Xₘ_ξ)_ζ + @turbo for k in eachnode(basis), j in eachnode(basis), i in eachnode(basis) + result = zero(eltype(contravariant_vectors)) + + for ii in eachnode(basis) + # Multiply derivative_matrix to k-dimension to differentiate wrt ζ + result += derivative_matrix[k, ii] * + (-node_coordinates[l, i, j, ii, element] * + jacobian_matrix[m, 1, i, j, ii, element]) + end + + contravariant_vectors[n, 2, i, j, k, element] = result + end + + # Second summand -0.5 * (Xₘ Xₗ_ζ - Xₗ Xₘ_ζ)_ξ + @turbo for k in eachnode(basis), j in eachnode(basis), i in eachnode(basis) + result = zero(eltype(contravariant_vectors)) + + for ii in eachnode(basis) + # Multiply derivative_matrix to i-dimension to differentiate wrt ξ + result += derivative_matrix[i, ii] * + (-node_coordinates[l, ii, j, k, element] * + jacobian_matrix[m, 3, ii, j, k, element]) + end + + contravariant_vectors[n, 2, i, j, k, element] -= result + end + + # Calculate Ja³ₙ = 0.5 * [ (Xₘ Xₗ_η - Xₗ Xₘ_η)_ξ - (Xₘ Xₗ_ξ - Xₗ Xₘ_ξ)_η ] + + # First summand 0.5 * (Xₘ Xₗ_η - Xₗ Xₘ_η)_ξ + @turbo for k in eachnode(basis), j in eachnode(basis), i in eachnode(basis) + result = zero(eltype(contravariant_vectors)) + + for ii in eachnode(basis) + # Multiply derivative_matrix to i-dimension to differentiate wrt ξ + result += derivative_matrix[i, ii] * + (-node_coordinates[l, ii, j, k, element] * + jacobian_matrix[m, 2, ii, j, k, element]) + end + + contravariant_vectors[n, 3, i, j, k, element] = result + end + + # Second summand -0.5 * (Xₘ Xₗ_ξ - Xₗ Xₘ_ξ)_η + @turbo for k in eachnode(basis), j in eachnode(basis), i in eachnode(basis) + result = zero(eltype(contravariant_vectors)) + + for ii in eachnode(basis) + # Multiply derivative_matrix to j-dimension to differentiate wrt η + result += derivative_matrix[j, ii] * + (-node_coordinates[l, i, ii, k, element] * + jacobian_matrix[m, 1, i, ii, k, element]) + end + + contravariant_vectors[n, 3, i, j, k, element] -= result + end + end + + return contravariant_vectors +end + """ New function to compute contravariant vectors """ From ee9a0b23a57358272ef3d122f1001190b61e9813 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Fri, 21 Jul 2023 16:34:08 +0200 Subject: [PATCH 11/19] First working version of the mimetic metrics and exact_jacobian with the p4estMesh FOR CONFORMING STRUCTURED MESHES -> Single and multi-element cases work -> The exact Jacobian and mimetic metrics are only implemented for the warped mesh test -> Since this is a structured mesh, it must be initialized with a single tree in each direction! -> Uniform refinement works (conforming meshes) but non-conforming meshes do not work yet! --- .../elixir_test_mimetic_metrics.jl | 96 +++++++++++++++++++ src/meshes/p4est_mesh.jl | 14 ++- src/solvers/dgsem_p4est/containers_3d.jl | 53 ++++++++-- src/solvers/dgsem_structured/containers_3d.jl | 30 ++---- 4 files changed, 161 insertions(+), 32 deletions(-) create mode 100644 examples/p4est_3d_dgsem/elixir_test_mimetic_metrics.jl diff --git a/examples/p4est_3d_dgsem/elixir_test_mimetic_metrics.jl b/examples/p4est_3d_dgsem/elixir_test_mimetic_metrics.jl new file mode 100644 index 00000000000..76d947173f0 --- /dev/null +++ b/examples/p4est_3d_dgsem/elixir_test_mimetic_metrics.jl @@ -0,0 +1,96 @@ +using LinearAlgebra +using OrdinaryDiffEq +using Trixi +using StaticArrays +using Plots + +############################################################################### +# semidiscretization of the linear advection equation + +advection_velocity = (0.2, -0.7, 0.5) +equations = LinearScalarAdvectionEquation3D(advection_velocity) + +# Mapping as described in https://arxiv.org/abs/2012.12040 +function mapping(xi, eta, zeta) + # Transform input variables between -1 and 1 onto our crazy domain + + y = eta + 0.1 * (cos(pi * xi) * + cos(pi * eta) * + cos(pi * zeta)) + + x = xi + 0.1 * (cos(pi * xi) * + cos(pi * eta) * + cos(pi * zeta)) + + z = zeta + 0.1 * (cos(pi * xi) * + cos(pi * eta) * + cos(pi * zeta)) + + #= x = xi + y = eta + z = zeta =# + return SVector(x, y, z) +end + +cells_per_dimension = (1, 1, 1) # The p4est implementation works with one tree per direction only + +polydeg_geo = 10 +polydeg = 5 +exact_jacobian = true +final_time = 1.0 + +initial_condition = initial_condition_constant + +solver = DGSEM(polydeg, flux_lax_friedrichs) + +# Create curved mesh with 8 x 8 x 8 elements +#mesh = StructuredMesh(cells_per_dimension, mapping; mimetic = false, exact_jacobian = exact_jacobian) +mesh = P4estMesh(cells_per_dimension; polydeg = polydeg_geo, mapping = mapping, mimetic = true, exact_jacobian = exact_jacobian, initial_refinement_level = 3) + +# A semidiscre tization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +# Create ODE problem with time span from 0.0 to 1.0 +ode = semidiscretize(semi, (0.0, final_time)); + +summary_callback = SummaryCallback() + +# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results +analysis_callback = AnalysisCallback(semi, interval=100) + +# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step +stepsize_callback = StepsizeCallback(cfl=0.1) + +#= amr_indicator = IndicatorHennemannGassner(semi, + alpha_max=1.0, + alpha_min=0.0001, + alpha_smooth=false, + variable=Trixi.energy_total) + +amr_controller = ControllerThreeLevel(semi, amr_indicator, + base_level=4, + max_level=6, max_threshold=0.01) + +#= amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable=first), + base_level=4, + med_level=5, med_threshold=0.1, + max_level=6, max_threshold=0.6) =# +amr_callback = AMRCallback(semi, amr_controller, + interval=5, + adapt_initial_condition=true, + adapt_initial_condition_only_refine=true) =# + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +callbacks = CallbackSet(summary_callback, analysis_callback, stepsize_callback) + + +############################################################################### +# run the simulation + +# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks +sol = solve(ode, Euler(), #CarpenterKennedy2N54(williamson_condition=false), + dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep=false, callback=callbacks); + + +summary_callback() \ No newline at end of file diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 60db285e04f..033dc3d0c83 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -24,10 +24,12 @@ mutable struct P4estMesh{NDIMS, RealT <: Real, IsParallel, P, Ghost, NDIMSP2, NN current_filename::String unsaved_changes::Bool p4est_partition_allow_for_coarsening::Bool + mimetic::Bool + exact_jacobian::Bool function P4estMesh{NDIMS}(p4est, tree_node_coordinates, nodes, boundary_names, current_filename, unsaved_changes, - p4est_partition_allow_for_coarsening) where {NDIMS} + p4est_partition_allow_for_coarsening, mimetic = false, exact_jacobian = false) where {NDIMS} if NDIMS == 2 @assert p4est isa Ptr{p4est_t} elseif NDIMS == 3 @@ -57,7 +59,9 @@ mutable struct P4estMesh{NDIMS, RealT <: Real, IsParallel, P, Ghost, NDIMSP2, NN boundary_names, current_filename, unsaved_changes, - p4est_partition_allow_for_coarsening) + p4est_partition_allow_for_coarsening, + mimetic, + exact_jacobian) # Destroy `p4est` structs when the mesh is garbage collected finalizer(destroy_mesh, mesh) @@ -166,7 +170,9 @@ function P4estMesh(trees_per_dimension; polydeg, coordinates_max = nothing, RealT = Float64, initial_refinement_level = 0, periodicity = true, unsaved_changes = true, - p4est_partition_allow_for_coarsening = true) + p4est_partition_allow_for_coarsening = true, + mimetic = false, + exact_jacobian = false) @assert ((coordinates_min === nothing)===(coordinates_max === nothing)) "Either both or none of coordinates_min and coordinates_max must be specified" @assert count(i -> i !== nothing, @@ -215,7 +221,7 @@ function P4estMesh(trees_per_dimension; polydeg, return P4estMesh{NDIMS}(p4est, tree_node_coordinates, nodes, boundary_names, "", unsaved_changes, - p4est_partition_allow_for_coarsening) + p4est_partition_allow_for_coarsening, mimetic, exact_jacobian) end # 2D version diff --git a/src/solvers/dgsem_p4est/containers_3d.jl b/src/solvers/dgsem_p4est/containers_3d.jl index e9994fe4569..8c1f6f98748 100644 --- a/src/solvers/dgsem_p4est/containers_3d.jl +++ b/src/solvers/dgsem_p4est/containers_3d.jl @@ -12,13 +12,43 @@ function init_elements!(elements, mesh::P4estMesh{3}, basis::LobattoLegendreBasi calc_node_coordinates!(node_coordinates, mesh, basis) - for element in 1:ncells(mesh) - calc_jacobian_matrix!(jacobian_matrix, element, node_coordinates, basis) + # Macros from `p4est` + p4est_root_len = 1 << P4EST_MAXLEVEL + p4est_quadrant_len(l) = 1 << (P4EST_MAXLEVEL - l) + + trees = unsafe_wrap_sc(p8est_tree_t, mesh.p4est.trees) - calc_contravariant_vectors!(contravariant_vectors, element, jacobian_matrix, - node_coordinates, basis) + node_coordinates_comp = zeros(3, nnodes(basis)) + + for tree in eachindex(trees) + offset = trees[tree].quadrants_offset + quadrants = unsafe_wrap_sc(p8est_quadrant_t, trees[tree].quadrants) - calc_inverse_jacobian!(inverse_jacobian, element, jacobian_matrix, basis) + for i in eachindex(quadrants) + element = offset + i + + if mesh.exact_jacobian || mesh.mimetic + quad = quadrants[i] + quad_length = p4est_quadrant_len(quad.level) / p4est_root_len + calc_node_coordinates_computational!(node_coordinates_comp, quad_length, p4est_root_len, quad, mesh, basis) + end + + if mesh.exact_jacobian + calc_jacobian_matrix_exact!(jacobian_matrix, element, node_coordinates, basis, node_coordinates_comp) + else + calc_jacobian_matrix!(jacobian_matrix, element, node_coordinates, basis) + end + + if mesh.mimetic + calc_contravariant_vectors_mimetic!(contravariant_vectors, element, jacobian_matrix, + node_coordinates, basis, node_coordinates_comp) + else + calc_contravariant_vectors!(contravariant_vectors, element, jacobian_matrix, + node_coordinates, basis) + end + + calc_inverse_jacobian!(inverse_jacobian, element, jacobian_matrix, basis) + end end return nothing @@ -30,7 +60,7 @@ function calc_node_coordinates!(node_coordinates, basis::LobattoLegendreBasis) # Hanging nodes will cause holes in the mesh if its polydeg is higher # than the polydeg of the solver. - @assert length(basis.nodes)>=length(mesh.nodes) "The solver can't have a lower polydeg than the mesh" + #@assert length(basis.nodes)>=length(mesh.nodes) "The solver can't have a lower polydeg than the mesh" calc_node_coordinates!(node_coordinates, mesh, basis.nodes) end @@ -75,6 +105,17 @@ function calc_node_coordinates!(node_coordinates, return node_coordinates end +function calc_node_coordinates_computational!(node_coordinates_comp, quad_length, p4est_root_len, quad, mesh, basis) + @unpack nodes = basis + + node_coordinates_comp[1,:] = 2 * (quad_length * 1 / 2 * (nodes .+ 1) .+ + quad.x / p4est_root_len) .- 1 + node_coordinates_comp[2,:] = 2 * (quad_length * 1 / 2 * (nodes .+ 1) .+ + quad.y / p4est_root_len) .- 1 + node_coordinates_comp[3,:] = 2 * (quad_length * 1 / 2 * (nodes .+ 1) .+ + quad.z / p4est_root_len) .- 1 +end + # Initialize node_indices of interface container @inline function init_interface_node_indices!(interfaces::P4estInterfaceContainer{3}, faces, orientation, interface_id) diff --git a/src/solvers/dgsem_structured/containers_3d.jl b/src/solvers/dgsem_structured/containers_3d.jl index f7d7b4477e6..9b6957f8323 100644 --- a/src/solvers/dgsem_structured/containers_3d.jl +++ b/src/solvers/dgsem_structured/containers_3d.jl @@ -21,19 +21,21 @@ function init_elements!(elements, mesh::StructuredMesh{3}, basis::LobattoLegendr calc_node_coordinates!(node_coordinates, element, cell_x, cell_y, cell_z, mesh.mapping, mesh, basis) - calc_node_coordinates_computational!(node_coordinates_comp, cell_x, cell_y, cell_z, mesh, basis) - + if mesh.exact_jacobian || mesh.mimetic + calc_node_coordinates_computational!(node_coordinates_comp, cell_x, cell_y, cell_z, mesh, basis) + end + if mesh.exact_jacobian calc_jacobian_matrix_exact!(jacobian_matrix, element, node_coordinates, basis, node_coordinates_comp) else - calc_jacobian_matrix_standard!(jacobian_matrix, element, node_coordinates, basis) + calc_jacobian_matrix!(jacobian_matrix, element, node_coordinates, basis) end if mesh.mimetic calc_contravariant_vectors_mimetic!(contravariant_vectors, element, jacobian_matrix, node_coordinates, basis, node_coordinates_comp) else - calc_contravariant_vectors_standard!(contravariant_vectors, element, jacobian_matrix, + calc_contravariant_vectors!(contravariant_vectors, element, jacobian_matrix, node_coordinates, basis) end @@ -122,27 +124,11 @@ function calc_jacobian_matrix_exact!(jacobian_matrix::AbstractArray{<:Any, 6}, e jacobian_matrix[1, 3, i, j, k, element] = theta_der3(nodes[1,i], nodes[2,j], nodes[3,k]) jacobian_matrix[2, 3, i, j, k, element] = theta_der3(nodes[1,i], nodes[2,j], nodes[3,k]) jacobian_matrix[3, 3, i, j, k, element] = 1 + theta_der3(nodes[1,i], nodes[2,j], nodes[3,k]) - - #= jacobian_matrix[1, 1, i, j, k, element] = 1 + theta_der1(nodes[1,i], nodes[2,j], nodes[3,k]) - jacobian_matrix[1, 2, i, j, k, element] = theta_der1(nodes[1,i], nodes[2,j], nodes[3,k]) - jacobian_matrix[1, 3, i, j, k, element] = theta_der1(nodes[1,i], nodes[2,j], nodes[3,k]) - - jacobian_matrix[2, 1, i, j, k, element] = theta_der2(nodes[1,i], nodes[2,j], nodes[3,k]) - jacobian_matrix[2, 2, i, j, k, element] = 1 + theta_der2(nodes[1,i], nodes[2,j], nodes[3,k]) - jacobian_matrix[2, 3, i, j, k, element] = theta_der2(nodes[1,i], nodes[2,j], nodes[3,k]) - - jacobian_matrix[3, 1, i, j, k, element] = theta_der3(nodes[1,i], nodes[2,j], nodes[3,k]) - jacobian_matrix[3, 2, i, j, k, element] = theta_der3(nodes[1,i], nodes[2,j], nodes[3,k]) - jacobian_matrix[3, 3, i, j, k, element] = 1 + theta_der3(nodes[1,i], nodes[2,j], nodes[3,k]) =# - - #jacobian_matrix[:, :, i, j, k, element] = transpose(jacobian_matrix[:, :, i, j, k, element]) end - - end # Calculate Jacobian matrix of the mapping from the reference element to the element in the physical domain -function calc_jacobian_matrix_standard!(jacobian_matrix::AbstractArray{<:Any, 6}, element, +function calc_jacobian_matrix!(jacobian_matrix::AbstractArray{<:Any, 6}, element, node_coordinates, basis) # The code below is equivalent to the following matrix multiplications but much faster. # @@ -200,7 +186,7 @@ end # Calculate contravariant vectors, multiplied by the Jacobian determinant J of the transformation mapping, # using the invariant curl form. # These are called Ja^i in Kopriva's blue book. -function calc_contravariant_vectors_standard!(contravariant_vectors::AbstractArray{<:Any, 6}, +function calc_contravariant_vectors!(contravariant_vectors::AbstractArray{<:Any, 6}, element, jacobian_matrix, node_coordinates, basis::LobattoLegendreBasis) From d18fa7ceaca652459c90c65f777b43b8b568f509 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Fri, 21 Jul 2023 17:21:56 +0200 Subject: [PATCH 12/19] p4est example uses a non-conforming grid --- .../elixir_test_mimetic_metrics.jl | 29 ++++++++++++++++--- 1 file changed, 25 insertions(+), 4 deletions(-) diff --git a/examples/p4est_3d_dgsem/elixir_test_mimetic_metrics.jl b/examples/p4est_3d_dgsem/elixir_test_mimetic_metrics.jl index 76d947173f0..c8a4d00de89 100644 --- a/examples/p4est_3d_dgsem/elixir_test_mimetic_metrics.jl +++ b/examples/p4est_3d_dgsem/elixir_test_mimetic_metrics.jl @@ -34,9 +34,9 @@ end cells_per_dimension = (1, 1, 1) # The p4est implementation works with one tree per direction only -polydeg_geo = 10 +polydeg_geo = 5 polydeg = 5 -exact_jacobian = true +exact_jacobian = false final_time = 1.0 initial_condition = initial_condition_constant @@ -45,7 +45,24 @@ solver = DGSEM(polydeg, flux_lax_friedrichs) # Create curved mesh with 8 x 8 x 8 elements #mesh = StructuredMesh(cells_per_dimension, mapping; mimetic = false, exact_jacobian = exact_jacobian) -mesh = P4estMesh(cells_per_dimension; polydeg = polydeg_geo, mapping = mapping, mimetic = true, exact_jacobian = exact_jacobian, initial_refinement_level = 3) +mesh = P4estMesh(cells_per_dimension; polydeg = polydeg_geo, mapping = mapping, mimetic = false, exact_jacobian = exact_jacobian, initial_refinement_level = 0) + +# Refine bottom left quadrant of each tree to level 3 +function refine_fn(p8est, which_tree, quadrant) + quadrant_obj = unsafe_load(quadrant) + if quadrant_obj.x == 0 && quadrant_obj.y == 0 && quadrant_obj.z == 0 && quadrant_obj.level < 2 + # return true (refine) + return Cint(1) + else + # return false (don't refine) + return Cint(0) + end + end + + # Refine recursively until each bottom left quadrant of a tree has level 3 + # The mesh will be rebalanced before the simulation starts + refine_fn_c = @cfunction(refine_fn, Cint, (Ptr{Trixi.p8est_t}, Ptr{Trixi.p4est_topidx_t}, Ptr{Trixi.p8est_quadrant_t})) + Trixi.refine_p4est!(mesh.p4est, true, refine_fn_c, C_NULL) # A semidiscre tization collects data structures and functions for the spatial discretization semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) @@ -61,6 +78,10 @@ analysis_callback = AnalysisCallback(semi, interval=100) # The StepsizeCallback handles the re-calculation of the maximum Δt after each time step stepsize_callback = StepsizeCallback(cfl=0.1) +# The SaveSolutionCallback allows to save the solution to a file in regular intervals +save_solution = SaveSolutionCallback(interval=100, + solution_variables=cons2prim) + #= amr_indicator = IndicatorHennemannGassner(semi, alpha_max=1.0, alpha_min=0.0001, @@ -81,7 +102,7 @@ amr_callback = AMRCallback(semi, amr_controller, adapt_initial_condition_only_refine=true) =# # Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver -callbacks = CallbackSet(summary_callback, analysis_callback, stepsize_callback) +callbacks = CallbackSet(summary_callback, analysis_callback, stepsize_callback, save_solution) ############################################################################### From 01289339f8274f8c2ffc55999cab783c73e313d8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Mon, 24 Jul 2023 15:50:36 +0200 Subject: [PATCH 13/19] some modifications to the mimetic p4est test --- .../elixir_test_mimetic_metrics.jl | 187 +++++++++++------- src/solvers/dgsem_p4est/containers_3d.jl | 2 +- 2 files changed, 115 insertions(+), 74 deletions(-) diff --git a/examples/p4est_3d_dgsem/elixir_test_mimetic_metrics.jl b/examples/p4est_3d_dgsem/elixir_test_mimetic_metrics.jl index c8a4d00de89..fdbd8d8baeb 100644 --- a/examples/p4est_3d_dgsem/elixir_test_mimetic_metrics.jl +++ b/examples/p4est_3d_dgsem/elixir_test_mimetic_metrics.jl @@ -34,84 +34,125 @@ end cells_per_dimension = (1, 1, 1) # The p4est implementation works with one tree per direction only -polydeg_geo = 5 -polydeg = 5 + exact_jacobian = false -final_time = 1.0 +mimetic = false +final_time = 1e-3 initial_condition = initial_condition_constant -solver = DGSEM(polydeg, flux_lax_friedrichs) - -# Create curved mesh with 8 x 8 x 8 elements -#mesh = StructuredMesh(cells_per_dimension, mapping; mimetic = false, exact_jacobian = exact_jacobian) -mesh = P4estMesh(cells_per_dimension; polydeg = polydeg_geo, mapping = mapping, mimetic = false, exact_jacobian = exact_jacobian, initial_refinement_level = 0) - -# Refine bottom left quadrant of each tree to level 3 -function refine_fn(p8est, which_tree, quadrant) - quadrant_obj = unsafe_load(quadrant) - if quadrant_obj.x == 0 && quadrant_obj.y == 0 && quadrant_obj.z == 0 && quadrant_obj.level < 2 - # return true (refine) - return Cint(1) - else - # return false (don't refine) - return Cint(0) - end - end - - # Refine recursively until each bottom left quadrant of a tree has level 3 - # The mesh will be rebalanced before the simulation starts - refine_fn_c = @cfunction(refine_fn, Cint, (Ptr{Trixi.p8est_t}, Ptr{Trixi.p4est_topidx_t}, Ptr{Trixi.p8est_quadrant_t})) - Trixi.refine_p4est!(mesh.p4est, true, refine_fn_c, C_NULL) - -# A semidiscre tization collects data structures and functions for the spatial discretization -semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) - -# Create ODE problem with time span from 0.0 to 1.0 -ode = semidiscretize(semi, (0.0, final_time)); - -summary_callback = SummaryCallback() - -# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results -analysis_callback = AnalysisCallback(semi, interval=100) - -# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step -stepsize_callback = StepsizeCallback(cfl=0.1) - -# The SaveSolutionCallback allows to save the solution to a file in regular intervals -save_solution = SaveSolutionCallback(interval=100, - solution_variables=cons2prim) - -#= amr_indicator = IndicatorHennemannGassner(semi, - alpha_max=1.0, - alpha_min=0.0001, - alpha_smooth=false, - variable=Trixi.energy_total) +#= polydeg_geo = 15 +polydeg = 10 =# + +max_polydeg = 20 +n_polydeg_geo = 4 +errors_sol_inf = zeros(max_polydeg,n_polydeg_geo) +errors_sol_L2 = zeros(max_polydeg,n_polydeg_geo) +polydeg_geos = [2, 3, 5, 10] +#polydeg_geos = [30] +for i in 1:n_polydeg_geo + polydeg_geo = polydeg_geos[i] + for polydeg in 1:max_polydeg + + solver = DGSEM(polydeg, flux_lax_friedrichs) + + # Create curved mesh with 8 x 8 x 8 elements + boundary_condition = BoundaryConditionDirichlet(initial_condition) + boundary_conditions = Dict( + :x_neg => boundary_condition, + :x_pos => boundary_condition, + :y_neg => boundary_condition, + :y_pos => boundary_condition, + :z_neg => boundary_condition, + :z_pos => boundary_condition + ) + println("polydeg_geo: ", polydeg_geo) + #mesh = P4estMesh(cells_per_dimension; polydeg = polydeg_geo, mapping = mapping, mimetic = mimetic, exact_jacobian = exact_jacobian, initial_refinement_level = 0, periodicity = true) + mesh = P4estMesh(cells_per_dimension; polydeg = polydeg_geo, mapping = mapping, mimetic = mimetic, exact_jacobian = exact_jacobian, initial_refinement_level = 0, periodicity = false) + + # Refine bottom left quadrant of each tree to level 3 + function refine_fn(p8est, which_tree, quadrant) + quadrant_obj = unsafe_load(quadrant) + if quadrant_obj.x == 0 && quadrant_obj.y == 0 && quadrant_obj.z == 0 && quadrant_obj.level < 2 + # return true (refine) + return Cint(1) + else + # return false (don't refine) + return Cint(0) + end + end + + # Refine recursively until each bottom left quadrant of a tree has level 3 + # The mesh will be rebalanced before the simulation starts + refine_fn_c = @cfunction(refine_fn, Cint, (Ptr{Trixi.p8est_t}, Ptr{Trixi.p4est_topidx_t}, Ptr{Trixi.p8est_quadrant_t})) + Trixi.refine_p4est!(mesh.p4est, true, refine_fn_c, C_NULL) + + # A semidiscre tization collects data structures and functions for the spatial discretization + #semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions = boundary_conditions) + + # Create ODE problem with time span from 0.0 to 1.0 + ode = semidiscretize(semi, (0.0, final_time)); + + summary_callback = SummaryCallback() + + # The AnalysisCallback allows to analyse the solution in regular intervals and prints the results + analysis_callback = AnalysisCallback(semi, interval=100) + + # The StepsizeCallback handles the re-calculation of the maximum Δt after each time step + stepsize_callback = StepsizeCallback(cfl=0.1) + + # The SaveSolutionCallback allows to save the solution to a file in regular intervals + save_solution = SaveSolutionCallback(interval=100, + solution_variables=cons2prim) + + #= amr_indicator = IndicatorHennemannGassner(semi, + alpha_max=1.0, + alpha_min=0.0001, + alpha_smooth=false, + variable=Trixi.energy_total) + + amr_controller = ControllerThreeLevel(semi, amr_indicator, + base_level=4, + max_level=6, max_threshold=0.01) + + #= amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable=first), + base_level=4, + med_level=5, med_threshold=0.1, + max_level=6, max_threshold=0.6) =# + amr_callback = AMRCallback(semi, amr_controller, + interval=5, + adapt_initial_condition=true, + adapt_initial_condition_only_refine=true) =# + + # Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver + callbacks = CallbackSet(summary_callback, analysis_callback, stepsize_callback, save_solution) + + + ############################################################################### + # run the simulation + + # OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks + sol = solve(ode, Euler(), #CarpenterKennedy2N54(williamson_condition=false), + dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep=false, callback=callbacks); + + + summary_callback() + errors = analysis_callback(sol) + + errors_sol_L2[polydeg, i] = errors.l2[1] + errors_sol_inf[polydeg, i] = errors.linf[1] -amr_controller = ControllerThreeLevel(semi, amr_indicator, - base_level=4, - max_level=6, max_threshold=0.01) - -#= amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable=first), - base_level=4, - med_level=5, med_threshold=0.1, - max_level=6, max_threshold=0.6) =# -amr_callback = AMRCallback(semi, amr_controller, - interval=5, - adapt_initial_condition=true, - adapt_initial_condition_only_refine=true) =# - -# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver -callbacks = CallbackSet(summary_callback, analysis_callback, stepsize_callback, save_solution) - - -############################################################################### -# run the simulation + end +end -# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks -sol = solve(ode, Euler(), #CarpenterKennedy2N54(williamson_condition=false), - dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback - save_everystep=false, callback=callbacks); +for i in 1:n_polydeg_geo + plot!(errors_sol_inf[:,i], xaxis=:log, yaxis=:log, label = "polydeg_geo="*string(polydeg_geos[i]), linewidth=2, thickness_scaling = 1) +end +plot!(title = "mimetic="*string(mimetic)*", exact_jacobian="*string(exact_jacobian)) +plot!(xlabel = "polydeg", ylabel = "|u_ex - u_disc|_inf") +plot!(ylims=(1e-15,1e-1)) -summary_callback() \ No newline at end of file +plot!(xticks=([2, 4, 8, 16], ["2", "4", "8", "16"])) \ No newline at end of file diff --git a/src/solvers/dgsem_p4est/containers_3d.jl b/src/solvers/dgsem_p4est/containers_3d.jl index 8c1f6f98748..5fb984b82bb 100644 --- a/src/solvers/dgsem_p4est/containers_3d.jl +++ b/src/solvers/dgsem_p4est/containers_3d.jl @@ -105,7 +105,7 @@ function calc_node_coordinates!(node_coordinates, return node_coordinates end -function calc_node_coordinates_computational!(node_coordinates_comp, quad_length, p4est_root_len, quad, mesh, basis) +function calc_node_coordinates_computational!(node_coordinates_comp, quad_length, p4est_root_len, quad, mesh::P4estMesh{3}, basis) @unpack nodes = basis node_coordinates_comp[1,:] = 2 * (quad_length * 1 / 2 * (nodes .+ 1) .+ From e652a266b6a1852a8f4c645adf08f04edbf8d785 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Thu, 27 Jul 2023 10:05:28 +0200 Subject: [PATCH 14/19] First implementation of inherited metric terms. ONLY WORKS FOR polydeg_geo = polydeg_parent_metrics WHEN mimetic = true OR exact_jacobian = true! --- .../elixir_test_mimetic_metrics_parent.jl | 158 ++++++++++++++++++ src/meshes/p4est_mesh.jl | 58 ++++++- src/solvers/dgsem_p4est/containers_3d.jl | 62 +++++-- 3 files changed, 262 insertions(+), 16 deletions(-) create mode 100644 examples/p4est_3d_dgsem/elixir_test_mimetic_metrics_parent.jl diff --git a/examples/p4est_3d_dgsem/elixir_test_mimetic_metrics_parent.jl b/examples/p4est_3d_dgsem/elixir_test_mimetic_metrics_parent.jl new file mode 100644 index 00000000000..be7faf6c695 --- /dev/null +++ b/examples/p4est_3d_dgsem/elixir_test_mimetic_metrics_parent.jl @@ -0,0 +1,158 @@ +#using LinearAlgebra +using OrdinaryDiffEq +using Trixi +#using StaticArrays +#using Plots + +############################################################################### +# semidiscretization of the linear advection equation + +advection_velocity = (0.2, -0.7, 0.5) +equations = LinearScalarAdvectionEquation3D(advection_velocity) + +# Mapping as described in https://arxiv.org/abs/2012.12040 +function mapping(xi, eta, zeta) + # Transform input variables between -1 and 1 onto our crazy domain + + y = eta + 0.1 * (cos(pi * xi) * + cos(pi * eta) * + cos(pi * zeta)) + + x = xi + 0.1 * (cos(pi * xi) * + cos(pi * y) * + cos(pi * zeta)) + + z = zeta + 0.1 * (cos(pi * x) * + cos(pi * y) * + cos(pi * zeta)) + + #= x = xi + y = eta + z = zeta =# + return SVector(x, y, z) +end + +cells_per_dimension = (1, 1, 1) # The p4est implementation works with one tree per direction only + + +exact_jacobian = false +mimetic = false +final_time = 1e-3 + +initial_condition = initial_condition_constant + +max_polydeg = 1 +n_polydeg_geo = 4 +errors_sol_inf = zeros(max_polydeg,n_polydeg_geo) +errors_sol_L2 = zeros(max_polydeg,n_polydeg_geo) +polydeg_geos = [2, 3, 5, 10] +#polydeg_geos = [5] + +polydeg = 5 + +# Refine bottom left quadrant of each tree to level 3 +function refine_fn(p8est, which_tree, quadrant) + quadrant_obj = unsafe_load(quadrant) + if quadrant_obj.x == 0 && quadrant_obj.y == 0 && quadrant_obj.z == 0 && quadrant_obj.level < 2 + # return true (refine) + return Cint(1) + else + # return false (don't refine) + return Cint(0) + end +end + +for i in 1:n_polydeg_geo + for polydeg in 1:max_polydeg + polydeg_geo = polydeg_geos[i] + + solver = DGSEM(polydeg, flux_lax_friedrichs) + + # Create curved mesh with 8 x 8 x 8 elements + boundary_condition = BoundaryConditionDirichlet(initial_condition) + boundary_conditions = Dict( + :x_neg => boundary_condition, + :x_pos => boundary_condition, + :y_neg => boundary_condition, + :y_pos => boundary_condition, + :z_neg => boundary_condition, + :z_pos => boundary_condition + ) + println("polydeg_geo: ", polydeg_geo) + #mesh = P4estMesh(cells_per_dimension; polydeg = polydeg_geo, mapping = mapping, mimetic = mimetic, exact_jacobian = exact_jacobian, initial_refinement_level = 0, periodicity = true) + mesh = P4estMesh(cells_per_dimension; polydeg = polydeg_geo, mapping = mapping, mimetic = mimetic, exact_jacobian = exact_jacobian, initial_refinement_level = 0, periodicity = false, polydeg_parent_metrics = polydeg_geo) + + # Refine recursively until each bottom left quadrant of a tree has level 3 + # The mesh will be rebalanced before the simulation starts + refine_fn_c = @cfunction(refine_fn, Cint, (Ptr{Trixi.p8est_t}, Ptr{Trixi.p4est_topidx_t}, Ptr{Trixi.p8est_quadrant_t})) + Trixi.refine_p4est!(mesh.p4est, true, refine_fn_c, C_NULL) + + # A semidiscre tization collects data structures and functions for the spatial discretization + #semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions = boundary_conditions) + + # Create ODE problem with time span from 0.0 to 1.0 + ode = semidiscretize(semi, (0.0, final_time)); + + summary_callback = SummaryCallback() + + # The AnalysisCallback allows to analyse the solution in regular intervals and prints the results + analysis_callback = AnalysisCallback(semi, interval=100) + + # The StepsizeCallback handles the re-calculation of the maximum Δt after each time step + stepsize_callback = StepsizeCallback(cfl=0.1) + + # The SaveSolutionCallback allows to save the solution to a file in regular intervals + save_solution = SaveSolutionCallback(interval=100, + solution_variables=cons2prim) + + #= amr_indicator = IndicatorHennemannGassner(semi, + alpha_max=1.0, + alpha_min=0.0001, + alpha_smooth=false, + variable=Trixi.energy_total) + + amr_controller = ControllerThreeLevel(semi, amr_indicator, + base_level=4, + max_level=6, max_threshold=0.01) + + #= amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable=first), + base_level=4, + med_level=5, med_threshold=0.1, + max_level=6, max_threshold=0.6) =# + amr_callback = AMRCallback(semi, amr_controller, + interval=5, + adapt_initial_condition=true, + adapt_initial_condition_only_refine=true) =# + + # Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver + callbacks = CallbackSet(summary_callback, analysis_callback, stepsize_callback, save_solution) + + + ############################################################################### + # run the simulation + + # OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks + sol = solve(ode, Euler(), #CarpenterKennedy2N54(williamson_condition=false), + dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep=false, callback=callbacks); + + + summary_callback() + errors = analysis_callback(sol) + + errors_sol_L2[polydeg, i] = errors.l2[1] + errors_sol_inf[polydeg, i] = errors.linf[1] + end +end +#=end + +for i in 1:n_polydeg_geo + plot!(errors_sol_inf[:,i], xaxis=:log, yaxis=:log, label = "polydeg_geo="*string(polydeg_geos[i]), linewidth=2, thickness_scaling = 1) +end +plot!(title = "mimetic="*string(mimetic)*", exact_jacobian="*string(exact_jacobian)) +plot!(xlabel = "polydeg", ylabel = "|u_ex - u_disc|_inf") + +plot!(ylims=(1e-15,1e-1)) + +plot!(xticks=([2, 4, 8, 16], ["2", "4", "8", "16"])) =# \ No newline at end of file diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 033dc3d0c83..b34057a9d2f 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -11,7 +11,7 @@ An unstructured curved mesh based on trees that uses the C library `p4est` to manage trees and mesh refinement. """ -mutable struct P4estMesh{NDIMS, RealT <: Real, IsParallel, P, Ghost, NDIMSP2, NNODES} <: +mutable struct P4estMesh{NDIMS, RealT <: Real, IsParallel, P, Ghost, NDIMSP2, NDIMSP3, NNODES} <: AbstractMesh{NDIMS} p4est :: P # Either PointerWrapper{p4est_t} or PointerWrapper{p8est_t} is_parallel :: IsParallel @@ -26,10 +26,12 @@ mutable struct P4estMesh{NDIMS, RealT <: Real, IsParallel, P, Ghost, NDIMSP2, NN p4est_partition_allow_for_coarsening::Bool mimetic::Bool exact_jacobian::Bool + polydeg_parent_metrics::Int + tree_contravariant_vectors::Array{RealT, NDIMSP3} # [dimension, i, j, k, tree] function P4estMesh{NDIMS}(p4est, tree_node_coordinates, nodes, boundary_names, current_filename, unsaved_changes, - p4est_partition_allow_for_coarsening, mimetic = false, exact_jacobian = false) where {NDIMS} + p4est_partition_allow_for_coarsening, mimetic = false, exact_jacobian = false, polydeg_parent_metrics = false, tree_contravariant_vectors = Array{eltype(tree_node_coordinates), NDIMS + 3}(undef, 3, ntuple(_ -> 0, NDIMS)..., 1)) where {NDIMS} if NDIMS == 2 @assert p4est isa Ptr{p4est_t} elseif NDIMS == 3 @@ -51,7 +53,7 @@ mutable struct P4estMesh{NDIMS, RealT <: Real, IsParallel, P, Ghost, NDIMSP2, NN ghost_pw = PointerWrapper(ghost) mesh = new{NDIMS, eltype(tree_node_coordinates), typeof(is_parallel), - typeof(p4est_pw), typeof(ghost_pw), NDIMS + 2, length(nodes)}(p4est_pw, + typeof(p4est_pw), typeof(ghost_pw), NDIMS + 2, NDIMS + 3, length(nodes)}(p4est_pw, is_parallel, ghost_pw, tree_node_coordinates, @@ -61,7 +63,9 @@ mutable struct P4estMesh{NDIMS, RealT <: Real, IsParallel, P, Ghost, NDIMSP2, NN unsaved_changes, p4est_partition_allow_for_coarsening, mimetic, - exact_jacobian) + exact_jacobian, + polydeg_parent_metrics, + tree_contravariant_vectors) # Destroy `p4est` structs when the mesh is garbage collected finalizer(destroy_mesh, mesh) @@ -172,7 +176,8 @@ function P4estMesh(trees_per_dimension; polydeg, unsaved_changes = true, p4est_partition_allow_for_coarsening = true, mimetic = false, - exact_jacobian = false) + exact_jacobian = false, + polydeg_parent_metrics = 0) @assert ((coordinates_min === nothing)===(coordinates_max === nothing)) "Either both or none of coordinates_min and coordinates_max must be specified" @assert count(i -> i !== nothing, @@ -209,6 +214,47 @@ function P4estMesh(trees_per_dimension; polydeg, calc_tree_node_coordinates!(tree_node_coordinates, nodes, mapping, trees_per_dimension) + tree_contravariant_vectors = Array{RealT, NDIMS + 3}(undef, NDIMS, NDIMS, + ntuple(_ -> polydeg_parent_metrics + 1, + NDIMS)..., + prod(trees_per_dimension)) + + if polydeg_parent_metrics > 0 # Only available for one tree + basis_parent_metrics = LobattoLegendreBasis(RealT, polydeg_parent_metrics) + nodes_parent_metrics = basis_parent_metrics.nodes + + jacobian_matrix = zeros(3,3,ntuple(_ -> polydeg_parent_metrics+1,NDIMS)..., 1) + inverse_jacobian = zeros(ntuple(_ -> polydeg_parent_metrics+1,NDIMS)..., 1) + node_coordinates_comp = zeros(3, polydeg_parent_metrics+1) + node_coordinates_comp[1,:] = nodes_parent_metrics + node_coordinates_comp[2,:] = nodes_parent_metrics + node_coordinates_comp[3,:] = nodes_parent_metrics + + if exact_jacobian + calc_jacobian_matrix_exact!(jacobian_matrix, 1, tree_node_coordinates, basis_parent_metrics, node_coordinates_comp) + else + # TODO: This makes sense only for polydeg_geo = polydeg_parent_metrics + calc_jacobian_matrix!(jacobian_matrix, 1, tree_node_coordinates, basis_parent_metrics) + end + + if mimetic + calc_contravariant_vectors_mimetic!(tree_contravariant_vectors, 1, jacobian_matrix, + tree_node_coordinates, basis_parent_metrics, node_coordinates_comp) + else + # TODO: This makes sense only for polydeg_geo = polydeg_parent_metrics + calc_contravariant_vectors!(tree_contravariant_vectors, 1, jacobian_matrix, + tree_node_coordinates, basis_parent_metrics) + end + + calc_inverse_jacobian!(inverse_jacobian, 1, jacobian_matrix, basis_parent_metrics) + + #= @turbo for k in eachnode(basis_parent_metrics), j in eachnode(basis_parent_metrics), i in eachnode(basis_parent_metrics) + for r in 1:3, s in 1:3 + tree_contravariant_vectors[s,r,i,j,k,1] *= inverse_jacobian[i,j,k,1] + end + end =# + end + # p4est_connectivity_new_brick has trees in Z-order, so use our own function for this connectivity = connectivity_structured(trees_per_dimension..., periodicity) @@ -221,7 +267,7 @@ function P4estMesh(trees_per_dimension; polydeg, return P4estMesh{NDIMS}(p4est, tree_node_coordinates, nodes, boundary_names, "", unsaved_changes, - p4est_partition_allow_for_coarsening, mimetic, exact_jacobian) + p4est_partition_allow_for_coarsening, mimetic, exact_jacobian, polydeg_parent_metrics, tree_contravariant_vectors) end # 2D version diff --git a/src/solvers/dgsem_p4est/containers_3d.jl b/src/solvers/dgsem_p4est/containers_3d.jl index 5fb984b82bb..a2dda9edfa3 100644 --- a/src/solvers/dgsem_p4est/containers_3d.jl +++ b/src/solvers/dgsem_p4est/containers_3d.jl @@ -26,10 +26,10 @@ function init_elements!(elements, mesh::P4estMesh{3}, basis::LobattoLegendreBasi for i in eachindex(quadrants) element = offset + i + quad = quadrants[i] + quad_length = p4est_quadrant_len(quad.level) / p4est_root_len if mesh.exact_jacobian || mesh.mimetic - quad = quadrants[i] - quad_length = p4est_quadrant_len(quad.level) / p4est_root_len calc_node_coordinates_computational!(node_coordinates_comp, quad_length, p4est_root_len, quad, mesh, basis) end @@ -39,15 +39,19 @@ function init_elements!(elements, mesh::P4estMesh{3}, basis::LobattoLegendreBasi calc_jacobian_matrix!(jacobian_matrix, element, node_coordinates, basis) end - if mesh.mimetic - calc_contravariant_vectors_mimetic!(contravariant_vectors, element, jacobian_matrix, - node_coordinates, basis, node_coordinates_comp) - else - calc_contravariant_vectors!(contravariant_vectors, element, jacobian_matrix, - node_coordinates, basis) - end - calc_inverse_jacobian!(inverse_jacobian, element, jacobian_matrix, basis) + + if mesh.polydeg_parent_metrics > 0 + calc_contravariant_vectors_interpolate!(contravariant_vectors, element, quad, quad_length, tree, mesh, p4est_root_len, basis, inverse_jacobian) + else + if mesh.mimetic + calc_contravariant_vectors_mimetic!(contravariant_vectors, element, jacobian_matrix, + node_coordinates, basis, node_coordinates_comp) + else + calc_contravariant_vectors!(contravariant_vectors, element, jacobian_matrix, + node_coordinates, basis) + end + end end end @@ -116,6 +120,44 @@ function calc_node_coordinates_computational!(node_coordinates_comp, quad_length quad.z / p4est_root_len) .- 1 end +function calc_contravariant_vectors_interpolate!(contravariant_vectors, element, quad, quad_length, tree, mesh::P4estMesh{3}, p4est_root_len, basis, inverse_jacobian) + @unpack nodes = basis + + basis_parent_metrics = LobattoLegendreBasis(eltype(contravariant_vectors), mesh.polydeg_parent_metrics) + nodes_parent_metrics = basis_parent_metrics.nodes + + nodes_out_x = 2 * (quad_length * 1 / 2 * (nodes .+ 1) .+ + quad.x / p4est_root_len) .- 1 + nodes_out_y = 2 * (quad_length * 1 / 2 * (nodes .+ 1) .+ + quad.y / p4est_root_len) .- 1 + nodes_out_z = 2 * (quad_length * 1 / 2 * (nodes .+ 1) .+ + quad.z / p4est_root_len) .- 1 + + matrix1 = polynomial_interpolation_matrix(nodes_parent_metrics, nodes_out_x) + matrix2 = polynomial_interpolation_matrix(nodes_parent_metrics, nodes_out_y) + matrix3 = polynomial_interpolation_matrix(nodes_parent_metrics, nodes_out_z) + + multiply_dimensionwise!(view(contravariant_vectors, 1, :, :, :, :, element), + matrix1, matrix2, matrix3, + view(mesh.tree_contravariant_vectors, 1, :, :, :, :, tree)) + + multiply_dimensionwise!(view(contravariant_vectors, 2, :, :, :, :, element), + matrix1, matrix2, matrix3, + view(mesh.tree_contravariant_vectors, 2, :, :, :, :, tree)) + + multiply_dimensionwise!(view(contravariant_vectors, 3, :, :, :, :, element), + matrix1, matrix2, matrix3, + view(mesh.tree_contravariant_vectors, 3, :, :, :, :, tree)) + + #= @turbo for k in eachnode(basis), j in eachnode(basis), i in eachnode(basis) + for r in 1:3, s in 1:3 + contravariant_vectors[s,r,i,j,k,element] /= inverse_jacobian[i,j,k,element] + end + end =# + + contravariant_vectors[:,:,:,:,:,element] /= 2^(2 * quad.level) +end + # Initialize node_indices of interface container @inline function init_interface_node_indices!(interfaces::P4estInterfaceContainer{3}, faces, orientation, interface_id) From 24beec3dbe6f20790a9175e3f4974c5d8956beab Mon Sep 17 00:00:00 2001 From: Daniel Bach Date: Mon, 16 Sep 2024 14:58:50 +0200 Subject: [PATCH 15/19] Added the option to run multiple element simulations with the mimetic metrics --- Project.toml | 4 +- ...r_advection_free_stream_mimetic_metrics.jl | 121 +++++++++++++----- ...n_free_stream_mimetic_metrics_dirichlet.jl | 4 +- 3 files changed, 93 insertions(+), 36 deletions(-) diff --git a/Project.toml b/Project.toml index 828f4778f74..2f00e66c06f 100644 --- a/Project.toml +++ b/Project.toml @@ -19,7 +19,9 @@ MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" Octavian = "6fd5a793-0b7e-452c-907f-f8bfe9c57db4" OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" +OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" P4est = "7d669430-f675-4ae7-b43e-fab78ec5a902" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" @@ -59,8 +61,8 @@ HDF5 = "0.14, 0.15, 0.16" IfElse = "0.1" LinearMaps = "2.7, 3.0" LoopVectorization = "0.12.118" -Makie = "0.19" MPI = "0.20" +Makie = "0.19" MuladdMacro = "0.2.2" Octavian = "0.3.5" OffsetArrays = "1.3" diff --git a/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl b/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl index 528eb4c6f90..caa44644487 100644 --- a/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl +++ b/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl @@ -32,48 +32,62 @@ function mapping(xi, eta, zeta) return SVector(x, y, z) end -function exact_contravariant_vectors!(Ja, xi, eta, zeta) +function exact_contravariant_vectors!(Ja, xi, eta, zeta, dxi, deta, dzeta) theta_xi = theta_der1(xi, eta, zeta) theta_eta = theta_der2(xi, eta, zeta) theta_zeta = theta_der3(xi, eta, zeta) - Ja[1,1] = 1 + theta_eta + theta_zeta - Ja[1,2] = -theta_xi - Ja[1,3] = -theta_xi - Ja[2,1] = -theta_eta - Ja[2,2] = 1 + theta_xi + theta_zeta - Ja[2,3] = -theta_eta - Ja[3,1] = -theta_zeta - Ja[3,2] = -theta_zeta - Ja[3,3] = 1 + theta_xi + theta_eta + Ja[1,1] = deta*dzeta*(1 + theta_eta + theta_zeta) + Ja[1,2] = dxi*dzeta*(-theta_xi) + Ja[1,3] = dxi*deta*(-theta_xi) + Ja[2,1] = deta*dzeta*(-theta_eta) + Ja[2,2] = dxi*dzeta*(1 + theta_xi + theta_zeta) + Ja[2,3] = dxi*deta*(-theta_eta) + Ja[3,1] = deta*dzeta*(-theta_zeta) + Ja[3,2] = dxi*dzeta*(-theta_zeta) + Ja[3,3] = dxi*deta*(1 + theta_xi + theta_eta) end -function compute_error(solver, semi) +function compute_error(solver, semi, cells_per_dimension) @unpack nodes, weights = solver.basis exact_Ja = zero(MMatrix{3, 3, Float64}) error = zero(Float64) error_L2 = zero(Float64) - for k in eachnode(solver.basis) - for j in eachnode(solver.basis) - for i in eachnode(solver.basis) - exact_contravariant_vectors!(exact_Ja, nodes[i], nodes[j], nodes[k]) - error = max(error, maximum(abs.(semi.cache.elements.contravariant_vectors[:,:,i,j,k,1] - exact_Ja))) - error_L2 += norm(semi.cache.elements.contravariant_vectors[:,1,i,j,k,1] - exact_Ja[:,1]) * weights[i] * weights[j] * weights[k] + linear_indices = LinearIndices(size(semi.mesh)) + xi_scale = 1/cells_per_dimension[1] + eta_scale = 1/cells_per_dimension[2] + zeta_scale = 1/cells_per_dimension[3] + + for d3 in 1:cells_per_dimension[3] + for d2 in 1:cells_per_dimension[2] + for d1 in 1:cells_per_dimension[1] + node_coordinates_comp = zeros(3, nnodes(solver.basis)) + Trixi.calc_node_coordinates_computational!(node_coordinates_comp, d1, d2, d3, semi.mesh, solver.basis) + element = linear_indices[d1,d2,d3] + for k in eachnode(solver.basis) + for j in eachnode(solver.basis) + for i in eachnode(solver.basis) + exact_contravariant_vectors!(exact_Ja, node_coordinates_comp[1,i], node_coordinates_comp[2,j], node_coordinates_comp[3,k], xi_scale, eta_scale, zeta_scale) + error = max(error, maximum(abs.(semi.cache.elements.contravariant_vectors[:,:,i,j,k,element] - exact_Ja))) + error_L2 += norm(semi.cache.elements.contravariant_vectors[:,:,i,j,k,element] - exact_Ja) * weights[i] * weights[j] * weights[k] + end + end + end end end end - return error, error_L2 / 8 + return error, error_L2 / (8 * prod(cells_per_dimension)) end -cells_per_dimension = (1,1,1) +cells_per_dimension = (2,2,2) max_polydeg = 25 -errors_normals_inf = zeros(max_polydeg,2) -errors_normals_L2 = zeros(max_polydeg,2) -errors_sol_inf = zeros(max_polydeg,2) -errors_sol_L2 = zeros(max_polydeg,2) +errors_normals_inf = zeros(max_polydeg,3) +errors_normals_L2 = zeros(max_polydeg,3) +errors_sol_inf = zeros(max_polydeg,3) +errors_sol_L2 = zeros(max_polydeg,3) exact_jacobian = true -final_time = 3.00636132e-03 +final_time = 1e0 initial_condition = initial_condition_constant for polydeg in 1:max_polydeg println("Computing polydeg = ", polydeg) @@ -81,12 +95,12 @@ for polydeg in 1:max_polydeg solver = DGSEM(polydeg, flux_lax_friedrichs) # Create curved mesh with 8 x 8 x 8 elements - mesh = StructuredMesh(cells_per_dimension, mapping; mimetic = false, exact_jacobian = exact_jacobian) + mesh = StructuredMesh(cells_per_dimension, mapping; mimetic = false, exact_jacobian = false) # A semidiscre tization collects data structures and functions for the spatial discretization semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) - error_inf, error_L2 = compute_error(solver, semi) + error_inf, error_L2 = compute_error(solver, semi, cells_per_dimension) errors_normals_inf[polydeg,1] = error_inf errors_normals_L2[polydeg,1] = error_L2 @@ -117,12 +131,12 @@ for polydeg in 1:max_polydeg errors_sol_inf[polydeg, 1] = errors.linf[1] # Create curved mesh with 8 x 8 x 8 elements - mesh = StructuredMesh(cells_per_dimension, mapping; mimetic = true, exact_jacobian = exact_jacobian) + mesh = StructuredMesh(cells_per_dimension, mapping; mimetic = true, exact_jacobian = true) # A semidiscretization collects data structures and functions for the spatial discretization semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) - error_inf, error_L2 = compute_error(solver, semi) + error_inf, error_L2 = compute_error(solver, semi, cells_per_dimension) errors_normals_inf[polydeg,2] = error_inf errors_normals_L2[polydeg,2] = error_L2 @@ -151,13 +165,54 @@ for polydeg in 1:max_polydeg errors_sol_L2[polydeg, 2] = errors.l2[1] errors_sol_inf[polydeg, 2] = errors.linf[1] -end -#= plot(errors_normals_inf[:,1], yaxis=:log, label = "standard", linewidth=2, thickness_scaling = 1) -plot!(errors_normals_inf[:,2], yaxis=:log, label = "mimetic", linewidth=2, thickness_scaling = 1) =# + #= + # Create curved mesh with 8 x 8 x 8 elements + mesh = StructuredMesh(cells_per_dimension, mapping; mimetic = true, exact_jacobian = false) + + # A semidiscretization collects data structures and functions for the spatial discretization + semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + + error_inf, error_L2 = compute_error(solver, semi, cells_per_dimension) + errors_normals_inf[polydeg,3] = error_inf + errors_normals_L2[polydeg,3] = error_L2 + + # Create ODE problem with time span from 0.0 to 1.0 + ode = semidiscretize(semi, (0.0, final_time)); + + # The AnalysisCallback allows to analyse the solution in regular intervals and prints the results + analysis_callback = AnalysisCallback(semi, interval=100) + + # The StepsizeCallback handles the re-calculation of the maximum Δt after each time step + stepsize_callback = StepsizeCallback(cfl=0.1) -plot(3:max_polydeg,errors_sol_inf[3:end,1], xaxis=:log, yaxis=:log, label = "standard", linewidth=2, thickness_scaling = 1) -plot!(3:max_polydeg,errors_sol_inf[3:end,2], xaxis=:log, yaxis=:log, label = "mimetic", linewidth=2, thickness_scaling = 1) + # Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver + callbacks = CallbackSet(analysis_callback, stepsize_callback) + + + ############################################################################### + # run the simulation + + # OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks + sol = solve(ode, Euler(), #, CarpenterKennedy2N54(williamson_condition=false), + dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep=false, callback=callbacks); + + errors = analysis_callback(sol) + + errors_sol_L2[polydeg, 3] = errors.l2[1] + errors_sol_inf[polydeg, 3] = errors.linf[1] + =# +end +plot(errors_normals_L2[2:end,1], yaxis=:log, ylabel = "discrete L2 norm", xlabel = "polynomial degree", label = "standard", linewidth=2, thickness_scaling = 1) +plot!(errors_normals_L2[2:end,2], yaxis=:log, ylabel = "discrete L2 norm", xlabel = "polynomial degree", label = "mimetic", linewidth=2, thickness_scaling = 1) +#plot!(errors_normals_L2[2:end,3], yaxis=:log, label = "interpolated mapping mimetic", linewidth=2, thickness_scaling = 1) + +#= +plot(3:max_polydeg,errors_sol_inf[3:end,1], yaxis=:log, label = "standard", linewidth=2, thickness_scaling = 1) +plot!(3:max_polydeg,errors_sol_inf[3:end,2], yaxis=:log, label = "mimetic", linewidth=2, thickness_scaling = 1) +plot!(3:max_polydeg,errors_sol_inf[3:end,3], yaxis=:log, label = "interpolated mapping mimetic", linewidth=2, thickness_scaling = 1) +=# #= ############################################################################### # ODE solvers, callbacks etc. diff --git a/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics_dirichlet.jl b/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics_dirichlet.jl index 16e924e2db6..f2d49fcaaab 100644 --- a/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics_dirichlet.jl +++ b/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics_dirichlet.jl @@ -156,8 +156,8 @@ end plot!(errors_normals_inf[:,2], yaxis=:log, label = "mimetic", linewidth=2, thickness_scaling = 1) =# -plot(3:max_polydeg,errors_sol_inf[3:end,1], xaxis=:log, yaxis=:log, label = "standard", linewidth=2, thickness_scaling = 1) -plot!(3:max_polydeg,errors_sol_inf[3:end,2], xaxis=:log, yaxis=:log, label = "mimetic", linewidth=2, thickness_scaling = 1) +plot(3:max_polydeg,errors_normals_inf[3:end,1], xaxis=:log, yaxis=:log, label = "standard", linewidth=2, thickness_scaling = 1) +plot!(3:max_polydeg,errors_normals_inf[3:end,2], xaxis=:log, yaxis=:log, label = "mimetic", linewidth=2, thickness_scaling = 1) #= ############################################################################### # ODE solvers, callbacks etc. From 69f5c9ed15cbffef88dfdbe22a62371708ccf750 Mon Sep 17 00:00:00 2001 From: Daniel Bach Date: Wed, 18 Sep 2024 15:56:47 +0200 Subject: [PATCH 16/19] Sine wave test --- ...r_advection_free_stream_mimetic_metrics.jl | 41 ++++++++++--------- 1 file changed, 22 insertions(+), 19 deletions(-) diff --git a/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl b/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl index caa44644487..076aed6a97a 100644 --- a/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl +++ b/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl @@ -78,7 +78,9 @@ function compute_error(solver, semi, cells_per_dimension) return error, error_L2 / (8 * prod(cells_per_dimension)) end -cells_per_dimension = (2,2,2) +f(x,t,equations::LinearScalarAdvectionEquation3D) = SVector(sin(pi*(x[1]+x[2]+x[3]))) + +cells_per_dimension = (4,4,4) max_polydeg = 25 @@ -88,7 +90,8 @@ errors_sol_inf = zeros(max_polydeg,3) errors_sol_L2 = zeros(max_polydeg,3) exact_jacobian = true final_time = 1e0 -initial_condition = initial_condition_constant +initial_condition = f + for polydeg in 1:max_polydeg println("Computing polydeg = ", polydeg) # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux @@ -108,7 +111,7 @@ for polydeg in 1:max_polydeg ode = semidiscretize(semi, (0.0, final_time)); # The AnalysisCallback allows to analyse the solution in regular intervals and prints the results - analysis_callback = AnalysisCallback(semi, interval=100) + analysis_callback = AnalysisCallback(semi, interval=100, analysis_polydeg = 50) # The StepsizeCallback handles the re-calculation of the maximum Δt after each time step stepsize_callback = StepsizeCallback(cfl=0.1) @@ -125,11 +128,12 @@ for polydeg in 1:max_polydeg dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback save_everystep=false, callback=callbacks); + errors = analysis_callback(sol) errors_sol_L2[polydeg, 1] = errors.l2[1] errors_sol_inf[polydeg, 1] = errors.linf[1] - + #= # Create curved mesh with 8 x 8 x 8 elements mesh = StructuredMesh(cells_per_dimension, mapping; mimetic = true, exact_jacobian = true) @@ -165,8 +169,8 @@ for polydeg in 1:max_polydeg errors_sol_L2[polydeg, 2] = errors.l2[1] errors_sol_inf[polydeg, 2] = errors.linf[1] - - #= + =# + # Create curved mesh with 8 x 8 x 8 elements mesh = StructuredMesh(cells_per_dimension, mapping; mimetic = true, exact_jacobian = false) @@ -181,7 +185,7 @@ for polydeg in 1:max_polydeg ode = semidiscretize(semi, (0.0, final_time)); # The AnalysisCallback allows to analyse the solution in regular intervals and prints the results - analysis_callback = AnalysisCallback(semi, interval=100) + analysis_callback = AnalysisCallback(semi, interval=100, analysis_polydeg = 50) # The StepsizeCallback handles the re-calculation of the maximum Δt after each time step stepsize_callback = StepsizeCallback(cfl=0.1) @@ -200,19 +204,18 @@ for polydeg in 1:max_polydeg errors = analysis_callback(sol) - errors_sol_L2[polydeg, 3] = errors.l2[1] - errors_sol_inf[polydeg, 3] = errors.linf[1] - =# + errors_sol_L2[polydeg, 2] = errors.l2[1] + errors_sol_inf[polydeg, 2] = errors.linf[1] + end -plot(errors_normals_L2[2:end,1], yaxis=:log, ylabel = "discrete L2 norm", xlabel = "polynomial degree", label = "standard", linewidth=2, thickness_scaling = 1) -plot!(errors_normals_L2[2:end,2], yaxis=:log, ylabel = "discrete L2 norm", xlabel = "polynomial degree", label = "mimetic", linewidth=2, thickness_scaling = 1) -#plot!(errors_normals_L2[2:end,3], yaxis=:log, label = "interpolated mapping mimetic", linewidth=2, thickness_scaling = 1) - -#= -plot(3:max_polydeg,errors_sol_inf[3:end,1], yaxis=:log, label = "standard", linewidth=2, thickness_scaling = 1) -plot!(3:max_polydeg,errors_sol_inf[3:end,2], yaxis=:log, label = "mimetic", linewidth=2, thickness_scaling = 1) -plot!(3:max_polydeg,errors_sol_inf[3:end,3], yaxis=:log, label = "interpolated mapping mimetic", linewidth=2, thickness_scaling = 1) -=# +#plot(errors_normals_L2[2:end,1], yaxis=:log, ylabel = "discrete L2 norm", xlabel = "polynomial degree", label = "standard", linewidth=2, thickness_scaling = 1) +#plot!(errors_normals_L2[2:end,2], yaxis=:log, ylabel = "discrete L2 norm", xlabel = "polynomial degree", label = "mimetic", linewidth=2, thickness_scaling = 1) + + +plot(3:max_polydeg,errors_sol_L2[3:end,1], yaxis=:log, label = "standard", linewidth=2, thickness_scaling = 1) +plot!(3:max_polydeg,errors_sol_L2[3:end,2], yaxis=:log, label = "mimetic", linewidth=2, thickness_scaling = 1) + + #= ############################################################################### # ODE solvers, callbacks etc. From 4cfb43bbb95c7a88eb0b2986ba74f39bf03f9349 Mon Sep 17 00:00:00 2001 From: Daniel Bach Date: Wed, 18 Sep 2024 22:18:07 +0200 Subject: [PATCH 17/19] Implemented new error calculation for metric terms --- ...r_advection_free_stream_mimetic_metrics.jl | 31 ++++++++++++------- 1 file changed, 20 insertions(+), 11 deletions(-) diff --git a/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl b/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl index 076aed6a97a..dd67910ff21 100644 --- a/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl +++ b/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl @@ -57,18 +57,26 @@ function compute_error(solver, semi, cells_per_dimension) eta_scale = 1/cells_per_dimension[2] zeta_scale = 1/cells_per_dimension[3] + int_basis = Trixi.LobattoLegendreBasis(50) + int_nodes, int_weights = Trixi.gauss_lobatto_nodes_weights(51) + vandermonde = Trixi.polynomial_interpolation_matrix(nodes, int_nodes) + for d3 in 1:cells_per_dimension[3] for d2 in 1:cells_per_dimension[2] for d1 in 1:cells_per_dimension[1] - node_coordinates_comp = zeros(3, nnodes(solver.basis)) - Trixi.calc_node_coordinates_computational!(node_coordinates_comp, d1, d2, d3, semi.mesh, solver.basis) + node_coordinates_comp = zeros(3, nnodes(int_basis)) + Trixi.calc_node_coordinates_computational!(node_coordinates_comp, d1, d2, d3, semi.mesh, int_basis) element = linear_indices[d1,d2,d3] - for k in eachnode(solver.basis) - for j in eachnode(solver.basis) - for i in eachnode(solver.basis) + înterpolated_metric_values = zeros(3,3,51,51,51) + for j in 1:3 + Trixi.multiply_dimensionwise!(înterpolated_metric_values[j,:,:,:,:], vandermonde, semi.cache.elements.contravariant_vectors[j,:,:,:,:,element]) + end + for k in eachnode(int_basis) + for j in eachnode(int_basis) + for i in eachnode(int_basis) exact_contravariant_vectors!(exact_Ja, node_coordinates_comp[1,i], node_coordinates_comp[2,j], node_coordinates_comp[3,k], xi_scale, eta_scale, zeta_scale) - error = max(error, maximum(abs.(semi.cache.elements.contravariant_vectors[:,:,i,j,k,element] - exact_Ja))) - error_L2 += norm(semi.cache.elements.contravariant_vectors[:,:,i,j,k,element] - exact_Ja) * weights[i] * weights[j] * weights[k] + error = max(error, maximum(abs.(înterpolated_metric_values[:,:,i,j,k] - exact_Ja))) + error_L2 += norm(înterpolated_metric_values[:,:,i,j,k] - exact_Ja) * int_weights[i] * int_weights[j] * int_weights[k] end end end @@ -78,7 +86,7 @@ function compute_error(solver, semi, cells_per_dimension) return error, error_L2 / (8 * prod(cells_per_dimension)) end -f(x,t,equations::LinearScalarAdvectionEquation3D) = SVector(sin(pi*(x[1]+x[2]+x[3]))) +f(x,t,equations::LinearScalarAdvectionEquation3D) = SVector(sin(2*pi*(x[1]+x[2]+x[3]))) cells_per_dimension = (4,4,4) @@ -90,10 +98,11 @@ errors_sol_inf = zeros(max_polydeg,3) errors_sol_L2 = zeros(max_polydeg,3) exact_jacobian = true final_time = 1e0 -initial_condition = f +initial_condition = initial_condition_constant for polydeg in 1:max_polydeg println("Computing polydeg = ", polydeg) + #cells_per_dimension = (cld(50,polydeg),cld(50,polydeg),cld(50,polydeg)) # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux solver = DGSEM(polydeg, flux_lax_friedrichs) @@ -124,7 +133,7 @@ for polydeg in 1:max_polydeg # run the simulation # OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks - sol = solve(ode, Euler(), #CarpenterKennedy2N54(williamson_condition=false), + sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback save_everystep=false, callback=callbacks); @@ -198,7 +207,7 @@ for polydeg in 1:max_polydeg # run the simulation # OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks - sol = solve(ode, Euler(), #, CarpenterKennedy2N54(williamson_condition=false), + sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback save_everystep=false, callback=callbacks); From afabb2ff14a6f7cdb45600e693ad374f088be7fb Mon Sep 17 00:00:00 2001 From: Daniel Bach Date: Fri, 20 Sep 2024 11:37:44 +0200 Subject: [PATCH 18/19] fixed a bug in error norm calculations for metric terms --- ...r_advection_free_stream_mimetic_metrics.jl | 28 +++++++++---------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl b/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl index dd67910ff21..6d4df04e339 100644 --- a/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl +++ b/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl @@ -47,7 +47,7 @@ function exact_contravariant_vectors!(Ja, xi, eta, zeta, dxi, deta, dzeta) Ja[3,3] = dxi*deta*(1 + theta_xi + theta_eta) end -function compute_error(solver, semi, cells_per_dimension) +function compute_error(solver, semi, degree, cells_per_dimension) @unpack nodes, weights = solver.basis exact_Ja = zero(MMatrix{3, 3, Float64}) error = zero(Float64) @@ -57,8 +57,8 @@ function compute_error(solver, semi, cells_per_dimension) eta_scale = 1/cells_per_dimension[2] zeta_scale = 1/cells_per_dimension[3] - int_basis = Trixi.LobattoLegendreBasis(50) - int_nodes, int_weights = Trixi.gauss_lobatto_nodes_weights(51) + int_basis = Trixi.LobattoLegendreBasis(degree) + int_nodes, int_weights = Trixi.gauss_lobatto_nodes_weights(degree+1) vandermonde = Trixi.polynomial_interpolation_matrix(nodes, int_nodes) for d3 in 1:cells_per_dimension[3] @@ -67,16 +67,16 @@ function compute_error(solver, semi, cells_per_dimension) node_coordinates_comp = zeros(3, nnodes(int_basis)) Trixi.calc_node_coordinates_computational!(node_coordinates_comp, d1, d2, d3, semi.mesh, int_basis) element = linear_indices[d1,d2,d3] - înterpolated_metric_values = zeros(3,3,51,51,51) + interpolated_metric_values = zeros(3,3,degree+1,degree+1,degree+1) for j in 1:3 - Trixi.multiply_dimensionwise!(înterpolated_metric_values[j,:,:,:,:], vandermonde, semi.cache.elements.contravariant_vectors[j,:,:,:,:,element]) + Trixi.multiply_dimensionwise!(view(interpolated_metric_values,j,:,:,:,:), vandermonde, view(semi.cache.elements.contravariant_vectors,j,:,:,:,:,element)) end for k in eachnode(int_basis) for j in eachnode(int_basis) for i in eachnode(int_basis) exact_contravariant_vectors!(exact_Ja, node_coordinates_comp[1,i], node_coordinates_comp[2,j], node_coordinates_comp[3,k], xi_scale, eta_scale, zeta_scale) - error = max(error, maximum(abs.(înterpolated_metric_values[:,:,i,j,k] - exact_Ja))) - error_L2 += norm(înterpolated_metric_values[:,:,i,j,k] - exact_Ja) * int_weights[i] * int_weights[j] * int_weights[k] + error = max(error, maximum(abs.(interpolated_metric_values[:,:,i,j,k] - exact_Ja))) + error_L2 += norm(interpolated_metric_values[:,:,i,j,k] - exact_Ja) * int_weights[i] * int_weights[j] * int_weights[k] end end end @@ -86,9 +86,9 @@ function compute_error(solver, semi, cells_per_dimension) return error, error_L2 / (8 * prod(cells_per_dimension)) end -f(x,t,equations::LinearScalarAdvectionEquation3D) = SVector(sin(2*pi*(x[1]+x[2]+x[3]))) +f(x,t,equations::LinearScalarAdvectionEquation3D) = SVector(sin(pi*(x[1]+x[2]+x[3]))) -cells_per_dimension = (4,4,4) +cells_per_dimension = (2,2,2) max_polydeg = 25 @@ -102,7 +102,7 @@ initial_condition = initial_condition_constant for polydeg in 1:max_polydeg println("Computing polydeg = ", polydeg) - #cells_per_dimension = (cld(50,polydeg),cld(50,polydeg),cld(50,polydeg)) + # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux solver = DGSEM(polydeg, flux_lax_friedrichs) @@ -112,7 +112,7 @@ for polydeg in 1:max_polydeg # A semidiscre tization collects data structures and functions for the spatial discretization semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) - error_inf, error_L2 = compute_error(solver, semi, cells_per_dimension) + error_inf, error_L2 = compute_error(solver, semi, 50, cells_per_dimension) errors_normals_inf[polydeg,1] = error_inf errors_normals_L2[polydeg,1] = error_L2 @@ -186,9 +186,9 @@ for polydeg in 1:max_polydeg # A semidiscretization collects data structures and functions for the spatial discretization semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) - error_inf, error_L2 = compute_error(solver, semi, cells_per_dimension) - errors_normals_inf[polydeg,3] = error_inf - errors_normals_L2[polydeg,3] = error_L2 + error_inf, error_L2 = compute_error(solver, semi, 50, cells_per_dimension) + errors_normals_inf[polydeg,2] = error_inf + errors_normals_L2[polydeg,2] = error_L2 # Create ODE problem with time span from 0.0 to 1.0 ode = semidiscretize(semi, (0.0, final_time)); From 95603964561b2d1f9cfa83fcbd1f03935e7cc456 Mon Sep 17 00:00:00 2001 From: Daniel Bach Date: Wed, 9 Oct 2024 19:25:03 +0200 Subject: [PATCH 19/19] removed sine initial condition --- .../elixir_advection_free_stream_mimetic_metrics.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl b/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl index 6d4df04e339..830175fc197 100644 --- a/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl +++ b/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl @@ -86,8 +86,6 @@ function compute_error(solver, semi, degree, cells_per_dimension) return error, error_L2 / (8 * prod(cells_per_dimension)) end -f(x,t,equations::LinearScalarAdvectionEquation3D) = SVector(sin(pi*(x[1]+x[2]+x[3]))) - cells_per_dimension = (2,2,2) max_polydeg = 25