|
| 1 | +#= |
| 2 | +
|
| 3 | +# 264 : Stokes+Darcy |
| 4 | +([source code](@__SOURCE_URL__)) |
| 5 | +
|
| 6 | +This example solves the coupled Stokes-Darcy problem in a domain |
| 7 | +``\Omega := \Omega_\text{FF} \cup \Omega_\text{PM}`` that has a free flow region |
| 8 | +``\Omega_\text{FF}`` and a porous media regions ``\Omega_\text{PM}``. |
| 9 | +In the free flow region a Stokes problem is solved that seeks a velocity |
| 10 | +``\mathbf{u}_\text{FF}`` and a pressure ``\mathbf{p}_\text{FF}`` |
| 11 | +such that |
| 12 | +```math |
| 13 | +\begin{aligned} |
| 14 | +- 2\mu \mathrm{div}(\epsilon(\mathbf{u}_\text{FF}) - p_\text{FF}I) & = \mathbf{f}_\text{FF}\\ |
| 15 | +\mathrm{div}(\mathbf{u}_\text{FF}) & = 0. |
| 16 | +\end{aligned} |
| 17 | +``` |
| 18 | +In the porous media region the Darcy problem is solved that seeks a velocity ``\mathbf{u}_\text{PM}`` |
| 19 | +and a pressure ``\mathbf{p}_\text{PM}`` such that |
| 20 | +```math |
| 21 | +\begin{aligned} |
| 22 | +\mathbf{u}_\text{PM} + k \nabla p_\text{PM} & = 0\\ |
| 23 | +\mathrm{div}(\mathbf{u}_\text{PM}) & = f_\text{PM}. |
| 24 | +\end{aligned} |
| 25 | +``` |
| 26 | +On the interface ``\Gamma := \partial \Omega_\text{FF} \cap \partial \Omega_\text{PM}`` |
| 27 | +the two velocities are coupled via several conditions, i.e., the conservation of mass |
| 28 | +```math |
| 29 | +\mathbf{u}_\text{FF} \cdot \mathbf{n} = \mathbf{u}_\text{PM} \cdot \mathbf{n} \quad \text{on } \Gamma, |
| 30 | +``` |
| 31 | +the balance of normal forces |
| 32 | +```math |
| 33 | +p_\text{FF} - 2\mu \epsilon(\mathbf{u}_\text{FF}) \mathbf{n}_\text{FF} \cdot \mathbf{n}_\text{FF} = p_\text{PM}, |
| 34 | +``` |
| 35 | +and the Beavers-Joseph-Saffman condition |
| 36 | +```math |
| 37 | +\mathbf{u}_\text{FF} \cdot \mathbf{\tau} |
| 38 | += -\frac{\sqrt{\mu k}}{\mu \alpha} 2 \epsilon(\mathbf{u}_\text{FF}) \mathbf{n}_\text{FF} \cdot \mathbf{\tau}. |
| 39 | +``` |
| 40 | +The interface condition for the normal fluxes is realized weakly via a Lagrange multiplier ``\lambda`` |
| 41 | +that only lives on the interface. |
| 42 | +
|
| 43 | +The weak formulation leads to the problem: |
| 44 | +seek ``(\mathbf{u}_\text{FF}, \mathbf{u}_\text{PM}, p_\text{FF}, p_\text{PM}, \lambda)`` |
| 45 | +such that, for all ``(\mathbf{v}_\text{FF}, \mathbf{v}_\text{PM}, q_\text{FF}, q_\text{PM}, \chi)``, |
| 46 | +```math |
| 47 | +\begin{aligned} |
| 48 | +a_1(\mathbf{u}_\text{FF},\mathbf{v}_\text{FF}) + a_2(\mathbf{u}_\text{PM}, \mathbf{v}_\text{PM}) |
| 49 | ++ b_1(p_\text{FF},\mathbf{v}_\text{FF}) + b_2(p_\text{PM},\mathbf{v}_\text{PM}) |
| 50 | ++ b_{\Gamma}(\mathbf{v}_\text{FF} - \mathbf{v}_\text{PM}, \lambda) |
| 51 | +& = (\mathbf{f}_\text{FF}, \mathbf{v}_\text{FF})_{L^2(\Omega_\text{FF})}\\ |
| 52 | +b_1(q_\text{FF},\mathbf{u}_\text{FF}) + b_2(q_\text{PM},\mathbf{u}_\text{PM}) & = (f_\text{PM}, q_\text{PM})_{L^2(\Omega_\text{PM})}\\ |
| 53 | +b_{\Gamma}(\mathbf{u}_\text{FF} - \mathbf{u}_\text{PM}, \chi) & = 0. |
| 54 | +\end{aligned} |
| 55 | +``` |
| 56 | +
|
| 57 | +The bilinearforms read |
| 58 | +```math |
| 59 | +\begin{aligned} |
| 60 | +a_1(\mathbf{u}_\text{FF},\mathbf{v}_\text{FF}) & := 2\mu (\epsilon(\mathbf{u}_\text{FF}), \epsilon(\mathbf{v}_\text{FF}))_{L^2(\Omega_\text{FF})} |
| 61 | ++ \frac{αμ}{\sqrt{μk}} (\mathbf{u}_\text{FF} \cdot \mathbf{\tau},\mathbf{v}_\text{FF} \cdot \mathbf{\tau})_{L^2(\Gamma)}\\ |
| 62 | +a_2(\mathbf{u}_\text{PM}, \mathbf{v}_\text{PM}) & := (\mathbf{u}_\text{PM}, \mathbf{v}_\text{PM})_{L^2(\Omega_\text{PM})}\\ |
| 63 | +b_1(q_\text{FF},\mathbf{v}_\text{FF}) & := -(\mathrm{div} \mathbf{v}_\text{FF}, q_\text{FF})_{L^2(\Omega_\text{FF})}\\ |
| 64 | +b_2(q_\text{PM},\mathbf{v}_\text{PM}) & := -(\mathrm{div} \mathbf{v}_\text{PM}, q_\text{PM})_{L^2(\Omega_\text{PM})}\\ |
| 65 | +b_{\Gamma}(\mathbf{v}, \lambda) & := (\mathbf{v} \cdot \mathbf{n}, \lambda)_{L^2(\Gamma)} |
| 66 | +\end{aligned} |
| 67 | +``` |
| 68 | +
|
| 69 | +Details on the model can be found e.g. in the reference below. |
| 70 | +
|
| 71 | +!!! reference |
| 72 | +
|
| 73 | + ''Coupling Fluid Flow with Porous Media Flow'' |
| 74 | + SIAM Journal on Numerical Analysis 2002 40:6, 2195-2218 |
| 75 | + [>Link<](https://doi.org/10.1137/S0036142901392766) |
| 76 | +
|
| 77 | +In this example an analytic benchmark problem is solved with the Taylor--Hood FEM in the free flow |
| 78 | +domain and the Raviart--Thomas FEM in the porous media domain. |
| 79 | +
|
| 80 | +The computed solution for the default parameters looks like this: |
| 81 | +
|
| 82 | + |
| 83 | +
|
| 84 | +=# |
| 85 | + |
| 86 | +module Example264_StokesDarcy |
| 87 | + |
| 88 | +using ExtendableFEM |
| 89 | +using ExtendableFEMBase |
| 90 | +using ExtendableGrids |
| 91 | +using GridVisualize |
| 92 | +using SimplexGridFactory |
| 93 | +using Triangulate |
| 94 | +using Metis |
| 95 | +using Test #hide |
| 96 | + |
| 97 | +## exact solution and data functions |
| 98 | +function u!(result, qpinfo) |
| 99 | + x = qpinfo.x |
| 100 | + result[1] = -cos(pi * x[1]) * sin(pi * x[2]) |
| 101 | + result[2] = sin(pi * x[1]) * cos(pi * x[2]) |
| 102 | + return nothing |
| 103 | +end |
| 104 | +function f_FF!(result, qpinfo) |
| 105 | + x = qpinfo.x |
| 106 | + result[1] = -(4 * pi^2 + 1) * cos(pi * x[1]) * sin(pi * x[2]) |
| 107 | + return result[2] = - sin(pi * x[1]) * cos(pi * x[2]) |
| 108 | +end |
| 109 | +function p!(result, qpinfo) |
| 110 | + x = qpinfo.x |
| 111 | + result[1] = -sin(pi * x[1])sin(pi * x[2]) * (1 / pi + 2 * pi) + 4 / pi |
| 112 | + return nothing |
| 113 | +end |
| 114 | +function q!(result, qpinfo) |
| 115 | + k = qpinfo.params[1] |
| 116 | + x = qpinfo.x |
| 117 | + result[1] = -sin(pi * x[1])sin(pi * x[2]) / pi + 4 / pi |
| 118 | + return nothing |
| 119 | +end |
| 120 | +function v!(result, qpinfo) # = -K∇q |
| 121 | + k = qpinfo.params[1] |
| 122 | + x = qpinfo.x |
| 123 | + result[1] = cos(pi * x[1]) * sin(pi * x[2]) |
| 124 | + result[2] = sin(pi * x[1]) * cos(pi * x[2]) |
| 125 | + return nothing |
| 126 | +end |
| 127 | +function f_PM!(result, qpinfo) # = div(v) |
| 128 | + x = qpinfo.x |
| 129 | + return result[1] = -2 * pi * sin(pi * x[1]) * sin(pi * x[2]) |
| 130 | +end |
| 131 | + |
| 132 | +## kernel functions for operators |
| 133 | +function kernel_stokes_standard!(result, u_ops, qpinfo) |
| 134 | + ∇u, p = view(u_ops, 1:4), view(u_ops, 5) |
| 135 | + μ = 2 * qpinfo.params[1] |
| 136 | + result[1] = μ * ∇u[1] - p[1] |
| 137 | + result[2] = μ * (∇u[2] + ∇u[3]) |
| 138 | + result[3] = μ * (∇u[2] + ∇u[3]) |
| 139 | + result[4] = μ * ∇u[4] - p[1] |
| 140 | + result[5] = -(∇u[1] + ∇u[4]) |
| 141 | + return nothing |
| 142 | +end |
| 143 | + |
| 144 | +function coupling_normal!(result, input, qpinfo) |
| 145 | + return result[1] = input[2] - input[1] |
| 146 | +end |
| 147 | + |
| 148 | +function exact_error!(f!::Function) |
| 149 | + return function closure(result, u, qpinfo) |
| 150 | + f!(result, qpinfo) |
| 151 | + result .-= u |
| 152 | + return result .= result .^ 2 |
| 153 | + end |
| 154 | +end |
| 155 | + |
| 156 | +## everything is wrapped in a main function |
| 157 | +function main(; |
| 158 | + μ = 1, # viscosity |
| 159 | + k = 1, # permeability |
| 160 | + α = 1, # parameter in interface condition |
| 161 | + coupled = true, # solve coupled problem with interface conditions or decoupled problems on subgrids ? |
| 162 | + order_u = 2, # polynomial order for free flow velocity |
| 163 | + order_p = order_u - 1, # polynomial order for free flow pressure |
| 164 | + order_v = order_u - 1, # polynomial order for porous media velocity |
| 165 | + nrefs = 4, # number of mesh refinements |
| 166 | + Plotter = nothing, # backend for Plotting (e.g. GLMakie) |
| 167 | + parallel = false, # do parallel assembly? |
| 168 | + npart = 8, # number of partitions for grid coloring (if parallel = true) |
| 169 | + kwargs... |
| 170 | + ) |
| 171 | + |
| 172 | + ## region numbers |
| 173 | + rFF = 1 # free flow region |
| 174 | + rPM = 2 # porous media region |
| 175 | + |
| 176 | + ## load mesh and refine |
| 177 | + xgrid = |
| 178 | + uniform_refine( |
| 179 | + simplexgrid( |
| 180 | + Triangulate; |
| 181 | + points = [0 0; 0 -1; 1 -1; 1 0; 1 1; 0 1]', |
| 182 | + bfaces = [1 2; 2 3; 3 4; 4 5; 5 6; 6 1; 1 4]', |
| 183 | + bfaceregions = [1, 2, 3, 4, 5, 6, 7], |
| 184 | + regionpoints = [0.5 0.5; 0.5 -0.5]', |
| 185 | + regionnumbers = [1, 2], |
| 186 | + regionvolumes = [1.0, 1.0] |
| 187 | + ), nrefs |
| 188 | + ) |
| 189 | + |
| 190 | + if parallel |
| 191 | + xgrid = partition(xgrid, RecursiveMetisPartitioning(npart = npart)) |
| 192 | + end |
| 193 | + |
| 194 | + ## define unknowns |
| 195 | + u = Unknown("u"; name = "velocity FF", dim = 2) |
| 196 | + p = Unknown("p"; name = "pressure FF", dim = 1) |
| 197 | + v = Unknown("q"; name = "velocity PM", dim = 2) |
| 198 | + q = Unknown("q"; name = "pressure PM", dim = 2) |
| 199 | + λ = Unknown("λ"; name = "interface LM", dim = 1) |
| 200 | + |
| 201 | + ## problem description |
| 202 | + PD = ProblemDescription("Stokes-Darcy problem") |
| 203 | + assign_unknown!(PD, u) |
| 204 | + assign_unknown!(PD, p) |
| 205 | + assign_unknown!(PD, v) |
| 206 | + assign_unknown!(PD, q) |
| 207 | + assign_unknown!(PD, λ) |
| 208 | + if coupled |
| 209 | + ## when coupled no boundary conditions on boundary region 7 |
| 210 | + bnd_stokes = [4, 5, 6] |
| 211 | + bnd_darcy = [1, 2, 3] |
| 212 | + else |
| 213 | + bnd_stokes = [4, 5, 6, 7] |
| 214 | + bnd_darcy = [1, 2, 3, 7] |
| 215 | + end |
| 216 | + |
| 217 | + ## define free flow problem: Stokes equations to solve for velocity u and pressure p |
| 218 | + assign_operator!(PD, BilinearOperator(kernel_stokes_standard!, [grad(u), id(p)]; params = [μ], parallel = parallel, regions = [rFF], kwargs...)) |
| 219 | + assign_operator!(PD, LinearOperator(f_FF!, [id(u)]; bonus_quadorder = 4, parallel = parallel, regions = [rFF], kwargs...)) |
| 220 | + assign_operator!(PD, InterpolateBoundaryData(u, u!; bonus_quadorder = 4, regions = bnd_stokes, kwargs...)) |
| 221 | + |
| 222 | + ## define porous media flow: Darcy problem to solve for pressure q |
| 223 | + assign_operator!(PD, BilinearOperator([id(v)]; parallel = parallel, regions = [rPM], kwargs...)) |
| 224 | + assign_operator!(PD, BilinearOperator([div(v)], [id(q)]; factor = -1, transposed_copy = -1, parallel = parallel, regions = [rPM], kwargs...)) |
| 225 | + assign_operator!(PD, LinearOperator(f_PM!, [id(q)]; bonus_quadorder = 4, parallel = parallel, regions = [rPM], kwargs...)) |
| 226 | + |
| 227 | + ## coupling conditions on interface (bregion 7) |
| 228 | + if coupled |
| 229 | + ## interface condition for tangential part |
| 230 | + assign_operator!(PD, BilinearOperator([tangentialflux(u)]; factor = α * μ / sqrt(μ * k), entities = ON_FACES, regions = [7])) |
| 231 | + ## interface condition for normal part via Lagrange multiplier |
| 232 | + assign_operator!(PD, BilinearOperator(coupling_normal!, [id(λ)], [normalflux(u), normalflux(v)]; transposed_copy = 1, entities = ON_FACES, regions = [7])) |
| 233 | + else |
| 234 | + ## dummy problem for Lagrange multiplier |
| 235 | + assign_operator!(PD, BilinearOperator([id(λ)]; entities = ON_FACES, parallel = parallel, regions = [7])) |
| 236 | + ## Dirichlet condition for Darcy problem with different factor (due to different direction of normal) |
| 237 | + assign_operator!(PD, LinearOperator(q!, [normalflux(v)]; bonus_quadorder = 4, factor = 1, entities = ON_FACES, params = [k], regions = [7], kwargs...)) |
| 238 | + end |
| 239 | + assign_operator!(PD, LinearOperator(q!, [normalflux(v)]; bonus_quadorder = 4, factor = -1, entities = ON_BFACES, params = [k], regions = setdiff(bnd_darcy, 7), kwargs...)) |
| 240 | + |
| 241 | + ## generate FESpaces and a solution vector for all 3 unknowns |
| 242 | + FEType_λ = order_v == 0 ? L2P0{1} : H1Pk{1, 1, min(order_v, order_u)} |
| 243 | + FES = [ |
| 244 | + FESpace{H1Pk{2, 2, order_u}}(xgrid; regions = [rFF]), |
| 245 | + FESpace{H1Pk{1, 2, order_p}}(xgrid; regions = [rFF]), |
| 246 | + FESpace{HDIVRTk{2, order_v}}(xgrid; regions = [rPM]), |
| 247 | + FESpace{H1Pk{1, 2, order_v}}(xgrid; regions = [rPM], broken = true), |
| 248 | + FESpace{FEType_λ, ON_FACES}(xgrid; regions = [7], broken = true), |
| 249 | + ] |
| 250 | + sol = FEVector(FES; tags = PD.unknowns) |
| 251 | + |
| 252 | + ## solve the coupled problem |
| 253 | + sol = solve(PD, [FES[1], FES[2], FES[3], FES[4], FES[5]]; kwargs...) |
| 254 | + |
| 255 | + ## correct integral mean of pressures |
| 256 | + pmean_exact = integrate(xgrid, ON_CELLS, p!, 1; quadorder = 10, regions = [rFF]) |
| 257 | + qmean_exact = integrate(xgrid, ON_CELLS, q!, 1; params = [k], quadorder = 10, regions = [rPM]) |
| 258 | + pintegrate = ItemIntegrator([id(p)]; regions = [rFF]) |
| 259 | + qintegrate = ItemIntegrator([id(q)]; regions = [rPM]) |
| 260 | + pmean = sum(evaluate(pintegrate, sol)) |
| 261 | + qmean = sum(evaluate(qintegrate, sol)) |
| 262 | + if coupled |
| 263 | + view(sol[p]) .-= (pmean + qmean) |
| 264 | + view(sol[q]) .-= (pmean + qmean) |
| 265 | + else |
| 266 | + view(sol[p]) .-= pmean - pmean_exact |
| 267 | + view(sol[q]) .-= qmean - qmean_exact |
| 268 | + end |
| 269 | + |
| 270 | + ## error calculation |
| 271 | + ErrorIntegratorExactU = ItemIntegrator(exact_error!(u!), [id(u)]; bonus_quadorder = 4, regions = [rFF], kwargs...) |
| 272 | + ErrorIntegratorExactV = ItemIntegrator(exact_error!(v!), [id(v)]; params = [k], bonus_quadorder = 4, regions = [rPM], kwargs...) |
| 273 | + ErrorIntegratorExactP = ItemIntegrator(exact_error!(p!), [id(p)]; bonus_quadorder = 4, regions = [rFF], kwargs...) |
| 274 | + ErrorIntegratorExactQ = ItemIntegrator(exact_error!(q!), [id(q)]; params = [k], bonus_quadorder = 4, regions = [rPM], kwargs...) |
| 275 | + errorU = evaluate(ErrorIntegratorExactU, sol) |
| 276 | + errorV = evaluate(ErrorIntegratorExactV, sol) |
| 277 | + errorP = evaluate(ErrorIntegratorExactP, sol) |
| 278 | + errorQ = evaluate(ErrorIntegratorExactQ, sol) |
| 279 | + L2errorU = sqrt(sum(view(errorU, 1, :)) + sum(view(errorU, 2, :))) |
| 280 | + L2errorV = sqrt(sum(view(errorV, 1, :)) + sum(view(errorV, 2, :))) |
| 281 | + L2errorP = sqrt(sum(view(errorP, 1, :))) |
| 282 | + L2errorQ = sqrt(sum(view(errorQ, 1, :))) |
| 283 | + @info "L2error(u) = $L2errorU" |
| 284 | + @info "L2error(v) = $L2errorV" |
| 285 | + @info "L2error(p) = $L2errorP" |
| 286 | + @info "L2error(q) = $L2errorQ" |
| 287 | + |
| 288 | + ## plot |
| 289 | + plt = plot([id(u), id(v), id(p), id(q), grid(u)], sol; Plotter = Plotter, ncols = 2, rasterpoints = 12, width = 800, height = 1200) |
| 290 | + |
| 291 | + return [L2errorU, L2errorV, L2errorP, L2errorQ], plt |
| 292 | +end |
| 293 | + |
| 294 | +generateplots = ExtendableFEM.default_generateplots(Example264_StokesDarcy, "example264.png") #hide |
| 295 | +function runtests() #hide |
| 296 | + errors, plt = main(; μ = 1, k = 1, σ = 1, nrefs = 3, coupled = true) #hide |
| 297 | + @info errors |
| 298 | + @test errors[1] <= 0.000804 #hide |
| 299 | + @test errors[2] <= 0.004468 #hide |
| 300 | + @test errors[3] <= 0.043904 #hide |
| 301 | + @test errors[4] <= 0.002966 #hide |
| 302 | + return nothing #hide |
| 303 | +end #hide |
| 304 | +end # module |
0 commit comments