Advanced Search
Article Contents

An Adaptive Nonhydrostatic Atmospheric Dynamical Core Using a Multi-Moment Constrained Finite Volume Method


doi: 10.1007/s00376-021-1185-9

  • An adaptive 2D nonhydrostatic dynamical core is proposed by using the multi-moment constrained finite-volume (MCV) scheme and the Berger-Oliger adaptive mesh refinement (AMR) algorithm. The MCV scheme takes several point-wise values within each computational cell as the predicted variables to build high-order schemes based on single-cell reconstruction. Two types of moments, such as the volume-integrated average (VIA) and point value (PV), are defined as constraint conditions to derive the updating formulations of the unknowns, and the constraint condition on VIA guarantees the rigorous conservation of the proposed model. In this study, the MCV scheme is implemented on a height-based, terrain-following grid with variable resolution to solve the nonhydrostatic governing equations of atmospheric dynamics. The AMR grid of Berger-Oliger consists of several groups of blocks with different resolutions, where the MCV model developed on a fixed structured mesh can be used directly. Numerical formulations are designed to implement the coarse-fine interpolation and the flux correction for properly exchanging the solution information among different blocks. Widely used benchmark tests are carried out to evaluate the proposed model. The numerical experiments on uniform and AMR grids indicate that the adaptive model has promising potential for improving computational efficiency without losing accuracy.
    摘要: 本文提出了一种基于多矩限制有限体积法和Berger-Oliger自适应加密网格的二维非静力大气动力框架。多矩限制有限体积法算法在每个计算单元内选择多个点值作为预报量,利用独立单元上的局地高阶重构构造高精度的数值格式。多矩限制有限体积法算法定义了两种类型的矩,即积分平均值矩和点值矩,它们作为限制条件约束预报量的更新,其中对积分平均值矩的限制条件保证了数值格式的严格守恒性。本文在基于高度和地形跟踪的变分辨率网格上采用多矩限制有限体积法算法求解非静力大气运动方程,其中所使用的Berger-Oliger自适应加密网格是由几组具有不同分辨率的均匀网格块构成的,在固定的结构网格上发展的多矩限制有限体积法算法可以直接在这些网格块上使用。我们还设计了粗细网格插值和通量修正的数值方法,使不同层网格块之间的信息得以双向交换。一些广泛使用的标准数值算例被用于检验本文所提出的模式,数值实验的结果表明,在保证数值精度的前提下,自适应模式对于提升计算效率具有可观的前景。
  • 加载中
  • Figure 1.  Distribution of DOFs for the (a) 1D case and (b) 2D Cartesian grid.

    Figure 2.  Diagram depicting the subcycling of AMR levels in time.

    Figure 3.  Flux correction along the interface between the coarse and fine levels.

    Figure 4.  Fine-to-coarse solution transfer in the overlap regions.

    Figure 5.  Coarse-fine interpolations based on multi-moments for evaluating the ghost cells.

    Figure 6.  Coarse-fine interpolations based on multi-moments for evaluating the flow variables in newly generated fine cells.

    Figure 7.  Contour plots of a potential temperature perturbation ($\theta '$) and its errors for a rising thermal bubble on $100 \times 50 \times 3 \times 2$ grid at (a) t = 0 s, (b) t = 250 s, (c) t = 500 s, (d) t = 750 s, (e) t = 1000 s, and (f) absolute errors.

    Figure 8.  Contour plots of the potential temperature perturbation ($\theta '$) and its errors for the density current on $132 \times 32 \times 2 \times 4$ grid at (a) t = 0 s, (b) t = 225 s, (c) t = 450 s, (d) t = 675 s, (e) t = 900 s, and (f) absolute errors.

    Figure 9.  Contour plots of potential temperature perturbation ($\theta '$) for the internal gravity waves on $75 \times 25 \times 3 \times 2$ grid at (a) t = 0 s, (b) t = 750 s, (c) t =1500 s, (d) t = 2250 s, and (e) t = 3000 s.

    Figure 10.  Error distribution of the potential temperature perturbation ($\theta '$) for the internal gravity waves on the (a) uniform grid and (b) $75 \times 25 \times 3 \times 2$ grid at t = 3000 s.

    Figure 11.  Contour plots of horizontal velocity ( u ) of the Schär mountain case at (a) t = 7200 s, (b) t = 16800 s, (c) t = 26400 s, (d) t = 36000 s, and (e) the shape and size of the Schär mountain.

    Table 1.  Normalized errors of the potential temperature perturbation ($\theta '$) and elapsed CPU time for the thermal bubble test running on different grids.

    Resolution (m)Grid configuration$l_1$$l_2$$ {l_\infty }$ CPU time (s)
    200 × 200100 × 50 × 1 × 10.230.210.281 (187.34)
    100 × 100200 × 100 × 1 × 16.40 × 10−26.86 × 10−20.118.45 (1584.21)
    100 × 50 × 2 × 26.44 × 10−26.87 × 10−20.112.84 (532.55)
    50 × 50400 × 200 × 1 × 18.82 × 10−31.09 × 10−22.43 × 10−275.00 (14051.62)
    200 × 100 × 2 × 29.34 × 10−31.15 × 10−22.43 × 10−210.17 (1905.17)
    100 × 50 × 2 × 48.90 × 10−31.07 × 10−22.53 × 10−219.54 (3660.69)
    100 × 50 × 3 × 29.28 × 10−31.13 × 10−22.35 × 10−213.63 (2544.14)
    25 × 25800 × 400 × 1 × 1631.39 (118285.12)
    DownLoad: CSV

    Table 2.  Normalized errors of the potential temperature perturbation ($\theta '$) and elapsed CPU time for the density current test running on different grids.

    Resolution (m)Grid configuration${l_1}$${l_2}$$ {l_\infty } $CPU time (s)
    200 × 200132 × 32 × 1 × 10.280.330.941 (235.36)
    100 × 100265 × 64 × 1 × 10.130.170.608.69 (2045.44)
    132 × 32 × 2 × 20.130.170.643.44 (809.19)
    50 × 50530 × 128 × 1 × 12.63 × 10−23.59 × 10−20.1569.93 (16458.58)
    265 × 64 × 2 × 22.72 × 10−23.67 × 10−20.1711.10 (2613.85)
    132 × 32 × 2 × 42.64 × 10−23.59 × 10−20.1525.99 (6117.59)
    132 × 32 × 3 × 22.68 × 10−23.66 × 10−20.1820.00 (4047.29)
    25 × 251060 × 256 × 1 × 1640.80 (150820.42)
    DownLoad: CSV

    Table 3.  Normalized errors of potential temperature perturbation ($\theta '$) and elapsed CPU time for internal gravity waves test running on different grids.

    Resolution (m)Grid configuration${l_1}$${l_2}$$ {l_\infty } $CPU time (s)
    4000 × 40075 × 25 × 1 × 10.250.290.321 (48.60)
    2000 × 200150 × 50 × 1 × 18.20 × 10−20.110.159.77 (474.99)
    75 × 25 × 2 × 27.95 × 10−20.110.164.13 (200.48)
    1000 × 100300 × 100 × 1 × 11.48 × 10−22.12 × 10−23.58 × 10−278.87 (3833.06)
    150 × 50 × 2 × 22.77 × 10−23.60 × 10−26.09 × 10−231.37 (1524.60)
    75 × 25 × 2 × 41.62 × 10−22.14 × 10−23.57 × 10−227.58 (1340.39)
    75 × 25 × 3 × 22.84 × 10−23.61 × 10−26.08 × 10−226.33 (1279.57)
    250 × 251200 × 400 × 1 × 16375.50 (309849.29)
    DownLoad: CSV

    Table 4.  Maximum and minimum of vertical velocity, $w$, and potential temperature perturbation, $\theta '$ , for the internal gravity waves test on the AMR grids with the finest resolution, $\Delta z = 100{\text{ m}}$ after 3000 s.

    Gridwmax (m s−1)wmin (m s−1)$\theta {'_{\max }}$ (K)$ \theta {'_{\min }} $ (K)
    75 × 25 × 2 × 42.47 × 10−3−2.42 × 10−32.80 × 10−3−1.52 × 10−3
    75 × 25 × 3 × 22.47 × 10−3−2.42 × 10−32.80 × 10−3−1.52 × 10−3
    150 × 50 × 2 × 22.47 × 10−3−2.43 × 10−32.80 × 10−3−1.52 × 10−3
    DownLoad: CSV

    Table 5.  Normalized errors of horizontal velocity ($u$) and elapsed CPU time for the Schär mountain test running on different grids.

    Resolution (m)Grid configuration${l_1}$${l_2}$$ {l_\infty } $CPU time (s)
    500 × 420100 × 50 × 1 × 12.97 × 10−35.66 × 10−32.19 × 10−21 (4028.53)
    250 × 210200 × 100 × 1 × 14.31 × 10−48.66 × 10−49.41 × 10−38.66 (34899.12)
    100 × 50 × 2 × 28.61 × 10−41.28 × 10−31.03 × 10−23.83 (15429.89)
    125 × 105400 × 200 × 1 × 183.10 (334779.90)
    DownLoad: CSV
  • Ahmad, N., and J. Lindeman, 2007: Euler solutions using flux-based wave decomposition. International Journal for Numerical Methods in Fluids, 54, 47−72, https://doi.org/10.1002/fld.1392.
    Bacon, D. P., and Coauthors, 2000: A dynamically adapting weather and dispersion model: The operational multiscale environment model with grid adaptivity (OMEGA). Mon. Wea. Rev., 128, 2044−2076, https://doi.org/10.1175/1520-0493(2000)128<2044:ADAWAD>2.0.CO;2.
    Behrens, J., 1996: An adaptive semi-lagrangian advection scheme and its parallelization. Mon. Wea. Rev., 124, 2386−2395, https://doi.org/10.1175/1520-0493(1996)124<2386:AASLAS>2.0.CO;2.
    Berger, M., and I. Rigoutsos, 1991: An algorithm for point clustering and grid generation. IEEE Transactions on Systems, Man, and Cybernetics, 21, 1278−1286, https://doi.org/10.1109/21.120081.
    Berger, M. J., and J. Oliger, 1984: Adaptive mesh refinement for hyperbolic partial differential equations. J. Comput. Phys., 53, 484−512, https://doi.org/10.1016/0021-9991(84)90073-1.
    Berger, M. J., and P. Colella, 1989: Local adaptive mesh refinement for shock hydrodynamics. J. Comput. Phys., 82, 64−84, https://doi.org/10.1016/0021-9991(89)90035-1.
    Berger, M. J., and R. J. LeVeque, 1998: Adaptive mesh refinement using wave-propagation algorithms for hyperbolic systems. SIAM Journal on Numerical Analysis, 35, 2298−2316, https://doi.org/10.1137/S0036142997315974.
    Blaise, S., and A. St-Cyr, 2012: A dynamic hp-adaptive discontinuous Galerkin method for shallow-water flows on the sphere with application to a global tsunami simulation. Mon. Wea. Rev., 140, 978−996, https://doi.org/10.1175/MWR-D-11-00038.1.
    Blayo, E., and L. Debreu, 1999: Adaptive mesh refinement for finite-difference ocean models: First experiments. J. Phys. Oceanogr., 29, 1239−1250, https://doi.org/10.1175/1520-0485(1999)029<1239:AMRFFD>2.0.CO;2.
    Carpenter, R. L. Jr., K. K. Droegemeier, P. R. Woodward, and C. E. Hane, 1990: Application of the piecewise parabolic method (PPM) to meteorological modeling. Mon. Wea. Rev., 118, 586−612, https://doi.org/10.1175/1520-0493(1990)118<0586:AOTPPM>2.0.CO;2.
    Chen, C. G., and F. Xiao, 2008: Shallow water model on cubed-sphere by multi-moment finite volume method. J. Comput. Phys., 227, 5019−5044, https://doi.org/10.1016/j.jcp.2008.01.033.
    Chen, C. G., F. Xiao, and X. L. Li, 2011: An adaptive multimoment global model on a cubed sphere. Mon. Wea. Rev., 139, 523−548, https://doi.org/10.1175/2010MWR3365.1.
    Giraldo, F. X., and T. E. Rosmond, 2004: A scalable spectral element Eulerian atmospheric model (SEE-AM) for NWP: Dynamical core tests. Mon. Wea. Rev., 132, 133−153, https://doi.org/10.1175/1520-0493(2004)132<0133:ASSEEA>2.0.CO;2.
    Giraldo, F. X., and M. Restelli, 2008: A study of spectral element and discontinuous Galerkin methods for the Navier-Stokes equations in nonhydrostatic mesoscale atmospheric modeling: Equation sets and test cases. J. Comput. Phys., 227, 3849−3877, https://doi.org/10.1016/j.jcp.2007.12.009.
    Hubbard, M. E., and N. Nikiforakis, 2003: A three-dimensional, adaptive, Godunov-type model for global atmospheric flows. Mon. Wea. Rev., 131, 1848−1864, https://doi.org/10.1175//2568.1.
    Ii, S., and F. Xiao, 2007: CIP/multi-moment finite volume method for Euler equations: A semi-Lagrangian characteristic formulation. J. Comput. Phys., 222, 849−871, https://doi.org/10.1016/j.jcp.2006.08.015.
    Ii, S., and F. Xiao, 2009: High order multi-moment constrained finite volume method. Part I: Basic formulation. J. Comput. Phys., 228, 3669−3707, https://doi.org/10.1016/j.jcp.2009.02.009.
    Ii, S., M. Shimuta, and F. Xiao, 2005: A 4th-order and single-cell-based advection scheme on unstructured grids using multi-moments. Computer Physics Communications, 173, 17−33, https://doi.org/10.1016/j.cpc.2005.07.003.
    Iskandarani, M., D. B. Haidvogel, J. C. Levin, E. Curchitser, and C. A. Edwards, 2002: Multi-scale geophysical modeling using the spectral element method. Computing in Science & Engineering, 4, 42−48, https://doi.org/10.1109/MCISE.2002.1032428.
    Jablonowski, C., M. Herzog, J. E. Penner, R. C. Oehmke, Q. F. Stout, B. van Leer, and K. G. Powell, 2006: Block-structured adaptive grids on the sphere: Advection experiments. Mon. Wea. Rev., 134, 3691−3713, https://doi.org/10.1175/MWR3223.1.
    Kessler, M., 1999: Development and analysis of an adaptive transport scheme. Atmos. Environ., 33, 2347−2360, https://doi.org/10.1016/S1352-2310(98)00415-4.
    Levy, M. N., R. D. Nair, and H. M. Tufo, 2007: High-order Galerkin methods for scalable global atmospheric models. Computers & Geosciences, 33, 1022−1035, https://doi.org/10.1016/j.cageo.2006.12.004.
    Li, X. L., C. G., Chen, X. S., Shen and F., Xiao, 2013: A Multimoment constrained finite-volume model for nonhydrostatic atmospheric dynamics. Mon. Wea. Rev, 141, 1216−1240, https://doi.org/10.1175/MWR-D-12-00144.1.
    Nair, R. D., H. W. Choi, and H. M. Tufo, 2009: Computational aspects of a scalable high-order discontinuous Galerkin atmospheric dynamical core. Computers & Fluids, 38, 309−319, https://doi.org/10.1016/j.compfluid.2008.04.006.
    Norman, M. R., R. D. Nair, and F. H. M. Semazzi, 2011: A low communication and large time step explicit finite-volume solver for non-hydrostatic atmospheric dynamics. J. Comput. Phys., 230, 1567−1584, https://doi.org/10.1016/j.jcp.2010.11.022.
    Pielke, R. A., and Coauthors, 1992: A comprehensive meteorological modeling system-RAMS. Meteorol. Atmos. Phys., 49, 69−91, https://doi.org/10.1007/BF01025401.
    Qin, Q. C., X. S. Shen, C. G. Chen, F. Xiao, Y. J. Dai, and X. L. Li, 2019: A 3D nonhydrostatic compressible atmospheric dynamic core by multi-moment constrained finite volume method. Adv. Atmos. Sci., 36, 1129−1142, https://doi.org/10.1007/s00376-019-9002-4.
    Schär, C., D. Leuenberger, O. Fuhrer, D. Lüthi, and C. Girard, 2002: A new terrain-following vertical coordinate formulation for atmospheric prediction models. Mon. Wea. Rev., 130, 2459−2480, https://doi.org/10.1175/1520-0493(2002)130<2459:ANTFVC>2.0.CO;2.
    Shu, C. W., 1988: Total-variation-diminishing time discretizations. SIAM Journal on Scientific and Statistical Computing, 9, 1073−1084, https://doi.org/10.1137/0909073.
    Skamarock, W. C., and J. B. Klemp, 1993: Adaptive grid refinement for two-dimensional and three-dimensional nonhydrostatic atmospheric flow. Mon. Wea. Rev., 121, 788−804, https://doi.org/10.1175/1520-0493(1993)121<0788:AGRFTD>2.0.CO;2.
    Skamarock, W. C., and J. B. Klemp, 1994: Efficiency and accuracy of the Klemp-Wilhelmson time-splitting technique. Mon. Wea. Rev., 122, 2623−2630, https://doi.org/10.1175/1520-0493(1994)122<2623:EAAOTK>2.0.CO;2.
    Skamarock, W., J. Oliger, and R. L. Street, 1989: Adaptive grid refinement for numerical weather prediction. J. Comput. Phys., 80, 27−60, https://doi.org/10.1016/0021-9991(89)90089-2.
    St-Cyr, A., C. Jablonowski, J. M. Dennis, H. M. Tufo, and S. J. Thomas, 2008: A comparison of two shallow-water models with nonconforming adaptive grids. Mon. Wea. Rev., 136, 1898−1922, https://doi.org/10.1175/2007MWR2108.1.
    Stevens, D. E., and S. Bretherton, 1996: A forward-in-time advection scheme and adaptive multilevel flow solver for nearly incompressible atmospheric flow. J. Comput. Phys., 129, 284−295, https://doi.org/10.1006/jcph.1996.0250.
    Straka, J. M., R. B. Wilhelmson, L. J. Wicker, J. R. Anderson, and K. K. Droegemeier, 1993: Numerical solutions of a non-linear density current: A benchmark solution and comparisons. International Journal for Numerical Methods in Fluids, 17, 1−22, https://doi.org/10.1002/fld.1650170103.
    Thomas, S. J., and R. D. Loft, 2000: Parallel semi-implicit spectral element methods for atmospheric general circulation models. J. Comput. Phys., 15, 499−518, https://doi.org/10.1023/A:1011188832645.
    Tomlin, A., M. Berzins, J. Ware, J. Smith, and M. J. Pilling, 1997: On the use of adaptive gridding methods for modelling chemical transport from multi-scale sources. Atmos. Environ., 31, 2945−2959, https://doi.org/10.1016/S1352-2310(97)00120-9.
    Wicker, L. J., and W. C. Skamarock, 1998: A time-splitting scheme for the elastic equations incorporating second-order Runge-Kutta time differencing. Mon. Wea. Rev., 126, 1992−1999, https://doi.org/10.1175/1520-0493(1998)126<1992:ATSSFT>2.0.CO;2.
    Williamson, D. L., J. B. Drake, J. J. Hack, R. Jakob, and P. N. Swarztrauber, 1992: A standard test set for numerical approximations to the shallow water equations in spherical geometry. J. Comput. Phys., 102, 211−224, https://doi.org/10.1016/S0021-9991(05)80016-6.
    Xiao, F., R. Akoh, and S. Ii, 2006: Unified formulation for compressible and incompressible flows by using multi-integrated moments II: Multi-dimensional version for compressible and incompressible flows. J. Comput. Phys., 213, 31−56, https://doi.org/10.1016/j.jcp.2005.08.002.
    Yabe, T., and T. Aoki, 1991: A universal solver for hyperbolic equations by cubic-polynomial interpolation I. One-dimensional solver. Computer Physics Communications, 66, 219−232, https://doi.org/10.1016/0010-4655(91)90071-R.
    Yabe, T., F. Xiao, and T. Utsumi, 2001: The constrained interpolation profile method for multiphase analysis. J. Comput. Phys., 169, 556−593, https://doi.org/10.1006/jcph.2000.6625.
    Yessad, K., and P. Bénard, 1996: Introduction of a local mapping factor in the spectral part of the Météo-France global variable mesh numerical forecast model. Quart. J. Roy. Meteor. Soc., 122, 1701−1719, https://doi.org/10.1002/QJ.49712253511.
  • [1] Qingchang QIN, Xueshun SHEN, Chungang CHEN, Feng XIAO, Yongjiu DAI, Xingliang LI, 2019: A 3D Nonhydrostatic Compressible Atmospheric Dynamic Core by Multi-moment Constrained Finite Volume Method, ADVANCES IN ATMOSPHERIC SCIENCES, 36, 1129-1142.  doi: 10.1007/s00376-019-9002-4
    [2] Yang XIA, Bin WANG, Lijuan LI, Li LIU, Jianghao LI, Li DONG, Shiming XU, Yiyuan LI, Wenwen XIA, Wenyu HUANG, Juanjuan LIU, Yong WANG, Hongbo LIU, Ye PU, Yujun HE, Kun XIA, 2024: A Neural-network-based Alternative Scheme to Include Nonhydrostatic Processes in an Atmospheric Dynamical Core, ADVANCES IN ATMOSPHERIC SCIENCES, 41, 1083-1099.  doi: 10.1007/s00376-023-3119-1
    [3] Sun Litan, Huang Meiyuan, 1994: Improving the Vorticity-Streamfunction Method to Solve Two-Dimensional Anelastic and Nonhydrostatic Model, ADVANCES IN ATMOSPHERIC SCIENCES, 11, 247-249.  doi: 10.1007/BF02666551
    [4] Daosheng XU, Dehui CHEN, Kaixin WU, 2021: Properties of High-Order Finite Difference Schemes and Idealized Numerical Testing, ADVANCES IN ATMOSPHERIC SCIENCES, 38, 615-626.  doi: 10.1007/s00376-020-0130-7
    [5] SHEN Xueshun, Akimasa SUMI, 2005: A High Resolution Nonhydrostatic Tropical Atmospheric Model and Its Performance, ADVANCES IN ATMOSPHERIC SCIENCES, 22, 30-38.  doi: 10.1007/BF02930867
    [6] LI Xingliang, SHEN Xueshun, PENG Xindong, XIAO Feng, ZHUANG Zhaorong, CHEN Chungang, 2013: An Accurate Multimoment Constrained Finite Volume Transport Model on Yin-Yang Grids, ADVANCES IN ATMOSPHERIC SCIENCES, 30, 1320-1330.  doi: 10.1007/s00376-013-2217-x
    [7] Ji Zhongzhen, Li Jing, 1998: A High-Order Compact Scheme with Square-Conservativity, ADVANCES IN ATMOSPHERIC SCIENCES, 15, 580-584.  doi: 10.1007/s00376-998-0034-4
    [8] Jie TANG, Chungang CHEN, Xueshun SHEN, Feng XIAO, Xingliang LI, 2021: A Positivity-preserving Conservative Semi-Lagrangian Multi-moment Global Transport Model on the Cubed Sphere, ADVANCES IN ATMOSPHERIC SCIENCES, 38, 1460-1473.  doi: 10.1007/s00376-021-0393-7
    [9] ZHAO Bin, ZHONG Qing, 2010: The Development of a Nonhydrostatic Global Spectral Model, ADVANCES IN ATMOSPHERIC SCIENCES, 27, 676-684.  doi: 10.1007/s00376-009-9080-9
    [10] Shaoying LI, Jun PENG, Weimin ZHANG, Jianping WU, Qiang YAO, Xiangrong YANG, Tengling LUO, 2023: Effects of a Dry-Mass Conserving Dynamical Core on the Simulation of Tropical Cyclones, ADVANCES IN ATMOSPHERIC SCIENCES, 40, 464-482.  doi: 10.1007/s00376-022-2085-3
    [11] Rui LYU, Fei HU, Lei LIU, Jingjing XU, Xueling CHENG, 2018: High-Order Statistics of Temperature Fluctuations in an Unstable Atmospheric Surface Layer over Grassland, ADVANCES IN ATMOSPHERIC SCIENCES, 35, 1265-1276.  doi: 10.1007/s00376-018-7248-x
    [12] Pengfei WANG, 2017: A High-Order Spatiotemporal Precision-Matching Taylor-Li Scheme for Time-Dependent Problems, ADVANCES IN ATMOSPHERIC SCIENCES, 34, 1461-1471.  doi: 10.1007/s00376-017-7018-1
    [13] Wang Weiguo, Xie Yingqi, Qiu Jinhuan, Liu Qing, 1998: The Regional Dynamical Model of the Atmospheric Ozonosphere, ADVANCES IN ATMOSPHERIC SCIENCES, 15, 74-82.  doi: 10.1007/s00376-998-0019-3
    [14] Fei Jianfang, Zhang Ming, Lu Hancheng, Yang Guoxiang, 1992: A Fine-Mesh Numerical Model with Detailed Boundary Layer Parameterization, ADVANCES IN ATMOSPHERIC SCIENCES, 9, 465-476.  doi: 10.1007/BF02677079
    [15] Wang Guomin, Wang Shiwen, Li Jianjun, 1996: A Bogus Typhoon Scheme and Its Application to a Movable Nested Mesh Model, ADVANCES IN ATMOSPHERIC SCIENCES, 13, 103-114.  doi: 10.1007/BF02657031
    [16] Ye Weizuo, 1990: A Coupled Dynamical-Radiational Model of Stratocumulus, ADVANCES IN ATMOSPHERIC SCIENCES, 7, 197-210.  doi: 10.1007/BF02919158
    [17] Ye Weizuo, 1990: Application of the Coupled Dynamical-Radiational Model of Stratocumulus, ADVANCES IN ATMOSPHERIC SCIENCES, 7, 331-346.  doi: 10.1007/BF03179765
    [18] XIANG Jie, LIAO Qianfeng, HUANG Sixun, LAN Weiren, FENG Qiang, ZHOU Fengcai, 2006: An Application of the Adjoint Method to a Statistical-Dynamical Tropical-Cyclone Prediction Model (SD–90) II: Real Tropical Cyclone Cases, ADVANCES IN ATMOSPHERIC SCIENCES, 23, 118-126.  doi: 10.1007/s00376-006-0012-7
    [19] Liu Yudi, Ji Zhongzhen, Wang Bin, 2002: Study on Computational Properties of Several Vertical Grids with a Nonhydrostatic Model in Comparison to Analytical Solutions, ADVANCES IN ATMOSPHERIC SCIENCES, 19, 528-543.  doi: 10.1007/s00376-002-0084-y
    [20] Yu ZHANG, Yuanfu XIE, Hongli WANG, Dehui CHEN, Zoltan TOTH, 2016: Ensemble Transform Sensitivity Method for Adaptive Observations, ADVANCES IN ATMOSPHERIC SCIENCES, 33, 10-20.  doi: 10.1007/s00376-015-5031-9

Get Citation+

Export:  

Share Article

Manuscript History

Manuscript received: 21 May 2021
Manuscript revised: 14 August 2021
Manuscript accepted: 26 August 2021
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

An Adaptive Nonhydrostatic Atmospheric Dynamical Core Using a Multi-Moment Constrained Finite Volume Method

    Corresponding author: Chungang CHEN, cgchen@xjtu.edu.cn
  • 1. State Key Laboratory for Strength and Vibration of Mechanical Structures, Xi’an Jiaotong University, Xi’an 710049, China
  • 2. Center of Numerical Weather Prediction, China Meteorological Administration, Beijing 100081, China
  • 3. Department of Mechanical Engineering, Tokyo Institute of Technology, Tokyo 226-8502, Japan

Abstract: An adaptive 2D nonhydrostatic dynamical core is proposed by using the multi-moment constrained finite-volume (MCV) scheme and the Berger-Oliger adaptive mesh refinement (AMR) algorithm. The MCV scheme takes several point-wise values within each computational cell as the predicted variables to build high-order schemes based on single-cell reconstruction. Two types of moments, such as the volume-integrated average (VIA) and point value (PV), are defined as constraint conditions to derive the updating formulations of the unknowns, and the constraint condition on VIA guarantees the rigorous conservation of the proposed model. In this study, the MCV scheme is implemented on a height-based, terrain-following grid with variable resolution to solve the nonhydrostatic governing equations of atmospheric dynamics. The AMR grid of Berger-Oliger consists of several groups of blocks with different resolutions, where the MCV model developed on a fixed structured mesh can be used directly. Numerical formulations are designed to implement the coarse-fine interpolation and the flux correction for properly exchanging the solution information among different blocks. Widely used benchmark tests are carried out to evaluate the proposed model. The numerical experiments on uniform and AMR grids indicate that the adaptive model has promising potential for improving computational efficiency without losing accuracy.

摘要: 本文提出了一种基于多矩限制有限体积法和Berger-Oliger自适应加密网格的二维非静力大气动力框架。多矩限制有限体积法算法在每个计算单元内选择多个点值作为预报量,利用独立单元上的局地高阶重构构造高精度的数值格式。多矩限制有限体积法算法定义了两种类型的矩,即积分平均值矩和点值矩,它们作为限制条件约束预报量的更新,其中对积分平均值矩的限制条件保证了数值格式的严格守恒性。本文在基于高度和地形跟踪的变分辨率网格上采用多矩限制有限体积法算法求解非静力大气运动方程,其中所使用的Berger-Oliger自适应加密网格是由几组具有不同分辨率的均匀网格块构成的,在固定的结构网格上发展的多矩限制有限体积法算法可以直接在这些网格块上使用。我们还设计了粗细网格插值和通量修正的数值方法,使不同层网格块之间的信息得以双向交换。一些广泛使用的标准数值算例被用于检验本文所提出的模式,数值实验的结果表明,在保证数值精度的前提下,自适应模式对于提升计算效率具有可观的前景。

    • During the last few decades, increased attention has been paid toward improving numerical weather prediction (NWP) models for better resolving small-scale processes at high resolution. In terms of a dynamical core, much effort has been put forth to develop new numerical schemes and computational meshes with the intent of attaining higher accuracy and efficiency. Aside from the traditional methods, the high-order numerical schemes that are based on “local” or single-cell based spatial reconstruction are becoming more popular recently in developing the atmospheric dynamical cores because of their high-order accuracy, geometric flexibility, and their scalability on massive, parallel computer systems. The representative schemes are the spectral element method (Thomas and Loft, 2000; Iskandarani et al., 2002; Giraldo and Rosmond, 2004) and the discontinuous Galerkin method (Levy et al., 2007; Giraldo and Restelli, 2008; Nair et al., 2009; Blaise and St-Cyr, 2012) among others.

      The multi-moment constrained finite-volume (MCV) method (Ii and Xiao, 2009) was recently proposed as a new framework to develop high-order schemes. Several point-wise values at equidistantly distributed solution points are defined as predicted variables (unknowns) within each computational cell. The formulations for the spatial discretization of unknowns are derived through the constraint conditions provided by introducing different types of moments (Yabe and Aoki, 1991; Yabe et al., 2001; Ii et al., 2005; Xiao et al., 2006; Ii and Xiao, 2007; Chen and Xiao, 2008), such as point values, volume-integrated averages, derivatives of different orders, and so on. Using the MCV method, 2D and 3D nonhydrostatic dynamical cores have been proposed in Li et al. (2013) and Qin et al. (2019).

      Though computational power has made rapid gains recently, atmospheric modeling is still highly time-consuming and computationally expensive work. Considering the multi-scale phenomena at operating in the atmosphere, the development of adaptive numerical models, which are capable of running on the computational meshes with variable spatial resolution, is a very effective way to save on computational costs. Some existing methods can conduct multi-resolution simulations in atmospheric modeling, such as the nested grid (Pielke et al., 1992), the stretched grid (Yessad and Bénard, 1996), and the AMR grid (Berger and Oliger, 1984; Behrens, 1996; Stevens and Bretherton, 1996; Tomlin et al., 1997; Kessler, 1999). Using an AMR grid is very attractive among these techniques for its flexibility and efficiency in dynamic adjusting a computational mesh according to the time-evolution of predicted variables. AMR techniques have been applied in some atmospheric models to improve the computational efficiency, examples include the studies described in Skamarock et al. (1989), Skamarock and Klemp (1993), Blayo and Debreu (1999), Bacon et al. (2000), Hubbard and Nikiforakis (2003), Jablonowski et al. (2006), St-Cyr et al. (2008), and Chen et al. (2011). In this study, we will develop an adaptive MCV nonhydrostatic model with the applications of Berger and Oliger’s AMR algorithm based on the 2D model proposed in Li et al. (2013).

      The remainder of this paper is organized as follows. A 2D nonhydrostatic dynamical core using a 3-point MCV scheme is briefly described in section 2. The setup and dynamic adjustment of an AMR mesh, along with numerical operations to exchange the solution information among different blocks under the framework of the MCV scheme, are described in section 3. Numerical results of the widely used benchmark tests are given in section 4, and a short summary is given in section 5.

    2.   The two-dimensional (2D) compressible nonhydrostatic MCV dynamical core
    • A 2D compressible nonhydrostatic atmospheric dynamical core was proposed in Li et al. (2013) using a multi-moment constrained finite volume method. We implement this MCV dynamical core on the AMR grid in this study. The governing equations for the compressible nonhydrostatic dynamical core are written in conservative flux-form as

      where the dependent variables are the density ($\rho $), the momentum vector $\left( {\rho {\boldsymbol{u}},\rho {\boldsymbol{w}}} \right)$, and the product of density and potential temperature $\rho \theta $; p is pressure and g is the acceleration of gravity. The effects of rotation are not considered in this implementation.

      To deal with bottom topography, the height-based, terrain-following coordinates $\left( {x,\zeta } \right)$ proposed in Schär et al. (2002) are adopted. The Jacobians of the transformation are $\sqrt G = \dfrac{{\partial z}}{{\partial \zeta }}$ and ${G_{13}} = \dfrac{{\partial \zeta }}{{\partial x}}$. In computational space, the uniform grid $\left( {x,\zeta } \right)$ is related to the Cartesian coordinate as

      where ${h_s}$ is the height of the bottom of the mountain, H is the model top, and the parameter S = 3000 m is adopted. $\widetilde w$ is the contravariant vertical velocity component in the transformed coordinate defined as

      where u and w are velocity components for the Cartesian coordinate system $\left( {x,z} \right)$

      Splitting of the hydrostatic reference state is usually applied in atmospheric modeling to improve accuracy. As a result, the thermodynamical variables are written as

    • A 3-point MCV scheme of third-order accuracy (Ii and Xiao, 2009) is adopted to solve the governing equations, Eqs. (1)–(4) in this study. Hereafter, we briefly describe the MCV formulations by considering the 1D hyperbolic conservation law having the form of

      where ${\boldsymbol{q}}$ is the vector of dependent variables and ${\boldsymbol{f}}$ is the vector of flux function.

      Within each cell, three local degrees of freedom (DOFs), in terms of point-wise values ${{\boldsymbol{q}}_{i,m}}$ (m = 1 to 3), are defined at two endpoints ${x_{i \pm \frac{1}{2}}}$ and the center ${x_i}$ as shown in Fig. 1a for the cell ${C_i}$.

      Figure 1.  Distribution of DOFs for the (a) 1D case and (b) 2D Cartesian grid.

      With known values of local DOFs, a quadratic Lagrangian polynomial ${{\boldsymbol{Q}}_i}\left( x \right)$ can be constructed as

      where the Lagrangian basis function is

      To build the MCV scheme, two kinds of moments including the PV (point value) moment ${\overline {{}^P{\boldsymbol{q}}} _{i - \frac{1}{2}}}$ at cell interface and the VIA (volume-integrated average) moment ${\overline {{}^V{\boldsymbol{q}}} _i}$ over each cell are defined to derive the spatial discretization formulations (Ii and Xiao, 2009). At the cell interface, the PV moment (same as the DOF defined at the cell’s endpoint) is updated using the differential-form governing equations through solving the derivative Riemann problem. The VIA moment is updated using flux-form finite volume formulations and assures the rigorous numerical conservation. The updating formulations for the DOF defined at the cell center are then derived through a constraint condition on the VIA moment. The details of the numerical procedure to derive the updating formulations for the 3-point MCV scheme can be referred to Ii and Xiao (2009). For the sake of brevity, we directly give the semi-discrete formulations for three different DOFs as

      The numerical fluxes ${\widehat {\boldsymbol{f}}_{i \pm \frac{1}{2}}}$ at cell interfaces are C0-continuous (i.e., the zeroth derivative of numerical fluxes at cell interfaces are continuous) and can be directly computed from the values of known DOFs. The first-order derivatives of flux functions ${\widehat {\boldsymbol{f}}_{x,i \pm \frac{1}{2}}}$ are generally discontinuous and determined by solving derivative Riemann problems. In this study, a local Lax-Friedrich (LLF) solver (Shu, 1988) is adopted and the derivatives of flux functions at the cell interface $x = {x_{i + \frac{1}{2}}}$ are evaluated by

      where a is the maximum absolute value of the eigenvalues of the Jacobian matrix ${{\partial {\boldsymbol{f}}}}/{{\partial {\boldsymbol{q}}}}$. The derivatives of dependent variables and flux functions are evaluated using the Lagrangian interpolation polynomials having the form of Eq. (10), and the subscripts L and R denote the cells ${C_i}$ and ${C_{i + 1}}$.

      As shown in Ii and Xiao (2009), the 3-point MCV scheme shown above is of 3rd-order accuracy and preserves the rigorous numerical conservation in terms of the VIA moment given by

      The MCV scheme can be easily extended to the multi-dimensional models by applying the above one-dimensional formulations in different directions one-by-one (Ii and Xiao, 2009). As shown in Fig. 1b, nine local DOFs in terms of point-wise values are defined within the cell ${C_{i,j}}$ to build the two-dimensional numerical models using the 3-point MCV scheme. Similarly, the numerical conservation is assured by the VIA moment ${\overline {{}^V{\boldsymbol{q}}} _{i,j}}$, which is determined by the two-dimensional Simpson’s rule as

      where the indices i and j are omitted on the right-hand side of Eq. (17) for the sake of brevity.

      For time-stepping, a third-order Runge-Kutta scheme is adopted to update the following semi-discrete ordinary differential equation (ODE) obtained from the above spatial discretization:

      The solution is advanced from ${{\boldsymbol{q}}_{(n)}}$ at time ${t_{(n)}}$ to ${{\boldsymbol{q}}_{(n + 1)}}$ at time ${t_{(n + 1)}}$ by

      with $\Delta t = {t_{(n + 1)}} - {t_{(n)}}$ as the time step. The time integration strategy for the AMR grids will be demonstrated in later contents.

      The details of the numerical formulations for a two-dimensional MCV nonhydrostatic dynamical core are described in Li et al. (2013).

    3.   Adaptive mesh refinement algorithm
    • The AMR algorithm proposed by Berger and Oliger (1984) uses a set of blocks to cover the whole computational domain. On each block, the structured computational mesh is constructed to implement the spatial discretization. These blocks are grouped into different levels according to their spatial resolutions and the blocks with finer resolution always overlay the blocks belonging to the next coarser levels. The governing equations can be solved on each block independently if enough ghost cells, which are either copied from the adjacent blocks of the same level or evaluated by coarse-fine interpolations from blocks of the next coarser level, are supplemented. Additionally, the AMR algorithm also implements the dynamic adjustment of blocks according to the time evolution of the predicted fields. Since the high-order discretization is implemented using the single-cell based stencils in MCV models, the coarse-fine interpolations are easily accomplished and extra numerical errors due to the non-uniform resolution can be more effectively suppressed in comparison with existing variable-resolution models, as shown in our previous study (Chen et al., 2011). The operations to generate and adjust the AMR grid along with some additional numerical procedures to implement the MCV model on an AMR grid are briefly introduced in this section. Without losing generality, here, we consider the uniform Cartesian grid $\left( {x,y} \right)$ on each block, i.e., the grid spacing is ${h_k} = \Delta {x_k} = \Delta {y_k}$ on level k.

    • In the beginning, a base block (level 0) is first built, which covers the whole computational domain and has the coarsest resolution. Then, finer levels are generated one by one following a similar procedure. Given the kth-level blocks, the next finer level is generated as follows.

    • All cells belonging to level k are checked by some refinement criterion to find those cells where the physical fields should be solved using finer resolution to improve the accuracy. Berger and Oliver (1984) introduced a criterion based on the error estimation using Richardson extrapolation. Though it is a general methodology for different types of applications, a more efficient criterion is to check the known flow structures of some representative physical variable in atmospheric modeling. In this study, we flag the cell to be refined if the variation of some indicative variable across the cell is larger than a prescribed threshold.

    • The simplest way is to use just one block, which is the smallest rectangle (in 2D computational space) is to cover all flagged cells. However, it is inefficient since, in general, the domain consisting of the flagged cells is highly irregular, and using a single structured block will result in many un-flagged cells being involved at the finer level $k{\text{ + 1}}{\text{.}}$ To improve the efficiency, we divide one big block into two small ones and it is expected that the area of each small one then can be reduced to cover the flagged cells. The division of an existing block will be conducted continuously until the ratio between the area covered by un-flagged cells and the total area on this block is smaller than a prescribed value. In this study, we adopted a division method based on pattern recognition (Berger and Rigoutsos, 1991) instead of a simple bisection to speed up this operation.

      The above procedure is repeated until the finest level has been constructed. During the simulations, the computational meshes are completely or partially adjusted every several time steps to track the evolution of predicted fields. This adjustment starts from the fine level and most of the above operations can be applied with some modifications. When the grid adjustment is finished, values of predicted variables in the newly generated fine cells are interpolated from the coarser level. Coarse-fine interpolation operations based on a single-cell stencil are described later in this section.

    • In this study, a time marching step is chosen to keep the same CFL (Courant-Friedrichs-Lewy) number on different levels. In two-dimensional case, the CFL number on level k is $\frac{{u\Delta {t_k}}}{{\Delta {x_k}}} + \frac{{v\Delta {t_k}}}{{\Delta {y_k}}}$. Considering that the refinement ratio of grid resolution is an integer, $\lambda $, between levels $k$and $k + 1$, the corresponding grid spacings and time steps on these two adjacent levels are ${h_k} = \lambda {h_{k + 1}}$ and $\Delta {t_k} = \lambda \Delta {t_{k + 1}}$. With variable time steps on different levels, the governing equations are advanced for $\lambda $ steps on level $k + 1$ during one step on level $k$. In Fig. 2, the time marching procedure from ${t_n}$ to ${t_{n + 1}} = {t_n} + \Delta t$ is illustrated for a 3-level AMR grid. Each level is advanced through Eq. (19) with its own time integration interval. For the sake of simplicity, the refinement ratio is chosen as $\lambda = 2$ hereafter in all illustrative examples. A special recursive procedure is adopted for time marching on an AMR grid, i.e., any coarse level cannot be advanced to the next time step until all finer levels have reached the current time instant. Additionally, the coarse-fine interpolations are conducted both in time and space to evaluate the ghost cells. In this study, the temporal interpolation is accomplished by a linear reconstruction of the predicted fields at two adjacent time steps.

      Figure 2.  Diagram depicting the subcycling of AMR levels in time.

      To guarantee the numerical conservation, flux correction should be conducted along the interfaces between the coarse and fine levels (Berger, 1989, 1998). As shown in Fig. 3, $F_{i + 1,j,k}^{}$ represents the time-averaged numerical flux during one time step on level $ k $, and is used to update $ \overline {{}^Vq} _{i,j,k}^{} $, the VIA moment of the coarse cell $C_{i,j,k}^{}$. During the same time interval, VIA moments of the adjacent fine cells $ C_{2i + 1,2j,k + 1}^{} $ and $C_{2i + 1,2j - 1,k + 1}^{}$ are advanced by time-averaged numerical fluxes $f_{2i + 1,2j,k + 1}^{}$ and $f_{2i + 1,2j - 1,k + 1}^{}$. Since the relation $F_{i + 1,j,k}^{} = {(f_{2i + 1,2j - 1,k + 1}^{} + f_{2i + 1,2j,k + 1})}/{2}$ is not assured by the raw scheme, the AMR model is not conservative without a flux correction operation.

      Figure 3.  Flux correction along the interface between the coarse and fine levels.

      In the proposed model, the VIA moment of a coarse cell is corrected to assure the numerical conservation by:

      where $A_{i,j,k}^{}$ is the area of the coarse cell, $C_{i,j,k}^{}$.

    • The solution transfer between coarse and fine levels plays an important role to accomplish the numerical modeling on an AMR grid. In the proposed model, numerical solutions on level $k$ are only exchanged with the two adjacent levels, $k - 1$ and $k + 1$.

      Since the numerical solutions with finer resolution are of higher accuracy, they are used to replace the solutions in the overlap regions on the next coarser level through a fine-to-coarse solution transfer. As shown in Fig. 4, at the solution points defined along the edges of the coarse cell (denoted by solid circles), the values of predicted variables are directly copied from the corresponding solution points (denoted by hollow circles) on the finer level. At the center of the coarse cell, the point-wise values (denoted by solid triangles) are evaluated to preserve the numerical conservation, i.e., it is calculated through Eq. (17) with the VIA obtained through the relation

      Figure 4.  Fine-to-coarse solution transfer in the overlap regions.

      The coarse-to-fine solution transfer is adopted to supplement the ghost cells and to evaluate the solutions within the newly generated fine cells.

      In Fig. 5, the coarse-to-fine solution transfer to supplement the ghost cells is illustrated. A coarse cell $C_{i,j,k}^{}$ is adjacent to the fine cells $C_{2i + 1,2j,k + 1}^{}$ and $C_{2i + 1,2j - 1,k + 1}^{}$ on level $k + 1$. Fine cells $C_{2i,2j,k + 1}^{}$ and $C_{2i,2j - 1,k + 1}^{}$ are ghost cells, where the solutions should be interpolated from the coarse cell $C_{i,j,k}^{}$ to provide the boundary conditions. Ghost points denoted by solid circles coincide with solution points within $C_{i,j,k}^{}$, and the values of predicted variables can be directly copied from the known solutions. Ghost points, denoted by solid squares, are not solution points of level $k$, and values of physical fields have to be evaluated through single-cell ($C_{i,j,k}^{}$) based Lagrangian interpolation polynomials, which are 2D quadratic polynomials. For some variable, $Q$, these take the form of

      Figure 5.  Coarse-fine interpolations based on multi-moments for evaluating the ghost cells.

      where the coefficients are decided by nine constraint conditions provided by known DOFs.

      The coarse-fine interpolations to evaluate predicted variables in newly generated fine cells are more complicated, as the rigorous numerical conservation should be preserved. As shown in Fig. 6, the coarse cell $C_{i,j,k}^{}$ is divided into four fine cells $C_{2i - 1,2j - 1,k + 1}^{}$, $C_{2i - 1,2j,k + 1}^{}$, $C_{2i,2j - 1,k + 1}^{}$, and $C_{2i,2j,k + 1}^{}$. The predicted variables should be evaluated at all solution points on level $k + 1$ (denoted by solid shapes). At solution points denoted by solid circles and squares, the same algorithm is adopted either by copy or interpolation, as mentioned above, to set up the fine ghost cells. However, at internal points denoted by solid triangles, the point-wise values should be determined to preserve the numerical conservation. First, the VIAs of some variable $q$ within four fine cells are evaluated by

      Figure 6.  Coarse-fine interpolations based on multi-moments for evaluating the flow variables in newly generated fine cells.

      Then the point-wise value at the cell center denoted by the solid triangle is obtained through the two-dimensional Simpson’s rule [see Eq. (17)] in each newly generated fine cell.

    4.   Numerical tests and results
    • To check the performance of the adaptive model in improving the computational efficiency, the widely used benchmark test cases are carried out in this section. The performance of the MCV dynamical core on a fixed grid has been validated in Li et al. (2013) and the numerical results are comparable to those of existing advanced models. In this study, we run the MCV dynamical core on adaptive grids with different configurations. The adaptive grid is denoted by $N \times M \times l \times \lambda $. N and M represent the number of cells along the horizontal and vertical directions respectively, which means the base block (the coarsest level) consists of $N \times M$ computational cells, the finest level is level $l$ and the refinement ratio between two adjacent levels is $\lambda $. It is noted that the minimum grid spacing is allowed by the finest level l. The normalized ${l_1}$, ${l_2}$ and ${l_\infty }$ errors (Williamson et al., 1992) and computational costs are examined. Definitions of the normalized errors are given as follows:

      where $\varOmega$ is the computational domain, $\overline q $ and $ {\overline q _{\text{e}}} $ are the numerical and reference solutions in terms of cell-integrated averages. All numerical tests in this study are carried out on the AMD EYPC 7702 CPU (single processor).

      The initial hydrostatic state is specified as

      where the Exner pressure is given by $\Pi = {\left( {\frac{p}{{{p_0}}}} \right)^{\frac{{{R_d}}}{{{c_p}}}}}$. The constants are $ {c_p} = 1004.5{\text{ J k}}{{\text{g}}^{ - 1}}{\text{ }}{{\text{K}}^{ - 1}} $ and ${R_d} = 287{\text{ J k}}{{\text{g}}^{ - 1}}{\text{ }}{{\text{K}}^{ - 1}}$.

      The basic state of potential temperature, $\theta $, is chosen to be a constant

      or specified with a constant Brunt-Väisälä frequency as

      where $N_0^2 = g {{{\rm{d}}{\text{ln}}\overline \theta }}/{{{\rm{d}}z}}$. The hydrostatic density $\overline \rho $ is then computed by:

      where $\overline T = \frac{{\overline \theta }}{{\overline \Pi }}$.

      The refinement criterion is chosen as the variation of the potential temperature perturbation for the rising thermal bubble, density current flow, and internal gravity waves, as well the variation of horizontal velocity for the Schär mountain case. Considering a two-dimensional cell $C_{i,j,k}^{}$ on level $k$, it is flagged to be refined if

      where $\delta $ is a prescribed threshold, and $q$ can be the potential temperature perturbation, $\theta '$, or the horizontal velocity, $u$.

    • Atmospheric motions caused by thermodynamic effects (Carpenter et al., 1990; Wicker and Skamarock, 1998; Ahmad and Lindeman, 2007; Norman et al., 2011) are common phenomena in nature. They are extensively used to verify the performance of dynamical cores. In this test, a thermal bubble that is warmer than the ambient air is initialized by specifying a positive potential temperature perturbation as

      where $r = \sqrt {{{\left( {x - {x_0}} \right)}^2} + {{\left( {z - {z_0}} \right)}^2}} $ and $ R = 2000 $ m (the radius of the bubble). The initial center of the thermal bubble is located at $ \left( {{x_0},{\text{ }}{z_0}} \right) = \left( {10000{\text{ m}},{\text{ }}2000{\text{ m}}} \right) $.

      In this case, the uniform background potential temperature is specified as $\overline \theta = 300{\text{ }}K$, the maximal perturbation is $\Delta \theta = 2{\text{ }}K$ and the computational domain is [0, 20] km $ \times $ [0, 10] km. The simulation runs for 1000 s and all boundaries are slippery walls. The explicit dissipation filter (Li et al., 2013) with a viscosity coefficient, $\mu = 10{\text{ }}{{\text{m}}^2}{\text{ }}{{\text{s}}^{ - 1}}$, is used to eliminate numerical noise. The threshold for refinement is set as $\delta = 0.04$. During the simulation, the thermal bubble rises and finally attains a mushroom-like shape. The numerical results on a uniform grid, with a resolution of 25 m $ \times $ 25 m, are used as a reference solution to calculate the normalized errors.

      Normalized errors and elapsed CPU time for different grids are given in Table 1, where the results are grouped according to the finest resolution on an adaptive grid. The normalized CPU time is also computed by dividing the CPU time on the coarsest uniform grid. From Table 1, it can be found that the normalized errors of the numerical results on the uniform and AMR grids within each group are quite similar, noting that much less CPU time is consumed by the AMR model. Contour plots of the potential temperature perturbation $\left( {\theta '} \right)$ and the absolute errors on a $100 \times 50 \times 3 \times 2$ grid at different times are shown in Fig. 7. The solid thick lines represent the interfaces between different levels. The symmetric distribution of $\theta '$ is perfectly reproduced on the adaptive grid as in our previous studies (Li et al., 2013). It is observed that the fine levels are dynamically adjusted to fit the change of flow field. AMR grids properly locate the disturbed region and put the fine blocks there to assure numerical accuracy. In undisturbed regions, the coarse blocks are applied to save computational costs. The relative total mass error along the simulation is also recorded, which supports that the numerical conservation property is well achieved through the flux correction procedure.

      Resolution (m)Grid configuration$l_1$$l_2$$ {l_\infty }$ CPU time (s)
      200 × 200100 × 50 × 1 × 10.230.210.281 (187.34)
      100 × 100200 × 100 × 1 × 16.40 × 10−26.86 × 10−20.118.45 (1584.21)
      100 × 50 × 2 × 26.44 × 10−26.87 × 10−20.112.84 (532.55)
      50 × 50400 × 200 × 1 × 18.82 × 10−31.09 × 10−22.43 × 10−275.00 (14051.62)
      200 × 100 × 2 × 29.34 × 10−31.15 × 10−22.43 × 10−210.17 (1905.17)
      100 × 50 × 2 × 48.90 × 10−31.07 × 10−22.53 × 10−219.54 (3660.69)
      100 × 50 × 3 × 29.28 × 10−31.13 × 10−22.35 × 10−213.63 (2544.14)
      25 × 25800 × 400 × 1 × 1631.39 (118285.12)

      Table 1.  Normalized errors of the potential temperature perturbation ($\theta '$) and elapsed CPU time for the thermal bubble test running on different grids.

      Figure 7.  Contour plots of a potential temperature perturbation ($\theta '$) and its errors for a rising thermal bubble on $100 \times 50 \times 3 \times 2$ grid at (a) t = 0 s, (b) t = 250 s, (c) t = 500 s, (d) t = 750 s, (e) t = 1000 s, and (f) absolute errors.

    • In the density current test case (Straka et al., 1993; Giraldo and Restelli, 2008), a cold bubble is put in a neutrally stratified atmosphere by specifying a negative potential temperature perturbation as

      where

      $ \overline \theta = 300{\text{ K}} $, $\Delta \theta = - 15{\text{ K}}$, $\left( {{x_0},{\text{ }}{z_0}} \right) = \left( {0 , 3000 } \right)\;{\rm{m}}$ and $\left( {{x_r},{\text{ }}{z_r}} \right) = \left( {4000 , 2000 } \right)\;{\rm{m}}$.

      During the simulation, the cold bubble drops to the ground and spreads out in the horizontal direction to form three Kelvin-Helmholtz shear instability rotors along the cold frontal surface. The computational domain of this case is [–26.5, 26.5] km $ \times $ [0, 6.4] km and the simulation time is 900 s. All of the boundaries are slippery walls. A dissipation filter with a coefficient of $\mu = 10{\text{ }}{{\text{m}}^2}{\text{ }}{{\text{s}}^{ - 1}}$ is used to meet the requirement of a physical process. The threshold of the refinement criterion of this case is $\delta = 0.2$. Again, the numerical results on a uniform grid with a resolution of 25 m × 25 m are adopted as reference solutions.

      Normalized errors of $\theta '$ and elapsed CPU time on different grids are given in Table 2. It can be observed that the AMR model consumes much less computational costs at the price of slightly larger normalized errors. Contour plots of $\theta '$ and the absolute errors on a $132 \times 32 \times 2 \times 4$ grid are shown in Fig. 8, which have the finest resolution of 50 m $ \times $ 50 m. Only half of the computational domain is depicted due to the symmetry of the flow field. The numerical result agrees well with those on the uniform grid with the same resolution given in Li et al., (2013).

      Resolution (m)Grid configuration${l_1}$${l_2}$$ {l_\infty } $CPU time (s)
      200 × 200132 × 32 × 1 × 10.280.330.941 (235.36)
      100 × 100265 × 64 × 1 × 10.130.170.608.69 (2045.44)
      132 × 32 × 2 × 20.130.170.643.44 (809.19)
      50 × 50530 × 128 × 1 × 12.63 × 10−23.59 × 10−20.1569.93 (16458.58)
      265 × 64 × 2 × 22.72 × 10−23.67 × 10−20.1711.10 (2613.85)
      132 × 32 × 2 × 42.64 × 10−23.59 × 10−20.1525.99 (6117.59)
      132 × 32 × 3 × 22.68 × 10−23.66 × 10−20.1820.00 (4047.29)
      25 × 251060 × 256 × 1 × 1640.80 (150820.42)

      Table 2.  Normalized errors of the potential temperature perturbation ($\theta '$) and elapsed CPU time for the density current test running on different grids.

      Figure 8.  Contour plots of the potential temperature perturbation ($\theta '$) and its errors for the density current on $132 \times 32 \times 2 \times 4$ grid at (a) t = 0 s, (b) t = 225 s, (c) t = 450 s, (d) t = 675 s, (e) t = 900 s, and (f) absolute errors.

    • The internal gravity waves test involves the evolution of a potential temperature perturbation in a channel with periodic boundary conditions in the horizontal direction. Initial conditions used in Skamarock and Klemp (1994) are adopted. The potential temperature field is initialized as

      where $ {\theta _0} = 300{\text{ K}} $, $H = 10{\text{ km}}$, $\Delta \theta = 0.01{\text{ K}}$, ${x_0} = 100{\text{ km}}$, and $a = 5{\text{ km}}$.

      The background state of potential temperature $ \overline \theta \left( z \right) $ is obtained by using a constant Brunt-Väisälä frequency. A constant mean flow of $\overline u = 20{\text{ m }}{{\text{s}}^{ - 1}}$ for the advection of the internal gravity waves is adopted. The bottom and top boundaries are slippery walls. The computational domain of this case is [0, 300] km $ \times $ [0, 10] km, and the simulation time is 3000 s. $\delta = 1.8 \times {10^{ - 4}}$ is chosen for grid refinement. Numerical results on a 250 m $ \times $ 25 m grid are taken as the reference solutions. It is noted that the computational cells are no longer square, an aspect ratio of grid spacing $\Delta x = 10\Delta z$ is adopted in this case.

      Normalized errors of $\theta '$ and elapsed CPU times on various grids are given in Table 3. Comparing to the first two test cases, the effect of AMR grids in saving computational costs is less obvious. The reason is that the internal gravity waves are spreading in the horizontal direction during the simulation and the disturbed regions are continuously growing as shown in Fig. 9. Maximum and minimum values of vertical velocity and potential temperature perturbation on the AMR grids with the finest resolution of $\Delta z = 100{\text{ m}}$after 3000 s are given in Table 4, which is the same as that obtained by Li et al (2013). The distribution of the absolute errors of $\theta '$ on a uniform grid and a 75 × 25 × 3 × 2 grid at $t = 3000{\text{ s}}$ is depicted in Figs. 10a, b. No obvious increase of errors is found around the boundaries between different levels, which reveals that the computational mode due to the change of grid resolution is effectively suppressed in the proposed AMR model. The general errors in the AMR model are also affected by the prescribed refinement threshold. More accurate solutions are obtained when a larger area of the computational domain is refined with a smaller refinement threshold.

      Resolution (m)Grid configuration${l_1}$${l_2}$$ {l_\infty } $CPU time (s)
      4000 × 40075 × 25 × 1 × 10.250.290.321 (48.60)
      2000 × 200150 × 50 × 1 × 18.20 × 10−20.110.159.77 (474.99)
      75 × 25 × 2 × 27.95 × 10−20.110.164.13 (200.48)
      1000 × 100300 × 100 × 1 × 11.48 × 10−22.12 × 10−23.58 × 10−278.87 (3833.06)
      150 × 50 × 2 × 22.77 × 10−23.60 × 10−26.09 × 10−231.37 (1524.60)
      75 × 25 × 2 × 41.62 × 10−22.14 × 10−23.57 × 10−227.58 (1340.39)
      75 × 25 × 3 × 22.84 × 10−23.61 × 10−26.08 × 10−226.33 (1279.57)
      250 × 251200 × 400 × 1 × 16375.50 (309849.29)

      Table 3.  Normalized errors of potential temperature perturbation ($\theta '$) and elapsed CPU time for internal gravity waves test running on different grids.

      Figure 9.  Contour plots of potential temperature perturbation ($\theta '$) for the internal gravity waves on $75 \times 25 \times 3 \times 2$ grid at (a) t = 0 s, (b) t = 750 s, (c) t =1500 s, (d) t = 2250 s, and (e) t = 3000 s.

      Gridwmax (m s−1)wmin (m s−1)$\theta {'_{\max }}$ (K)$ \theta {'_{\min }} $ (K)
      75 × 25 × 2 × 42.47 × 10−3−2.42 × 10−32.80 × 10−3−1.52 × 10−3
      75 × 25 × 3 × 22.47 × 10−3−2.42 × 10−32.80 × 10−3−1.52 × 10−3
      150 × 50 × 2 × 22.47 × 10−3−2.43 × 10−32.80 × 10−3−1.52 × 10−3

      Table 4.  Maximum and minimum of vertical velocity, $w$, and potential temperature perturbation, $\theta '$ , for the internal gravity waves test on the AMR grids with the finest resolution, $\Delta z = 100{\text{ m}}$ after 3000 s.

      Figure 10.  Error distribution of the potential temperature perturbation ($\theta '$) for the internal gravity waves on the (a) uniform grid and (b) $75 \times 25 \times 3 \times 2$ grid at t = 3000 s.

    • The Schär mountain case (Schär et al., 2002) is used to test the ability of a model to deal with the effects of complex terrain. A mountain with five peaks is used as topography which is specifically defined as

      where ${h_0} = 250{\text{ m}}$, ${a_0} = 5000{\text{ m}}$, and ${\lambda _0} = 4000{\text{ m}}$. An initial background state of potential temperature is obtained by using a constant Brunt-Väisälä frequency of ${N_0} = {10^{ - 2}}{\text{ }}{{\text{s}}^{ - 1}}$ and ${\theta _0} = 280{\text{ K}}$. A constant mean flow of $\overline u = 10{\text{ m }}{{\text{s}}^{ - 1}}$along the horizontal is also initialized. This case is running on a domain of [–25, 25] km $ \times $ [0, 21] km for 10 hours. Non-reflecting boundary conditions are implemented by putting sponge layers in the area of $z \geqslant 9$ km for the top boundary and $ \left| x \right| \geqslant 15 $ km for the lateral inflow and outflow boundaries. The bottom boundary is a slippery wall. A refinement threshold, $\delta = 0.3$, is applied to this case.

      The profile of the Schär mountain is depicted in Fig. 11e, noting that only part of the computational domain is displayed. Contour plots of the horizontal velocity, $u$, at various stages in the simulations are displayed in Figs. 11ad. Since the disturbed regions quickly extended to cover the entire computational domain, it is not expected that the computational costs can be significantly saved by AMR models. Numerical results on a uniform fine grid with a resolution of 125 m × 105 m are used as reference solutions and the corresponding normalized errors are given in Table 5. The elapsed CPU time of the AMR model is about 44% of that on the uniform grid. The AMR grid does not affect accuracy since the numerical errors after 10 hours are almost the same as those on the uniform grid with the same resolution.

      Figure 11.  Contour plots of horizontal velocity ( u ) of the Schär mountain case at (a) t = 7200 s, (b) t = 16800 s, (c) t = 26400 s, (d) t = 36000 s, and (e) the shape and size of the Schär mountain.

      Resolution (m)Grid configuration${l_1}$${l_2}$$ {l_\infty } $CPU time (s)
      500 × 420100 × 50 × 1 × 12.97 × 10−35.66 × 10−32.19 × 10−21 (4028.53)
      250 × 210200 × 100 × 1 × 14.31 × 10−48.66 × 10−49.41 × 10−38.66 (34899.12)
      100 × 50 × 2 × 28.61 × 10−41.28 × 10−31.03 × 10−23.83 (15429.89)
      125 × 105400 × 200 × 1 × 183.10 (334779.90)

      Table 5.  Normalized errors of horizontal velocity ($u$) and elapsed CPU time for the Schär mountain test running on different grids.

    5.   Conclusion
    • In this paper, an adaptive nonhydrostatic atmospheric dynamical core is proposed by using the multi-moment constrained, finite-volume method and Berger-Oliger’s adaptive mesh refinement algorithm. The high-order scheme is constructed based on two kinds of local degrees of freedom. As a result, the spatial reconstruction for a 3rd-order scheme is accomplished based on a single-cell-based stencil. The compact spatial stencil is beneficial for the efficient implementation of coarse-fine interpolation and reduction of the extra numerical errors due to the non-uniform resolution of an AMR grid. The proposed model runs on the adaptive grid consisting of blocks with different resolutions. The blocks with fine resolutions are placed in those regions that contain large variations of physical fields in order to improve the numerical accuracy, while the blocks with coarse resolutions are placed in the remainder of the area to save computational costs. With the evolution of flow fields, the distributions of blocks are dynamically adjusted to guarantee that the regions with complex flow structures are always resolved with fine resolutions. Numerical results of widely used benchmark tests reveal the effectiveness of the AMR technique in saving computational costs. It is observed that the use of an adaptive grid has only a very slight negative impact on computational accuracy in comparison with results on fixed grids with the same resolution.

      Today, finer predictions of small-scale weather and climate phenomena and their interactions with large-scale ones gain more and more attention. Though the supercomputing capacity has been greatly increased recently, the use of an AMR grid would be beneficial in reducing unnecessary computational cells and make it more affordable to accurately track spatially complex and time-dependent small-scale phenomena. However, using an AMR grid has the undesirable side-effect of making it more difficult to achieve high scalability, since it is not easy to keep satisfying load balance among different CPUs on time-varying computational meshes. It is a major challenge to construct a practical atmospheric model using adaptive mesh.

      Promising results by the proposed adaptive model have shown great potential in applying the proposed numerical framework to develop a more efficient atmospheric dynamical core. Future research will mainly focus on the extension of an adaptive model to the three-dimensional regional and global dynamical cores.

      Acknowledgements. This work is supported by The National Key Research and Development Program of China (Grants Nos. 2017YFA0603901 and 2017YFC1501901) and The National Natural Science Foundation of China (Grant No. 41522504).

Reference

Catalog

    /

    DownLoad:  Full-Size Img  PowerPoint
    Return
    Return