Stéphane Abide
email: stephane.abide@univ-perp.fr

Stéphane Viazzo

Isabelle Raspo

Anthony Randriamampianina

Higher-order compact scheme for high-performance computing of stratified rotating flows

Keywords:

a b s t r a c t

To take advantage of modern generation computing hardware, a scalable numerical method, based on higher-order compact scheme, is described to solve rotating stratified flows in cylindrical annular domains. An original approach combining 2d-pencil decomposition and reduced Parallel Diagonal Dominant is proposed to improve the parallelization performance during the computation of Poisson/Helmholtz solvers and time explicit terms. The developed technique is validated with respect to analytical solutions, using the method of manufactured solutions, and available data for two specific configurations. The purpose is to demonstrate its ability to correctly capture the flow characteristics in strato-rotational instability and in baroclinic instability with associated small-scale features. Moreover, this code is found to drastically reduce the huge execution times often preventing detailed numerical investigations of these complex phenomena. Strong scaling test is carried out to assess the performance for up to 1024 cores using grid up to 128 × 568 × 568 in radial, axial and azimuthal directions.

Introduction

Current trends in computational science stem from the increasing prominence of modern multi-and many-core hardware such as graphics processing units (GPUs) in high-performance computing (HPC) infrastructure. Improved multi-core central processing units (CPUs) associated with GPU-accelerated computing can greatly accelerate scientific studies especially in large-scale simulations [START_REF] Yakovlev | To CG or to HDG: a comparative study in 3d[END_REF][START_REF] Witherden | PYFR: an open source framework for solving advection-diffusion type problems on streaming architectures using the flux reconstruction approach[END_REF][START_REF] King | Exploiting batch processing on streaming architectures to solve 2d elliptic finite element problems: a hybridized discontinuous galerkin (HDG) case study[END_REF] . However, the reduction of execution times strongly depends on the optimization of the performance of discretization methods used on each node while good computational scaling results from small communication footprint [START_REF] Moxey | Optimising the performance of the spectral/hp element method with collective linear algebra operations[END_REF] .

Even though finite element methods are becoming commonplace for simulations of most engineering and biomedical applications on complex geometries, low-order discretizations suffer from strong mesh refinement required to capture certain specific complex solution features. In contrast, high-order techniques, such as spectral approaches, provide improved numerical characteristics at reduced computational cost for a given number of degrees of freedom but they usually remain restricted to regular configurations. Recently, the spectral/hp element methods [START_REF] Karniadakis | Spectral/hp element methods for computational fluid dynamics[END_REF][START_REF] Cantwell | Nektar ++ : an open-source spectral/hp element framework[END_REF][START_REF] Fischer | Nek50 0 0: Open Source Spectral Element CFD Solver[END_REF] combine high accuracy, known as the spectral convergence, and geometric flex-ibility to tackle challenging aeronautical flow simulations using between O (10 4 ) and O (10 5 ) nodes [START_REF] Lombard | Implicit large-eddy simulation of a wingtip vortex[END_REF] , and even up to O (10 6 ) nodes [START_REF] Kerkemeier | scalability of the nek50 0 0 spectral element code[END_REF] .

In the present study, a higher-order compact scheme solver [START_REF] Abide | An efficient parallel high-order compact scheme for the 3d incompressible Navier-Stokes equations[END_REF] is tailored for solving complex rotating stratified flows in annular domains to take advantage of modern generation computing hardware. Our goal is to reduce the huge total time due to the very long transient stage resulting from the thermal stratification associated with the rotation. We have chosen to treat two specific idealized configurations: a Taylor-Couette setup submitted to axial stratification for the strato-rotational instability (SRI) [START_REF] Gellert | Stratorotational instability in Taylor-Couette flow heated from above[END_REF] and a uniformly rotating cavity submitted to radial thermal gradient for the baroclinic instability [START_REF] Hide | An experimental study of thermal convection in a rotating liquid[END_REF] . In the first case, the geometry is defined by a large aspect ratio while the second one is characterized by the coexistence of very different spatio-temporal scales [START_REF] Randriamampianina | Inertia-gravity waves in a liquidfilled, differentially heated, rotating annulus[END_REF] . Before considering more deeply actual physics arising in these two configurations, this work reports on the validation of the developed numerical tool with respect to data available in the literature. We have restricted our simulations to remain in the domain of the Boussinesq approximation: weak density stratification and weak rotation rate to neglect the centrifugal acceleration.

Most of common discretization such as finite volumes, finite differences, finite elements or lattice Boltzmann method, is efficiently implemented in parallel computing environment, even though massively parallel linear solvers still remain a challenging work [START_REF] Schornbaum | Massively parallel algorithms for the lattice Boltzmann method on nonuniform grids[END_REF][START_REF] Vàzquez | Alya: multiphysics engineering simulation toward exascale[END_REF][START_REF] Baker | Scaling algebraic multigrid solvers: on the road to exascale[END_REF][START_REF] Notay | An aggregation-based algebraic multigrid method[END_REF] . Concerning high accuracy discretization, the picture is more mixed. The computational stencils of such discretizations are larger than those of the common lower-order formulations, therefore inducing large parallel communications. The 2d-pencil decomposition became the common strategy to deal with large computational stencils [START_REF] Laizet | A numerical strategy to combine high-order schemes, complex geometry and parallel computing for high resolution DNS of fractal generated turbulence[END_REF][START_REF] He | A high order finite difference solver for massively parallel simulations of stably stratified turbulent channel flows[END_REF][START_REF] Thais | A massively parallel hybrid scheme for direct numerical simulation of turbulent viscoelastic channel flow[END_REF][START_REF] Chatterjee | Scaling of a fast fourier transform and a pseudo-spectral fluid solver up to 196608 cores[END_REF][START_REF] Shi | A hybrid mpi-openmp parallel implementation for pseudospectral simulations with application to taylorcouette flow[END_REF] . The successful example of the massively parallel FFT computations with 1d then 2d pencil decomposition [START_REF] Chatterjee | Scaling of a fast fourier transform and a pseudo-spectral fluid solver up to 196608 cores[END_REF] has been tailoring to the Navier-Stokes solvers based on higher-order discretizations [START_REF] Laizet | A numerical strategy to combine high-order schemes, complex geometry and parallel computing for high resolution DNS of fractal generated turbulence[END_REF][START_REF] Mininni | A hybrid MPI-Openmp scheme for scalable parallel pseudospectral computations for fluid turbulence[END_REF] . Nowadays, this strategy of parallelization offers a significant reduction of the wall-time and allows the study of highly turbulent flows. Because of the expensive cost of the communications inherent to 2dpencil decompositions, some choices in the design of the numerical method must be made to limit communications. For instance, one can choose time-explicit advancement schemes for the diffusive terms, Fourier series expansions, or a collocated variable layout [START_REF] Laizet | A numerical strategy to combine high-order schemes, complex geometry and parallel computing for high resolution DNS of fractal generated turbulence[END_REF] . Compact schemes are finite differences based on implicit relations [START_REF] Lele | Compact finite difference schemes with spectral-like resolution[END_REF] which rely on solutions of tridiagonal or pentadiagonal linear systems. Therefore, higher-order compact schemes can be considered as an interesting alternative to spectral methods, associating the robustness of finite difference with improved accuracy. In a recent work, Abide et al. [START_REF] Abide | An efficient parallel high-order compact scheme for the 3d incompressible Navier-Stokes equations[END_REF] proposed an approximate parallel tridiagonal solver [START_REF] Sun | Application and accuracy of the parallel diagonal dominant algorithm[END_REF] to avoid the use of 2d-pencil decomposition for the computation of compact derivatives and interpolations. In this case, the communications between neighbor process rely on halo exchange algorithm. Nevertheless, the 2d-pencil decomposition is still considered to tailor a parallel version of the full diagonalization method [START_REF] Haldenwang | Chebyshev 3d spectral and 2d pseudospectral solvers for the Helmholtz equation[END_REF] , despite the cost of the communications involved by the data transpositions. The combination of the two parallelization strategies allows the authors to successfully simulate incompressible flow in truly three-dimensional cartesian domain. Moreover, the basic numerical method features for turbulent flow simulations such as the staggered grid and implicit viscous terms are conserved, while exhibiting a good strong scaling.

The proposed numerical method is based on Fourier/Fourthorder compact scheme discretization defined on staggered grid. The second order Adams-Bashforth/Backward Euler time scheme combined with an efficient projection scheme has been considered: this approach relaxes the stringent time step constraints inherent to the diffusive terms [START_REF] Hugues | An improved projection scheme applied to pseudospectral methods for the incompressible Navier-Stokes equations[END_REF] . It corresponds to an extension of a previous numerical code, developed in cartesian coordinates [START_REF] Abide | An efficient parallel high-order compact scheme for the 3d incompressible Navier-Stokes equations[END_REF] , for solving rotating stratified flows in annular domains. Moreover, the present parallelization strategy is concerned with a 2d-pencil decomposition and reduced Parallel Diagonal Dominant (rPDD) to address the numerical solutions of Helmholtz/Poisson problems and the computations of the compact derivatives and interpolations respectively. The implementation of the combination of these two techniques brings novelty and originality to the present approach in comparison with previous works based on higher-order codes [START_REF] Laizet | A numerical strategy to combine high-order schemes, complex geometry and parallel computing for high resolution DNS of fractal generated turbulence[END_REF] .

The rotating stratified flows concern a wide spectrum of applications as well as in geophysics such as large scale circulations in atmosphere and oceans, astrophysics such as accretion disk, than in disk storage, turbomachinery … The Strato-Rotational Instability (SRI) is a purely hydrodynamic instability with distinctive local features and may be studied from specifically designed laboratory experiments and numerical simulations in an axially-stratified Taylor-Couette setup [START_REF] Gellert | Stratorotational instability in Taylor-Couette flow heated from above[END_REF][START_REF] Bars | Experimental analysis of the stratorotational instability in a cylindrical Couette flow[END_REF][START_REF] Rüdiger | Magnetic processes in astrophysics: theory, simulations, experiments[END_REF] . Here axial stratification is obtained by cooling the bottom and heating the upper lid [START_REF] Gellert | Stratorotational instability in Taylor-Couette flow heated from above[END_REF] . The temperature stratification represents the axial stratification of an accretion disk centered around a hot star and the rotation of the cylinders mimics the astrophysical rotation laws with rotation rate decreasing outwards. Understanding the mechanisms that can result in an outward angular momentum transport is the central problem of planet formation, particularly in the theory of accretion disks. When a planet forms in a disk, angular momentum has to be carried away from the planet otherwise its rotation speed would be far too large. Only turbulence can achieve such a large angular momentum transport. Recent studies mentioned that, at fixed value of the Froude number, F r = i /N, with N = gα∂ T /∂ z the buoyancy frequency and i the angular velocity of the inner cylinder, the system is more stable when increasing the Reynolds number [START_REF] Ibanez | Observations of the stratorotational instability in rotating concentric cylinders[END_REF][START_REF] Rüdiger | The stratorotational instability of Taylor-Couette flows with moderate Reynolds numbers[END_REF] . The final objective is to eventually explore the not yet understood turbulent flows, to get closer to actual astrophysical applications, by taking advantage of optimized HPC.

Re = i R i (R o -R i ) /ν
Baroclinic instability is recognized to be one of the dominant energetic processes in the large-scale atmospheres of terrestrial planets, such as the Earth and Mars [START_REF] Pierrehumbert | Baroclinic instability[END_REF] , and in the oceans [START_REF] Orlanski | Baroclinic instability in ocean currents[END_REF] . Its time-dependent behaviour exerts a dominant influence on the intrinsic predictability of the atmospheric circulation and the degree of chaotic variability in its large-scale meteorology [START_REF] Pierrehumbert | Baroclinic instability[END_REF][START_REF] Read | Wave interactions and baroclinic chaos: a paradigm for long timescale variability in planetary atmospheres[END_REF][START_REF] Read | Transition to geostrophic turbulence in the laboratory, and as a paradigm in atmospheres and oceans[END_REF] . Inertia gravity waves (IGWs) are ubiquitous in the atmosphere and the oceans and are known to play a fundamental role in a wide variety of processes. The contributions mainly concern the transport of a significant amount of energy and momentum, the initiation and organization of convection, the induction and modulation of turbulence, as well as the modification of the mean circulation and thermal structure of atmospheric and oceanic motions [START_REF] Fritts | Gravity wave dynamics and effects in the middle atmosphere[END_REF] . Observations and simulations have revealed their spontaneous occurrence during the development of baroclinic instability. In spite of intensive research activities carried out over the last decades, the generation mechanism and the propagation of IGWs, as well as their interaction with large-scale structures, remain poorly understood. A better understanding of these phenomena is therefore mandatory for the improvement of the IGW parameterization schemes actually required to upgrade numerical global weather predictions. Since the pioneering works of Hide [START_REF] Hide | An experimental study of thermal convection in a rotating liquid[END_REF] , the differentially heated, rotating cylindrical annulus has been an archetypal means of studying the properties of fully-developed baroclinic instability in the laboratory. Recent works mainly based on direct numerical simulations reported the occurrence of IGWs in water-filled baroclinic cavity along the inner cold cylinder [START_REF] Borchert | Gravity wave emission in an atmospherelike configuration of the differentially heated rotating annulus experiment[END_REF][START_REF] Hien | Spontaneous inertia-gravity wave emission in the differentially heated rotating annulus experiment[END_REF][START_REF] Larcher | Instabilities and small-scale waves within the Stewartson layers of a thermally driven rotating annulus[END_REF] . Moreover, Von Larcher et al. [START_REF] Larcher | Instabilities and small-scale waves within the Stewartson layers of a thermally driven rotating annulus[END_REF] observed the presence of additional small-scale ripples resulting from hydrodynamical instability along the hot wall.

The paper is organized as follows. First, the physical and mathematical models of the SRI and Baroclinic configurations are detailed. The numerical methods to solve the coupled Navier-Stokes and energy equations in annuli are introduced before describing the parallelizations strategies. The last section is devoted to the validation of the new numerical solver. Accuracy will be checked, validations with existing findings for strato-rotational and baroclinic instabilities are performed including a strong scaling test.

Physical and mathematical models

Physical model

This study focuses on flows in rotating annular cavities filled with an incompressible fluid and submitted to either a radial or a vertical thermal gradient. In these configurations, the flow results from the competition between the rotation of the cylinders and the buoyancy force arising from the prescribed temperature gradient.

The two configurations considered together with the notations are depicted in Fig. 1 .

The first configuration corresponds to a Taylor-Couette system ( Fig. 1 a) where the fluid is submitted to a stable vertical thermal gradient, which means that the top boundary is heated whereas the bottom one is cooled. The inner cylinder rotates faster than the outer one. To reduce the edge effects, stress-free conditions are prescribed for velocity on the two horizontal disks. Despite the stabilizing effect of the vertical stratification induced by the stable thermal gradient, an hydrodynamic instability, named the strato-rotational instability, may develop for specific values of the

μ = o i , Re = i dR i ν , F r = i N (1)
with d = R o -R i , ν the kinematic viscosity, i and o the angular speeds of the inner and outer cylinders respectively, and N the Brunt-Väisälä frequency which is a measure of the vertical stratification scale. In the framework of the Boussinesq approximation (see section 2.2), N is defined by:

N 2 = αg T H (2)
with α the thermal expansion coefficient.

In the second configuration ( Fig. 1 b), a radial temperature gradient is prescribed by heating the outer cylinder and cooling the inner one. Unlike the previous case, the two vertical cylinders rotate with the same angular velocity . The top boundary is an open free surface. This configuration was found to be a suitable testbed for studying mid-latitude atmospheric flows where a baroclinic instability together with inertia gravity waves may develop. The classical control parameters are the Taylor number Ta and the Rossby number Ro [START_REF] Hide | An experimental study of thermal convection in a rotating liquid[END_REF] , defined by:

T a = 4 2 d 5 ν 2 H , Ro = α T gH 2 d 2 (3)

Mathematical model

In both configurations, the flow and heat transfer are governed by the incompressible Navier-Stokes equations coupled with energy equation through the Boussinesq approximation. With these assumptions, the set of equations reads:

∇. u = 0 in D ∂ t u + 1 2 [ (u . ∇ ) u + ∇. (uu ) ] = -∇ p + ν u + F in D ∂ t T + 1 2 [ (u . ∇ ) T + ∇. (u T ) ] = κ∇ 2 T in D (4)
D being the computational domain, κ the fluid thermal diffusivity and u = (u r , u θ , u z ) , p and T standing for the velocity, pressure and temperature fields, respectively. The vectorial Laplacian

Table 1

Source terms F and boundary conditions in the SRI and baroclinic flows ( T h and T c stand for the hot and cold temperatures).

Vertical ∇T Radial ∇T F α(T -T 0 ) g α(T -T 0 ) g + 2 × u + × ( × r )
u = ∇ 2 (u r ) -u r /r 2 -(2 /r 2 ) ∂ θ u θ ∇ 2 (u θ ) -u θ /r 2 + (2 /r 2 ) ∂ θ u r ∇ 2 (u z ) (5)
where ∇ 2 (.) stands for the scalar Laplacian defined by:

∂ 2 r + (1 /r) ∂ r + (1 /r 2 ) ∂ 2 θ + ∂ 2 z ( 6
)
In Eq. ( 4) , the source term F accounts for external forces, mainly the buoyancy force. For the first case corresponding to a Taylor-Couette configuration ( Fig. 1 a), the equations are expressed in a fixed frame. For the second configuration with a radial thermal gradient ( Fig. 1 b), the equations are solved in a rotating reference frame and, as a consequence, in addition to the buoyancy force, Coriolis and centrifugal forces are also included in F . The expression of F and the boundary conditions are given in Table 1 for the two configurations.

Numerical method

Time discretization and projection scheme

The second order semi-implicit Adams-Bashforth/Backward-Euler scheme is used for the time advancement of Eq. ( 4) . The diffusive terms are implicitly treated except for the crossed terms of the vectorial Laplacian, resulting from cylindrical coordinates, which are explicitly evaluated. The temporal scheme reads:

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ ∇ • u n +1 = 0 1 2 t 3 u n +1 -4 u n + u n -1 + 2 H n u -H n -1 u = -∇ p n +1 + ν u n +1 + F n + 1 1 2 t 3 T n +1 -4 T n + T n -1 + 2 H n T -H n -1 T = κ∇ 2 T n +1 (7)
where H u and H T account for non-linear convective terms. For the momentum equations, the crossed contributions in the diffusive term are also included in H u . The velocity-pressure coupling is solved by a projection algorithm. Specifically, we used the improved method proposed by Hugues and Randriamampianina [START_REF] Hugues | An improved projection scheme applied to pseudospectral methods for the incompressible Navier-Stokes equations[END_REF] which allows a temporal evolution of the normal pressure gradient at boundaries. This projection algorithm is detailed in [START_REF] Hugues | An improved projection scheme applied to pseudospectral methods for the incompressible Navier-Stokes equations[END_REF][START_REF] Raspo | A spectral projection method for the simulation of complex three -dimensional rotating flows[END_REF] . Thus, the discretized Boussinesq-Navier-Stokes equations ( 7) reduce to two Poisson problems for a preliminary pressure and a pressure correction respectively, and four Helmholtz problems arising from the momentum and energy equations. The accuracy and the computational efficiency rely on the choice of the space discretization. In this work, we choose the spectral discretization in the azimuthal direction, using Fourier series, and the compact finite differences scheme in the two other directions. A short description is given hereafter.

Compact scheme discretization

The spatial approximation in the non-periodic r and z directions of the different terms of Eq. ( 4) is based on a fourth-order compact scheme discretization [START_REF] Lele | Compact finite difference schemes with spectral-like resolution[END_REF][START_REF] Abide | A 2d compact fourth-order projection decomposition method[END_REF][START_REF] Knikker | Study of a staggered fourth-order compact scheme for unsteady incompressible viscous flows[END_REF] . Compact schemes have a better resolution at high wavenumber than explicit finite differences. However, compact schemes induce a computational overhead due to the solution of tridiagonal/pentadiagonal linear systems. This overhead remains acceptable since the algorithm complexity is linear. A staggered grid arrangement of the velocity and the pressure is considered to ensure some conservation properties [START_REF] Knikker | Study of a staggered fourth-order compact scheme for unsteady incompressible viscous flows[END_REF][START_REF] Schiestel | A hermitian-fourier numerical method for solving the incompressible Navier-Stokes equations[END_REF] . For instance, in the radial direction, the grid in the interval ( R i , R o ) is composed of the pressure nodes, r i +1 / 2 = R i + ( i -1 / 2 ) h r with 0 ≤ i ≤ n r + 1 , and of the velocity nodes r i = ih r with 0 ≤ i ≤ n r . The uniform space step is h r = (R o -R i ) /n r . Non-uniform grids are obtained by means of an analytical mesh transformation, introducing the metric coefficients which are analytically computed [START_REF] Abide | A 2d compact fourth-order projection decomposition method[END_REF] . The staggered fourth-order compact scheme approximations from the pressure to velocity nodes, reads for the staggered derivation: [START_REF] Lombard | Implicit large-eddy simulation of a wingtip vortex[END_REF] and for the staggered interpolation:

α f i -1 + f i + α f i +1 = a f i +1 / 2 -f i -1 / 2 /h r
α f i -1 + f i + α f i +1 = a f i +1 / 2 + f i -1 / 2 / 2 (9)
where the coefficients α and a are calculated from accuracy requirements [START_REF] Lele | Compact finite difference schemes with spectral-like resolution[END_REF] . The boundary condition relations are: [START_REF] Abide | An efficient parallel high-order compact scheme for the 3d incompressible Navier-Stokes equations[END_REF] The implicit nature of this finite difference scheme is expressed by the equivalent matrix formulation M f (p) = B f of Eqs. ( 9) and [START_REF] Abide | An efficient parallel high-order compact scheme for the 3d incompressible Navier-Stokes equations[END_REF] where superscript p refers to either derivation ( ) p = 1 or interpolation p = 0 , and M refers to a square tridiagonal matrix. All the coefficients defining the fourth-order accurate compact scheme operators used in this study are detailed in Table 2 .

f (p) 1 + α f (p) 2 = a f -1 / 2 + b f 1 / 2 + c f 3 / 2 f (p) n + α f (p) n -1 = a f n +1 / 2 + b f n -1 / 2 + c f n -3 / 2

The full diagonalization method

This section presents the full diagonalization method used for solving the Helmholtz/Poisson problems ( Eq. ( 7) ) which is an efficient direct solver for linear systems arising from higher-order discretizations [START_REF] Haldenwang | Chebyshev 3d spectral and 2d pseudospectral solvers for the Helmholtz equation[END_REF] . The main advantage of this approach is to eliminate the bad condition number of the operator matrices responsible for high computational costs with iterative solvers [START_REF] Abide | A fourth-order iterative solver for the singular poisson equation[END_REF] . However, it is worth mentioning the three major drawbacks of the full diagonalization method which are the high algorithm complexity, the limitation to constant coefficients, and the required tensorial framework. For the sake of clarity of the following sections, this linear solver is detailed in the case only for a Poisson equation in cylindrical coordinates. With the tensorial notations and using Fourier series expansion, we have to compute for each Fourier mode i the following 2D solution:

D 2 ir I z + I r D 2 z ˆ φ i = ˆ s i , 0 ≤ i ≤ n θ / 2 ( 11
)
where the operators D 2 ir and D 2 z are defined by:

D 2 ir = 1 r δ f c r r δ c f r - i 2 r 2 , D 2 z = δ f c z δ c f z ( 12
)
and ˆ s i represents the i th mode of the 1D Fourier transform in the azimuthal direction of the source term s denoted by F θ (s ) .

It should be noted that these operators are modified to include boundary conditions [START_REF] Abide | A 2d compact fourth-order projection decomposition method[END_REF] and that the radial operator depends on the wavenumber i . The full diagonalization method consists of a two step method. Consequently, in the first step named preprocessing stage, the eigenvalues ( ) and eigenvectors (P) of the operators are computed such as D 2 ir = P -1 ir ir P ir and D 2 z = P -1 z z P z . This expensive step is performed only once during the preprocessing. During the processing stage the solution is computed at each time step. It involves fast Fourier transforms and tensorial/pointwise products for each wavenumber of the Fourier transform as follows:

1. compute the forward FFT of the source term: ˆ S = F θ (S) .

perform tensorial product

F i = P -1 ir P -1 z ˆ S i 3. compute pointwise product F i jk = F i jk / λ ir, j + λ z,k 4. perform tensorial product F i = ( P ir P z ) F i 5. compute the backward FFT to get the solution φ = F -1 θ F
The second step involves 2 Fast Fourier Transforms and 4 tensorial products. The algorithm complexity of this linear solver is

O (n 2 r n z + n r n 2 z + n r n z log n θ ) n θ .
Thus, the Fourier series expansion in the azimuthal direction contributes to reduce computational resources, otherwise a full 3D solver has a higher linear complexity O n θ n 2 r n z + n θ n r n 2 z + n r n z n 2 θ .

Code parallelization

The parallelization of the present code relies on two distinct strategies: to evaluate the compact derivative/interpolation and to get solutions of Poisson and equation problem [START_REF] Abide | An efficient parallel high-order compact scheme for the 3d incompressible Navier-Stokes equations[END_REF] . Both approaches are described in the following sections.

Helmholtz/Poisson solvers

First, the parallelization of the full diagonalization method still remains a tricky task. Indeed, the core of the method is mainly based on tensorial products to cross between the physical and eigenvectors spaces, and the lack of sparsity in the eigenvector matrices involves expensive communications. Here, the 2d-pencil decomposition is considered to get a parallel version of the full diagonalization method [START_REF] Abide | An efficient parallel high-order compact scheme for the 3d incompressible Navier-Stokes equations[END_REF] . The guideline of the 2d-pencil decomposition is to provide a 2D block decomposition where processors are filled into a 2D grid. Two of the dimensions of 3D data arrays are distributed on the processor grid. Such a decomposition leads to entirely maintain one dimension of the grid within the local memory. On this former direction, FFT or tensorial products can be performed in a serial way. An example of such an implementation is provided by the libraries P3DFFT or 2DECOMPFFT [START_REF] Pekurovsky | P3DFFT: a framework for parallel computations of fourier transforms in three dimensions[END_REF][START_REF] Li | 2DECOMP&FFT-a highly scalable 2d decomposition library and FFT interface[END_REF] . Fig. 2 illustrates the 2d-pencil decomposition in the cylindrical configuration, and the three pencil states.

Global data transposes are performed to switch between the three states inducing expensive MPI-all-to-all communications. The parallel algorithm of the diagonalization method can be summarized as follows:

1. compute ˆ S = F θ (S) on the θ -pencil state, 2. transpose ˆ S from θ to r -pencil then compute

F i = P -1 ir I z ˆ S i
3. transpose ˆ F from r to z -pencil then compute G i = I r P -1 z F i 4. compute F i jk = G i jk / λ ir, j + λ z,k 5. compute G i = ( I r P z ) F i then transpose G from z to r -pencil, 6. compute F i = ( P ir I z ) G i then transpose F from r to θ -pencil, 7. compute the final solution using a backward FTT: φ =

F -1 θ (F )
The previous algorithm involves 4 all-to-all communications during steps (2-3) while [START_REF] Karniadakis | Spectral/hp element methods for computational fluid dynamics[END_REF][START_REF] Cantwell | Nektar ++ : an open-source spectral/hp element framework[END_REF] are provided by a 2d-pencil decomposition library. Due to the implicit nature of the boundary conditions [START_REF] Abide | A 2d compact fourth-order projection decomposition method[END_REF] two additional transpositions are required. Therefore, the total number of all-to-all communications necessary to solve each Helmholtz/Poisson problem is 6. Although the all-to-all communications are expensive, this approach exhibits a correct strong scaling [START_REF] Abide | An efficient parallel high-order compact scheme for the 3d incompressible Navier-Stokes equations[END_REF] and contributes to drastically reduce the real execution time of a simulation. However, it should be stressed that a poor weak scaling is expected due to the O(n 4 ) algorithmic complexity of the diagonalization method. This is inconsistent with the concept of weak scaling, unlike methods with O(n ) algorithmic complexity as the multigrid methods.

Derivative and interpolation computations

The spectral evaluation of derivatives and interpolations in the periodic direction θ involves real-to-complex and complex-to-real FFT. Computations are made in a serial manner since the default pencil-state is the θ -pencil ( Fig. 2 ). The reduced Partial Diagonal Dominant (rPDD) algorithm of Sun [START_REF] Sun | Application and accuracy of the parallel diagonal dominant algorithm[END_REF] underlies the computation of the compact derivatives and interpolations in the present approach. This method of solving tridiagonal linear systems shows interesting parallel performances [START_REF] Abide | An efficient parallel high-order compact scheme for the 3d incompressible Navier-Stokes equations[END_REF][START_REF] Sun | Application and accuracy of the parallel diagonal dominant algorithm[END_REF] . The diagonal dominance is exploited to derive an approximate solver involving only neighborto-neighbor communications, greatly improving the parallel efficiency. Specifically, 2 × 2 linear systems involving interface nodes (reduced system) is derived by dropping terms through the diagonal dominance property [START_REF] Sun | Application and accuracy of the parallel diagonal dominant algorithm[END_REF] . The analysis and the reliability of the rPPD algorithm in the specific framework of the incompressible Navier-Stokes equations in cartesian coordinates is detailed by [START_REF] Abide | An efficient parallel high-order compact scheme for the 3d incompressible Navier-Stokes equations[END_REF] : a description is briefly recalled hereafter. The calculation of a compact scheme derivative like Eqs. ( 9) and ( 10) using rPDD , consists of a three-step algorithm involving mainly the solution of the distributed tridiagonal linear systems M p on P processors. These steps are summarized hereafter:

1. compute an approximate compact scheme derivative ˜ u p = M -1 p B p u p , p standing for a subdomain partition, 2. exchange boundary informations (neighbor-to-neighbor communications) and compute solutions of the reduced systems, 3. update the solution ˜ u p from solutions of the reduced systems independently on each processor.

The accuracy of this approximation depends on the diagonal dominance of the linear system [START_REF] Sun | Application and accuracy of the parallel diagonal dominant algorithm[END_REF] . For the fourth-order compact schemes ( Table 2 ), it is shown that using about 16 nodes per subdomain does not bring significant deterioration of either the accuracy or conservation properties. A good strong scaling has been observed for up to 10 4 processors [START_REF] Abide | An efficient parallel high-order compact scheme for the 3d incompressible Navier-Stokes equations[END_REF] . It is worth mentioning that the staggered grid layout involves about 20 nodes for the convective terms which would require naively about 80 all-to-all mpicommunications with the pencil-decomposition strategy. In conclusion, the main advantage of the reduced PDD algorithm is to address a less intensive parallel communication than the pencil decomposition.

Results and discussions

This section investigates the effects of the rPDD algorithm within the context of rotating flows, and estimates its parallel performances. First, the accuracy of the temporal and spatial schemes is demonstrated over an analytical solution of the Navier-Stokes equations. Then, results for the two configurations described in Section 2.1 are presented to assess the rPDD algorithm performance and its parallel efficiency.

Numerical accuracy

Both spatial and temporal accuracies are determined with the Method of Manufactured Solutions (MMS) [START_REF] Roache | Code verification by the method of manufactured solutions[END_REF] . A fluid of viscosity ν = 0 . 1 is considered in an annular cavity defined by η

= R i /R o = 0 . 5 , d = R o -R i = 1 and = H/d = 2 .
According to the MMS, a source term is derived from an analytical divergence-free velocity and associated pressure fields satisfying the Navier-Stokes equations. Here, the following divergence-free velocity and pressure fields have been considered:

u r = + 1 2 π sin 2 ( π ( r -r i ) /d ) cos ( θ ) sin ( 2 π z/H ) ( 1 + β cos ( 2 πt ) ) ( 13
)
u θ = - 1 2 π sin 2 ( π ( r -r i ) /d ) sin ( θ ) sin ( 2 π z/H ) ( 1 + β cos ( 2 πt ) ) ( 14
)
u z = - H 2 π d sin 2 ( 2 π ( r -r i ) /d ) cos ( θ ) sin ( π z/H ) ( 1 + β cos ( 2 πt ) ) (15) p = ( sin ( π ( r -r i ) /d ) + sin ( π z/H ) ) cos ( θ ) ( 1 + β cos ( 2 πt ) ) ( 16
)
The steady solution is obtained for β = 0 , whereas β = 1 leads to unsteady solution. The discussion on the numerical accuracy is based on the maximum relative error:

φ (t) = φ num -φ ex ∞ φ ex ∞ ( 17
)
where φ num and φ ex stand respectively for the numerical and analytical fields of the velocity components and the pressure. First, the space accuracy is estimated from error (17) using various grids with increasing sizes. Tests are performed with the steady solution ( β = 0 ), and the numerical error for each variable is reported in Fig. 3 a.

The numerical error decreases with a -4 slope corresponding to the expected accuracy of the present compact scheme discretization.

The temporal accuracy is then evaluated by considering the unsteady solution ( β = 1 ) leading to a time-dependent maximum relative error. The time-accuracy analysis is performed on the maximum value of this error during several periods. Fig. 4 shows the history of the numerical error and its maximum values after several periods versus the time step. Fig. 4 exhibits a -2 slope for each variable, in agreement with the expected order 2 of accuracy. These findings confirm that the accuracy requirements are well fulfilled by the present code. The influence of the rPDD algorithm over the computed solutions is now investigated. For this purpose, node-wise comparisons between simulations using the rPDD and the classical ones are performed. It is recalled in Section 3 that the accuracy of the rPDD algorithm depends on the number of nodes involved in a partition. This parameter is denoted by n / p . Thus, the deviation introduced by the rPDD on the numerical solution is measured as:

φ (n/p) = φ rpdd -φ std ∞ φ std ∞ ( 18
)
where std refers to the standard evaluation of compact scheme, or without error. In the following tests, the size number of the partition is kept identical in each non-periodic directions ( n/p = n r /p r = n z /p z ). The maximum difference φ versus the partition size n / p is plotted in Fig. 3 b. It is observed that the maximum difference drastically decreases with respect to the parameter n / p , in accordance with previous analysis [START_REF] Abide | An efficient parallel high-order compact scheme for the 3d incompressible Navier-Stokes equations[END_REF][START_REF] Sun | Application and accuracy of the parallel diagonal dominant algorithm[END_REF] . This behaviour indicates that above a partition size n/p = 15 no significant difference is observed, and consequently that the accuracy is not deteriorated.

Taylor-Couette flow with an axial temperature gradient

Increased attention has recently been paid to the Taylor-Couette flow with a vertical stably stratified temperature field. Indeed, despite the stabilizing differential rotation and thermal stratification, the flow may become unstable at specific values of the Reynolds and Froude numbers. This instability, named Stratorotational Instability (SRI), is suspected to be responsible for the outward angular momentum flux which is a necessary condition for materials to be transported inwards towards the central body in accretion disks. The major difficulty encountered in the numerical study of this configuration stems from the very long transient regime before the occurrence of the instability. As a consequence, simulation times to reach converged solutions are extremely lengthy, leading to huge computational times when using a sequential code. Moreover, the large aspect ratio in the axial direction, typical in this specific configuration, requires consequent number of modes, further increasing the simulation times. Thus, the present algorithm is a promising strategy to strongly reduce wall-clock times, without nevertheless reducing the algorithm complexity of the diagonalization method but excluding ideal weak scaling. For instance, the Fig. 5 illustrates the occurrence of a Strato-rotational instability by probing the velocity components at mid-height and radius of the cavity.

The converged solution requires about 30 0 0 0 0 temporal iterations. In this section, several simulations are carried out to check the effects of the rPDD algorithm on the variation of the angular momentum transport.

We consider the flow in a tall annular cavity with η

= R i /R o = 0 . 5 , d = R o -R i = 1 and = H/ (R o -R i ) = 6 .
The flow and heat transfer are governed by the dimensionless parameters Re, Fr and μ (see Section 2.1 ). The angular momentum is related to the normalized Reynolds stress, represented by the β viscosity [START_REF] Gellert | Stratorotational instability in Taylor-Couette flow heated from above[END_REF] :

β = u r u θ R 2 i 2 i ( 19
)
where fluctuating velocities u r and u θ correspond to deviations from the mean value in the azimuthal direction.

First, the β-viscosity is computed for Re = 40 0 , 60 0, 70 0 and 800 with fixed values of parameters (P r, F r, μ) = (1 , 1 . 4 , 0 . 354) . The mean β-viscosity along the z -axis is considered for comparison with the simulations of Gellert and Rüdiger [START_REF] Gellert | Stratorotational instability in Taylor-Couette flow heated from above[END_REF] at mid-gap r = 1 . 5 d. Computations were carried out on a grid size 48 × 96 × 96 in ( r, θ , z ), found to be sufficient for describing the angular momentum. Fig. 6 a shows the Reynolds stress pattern in a rz plane for Re = 600 . The evolution of the β-viscosity versus the Reynolds number is displayed in Fig. 6 b. Despite a discrepancy from the profiles reported by Gellert and Rüdiger [START_REF] Gellert | Stratorotational instability in Taylor-Couette flow heated from above[END_REF] , which can be explained by the very different numerical techniques used including grid resolution, the linear dependence on the Reynolds number is recovered for Re ≥ 600.

To investigate the effects of the approximation introduced by the rPDD algorithm (section 3.4.2), the β-viscosity is computed for several node densities n / p at Reynolds number Re = 600 . Table 3 gathers the differences between the solutions obtained with and without the rPDD algorithm. The rPDD algorithm is used for different partition size ranging from 48 nodes to 6 nodes per subdomain. The simulation performed with n/p = 96 corresponds to the standard evaluation of the compact scheme derivatives. Table 3 reveals noticeable discrepancies on the velocity field and divergence below n/p = 12 nodes per subdomain. A significant difference is obtained for 6 nodes per subdomain, where a difference of about 10 -4 is observed. It should be mentioned that the divergence of the flow is around 5 10 -12 , and according to our experience this value corresponds to the threshold of the round-off error of the present simulation.

Uniformly rotating flow with a radial temperature gradient: The baroclinic cavity

The aim of this section is to assess the parallel performance of the present solver in the baroclinic cavity, in particular its ability to capture the different scales reported by Von Larcher et al. [START_REF] Larcher | Instabilities and small-scale waves within the Stewartson layers of a thermally driven rotating annulus[END_REF] within the same configuration. For this purpose, the flow in a differentially heated rotating annulus is considered, and more specifically the regime involving baroclinic instability. This flow configuration was extensively studied in the literature. Among others we can refer to the experimental work of Harlander et al. [START_REF] Harlander | PIV-and LDV-measurements of baroclinic wave interactions in a thermally driven rotating annulus[END_REF] , the numerical simulations of von Larcher and Dörnbrack [START_REF] Von Larcher | Numerical simulations of baroclinic driven flows in a thermally driven rotating annulus using the immersed boundary method[END_REF] and the recent benchmark of Vincze et al. [START_REF] Vincze | Benchmarking in a rotating annulus: a comparative experimental and numerical study of baroclinic wave dynamics[END_REF] . The test case to validate the code is for a 6 rpm rotation rate and with a radial temperature difference of T = 8 K. The governing dimensionless parameters are the Taylor and the Rossby numbers defined by Eq. ( 3) . In the following validation tests T a = 1 . 69 × 10 6 and Ro = 0 . 2 . It is worth mentioning that the initial state to start the present simulation is an under-resolved solution (64 × 64 × 64 in radial, azimuthal and axial directions respectively) which was interpolated ( t = 0 ) to a finer grid 256 × 128 × 150. Both plots of Fig. 7 qualitatively show good agreement with previous results [START_REF] Vincze | Benchmarking in a rotating annulus: a comparative experimental and numerical study of baroclinic wave dynamics[END_REF] .

The Fig. 8 shows instantaneous isocontours in ( θ , z ) plane of the horizontal divergence of the velocity ∇ h . u ≡ -∂w ∂z along the inner cold wall ( Fig. 8 a) and along the outer hot cylinder ( Fig. 8 b). This variable is introduced to exhibit the occurrence of small-scale features simultaneously with the large-scale baroclinic waves. In the first plot, the small-scale structures, developing towards the bottom wall, have been identified as inertia gravity waves (IGWs) by different authors [START_REF] Borchert | Gravity wave emission in an atmospherelike configuration of the differentially heated rotating annulus experiment[END_REF][START_REF] Hien | Spontaneous inertia-gravity wave emission in the differentially heated rotating annulus experiment[END_REF][START_REF] Larcher | Instabilities and small-scale waves within the Stewartson layers of a thermally driven rotating annulus[END_REF] in similar water-filled cavities. The present observed features recall such IGWs reported by these authors. Recently, for the same present configuration, Von Larcher et al. [START_REF] Larcher | Instabilities and small-scale waves within the Stewartson layers of a thermally driven rotating annulus[END_REF] mentioned the presence of ripples resulting from hydrodynamical instability along the outer hot cylinder. We also capture the same phenomenon, as illustrated by the isocontours of the horizontal divergence.

A quantitative validation is obtained from the comparison of the values of drift rate of the large-scale baroclinic waves. The drift rate c m ( t ) of a wave mode m is calculated from the phase shift φ m ( t ) at each time t by : 0 . 0115 rad.s -1 , which favorably compares with the values previously obtained experimentally and numerically (see Fig. [START_REF] Abide | An efficient parallel high-order compact scheme for the 3d incompressible Navier-Stokes equations[END_REF] in [START_REF] Vincze | Benchmarking in a rotating annulus: a comparative experimental and numerical study of baroclinic wave dynamics[END_REF] ). The code being validated on the prediction of the baroclinic instability and associated small-scale features, we now focus on the parallel performance of the proposed method.

c m (t) = 1 m d φ m d t (20)
A strong scaling test is performed to assess the computational performance. It consists in measuring the computational time versus the number of cores. Three meshes are considered: 128 × 128 × 128, 128 × 264 × 264 and 128 × 568 × 568. The mesh in the r and z-directions is distributed on a square process grid of size P 2 ranging from 4 to 1024 cores. The computational costs to calculate the convective terms, t conv , to solve a Poisson equation, t diag , and to simulate one time step, t it , are measured. These values are presented in order to compare the performances of the two parallel strategies: the rPDD (explicit terms) and the 2d-pencil decomposition (implicit terms). The tests were performed on the CALMIPs EOS cluster facilities which consists of 612 computation nodes of 2 Intel(r) IVYBRIDGE 2,8 GHz 10-cores. Fig. 10 shows the computational times versus the number of used cores.

For each mesh size, a similar behaviour is observed. The computational cost decreases with respect to the core number ( Fig. 10 a). For small core numbers, the scalability is satisfactory, but for the largest ones, a loss of scalability is observed with the coarser meshes. It is reasonable to speculate that this loss of performance will also occur with finer mesh for larger core numbers. This behaviour can be attributed to the ratio between the volume of computations per core and the data exchange involved in the MPI communications. global communication has to be tuned according to the cluster to get better efficiency. Fig. 10 c shows the strong scaling test during the computation of convective terms. This computation relies only on the calculation of the compact interpolations and derivatives using the rPDD approach. In this case, a good scaling is noted for the two finest grids up to 10 3 cores. In comparison with the 2 dpencil decomposition, the rPDD involves only neighborhood communications which are less expensive. However, with the coarsest grid, the scalability from 10 2 cores is deteriorated, which can be also explained by the unbalance between volume of computations and communications. Finally, it must be noted that the scalability and the computational time are dominated by the solution of Poisson/Helmholtz equations as shown by the similarities of the curves in Fig. (10 a) and ( 10 b). It should be noted that better performances are observed in Fig. 10 b with the largest grid, wich may be explained by a memory cache effect. Even if a specific tuning can further improve performances, the present method exhibits a wall time reduction from 25 s with 4 cores to 0.35 s with 1024 cores when using a standard mesh of 128 × 264 × 264.

Conclusion

In this manuscript, we have extended a cartesian parallel higher-order compact scheme solver to cylindrical coordinates. The azimuthal direction has been discretized using Fourier series expansion to benefit from the natural periodicity, and the favorable complexity algorithm of the Fast Fourier Transform. The space discretization in the two wall-bounded directions relies on the fourthorder compact scheme approximations. The guideline concerns the simulation of two specific configurations of rotating stratified flows: the Taylor-Couette setup under an axial thermal stratification and the baroclinic cavity. The code parallelization strategy combines two approaches. The first one is the 2 d -pencil decomposition to address the parallel solution of the implicit viscous terms and the pressure-like equations based on the diagonalization method. The second strategy of parallelization consists in the calculation of the compact derivatives/interpolations. The approximate tridiagonal solver named reduced Partial Diagonal Dominant (rPDD) algorithm is used to evaluate compact derivatives and interpolations. According to our knowledge, the proposed parallelization strategy is a totally new method.

The present code is successfully validated against the Method of Manufactured Solution, and in the two test configurations. The effects of the approximation introduced by the rPDD on the solution accuracy is evaluated. Precisely, we show that with partition domains beyond 16 nodes the rPDD does not introduce noticeable error. Finally, a strong scaling test is performed to outline the potentiality of the present parallel strategy to simulate "specific rotating stratified flows". In particular, the present tool has allowed to drastically reduce the execution time for the two configurations, in comparison with sequential approach. Moreover, the proposed technique is able to correctly capture all the characteristic solution features : the SRI, the baroclinic instability and the associated small-scale structures (IGWs and outer ripples). Further exploration of these complex phenomena is now considered. It particularly concerns DNS of the turbulence in SRI. For the baroclinic configuration, the next studies should be dedicated to the large-tank facilities available at BTU, Cottbus-Senftenberg Germany (team of Professor U. Harlander). This setup allows for investigating IGWs in "atmosphere-like" situation with rotation rate smaller than thermal stratification. Nevertheless, the scalability test puts forward the pitfall introduced by the large communication scheme in the 2d-pencil decomposition. Performances could benefit from the hybrid OpenMP/MPI parallelization to enhance the computation/communication balance. Recherche". The authors acknowledge financial support from "Ministère des Affaires Etrangères", France, in the framework of PRO-COPE projects ( PHC 28369XH and PHC 35299PM ). The authors are grateful to Pr U. Harlander (BTU, Cottbus-Senftenberg, Germany) for fruitful discussions.

Fig. 1 .

 1 Fig. 1. Rotating annular cavities configurations.

Fig. 2 .

 2 Fig. 2. The 3 pencil decomposition states defined on 3 × 4 processor grids.

Fig. 3 .

 3 Fig. 3. Accuracy assessment of the space discretization scheme (a) the numerical errors versus size mesh, (b) difference error versus the partition size.

Fig. 4 .

 4 Fig. 4. Accuracy assessment of the temporal scheme (a) time history of the numerical errors, (b) maximum of the numerical error versus time step.

Fig. 5 .

 5 Fig. 5. Occurrence of the SRI: time history of the sampled velocity at mid-cavity for (Re, Pr, F r, μ) = (600 , 1 , 1 . 4 , 0 . 354) (skipped symbol 20).

Fig. 6 .

 6 Fig. 6. The β-viscosity for F r = 1 . 4 , η = 0 . 5 and μ = 0 . 354 (a) Reynolds stress pattern for Re = 600 , (b) Variation of β-viscosity with respect to the Reynolds number Re at mid-gap.

Fig. 7 a

 7 displays a snapshot of the free surface temperature showing the occurrence of the dominant azimuthal mode m = 3 on the baroclinic instability. This wave-3 structure is confirmed by the corresponding Hövmoller azimuth-time diagram at mid-radius and height, presented in Fig. 7 b.

Fig. 9

 9 Fig. 9 shows the time history of the phase shift divided by the azimuthal wavenumber m . The uniform drift for the dominant mode m = 3 and the associated slave mode m = 6 can thus be obtained. The predicted drift associated with the dominant mode m = 3 is c 3 (t) =

Fig. 7 .

 7 Fig. 7. Temperature at upper free-surface: (a) Snapshot (b) Hövmoller diagram at mid-radius and height in azimuth-time.

Fig. 8 .

 8 Fig. 8. Instantaneous isocontours in ( θ, z ) plane of the horizontal divergence of the velocity: at radius close to inner wall (a), at radius close to outer wall (b).

  Fig. 10 b concerns the strong scaling test for the solution of a Poisson equation based on the 2 d -pencil decomposition. The trend observed in Fig. 10 a still holds. This can be explained by the large communication volume in the switch between two states of the 2d-pencil decomposition [18] . Indeed, a

Fig. 9 .

 9 Fig. 9. Time history of the azimuthal distances m / m of wave modes m = 3 and m = 6 .

Fig. 10 .

 10 Fig. 10. Strong scaling test (a) temporal iteration (b) solution of Poisson equation (c) convective terms.

Table 2

 2 Compact scheme coefficients.

		Inner nodes	Boundary nodes
	Operators	α	a	α	a	b	c
	D c f ξ D f c ξ I c f ξ I f c ξ	1/22 1/22 1/6 1/6	12/11 12/11 2/3 2/3	-1 23 5 1	-1 -25 15/4 1/4	2 26 5/2 3/2	-1 -1 -1 / 4 1/4

Table 3

 3 Effects of the approximate rPDD algorithm at Re = 600 for different node densities n / p ( n/p = 96 corresponds to the simulation without using PDD algorithm and x(m) denotes x × 10 m ).

	n / p	u r	u φ	u z	T	∇ • u	β
	96	-	-	-	-	-	9.72( -12)
	48	3.60( -11) 3.20( -12) 2.52( -11) 3.14( -13)	5.51( -13) 1.19( -11)
	24	3.51( -11)	2.95( -12) 2.63( -11) 3.20( -13) 5.20( -13) 1.00( -11)
	12	3.13( -09)	3.34( -10) 1.82( -09) 8.09( -11) 1.72( -10)	1.00( -11)
	6	2.98( -04) 1.23( -04) 2.67( -04) 1.75( -05) 3.22( -05) 3.66( -06)

Acknowledgements

This work was granted access to the HPC resources of CALMIP supercomputing center under the allocation 2017P17027 and to the HPC resources of Aix-Marseille Université financed by the project Equip@Meso (ANR-10-EQPX-29-01) of the program "Investissements d'Avenir" supervised by the "Agence Nationale de la


Shan Zhao

Guifen Chen

Siwei

Enze Xiao Fu

Clc

Study on Precision Fertilization Model Based on Fusion Algorithm of Cluster and RBF Neural Network

Keywords: fusion optimization algorithm, FCM, OLS, RBF, precision fertili zation model

Precision fertilization is the core content of precision agriculture technology, There was a complex non-linear relationship between crop optimal fertilization and soil fertility, The single model is difficult to accurately describe its complex relationship and change law, making the crop accurate fertilizer is difficult to determine; Neural network technology to solve this problem provides a new way of thinking, But a single radial basis function network (RBF) neural network fertilization model is too dependent on the selection of the hidden layer data center, Therefore, this paper proposes a decision-making technique based on fuzzy C-means (FCM) clustering and RBF neural network fusion algorithm. The fusion algorithm first uses the FCM algorithm to select multiple RBF networks in the training samples. Based on this, the least squares (OLS) training network is used to optimize the data center, Finally, an improved RBF neural network model is established. In this paper, the model is applied to the maize precision operation demonstration base, soil nutrient and maize yield as the the input of neural network, using the precision fertilization amount of maize as output, The model of precision fertilization of maize was established, And the model was used to make the precision fertilization decision of maize. Experimental results show: The improved RBF neural network is compared with the traditional BP network to reduce the error by 0.47, Compared with the model, the error of the RBF neural network method is reduced by 0.045, Significantly improve the prediction accuracy, reduce the calculation time, Can effectively guide the precise fertilization.

Introduction

Maize production in Jilin Province, agricultural production occupies an extremely important position, is the largest planting area in Jilin Province, food crops, production accounted for 12.5% of the country, to 27.88 million t. corn production of Jilin Province is essential to China's food security. Reasonable application of chemical fertilizers is one of the important guarantees of high and stable yield of maize.

According to the soil type, soil fertility and maize fertilizer characteristics to fertilize scientifically and reasonably is particularly important. Variable fertilization technology is an important part of precision agriculture, It is based on the actual needs of crops, based on scientific fertilization methods, Identify variables inputs on crop To a certain extent, increased corn production. In the case of reasonable control of the amount of fertilizer, you can achieve both the purpose of reducing investment and increase production.

Because of the high nonlinearity between yield and soil nutrient and fertilization, neural network is considered to be a powerful tool to solve nonlinear problems, Therefore, a number of different forms of Guifen Chen(1956-), College of Information Technology, Jilin Agricultural University,ChangChun,China network in terms of approaching ability and learning speed, and can simplify the structure and learning rules of the network. Therefore, this paper uses FCM algorithm and OLS fusion algorithm to optimize RBF neural network, and constructs a RBF neural network prediction model based on fusion algorithm, To adapt the content of soil nutrients, fertilizer and food production randomness, complexity and nonlinear variation, Through the different soil types and yield prediction fertilization, to explore the realization of high yield of corn nitrogen, phosphorus, potassium fertilizer the best amount of application, For the rational use of fertilizer on the production, improve fertilizer utilization to provide a reference. To guide the variable fertilization in production of maize in central Jilin Province, in order to improve food productivity, to achieve sustainable development of productivity.

Technology and methods

Techniques employed herein route is divided into three parts: first, Normalizing treatment the Data samples, and using the FCM algorithm select primary data center; the second usage OLS to train RBF neural network to generate a regression matrix, preferably the data center; Thirdly, using the RBF neural network to construct the precision fertilization model. The specific description shown in Figure 1:

RBF neural network prediction algorithm

The RBF (Radial Basis Function) neural network model was proposed by Moody and Darken in 1988, it is a three-layer feedforward neural network, the first layer is the input layer, composed of the source node; the second layer for the hidden layer, the activation function using radial basis function (this article is Gaussian function),and Usually defined as any data sample to a data center Euclidean distance between the monotonic function; The third layer is the output layer, the hidden layer to the output layer using linear activation function. In this paper, four factors affecting the amount of maize fertilizer were obtained, so the input vector n = 4, h hidden layer, m output vector is the amount of soil fertilization of Alkaline N (N) Available phosphorus (P) Available potassium (K).so m = 3 By RBFNN learning process and its network mapping relationship, RBF performance depends on the strong choice of the center number, the center value and its location. Therefore, if all the input as a candidate space, it will bring a very large amount of calculation, reducing network performance. In order to overcome this deficiency, FCM algorithm can be used to narrow the candidate space and reduce the network size to improve the computational efficiency.

FCM clustering algorithm primary network center

     c ij n j u 1 i ,..., 1 , 1 (1)
The objective function of FCM algorithm is generally as follows: (2):

       c i n j ij k ij c i i c d u J c c U J 1 2 1 1 ) ,..., , ( (2)
In the formula, uij Is the degree of membership; dij=‖ci-xj‖Represents the Euclidean distance between the i-th cluster center and the j-th data point; ci is the clustering center of each fuzzy group; k∈[1,∞)Is a weighted power exponent. In order to achieve the minimum target function, We introduce the Lagrangian multiplier λj(j= [1,n]) The constructor is shown in Equation 3:

                n j c i n j n i c i ij j ij k ij c i ij j c n c i u d u u c c U J c c U J 1 1 1 1 2 1 1 ) 1 ( ) 1 ( ) ,..., , ( ) , ,..., , (    (3)
Derivating the parameters for each parameter, the requirement of Making Formula 3 Minimal is

     n i k ij n j j k ij i u x u 1 1 c (4) And C1 C2 C3 Σ Σ ... ... ...     c k k kj ij d d 1 ) 1 ( 2 ij ) ( 1 u (5)
The specific FCM algorithm steps are as follows:

Step 1: Given the clustering category c, 2≤c≤N, N is the number of initial data samples, Initial membership matrix is U(0), And set the thresholdε;

Step 2: initialize the cluster center ci (t), i = 1,2, ..., n;

Step 2: calculate the membership degree uij (t) with the formula (5), i = 1,2, ..., c; j = 1,2, ..., n;

Step 4: Use the obtained membership degree to correct the cluster center ci (t + 1), i = 1,2, ..., c according to (4);

Step 5: Calculate the error e=

    c 1 i 2 ) ( ) 1 ( ‖ ‖ t c t c i i
If e <Emax the operation is over, otherwise let t + 1 = t, return to the third step.

Orthogonal least squares optimization center

RBF network center selection orthogonal least squares (orthogoral least squares algorithm OLS) algorithm, Which solves the key problems in the construction of radial basis function network. OLS algorithm uses the orthogonalization method to independently calculate the contribution of the regression operator to the output, so the center selection step is simple and effective. OLS algorithm is the basic idea:

The center of the radial basis function is chosen a subset of the training pattern, Choose one sample at a time, Through the orthogonal regression matrix, the regression operator with large error compression ratio is selected, And the number of regression operators is determined by the selected tolerance, And then find the network weight. The RBF neural network is trained with OLS, and the matrix representation of the expected output response is:

E W E Y Y       (6)
Where Y∈R n * 3 is the desired output matrix; Y ̂∈R n * 3 is the output matrix of the neural network model; Ф∈R n * h is the regression matrix; E∈R n * 3 is the output error matrix. In this paper, Schmidt orthogonalization method to reduce the output error. OLS algorithm can save storage space and improve operation efficiency, and is very suitable for training RBF network. The RBF neural network trained by OLS can obtain a more optimized data center, improve the accuracy of the network, greatly reduce the computational complexity of the network and improve the training speed.

Application of RBF Neural Network Fusion Algorithm in Precision Fertilization Decision of Maize

Experimental data

Experimental data from the "National Spark Program" in Jilin Province, Nong'an County, Helong Town, Chen Jiadian Village and Kai'an town experimental field. In order to achieve accurate fertilization, the grid was divided into grid, the grid size is 40m × 40m.Figure 2

Pretreatment of training samples

In the process of network model training will appear a larger input variable to suppress the phenomenon of smaller input variables, resulting in increased network training time, and even the network can not converge. If the data dimension is inconsistent, it is easy to cause the sickness of the model. In order to remove the data in the dimension and the order of magnitude to establish the impact of network accuracy, After obtaining the sample vector, this paper takes the normalization of the data samples to deal with, namely:

} min{ } max{ } min{ ' i i i pi pi x x x x x    (7)
In the formula: Xpi is the raw data of the i-th variable of the pth sample, Xpi'is the normalized processing data for the i-th variable of the pth sample. The preprocessed sample data range is (0,1).

Construction of fertilization model

In order to achieve the target yield of corn, this paper establishes a RBF neural network model based on fusion algorithm, Study the best fertilization model of maize to achieve high yield of maize. Taking the soil nutrient content of Chenjiadian experiment field and the input as the neural network, the amount of fertilizer applied as the output of the neural network, In this paper, based on the RBF neural network fertilization model shown in Figure 4: 3: It can be seen from Table 3 that from the prediction accuracy, the FCM-OLS-RBF model constructed in this paper has a prediction error of 0.2, which is obviously better than RBF network and BP network.

From the processing speed, FCM-OLS-RBF network is 1.93,1.73 times faster than single RBF network and BP neural network in predicting 10 samples. Thus, the RBF network after the improved algorithm has higher accuracy and faster speed than the single RBF network and the traditional BP network, and the practicability is stronger.

( It can be seen from Figure 5 that the predictive curve of FCM-OLS-RBF network is in good agreement with the measured curve, especially the predicted value and measured value of N and P, and the deviation is very small, Showing its higher prediction accuracy.

Figure 6: FCM-OLS-RBF network prediction error

It can be seen from Fig. 6 that the relative error of nitrogen, phosphorus and potassium fertilizers is about 0.1%, and the maximum error is below 0.25%, which can meet the requirements of actual production.

Conclusion

This paper focuses on the RBF neural network algorithm based on FCM clustering and OLS, and compares it with the traditional RBF and BP neural network. Finally, the method is applied to the precision fertilization problem. According to the experimental test, the following conclusions are drawn:

(1) FCM-OLS-RBF fertilization prediction model training speed, good convergence, In the case of large quantities of non-linear and accurate fertilization applications, Compared with the traditional RBF and BP network prediction results, it also shows its better evaluation accuracy. Both to ensure the accuracy of the network while simplifying the algorithm and network structure;

(2) The fertilization model based on FCM-OLS-RBF neural network can simulate the nonlinear relationship between soil test and fertilization, which is superior to the traditional fertilization model and can be used as a quantitative model to guide the accurate fertilization;

(3) A large number of training samples are needed to predict with a neural network, or the accuracy of prediction can not be accurately fitted due to the high degree of non-linear relationship between the factors and fertilization.

And because the factors affecting the amount of corn fertilization varied and intricate, the factors that affect the yield is not only the amount of fertilizer problems. In the future, with the deepening of the practice of variable fertilization, the accumulation of data is rich, can provide a broader basis for improving production.

  neural networks are used for precise fertilization decisions. In 2001, Pokrajac et al. Presented a neural network-based decision support system to provide precise fertilization recommendations; In 2002, Teng Qingfang and others tried to use the neural network to establish the soil fertilization model; In 2004, Ma Weilin et al. Proposed a method of variable fertilization decision-making based on data envelopment analysis and artificial neural network. In 2007, Lan Weijuan et al. Applied radial basis function network to variable fertilization decision. The use of RBF network to study the variable fertilization decision has the advantages of simple structure, rapid training process and good promotion ability. It is superior to the BP ShanZhao (1995-),College of Information Technology, Jilin Agricultural University,ChangChun,China

Figure 1 :

 1 Figure 1: the technical roadmap for this article

Figure 2 :

 2 Figure 2:RBF network structure

  Fuzzy C-means (FCM) is proposed by Bezdek et al. FCM compared with the traditional clustering analysis algorithm,It introduced The fuzzy mathematics theory, which is an unsupervised clustering algorithm based on objective function. Corn soil nutrient, yield and fertilizer data samples usually have serious non-linear characteristics, Single model is difficult to describe, so it is necessary to divide the data sample set into different subclasses. The In order to reduce the number of RBFNN centers, FCM clustering is used to initialize the sample training set.And can improve the speed of network computing efficiency, so that the forecast is more accurate. The basic idea of FCM is to divide all data samples into c fuzzy groups, The cluster center of each fuzzy group is determined by minimizing the objective function value. And the degree of membership in the [0,1] is used to determine the degree to which each given data belongs to each fuzzy group. Its membership matrix U must satisfy equation (1):

  is Chenjiadian village and Kai'an town experimental field grid.

  (a)Chen Jiadian village plot grid (b)The kaian town grid diagram

Figure 3 :

 3 Figure 3: Experimental area gridIn the year of 2015, the experimental field was tested, fertilized and tested, and finally 64 samples were obtained. Each sample consisted of seven attributes: soil N, P, K content, actual yield, N, P, K fertilizer The amount of application, part of the experimental data shown in Table1 and Table 2:

Figure 4 :

 4 Figure 4: RBF neural network fertilization model (1) different algorithms to predict performance comparison In the training RBF network, the FCM algorithm, take c = 20, in all the samples, through clustering, the first selected cluster center number 20 data centers; use OLS training RBF, take the error threshold of 0.01 to get regression Matrix information, On the basis of the initial selection of 20 centers, and then enter the choice, and finally get 10 effective data center. Using BP network, single RBF network and post-training RBF neural network (FCM-OLS-RBF) to test the sample, Test results are shown in Table3:

  ) FCM-OLS-RBF fusion algorithm simulation experiment In this experiment, 64 samples of soil nutrient content, yield and actual fertilization amount in Chenjiadian and Kai'an were divided into three parts, Of which 44 samples were used as training set, 10 as validation sets, and the remaining 10 samples were used as test sets. Before training RBF network, we first use FCM algorithm to select 20 centers in all training samples, then train the network with OLS, and finally select 10 centers. Finally, the trained RBF neural network was used to predict the soil fertilization amount of maize in the field of maize, Simulation results are shown in Figure 5:

Figure 5 :

 5 Figure 5: FCM-OLS-RBF network simulation results

Table

Table 2 :Table 1 :

 21 1 and Chenjiadian Village Experimental Field

	Plots N	P	K	N	P	K	Measured
		conten	content	content	amount	amount	amount	yield
		tmg/kg	mg/kg	mg/kg	kg/hm	kg/hm	kg/hm	kg/hm
	a1	115	14.8	159	75.6	32.3	64.8	8303
	a2	115	14.8	159	75.6	32.3	66.1	8470
	a3	136	17.55	159	76.9	31.4	66.8	8556
	a4	136	17.55	159	76.9	31.4	66.1	8461
	a5	133	14.6	182	76.7	32.4	66.3	8684
	a6	133	14.6	182	76.7	32.4	56.2	7367
	a7	129	17.55	182	76.5	31.4	61.1	8002
	a8	129	17.55	182	76.5	31.4	65.4	8564
	a9	118	16.3	165	75.3	32.3	60.1	7745
	a10	125	15.8	179	77.1	32.3	62.7	8192

Table 2 :

 2 Kai'an Town Experimental Field

	Plots N	P	K	N	P	K	Measured
		conten	content	content	amount	amount	amount	yield
		tmg/kg	mg/kg	mg/kg	kg/hm	kg/hm	kg/hm	kg/hm
	a6	133	26.7	240	76.7	28.80	70.0	7696.5
	a7	108.5	38.6	270	75.2	26.20	61.6	7592.25
	a8	161	30.3	180	78.7	27.96	72.2	8180.25
	b5	108.5	35.8	230	75.2	26.76	65.7	7805.25
	b6	101.5	35.2	230	74.8	26.89	66.0	8817
	b7	140	30.7	160	77.2	27.87	73.4	7587.75
	b8	175	36.2	200	79.8	26.68	67.5	7221.75
	c4	98	21.4	220	74.6	30.22	74.9	6804
	c5	112	32.45	280	75.4	27.46	63.9	7360.5
	c6	168	33.2	240	79.2	27.31	66.3	7576.5

Table 3 :

 3 Comparison of different algorithms for prediction performance

		BP	RBF	FCM-OLS-RBF
	TIME(s)	2.008	2.234	1.158
	MSE	0.67115	0.243501	0.2

Acknowledgments

This work was funded by National Spark Program "Based on the Internet of Things precise operation of corn technology integration and demonstration" (2015GA660004) 2015-2017


Yinglun Li

Guifen Chen
email: guifchen@163.com

Dongxue Wang

Research and Application of 3D Visualization Plug-in Integration with ArcGIS

Keywords: ArcGIS, 3D visualization, auxiliary plug ins, farmland monitoring 1

With the development of agricultural production data, more and more three-dimensional visualization of monitoring data has become the focus of agricultural information technology research. The application of ArcGIS visualization technology, 3D visualization of auxiliary plug-in, and the national Spark Program demonstration area of farmland monitoring data obtained in Jilin province Nong'an County Helong town of application. In this study, the use of C++ language, the introduction of three-dimensional visualization library SDK, compiled auxiliary plug-in ArcGIS software, optimize the joint jiulongzhen farmland monitoring data 3D visualization. The experimental results show that the preparation of the postprocessing plug-ins integrated into ArcGis software can optimize the 3D visualization of monitoring data, the two-dimensional spatial variability of thematic map than the traditional representation of 3D data more natural, more clear, more intuitive, is of great significance for promoting the development of modern agriculture.

1. visual data is discrete data. In order to improve accuracy, difference method is needed.

2. visual data can be regular or irregular. Regular data can save storage space, and irregular data can give us greater freedom of expression.

3., the visual data must have an extension dimension, that is, the 5 data types mentioned earlier.

After identifying the issues related to data attributes, you need to prepare for dealing with large amounts of data. Data is produced at a much faster rate than it is managed and analyzed. The cardinality of the columns that need to be visualized is also a factor that should be taken into account. The high cardinality means that the column has a large number of different values, while the lower cardinality indicates that the column has a large number of duplicate values, such as [START_REF] Bo | Research on 3D graphic planning system of GIS and VR technology[END_REF] .

data preprocessing

The visual module inputs analog data and experimental / measurement data from the computer, and includes data calls in the database. The objects of visual processing include symbols, structures, images, TIN is a continuous blanket of topographic surfaces consisting of many adjacent but not overlapping triangles. Different types of TIN have different generation methods. According to the geometric structure of the spatial distribution, there can be a general triangulation and Delaunay triangulation.

In contrast to TIN, Grid is a composite surface composed of several consecutive squares of the same shape in space. The data structure of Grid is very regular and simple. It is possible to construct continuous, smooth and smooth surfaces through these data, and the structure of the mathematical model is simple.

Thus, Grid has the advantage of constructing a visual representation of good ground models [START_REF] Xie Zhongjun | A summary of three dimensional GIS aided agricultural planning[END_REF] .

Many scholars at home and abroad have made extensive research in the geological level generation methods, some scholars have used Grid triangulation into triangular data to generate geological level [START_REF] Penggen | Research and application of 3D geological model construction method[END_REF] .

The triangulation is the simplest for Grid, which splits into the triangular mesh structure and gets the raw data in a regular and systematic way. From the rules of the generation of TIN is mainly in the grid decomposition can be combined in the triangle, regular grid case, with one or two days to form diagonal grid decomposition of a series of rules of the triangle. After data processing, the interpolation of geological space is carried out, and the purpose is to grid the data [START_REF] Yingbiao | Three dimensional visualization management system of mine geological exploration and Research on the modeling technology Tian[END_REF] . Based on the grid data obtained, geological layers are plotted by using VC++ and OpenGL.

The main idea of the establishment of this method based on TIN multilayer geological body is, according to the geological layers with the elevation information to establish the TIN into a boundary between two layers of TIN, using the boundary and TIN were built erected under two TIN layer, the vertical TIN to simulate the geological body solid part.

(1) data arrangement;

The main tasks of data collation are drilling data extraction and data slicing. Thus, the TIN formed by the geological point is the surface, and the geological body with the thickness of the TIN as the lateral body is formed.

In support of the national Spark Program, the experiment and research in the field of monitoring data obtained in Nongan County, Jilin Province town of helong. The total area of the experimental plot is about 375 mu, and the grid size is set up [START_REF] Yingbiao | Bell ear research mine geological exploration and visualization management system modeling[END_REF] . There are several sampling points in the plot. We chose the moisture data for comprehensive analysis. After finishing, we drew a simple 3D block rendering by drawing the plots and the marked TIN points in ArcScene. As shown in figure 4. The first parameter is a selector whose value indicates the execution phase the plug-in is about to enter.

The second parameter is a callback function that points to the large data structure that contains the host application status information and hooks, which is a primary parameter.

The third parameter is a pointer to a structure that communicates its needs to the host as a plug-in. For example, error messages can be sent to the host via mail.

The structure of pf_layerdef output parameter point, which contains the output image information (including a pointer to the data): its size, its rowbytes value, and extent_hint, which is directly to the border area, really need to render [START_REF] Wang Minghua | Three-dimensional geological modeling of clouds present situation and development trend geotextile base[END_REF] .

Finally, the extra pointer is a special use parameter that is largely unused to the program. It points to different data structures at different times. If you decide to implement a custom user interface element, or a complex custom data type, you will only process this parameter [START_REF] Zhang Kai Election | and other three-dimensional geological ArcGIS-based visualization[END_REF] .

ArcGIS camera plug-in function

In the general case (without opening the 3 dimension layer attributes) camera has no effect, open the layer 3D attribute or some 3D effects can play a role, it can be simulated in 3D scene picture (picture Z attribute axis), allowing multiple layers to produce perspective, and also to support the light and material, projection / shadow. In software, the camera is virtual and will not have a track or a three stand for you to use. However, by simulating real-life camera movements, we can migrate our understanding to virtual cameras, so that we can better understand and apply it. In practice, we will find that the simple on camera animation, it is very difficult to achieve the desired effect, not our understanding is not enough, but the software itself is designed so that we need to find some way to simulate the camera in real life, to achieve our purpose.

ArcGIS camera plug-in preparation

First, define the camera interface; At this point, our plug-in program is complete, the plug-in code generation plug-in format files, import into the ArcGis software database, open the software can be invoked. Figure 5 is the effect of observing the block data in the software after calling the camera plug-in.

#
Figure 5 Integrated camera plug-in after the soil moisture data observation effect diagram

Conclusion

With the help of data visualization is the main graphical means to communicate clearly and effectively communicate with the information in order to effectively convey the concept, through the convey and key characteristics, in order to achieve in-depth insight into the rather sparse and complex data sets.

Research and experimental results show that

(1) based on the ArcGis software SDK source code to write the auxiliary plug-in and ArcGis software integration, to make up for the ArcGis software auxiliary function is not perfect, and optimize the visualization effect of agricultural monitoring data is feasible.

(2) integrated plug-in ArcGis software in the field monitoring data, obvious optimization effect of 3D visualization of monitoring data, the two-dimensional spatial variability of thematic map than the traditional representation of 3D data more natural, more clear, more intuitive. The results of the study are of great significance to the development of modern agriculture.

  and signals". The generation and acquisition of primitive symbols, structures, images and signals are the tasks of scientists or engineers engaged in scientific and engineering calculations and experiments. How to transform these data into visual graphics (images) is the task of scientists and engineers engaged in scientific computing visualization applications. They combine with scientists engaged in basic research or engineering research, and use visual tools to implement and complete the process of converting invisible objects into visible images. Usually, the processing of these raw data should include the following two aspects: Raw data preprocessing and storage are usually used to visualize raw data, and the preprocessing section involves the following operations: 1) data format and standardization; 2) data transformation technology; 3) data compression and decompression. The data processed through the above may have larger data capacity. It is usually stored in a standardized data interchange format. The NASA and other units of the National Supercomputing Center study the format and structure of scientific data storage and exchange, such as HDF, CDF, and NetCDF. These data formats are typically self describing, easy to store and read, and are widely used in large-scale numerical computations. Through the sensor Nong'an County Joint jiulongzhen soil moisture real-time data as shown in figure 1.

Figure 1 moisture

 1 Figure 1 moisture Nong'an County Helong town farmland real time data

( 2 )

 2 data format conversion Using the development program of AutoCAD2000 software, we can draw the above hierarchical data files, eliminate the error points and save the.Dwg format. (3) loading elevation point data; Open Arcsence and add the "point" layer of the saved.Dwg data to the layer. Since the elevation of each stratum is different, the multilayer data is stacked up and down after loading the multi layer data Point layer. (4) establishing TIN; Using create, Analyst, create/modify, TIN, TIN, from, and features under 3D, choose to load the point layer to build the TIN and build the TIN (Elevation) to build the TIN. As shown in figure 2.

Figure 2

 2 Figure 2 the establishment of Nong'an County Helong town experimental field plots TIN model

Figure 3

 3 Figure 3 markers in Nong'an County Helong town experimental field plots TIN boundary (6) build DEM stretch layer When you open the ArcScene and load the inner boundary, you can see that it is at the height of the two layer point, adjusted properly until the inner edge of the interface is in the middle of the two layer. Again, open the create TIN from Elevation create/modify Analyst under 3D TIN features, load the geological points and the inner edge interface, select a layer, set the Height source to the other attribute in the default format. Then select its corresponding inner edge interface, set its Triangulate attribute to hard erase, and the other two items to None. Select the saved path and click OK to build an upright TIN.

Figure 4

 4 Figure 4 Visual expression of soil moisture in farmland of Nong'an County Helong town

  The borehole coordinate results table.Xls table data and coal quality testing results table.Xls table finishing, let it become a format for the "number, X, Y, Z for the superposition of the formation of documents. Converts to text files for subsequent program processing.

Acknowledgments

This work was funded by the China Spark Program.2015GA660004."Integration and demonstration of corn precise operation technology based on Internet of things"