From 1b9b0ebce2b23baf82a2bfe40e3f9b554a5d2371 Mon Sep 17 00:00:00 2001 From: cibinjoseph Date: Fri, 20 Oct 2023 14:58:29 -0600 Subject: [PATCH 1/3] Add sharpTE provision --- src/FLOWPanel_postprocess.jl | 120 ++++++++++++++++++++++++++++++++--- 1 file changed, 110 insertions(+), 10 deletions(-) diff --git a/src/FLOWPanel_postprocess.jl b/src/FLOWPanel_postprocess.jl index 1295976..6be3e0b 100644 --- a/src/FLOWPanel_postprocess.jl +++ b/src/FLOWPanel_postprocess.jl @@ -619,7 +619,8 @@ end function calcfield_Ugradmu!(out::AbstractMatrix, body::AbstractBody, areas::AbstractVector, normals::AbstractMatrix, controlpoints::AbstractMatrix; - fieldname="Ugradmu", addfield=true, Gammai=1) + fieldname="Ugradmu", addfield=true, Gammai=1, + sharpTE=false) # Error cases @assert size(out, 1)==3 && size(out, 2)==body.ncells ""* @@ -680,32 +681,131 @@ function calcfield_Ugradmu!(out::AbstractMatrix, body::AbstractBody, A[3, 3] = t3[2]-t0[2] b .= nodal_data[vtx] + res = A\b dx = res[2] dy = res[3] # Transform slopes back to global coordinate system - out[:, i] = @. -0.5 * (dx*e1 + dy*e2) + out[:, i] = @. (dx*e1 + dy*e2) + end + + # If it's a sharp TE, do not use contribution from the other side of the mesh + # while converting from cell to nodal data + if sharpTE && body.grid.orggrid.loop_dim == 1 + if body.grid.dimsplit == 1 + # Compute TE node indices + lin_node = LinearIndices(body.grid._ndivsnodes) + TE_idx = lin_node[1, :, 1] + + # Compute trailing cell indices that share vertices with TE nodes + lin_cell = LinearIndices(body.grid._ndivscells[1:2]) + cells_U = vec(lin_cell[end-1:end, :]) + cells_L = vec(lin_cell[1:2, :]) + + nodal_data_U, nodal_data_L = gt.get_nodal_data_TEcells(body.grid, Gammas, + TE_idx, + cells_U, cells_L; + areas=areas) + + # Overwrite TE node values for cells on upper side + nodal_data[TE_idx] .= nodal_data_U + for i in cells_U + # Convert cell vertices to a local x,y coordinate frame + vtx = gt.get_cell(body.grid, i) + gt.project_3d_2d!(t2, t3, e1, e2, + nodes[:, vtx[1]], + nodes[:, vtx[2]], + nodes[:, vtx[3]]) + + # The (x, y) coordinate of t1 is always at origin + t0 = @. (t1 + t2 + t3) / 3.0 + + # Find slope of 'plane' created by field values at vertices + # in the local coordinate system + A[1, 1] = 1.0 + A[1, 2] = t1[1]-t0[1] + A[1, 3] = t1[2]-t0[2] + + A[2, 1] = 1.0 + A[2, 2] = t2[1]-t0[1] + A[2, 3] = t2[2]-t0[2] + + A[3, 1] = 1.0 + A[3, 2] = t3[1]-t0[1] + A[3, 3] = t3[2]-t0[2] + + b .= nodal_data[vtx] + + res = A\b + dx = res[2] + dy = res[3] + + # Transform slopes back to global coordinate system + out[:, i] = @. (dx*e1 + dy*e2) + end + + # Overwrite TE node values for cells on lower side + nodal_data[TE_idx] .= nodal_data_L + for i in cells_L + # Convert cell vertices to a local x,y coordinate frame + vtx = gt.get_cell(body.grid, i) + gt.project_3d_2d!(t2, t3, e1, e2, + nodes[:, vtx[1]], + nodes[:, vtx[2]], + nodes[:, vtx[3]]) + + # The (x, y) coordinate of t1 is always at origin + t0 = @. (t1 + t2 + t3) / 3.0 + + # Find slope of 'plane' created by field values at vertices + # in the local coordinate system + A[1, 1] = 1.0 + A[1, 2] = t1[1]-t0[1] + A[1, 3] = t1[2]-t0[2] + + A[2, 1] = 1.0 + A[2, 2] = t2[1]-t0[1] + A[2, 3] = t2[2]-t0[2] + + A[3, 1] = 1.0 + A[3, 2] = t3[1]-t0[1] + A[3, 3] = t3[2]-t0[2] + + b .= nodal_data[vtx] + + res = A\b + dx = res[2] + dy = res[3] + + # Transform slopes back to global coordinate system + out[:, i] = @. (dx*e1 + dy*e2) + end + end end + # Gamma / 2 + @. out *= -0.5 + # Save field in body if addfield + add_field(body, "gamma_node", "scalar", nodal_data, "node") add_field(body, fieldname, "vector", eachcol(out), "cell") end end function calcfield_Ugradmu!(out::AbstractMatrix, mbody::MultiBody, - areas::AbstractVector, - normals::AbstractMatrix, - controlpoints::AbstractMatrix, args...; - fieldname="Ugradmu", addfield=true, - optargs... - ) + areas::AbstractVector, + normals::AbstractMatrix, + controlpoints::AbstractMatrix, args...; + fieldname="Ugradmu", addfield=true, + optargs... + ) # Error cases @assert size(out, 1)==3 && size(out, 2)==mbody.ncells ""* - "Invalid `out` matrix."* - " Expected size $((3, mbody.ncells)); got $(size(out))." + "Invalid `out` matrix."* + " Expected size $((3, mbody.ncells)); got $(size(out))." @assert length(areas)==mbody.ncells ""* "Invalid `areas` vector."* " Expected length $(mbody.ncells); got $(length(areas))." From 11aeaac70dad0642e5fa6635f10b895e76824b86 Mon Sep 17 00:00:00 2001 From: cibinjoseph Date: Fri, 20 Oct 2023 16:09:04 -0600 Subject: [PATCH 2/3] Add get_tri_gradient!() in GT and use that instead of redundant code --- src/FLOWPanel_postprocess.jl | 72 ++++++------------------------------ 1 file changed, 11 insertions(+), 61 deletions(-) diff --git a/src/FLOWPanel_postprocess.jl b/src/FLOWPanel_postprocess.jl index 6be3e0b..a0b2ec9 100644 --- a/src/FLOWPanel_postprocess.jl +++ b/src/FLOWPanel_postprocess.jl @@ -647,12 +647,13 @@ function calcfield_Ugradmu!(out::AbstractMatrix, body::AbstractBody, # Pre-allocate required arrays A = Array{Float64}(undef, 3, 3) - b = zeros(3) + t0 = zeros(2) t1 = zeros(2) t2 = zeros(2) t3 = zeros(2) e1 = zeros(3) e2 = zeros(3) + grad = zeros(3) # Compute cell-based gradient for each cell for i = 1:prod(body.grid._ndivscells[1:2]) @@ -666,28 +667,11 @@ function calcfield_Ugradmu!(out::AbstractMatrix, body::AbstractBody, # The (x, y) coordinate of t1 is always at origin t0 = @. (t1 + t2 + t3) / 3.0 - # Find slope of 'plane' created by field values at vertices - # in the local coordinate system - A[1, 1] = 1.0 - A[1, 2] = t1[1]-t0[1] - A[1, 3] = t1[2]-t0[2] - - A[2, 1] = 1.0 - A[2, 2] = t2[1]-t0[1] - A[2, 3] = t2[2]-t0[2] - - A[3, 1] = 1.0 - A[3, 2] = t3[1]-t0[1] - A[3, 3] = t3[2]-t0[2] - - b .= nodal_data[vtx] - - res = A\b - dx = res[2] - dy = res[3] + # Get gradient of plane formed by three scalar values at vertices + gt.get_tri_gradient!(grad, t1, t2, t3, t0, e1, e2, A, nodal_data[vtx]) # Transform slopes back to global coordinate system - out[:, i] = @. (dx*e1 + dy*e2) + out[:, i] = @. (grad[2]*e1 + grad[3]*e2) end # If it's a sharp TE, do not use contribution from the other side of the mesh @@ -721,28 +705,11 @@ function calcfield_Ugradmu!(out::AbstractMatrix, body::AbstractBody, # The (x, y) coordinate of t1 is always at origin t0 = @. (t1 + t2 + t3) / 3.0 - # Find slope of 'plane' created by field values at vertices - # in the local coordinate system - A[1, 1] = 1.0 - A[1, 2] = t1[1]-t0[1] - A[1, 3] = t1[2]-t0[2] - - A[2, 1] = 1.0 - A[2, 2] = t2[1]-t0[1] - A[2, 3] = t2[2]-t0[2] - - A[3, 1] = 1.0 - A[3, 2] = t3[1]-t0[1] - A[3, 3] = t3[2]-t0[2] - - b .= nodal_data[vtx] - - res = A\b - dx = res[2] - dy = res[3] + # Get gradient of plane formed by three scalar values at vertices + gt.get_tri_gradient!(grad, t1, t2, t3, t0, e1, e2, A, nodal_data[vtx]) # Transform slopes back to global coordinate system - out[:, i] = @. (dx*e1 + dy*e2) + out[:, i] = @. (grad[2]*e1 + grad[3]*e2) end # Overwrite TE node values for cells on lower side @@ -758,28 +725,11 @@ function calcfield_Ugradmu!(out::AbstractMatrix, body::AbstractBody, # The (x, y) coordinate of t1 is always at origin t0 = @. (t1 + t2 + t3) / 3.0 - # Find slope of 'plane' created by field values at vertices - # in the local coordinate system - A[1, 1] = 1.0 - A[1, 2] = t1[1]-t0[1] - A[1, 3] = t1[2]-t0[2] - - A[2, 1] = 1.0 - A[2, 2] = t2[1]-t0[1] - A[2, 3] = t2[2]-t0[2] - - A[3, 1] = 1.0 - A[3, 2] = t3[1]-t0[1] - A[3, 3] = t3[2]-t0[2] - - b .= nodal_data[vtx] - - res = A\b - dx = res[2] - dy = res[3] + # Get gradient of plane formed by three scalar values at vertices + gt.get_tri_gradient!(grad, t1, t2, t3, t0, e1, e2, A, nodal_data[vtx]) # Transform slopes back to global coordinate system - out[:, i] = @. (dx*e1 + dy*e2) + out[:, i] = @. (grad[2]*e1 + grad[3]*e2) end end end From e19d9b35c3c5a7faea03ee28411c4e391ca9688a Mon Sep 17 00:00:00 2001 From: cibinjoseph Date: Fri, 20 Oct 2023 16:11:13 -0600 Subject: [PATCH 3/3] Use sharpTE option when computing Ugradmu to avoid peak at TE --- examples/sweptwing.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/sweptwing.jl b/examples/sweptwing.jl index 1f3cf75..b18832a 100644 --- a/examples/sweptwing.jl +++ b/examples/sweptwing.jl @@ -154,7 +154,7 @@ println("Post processing...") # doublet strength to get an accurate surface velocity # Calculate surface velocity U_∇μ due to the gradient of the doublet strength -UDeltaGamma = pnl.calcfield_Ugradmu(body) +UDeltaGamma = pnl.calcfield_Ugradmu(body; sharpTE=true) # Add both velocities together pnl.addfields(body, "Ugradmu", "U")