diff --git a/examples/Example240_FiniteElementVelocities.jl b/examples/Example240_FiniteElementVelocities.jl index 871bc7359..ecf484263 100644 --- a/examples/Example240_FiniteElementVelocities.jl +++ b/examples/Example240_FiniteElementVelocities.jl @@ -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$. =# @@ -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) @@ -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 @@ -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) @@ -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 @@ -214,55 +187,72 @@ 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] @@ -270,8 +260,8 @@ 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] diff --git a/ext/VoronoiFVMExtendableFEMBaseExt.jl b/ext/VoronoiFVMExtendableFEMBaseExt.jl index 67b51d3c8..ba0e2c282 100644 --- a/ext/VoronoiFVMExtendableFEMBaseExt.jl +++ b/ext/VoronoiFVMExtendableFEMBaseExt.jl @@ -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)