Skip to content

Commit

Permalink
corrected edgevelocities computation; restructured and cleaned up Exa…
Browse files Browse the repository at this point in the history 
…mple script
  • Loading branch information
@Da-Be-Ru
Da-Be-Ru committed Oct 21, 2024
1 parent 4c2922d commit be921f5
Show file tree
Hide file tree
Showing 2 changed files with 97 additions and 107 deletions.
202 changes: 96 additions & 106 deletions examples/Example240_FiniteElementVelocities.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,19 +6,20 @@
Solve the equation
```math
-\nabla ( D \nabla u - v u) = 0
-\nabla \cdot ( D \nabla u - \mathbf{v} u) = 0
```
in $\Omega=(0,L)\times (0,H)$ with a homogeneous Neumann boundary condition
at $x=0$, an outflow boundary condition at $x=L$, a Dirichlet inflow
condition at $y=H$, and a homogeneous Dirichlet boundary condition
on part of $y=0$.
The analytical expression for the velocity field is $v(x,y)=(x^2,-2xy)$ in
cartesian coordinates and $v(r,z)=(r^2,-3rz)$ in cylindrical coordinates, i.e.
The analytical expression for the (quadratic variant of the) velocity field
is $\mathbf{v}(x,y)=(x^2,-2xy)$ in cartesian coordinates and (for the linear variant)
$\mathbf{v}(r,z)=(r,-2z)$ in cylindrical coordinates, i.e.
where the system is solved on $\Omega$ to represent a solution on the solid
of revolution arising from rotating $\Omega$ around $x=0$.
We compute the solution $u$ in both coordinate systems where $v$ is given
We compute the solution $u$ in both coordinate systems where $\mathbf{v}$ is given
as an analytical expression and as a finite element interpolation onto
the grid of $\Omega$.
=#
Expand All @@ -36,11 +37,12 @@ function stagnation_flow_cartesian(x, y)
return (x^2, -2x * y)
end

# in cylindrical case: since the reconstruction space HDIVBDM2
# is only quadratic, but we have to reconstruct r*v for a
# div-free solution, so we can only resolve at most the linear case exactly
# In the cylindrical case: since the reconstruction space $\mathtt{HDIVBDM2}$
# is only quadratic, but we have to reconstruct $r \, \mathbf{v}(r,z)$ for a
# (reconstructed) divergence-free solution,
# we can only resolve at most the linear case exactly.
function stagnation_flow_cylindrical(r, z)
return (r, -2*z)
return (r, -2 * z)
end

function inflow_cylindrical(u, qpinfo)
Expand All @@ -61,17 +63,17 @@ function flux!(f, u, edge, data)
end

function bconditions!(f, u, node, data)
# catalytic Dirichlet condition at y=0
## catalytic Dirichlet condition at y=0
if node.region == 5
boundary_dirichlet!(f, u, node, 1, node.region, 0.0)
end

# outflow condition at x = L
## outflow condition at x = L
if node.region == 2
f[1] = data.bfvelo[node.ibnode, node.ibface] * u[1]
end

# inflow condition at y = H
## inflow condition at y = H
if node.region == 3
boundary_dirichlet!(f, u, node, 1, node.region, data.cin)
end
Expand All @@ -86,81 +88,43 @@ mutable struct Data
Data() = new()
end

function main(; coord_system = Cartesian2D, usefem = false, nref = 0, Plotter = nothing,
D = 0.01, cin = 1.0, assembly = :edgewise)
if coord_system == Cylindrical2D
analytical_velocity = stagnation_flow_cylindrical
else
analytical_velocity = stagnation_flow_cartesian
end

data = Data()
data.D = D
data.cin = cin

#=
Calculate the analytical or FEM solution to the stagnation point flow field $\mathbf{v}$ and
use this as input to solve for the species concentration $u$ of the corresponding
convection-diffusion system.
The passed flags regulate the following behavior:
- `cylindrical_grid`: if true, calculates both the velocity field $\mathbf{v}(r,z)$
and the species concentration $u(r,z)$ in cylindrical coordinates, assuming
rotationally symmetric solutions for both.
- `usefem`: if true, calculates the velocity field $\mathbf{v}$ using the
finite element method provided by [`ExtendableFEM`](https://github.com/chmerdon/ExtendableFEM.jl).
- `reconst`: if true, interpolates the FEM-calculated
velocity field onto a "reconstruction" finite element space that
provides an exactly divergence-free solution. In a cylindrical grid,
this returns not $\mathbf{v}(r,z)$, but $r \, \mathbf{v}(r,z)$ as the velocity.
- `use_different_grids`: if true, calculates the FEM solution of the
velocity field on a uniform `flowgrid` and the species concentration
on an adaptively refined `chemgrid` while still interpolating the
calculated velocity correctly onto the `chemgrid`.
=#
function main(;
cylindrical_grid = false, usefem = true, reconst = cylindrical_grid, use_different_grids = false, nref = 1,
Plotter = nothing, μ = 1.0e-02, D = 0.01, cin = 1.0, assembly = :edgewise,
interpolation_eps = 1.0e-09)
H = 1.0
L = 5.0
grid = simplexgrid(range(0, L; length = 20 * 2^nref),
range(0, H; length = 5 * 2^nref))
bfacemask!(grid, [L / 4, 0.0], [3 * L / 4, 0.0], 5)

if usefem
function interpolate_vel!(result, qpinfo)
x = qpinfo.x
result .= analytical_velocity(x[1], x[2])
end
FES = FESpace{H1P2{2, 2}}(grid)
fem_velocity = FEVector(FES)[1]
interpolate!(fem_velocity, interpolate_vel!)

evelo = edgevelocities(grid, fem_velocity)
bfvelo = bfacevelocities(grid, fem_velocity)
if cylindrical_grid
coord_system = Cylindrical2D
else
evelo = edgevelocities(grid, analytical_velocity)
bfvelo = bfacevelocities(grid, analytical_velocity)
coord_system = Cartesian2D
end

data.evelo = evelo
data.bfvelo = bfvelo

physics = VoronoiFVM.Physics(; flux = flux!, breaction = bconditions!, data)
sys = VoronoiFVM.System(grid, physics; assembly = assembly)
enable_species!(sys, 1, [1])

sol = solve(sys; inival = 0.0)

vis = GridVisualizer(; Plotter = Plotter)

scalarplot!(vis[1, 1], grid, sol[1, :]; flimits = (0, cin + 1.0e-5),
show = true)
sol, evelo, bfvelo
end

using Test
function runtests()
for coord_system in [Cartesian2D, Cylindrical2D]
sol_analytical, evelo_analytical, bfvelo_analytical = main(;
coord_system, usefem = false)
sol_fem, evelo_fem, bfvelo_fem = main(; coord_system, usefem = true)
@test norm(evelo_analytical .- evelo_fem, Inf) ≤ 1.0e-11
@test norm(bfvelo_analytical .- bfvelo_fem, Inf) ≤ 1.0e-09
@test norm(sol_analytical .- sol_fem, Inf) ≤ 1.0e-10
end
end

# test analytical distance and calc_divergences
# TODO: make this properly work for cylindrical case (where this is still inexplicably broken)
function full_fem_demo(;
coord_system = Cartesian2D, nref = 1, usefem = true, usedifferentgrids = false,
Plotter = nothing, μ = 1.0e-02, D = 0.01, cin = 1.0, assembly = :edgewise,
interpolation_eps = 1.0e-09)
H = 1.0
L = 5.0

flowgrid = simplexgrid(range(0, L; length = 20 * 2^nref),
range(0, H; length = 5 * 2^nref))

if usedifferentgrids
if use_different_grids
h_fine = 1.0e-01
X_bottom = geomspace(0.0, L / 2, 5.0e-01, h_fine)
X_cat = range(L / 2, L; step = h_fine)
Expand All @@ -173,26 +137,35 @@ function full_fem_demo(;
end

if usefem
velocity = compute_velocity(flowgrid, coord_system, μ; interpolation_eps)
velocity = compute_velocity(
flowgrid, cylindrical_grid, reconst, μ; interpolation_eps)
DivIntegrator = L2NormIntegrator([div(1)]; quadorder = 2 * 2, resultdim = 1)
div_v = sqrt(sum(evaluate(DivIntegrator, [velocity])))
@info "||div(R(v))||_2 = $(div_v)"
else
if coord_system == Cartesian2D
velocity = stagnation_flow_cartesian
elseif coord_system == Cylindrical2D
if cylindrical_grid
velocity = stagnation_flow_cylindrical
else
velocity = stagnation_flow_cartesian
end
end

chemgrid[CoordinateSystem] = coord_system
if coord_system == Cylindrical2D
analytical_velocity = stagnation_flow_cylindrical
else
analytical_velocity = stagnation_flow_cartesian
end

if cylindrical_grid
chemgrid[CoordinateSystem] = Cylindrical2D
end

data = Data()
data.D = D
data.cin = cin

evelo = edgevelocities(chemgrid, velocity)
bfvelo = bfacevelocities(chemgrid, velocity)
evelo = edgevelocities(chemgrid, velocity; reconst)
bfvelo = bfacevelocities(chemgrid, velocity; reconst)

data.evelo = evelo
data.bfvelo = bfvelo
Expand All @@ -214,64 +187,81 @@ function full_fem_demo(;
minmax = extrema(sol)
@info "Minimal/maximal values of concentration: $(minmax)"

return sys, velocity, sol, evelo, bfvelo, fvm_divs
#return velocity
return sol, evelo, bfvelo, minmax
end

function compute_velocity(flowgrid, coord_system, μ = 1.0e-02; interpolation_eps = 1.0e-10)
axisymmetric = coord_system == Cylindrical2D ? true : false
using Test
function runtests()
cin = 1.0
for cylindrical_grid in [false, true]
sol_analytical, evelo_analytical, bfvelo_analytical, minmax_analytical = main(;
cylindrical_grid, cin, usefem = false)
sol_fem, evelo_fem, bfvelo_fem, minmax_fem = main(;
cylindrical_grid, cin, usefem = true)
@test norm(evelo_analytical .- evelo_fem, Inf) ≤ 1.0e-11
@test norm(bfvelo_analytical .- bfvelo_fem, Inf) ≤ 1.0e-09
@test norm(sol_analytical .- sol_fem, Inf) ≤ 1.0e-10
@test norm(minmax_analytical .- [0.,cin], Inf) ≤ 1.0e-15
@test norm(minmax_fem .- [0.,cin], Inf) ≤ 1.0e-11
end
end

# define finite element spaces
function compute_velocity(
flowgrid, cylindrical_grid, reconst, μ = 1.0e-02; interpolation_eps = 1.0e-10)
## define finite element spaces
FE_v, FE_p = H1P2B{2, 2}, L2P1{1}
reconst_FEType = HDIVBDM2{2}
FES = [FESpace{FE_v}(flowgrid), FESpace{FE_p}(flowgrid; broken = true)]

# describe problem
## describe problem
Problem = ProblemDescription("incompressible Stokes problem")
v = Unknown("v"; name = "velocity")
p = Unknown("p"; name = "pressure")
assign_unknown!(Problem, v)
assign_unknown!(Problem, p)

# assign stokes operator
## assign stokes operator
assign_operator!(Problem,
BilinearOperator(
kernel_stokes!, axisymmetric ? [id(v), grad(v), id(p)] : [grad(v), id(p)];
kernel_stokes!, cylindrical_grid ? [id(v), grad(v), id(p)] : [grad(v), id(p)];
bonus_quadorder = 2, store = false,
params = [μ, axisymmetric], verbosity = 2))
params = [μ, cylindrical_grid]))

# assign boundary condition operators
if axisymmetric
# inflow
## assign Dirichlet boundary conditions on all boundary regions to
## enforce match with analytical solution
if cylindrical_grid
assign_operator!(
Problem, InterpolateBoundaryData(v, inflow_cylindrical; regions = [1, 2, 3, 4]))
else
# inflow and outflow
assign_operator!(
Problem, InterpolateBoundaryData(v, inflow_cartesian; regions = [1, 2, 3, 4]))
end

velocity_solution = solve(Problem, FES)

# ensure divergence free solution by projecting onto reconstruction spaces
## ensure divergence free solution by projecting onto reconstruction spaces
FES_reconst = FESpace{reconst_FEType}(flowgrid)
R = FEVector(FES_reconst)
if axisymmetric
lazy_interpolate!(R[1], velocity_solution, [id(v)]; postprocess = multiply_r,
bonus_quadorder = 2, eps = interpolation_eps)
if reconst
if cylindrical_grid
lazy_interpolate!(R[1], velocity_solution, [id(v)]; postprocess = multiply_r,
bonus_quadorder = 2, eps = interpolation_eps)
else
lazy_interpolate!(
R[1], velocity_solution, [id(v)];
bonus_quadorder = 2, eps = interpolation_eps)
end
else
lazy_interpolate!(
R[1], velocity_solution, [id(v)];
bonus_quadorder = 2, eps = interpolation_eps)
return velocity_solution[1]
end

return R[1]
end

function kernel_stokes!(result, u_ops, qpinfo)
μ = qpinfo.params[1]
axisymmetric = qpinfo.params[2]
if axisymmetric > 0
cylindrical_grid = qpinfo.params[2]
if cylindrical_grid > 0
r = qpinfo.x[1]
u, ∇u, p = view(u_ops, 1:2), view(u_ops, 3:6), view(u_ops, 7)
result[1] = μ / r * u[1] - p[1]
Expand Down
2 changes: 1 addition & 1 deletion ext/VoronoiFVMExtendableFEMBaseExt.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -125,7 +125,7 @@ function _integrate_along_segments(p1, p2, hnormal, aug_vec_block::AugmentedFEVe

edge_length = norm(p1 - p2, 2)
avg_r = (p1[1] + p2[1]) / 2
if avg_r < eps()
if axisymmetric && avg_r < eps()
return 0
end
bp1 = zeros(3)
Expand Down

0 comments on commit be921f5

Please sign in to comment.