Skip to content

Commit

Permalink
Reduction for linf_divb analysis (#2206)
Browse files Browse the repository at this point in the history
  • Loading branch information
@DanielDoehring
DanielDoehring authored Dec 12, 2024
1 parent ab8ecad commit 07cbc30
Show file tree
Hide file tree
Showing 4 changed files with 26 additions and 26 deletions.
42 changes: 21 additions & 21 deletions src/callbacks_step/analysis_dg1d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,81 +40,81 @@ function create_cache_analysis(analyzer, mesh::StructuredMesh{1},
end

function calc_error_norms(func, u, t, analyzer,
mesh::StructuredMesh{1}, equations, initial_condition,
mesh::TreeMesh{1}, equations, initial_condition,
dg::DGSEM, cache, cache_analysis)
@unpack vandermonde, weights = analyzer
@unpack node_coordinates, inverse_jacobian = cache.elements
@unpack u_local, x_local, jacobian_local = cache_analysis
@unpack node_coordinates = cache.elements
@unpack u_local, x_local = cache_analysis

# Set up data structures
l2_error = zero(func(get_node_vars(u, equations, dg, 1, 1), equations))
linf_error = copy(l2_error)
total_volume = zero(real(mesh))

# Iterate over all elements for error calculations
for element in eachelement(dg, cache)
# Interpolate solution and node locations to analysis nodes
multiply_dimensionwise!(u_local, vandermonde, view(u, :, :, element))
multiply_dimensionwise!(x_local, vandermonde,
view(node_coordinates, :, :, element))
multiply_scalar_dimensionwise!(jacobian_local, vandermonde,
inv.(view(inverse_jacobian, :, element)))

# Calculate errors at each analysis node
volume_jacobian_ = volume_jacobian(element, mesh, cache)

for i in eachnode(analyzer)
u_exact = initial_condition(get_node_coords(x_local, equations, dg, i), t,
equations)
diff = func(u_exact, equations) -
func(get_node_vars(u_local, equations, dg, i), equations)
# We take absolute value as we need the Jacobian here for the volume calculation
abs_jacobian_local_i = abs(jacobian_local[i])

l2_error += diff .^ 2 * (weights[i] * abs_jacobian_local_i)
l2_error += diff .^ 2 * (weights[i] * volume_jacobian_)
linf_error = @. max(linf_error, abs(diff))
total_volume += weights[i] * abs_jacobian_local_i
end
end

# For L2 error, divide by total volume
l2_error = @. sqrt(l2_error / total_volume)
total_volume_ = total_volume(mesh)
l2_error = @. sqrt(l2_error / total_volume_)

return l2_error, linf_error
end

function calc_error_norms(func, u, t, analyzer,
mesh::TreeMesh{1}, equations, initial_condition,
mesh::StructuredMesh{1}, equations, initial_condition,
dg::DGSEM, cache, cache_analysis)
@unpack vandermonde, weights = analyzer
@unpack node_coordinates = cache.elements
@unpack u_local, x_local = cache_analysis
@unpack node_coordinates, inverse_jacobian = cache.elements
@unpack u_local, x_local, jacobian_local = cache_analysis

# Set up data structures
l2_error = zero(func(get_node_vars(u, equations, dg, 1, 1), equations))
linf_error = copy(l2_error)
total_volume = zero(real(mesh))

# Iterate over all elements for error calculations
for element in eachelement(dg, cache)
# Interpolate solution and node locations to analysis nodes
multiply_dimensionwise!(u_local, vandermonde, view(u, :, :, element))
multiply_dimensionwise!(x_local, vandermonde,
view(node_coordinates, :, :, element))
multiply_scalar_dimensionwise!(jacobian_local, vandermonde,
inv.(view(inverse_jacobian, :, element)))

# Calculate errors at each analysis node
volume_jacobian_ = volume_jacobian(element, mesh, cache)

for i in eachnode(analyzer)
u_exact = initial_condition(get_node_coords(x_local, equations, dg, i), t,
equations)
diff = func(u_exact, equations) -
func(get_node_vars(u_local, equations, dg, i), equations)
l2_error += diff .^ 2 * (weights[i] * volume_jacobian_)
# We take absolute value as we need the Jacobian here for the volume calculation
abs_jacobian_local_i = abs(jacobian_local[i])

l2_error += diff .^ 2 * (weights[i] * abs_jacobian_local_i)
linf_error = @. max(linf_error, abs(diff))
total_volume += weights[i] * abs_jacobian_local_i
end
end

# For L2 error, divide by total volume
total_volume_ = total_volume(mesh)
l2_error = @. sqrt(l2_error / total_volume_)
l2_error = @. sqrt(l2_error / total_volume)

return l2_error, linf_error
end
Expand Down Expand Up @@ -214,7 +214,7 @@ function analyze(::Val{:linf_divb}, du, u, t,

# integrate over all elements to get the divergence-free condition errors
linf_divb = zero(eltype(u))
for element in eachelement(dg, cache)
@batch reduction=(max, linf_divb) for element in eachelement(dg, cache)
for i in eachnode(dg)
divb = zero(eltype(u))
for k in eachnode(dg)
Expand Down
4 changes: 2 additions & 2 deletions src/callbacks_step/analysis_dg2d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -338,7 +338,7 @@ function analyze(::Val{:linf_divb}, du, u, t,

# integrate over all elements to get the divergence-free condition errors
linf_divb = zero(eltype(u))
for element in eachelement(dg, cache)
@batch reduction=(max, linf_divb) for element in eachelement(dg, cache)
for j in eachnode(dg), i in eachnode(dg)
divb = zero(eltype(u))
for k in eachnode(dg)
Expand Down Expand Up @@ -368,7 +368,7 @@ function analyze(::Val{:linf_divb}, du, u, t,

# integrate over all elements to get the divergence-free condition errors
linf_divb = zero(eltype(u))
for element in eachelement(dg, cache)
@batch reduction=(max, linf_divb) for element in eachelement(dg, cache)
for j in eachnode(dg), i in eachnode(dg)
divb = zero(eltype(u))
# Get the contravariant vectors Ja^1 and Ja^2
Expand Down
4 changes: 2 additions & 2 deletions src/callbacks_step/analysis_dg3d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -378,7 +378,7 @@ function analyze(::Val{:linf_divb}, du, u, t,

# integrate over all elements to get the divergence-free condition errors
linf_divb = zero(eltype(u))
for element in eachelement(dg, cache)
@batch reduction=(max, linf_divb) for element in eachelement(dg, cache)
for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
divb = zero(eltype(u))
for l in eachnode(dg)
Expand Down Expand Up @@ -416,7 +416,7 @@ function analyze(::Val{:linf_divb}, du, u, t,

# integrate over all elements to get the divergence-free condition errors
linf_divb = zero(eltype(u))
for element in eachelement(dg, cache)
@batch reduction=(max, linf_divb) for element in eachelement(dg, cache)
for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
divb = zero(eltype(u))
# Get the contravariant vectors Ja^1, Ja^2, and Ja^3
Expand Down
2 changes: 1 addition & 1 deletion src/callbacks_step/analysis_dgmulti.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -125,7 +125,7 @@ function analyze(::Val{:linf_divb}, du, u, t,
uEltype = eltype(B1)
linf_divB = zero(uEltype)
local_divB = zeros(uEltype, size(B1, 1))
for e in eachelement(mesh, dg, cache)
@batch reduction=(max, linf_divb) for e in eachelement(mesh, dg, cache)
compute_local_divergence!(local_divB, e, view.(B, :, e), mesh, dg, cache)

# compute maximum norm
Expand Down

0 comments on commit 07cbc30

Please sign in to comment.