Some multilevel decoupled algorithms for a mixed Navier-Stokes/Darcy model

TL;DR

This paper introduces multilevel decoupled algorithms for a mixed Navier-Stokes/Darcy model, demonstrating their accuracy and efficiency through error analysis and numerical experiments with different discretizations.

Contribution

It proposes new multilevel decoupled algorithms for the Navier-Stokes/Darcy model, with proven error estimates and validated effectiveness.

Findings

Algorithms achieve accurate approximations as shown by error estimates.

Numerical experiments confirm the effectiveness of the proposed methods.

Both first and second order discretizations are successfully implemented.

Abstract

In this work, several multilevel decoupled algorithms are proposed for a mixed Navier-Stokes/Darcy model. These algorithms are based on either successively or parallelly solving two linear subdomain problems after solving a coupled nonlinear coarse grid problem. Error estimates are given to demonstrate the approximation accuracy of the algorithms. Experiments based on both the first order and the second order discretizations are presented to show the effectiveness of the decoupled algorithms.

Tables

Table 1. Table 1: Tests for the two-level algorithms under the Mini/ P 1 subscript 𝑃 1 P_{1} element discretizations. The FE solution errors, the intermediate-step two-level solution errors, and the final two-level solution errors with the scaling between the two-level sizes being set as h = H 3 ℎ superscript 𝐻 3 h=H^{3} .

$h=H^3$	$e_0^{{\varphi }_h}$	$e_1^{{\varphi }_h}$	$e_0^{u_h}$	$e_1^{u_h}$	$e_0^{v_h}$	$e_1^{v_h}$	$e_0^{p_h}$
$1/8$	1.736E-3	6.134E-2	3.685E-3	1.263E-1	2.588E-3	1.066E-1	7.420E-2
$1/27$	1.552E-4	1.823E-2	3.213E-4	3.714E-2	2.251E-4	3.070E-2	9.113E-3
$1/64$	2.766E-5	7.693E-3	5.697E-5	1.564E-2	3.996E-5	1.289E-2	2.255E-3
$1/125$	7.253E-6	3.939E-3	1.491E-5	8.000E-3	1.046E-5	6.587E-3	7.906E-4
$h=H^3$	$e_{0,A}^{{\varphi }_h^{\ast }}$	$e_{1,A}^{{\varphi }_h^{\ast }}$	$e_{0,A}^{u_h^{\ast }}$	$e_{1,A}^{u_h^{\ast }}$	$e_{0,A}^{v_h^{\ast }}$	$e_{1,A}^{v_h^{\ast }}$	$e_{0,A}^{p_h^{\ast }}$
$1/8$	7.649E-3	6.848E-2	3.830E-3	1.263E-1	2.573E-3	1.067E-1	7.541E-2
$1/27$	3.051E-3	2.276E-2	4.089E-4	3.718E-2	2.392E-4	3.072E-2	1.355E-2
$1/64$	1.646E-3	1.073E-2	9.572E-5	1.567E-2	6.792E-5	1.290E-2	5.606E-3
$1/125$	1.049E-3	6.206E-3	4.773E-5	8.024E-3	3.034E-5	6.600E-3	3.026E-3
$h=H^3$	$e_{0,A}^{{\varphi }^h}$	$e_{1,A}^{{\varphi }^h}$	$e_{0,A}^{u^h}$	$e_{1,A}^{u^h}$	$e_{0,A}^{v^h}$	$e_{1,A}^{v^h}$	$e_{0,A}^{p^h}$
$1/8$	1.741E-3	6.134E-2	3.685E-3	1.263E-1	2.588E-3	1.066E-1	7.421E-2
$1/27$	1.580E-4	1.823E-2	2.891E-4	3.715E-2	2.251E-4	3.070E-2	9.194E-3
$1/64$	2.766E-5	7.693E-3	5.156E-5	1.564E-2	3.988E-5	1.289E-2	2.262E-3
$1/125$	6.686E-6	3.939E-3	1.361E-5	8.000E-3	1.044E-5	6.587E-3	7.919E-4
$h=H^3$	$e_{0,B}^{{\varphi }_h^{\ast }}$	$e_{1,B}^{{\varphi }_h^{\ast }}$	$e_{0,B}^{u_h^{\ast }}$	$e_{1,B}^{u_h^{\ast }}$	$e_{0,B}^{v_h^{\ast }}$	$e_{1,B}^{v_h^{\ast }}$	$e_{0,B}^{p_h^{\ast }}$
$1/8$	1.725E-3	6.134E-2	3.692E-3	1.263E-1	2.573E-3	1.067E-1	9.835E-2
$1/27$	1.476E-4	1.823E-2	4.163E-4	3.716E-2	2.633E-4	3.075E-2	3.286E-2
$1/64$	2.195E-5	7.693E-3	1.206E-4	1.567E-2	1.143E-4	1.293E-2	1.766E-2
$1/125$	6.615E-6	3.939E-3	7.639E-5	8.031E-3	8.135E-5	6.634E-3	1.111E-3
$h=H^3$	$e_{0,B}^{{\varphi }^h}$	$e_{1,B}^{{\varphi }^h}$	$e_{0,B}^{u^h}$	$e_{1,B}^{u^h}$	$e_{0,B}^{v^h}$	$e_{1,B}^{v^h}$	$e_{0,B}^{p^h}$
$1/8$	1.736E-3	6.134E-2	3.685E-3	1.263E-1	2.588E-3	1.066E-1	7.418E-2
$1/27$	1.553E-4	1.823E-2	3.225E-4	3.714E-2	2.252E-4	3.070E-2	9.110E-3
$1/64$	2.773E-5	7.693E-3	5.767E-5	1.564E-2	3.998E-5	1.289E-2	2.248E-3
$1/125$	7.319E-6	3.939E-3	1.363E-5	8.000E-3	1.052E-5	6.587E-3	7.879E-4
$h=H^3$	$e_{0,C}^{{\varphi }^h}$	$e_{1,C}^{{\varphi }^h}$	$e_{0,C}^{u^h}$	$e_{1,C}^{u^h}$	$e_{0,C}^{v^h}$	$e_{1,C}^{v^h}$	$e_{0,B}^{p^h}$
$1/8$	7.649E-3	6.848E-2	3.692E-3	1.263E-1	2.573E-3	1.067E-1	9.835E-2
$1/27$	3.051E-3	2.276E-2	4.163E-4	3.716E-2	2.633E-4	3.075E-2	3.286E-2
$1/64$	1.646E-3	1.073E-2	1.206E-4	1.567E-2	1.143E-4	1.293E-2	1.766E-2
$1/125$	1.049E-3	6.206E-3	7.639E-5	8.031E-3	8.135E-5	6.634E-3	1.111E-3

Table 2. Table 2: Tests for the two-level algorithms under the Taylor-Hood/ P 2 subscript 𝑃 2 P_{2} element discretizations. The FE solution errors, the intermediate-step two-level solution errors, and the final two-level solution errors with the scaling between the two-level mesh sizes being set as h ≈ H 5 / 2 ℎ superscript 𝐻 5 2 h\approx H^{5/2} .

$h\approx H^{5/2},H$	$e_0^{{\varphi }_h}$	$e_1^{{\varphi }_h}$	$e_0^{u_h}$	$e_1^{u_h}$	$e_0^{v_h}$	$e_1^{v_h}$	$e_0^{p_h}$
$1/6,1/2$	1.056E-4	2.522E-3	4.017E-4	8.270E-3	2.266E-4	4.882E-3	2.837E-3
$1/16,1/3$	5.584E-6	3.648E-4	2.221E-5	1.156E-3	1.161E-5	6.732E-4	2.930E-4
$1/32,1/4$	7.000E-7	9.201E-5	2.797E-6	2.888E-4	1.446E-6	1.678E-4	7.002E-5
$1/56,1/5$	1.308E-7	3.016E-5	5.228E-7	9.432E-5	2.696E-7	5.477E-5	2.263E-5
$h\approx H^{5/2},H$	$e_{0,A}^{{\varphi }_h^{\ast }}$	$e_{1,A}^{{\varphi }_h^{\ast }}$	$e_{0,A}^{u_h^{\ast }}$	$e_{1,A}^{u_h^{\ast }}$	$e_{0,A}^{v_h^{\ast }}$	$e_{1,A}^{v_h^{\ast }}$	$e_{0,A}^{p_h^{\ast }}$
$1/6,1/2$	1.808E-4	2.830E-3	4.006E-4	8.270E-3	2.275E-4	4.883E-3	2.849E-3
$1/16,1/3$	3.077E-5	4.666E-4	2.233E-5	1.156E-3	1.212E-5	6.736E-4	2.947E-4
$1/32,1/4$	8.907E-6	1.317E-4	2.955E-6	2.889E-4	1.794E-6	1.680E-4	7.060E-5
$1/56,1/5$	3.484E-6	5.018E-5	6.489E-7	9.433E-5	4.811E-7	5.485E-5	2.300E-5
$h\approx H^{5/2},H$	$e_{0,A}^{{\varphi }^h}$	$e_{1,A}^{{\varphi }^h}$	$e_{0,A}^{u^h}$	$e_{1,A}^{u^h}$	$e_{0,A}^{v^h}$	$e_{1,A}^{v^h}$	$e_{0,A}^{p^h}$
$1/6,1/2$	1.055E-4	2.522E-3	4.018E-4	8.270E-3	2.266E-4	4.882E-3	2.837E-3
$1/16,1/3$	5.587E-5	3.648E-4	2.221E-5	1.156E-3	1.161E-5	6.732E-4	2.930E-4
$1/32,1/4$	7.087E-7	9.202E-5	2.797E-6	2.888E-4	1.446E-6	1.678E-4	7.002E-5
$1/56,1/5$	1.380E-7	3.016E-5	5.229E-7	9.432E-5	2.696E-7	5.477E-5	2.263E-5
$h\approx H^{5/2},H$	$e_{0,B}^{{\varphi }_h^{\ast }}$	$e_{1,B}^{{\varphi }_h^{\ast }}$	$e_{0,B}^{u_h^{\ast }}$	$e_{1,B}^{u_h^{\ast }}$	$e_{0,B}^{v_h^{\ast }}$	$e_{1,B}^{v_h^{\ast }}$	$e_{0,B}^{p_h^{\ast }}$
$1/6,1/2$	1.052E-4	2.523E-3	4.038E-4	8.277E-3	2.530E-4	4.949E-3	3.018E-3
$1/16,1/3$	6.126E-6	3.656E-4	3.195E-5	1.160E-3	2.663E-5	6.978E-4	3.472E-4
$1/32,1/4$	1.020E-6	9.227E-5	7.704E-6	2.907E-4	7.129E-6	1.780E-4	9.291E-5
$1/56,1/5$	3.095E-7	3.027E-5	2.916E-6	9.535E-5	2.697E-6	6.035E-5	3.489E-5
$h\approx H^{5/2},H$	$e_{0,B}^{{\varphi }^h}$	$e_{1,B}^{{\varphi }^h}$	$e_{0,B}^{u^h}$	$e_{1,B}^{u^h}$	$e_{0,B}^{v^h}$	$e_{1,B}^{v^h}$	$e_{0,B}^{p^h}$
$1/6,1/2$	1.057E-4	2.522E-3	4.019E-4	8.270E-3	2.266E-4	4.882E-3	2.836E-3
$1/16,1/3$	5.585E-6	3.648E-4	2.223E-5	1.156E-3	1.162E-5	6.732E-4	2.930E-4
$1/32,1/4$	7.002E-7	9.201E-5	2.800E-6	2.888E-4	1.449E-6	1.678E-4	7.003E-5
$1/56,1/5$	1.309E-7	3.016E-5	5.242E-7	9.432E-5	2.717E-7	5.477E-5	2.264E-5
$h\approx H^{5/2},H$	$e_{0,C}^{{\varphi }^h}$	$e_{1,C}^{{\varphi }^h}$	$e_{0,C}^{u^h}$	$e_{1,C}^{u^h}$	$e_{0,C}^{v^h}$	$e_{1,C}^{v^h}$	$e_{0,C}^{p^h}$
$1/6,1/2$	1.808E-4	2.830E-3	4.038E-3	8.277E-3	2.530E-4	4.949E-3	3.018E-3
$1/16,1/3$	3.077E-5	4.666E-4	3.195E-5	1.160E-3	2.663E-5	6.978E-4	3.472E-4
$1/32,1/4$	8.907E-6	1.317E-4	7.704E-6	2.907E-4	7.129E-6	1.780E-4	9.291E-5
$1/56,1/5$	3.484E-6	5.018E-5	2.916E-6	9.535E-5	2.697E-6	6.035E-5	3.489E-5

Table 3. Table 3: Tests for the multilevel algorithms under the Mini/ P 1 subscript 𝑃 1 P_{1} element discretization. Comparisons of Algorithm A , Algorithm C and Algorithm D .

$h_l=h_{l-1}^2$	$e_0^{{\varphi }_h}$	$e_1^{{\varphi }_h}$	$e_0^{u_h}$	$e_1^{u_h}$	$e_0^{v_h}$	$e_1^{v_h}$	$e_0^{p_h}$
$2^{-1}$	2.153E-2	2.351E-1	5.524E-2	5.087E-1	4.295E-2	5.252E-1	1.024E-0
$2^{-2}$	6.566E-3	1.215E-1	1.467E-2	2.539E-1	1.049E-2	2.261E-1	2.649E-1
$2^{-4}$	4.406E-4	3.075E-2	9.178E-4	6.282E-2	6.429E-4	5.217E-3	2.201E-2
$2^{-8}$	1.729E-6	1.923E-3	3.551E-6	3.905E-3	2.493E-6	3.214E-3	2.628E-4
$h_l=h_{l-1}^2$	$e_{0,A}^{{\varphi }_h^{\ast }}$	$e_{1,A}^{{\varphi }_h^{\ast }}$	$e_{0,A}^{u_h^{\ast }}$	$e_{1,A}^{u_h^{\ast }}$	$e_{0,A}^{v_h^{\ast }}$	$e_{1,A}^{v_h^{\ast }}$	$e_{0,A}^{p_h^{\ast }}$
$2^{-2}$	1.049E-2	1.247E-1	1.474E-2	2.538E-1	1.047E-2	2.261E-1	2.612E-1
$2^{-4}$	1.868E-3	3.163E-2	9.574E-4	6.283E-2	6.481E-4	5.217E-3	2.213E-2
$2^{-8}$	1.031E-4	1.979E-3	3.112E-6	3.905E-3	3.299E-6	3.214E-3	3.692E-4
$h_l=h_{l-1}^2$	$e_{0,A}^{{\varphi }^h}$	$e_{1,A}^{{\varphi }^h}$	$e_{0,A}^{u^h}$	$e_{1,A}^{u^h}$	$e_{0,A}^{v^h}$	$e_{1,A}^{v^h}$	$e_{0,A}^{p^h}$
$2^{-2}$	6.570E-3	1.215E-1	1.467E-2	2.539E-1	1.049E-2	2.261E-1	2.649E-1
$2^{-4}$	4.424E-4	3.075E-2	9.184E-4	6.282E-2	6.429E-4	5.217E-3	2.201E-2
$2^{-8}$	1.657E-6	1.923E-3	3.818E-6	3.905E-3	2.491E-6	3.214E-3	2.692E-4
$h_l=h_{l-1}^2$	$e_{0,D}^{{\varphi }_h^{\ast }}$	$e_{1,D}^{{\varphi }_h^{\ast }}$	$e_{0,D}^{u_h^{\ast }}$	$e_{1,D}^{u_h^{\ast }}$	$e_{0,D}^{v_h^{\ast }}$	$e_{1,D}^{v_h^{\ast }}$	$e_{0,D}^{p_h^{\ast }}$
$2^{-2}$	1.049E-2	1.247E-1	1.474E-2	2.538E-1	1.047E-2	2.261E-1	2.612E-1
$2^{-4}$	1.865E-3	3.163E-2	9.571E-4	6.283E-2	6.481E-4	5.217E-3	2.212E-2
$2^{-8}$	9.986E-5	1.976E-3	3.177E-6	3.905E-3	3.332E-6	3.214E-3	3.507E-4
$h_l=h_{l-1}^2$	$e_{0,D}^{{\varphi }^h}$	$e_{1,D}^{{\varphi }^h}$	$e_{0,D}^{u^h}$	$e_{1,D}^{u^h}$	$e_{0,D}^{v^h}$	$e_{1,D}^{v^h}$	$e_{0,D}^{p^h}$
$2^{-2}$			1.466E-2	2.539E-1	1.049E-2	2.261E-1	2.706E-1
$2^{-4}$			9.101E-4	6.282E-2	6.404E-4	5.217E-3	2.357E-2
$2^{-8}$			3.318E-6	3.905E-3	2.746E-6	3.214E-3	3.933E-4
$h_l=h_{l-1}^2$	$e_{0,C}^{{\varphi }^h}$	$e_{1,C}^{{\varphi }^h}$	$e_{0,C}^{u^h}$	$e_{1,C}^{u^h}$	$e_{0,C}^{v^h}$	$e_{1,C}^{v^h}$	$e_{0,C}^{p^h}$
$2^{-2}$	1.049E-2	1.247E-1	1.466E-2	2.539E-1	1.049E-2	2.261E-1	2.752E-1
$2^{-4}$	1.868E-3	3.163E-2	8.931E-4	6.284E-2	6.520E-4	5.220E-3	3.974E-2
$2^{-8}$	8.642E-5	1.967E-3	2.998E-5	3.910E-3	3.463E-6	3.233E-3	5.265E-4

Table 4. Table 4: Tests for the multilevel algorithms under the Taylor-Hood/ P 2 subscript 𝑃 2 P_{2} element discretization. Comparisons of Algorithm A , Algorithm C and Algorithm D .

$h_l=h_{l-1}^2$	$e_0^{{\varphi }_h}$	$e_1^{{\varphi }_h}$	$e_0^{u_h}$	$e_1^{u_h}$	$e_0^{v_h}$	$e_1^{v_h}$	$e_0^{p_h}$
$2^{-1}$	3.009E-3	2.100E-2	8.877E-3	7.608E-2	6.391E-3	4.943E-2	7.601E-2
$2^{-2}$	3.587E-4	5.552E-3	1.288E-3	1.874E-2	7.934E-4	1.125E-2	8.659E-3
$2^{-4}$	5.584E-6	3.648E-4	2.221E-5	1.156E-3	1.161E-5	6.732E-4	2.930E-4
$2^{-8}$	1.372E-9	1.449E-6	5.482E-9	4.514E-6	2.827E-9	2.620E-6	1.078E-6
$h_l=h_{l-1}^2$	$e_{0,A}^{{\varphi }_h^{\ast }}$	$e_{1,A}^{{\varphi }_h^{\ast }}$	$e_{0,A}^{u_h^{\ast }}$	$e_{1,A}^{u_h^{\ast }}$	$e_{0,A}^{v_h^{\ast }}$	$e_{1,A}^{v_h^{\ast }}$	$e_{0,A}^{p_h^{\ast }}$
$2^{-2}$	3.968E-4	5.698E-3	1.286E-3	1.873E-2	7.940E-4	1.126E-2	8.676E-3
$2^{-4}$	9.394E-6	3.743E-4	2.219E-5	1.156E-3	1.165E-5	6.732E-4	2.932E-4
$2^{-8}$	7.161E-8	1.632E-6	9.735E-9	4.514E-6	9.342E-9	2.621E-6	1.080E-6
$h_l=h_{l-1}^2$	$e_{0,A}^{{\varphi }^h}$	$e_{1,A}^{{\varphi }^h}$	$e_{0,A}^{u^h}$	$e_{1,A}^{u^h}$	$e_{0,A}^{v^h}$	$e_{1,A}^{v^h}$	$e_{0,A}^{p^h}$
$2^{-2}$	3.585E-4	5.552E-3	1.288E-3	1.874E-2	7.934E-4	1.125E-2	8.659E-3
$2^{-4}$	5.582E-6	3.648E-4	2.221E-5	1.156E-3	1.161E-5	6.732E-4	2.930E-4
$2^{-8}$	6.394E-9	1.451E-6	5.487E-9	4.514E-6	2.842E-9	2.620E-6	1.078E-6
$h_l=h_{l-1}^2$	$e_{0,D}^{{\varphi }_h^{\ast }}$	$e_{1,D}^{{\varphi }_h^{\ast }}$	$e_{0,D}^{u_h^{\ast }}$	$e_{1,D}^{u_h^{\ast }}$	$e_{0,D}^{v_h^{\ast }}$	$e_{1,D}^{v_h^{\ast }}$	$e_{0,D}^{p_h^{\ast }}$
$2^{-2}$	3.968E-4	5.698E-3	1.286E-3	1.873E-2	7.940E-4	1.126E-2	8.676E-3
$2^{-4}$	9.394E-6	3.743E-4	2.219E-5	1.156E-3	1.165E-5	6.732E-4	2.932E-4
$2^{-8}$	7.161E-8	1.632E-6	9.735E-9	4.514E-6	9.342E-9	2.621E-6	1.080E-6
$h_l=h_{l-1}^2$	$e_{0,D}^{{\varphi }^h}$	$e_{1,D}^{{\varphi }^h}$	$e_{0,D}^{u^h}$	$e_{1,D}^{u^h}$	$e_{0,D}^{v^h}$	$e_{1,D}^{v^h}$	$e_{0,D}^{p^h}$
$2^{-2}$			1.286E-3	1.873E-2	7.941E-4	1.125E-2	8.681E-3
$2^{-4}$			2.219E-5	1.156E-3	1.165E-5	6.732E-4	2.932E-4
$2^{-8}$			9.731E-9	4.514E-6	9.339E-9	2.621E-6	1.080E-6
$h_l=h_{l-1}^2$	$e_{0,C}^{{\varphi }^h}$	$e_{1,C}^{{\varphi }^h}$	$e_{0,C}^{u^h}$	$e_{1,C}^{u^h}$	$e_{0,C}^{v^h}$	$e_{1,C}^{v^h}$	$e_{0,C}^{p^h}$
$2^{-2}$	3.968E-4	5.698E-3	1.277E-3	1.873E-2	8.033E-4	1.128E-2	8.821E-3
$2^{-4}$	8.278E-6	3.711E-4	2.241E-5	1.156E-3	1.220E-5	6.747E-4	2.967E-4
$2^{-8}$	4.817E-7	4.613E-6	2.700E-7	4.797E-6	3.334E-7	3.887E-6	2.910E-6

Table 5. Table 5: Multilevel tests: Comparisons of Algorithm A , Algorithm C and Algorithm D using Taylor-Hood/P2 discretization. The errors between the exact solution and the solutions of the multilevel algorithms with h l = h l − 1 3 / 2 subscript ℎ 𝑙 superscript subscript ℎ 𝑙 1 3 2 h_{l}=h_{l-1}^{3/2} .

$h_l=h_{l-1}^{3/2}$	$e_0^{{\varphi }_h}$	$e_1^{{\varphi }_h}$	$e_0^{u_h}$	$e_1^{u_h}$	$e_0^{v_h}$	$e_1^{v_h}$	$e_0^{p_h}$
$\frac{1}{2}$	3.009E-3	2.100E-2	8.877E-3	7.608E-2	6.391E-3	4.943E-2	7.601E-2
$\frac{1}{3}$	8.609E-4	9.673E-3	2.903E-3	3.352E-2	1.951E-3	2.060E-2	2.066E-2
$\frac{1}{4}$	1.056E-4	2.522E-3	4.017E-4	8.270E-3	2.266E-4	4.882E-3	2.837E-3
$\frac{1}{16}$	6.774E-6	4.146E-4	2.692E-5	1.315E-3	1.410E-5	7.663E-4	3.361E-4
$\frac{1}{58}$	1.241E-7	2.911E-5	4.957E-7	9.104E-5	2.556E-7	5.287E-5	2.184E-5
$h_l=h_{l-1}^{3/2}$	$e_{0,D}^{{\varphi }_h^{\ast }}$	$e_{1,D}^{{\varphi }_h^{\ast }}$	$e_{0,D}^{u_h^{\ast }}$	$e_{1,D}^{u_h^{\ast }}$	$e_{0,D}^{v_h^{\ast }}$	$e_{1,D}^{v_h^{\ast }}$	$e_{0,D}^{p_h^{\ast }}$
$\frac{1}{3}$	9.054E-4	9.760E-3	2.902E-3	3.352E-2	1.951E-3	2.059E-2	2.066E-2
$\frac{1}{4}$	1.114E-4	2.538E-3	4.014E-4	8.270E-3	2.268E-4	4.882E-3	2.839E-3
$\frac{1}{16}$	6.848E-6	4.150E-4	2.692E-5	1.315E-3	1.411E-5	7.663E-4	3.361E-4
$\frac{1}{58}$	1.270E-7	2.912E-5	4.957E-7	9.104E-5	2.556E-7	5.287E-5	2.184E-5
$h_l=h_{l-1}^{3/2}$	$e_{0,D}^{{\varphi }^h}$	$e_{1,D}^{{\varphi }^h}$	$e_{0,D}^{u^h}$	$e_{1,D}^{u^h}$	$e_{0,D}^{v^h}$	$e_{1,D}^{v^h}$	$e_{0,D}^{p^h}$
$\frac{1}{3}$			2.901E-3	3.352E-2	1.952E-3	2.059E-2	2.069E-2
$\frac{1}{4}$			4.014E-4	8.270E-3	2.268E-4	4.882E-3	2.839E-3
$\frac{1}{16}$			2.692E-5	1.315E-3	1.411E-5	7.663E-4	3.361E-4
$\frac{1}{58}$			4.957E-7	9.104E-5	2.556E-7	5.287E-5	2.184E-5
$h_l=h_{l-1}^{3/2}$	$e_{0,C}^{{\varphi }^h}$	$e_{1,C}^{{\varphi }^h}$	$e_{0,C}^{u^h}$	$e_{1,C}^{u^h}$	$e_{0,C}^{v^h}$	$e_{1,C}^{v^h}$	$e_{0,C}^{p^h}$
$\frac{1}{3}$	9.054E-4	9.760E-3	2.888E-3	3.352E-2	1.957E-3	2.058E-2	2.079E-2
$\frac{1}{4}$	1.084E-4	2.530E-3	4.003E-4	8.270E-3	2.275E-4	4.883E-3	2.854E-3
$\frac{1}{16}$	6.814E-6	4.149E-4	2.691E-5	1.315E-3	1.411E-5	7.663E-4	3.364E-4
$\frac{1}{58}$	1.369E-7	2.912E-5	4.970E-7	9.104E-5	2.601E-7	5.287E-5	2.185E-5

Table 6. Table 6: 3-level test: comparisons of Algorithm A and Algorithm C using Mini/P1 discretization. The errors between the exact solution and the solutions of the multilevel algorithm with h 1 = h 0 3 subscript ℎ 1 superscript subscript ℎ 0 3 h_{1}=h_{0}^{3} and h 2 = h 1 2 subscript ℎ 2 superscript subscript ℎ 1 2 h_{2}=h_{1}^{2} .

$h_l$	$e_0^{{\varphi }_h}$	$e_1^{{\varphi }_h}$	$e_0^{u_h}$	$e_1^{u_h}$	$e_0^{v_h}$	$e_1^{v_h}$	$e_0^{p_h}$
$2^{-1}$	2.153E-2	2.351E-1	5.524E-2	5.087E-1	4.295E-2	5.252E-1	1.024E-0
$2^{-3}$	1.736E-3	6.134E-2	3.685E-3	1.263E-1	2.588E-3	1.066E-1	7.420E-2
$2^{-6}$	2.766E-5	7.693E-3	5.697E-5	1.564E-2	3.996E-5	1.289E-2	2.255E-3
$h_l$	$e_{0,A}^{{\varphi }_h^{\ast }}$	$e_{1,A}^{{\varphi }_h^{\ast }}$	$e_{0,A}^{u_h^{\ast }}$	$e_{1,A}^{u_h^{\ast }}$	$e_{0,A}^{v_h^{\ast }}$	$e_{1,A}^{v_h^{\ast }}$	$e_{0,A}^{p_h^{\ast }}$
$2^{-3}$	7.649E-3	6.848E-2	3.830E-3	1.263E-1	2.573E-3	1.067E-1	7.541E-2
$2^{-6}$	4.213E-4	7.915E-3	5.485E-5	1.564E-3	4.054E-5	1.289E-2	2.414E-3
$h_l$	$e_{0,A}^{{\varphi }^h}$	$e_{1,A}^{{\varphi }^h}$	$e_{0,A}^{u^h}$	$e_{1,A}^{u^h}$	$e_{0,A}^{v^h}$	$e_{1,A}^{v^h}$	$e_{0,A}^{p^h}$
$2^{-3}$	1.741E-3	6.134E-2	3.685E-3	1.263E-1	2.588E-3	1.066E-1	7.421E-2
$2^{-6}$	2.760E-5	7.693E-3	5.709E-5	1.564E-3	3.995E-5	1.289E-2	2.255E-3
$h_l$	$e_{0,C}^{{\varphi }^h}$	$e_{1,C}^{{\varphi }^h}$	$e_{0,C}^{u^h}$	$e_{1,C}^{u^h}$	$e_{0,C}^{v^h}$	$e_{1,C}^{v^h}$	$e_{0,C}^{p^h}$
$2^{-3}$	7.649E-3	6.848E-2	3.830E-4	1.263E-1	2.573E-3	1.067E-1	7.541E-2
$2^{-6}$	4.090E-4	7.903E-3	1.094E-4	1.566E-3	1.416E-4	1.297E-2	2.151E-2

Table 7. Table 7: 3-level tests: comparisons of Algorithm A and Algorithm C using Taylor-Hood/ P 2 subscript 𝑃 2 P_{2} discretization. The errors between the exact solution and the solutions of the multilevel algorithm with h 1 ≈ h 0 5 / 2 subscript ℎ 1 superscript subscript ℎ 0 5 2 h_{1}\approx h_{0}^{5/2} and h 2 = h 1 2 subscript ℎ 2 superscript subscript ℎ 1 2 h_{2}=h_{1}^{2} .

$h_l$	$e_0^{{\varphi }_h}$	$e_1^{{\varphi }_h}$	$e_0^{u_h}$	$e_1^{u_h}$	$e_0^{v_h}$	$e_1^{v_h}$	$e_0^{p_h}$
$1/2$	2.153E-2	2.351E-1	5.524E-2	5.087E-1	4.295E-2	5.252E-1	1.024E-0
$1/6$	1.056E-4	2.522E-3	4.017E-4	8.270E-3	2.266E-4	4.882E-3	2.837E-2
$1/36$	4.918E-7	7.277E-5	1.965E-6	2.282E-4	1.015E-6	1.326E-4	5.515E-5
$h_l$	$e_{0,A}^{{\varphi }_h^{\ast }}$	$e_{1,A}^{{\varphi }_h^{\ast }}$	$e_{0,A}^{u_h^{\ast }}$	$e_{1,A}^{u_h^{\ast }}$	$e_{0,A}^{v_h^{\ast }}$	$e_{1,A}^{v_h^{\ast }}$	$e_{0,A}^{p_h^{\ast }}$
$1/6$	1.808E-4	2.830E-3	4.006E-4	8.270E-3	2.275E-4	4.883E-3	2.849E-3
$1/36$	1.300E-6	7.463E-5	1.966E-6	2.282E-4	1.017E-6	1.326E-4	5.521E-5
$h_l$	$e_{0,A}^{{\varphi }^h}$	$e_{1,A}^{{\varphi }^h}$	$e_{0,A}^{u^h}$	$e_{1,A}^{u^h}$	$e_{0,A}^{v^h}$	$e_{1,A}^{v^h}$	$e_{0,A}^{p^h}$
$1/6$	1.055E-4	2.522E-3	4.018E-4	8.270E-3	2.266E-4	4.882E-3	2.837E-3
$1/36$	4.921E-7	7.277E-5	1.965E-6	2.282E-4	1.015E-6	1.326E-4	5.515E-5
$h_l$	$e_{0,C}^{{\varphi }^h}$	$e_{1,C}^{{\varphi }^h}$	$e_{0,C}^{u^h}$	$e_{1,C}^{u^h}$	$e_{0,C}^{v^h}$	$e_{1,C}^{v^h}$	$e_{0,C}^{p^h}$
$1/6$	1.808E-4	2.830E-3	4.006E-4	8.270E-3	2.275E-4	4.883E-3	2.849E-3
$1/36$	1.121E-5	1.225E-4	8.057E-6	2.310E-4	8.614E-6	1.470E-4	8.067E-5

Equations

$$\cases{-\nu\Delta{\bf u}+\nabla p+\rho({\bf u}\cdot\nabla){\bf u}={\bf f}&$\forall{\bf x}\in\Omega_{f}$,\cr{\rm div\,}{\bf u}=0&$\forall{\bf x}\in\Omega_{f}$,}$$

$$-\mathrm{div}\left(\frac{\mathbf{K}}{n}\nabla \varphi \right)=f_p\forall \mathbf{x}\in {\Omega }_p.$$

$${\mathbf{K}}_{ij}={\mathbf{K}}_{ji},\text{mbox}and{\alpha }_1(\mathbf{x},\mathbf{x})\le (\mathbf{K}\mathbf{x},\mathbf{x})\le {\alpha }_2(\mathbf{x},\mathbf{x})\forall \mathbf{x}\in {\Omega }_p.$$

$${\mathbf{u}}_p=-\frac{\mathbf{K}}{n}\nabla \varphi \text{mbox}and{p}_p=\rho g(\varphi -z).$$

$$\left\{\begin{array}[]{ll}{\bf u}\cdot{\bf n}_{f}={\bf u}_{p}\cdot{\bf n}_{f}=-\frac{{\bf K}}{n}\nabla\phi\cdot{\bf n}_{f},\\ -\nu(\nabla{\bf u}{\bf n}_{f})\cdot{\bf n}_{f}+p=\rho g\phi,\\ -\nu(\nabla{\bf u}{\bf n}_{f})\cdot{\boldsymbol{\tau}}_{i}=\frac{\nu\alpha_{{}_{BJS}}}{\sqrt{\nu{\boldsymbol{\tau}}_{i}\cdot{\bf K}{\boldsymbol{\tau}}_{i}}}{\bf u}\cdot{\boldsymbol{\tau}}_{i},\quad i=1,...,d-1.\end{array}\right.$$

$$\cases{{\bf u}={\bf 0}&on $\Gamma_{f}$,\cr\phi=0&on $\Gamma_{D}$.}$$

$${\mathbf{X}}_f$$

$$Q$$

$$X_p$$

$$\left\{\begin{array}[]{ll}a(\underline{u},\underline{v})+c({\bf u},{\bf u},{\bf v})+b({\bf v},p)=f(\underline{v})&\forall\underline{v}=({\bf v},\psi)\in\underline{X},\\ b({\bf u},q)=0&\forall q\in Q,\end{array}\right.$$

$$a(\underline{u},\underline{v})=a_f(\mathbf{u},\mathbf{v})+a_p(\varphi ,\psi )+a_{\Gamma }(\underline{u},\underline{v}),b(\mathbf{v},p)=-{\int }_{{\Omega }_f}p\nabla \cdot \mathbf{v},$$

$$c(\mathbf{u},\mathbf{v},\mathbf{w})=\rho {\int }_{{\Omega }_f}(\mathbf{u}\cdot \nabla )\mathbf{v}\cdot \mathbf{w},f(\underline{v})={\int }_{{\Omega }_f}{\mathbf{f}}_f\cdot \mathbf{v}+\rho g{\int }_{{\Omega }_p}{f}_p\psi$$

$$a_f(\mathbf{u},\mathbf{v})$$

$$a_p(\varphi ,\psi )$$

$$(a)|c(\mathbf{u},\mathbf{v},\mathbf{w})|\lesssim |\mathbf{u}{|}_{1,{\Omega }_f}|\mathbf{v}{|}_{1,{\Omega }_f}|\mathbf{w}{|}_{1,{\Omega }_f}\forall \mathbf{u},\mathbf{v},\mathbf{w}\in {\mathbf{X}}_f,$$

$$(b)|c(\mathbf{u},\mathbf{v},\mathbf{w})|\lesssim |\mathbf{u}{|}_{0,{\Omega }_f}|\mathbf{v}{|}_{1,{\Omega }_f}\Vert \mathbf{w}{\Vert }_{2,{\Omega }_f}\forall \mathbf{u},\mathbf{v}\in {\mathbf{X}}_f,\mathbf{w}\in {\mathbf{H}}^2({\Omega }_f),$$

$$(c)|c(\mathbf{u},\mathbf{v},\mathbf{w})|\lesssim |\mathbf{u}{|}_{1,{\Omega }_f}|\mathbf{v}{|}_{0,{\Omega }_f}\Vert \mathbf{w}{\Vert }_{2,{\Omega }_f}\forall \mathbf{u},\mathbf{v}\in {\mathbf{X}}_f,\mathbf{w}\in {\mathbf{H}}^2({\Omega }_f).$$

$$\left\{\begin{array}[]{ll}a(\underline{u}_{h},\underline{v}_{h})+c({\bf u}_{h},{\bf u}_{h},{\bf v}_{h})+b({\bf v}_{h},p_{h})=f(\underline{v}_{h})&\forall\underline{v}_{h}=({\bf v}_{h},\psi_{h})\in\underline{X}_{h},\\ b({\bf u}_{h},q_{h})=0&\forall q_{h}\in Q_{h}.\end{array}\right.$$

$$\underset{{\mathbf{v}}_h\in {\mathbf{X}}_{f,h}}{\sup }\frac{b({\mathbf{v}}_h,q_h)}{|{\mathbf{v}}_h{|}_{1,{\Omega }_f}}\ge \beta \Vert q_h{\Vert }_{0,{\Omega }_f}\forall q_h\in Q_h.$$

$$\inf_{{\bf v}_{h}\in{\bf X}_{f,h},q_{h}\in Q_{h}}\big{\{}h|{\bf u}-{\bf v}_{h}|_{1,\Omega_{f}}+\|{\bf u}-{\bf v}_{h}\|_{0,\Omega_{f}}+h\|p-q_{h}\|_{0,\Omega_{f}}\big{\}}\lesssim h^{k+1}(|{\bf u}|_{k+1,\Omega_{f}}+|p|_{k,\Omega_{f}});$$

$$\inf_{\psi_{h}\in X_{p,h}}\big{\{}h|\phi-\psi_{h}|_{1,\Omega_{p}}+\|\phi-\psi_{h}\|_{0,\Omega_{p}}\big{\}}\lesssim h^{k+1}|\phi|_{k+1,\Omega_{p}}.$$

$$|\mathbf{u}-{\mathbf{u}}_h{|}_{1,{\Omega }_f}+|\varphi -{\varphi }_h{|}_{1,{\Omega }_p}+\Vert p-p_h{\Vert }_{0,{\Omega }_f}\lesssim h^k.$$

$$\Vert \mathbf{u}-{\mathbf{u}}_h{\Vert }_{0,{\Omega }_f}+\Vert \varphi -{\varphi }_h{\Vert }_{0,{\Omega }_p}\lesssim h^{k+1}.$$

$$\cases{a(\underline{u}_{H},\underline{v}_{H})+c({\bf u}_{H},{\bf u}_{H},{\bf v}_{H})+b(v_{H},p_{H})=(f,v_{H})&$\forall v_{H}=({\bf v}_{H},\psi_{H})\in{\underline{X}_{H}}$,\cr b(u_{H},q_{H})=0&$\forall q_{H}\in Q_{H}$.}$$

$${\tilde{a}}_f(\mathbf{u},\mathbf{v},\mathbf{w})=a_f(\mathbf{u},\mathbf{v})+c(\mathbf{u},\mathbf{v},\mathbf{w})+c(\mathbf{v},\mathbf{u},\mathbf{w}).$$

$$a_p({\varphi }_{A,h_l}^{\ast },{\psi }_{h_l})=(f_p,{\psi }_{h_l})+\rho g({\psi }_{h_l},{\mathbf{u}}_A^{h_{l-1}}\cdot {\mathbf{n}}_f{)}_{\Gamma }\forall {\psi }_{h_l}\in X_{p,h_l}.$$

$$\text{cases}{\tilde{a}}_f({\mathbf{u}}_A^{h_{l-1}},{\mathbf{u}}_{A,h_l}^{\ast },{\mathbf{v}}_{h_l})+b({\mathbf{v}}_{h_l},p_{A,h_l}^{\ast })=(\tilde{\mathbf{f}},{\mathbf{v}}_{h_l})-\rho g({\varphi }_{A,h_l}^{\ast },{\mathbf{v}}_{h_l}\cdot {\mathbf{n}}_f{)}_{\Gamma }\text{cr}b({\mathbf{u}}_{A,h_l}^{\ast },q_{h_l})=0.$$

$$a_p({\varphi }_A^{h_l},{\psi }_{h_l})=(f_p,{\psi }_{h_l})+\rho g({\psi }_{h_l},{\mathbf{u}}_{A,h_l}^{\ast }\cdot {\mathbf{n}}_f{)}_{\Gamma }\forall {\psi }_{h_l}\in X_{p,h_l}.$$

$$\text{cases}{\tilde{a}}_f({\mathbf{u}}_A^{h_{l-1}},{\mathbf{u}}_A^{h_l},{\mathbf{v}}_{h_l})+b({\mathbf{v}}_{h_l},p_A^{h_l})=(\bar{\mathbf{f}},{\mathbf{v}}_{h_l})-\rho g({\varphi }_A^{h_l},{\mathbf{v}}_{h_l}\cdot {\mathbf{n}}_f{)}_{\Gamma }\text{cr}b({\mathbf{u}}_A^{h_l},q_{h_l})=0.$$

$$\text{cases}{\tilde{a}}_f({\mathbf{u}}_B^{h_{l-1}},{\mathbf{u}}_{B,h_l}^{\ast },{\mathbf{v}}_{h_l})+b({\mathbf{v}}_{h_l},p_{B,h_l}^{\ast })=(\tilde{\mathbf{f}},{\mathbf{v}}_{h_l})-\rho g({\varphi }_B^{h_{l-1}},{\mathbf{v}}_{h_l}\cdot {\mathbf{n}}_f{)}_{\Gamma }\text{cr}b({\mathbf{u}}_{B,h_l}^{\ast },q_{h_l})=0.$$

$$a_p({\varphi }_{B,h_l}^{\ast },{\psi }_{h_l})=(f_p,{\psi }_{h_l})+\rho g({\psi }_{h_l},{\mathbf{u}}_{B,h_j}^{\ast }\cdot {\mathbf{n}}_f{)}_{\Gamma }\forall {\psi }_{h_l}\in X_{p,h_l}.$$

$$\text{cases}{\tilde{a}}_f({\mathbf{u}}_B^{h_{l-1}};{\mathbf{u}}_B^{h_l},{\mathbf{v}}_{h_l})+b({\mathbf{v}}_{h_l},p_B^{h_l})=(\bar{\mathbf{f}},{\mathbf{v}}_{h_l})-\rho g({\varphi }_{B,h_l}^{\ast },{\mathbf{v}}_{h_l}\cdot {\mathbf{n}}_f{)}_{\Gamma }\text{cr}b({\mathbf{u}}_B^{h_l},q_{h_l})=0.$$

$$a_p({\varphi }_B^{h_l},{\psi }_{h_l})=(f_p,{\psi }_{h_l})+\rho g({\psi }_{h_l},{\mathbf{u}}_B^{h_l}\cdot {\mathbf{n}}_f{)}_{\Gamma }\forall {\psi }_{h_l}\in X_{p,h_l}.$$

Full text

{AMS}

65F08, 65F10, 65N30, 65N55

1 Introduction

The coupling of incompressible fluid flow with porous media flow is an interesting but challenging topic. For describing the interactions of the fluid flow with the porous media flow, a coupled Stokes/Darcy or Naiver-Stokes/Darcy system is typically used as a macro-scale sharp interface model [2, 3, 6, 8, 9, 12, 13, 14, 16, 17, 18, 19, 20, 21, 23, 24, 25, 34, 35, 37, 40, 41, 42]. The coupled Navier-Stokes/Darcy model is composed of a nonlinear Navier-Stokes equations for fluid flow, a Darcy law equation for porous media flow, plus certain interface conditions for describing the interactions of the different types of flows. Numerical methods for this model [8, 13, 23, 42] usually result in a coupled and nonlinear saddle point problem, for which numerical difficulties increase as the mesh size decreases.

Let us consider a domain $\Omega\subset R^{d}$ ($d=2$ or 3), consisting of a fluid region $\Omega_{f}$ and a porous media region $\Omega_{p}$ separated by an interface $\Gamma$. As shown in Fig. 1, $\Omega=\Omega_{f}\bigcup\Omega_{p}$ and $\Gamma=\overline{\Omega}_{f}\bigcap\overline{\Omega}_{p}$. The interface $\Gamma$ is assumed to be smooth enough [23].

The fluid flow in $\Omega_{f}$ is governed by the steady state Navier-Stokes equations:

[TABLE]

where $\rho$ is the density of the fluid flow, ${\bf u}$ is the velocity vector, $p$ is the pressure, ${\bf f}$ is the external force, $\nu>0$ is the viscosity coefficient.

In the porous media region $\Omega_{p}$, the governing equation becomes

[TABLE]

Here, $\phi$ is the piezometric head, $f_{p}$ is the source term due to injection or pump, $n$ is the volumetric porosity, ${\bf K}$ is the hydraulic conductivity tensor of the porous media satisfying

[TABLE]

Here, $\alpha_{1}$ and $\alpha_{2}$ are positive constants. Typically, ${\bf K}_{ij}$ is proportional to $\frac{\epsilon^{2}}{\nu}$ with $\epsilon$ being the characteristic length of the porous media. For simplicity, in this paper, we will assume that ${\bf K}=\frac{\epsilon^{2}}{\nu}{\bf I}$. In $\Omega_{p}$, the flow velocity and pressure can be calculated by

[TABLE]

Here, $z$, representing the elevation from a reference level, is assumed to be [math], $p_{p}$ is the pressure in $\Omega_{p}$, and $g$ is the gravity acceleration.

The key part of the coupled model is the transmission conditions at the interface, which describe the interaction mechanism of the two different types of flows. The following interface conditions have been extensively used and studied in the literature [4, 16, 17, 27, 34, 38]:

[TABLE]

Here, ${\bf n}_{f}$ is the unit outward normal directions on $\partial\Omega_{f}$ at $\Gamma$, $\{{\boldsymbol{\tau}}_{i}\}_{i=1}^{d-1}$ is the unit tangent vector on $\Gamma$, $\alpha_{{}_{BJS}}$ is a positive parameter depending on the properties of the porous medium. The first interface condition ensures the mass conservation across $\Gamma$. The second one is the balance of normal forces across the interface. The third condition is well known as Beavers-Joseph-Saffman’s law [4, 38], which states that the slip velocity is proportional to the shear stress along $\Gamma$.

For boundary conditions, without loss of generality, we impose homogeneous Dirichlet boundary conditions on $\Gamma_{f}=\partial\Omega_{f}/\Gamma$ and $\Gamma_{p}=\partial\Omega_{p}/\Gamma$:

[TABLE]

The Finite Element method (FEM) discretization of the coupled Navier-Stokes/Darcy model will result in a coupled nonlinear saddle point problem, which is very difficult to solve. In this work, we are interested in developing decoupled and linearized methods so that they not only allow for easy and efficient implementation and software reuse, but also are numerically effective and efficient. We propose and investigate four multilevel decoupled algorithms. In all these algorithms, the coupled nonlinear system only needs to be solved on a very coarse grid level. After that, decoupled linearized Navier-Stokes and Darcy subproblems are solved on all the subsequently refined meshes. In Algorithm A, we solve a Darcy subproblem firstly and using the coarse grid solution to provide its boundary condition at the interface, and then solve a linearized Navier-Stokes problem using the Darcy problem to provide its boundary condition at the interface, and finally, on the same fine grid level, correct both the Darcy problem and the Navier-Stokes equations using the most updated subproblems to supplement the boundary conditions to each other at the interface. Algorithm B is similar to Algorithm A. Compared with Algorithm A, we only exchange the order of solving the linearized Navier-Stokes equations and the Darcy problem in Algorithm B [40]. In Algorithm C, on all fine grid levels, we use the previous level solution to provide boundary conditions for each subproblems and solve them in parallel [35, 8, 40, 24]. In Algorithm D, on all fine grid levels, the correction step is only applied to the Navier-Stokes part, the boundary conditions of each subproblem are provided by using the most updated numerical solutions.

These multilevel algorithms are extended from the existing two-level algorithms [26, 8, 35, 15, 24, 41, 29, 30, 44]. However, the error estimates of the multilevel algorithms are much more difficult than those of two-level algorithms. In this paper, a theoretical analysis is given for Algorithm A. We apply mathematical induction method to give the estimates of the multilevel algorithm. Different from other existing papers, in which most of the researchers only analyze and test the first order discretization, our theory is valid not only for the first order discretization, but also valid for a general $k$-th order discretization. In particular, for both the first order and the second order discretizations, it is shown that if the mesh sizes of the two successive mesh levels are scaled with $h_{l}=h_{l-1}^{2}$, $l=1,2,\cdots,L$, then the energy norm errors in the final-step approximation are of optimal order. This means that the final approximation is of the same order of accuracy as the Finite Element approximation to $(\phi,{\bf u},p)$ obtained by solving exactly the coupled nonlinear system on the finest mesh. The results are similar to the so-called mesh independent principle justified for the multilevel algorithm for a single Navier-Stokes model by W. Layton [31, 32]. The advantages of these multilevel algorithms are: they are numerically efficient because they enable the application of the most efficient and optimized local linear solvers on the fine grid that have been well developed for the linearized Navier-Stokes and Darcy models. Furthermore, in this work, we are interested in not only the mathematical analysis, but also the comparisons of different algorithms. Extensive numerical experiments for both the first order and the second order discretizations are provided to compare the different multilevel algorithms and to illustrate the effectiveness of these algorithms. In our numerical experiments, we firstly compare the algorithms in the two-level cases, then careful tests are designed to verify the theoretical predictions; some three-level experiments are also conducted to highlight the possible improvements of the theoretical analysis and the numerical algorithms.

The rest of the paper is organized as follows. The weak problem and a coupled and nonlinear algorithm are introduced in Section 2. Some multilevel algorithms are proposed in Section 3. Numerical analysis for Algorithm A is conducted in Section 4 to show that the decoupled and linearized multilevel algorithm retains the same order of approximation accuracy as the coupled and nonlinear algorithms if the scalings between the successive mesh levels are properly selected. In Section 5, we first compare the proposed two-level algorithms and then investigate the multilevel algorithms.

2 Weak form and finite element approximations

We begin with some notations. Let

[TABLE]

be the functional spaces for ${\bf u}$, $p$ and $\phi$, respectively. We denote $\underline{X}={{\bf X}}_{f}\times X_{p}$. By multiplying test functions to (1) and (2), integrating by parts and plugging in the interface boundary conditions (3)-(4), the weak form of the coupled NS/Darcy model reads as: find $\underline{u}=({\bf u},\phi)\in\underline{X},p\in Q$ such that

[TABLE]

where

[TABLE]

with

[TABLE]

Here, $c({\bf u},{\bf v},{\bf w})$ corresponds to the nonlinear term, $a_{f}({\bf u},{\bf v})+a_{p}(\phi,\psi)+a_{\Gamma}(u,v)$ is the corresponding bilinear form to the linear coupled Stokes/Darcy problem [16, 35]. The following results have been well established: $a(\cdot,\cdot)$ is bounded and coercive; $b(\cdot,\cdot)$ is bounded and satisfies the inf-sup condition [16, 18, 36]; and the nonlinear term satisfies the following estimates [22, 29, 30, 33].

Lemma 2.1

Suppose that the boundary of the domain $\Omega_{f}$ satisfies the strong Lipschitz condition of Adams [1]. We have

[TABLE]

Here and thereafter, we will use $a\lesssim b$ to denote that there exists a generic constant $C$, such that $a\leq Cb$. For the wellposedness of the coupled NS/Darcy model, we refer to [3, 13, 18, 23, 42]. It is shown that the coupled NS/Darcy problem (5) is well-posed if the normal velocity across the interface is sufficiently small and the viscosity $\nu$ is sufficiently large. Moreover, there holds the a priori bound of the weak solution [23, 42].

Now, we discuss the Finite Element approximations of problem (5). For a subdomain $\Omega_{d}$, we denote $\|\cdot\|_{s,\Omega_{d}}$ and $|\cdot|_{s,\Omega_{d}}$ as the usual Sobolev norm and seminorm for $H^{s}(\Omega_{d})$, respectively [1]. $(\cdot,\cdot)_{\Omega_{d}}$ represents the $L^{2}$ inner product on $\Omega_{d}$, where $\Omega_{d}$ can be the interface $\Gamma$ or one of the subdomains. We partition $\Omega_{f}$ and $\Omega_{p}$ by quasi-uniform regular triangulations ${\cal T}_{f,h}$ and ${\cal T}_{p,h}$ with a characteristic meshsize $h$. Moreover, we assume that the two subdomain triangulations coincide at $\Gamma$. If a conventional conforming Finite Element method is applied to the model problem (5), the discrete problem reads as: Find $\underline{u}_{h}=({\bf u}_{h},\phi_{h})\in\underline{X}_{h}={{\bf X}}_{f,h}\times X_{p,h}\subset{{\bf X}}_{f}\times X_{p},p_{h}\in Q_{h}\subset Q$ such that

[TABLE]

Here, the FE pairs ${{\bf X}}_{f,h}\times Q_{h}$ needs to be stable [5, 22], i.e., there exists a positive constant $\beta$ such that

[TABLE]

We assume that the solution of (5) is smooth enough and the FE spaces have the following typical approximation properties: let $k$ be a natural number, for all $({\bf u},p)\in{\bf H}^{k+1}(\Omega_{f})\cap{\bf X}_{f}\times H^{k}(\Omega_{f})$ and $\phi\in H^{k+1}(\Omega_{p})\cap X_{p}$,

[TABLE]

[TABLE]

There are several well-known Finite Element spaces satisfying the discrete inf-sup condition and the approximation properties (7)-(9). For instance, if $k=1$, one can apply the Mini elements [5, 22] in $\Omega_{f}$ and the piecewise linear elements in $\Omega_{p}$. If $k\geq 2$, the $k$-th order Taylor-Hood elements [5, 22, 39] and $P_{k}$ elements can be applied in $\Omega_{f}$ and $\Omega_{p}$ respectively [5, 7, 26]. For simplicity, we will only consider the cases $k=1$ and $k=2$ for numerical experiments in this paper.

For the coupled discrete problem (6), the energy norm error estimates can be derived by using a fixed-point framework [8, 22], the $L^{2}$ error analysis can be obtained by using the Aubin-Nitsche duality argument [8]. In summary, we have

Lemma 2.2

Let $({\bf u},\phi,p)\in{\bf H}^{k+1}(\Omega_{f})\times H^{k+1}(\Omega_{p})\times H^{k}(\Omega_{f})$ be the solution of the NS/Darcy model (5) and $({\bf u}_{h},\phi_{h},p_{h})$ be the Finite Element solution of (6). Assuming that $\nu$ is sufficiently large and $h$ is sufficiently small, there holds the following energy norm estimate for the problem (6).

[TABLE]

Moreover, we have the following $L^{2}$ error estimate:

[TABLE]

Furthermore, one can derive the a priori bound of the FE solution under the assumption that the viscosity $\nu$ is sufficiently large [23, 42].

3 Multilevel decoupled algorithms

In this section, we introduce four multilevel decoupled algorithms for the coupled Navier-Stokes/Darcy model. In the first step of all these algorithms, we solve the coupled nonlinear problem $(\ref{Discrete-Pro})$ on a coarse mesh level: find $u_{H}=({\bf u}_{H},\phi_{H})\in{\underline{X}_{H}}\subset{\underline{X}},~{}p_{H}\in Q_{H}\subset Q$ such that

[TABLE]

In the following, for the ease of notations, we denote

[TABLE]

The first multi-level algorithm is actually an extension of the two-level algorithm developed in [26]. After solving the nonlinear coupled problem on a coarse grid level (cf. equation (12)), the fine-level steps read as:

Algorithm A

1:Set $h_{0}=H$, $\phi^{h_{0}}_{A}=\phi_{H}$, ${\bf u}^{h_{0}}_{A}={\bf u}_{H}$, and $p^{h_{0}}_{A}=p_{H}$.

2:For $l=1,2,...,L$, Step a: Solve a Darcy problem on a fine grid: find $\phi^{*}_{A,h_{l}}\in X_{p,{h_{l}}}\supset X_{p,h_{l-1}}$ such that

[TABLE]

Step b: Solve a linearized NS problem on a fine grid: find ${\bf u}^{*}_{A,h_{l}}\in{{\bf X}}_{f,h_{l}}\supset{{\bf X}}_{f,h_{l-1}},~{}p^{*}_{A,h_{l}}\in Q_{h_{l}}\supset Q_{h_{l-1}}$ such that $\forall{\bf v}_{h_{l}}\in{{\bf X}}_{f,h_{l}}$, $\forall q_{h_{l}}\in Q_{h_{l}}$,

[TABLE]

Here, $(\tilde{\bf f},{\bf v}_{h_{l}})=({\bf f},{\bf v}_{h_{l}})+c({\bf u}^{h_{l-1}}_{A},{\bf u}^{h_{l-1}}_{A},{\bf v}_{h_{l}})$. Step c: Correct the Darcy solution on the same fine grid: find $\phi^{h_{l}}_{A}\in X_{p,{h_{l}}}$ such that

[TABLE]

Step d: Correct the NS solution on the same fine grid: find ${\bf u}^{h_{l}}_{A}\in{{\bf X}}_{f,h_{l}}\supset{{\bf X}}_{f,h_{l-1}},p^{h_{l}}_{A}\in Q_{h_{l}}\supset Q_{h_{l-1}}$ such that $\forall{\bf v}_{h_{l}}\in{{\bf X}}_{f,h_{l}}$, $\forall q_{h_{l}}\in Q_{h_{l}}$,

[TABLE]

    Here, $(\bar{\bf f},{\bf v}_{h_{l}})=({\bf f},{\bf v}_{h_{l}})+c({\bf u}^{h_{l-1}}_{A},{\bf u}^{*}_{A,h_{l}},{\bf v}_{h_{l}})+c({\bf u}^{*}_{A,h_{l}},{\bf u}^{h_{l-1}}_{A}-{\bf u}^{*}_{A,h_{l}},{\bf v}_{h_{l}}).$

3:End

In the second multi-level algorithm, different from Algorithm A, we exchange the order of solving the two subproblems on fine grid levels [41]. Specifically, after solving the coupled nonlinear problem on a coarse grid level, the fine-level steps of the second multi-level algorithm read as:

Algorithm B

1:Set $h_{0}=H$, $\phi^{h_{0}}_{B}=\phi_{H}$, ${\bf u}^{h_{0}}_{B}={\bf u}_{H}$, and $p^{h_{0}}_{B}=p_{H}$.

2:For $l=1,...,L$. Step a: Solve a linearized NS problem on a fine grid: find ${\bf u}^{*}_{B,h_{l}}\in{{\bf X}}_{f,h_{l}}\supset{{\bf X}}_{f,h_{l-1}},~{}p^{*}_{B,h_{l}}\in Q_{h_{l}}\supset Q_{h_{l-1}}$ such that

[TABLE]

Here, $(\tilde{\bf f},{\bf v}_{h_{l}})=({\bf f},{\bf v}_{h_{l}})+c({\bf u}_{B}^{h_{l-1}},{\bf u}_{B}^{h_{l-1}},{\bf v}_{h_{l}})$. Step b: Solve a Darcy problem on a fine grid: find $\phi^{*}_{B,h_{l}}\in X_{p,{h_{l}}}$ such that

[TABLE]

Step c: Correct the NS solution on the same fine grid: find ${\bf u}^{h_{l}}_{B}\in{{\bf X}}_{f,h_{l}}\supset{{\bf X}}_{f,h_{l-1}},p^{h_{l}}_{B}\in Q_{h_{l}}\supset Q_{h_{l-1}}$ such that

[TABLE]

Here, $(\bar{\bf f},{\bf v}_{h_{l}})=({\bf f},{\bf v}_{h_{l}})+c({\bf u}^{h_{l-1}}_{B},{\bf u}^{*}_{B,h_{l}},{\bf v}_{h_{l}})+c({\bf u}^{*}_{B,h_{l}},{\bf u}^{h_{l-1}}_{B}-{\bf u}^{*}_{B,h_{l}},{\bf v}_{h_{l}}).$ Step d: Correct the Darcy solution on the same fine grid: find $\phi^{h_{l}}_{B}\in X_{p,{h_{l}}}$ such that

[TABLE]

3:End

In the third multilevel algorithm, after solving the coupled nonlinear problem on a coarse grid, we will solve the two subproblems in parallel on all fine grid levels. Specifically, the fine-level steps read as:

Algorithm C

1:Set $h_{0}=H$, $\phi^{h_{0}}_{C}=\phi_{H}$, ${\bf u}^{h_{0}}_{C}={\bf u}_{H}$, and $p^{h_{0}}_{C}=p_{H}$.

2:For $l=1,...,L$, Solve a linearized NS equation: find ${\bf u}^{h_{l}}_{C}\in{{\bf X}}_{f,h_{l}}\supset{{\bf X}}_{f,h_{l-1}},~{}p^{h_{l}}_{C}\in Q_{h_{l}}\supset Q_{h_{l-1}}$ such that $\forall{\bf v}_{h_{l}}\in{{\bf X}}_{f,h_{l}}$, $\forall q_{h_{l}}\in Q_{h_{l}}$,

[TABLE]

Solve the local Darcy problem on a fine grid: find $\phi^{h_{l}}_{C}\in X_{p,{h_{l}}}$ such that

[TABLE]

3:End

In the last multi-level algorithm, we skip the correction step for the Darcy problem. After solving the coupled nonlinear problem on a coarse grid level, the fine-level steps of the algorithm reads as:

Algorithm D

1:Set $h_{0}=H$, $\phi^{h_{0}}_{D}=\phi_{H}$, ${\bf u}^{h_{0}}_{D}={\bf u}_{H}$ and $p^{h_{0}}_{D}=p_{H}$.

2:For $l=1,...,L$, Step a: Solve the Darcy problem on a fine grid: find $\phi^{*}_{D,h_{l}}\in X_{p,{h_{l}}}\supset X_{p,h_{l-1}}$ such that

[TABLE]

Step b: Solve a linearized NS problem on a fine grid: find ${\bf u}^{*}_{D,h_{l}}\in{{\bf X}}_{f,h_{l}}\supset{{\bf X}}_{f,h_{l-1}},~{}p^{*}_{D,h_{l}}\in Q_{h_{l}}\supset Q_{h_{l-1}}$ such that $\forall{\bf v}_{h_{l}}\in{{\bf X}}_{f,h_{l}}$, $\forall q_{h_{l}}\in Q_{h_{l}}$,

[TABLE]

Here, $(\tilde{\bf f},{\bf v}_{h_{l}})=({\bf f},{\bf v}_{h_{l}})+c({\bf u}^{h_{l-1}}_{D},{\bf u}^{h_{l-1}}_{D},{\bf v}_{h_{l}})$. Step c: Correct the NS solution on the same fine grid: find ${\bf u}^{h_{l}}_{D}\in{{\bf X}}_{f,h_{l}}\supset{{\bf X}}_{f,h_{l-1}},p^{h_{l}}_{D}\in Q_{h_{l}}\supset Q_{h_{l-1}}$ such that $\forall{\bf v}_{h_{l}}\in{{\bf X}}_{f,h_{l}}$, $\forall q_{h_{l}}\in Q_{h_{l}}$,

[TABLE]

Here, $(\bar{\bf f},{\bf v}_{h_{l}})=({\bf f},{\bf v}_{h_{l}})+c({\bf u}^{h_{l-1}}_{D},{\bf u}^{*}_{D,h_{l}},{\bf v}_{h_{l}})+c({\bf u}^{*}_{D,h_{l}},{\bf u}^{h_{l-1}}_{A}-{\bf u}^{*}_{D,h_{l}},{\bf v}_{h_{l}})$.

3:End

We see that when $L=1$, Algorithm A is reduced to the two-level algorithm developed in [26], Algorithm C degenerates to the two-level algorithm proposed in [35, 8]. Algorithm B is an extension of the two-grid algorithm proposed in [40]. Algorithm D differs from Algorithm A in that there is no correction step for the Darcy problem. Intuitively, each of the above multilevel algorithms can be thought of as a recursive call of a certain two-level algorithm. Moreover, it is not difficult to see that Algorithm A and Algorithm B require more operation cost, while Algorithm C requires the least operation cost on every mesh level.

4 Theoretical Analysis

In this section, we only analyze the solution by the decoupled multilevel Algorithm A. For Algorithm A and Algorithm B, we will see that they produce almost the same accurate solution from our numerical experiments in Section 5. The analysis of Algorithm C in the linear case can be found in [10]. As previously pointed out, when $L=1$, all the above algorithms degenerate to the two level algorithms. We firstly present the results for $L=1$ case, and then provide the error analysis for analyzing the numerical solution on a general meshlevel $l$.

4.1 Results for the two level algorithms

For Algorithm A in the two-level case, we have the following results [26].

Lemma 4.1

Let $H$ and $h$ be the coarse grid size and the fine grid size, i.e., $h_{0}=H$ and $h_{1}=h$, and let $(\phi,{\bf u},p)$, $(\phi^{*}_{A,h},{\bf u}^{*}_{A,h},p^{*}_{A,h})$, and $(\phi^{h}_{A},{\bf u}^{h}_{A},p^{h}_{A})$ be defined by the problems (5), (13)-(14), and (15)-(16), respectively. Under the assumptions of Lemma 2.2, there holds

[TABLE]

[TABLE]

[TABLE]

[TABLE]

For Algorithm C in the two level case, the corresponding analysis for the linear case can be found in [24]. In short, there holds

[TABLE]

Remarks. For Algorithm A, we comment here that Lemma 4.1 indicates that when $h=H^{3}$ if $k=1$ ($h=H^{5/2}$ if $k=2$), the final-step solution of Algorithm A possesses the same order accuracy as the Finite Element solution in the energy norm. In comparison, for Algorithm C, the theoretical estimates of energy norm errors in (30) suggest that one needs to take the scaling $h=H^{2}$ if $k=1$ ($h=H^{3/2}$ if $k=2$). For Algorithm A, to ensure the final-step solutions have optimal $L^{2}$ norm errors, one has to take the scaling $h=H^{2}$ if $k=1,2$; To ensure the intermediate-step solutions have optimal energy convergence, one has to take $h=H^{\frac{k+1}{k}}$ ($h=H^{3/2}$ if $k=2$); To ensure the intermediate-step solutions have optimal $L^{2}$ norm errors, the scaling between the two grid sizes has to be taken as $h=\max\{H^{\frac{k+1}{k}},H^{\frac{2k+1}{k+1}}\}$ ($h=H^{3/2}$ if $k=1,2$).

4.2 Analysis of the multilevel algorithms

The main purpose in this part is to show that the multilevel decoupled and linearized Algorithm A, with a properly chosen scalings of the two successive meshlevel sizes, is of the same order of approximation accuracy as the coupled and nonlinear algorithm. Note that Algorithm A may be viewed as an approximation to the coupled Finite Element algorithm, we will analyze the difference between the solution by Algorithm A and the solution by using the nonlinear coupled algorithm.

To estimate the $L^{2}$-error of the intermediate-step solution of Algorithm A on the $l$-th ($l\geq 1$) mesh level, we will consider the following the dual problem of the linearized problem: given ${\bf g}\in{\bf L}^{2}(\Omega_{f})$, find $({\bf w},r)\in{\bf X}_{f}\times Q$ such that $\forall({\bf v},q)\in{\bf X}_{f}\times Q$

[TABLE]

If the solution of the linearized coupled NS/Darcy model has the regularity $({\bf u},\phi)\in(H^{2}(\Omega_{f}))^{d}\times H^{2}(\Omega_{p})$ as assumed in Lemma 2.2 and for $\nu$ sufficiently large, the two convection terms $c({\bf u},{\bf v},{\bf w})$ and $c({\bf v},{\bf u},{\bf w})$ in the linear dual problem (31) can be properly bounded, and thus we may assume that the solution of the problem (31) is locally smooth and has the regularity

[TABLE]

Theorem 4.1

Let $(\phi,{\bf u},p)$, $(\phi^{*}_{A,h_{l}},{\bf u}^{*}_{A,h_{l}},p^{*}_{A,h_{l}})$ and $(\phi^{*}_{A,h_{l}},{\bf u}^{*}_{A,h_{l}},p^{*}_{A,h_{l}},\phi^{h_{l}}_{A},{\bf u}^{h_{l}}_{A},p^{h_{l}}_{A})$ be defined by problem (5) and (13)-(16) (on a fine grid level with the grid size $h_{l}$), respectively. Under the assumptions of Lemma 2.2, the following error estimates hold:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

**Proof ** (i). The proof of (33) is very similar to the estimate of $\phi$ in the two-grid algorithms developed in [35, 8]. First, by taking $\underline{v}_{h_{l}}=({\bf 0},\phi_{h_{l}}-\phi^{*}_{A,h_{l}})$ in (6) and comparing with the discrete model (23), we have

[TABLE]

Let $\theta\in H^{1}(\Omega_{f})$ be the solution of the problem:

[TABLE]

$\theta$ is the harmonic extension of $\phi_{h_{l}}-\phi^{*}_{A,h_{l}}$ to the fluid flow region and satisfies the following estimate [35].

[TABLE]

Then, integrating by parts and noting that both ${\bf u}_{h_{l}}$ and ${{\bf u}}^{*}_{h_{l}}$ satisfy the discrete divergence-free property, we have, for any $q_{h_{l-1}}\in Q_{h_{l-1}}$,

[TABLE]

By applying the Cauchy-Schwarz inequality and the inequalities (11), there holds

[TABLE]

By applying triangle inequality and the estimate of the Finite Element solution, we see that (33) holds true.

(ii). We only provide a proof for the error estimate of ${\bf u}_{h_{l}}-{\bf u}^{*}_{A,h_{l}}$. Similar to the techniques used in [35, 8, 41, 40], the estimate for $p_{h_{l}}-p^{*}_{A,h_{l}}$ then follows from the discrete inf-sup condition and the estimate of ${\bf u}_{h_{l}}-{\bf u}^{*}_{A,h_{l}}$. To prove (34), we compare the coupled nonlinear discrete problem (6) with the linearized Navier-Stokes model (14). We see that

[TABLE]

Taking ${\bf v}_{h_{l}}={\bf u}_{h_{l}}-{\bf u}^{*}_{A,h_{l}}$, due to the discrete divergence-free property of ${\bf u}_{h_{l}}$ and ${\bf u}^{*}_{A,h_{l}}$, there holds $b({\bf u}_{h_{l}}-{\bf u}^{*}_{A,h_{l}},p_{h_{l}}-p^{*}_{A,h_{l}})=0$. The interface term in (4.2) can be controlled by $|\phi_{h_{l}}-\phi^{*}_{A,h_{l}}|_{1,\Omega_{p}}|{\bf u}_{h_{l}}-{\bf u}^{*}_{A,h_{l}}|_{1,\Omega_{f}}$. For the trilinear terms, it is easy to verify the following identity.

[TABLE]

Note that the viscosity $\nu$ is sufficiently large, then the weak solution, the Finite Element solution as well as the multilevel solution have the a-priori bounds [23, 42]. Then, roughly speaking, the following inequality holds true.

[TABLE]

Thus, by using (4.2), (4.2), and Lemma 2.1, we see that

[TABLE]

Hence, by applying the triangle inequality, $|{\bf u}-{\bf u}^{*}_{A,h_{l}}|_{1,\Omega_{f}}\leq|{\bf u}-{\bf u}_{h_{l}}|_{1,\Omega_{f}}+|{\bf u}_{h_{l}}-{\bf u}^{*}_{A,h_{l}}|_{1,\Omega_{f}}$, and the energy norm estimate of the Finite Element solution (10), we see that the inequality (34) hods true.

(iii). For estimating the $L^{2}$ error of the intermediate-step solution, we set ${\bf g}={\bf u}-{\bf u}^{*}_{A,h_{l}}$ and $({\bf v},q)=({\bf u}-{\bf u}^{*}_{A,h_{l}},p-p^{*}_{A,h_{l}})$ in (31), and then splitting the two trilinear terms into four terms, we obtain

[TABLE]

Taking $\underline{v}=({\bf v}_{h_{l}},0)$ in (5) and subtracting with (14), we obtain

[TABLE]

Subtracting (4.2) from (41), we have

[TABLE]

Here, in the last inequality we have used the Cauchy-Schwarz inequality and the estimate in Lemma 2.1 for the trilinear term, we have dropped those higher order terms as $k\geq 1$ and $h_{l}\leq h_{l-1}$. By the approximation error estimate (8)-(9), discarding the terms which are of the same order or higher order errors (for example, $|\phi-\phi^{*}_{A,h_{l}}|_{1,\Omega_{f}}$ and $\|p-p^{*}_{A,h_{l}}\|_{0,\Omega_{f}}$ are of the same order as $|{\bf u}-{\bf u}^{*}_{A,h_{l}}|_{1,\Omega_{f}}$, and $|{\bf u}-{\bf u}^{*}_{A,h_{l}}|_{1,\Omega_{f}}\leq|{\bf u}-{\bf u}^{h_{l-1}}_{A}|_{1,\Omega_{f}}$), using Lemma 2.2 and the estimate (32), it follows that the $L^{2}$ error estimate for ${\bf u}^{*}_{A,h_{l}}$ holds true.

(iv). The estimate of (36) is similar to the estimate of (33). Taking $\underline{v}_{h}=({\bf 0},\phi_{h_{l}}-\phi^{h_{l}}_{A})$ in (6) and comparing with the discrete model (15), we have

[TABLE]

Let $\theta\in H^{1}(\Omega_{f})$ be a harmonic extension of $\phi_{h_{l}}-\phi^{h_{l}}_{A}$ to the fluid flow region with the Dirichlet data at $\Gamma$ being equal to equal to $\phi_{h_{l}}-\phi^{h_{l}}_{A}$. Then, we have

[TABLE]

Note that for any $q_{h_{l}}\in Q_{h_{l}}$, there holds

[TABLE]

Here, in the last equality, we have used the discrete divergence-free property for ${\bf u}_{h_{l}}$ and ${{\bf u}}^{*}_{A,h_{l}}$. Therefore,

[TABLE]

We see that the estimate (36) holds true.

(v). Now, we estimate the error of ${\bf u}_{h_{l}}-{{\bf u}}^{h_{l}}_{A}$ in the energy norm. Taking $\underline{v}_{h_{l}}=({\bf v}_{h_{l}},0)$ in (6) and comparing with (16) on the $l$-th level mesh, and splitting the trilinear terms in the right hand side, we obtain

[TABLE]

Similar to the proof of (34), letting ${\bf v}_{h_{l}}={\bf u}_{h_{l}}-{{\bf u}}^{h_{l}}_{A}$, we have

[TABLE]

The right hand side of (45) is bounded by

[TABLE]

Applying the triangle inequality, the energy norm error estimate of Finite Element solution (cf. Lemma 2.2), then discarding the terms which are of the same order or higher order errors in (4.2) and (45), we see that (4.1) holds true.

(vi). For estimating the $L^{2}$-error of ${\bf u}-{\bf u}^{h_{l}}_{A}$, let us consider the dual problem of a linearized problem, which is similar to problem (31).

[TABLE]

Moreover, we assume that a regularity estimate which is similar to (32) holds. Setting ${\bf g}={\bf u}-{\bf u}^{h_{l}}_{A}$ and $({\bf v},q)=({\bf u}-{\bf u}^{h_{l}}_{A},p-p^{h_{l}}_{A})$ in (46), splitting the two trilinear terms in (46) into four terms, we have

[TABLE]

Taking $\underline{v}=({\bf v}_{h_{l}},0)$ in (5) and subtracting with (16), we obtain

[TABLE]

Combining (47) and (4.2), we arrive at

[TABLE]

Similar to the proof for $\|{\bf u}-{\bf u}^{*}_{A,h_{l}}\|^{2}_{0,\Omega_{f}}$, there holds

[TABLE]

Then, the estimates for the last three terms in (4.2) are:

[TABLE]

[TABLE]

The estimate of the last term is similar to that for the above term

[TABLE]

Putting all the terms in the right hand side of (4.2) together, and combining with the regularity estimate, discarding the terms which are of the same order or higher order errors, we see that (38) holds true. $\hfill\Box$

We comment here that Theorem 4.1 is valid for a general $k$-th order discretization. However, because of the complex forms of the error terms, it is not easy to identify the scaling relationship for the two adjacent mesh level sizes. As previously mentioned, we are particularly interested in the first order and the second order discretizations. The following theorem states that for the first and the second order discretizations, if $h_{l}=h_{l-1}^{2}$, the energy norm errors of the final-step solution and the $L^{2}$ norm of ${\bf u}^{h_{l}}_{A}$ are of the same orders as those of the FE solution.

Theorem 4.2

Let $(\phi,{\bf u},p)$, $(\phi^{*}_{A,h_{l}},{\bf u}^{*}_{A,h_{l}},p^{*}_{A,h_{l}})$ and $(\phi^{h_{l}}_{A},{\bf u}^{h_{l}}_{A},p^{h_{l}}_{A})$ be defined by the problems (5), (13)-(14), and (15)-(16), respectively. For the first order and the second discretizations, i.e., $k=1$ or $k=2$, if $h_{l}=h_{l-1}^{2}$, under the assumptions of Lemma 2.2, there hold

[TABLE]

[TABLE]

[TABLE]

[TABLE]

**Proof ** We apply mathematical induction to the meshlevel $l$. For proving the results for the solution on meshlevel $l$, we will assume that the conclusions for the solutions on meshlevel $l-1$ hold true. From Lemma 4.1, we know that $L=1$ (in (26)-(29) by changing $H$ to be $h_{0}$ and $h$ to be $h_{1}$), the error estimates for the intermediate-step solution, and the final-step solution hold true. For both $k=1$ and $k=2$, if $h_{1}=h^{2}_{0}$ the estimates (49)-(52) hold true (see also Remark 4.1 in [26]). We are going to prove the results for a general meshlevel $l$. We will discuss the two cases: $k=1$ and $k=2$ separately.

If $k=1$, by using (33)-(38), we see that the estimates for $(\phi^{*}_{A,h_{l}},{\bf u}^{*}_{A,h_{l}},p^{*}_{A,h_{l}})$ and $(\phi^{h_{l}}_{A},{\bf u}^{h_{l}}_{A},p^{h_{l}}_{A})$ are as follows.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

If $k=2$, we see that we see that the estimates for $(\phi^{*}_{A,h_{l}},{\bf u}^{*}_{A,h_{l}},p^{*}_{A,h_{l}})$ and $(\phi^{h_{l}}_{A},{\bf u}^{h_{l}}_{A},p^{h_{l}}_{A})$ are as follows.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

$\hfill\Box$

From Theorem 4.2, under the scaling $h_{l}=h^{2}_{l-1}$, it is shown in (50) that the intermediate-step solution does not have optimal $L^{2}$ errors. To ensure the intermediate-step solution has optimal $L^{2}$ errors, one usually requires a very stringent scaling between the meshsizes of the two subsequent mesh levels. In practice, we are not interested in making the intermediate-step solution has optimal $L^{2}$ error. The estimate of the $L^{2}$ errors is for the purpose of estimating the energy norm of the final-step solution.

We would comment that the theoretical analysis of Algorithm D can be done similar to that for Algorithm A. Noting that there is no correction step in Algorithm D, the scaling of the meshsizes between two adjacent meshlevels are more stringent than that for Algorithm A. In the next section, we provide numerical experiments showing that for the first order discretization, the final-step solution of Algorithm D is still optimal if $h_{l}=h_{l-1}^{2}$. However, for the second order discretization, one has to take $h_{l}=h_{l-1}^{3/2}$ to ensure the final-step solution in the energy norm is optimal (in particular, for the variable $\phi^{h_{l}}_{D}$).

5 Numerical Experiments

We now present numerical experiments to demonstrate the effectiveness and the accuracy of the multi-level approach. In order to make our experiments more solid, we first compare different two-level algorithms then give the numerical experiments for the multilevel cases.

The computational domain is $\Omega\subset\mathbb{R}^{2}$ with $\Omega_{f}=(0,1)\times(1,2)$, $\Omega_{p}=(0,1)\times(0,1)$ and the interface $\Gamma=(0,1)\times\{1\}$. The components of ${\bf u}$ are denoted by $(u,v)$. For simplicity, all the parameters in the coupled NS/Darcy model are set to $1$. The boundary conditions and right hand side functions of the coupled NS/Darcy model are chosen so that the exact solution $(u,v,p,\phi)$ is given by

[TABLE]

The coupled nonlinear FE problem is solved by the Picard iteration: given $(\underline{u}^{0},p^{0})\in\underline{X}_{h}\times Q_{h}$, for $m\geq 0$, find $(\underline{u}^{m+1},p^{m+1})\in\underline{X}_{h}\times Q_{h}$ such that

[TABLE]

The stopping criterion for the Picard iteration is $\|{\bf U}^{m+1}-{\bf U}^{m}\|_{l^{2}}<10^{-7}$, where ${\bf U}^{m}$ is the nodal-value vector for the $m$-th iterate. In all algorithms, the symmetric positive definite linear system of the fine-grid Darcy problems are solved by the PCG method with the incomplete Cholesky factorization as preconditioner. The stopping criterion of PCG is set to be $10^{-9}$ and the dropping tolerance of the incomplete Cholesky factorization is $10^{-3}$. For solving the fine-grid linearized Navier-Stokes problems and the coarse-grid linear system at each step of the Picard iteration, we employ the preconditioned GMRES method with the stopping criterion $\frac{\|{\bf r}_{q}\|_{l^{2}}}{\|{\bf r}_{0}\|_{l^{2}}}<10^{-9}$, where ${\bf r}_{q}$ is the residual at the $q$-th iteration of the GMRES method. The preconditioners of these saddle point problems are designed by using the Green function theory [28]. Interested readers are referred to [7, 11] for more details. All experiments were performed using personal desktop computer with the processor Intel Core i3 2130 (Operating speed 3.4 GHz). For the tests we presented in this paper, the average number of Picard iteration is 6, the number of GMRES iterations for the coupled model or linearized Navier-Stokes model with Green function theory based preconditioner is around 30, the iterations of PCG with incomplete Cholesky factorization preconditioning are less than 280 in all tests.

In the implementation of the two-level and multilevel algorithms, the key part is the Finite Element interpolations. FE interpolations are applied from coarse grid to fine grid or between different submodels. For example, when solving the Darcy problem (13), we need to compute the Neumann data at the quadrature points of the fine grid by using the coarse grid solution. Standard FE interpolation is applied to supplement the Neumann data: take the coarse grid NS solution, and use the coarse grid basis functions to calculate the Neuman data at the quadrature points when assembling the right hand side of (13).

For Algorithm A, the following notations are used to measure the solution errors, the intermediate-step two-level solution errors and the final two-level solution errors for $\phi$ in the energy norm and the $L^{2}$ norm.

[TABLE]

Similarly, the notations, $e^{u^{*}_{h}}_{0,A},e^{u^{*}_{h}}_{1,A},e^{u^{h}}_{0,A},e^{u^{h}}_{1,A},e^{v^{*}_{h}}_{0,A},e^{v^{*}_{h}}_{1,A},e^{v^{h}}_{0,A},e^{v^{h}}_{1,A},e^{p^{*}_{h}}_{0,A}$, and $e^{p^{h}}_{0,A}$ are used to denote the corresponding errors for the velocity components and the pressure variable with specified norms. For the coupled nonlinear FE algorithm, we use $e^{\phi_{h}}_{1},e^{u_{h}}_{1},e^{v_{h}}_{1},e^{p_{h}}_{0}$ to denote the corresponding finite element errors. Similarly, if the algorithm is changed to be Algorithm B, Algorithm C or Algorithm D, the subindex of the errors will be changed correspondingly. We keep $4$ valid digits when calculating all the errors in the following tests.

5.1 Comparisons of the two-level algorithms

We firstly conduct numerical experiments for comparing all algorithms under the two-level cases. According to Lemma 4.1, for the first (second) order discretization, if $h=H^{3}$ ($h=H^{5/2}$), then Algorithm A still gives optimal errors in the energy norm. However, it is not clear whether Algorithm B and Algorithm C also give optimal errors under such a scaling. It is also important to know the differences of the different algorithms so that we can have better understanding of their generalizations to the multilevel cases.

For the first order discretization, Mini elements and piecewise linear elements are used in $\Omega_{f}$ and $\Omega_{p}$, respectively. The scaling between $h$ and $H$ is set to be $h=H^{3}$. In Table 1, we report the numerical results for Algorithm A, Algorithm B and Algorithm C. For Algorithm C, we note that the results for $\phi$ actually are the same as those for $\phi$ in the intermediate-step of Algorithm A and the results for fluid variables actually are the same as those for fluid variable in the intermediate-step solution of Algorithm B. As observed from Table 1, the Finite Element solution errors confirm the theoretical predictions (cf. Lemma 2.2). The pressure FE error is of an order between $O(h)$ and $O(h^{2})$ as there is a bubble function in the Mini element discretization for the fluid part [8, 7, 14]. From Table 1, for both Algorithm A and Algorithm B, we see that the final two-level solution errors in the energy norms, i.e., $e^{\phi^{h}}_{1,A}$, $e^{u^{h}}_{1,A}$, $e^{v^{h}}_{1,A}$, $e^{p^{h}}_{0,A}$, $e^{\phi^{h}}_{1,B}$, $e^{u^{h}}_{1,B}$, $e^{v^{h}}_{1,B}$, and $e^{p^{h}}_{0,B}$ are comparable with those of the coupled algorithm with the same meshsizes; Moreover, the $L^{2}$ errors of $\phi$, $u$ and $v$ are almost of the same order as that for $e^{\phi_{h}}_{0}$, $e^{u_{h}}_{0}$ and $e^{v_{h}}_{0}$. This indicates that the computational $L^{2}$ errors for the velocity components seem to be better than the theoretical predictions (cf. Remark 4.1 in [26]). In comparison, the intermediate-step errors $e^{\phi^{*}_{h}}_{0,A},e^{\phi^{*}_{h}}_{1,A},e^{u^{*}_{h}}_{0,A}$ and $e^{v^{*}_{h}}_{0,A}$ (and also $e^{\phi^{*}_{h}}_{0,B},e^{\phi^{*}_{h}}_{1,B},e^{u^{*}_{h}}_{0,B}$ and $e^{v^{*}_{h}}_{0,B}$) for the intermediate-step two-level solution are not optimal; the errors $e^{u^{*}_{h}}_{1,A},e^{v^{*}_{h}}_{1,A}$ and $e^{p^{*}_{h}}_{0,A}$ (and also $e^{u^{*}_{h}}_{1,B},e^{v^{*}_{h}}_{1,B}$ and $e^{p^{*}_{h}}_{0,B}$) are slightly worse than those corresponding errors for $u_{h}$ or $u^{h}$, $v_{h}$ or $v^{h}$ and $p_{h}$ or $p^{h}$. From the digital comparisons in Table 1, we see that Algorithm B gives almost the same errors as Algorithm A in the final step. However, from Table 5.1 and Figure 2, Algorithm C does not give optimal error order under the same scaling setting for the two level meshsizes (in particular, for the pressure errors).

For second order discretization, Taylor-Hood elements are applied in $\Omega_{f}$ and piecewise quadratic elements are applied in $\Omega_{p}$. The scaling between $h$ and $H$ is set to be $h\approx H^{5/2}$. (Actually, except the case $h=1/32=(1/4)^{5/2}$, the fine grid sizes $h$ are selected even slightly smaller than $H^{5/2}$.) Numerical results are reported in Table 5.2 and the comparisons of Algorithm A and Algorithm C are presented in Figure 3. As observed from Table 5.2 and Figure 3, the Finite Element solution errors confirm the theoretical analysis of Lemma 2.2; The final-step solution errors of the two-level Algorithm A in the energy norms are almost the same as those of the coupled nonlinear FE algorithm with the same meshsizes; Again, Algorithm A and Algorithm B give almost the same numerical solution; The results of Algorithm C actually correspond to the intermediate-step solution of Algorithm A and Algorithm B. From Table 5.2 and Figure 3, Algorithm C does not give optimal numerical errors under the scaling $h\approx H^{5/2}$ (in particular, for pressure errors). By comparing the digital results of Algorithm A and Algorithm C in Table 5.2, one can get the same conclusions as we have drawn for the first order discretization.

5.2 Comparisons of the multilevel algorithms

From Table 1 and Table 5.2, we see that Algorithm A and Algorithm B actually give almost the same numerical accuracy. They only have some difference in the intermediate-step solution errors. This means that it doesn’t matter whether the NS problem or the Darcy problem is solved firstly. Therefore, in the multilevel tests, we will only compare Algorithm A with Algorithm C and Algorithm D.

In Table 5.3 and Table 5.4 we report the numerical results based on the Mini/$P_{1}$ element discretization and the Taylor-Hood/$P_{2}$ element discretization, respectively. Correspondingly, Figure 5.3 and Figure 5.4 plot the results from Table 5.3 and Table 5.4. For the both the first order and the second order discretizations, the scalings of the two successive meshsizes are all set as $h_{l}=h_{l-1}^{2}$. From Table 5.3 and Figure 5.3, we see that all the final-step solution errors based on the decoupled multilevel Algorithm A are almost the same as those based on the coupled nonlinear algorithm. This clearly shows the approximation properties of Algorithm A. For Algorithm C and Algorithm D, although they provide accurate $H^{1}-$ norm errors for the variables $\phi$, $u$, $v$, they can not give accurate pressure errors and the $L^{2}$- norm errors. To be more precisely, Algorithm C can not provide optimal $L^{2}$ norm errors for all variables; Algorithm D can not provide optimal $L^{2}$ norm errors for $\phi$.

For the second order discretization, from the results reported in Table 5.4 and Figure 5.4, we can draw the same conclusions for Algorithm A as those based on the first order discretization. For Algorithm C, both the energy norm errors and the $L^{2}$ norm errors are not accurate enough because the scalings for the two successive mesh sizes are set as $h_{l}=h_{l-1}^{2}$. For Algorithm D, we note that Algorithm D does provide accurate energy norm errors for fluid variables. However, it does not give the optimal errors for $\phi$ (for both the $H^{1}$ and $L^{2}$ norm errors). The reason is that one can not theoretically guarantee the optimal convergence rate of Algorithm D if $h_{l}=h_{l-1}^{2}$ for the second order discretization. Instead, one has to set $h_{l}=h_{l-1}^{3/2}$ for the second order discretization. To verify this, we report the numerical results in Table 5.5. From the results in Table 5.5, we note that both Algorithm C and Algorithm D give almost the same errors as the coupled nonlinear algorithm. The results confirm our theoretical predications, and most importantly, the results suggest that it is necessary to have the correction steps for both the Navier-Stokes subproblem and the Darcy subproblem.

5.3 Experiments for the multilevel algorithms using different scalings between different meshlevel sizes

From [26], we see that for the first two levels of the multilevel Algorithm A, one can take $h_{1}=h_{0}^{3}$ to guarantee the optimal convergence of the energy norm errors (for simplicity, we use the first order discretization for the discussion). However, the analysis in Section 4 shows that one should take $h_{l}=h_{l-1}^{2}$ to guarantee the solution errors are optimal in the energy norm on all mesh levels. This suggests us to test the multilevel algorithms using different scalings on different meshlevels. In this subsection, we test the multilevel algorithms under the three level cases using different scalings between $h_{l}$ and $h_{l-1}$ for two adjacent mesh levels.

For the first order discretization, we take $h_{1}=h_{0}^{3}$ while $h_{2}=h_{1}^{2}$. The corresponding numerical results are reported in Table 5.6. For the second order discretization, we set $h_{1}=h_{0}^{5/2}$ while $h_{2}=h_{1}^{2}$. The corresponding numerical results are reported in Table 5.7. From both Table 5.6 and Table 5.7, by comparing the FE errors with the multilevel algorithm errors, we see that Algorithm A still gives optimal energy norm errors and optimal $L^{2}$- norm errors for velocity. For Algorithm C, as the scaling between $h_{l}$ and $h_{l-1}$ is of very higher order, neither theoretical analysis nor numerical experiments guarantee it can give optimal energy norm or $L^{2}$ norm solution errors.

6 Conclusion

In conclusion, we have proposed some decoupled and linearized multilevel algorithms for the coupled NS/Darcy model. These algorithms are numerically efficient and also enables easy and efficient implementation and software reuse. Numerical analysis are presented to show that the decoupled and linearized Algorithm A retains the same order of approximation accuracy as the coupled and nonlinear algorithm if the scalings between two successive mesh level sizes are properly chosen. Extensive numerical experiments are provided to verify the theoretical predictions and to compare the different algorithms.