Skip to content

Commit

Permalink
add x coordinates and time as input arguments to the variable functio…
Browse files Browse the repository at this point in the history 
…n to be integrated along a boundary (#1902)
  • Loading branch information
@andrewwinters5000
andrewwinters5000 committed Apr 10, 2024
1 parent efbd5a7 commit b1a84a6
Showing 1 changed file with 30 additions and 12 deletions.
42 changes: 30 additions & 12 deletions src/callbacks_step/analysis_surface_integral_2d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,7 @@ struct DragCoefficientPressure{RealT <: Real}
force_state::ForceState{RealT}
end

# Abstract base type used for dispatch of `analyze` for quantities
# Abstract base type used for dispatch of `analyze` for quantities
# requiring gradients of the velocity field.
abstract type VariableViscous end

Expand Down Expand Up @@ -175,15 +175,17 @@ function DragCoefficientShearStress(aoa, rhoinf, uinf, linf)
return DragCoefficientShearStress(ForceState(psi_drag, rhoinf, uinf, linf))
end

function (lift_coefficient::LiftCoefficientPressure)(u, normal_direction, equations)
function (lift_coefficient::LiftCoefficientPressure)(u, normal_direction, x, t,
equations)
p = pressure(u, equations)
@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 * linf)
end

function (drag_coefficient::DragCoefficientPressure)(u, normal_direction, equations)
function (drag_coefficient::DragCoefficientPressure)(u, normal_direction, x, t,
equations)
p = pressure(u, equations)
@unpack psi, rhoinf, uinf, linf = drag_coefficient.force_state
# Normalize as `normal_direction` is not necessarily a unit vector
Expand All @@ -193,8 +195,8 @@ 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
# 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)
Expand Down Expand Up @@ -233,7 +235,7 @@ function viscous_stress_vector(u, normal_direction, equations_parabolic,
return (visc_stress_vector_1, visc_stress_vector_2)
end

function (lift_coefficient::LiftCoefficientShearStress)(u, normal_direction,
function (lift_coefficient::LiftCoefficientShearStress)(u, normal_direction, x, t,
equations_parabolic,
gradients_1, gradients_2)
visc_stress_vector = viscous_stress_vector(u, normal_direction, equations_parabolic,
Expand All @@ -243,7 +245,7 @@ function (lift_coefficient::LiftCoefficientShearStress)(u, normal_direction,
(0.5 * rhoinf * uinf^2 * linf)
end

function (drag_coefficient::DragCoefficientShearStress)(u, normal_direction,
function (drag_coefficient::DragCoefficientShearStress)(u, normal_direction, x, t,
equations_parabolic,
gradients_1, gradients_2)
visc_stress_vector = viscous_stress_vector(u, normal_direction, equations_parabolic,
Expand Down Expand Up @@ -283,10 +285,18 @@ function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t,
i_node, j_node,
element)

# Coordinates at a boundary node
x = get_node_coords(node_coordinates, equations, dg, i_node, j_node,
element)

# L2 norm of normal direction (contravariant_vector) is the surface element
dS = weights[node_index] * norm(normal_direction)
# Integral over entire boundary surface
surface_integral += variable(u_node, normal_direction, equations) * dS

# Integral over entire boundary surface. Note, it is assumed that the
# `normal_direction` is normalized to be a normal vector within the
# function `variable` and the division of the normal scaling factor
# `norm(normal_direction)` is then accounted for with the `dS` quantity.
surface_integral += variable(u_node, normal_direction, x, t, equations) * dS

i_node += i_node_step
j_node += j_node_step
Expand Down Expand Up @@ -328,11 +338,15 @@ function analyze(surface_variable::AnalysisSurfaceIntegral{Variable},
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.
# i.e., *into* the structure.
normal_direction = get_normal_direction(direction, contravariant_vectors,
i_node, j_node,
element)

# Coordinates at a boundary node
x = get_node_coords(node_coordinates, equations, dg, i_node, j_node,
element)

# L2 norm of normal direction (contravariant_vector) is the surface element
dS = weights[node_index] * norm(normal_direction)

Expand All @@ -341,8 +355,12 @@ function analyze(surface_variable::AnalysisSurfaceIntegral{Variable},
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,
# Integral over whole boundary surface. Note, it is assumed that the
# `normal_direction` is normalized to be a normal vector within the
# function `variable` and the division of the normal scaling factor
# `norm(normal_direction)` is then accounted for with the `dS` quantity.
surface_integral += variable(u_node, normal_direction, x, t,
equations_parabolic,
gradients_1, gradients_2) * dS

i_node += i_node_step
Expand Down

0 comments on commit b1a84a6

Please sign in to comment.