Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix enstrophy computation for DGMulti #2143

Open
DanielDoehring opened this issue Nov 5, 2024 · 0 comments
Open

Fix enstrophy computation for DGMulti #2143

DanielDoehring opened this issue Nov 5, 2024 · 0 comments
Labels
bug Something isn't working enhancement New feature or request

Comments

@DanielDoehring
Copy link
Contributor

The original implementation of the enstrophy analysis callback introduced in PR1239 does not work and was not tested.

It is removed in #2140 .

A starting point for a fix would be (note that there are already couple changes to the original version)

function integrate(func::typeof(enstrophy), u,
                   mesh::DGMultiMesh,
                   equations, equations_parabolic::CompressibleNavierStokesDiffusion3D,
                   dg::DGMulti,
                   cache, cache_parabolic; normalize=true)
  gradients_x, gradients_y, gradients_z = cache_parabolic.gradients
  # allocate local storage for gradients.
  # TODO: can we avoid allocating here?
  local_gradient_quadrature_values = ntuple(_ -> copy(cache_parabolic.local_u_values_threaded), 3)
  integral = zero(eltype(eltype(u)))
  for e in eachelement(mesh, dg)
    u_quadrature_values = cache_parabolic.local_u_values_threaded[Threads.threadid()]
    gradient_x_quadrature_values = local_gradient_quadrature_values[1][Threads.threadid()]
    gradient_y_quadrature_values = local_gradient_quadrature_values[2][Threads.threadid()]
    gradient_z_quadrature_values = local_gradient_quadrature_values[3][Threads.threadid()]
    # interpolate to quadrature on each element
    apply_to_each_field(mul_by(dg.basis.Vq), u_quadrature_values, view(u, :, e))
    apply_to_each_field(mul_by(dg.basis.Vq), gradient_x_quadrature_values, view(gradients_x, :, e))
    apply_to_each_field(mul_by(dg.basis.Vq), gradient_y_quadrature_values, view(gradients_y, :, e))
    apply_to_each_field(mul_by(dg.basis.Vq), gradient_z_quadrature_values, view(gradients_z, :, e))
    # integrate over the element
    for i in eachindex(u_quadrature_values)
      gradients_i = SVector(gradient_x_quadrature_values[i],
                            gradient_y_quadrature_values[i],
                            gradient_z_quadrature_values[i])
      integral += mesh.md.wJq[i, e] * func(u_quadrature_values[i], gradients_i, equations_parabolic)
    end
  end
  return integral
end
@DanielDoehring DanielDoehring added bug Something isn't working enhancement New feature or request labels Nov 5, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working enhancement New feature or request
Projects
None yet
Development

No branches or pull requests

1 participant