diff --git a/Project.toml b/Project.toml index 89052859e42..5e1a60b7723 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.7.7-pre" +version = "0.7.8-pre" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl index 5bfa21c5322..1b485913ab2 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl @@ -110,7 +110,7 @@ semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabol # ODE solvers # Run for a long time to reach a state where forces stabilize up to 3 digits -tspan = (0.0, 1.0) +tspan = (0.0, 10.0) ode = semidiscretize(semi, tspan) # Callbacks @@ -130,12 +130,26 @@ lift_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names, u_inf(equations), l_inf())) +drag_coefficient_shear_force = AnalysisSurfaceIntegral(semi, force_boundary_names, + DragCoefficientShearStress(aoa(), + rho_inf(), + u_inf(equations), + l_inf())) + +lift_coefficient_shear_force = AnalysisSurfaceIntegral(semi, force_boundary_names, + LiftCoefficientShearStress(aoa(), + rho_inf(), + u_inf(equations), + l_inf())) + analysis_callback = AnalysisCallback(semi, interval = analysis_interval, output_directory = "out", save_analysis = true, analysis_errors = Symbol[], analysis_integrals = (drag_coefficient, - lift_coefficient)) + lift_coefficient, + drag_coefficient_shear_force, + lift_coefficient_shear_force)) alive_callback = AliveCallback(analysis_interval = analysis_interval) diff --git a/examples/structured_2d_dgsem/elixir_shallowwater_source_terms.jl b/examples/structured_2d_dgsem/elixir_shallowwater_source_terms.jl index 48fe37b9996..532fe8dbe7d 100644 --- a/examples/structured_2d_dgsem/elixir_shallowwater_source_terms.jl +++ b/examples/structured_2d_dgsem/elixir_shallowwater_source_terms.jl @@ -42,7 +42,7 @@ save_solution = SaveSolutionCallback(interval = 100, save_final_solution = true, solution_variables = cons2prim) -stepsize_callback = StepsizeCallback(cfl = 2.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, diff --git a/examples/structured_2d_dgsem/elixir_shallowwater_well_balanced.jl b/examples/structured_2d_dgsem/elixir_shallowwater_well_balanced.jl index a6a56aa807c..09abdf33843 100644 --- a/examples/structured_2d_dgsem/elixir_shallowwater_well_balanced.jl +++ b/examples/structured_2d_dgsem/elixir_shallowwater_well_balanced.jl @@ -108,7 +108,7 @@ save_solution = SaveSolutionCallback(interval = 1000, save_initial_solution = true, save_final_solution = true) -stepsize_callback = StepsizeCallback(cfl = 3.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl b/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl index a3df37fb966..af0da5d1768 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl @@ -77,7 +77,7 @@ save_solution = SaveSolutionCallback(interval = 100, save_initial_solution = true, save_final_solution = true) -stepsize_callback = StepsizeCallback(cfl = 3.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_well_balanced.jl b/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_well_balanced.jl index d9f1a52b500..a4f4b0189ba 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_well_balanced.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_well_balanced.jl @@ -70,7 +70,7 @@ save_solution = SaveSolutionCallback(interval = 1000, save_initial_solution = true, save_final_solution = true) -stepsize_callback = StepsizeCallback(cfl = 3.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl index 649e5023f6d..5851530e230 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl @@ -73,7 +73,7 @@ save_solution = SaveSolutionCallback(interval = 1000, save_initial_solution = true, save_final_solution = true) -stepsize_callback = StepsizeCallback(cfl = 3.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback) diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_ec.jl b/examples/tree_2d_dgsem/elixir_shallowwater_ec.jl index bc528ae7756..8221dfebe39 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_ec.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_ec.jl @@ -105,7 +105,7 @@ save_solution = SaveSolutionCallback(dt = 0.2, save_initial_solution = true, save_final_solution = true) -stepsize_callback = StepsizeCallback(cfl = 3.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback) diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced.jl b/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced.jl index 13023dfaba2..22043392b2a 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced.jl @@ -105,7 +105,7 @@ save_solution = SaveSolutionCallback(interval = 1000, save_initial_solution = true, save_final_solution = true) -stepsize_callback = StepsizeCallback(cfl = 3.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback) diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wall.jl b/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wall.jl index f50bd4e4f65..19073b0504a 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wall.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wall.jl @@ -108,7 +108,7 @@ save_solution = SaveSolutionCallback(interval = 1000, save_initial_solution = true, save_final_solution = true) -stepsize_callback = StepsizeCallback(cfl = 3.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback) diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_ec.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_ec.jl index 9122fb8287d..1f4aa414905 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_ec.jl +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_ec.jl @@ -105,7 +105,7 @@ save_solution = SaveSolutionCallback(interval = 1000, save_initial_solution = true, save_final_solution = true) -stepsize_callback = StepsizeCallback(cfl = 3.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback) diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_ec_shockcapturing.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_ec_shockcapturing.jl index 98408db5a78..c4736e8b9a5 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_ec_shockcapturing.jl +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_ec_shockcapturing.jl @@ -110,7 +110,7 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval) alive_callback = AliveCallback(analysis_interval = analysis_interval) -stepsize_callback = StepsizeCallback(cfl = 3.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, stepsize_callback) diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_well_balanced.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_well_balanced.jl index 6bad3a77f03..6cefca853c1 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_well_balanced.jl +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_well_balanced.jl @@ -104,7 +104,7 @@ save_solution = SaveSolutionCallback(interval = 1000, save_initial_solution = true, save_final_solution = true) -stepsize_callback = StepsizeCallback(cfl = 3.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback) diff --git a/src/Trixi.jl b/src/Trixi.jl index 8c4a987ebd5..b127b42f7f5 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -262,7 +262,8 @@ export SummaryCallback, SteadyStateCallback, AnalysisCallback, AliveCallback, AMRCallback, StepsizeCallback, GlmSpeedCallback, LBMCollisionCallback, EulerAcousticsCouplingCallback, TrivialCallback, AnalysisCallbackCoupled, - AnalysisSurfaceIntegral, DragCoefficientPressure, LiftCoefficientPressure + AnalysisSurfaceIntegral, DragCoefficientPressure, LiftCoefficientPressure, + DragCoefficientShearStress, LiftCoefficientShearStress export load_mesh, load_time, load_timestep, load_timestep!, load_dt, load_adaptive_time_integrator! diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl index 62048853b53..8f89af755a2 100644 --- a/src/callbacks_step/analysis.jl +++ b/src/callbacks_step/analysis.jl @@ -9,11 +9,11 @@ # - analysis_interval part as PeriodicCallback called after a certain amount of simulation time """ AnalysisCallback(semi; interval=0, - save_analysis=false, - output_directory="out", - analysis_filename="analysis.dat", - extra_analysis_errors=Symbol[], - extra_analysis_integrals=()) + save_analysis=false, + output_directory="out", + analysis_filename="analysis.dat", + extra_analysis_errors=Symbol[], + extra_analysis_integrals=()) Analyze a numerical solution every `interval` time steps and print the results to the screen. If `save_analysis`, the results are also saved in @@ -634,9 +634,7 @@ pretty_form_utf(quantity) = get_name(quantity) pretty_form_ascii(quantity) = get_name(quantity) # Special analyze for `SemidiscretizationHyperbolicParabolic` such that -# precomputed gradients are available. For now only implemented for the `enstrophy` -#!!! warning "Experimental code" -# This code is experimental and may be changed or removed in any future release. +# precomputed gradients are available. function analyze(quantity::typeof(enstrophy), du, u, t, semi::SemidiscretizationHyperbolicParabolic) mesh, equations, solver, cache = mesh_equations_solver_cache(semi) @@ -695,3 +693,19 @@ include("analysis_surface_integral_2d.jl") include("analysis_dg2d_parallel.jl") include("analysis_dg3d.jl") include("analysis_dg3d_parallel.jl") + +# Special analyze for `SemidiscretizationHyperbolicParabolic` such that +# precomputed gradients are available. Required for `enstrophy` (see above) and viscous forces. +# Note that this needs to be included after `analysis_surface_integral_2d.jl` to +# have `VariableViscous` available. +function analyze(quantity::AnalysisSurfaceIntegral{Variable}, + du, u, t, + semi::SemidiscretizationHyperbolicParabolic) where { + Variable <: + VariableViscous} + mesh, equations, solver, cache = mesh_equations_solver_cache(semi) + equations_parabolic = semi.equations_parabolic + cache_parabolic = semi.cache_parabolic + analyze(quantity, du, u, t, mesh, equations, equations_parabolic, solver, cache, + cache_parabolic) +end diff --git a/src/callbacks_step/analysis_surface_integral_2d.jl b/src/callbacks_step/analysis_surface_integral_2d.jl index 902fb0b4e85..c5a3d885eb0 100644 --- a/src/callbacks_step/analysis_surface_integral_2d.jl +++ b/src/callbacks_step/analysis_surface_integral_2d.jl @@ -13,15 +13,19 @@ boundary_symbol_or_boundary_symbols, variable) - This struct is used to compute the surface integral of a quantity of interest `variable` alongside - the boundary/boundaries associated with particular name(s) given in `boundary_symbol` - or `boundary_symbols`. - For instance, this can be used to compute the lift [`LiftCoefficientPressure`](@ref) or - drag coefficient [`DragCoefficientPressure`](@ref) of e.g. an airfoil with the boundary - name `:Airfoil` in 2D. +This struct is used to compute the surface integral of a quantity of interest `variable` alongside +the boundary/boundaries associated with particular name(s) given in `boundary_symbol` +or `boundary_symbols`. +For instance, this can be used to compute the lift [`LiftCoefficientPressure`](@ref) or +drag coefficient [`DragCoefficientPressure`](@ref) of e.g. an airfoil with the boundary +name `:Airfoil` in 2D. + +- `semi::Semidiscretization`: Passed in to retrieve boundary condition information +- `boundary_symbol_or_boundary_symbols::Symbol|Vector{Symbol}`: Name(s) of the boundary/boundaries + where the quantity of interest is computed +- `variable::Variable`: Quantity of interest, like lift or drag """ -struct AnalysisSurfaceIntegral{Semidiscretization, Variable} - semi::Semidiscretization # passed in to retrieve boundary condition information +struct AnalysisSurfaceIntegral{Variable} indices::Vector{Int} # Indices in `boundary_condition_indices` where quantity of interest is computed variable::Variable # Quantity of interest, like lift or drag @@ -29,8 +33,7 @@ struct AnalysisSurfaceIntegral{Semidiscretization, Variable} @unpack boundary_symbol_indices = semi.boundary_conditions indices = boundary_symbol_indices[boundary_symbol] - return new{typeof(semi), typeof(variable)}(semi, indices, - variable) + return new{typeof(variable)}(indices, variable) end function AnalysisSurfaceIntegral(semi, boundary_symbols::Vector{Symbol}, variable) @@ -41,8 +44,7 @@ struct AnalysisSurfaceIntegral{Semidiscretization, Variable} end sort!(indices) - return new{typeof(semi), typeof(variable)}(semi, indices, - variable) + return new{typeof(variable)}(indices, variable) end end @@ -50,7 +52,7 @@ struct ForceState{RealT <: Real} psi::Tuple{RealT, RealT} # Unit vector normal or parallel to freestream rhoinf::RealT uinf::RealT - l_inf::RealT + linf::RealT end struct LiftCoefficientPressure{RealT <: Real} @@ -61,13 +63,25 @@ struct DragCoefficientPressure{RealT <: Real} force_state::ForceState{RealT} end +# Abstract base type used for dispatch of `analyze` for quantities +# requiring gradients of the velocity field. +abstract type VariableViscous end + +struct LiftCoefficientShearStress{RealT <: Real} <: VariableViscous + force_state::ForceState{RealT} +end + +struct DragCoefficientShearStress{RealT <: Real} <: VariableViscous + force_state::ForceState{RealT} +end + """ - LiftCoefficientPressure(aoa, rhoinf, uinf, l_inf) + LiftCoefficientPressure(aoa, rhoinf, uinf, linf) Compute the lift coefficient ```math C_{L,p} \\coloneqq \\frac{\\oint_{\\partial \\Omega} p \\boldsymbol n \\cdot \\psi_L \\, \\mathrm{d} S} - {0.5 \\cdot \\rho_{\\infty} \\cdot U_{\\infty}^2 \\cdot L_{\\infty}} + {0.5 \\rho_{\\infty} U_{\\infty}^2 L_{\\infty}} ``` based on the pressure distribution along a boundary. Supposed to be used in conjunction with [`AnalysisSurfaceIntegral`](@ref) @@ -76,25 +90,25 @@ which stores the boundary information and semidiscretization. - `aoa::Real`: Angle of attack in radians (for airfoils etc.) - `rhoinf::Real`: Free-stream density - `uinf::Real`: Free-stream velocity -- `l_inf::Real`: Reference length of geometry (e.g. airfoil chord length) +- `linf::Real`: Reference length of geometry (e.g. airfoil chord length) """ -function LiftCoefficientPressure(aoa, rhoinf, uinf, l_inf) +function LiftCoefficientPressure(aoa, rhoinf, uinf, linf) # psi_lift is the normal unit vector to the freestream direction. # Note: The choice of the normal vector psi_lift = (-sin(aoa), cos(aoa)) # leads to positive lift coefficients for positive angles of attack for airfoils. # One could also use psi_lift = (sin(aoa), -cos(aoa)) which results in the same # value, but with the opposite sign. psi_lift = (-sin(aoa), cos(aoa)) - return LiftCoefficientPressure(ForceState(psi_lift, rhoinf, uinf, l_inf)) + return LiftCoefficientPressure(ForceState(psi_lift, rhoinf, uinf, linf)) end """ - DragCoefficientPressure(aoa, rhoinf, uinf, l_inf) + DragCoefficientPressure(aoa, rhoinf, uinf, linf) Compute the drag coefficient ```math C_{D,p} \\coloneqq \\frac{\\oint_{\\partial \\Omega} p \\boldsymbol n \\cdot \\psi_D \\, \\mathrm{d} S} - {0.5 \\cdot \\rho_{\\infty} \\cdot U_{\\infty}^2 \\cdot L_{\\infty}} + {0.5 \\rho_{\\infty} U_{\\infty}^2 L_{\\infty}} ``` based on the pressure distribution along a boundary. Supposed to be used in conjunction with [`AnalysisSurfaceIntegral`](@ref) @@ -103,28 +117,140 @@ which stores the boundary information and semidiscretization. - `aoa::Real`: Angle of attack in radians (for airfoils etc.) - `rhoinf::Real`: Free-stream density - `uinf::Real`: Free-stream velocity -- `l_inf::Real`: Reference length of geometry (e.g. airfoil chord length) +- `linf::Real`: Reference length of geometry (e.g. airfoil chord length) +""" +function DragCoefficientPressure(aoa, rhoinf, uinf, linf) + # `psi_drag` is the unit vector tangent to the freestream direction + psi_drag = (cos(aoa), sin(aoa)) + return DragCoefficientPressure(ForceState(psi_drag, rhoinf, uinf, linf)) +end + +""" + LiftCoefficientShearStress(aoa, rhoinf, uinf, linf) + +Compute the lift coefficient +```math +C_{L,f} \\coloneqq \\frac{\\oint_{\\partial \\Omega} \\boldsymbol \\tau_w \\cdot \\psi_L \\, \\mathrm{d} S} + {0.5 \\rho_{\\infty} U_{\\infty}^2 L_{\\infty}} +``` +based on the wall shear stress vector ``\\tau_w`` along a boundary. +Supposed to be used in conjunction with [`AnalysisSurfaceIntegral`](@ref) +which stores the boundary information and semidiscretization. + +- `aoa::Real`: Angle of attack in radians (for airfoils etc.) +- `rhoinf::Real`: Free-stream density +- `uinf::Real`: Free-stream velocity +- `linf::Real`: Reference length of geometry (e.g. airfoil chord length) +""" +function LiftCoefficientShearStress(aoa, rhoinf, uinf, linf) + # psi_lift is the normal unit vector to the freestream direction. + # Note: The choice of the normal vector psi_lift = (-sin(aoa), cos(aoa)) + # leads to negative lift coefficients for airfoils. + # One could also use psi_lift = (sin(aoa), -cos(aoa)) which results in the same + # value, but with the opposite sign. + psi_lift = (-sin(aoa), cos(aoa)) + return LiftCoefficientShearStress(ForceState(psi_lift, rhoinf, uinf, linf)) +end + +""" + DragCoefficientShearStress(aoa, rhoinf, uinf, linf) + +Compute the drag coefficient +```math +C_{D,f} \\coloneqq \\frac{\\oint_{\\partial \\Omega} \\boldsymbol \\tau_w \\cdot \\psi_D \\, \\mathrm{d} S} + {0.5 \\rho_{\\infty} U_{\\infty}^2 L_{\\infty}} +``` +based on the wall shear stress vector ``\\tau_w`` along a boundary. +Supposed to be used in conjunction with [`AnalysisSurfaceIntegral`](@ref) +which stores the boundary information and semidiscretization. + +- `aoa::Real`: Angle of attack in radians (for airfoils etc.) +- `rhoinf::Real`: Free-stream density +- `uinf::Real`: Free-stream velocity +- `linf::Real`: Reference length of geometry (e.g. airfoil chord length) """ -function DragCoefficientPressure(aoa, rhoinf, uinf, l_inf) +function DragCoefficientShearStress(aoa, rhoinf, uinf, linf) # `psi_drag` is the unit vector tangent to the freestream direction psi_drag = (cos(aoa), sin(aoa)) - return DragCoefficientPressure(ForceState(psi_drag, rhoinf, uinf, l_inf)) + return DragCoefficientShearStress(ForceState(psi_drag, rhoinf, uinf, linf)) end function (lift_coefficient::LiftCoefficientPressure)(u, normal_direction, equations) p = pressure(u, equations) - @unpack psi, rhoinf, uinf, l_inf = lift_coefficient.force_state + @unpack psi, rhoinf, uinf, linf = lift_coefficient.force_state # Normalize as `normal_direction` is not necessarily a unit vector n = dot(normal_direction, psi) / norm(normal_direction) - return p * n / (0.5 * rhoinf * uinf^2 * l_inf) + return p * n / (0.5 * rhoinf * uinf^2 * linf) end function (drag_coefficient::DragCoefficientPressure)(u, normal_direction, equations) p = pressure(u, equations) - @unpack psi, rhoinf, uinf, l_inf = drag_coefficient.force_state + @unpack psi, rhoinf, uinf, linf = drag_coefficient.force_state # Normalize as `normal_direction` is not necessarily a unit vector n = dot(normal_direction, psi) / norm(normal_direction) - return p * n / (0.5 * rhoinf * uinf^2 * l_inf) + return p * n / (0.5 * rhoinf * uinf^2 * linf) +end + +# Compute the three components of the symmetric viscous stress tensor +# (tau_11, tau_12, tau_22) based on the gradients of the velocity field. +# This is required for drag and lift coefficients based on shear stress, +# as well as for the non-integrated quantities such as +# skin friction coefficient (to be added). +function viscous_stress_tensor(u, normal_direction, equations_parabolic, + gradients_1, gradients_2) + _, dv1dx, dv2dx, _ = convert_derivative_to_primitive(u, gradients_1, + equations_parabolic) + _, dv1dy, dv2dy, _ = convert_derivative_to_primitive(u, gradients_2, + equations_parabolic) + + # Components of viscous stress tensor + # (4/3 * (v1)_x - 2/3 * (v2)_y) + tau_11 = 4.0 / 3.0 * dv1dx - 2.0 / 3.0 * dv2dy + # ((v1)_y + (v2)_x) + # stress tensor is symmetric + tau_12 = dv1dy + dv2dx # = tau_21 + # (4/3 * (v2)_y - 2/3 * (v1)_x) + tau_22 = 4.0 / 3.0 * dv2dy - 2.0 / 3.0 * dv1dx + + mu = dynamic_viscosity(u, equations_parabolic) + + return mu .* (tau_11, tau_12, tau_22) +end + +function viscous_stress_vector(u, normal_direction, equations_parabolic, + gradients_1, gradients_2) + # Normalize normal direction, should point *into* the fluid => *(-1) + n_normal = -normal_direction / norm(normal_direction) + + tau_11, tau_12, tau_22 = viscous_stress_tensor(u, normal_direction, + equations_parabolic, + gradients_1, gradients_2) + + # Viscous stress vector: Stress tensor * normal vector + visc_stress_vector_1 = tau_11 * n_normal[1] + tau_12 * n_normal[2] + visc_stress_vector_2 = tau_12 * n_normal[1] + tau_22 * n_normal[2] + + return (visc_stress_vector_1, visc_stress_vector_2) +end + +function (lift_coefficient::LiftCoefficientShearStress)(u, normal_direction, + equations_parabolic, + gradients_1, gradients_2) + visc_stress_vector = viscous_stress_vector(u, normal_direction, equations_parabolic, + gradients_1, gradients_2) + @unpack psi, rhoinf, uinf, linf = lift_coefficient.force_state + return (visc_stress_vector[1] * psi[1] + visc_stress_vector[2] * psi[2]) / + (0.5 * rhoinf * uinf^2 * linf) +end + +function (drag_coefficient::DragCoefficientShearStress)(u, normal_direction, + equations_parabolic, + gradients_1, gradients_2) + visc_stress_vector = viscous_stress_vector(u, normal_direction, equations_parabolic, + gradients_1, gradients_2) + @unpack psi, rhoinf, uinf, linf = drag_coefficient.force_state + return (visc_stress_vector[1] * psi[1] + visc_stress_vector[2] * psi[2]) / + (0.5 * rhoinf * uinf^2 * linf) end function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t, @@ -169,21 +295,88 @@ function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t, return surface_integral end -function pretty_form_ascii(::AnalysisSurfaceIntegral{<:Any, - <:LiftCoefficientPressure{<:Any}}) +function analyze(surface_variable::AnalysisSurfaceIntegral{Variable}, + du, u, t, mesh::P4estMesh{2}, + equations, equations_parabolic, + dg::DGSEM, cache, + cache_parabolic) where {Variable <: VariableViscous} + @unpack boundaries = cache + @unpack surface_flux_values, node_coordinates, contravariant_vectors = cache.elements + @unpack weights = dg.basis + + @unpack indices, variable = surface_variable + + # Additions for parabolic + @unpack viscous_container = cache_parabolic + @unpack gradients = viscous_container + + gradients_x, gradients_y = gradients + + surface_integral = zero(eltype(u)) + index_range = eachnode(dg) + for boundary in indices + element = boundaries.neighbor_ids[boundary] + node_indices = boundaries.node_indices[boundary] + direction = indices2direction(node_indices) + + i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], index_range) + j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], index_range) + + i_node = i_node_start + j_node = j_node_start + for node_index in index_range + u_node = Trixi.get_node_vars(cache.boundaries.u, equations, dg, node_index, + boundary) + # Extract normal direction at nodes which points from the elements outwards, + # i.e., *into* the structure. + normal_direction = get_normal_direction(direction, contravariant_vectors, + i_node, j_node, + element) + + # L2 norm of normal direction (contravariant_vector) is the surface element + dS = weights[node_index] * norm(normal_direction) + + gradients_1 = get_node_vars(gradients_x, equations_parabolic, dg, i_node, + j_node, element) + gradients_2 = get_node_vars(gradients_y, equations_parabolic, dg, i_node, + j_node, element) + + # Integral over whole boundary surface + surface_integral += variable(u_node, normal_direction, equations_parabolic, + gradients_1, gradients_2) * dS + + i_node += i_node_step + j_node += j_node_step + end + end + return surface_integral +end + +function pretty_form_ascii(::AnalysisSurfaceIntegral{<:LiftCoefficientPressure{<:Any}}) "CL_p" end -function pretty_form_utf(::AnalysisSurfaceIntegral{<:Any, - <:LiftCoefficientPressure{<:Any}}) +function pretty_form_utf(::AnalysisSurfaceIntegral{<:LiftCoefficientPressure{<:Any}}) "CL_p" end -function pretty_form_ascii(::AnalysisSurfaceIntegral{<:Any, - <:DragCoefficientPressure{<:Any}}) +function pretty_form_ascii(::AnalysisSurfaceIntegral{<:DragCoefficientPressure{<:Any}}) "CD_p" end -function pretty_form_utf(::AnalysisSurfaceIntegral{<:Any, - <:DragCoefficientPressure{<:Any}}) +function pretty_form_utf(::AnalysisSurfaceIntegral{<:DragCoefficientPressure{<:Any}}) "CD_p" end + +function pretty_form_ascii(::AnalysisSurfaceIntegral{<:LiftCoefficientShearStress{<:Any}}) + "CL_f" +end +function pretty_form_utf(::AnalysisSurfaceIntegral{<:LiftCoefficientShearStress{<:Any}}) + "CL_f" +end + +function pretty_form_ascii(::AnalysisSurfaceIntegral{<:DragCoefficientShearStress{<:Any}}) + "CD_f" +end +function pretty_form_utf(::AnalysisSurfaceIntegral{<:DragCoefficientShearStress{<:Any}}) + "CD_f" +end end # muladd diff --git a/src/equations/linearized_euler_1d.jl b/src/equations/linearized_euler_1d.jl index 22f452219b5..856a53a3d5d 100644 --- a/src/equations/linearized_euler_1d.jl +++ b/src/equations/linearized_euler_1d.jl @@ -24,7 +24,7 @@ Linearized euler equations in one space dimension. The equations are given by 0 \\ 0 \\ 0 \end{pmatrix} ``` -The bar ``\bar{(\cdot)}`` indicates uniform mean flow variables and c is the speed of sound. +The bar ``\bar{(\cdot)}`` indicates uniform mean flow variables and ``c`` is the speed of sound. The unknowns are the perturbation quantities of the acoustic velocity ``v_1'``, the pressure ``p'`` and the density ``\rho'``. """ diff --git a/src/equations/linearized_euler_2d.jl b/src/equations/linearized_euler_2d.jl index d497762bf62..d9d0d44e553 100644 --- a/src/equations/linearized_euler_2d.jl +++ b/src/equations/linearized_euler_2d.jl @@ -29,7 +29,7 @@ Linearized euler equations in two space dimensions. The equations are given by 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} ``` -The bar ``\bar{(\cdot)}`` indicates uniform mean flow variables and c is the speed of sound. +The bar ``\bar{(\cdot)}`` indicates uniform mean flow variables and ``c`` is the speed of sound. The unknowns are the acoustic velocities ``v' = (v_1', v_2')``, the pressure ``p'`` and the density ``\rho'``. """ struct LinearizedEulerEquations2D{RealT <: Real} <: diff --git a/src/equations/shallow_water_1d.jl b/src/equations/shallow_water_1d.jl index e348ef946b7..7007bea887a 100644 --- a/src/equations/shallow_water_1d.jl +++ b/src/equations/shallow_water_1d.jl @@ -551,7 +551,7 @@ end h = waterheight(u, equations) v = velocity(u, equations) - c = equations.gravity * sqrt(h) + c = sqrt(equations.gravity * h) return (abs(v) + c,) end diff --git a/src/equations/shallow_water_2d.jl b/src/equations/shallow_water_2d.jl index 74a299a51e6..73fae89a0fa 100644 --- a/src/equations/shallow_water_2d.jl +++ b/src/equations/shallow_water_2d.jl @@ -927,7 +927,7 @@ end h = waterheight(u, equations) v1, v2 = velocity(u, equations) - c = equations.gravity * sqrt(h) + c = sqrt(equations.gravity * h) return abs(v1) + c, abs(v2) + c end diff --git a/src/equations/shallow_water_quasi_1d.jl b/src/equations/shallow_water_quasi_1d.jl index 51c360104a7..08620021254 100644 --- a/src/equations/shallow_water_quasi_1d.jl +++ b/src/equations/shallow_water_quasi_1d.jl @@ -248,7 +248,7 @@ end h = waterheight(u, equations) v = velocity(u, equations) - c = equations.gravity * sqrt(h) + c = sqrt(equations.gravity * h) return (abs(v) + c,) end diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index edbdfe5a84f..a18dfb353b7 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -719,7 +719,10 @@ end 1.199362305026636, 0.9077214424040279, 5.666071182328691], tspan=(0.0, 0.001), - initial_refinement_level=0) + initial_refinement_level=0, + # With the default `maxiters = 1` in coverage tests, + # there would be no time steps after the restart. + coverage_override=(maxiters = 10_000,)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let @@ -728,6 +731,30 @@ end du_ode = similar(u_ode) @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end + + u_ode = copy(sol.u[end]) + du_ode = zero(u_ode) # Just a placeholder in this case + + u = Trixi.wrap_array(u_ode, semi) + du = Trixi.wrap_array(du_ode, semi) + + drag_p = Trixi.analyze(drag_coefficient, du, u, tspan[2], mesh, equations, solver, + semi.cache) + lift_p = Trixi.analyze(lift_coefficient, du, u, tspan[2], mesh, equations, solver, + semi.cache) + + drag_f = Trixi.analyze(drag_coefficient_shear_force, du, u, tspan[2], mesh, + equations, equations_parabolic, solver, + semi.cache, semi.cache_parabolic) + lift_f = Trixi.analyze(lift_coefficient_shear_force, du, u, tspan[2], mesh, + equations, equations_parabolic, solver, + semi.cache, semi.cache_parabolic) + + @test isapprox(drag_p, 0.17963843913309516, atol = 1e-13) + @test isapprox(lift_p, 0.26462588007949367, atol = 1e-13) + + @test isapprox(drag_f, 1.5427441885921553, atol = 1e-13) + @test isapprox(lift_f, 0.005621910087395724, atol = 1e-13) end end diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 64a1faf05b8..85eb681ec78 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -849,15 +849,15 @@ end @trixi_testset "elixir_shallowwater_source_terms.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), l2=[ - 0.0017285599436729316, - 0.025584610912606776, - 0.028373834961180594, + 0.0017286908591070864, + 0.025585037307655684, + 0.028374244567802766, 6.274146767730866e-5, ], linf=[ - 0.012972309788264802, - 0.108283714215621, - 0.15831585777928936, + 0.012973752001194772, + 0.10829375385832263, + 0.15832858475438094, 0.00018196759554722775, ], tspan=(0.0, 0.05)) diff --git a/test/test_tree_1d_shallowwater.jl b/test/test_tree_1d_shallowwater.jl index 41ad5c32bbd..8fe3291a938 100644 --- a/test/test_tree_1d_shallowwater.jl +++ b/test/test_tree_1d_shallowwater.jl @@ -12,10 +12,14 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_shallowwater_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec.jl"), - l2=[0.244729018751225, 0.8583565222389505, 0.07330427577586297], + l2=[ + 0.24476140682560343, + 0.8587309324660326, + 0.07330427577586297, + ], linf=[ - 2.1635021283528504, - 3.8717508164234453, + 2.1636963952308372, + 3.8737770522883115, 1.7711213427919539, ], tspan=(0.0, 0.25)) @@ -32,13 +36,13 @@ end @trixi_testset "elixir_shallowwater_ec.jl with initial_condition_weak_blast_wave" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec.jl"), l2=[ - 0.39464782107209717, - 2.03880864210846, + 0.39472828074570576, + 2.0390687947320076, 4.1623084150546725e-10, ], linf=[ - 0.778905801278281, - 3.2409883402608273, + 0.7793741954662221, + 3.2411927977882096, 7.419800190922032e-10, ], initial_condition=initial_condition_weak_blast_wave, diff --git a/test/test_tree_2d_shallowwater.jl b/test/test_tree_2d_shallowwater.jl index 01742644736..9a3ba36c7d4 100644 --- a/test/test_tree_2d_shallowwater.jl +++ b/test/test_tree_2d_shallowwater.jl @@ -13,15 +13,15 @@ EXAMPLES_DIR = joinpath(examples_dir(), "tree_2d_dgsem") @trixi_testset "elixir_shallowwater_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec.jl"), l2=[ - 0.991181203601035, - 0.734130029040644, - 0.7447696147162621, + 0.9911802019934329, + 0.7340106828033273, + 0.7446338002084801, 0.5875351036989047, ], linf=[ - 2.0117744577945413, - 2.9962317608172127, - 2.6554999727293653, + 2.0120253138457564, + 2.991158989293406, + 2.6557412817714035, 3.0, ], tspan=(0.0, 0.25)) @@ -280,15 +280,15 @@ end @trixi_testset "elixir_shallowwater_wall.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_wall.jl"), l2=[ - 0.13517233723296504, - 0.20010876311162215, - 0.20010876311162223, + 0.1351723240085936, + 0.20010881416550014, + 0.2001088141654999, 2.719538414346464e-7, ], linf=[ - 0.5303607982988336, - 0.5080989745682338, - 0.5080989745682352, + 0.5303608302490757, + 0.5080987791967457, + 0.5080987791967506, 1.1301675764130437e-6, ], tspan=(0.0, 0.25)) diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index 6814250dd47..35a380e4f24 100644 --- a/test/test_unstructured_2d.jl +++ b/test/test_unstructured_2d.jl @@ -301,15 +301,15 @@ end @trixi_testset "elixir_shallowwater_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec.jl"), l2=[ - 0.6106939484178353, - 0.48586236867426724, - 0.48234490854514356, - 0.29467422718511727, + 0.6107326269462766, + 0.48666631722018877, + 0.48309775159067053, + 0.29467422718511704, ], linf=[ - 2.775979948281604, - 3.1721242154451548, - 3.5713448319601393, + 2.776782342826098, + 3.2158378644333707, + 3.652920889487258, 2.052861364219655, ], tspan=(0.0, 0.25)) @@ -408,16 +408,16 @@ end @trixi_testset "elixir_shallowwater_source_terms.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), l2=[ - 0.0011197623982310795, - 0.04456344888447023, - 0.014317376629669337, - 5.089218476758975e-6, + 0.001118134082248467, + 0.044560486817464634, + 0.01430926600634214, + 5.089218476759981e-6, ], linf=[ - 0.007835284004819698, - 0.3486891284278597, - 0.11242778979399048, - 2.6407324614119432e-5, + 0.007798727223654822, + 0.34782952734839157, + 0.11161614702628064, + 2.6407324614341476e-5, ], tspan=(0.0, 0.025)) # Ensure that we do not have excessive memory allocations @@ -433,15 +433,15 @@ end @trixi_testset "elixir_shallowwater_source_terms.jl with FluxHydrostaticReconstruction" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), l2=[ - 0.001119678684752799, - 0.015429108794630785, - 0.01708275441241111, - 5.089218476758271e-6, + 0.0011196838135485918, + 0.01542895635133927, + 0.017082803023121197, + 5.089218476759981e-6, ], linf=[ - 0.014299564388827513, - 0.12785126473870534, - 0.17626788561725526, + 0.014299541415654371, + 0.12783948113206955, + 0.17626489583921323, 2.6407324614341476e-5, ], surface_flux=(FluxHydrostaticReconstruction(flux_hll, @@ -461,15 +461,15 @@ end @trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_ersing_etal" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), l2=[ - 0.0011196687776346434, - 0.044562672453443995, - 0.014306265289763618, + 0.001118046975499805, + 0.04455969246244461, + 0.014298120235633432, 5.089218476759981e-6, ], linf=[ - 0.007825021762002393, - 0.348550815397918, - 0.1115517935018282, + 0.007776521213640031, + 0.34768318303226353, + 0.11075311228066198, 2.6407324614341476e-5, ], surface_flux=(flux_wintermeyer_etal, @@ -490,15 +490,15 @@ end @trixi_testset "elixir_shallowwater_source_terms.jl with flux_hll" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), l2=[ - 0.0011196786847528799, - 0.015429108794631075, - 0.017082754412411742, + 0.0011196838135486059, + 0.015428956351339451, + 0.017082803023120943, 5.089218476759981e-6, ], linf=[ - 0.014299564388830177, - 0.12785126473870667, - 0.17626788561728546, + 0.01429954141565526, + 0.12783948113205668, + 0.176264895839215, 2.6407324614341476e-5, ], surface_flux=(flux_hll, @@ -561,15 +561,15 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec_shockcapturing.jl"), l2=[ - 0.6124656312639043, - 0.504371951785709, - 0.49180896200746366, - 0.29467422718511727, + 0.612551520607341, + 0.5039173660221961, + 0.49136517934903523, + 0.29467422718511704, ], linf=[ - 2.7639232436274392, - 3.3985508653311767, - 3.3330308209196224, + 2.7636771472622197, + 3.236168963021072, + 3.3363936775653826, 2.052861364219655, ], tspan=(0.0, 0.25))