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/p4est_3d_dgsem/elixir_test_mimetic_metrics.jl b/examples/p4est_3d_dgsem/elixir_test_mimetic_metrics.jl new file mode 100644 index 00000000000..fdbd8d8baeb --- /dev/null +++ b/examples/p4est_3d_dgsem/elixir_test_mimetic_metrics.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 * 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 + + +exact_jacobian = false +mimetic = false +final_time = 1e-3 + +initial_condition = initial_condition_constant + +#= 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] + + 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/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/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..830175fc197 --- /dev/null +++ b/examples/structured_3d_dgsem/elixir_advection_free_stream_mimetic_metrics.jl @@ -0,0 +1,260 @@ +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, 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] = 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, degree, cells_per_dimension) + @unpack nodes, weights = solver.basis + exact_Ja = zero(MMatrix{3, 3, Float64}) + error = zero(Float64) + error_L2 = zero(Float64) + 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] + + 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] + for d2 in 1:cells_per_dimension[2] + for d1 in 1:cells_per_dimension[1] + 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] + interpolated_metric_values = zeros(3,3,degree+1,degree+1,degree+1) + for j in 1:3 + 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.(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 + end + end + end + return error, error_L2 / (8 * prod(cells_per_dimension)) +end + +cells_per_dimension = (2,2,2) + +max_polydeg = 25 + +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 = 1e0 +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 = 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, 50, cells_per_dimension) + 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, analysis_polydeg = 50) + + # 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, 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 = 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, 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)); + + # 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] + =# + + # 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, 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)); + + # The AnalysisCallback allows to analyse the solution in regular intervals and prints the results + 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) + + # 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, 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_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. + +# 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/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..f2d49fcaaab --- /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_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. + +# 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 66878f4b459..804730f3f46 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! @@ -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/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 60db285e04f..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 @@ -24,10 +24,14 @@ 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 + 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) 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 @@ -49,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, @@ -57,7 +61,11 @@ 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, + polydeg_parent_metrics, + tree_contravariant_vectors) # Destroy `p4est` structs when the mesh is garbage collected finalizer(destroy_mesh, mesh) @@ -166,7 +174,10 @@ 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, + 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, @@ -203,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) @@ -215,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) + p4est_partition_allow_for_coarsening, mimetic, exact_jacobian, polydeg_parent_metrics, tree_contravariant_vectors) end # 2D version diff --git a/src/meshes/structured_mesh.jl b/src/meshes/structured_mesh.jl index df067db833d..1824d8230bb 100644 --- a/src/meshes/structured_mesh.jl +++ b/src/meshes/structured_mesh.jl @@ -20,6 +20,8 @@ mutable struct StructuredMesh{NDIMS, RealT <: Real} <: AbstractMesh{NDIMS} periodicity::NTuple{NDIMS, Bool} current_filename::String unsaved_changes::Bool + mimetic::Bool + exact_jacobian::Bool end """ @@ -44,6 +46,8 @@ 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, + exact_jacobian = false, periodicity = true, unsaved_changes = true, mapping_as_string = mapping2string(mapping, length(cells_per_dimension))) @@ -63,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) + unsaved_changes, mimetic, exact_jacobian) end """ diff --git a/src/solvers/dgsem_p4est/containers_3d.jl b/src/solvers/dgsem_p4est/containers_3d.jl index e9994fe4569..a2dda9edfa3 100644 --- a/src/solvers/dgsem_p4est/containers_3d.jl +++ b/src/solvers/dgsem_p4est/containers_3d.jl @@ -12,13 +12,47 @@ 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) + + 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) + + 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 + 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 - calc_contravariant_vectors!(contravariant_vectors, element, jacobian_matrix, - node_coordinates, basis) + calc_inverse_jacobian!(inverse_jacobian, element, jacobian_matrix, basis) - 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 return nothing @@ -30,7 +64,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 +109,55 @@ 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::P4estMesh{3}, 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 + +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) diff --git a/src/solvers/dgsem_structured/containers_3d.jl b/src/solvers/dgsem_structured/containers_3d.jl index e843e869bf5..9b6957f8323 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,10 +21,23 @@ 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 || 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!(jacobian_matrix, element, node_coordinates, basis) + end - 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, 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 @@ -61,6 +76,57 @@ 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, 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] + # # ∂/∂η + # 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 + @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, 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, 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]) + 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, node_coordinates, basis) @@ -242,6 +308,293 @@ function calc_contravariant_vectors!(contravariant_vectors::AbstractArray{<:Any, 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 +""" +function calc_contravariant_vectors_mimetic!(contravariant_vectors::AbstractArray{<:Any, 6}, + element, + jacobian_matrix, node_coordinates, + 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)}) + 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 + + # Project the mapping potential \vec{g} \in H_{curl} to \vec{G} \in V_1 + 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) + 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[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[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 + + # Evaluate the mapping potential at the Lagrange points + 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) + 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, 3, 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 + + return contravariant_vectors +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)