An optimization method, device, medium and product for an automatic-differentiable high-precision numerical weather prediction method

Abstract

The application discloses an optimization method, equipment, medium and product of a high-precision numerical weather prediction method capable of automatic differentiation, relates to the technical field of weather prediction, and comprises the following steps: acquiring an initial atmospheric state, terrain data and physical parameters; generating a cubic spherical grid, defining internal units and shadow units; precalculating a linear mapping matrix from the internal units to the shadow units based on the grid geometry structure; generating an initial field at the nth moment and calculating the average value of a source term; reconstructing the boundary integral point value based on the average value of the internal units; solving a Riemann problem to obtain a boundary flux, and integrating to obtain an average flux; calculating a time trend in combination with the average value of the source term; time-integrating to obtain a predicted field of the internal units at the moment; and filling the shadow units by using the linear mapping matrix to obtain a predicted field at the moment. The application can significantly improve the calculation efficiency and has an automatic differentiation function.

Description

Technical Field

[0001] This application relates to the field of weather forecasting technology, and in particular to an optimization method, device, medium, and product for a high-precision numerical weather forecasting method with automatic differentiation. Background Technology

[0002] The core task of numerical weather prediction is to obtain the weather conditions at future times by discretely solving atmospheric dynamic equations. As resolution continues to improve, numerical models are increasingly demanding higher computational efficiency, accuracy, and scalability.

[0003] Existing numerical weather prediction models typically employ iterative interpolation to fill shadow cells when dealing with panel boundaries. This requires multiple iterative calculations and inter-node communication at each time step, resulting in low computational efficiency and hindering the full utilization of the performance advantages of heterogeneous computing platforms such as GPUs. Furthermore, traditional models are mostly developed using languages ​​like Fortran, lacking automatic differentiation capabilities and thus limiting deep integration with artificial intelligence models. This restricts the development of "numerical + AI" integrated forecasting systems.

[0004] In addition, existing high-order numerical weather prediction methods typically require multiple sampling points for Gaussian integration at cell boundaries and within cells, resulting in high computational and storage costs, making it difficult to meet the needs of ultra-high resolution numerical prediction scenarios. Summary of the Invention

[0005] The purpose of this application is to provide an optimization method, device, medium, and product for a high-precision numerical weather prediction method with automatic differentiation, which can significantly improve computational efficiency and has automatic differentiation function.

[0006] To achieve the above objectives, this application provides the following solution: In a first aspect, this application provides an optimization method for a high-precision numerical weather prediction method with automatic differentiation, comprising: Acquire initial atmospheric state data and topographic data, and configure physical parameters based on geophysical constants; Based on the physical parameters, a cubic sphere mesh is generated and the sphere is decomposed into six panels; on each panel, internal cells for solving the forecast quantity and shadow cells located outside the panel boundary for handling boundary conditions are defined. Based on the geometry of the cubic sphere mesh, a linear mapping matrix from the internal cells to the shadow cells is pre-calculated; the linear mapping matrix is ​​stored in a sparse matrix format. The initial atmospheric state data is mapped onto the cubic sphere grid to generate the initial field at time n, and the source term average value of each of the internal cells is calculated based on the terrain data and the physical parameters; the initial field at time n includes the internal cell mean value stored on each of the internal cells; Based on the mean of the internal units, spatial reconstruction is performed on each internal unit. The distribution of physical quantities inside the unit is constructed using tensor product polynomials. The integration point value on the boundary of the internal unit is obtained by using the midpoint rule as a single sampling point. Using the integration point values ​​on the boundary of the inner element as the left and right states of the Riemann problem, the Riemann problem is solved using a low Mach number approximate Riemann solver to obtain the flux value at each integration point on the boundary of the inner element. Numerical integration is then performed along each boundary of the inner element to obtain the average flux on each boundary of the inner element. Based on the average flux on each boundary of the internal unit and the average source term of the internal unit, the time tendency of the mean of the internal unit with time is calculated according to the finite volume method; the average source term is obtained by integration with the center point of the unit as a single sampling point through the midpoint rule. Based on the aforementioned time tendency, an explicit Runge-Kutta time integration scheme is used to perform time integration, and the internal cell prediction field at time n+1 is obtained. Using the linear mapping matrix, for the The internal cell prediction field at time 1 is subjected to sparse matrix-vector multiplication to calculate the filling values ​​of all shadow cells, thus obtaining a prediction field containing complete shadow cells. Time forecast field.

[0007] In a second aspect, this application provides a computer device, including: a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the steps of the optimization method for the automatically differentiable high-precision numerical weather forecasting method described in any one of the above-mentioned methods.

[0008] Thirdly, this application provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the steps of the optimization method for the automatically differentiable high-precision numerical weather forecasting method described above.

[0009] Fourthly, this application provides a computer program product, including a computer program that, when executed by a processor, implements the steps of the optimization method for the automatically differentiable high-precision numerical weather forecasting method described above.

[0010] According to the specific embodiments provided in this application, this application has the following technical effects: This application provides an optimization method, apparatus, medium, and product for a high-precision numerical weather prediction method with automatic differentiation. By utilizing a polynomial basis to construct a numerical integration scheme, this application overcomes the problem of excessive sampling points in traditional Gaussian integration, requiring only a single sampling point at both cell boundaries and centers to achieve high-precision integration. This optimization scheme significantly reduces the number of reconstruction points, integration points, and storage overhead while maintaining arbitrary high-order accuracy, oscillation-free operation, strict mass conservation, and the inherent automatic differentiation characteristics. This further improves GPU parallel computing efficiency and is more suitable for ultra-high-resolution numerical weather prediction scenarios. This application pre-computes the linear mapping matrix from internal cells to shadow cells, transforming the shadow cell filling calculation, which traditionally requires multiple iterations per time step, into an offline task completed in one go during the initialization phase. During real-time integration, shadow cell filling is simplified to a single sparse matrix-vector multiplication operation, greatly reducing computational load. Since iterative calculations are no longer required at each time step, the number of communication steps between distributed computing nodes is also reduced, significantly improving the computational efficiency and parallel scalability of the model, providing efficient algorithmic support for high-resolution, long-term numerical forecasting. This application utilizes cubic spherical grids for numerical discretization, transforming core operations such as spatial reconstruction, flux integration, and time progression into computational forms suitable for parallel execution on heterogeneous computing platforms like GPUs. By organizing forecast data into a data structure suitable for batch processing, the massively parallel computing capabilities of GPUs can be fully utilized, effectively addressing the computational bottleneck faced by traditional CPU platforms in high-resolution computing. This application is based on an automatic differentiable framework, possessing automatic differentiating capabilities across all models. This feature enables seamless integration with neural network models without the need to develop additional cross-language interface modules. On one hand, the automatic differentiating function allows for direct end-to-end optimization of physically parameterized models; on the other hand, it greatly simplifies the development of model adjoint operators in four-dimensional variational assimilation. This lays the technical foundation for building a "numerical + AI" integrated forecasting system, enabling models to maintain the advantages of traditional numerical methods while possessing the automatic learning capabilities of artificial intelligence. Attached Figure Description

[0011] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0012] Figure 1 This is an application environment diagram of an optimization method for an automatically differentiable high-precision numerical weather prediction method according to an embodiment of this application; Figure 2A flowchart illustrating an optimization method for a high-precision numerical weather prediction method with automatic differentiation, provided as an embodiment of this application; Figure 3 This is a schematic diagram of a cubic mesh and computational space provided in an embodiment of this application; Figure 4 A schematic diagram of the function point distribution of the unit and boundary provided in an embodiment of this application; Figure 5 is a schematic diagram of points and elements near the panel boundary provided in an embodiment of this application; wherein, Figure 5(a) is the flux point on the interface between panel 1 and panel 4; Figure 5(b) is the shadow element outside the boundary of panel 4; Figure 5(c) is the template required for 5th order shadow element interpolation; Figure 6 This is a schematic diagram of the function point distribution on a unit provided in an embodiment of this application; Figure 7 This is a schematic diagram of a reconstructed coordinate system and a polynomial reconstruction template provided in an embodiment of this application; Figure 8 This is a schematic diagram of a 5th-order WENO 2D template provided in an embodiment of this application; Figure 9 This is a schematic diagram of the HOPE mode algorithm provided in an embodiment of this application; Figure 10 This is a schematic diagram of the structure of a computer device provided in an embodiment of this application. Detailed Implementation

[0013] The technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, and not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0014] This application relates to the field of weather forecasting technology, specifically to the High Order Prediction Environment (HOPE), which provides a numerical weather prediction method that combines high accuracy, oscillation-free operation, scalability, mass conservation, and automatic differentiability. It can accurately solve the spherical shallow water equations, simulate barotropic atmospheric flow fields, fully utilize the computing performance of GPU heterogeneous platforms, and is easily integrated with AI algorithms to construct numerical weather prediction models. This solves the software integration problem between traditional numerical models and artificial intelligence models.

[0015] Numerical weather prediction is a crucial support for modern weather forecasting operations. Its core function is to discretely solve atmospheric dynamic equations to obtain future weather conditions. Shallow water equations are a simplified set of atmospheric motion equations used to describe fluid free-interface fluctuations, such as free flow in the upper atmosphere and storm surges in the ocean. Shallow water equations are nonlinear partial differential equations, and no direct analytical general solution can be given. The task of shallow water numerical models is to obtain accurate numerical solutions to these equations through discrete solutions. In atmospheric science, shallow water models are also widely used to verify key elements of the solution algorithm, such as accuracy, conservation, and computational efficiency, to evaluate the algorithm's potential. This application focuses on solving shallow water equations and proposes a new method for high-precision discrete solutions on a sphere.

[0016] To make the above-mentioned objectives, features and advantages of this application more apparent and understandable, the application will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0017] This application provides an optimization method for a high-precision numerical weather prediction method with automatic differentiation, which can be applied to, for example... Figure 1 The application environment shown is illustrated. Terminal 102 communicates with server 104 via a network. A data storage system can store the data that server 104 needs to process. The data storage system can be set up independently, integrated into server 104, or placed in the cloud or on another server.

[0018] The terminal 102 can be, but is not limited to, various desktop computers, laptops, smartphones, tablets, IoT devices, and portable wearable devices. IoT devices can include smart speakers, smart TVs, smart air conditioners, and smart in-vehicle devices. Portable wearable devices can include smartwatches, smart bracelets, and head-mounted devices. The server 104 can be implemented using a standalone server or a server cluster composed of multiple servers, or it can be a cloud server.

[0019] In one exemplary embodiment, such as Figure 2 As shown, an optimization method for a high-precision numerical weather prediction method with automatic differentiation is provided. This method is executed by a computer device, specifically by a terminal or server alone, or by both a terminal and a server. In this embodiment, the method is applied to... Figure 1 Taking server 104 as an example, the explanation includes the following steps S1 to S9. Wherein: S1. Acquire initial atmospheric state data and terrain data, and configure physical parameters based on geophysical constants; the physical parameters in this embodiment may include gravitational acceleration, Earth's rotational angular velocity, Coriolis parameter, and Earth's radius.

[0020] S2. Based on the physical parameters, generate a cubic sphere mesh and decompose the sphere into six panels; define internal cells for solving the forecast quantity and shadow cells located outside the panel boundary for handling boundary conditions on each panel.

[0021] S3. Based on the geometry of the cubic sphere mesh, pre-calculate the linear mapping matrix from the internal cells to the shadow cells; the linear mapping matrix is ​​stored in a sparse matrix format.

[0022] S4. Map the initial atmospheric state data to the cubic sphere grid to generate the initial field at time n, and calculate the source term average value of each internal cell based on the terrain data and the physical parameters; the initial field at time n includes the internal cell average value stored on each internal cell.

[0023] S5. Based on the mean of the internal units, spatial reconstruction is performed on each internal unit. The distribution of physical quantities inside the unit is constructed using tensor product polynomials. The integral point value on the boundary of the internal unit is obtained by using the midpoint rule as a single sampling point.

[0024] S6. Using the integration point values ​​on the boundary of the inner element as the left and right states of the Riemann problem, the Riemann problem is solved using a low Mach number approximate Riemann solver to obtain the flux value at each integration point on the boundary of the inner element. Numerical integration is then performed along each boundary of the inner element to obtain the average flux on each boundary of the inner element.

[0025] S7. Based on the average flux on each boundary of the internal unit and the average source term of the internal unit, calculate the time tendency of the mean of the internal unit with time according to the finite volume method; the average source term is obtained by integration with the center point of the unit as a single sampling point through the midpoint rule.

[0026] S8. Based on the aforementioned time tendency, time integration is performed using the explicit Runge-Kutta time integration scheme to obtain the internal cell prediction field at time n+1.

[0027] S9. Using the linear mapping matrix, for the... The internal cell prediction field at time 1 is subjected to sparse matrix-vector multiplication to calculate the filling values ​​of all shadow cells, thus obtaining a prediction field containing complete shadow cells. Time forecast field.

[0028] In another exemplary embodiment of this application, the linear mapping matrix is ​​used to... The internal cell prediction field at time 1 is subjected to sparse matrix-vector multiplication to calculate the filling values ​​of all shadow cells, thus obtaining a prediction field containing complete shadow cells. Following the time forecast field, it also includes: Based on the aforementioned complete shadow unit The system forecasts the field at any time and outputs the solutions of meteorological elements as numerical weather prediction results. The solutions of the meteorological elements include one or more of the following: wind field, height field, vorticity, divergence, potential vorticity, mass conservation quantity, and energy conservation quantity.

[0029] In another exemplary embodiment of this application, the linear mapping matrix is ​​used to... The internal cell prediction field at time 1 is subjected to sparse matrix-vector multiplication to calculate the filling values ​​of all shadow cells, thus obtaining a prediction field containing complete shadow cells. Following the time forecast field, it also includes: Based on the aforementioned complete shadow unit The system forecasts the field at any time and outputs the solutions of meteorological elements as numerical weather prediction results. The solutions of the meteorological elements include one or more of the following: wind field, height field, vorticity, divergence, potential vorticity, mass conservation quantity, and energy conservation quantity.

[0030] In another exemplary embodiment of this application, step S1 mainly includes front-end data preparation.

[0031] Step S11, Initial Atmospheric State Data: Obtaining the global geopotential height field ( Wind speed and zonal wind speed: obtained from reanalysis data or ideal experimental initial fields. and meridional wind speed This formula is used to calculate the core variables in meteorology. This represents the potential height, and its unit is meters (m). 2 / s 2 g represents the acceleration due to gravity, and its unit is m / s². 2 Standard value g = 9.80616 m / s 2 h represents the thickness of the atmosphere in the vertical direction, with the unit being meters (m).

[0032] Step S12, Terrain Data: Obtaining surface elevation ( ) is used to calculate the terrain forcing term ( These data typically come from the terrain settings in a digital elevation model (DEM) or ideal test. Since solving partial differential equations requires initial and boundary conditions to ensure the well-posedness of the solution, the reanalysis data and DEMs mentioned here only provide examples of initial and boundary conditions. These data can directly provide the initial and boundary values ​​required by the equations.

[0033] Step S13, Physical parameter configuration: Including gravitational acceleration (g), Earth's rotational angular velocity (ωt), and so on. Constants such as Coriolis parameter (f) and Earth radius (r) are also included.

[0034] In another exemplary embodiment of this application, the cubic sphere mesh generation in step S2 mainly includes the following steps.

[0035] Please see Figure 3 ,in, Figure 3 In this example, (a) represents the cubic sphere mesh, and (b) represents the computational space. An isoangular cubic sphere mesh is used in this embodiment, decomposing the sphere into six faces, each mapped to a rectangular computational space for subsequent discretization of the equations.

[0036] The isoangular cubic sphere mesh is a mesh generation method for spherical numerical computation. Its core idea is to divide the six faces of a cube into an isoangular grid, and then project this grid onto a circumscribed sphere through a central projection, thereby constructing a seamless computational grid on the sphere consisting of six identical regions.

[0037] Foundations of geometry: Initial definition: There exists a cube with side length 2a, whose eight vertices lie on a radius of... On the sphere (i.e., the circumsphere of the cube, and) ).

[0038] Mesh drawing: On each face of the cube, establish a coordinate system based on the isoangular center. The local computational grid, in which This is the face number.

[0039] Projection mapping: By projecting from the center, the grid points defined on each face of the cube are mapped. Mapped to the corresponding point on the circumscribed sphere.

[0040] The governing equations are a set of spherical shallow water motion equations in generalized curvilinear coordinate system form. The transformation matrix between the spherical momentum vector and the inverse momentum vector in generalized curvilinear coordinate system is: ; Based on the coordinate mapping relationship of the cubic sphere grid, it can be further written as: ; Where p is Figure 3 The corresponding face numbers of the cube sphere are 1 to 6, corresponding to the 6 panels of the cube sphere. They are longitude and latitude, respectively. This represents the longitude of the p-th cubic sphere. The transformation matrix parameters represent the 5th and 6th panels (pole panels) of the cubic sphere, which are used to correct the coordinate mapping specialness of the pole region and ensure the physical consistency of the momentum vector in this region during cross-panel transformation.

[0041] This represents the auxiliary parameter for the momentum vector transformation matrix of the 5th and 6th panels of the cube sphere (corresponding to the pole panels of the North and South Pole regions). It is based on the trigonometric function relationship between longitude λ and geocentric latitude θ, correcting the nonlinear deviation of the coordinate mapping in the pole region. They jointly constructed the momentum transformation matrix for this region. .

[0042] ; Covariant metric tensor ( Inverse metric tensor , measuring Jacobi .

[0043] In another exemplary embodiment of this application, the pre-calculation of the linear mapping matrix from the inner element to the shadow element based on the geometry of the cubic sphere mesh in step S3 can be replaced by the following step: S31. Initialize the value of all shadow units to zero.

[0044] S32. In each iteration, based on the currently known mean of the internal cells and the value of the shadow cell in the current iteration step, the tensor product polynomial interpolation method is used to calculate the value of all Gaussian integral points in each shadow cell.

[0045] S33. Perform Gaussian integration on the values ​​of all Gaussian integration points within each shadow unit to obtain the new average value of the shadow unit, thereby obtaining the updated values ​​of all shadow units.

[0046] S34. Continue until the norm of the difference between the values ​​of the shadow units obtained from the two iterations is less than the preset threshold, and the convergence state is reached.

[0047] S35. Treat the value of the shadow unit in the convergent state as a linear function of the mean of the inner unit, use the automatic differentiation framework to differentiate the linear function to obtain the Jacobian matrix from the inner unit to all shadow units, and use the Jacobian matrix as the linear mapping matrix.

[0048] Specifically, the linear mapping matrix is ​​pre-computed. This transforms the shadow cell filling calculation, which is traditionally required at each time step and may involve iterations, into an offline task that is completed in one go during the initialization phase. In the subsequent real-time integration, shadow cell filling is simplified to a single efficient (sparse matrix-vector multiplication) operation, thereby greatly improving computational efficiency, reducing the number of communications, and maintaining the algorithm's high accuracy.

[0049] In HOPE, the skewed single-sided shadow unit reconstruction scheme of previous studies is abandoned, solving the problem that the computational accuracy cannot exceed the 7th order (excluding 7). A central reconstruction scheme is adopted. In the mode initialization stage, the unknown unit information is first filled in through an iterative algorithm. Then, based on the linear characteristics of the mapping algorithm, the Jacobian matrix of the mapping process is obtained by utilizing the automatic differentiability property of the model. (This matrix is ​​the linear mapping corresponding to the aforementioned iterative process), and save the matrix as a CSR (Compressed Sparse Row) sparse matrix.

[0050] The matrix A deterministic linear relationship was established from the forecast values ​​of all internal cells to the filling values ​​of all shadow cells. In subsequent calculations, the entire filling process is directly expressed as the multiplication of the mapping matrix and the atmospheric state vector, eliminating the need for iterative calculations, greatly improving computational efficiency, and reducing the number of communications between distributed computing nodes.

[0051] This matrix can be pre-computed and stored as a sparse matrix in CSR (Compressed Sparse Rows) format during pattern initialization. This matrix depends only on the mesh geometry and the Gaussian integration scheme and does not change over time. Efficiently generated by the automatic differentiation framework (PyTorch), it represents a linear mapping relationship of "known internal element → shadow element".

[0052] Please see Figure 4 It represents the distribution of function points in the unit and boundary. Blue represents the halo region (shadow region), white represents the prediction region (internal calculation unit), purple points are Gaussian integral points (GQP) in the shadow unit, and green points are the center points (CP).

[0053] Combination Figure 4 The goal of this step is to use the known cell average value (green dots) of the forecast area (white) to calculate the value of all GQPs (purple dots) in the shadow cell (blue) through high-precision interpolation, and then obtain the filling value of the shadow cell through integration.

[0054] The calculation process is as follows: Shadow unit value initialization: Initialize the value of the shadow unit to zero: ,in The superscript indicates the field value of the shadow unit, and the superscript indicates the iteration number.

[0055] Tensor product polynomial (TPP) interpolation: This step interpolates Gaussian Quadrature Points (GQPs) within the shadow cell using known field values ​​from adjacent panels: for each GQP, a reconstruction template on the adjacent panel is selected, and the GQP value is reconstructed using the Tensor Product Polynomial (TPP). It should be noted that the GQP selection method involves first defining the coordinates of the GQP within the cell of the computational space, then mapping them to spherical coordinates using a cube sphere generation algorithm, and finally calculating their corresponding metric Jacobi and metric tensor.

[0056] Please refer to Figure 5, where Figure 5(a) shows the flux points at the interface between panel 1 and panel 4: the flux points (red dots) located at the interface between the two panels, and the flux values ​​of these points are determined by the Riemann solver. Figure 5(b) shows the shadow elements outside the boundary of panel 4: the shaded area represents the shadow elements outside the boundary of panel 4, and the green dots represent the integration points (GQPs) within these ghost elements. The values ​​of these points need to be obtained from the known field values ​​of the adjacent panels through interpolation. Figure 5(c) shows the template required for 5th-order shadow element interpolation: the red dots represent the target GQPs located within the shadow elements of panel 4. The blue shaded area is the Tensor Product Polynomial (TPP) reconstruction template on panel 1, used to interpolate the values ​​of these red GQPs. This template covers multiple elements of the adjacent panels to ensure high-precision interpolation.

[0057] Iterative update of shadow unit values: (1) Interpolate all Gaussian integral points (GQP): Using the current field value (known element value) and current shadow unit value Interpolate all GQPs.

[0058] (2) Calculate the new shadow unit value: The average value of the shadow unit is obtained by integrating the GQP value using the Gaussian integral formula, and then the updated shadow unit value is obtained. .

[0059] (3) Convergence judgment: Calculate residuals ,like (in double precision) If the iteration stops, the iteration continues; otherwise, the previous steps are repeated. The purpose of this sub-step is to determine whether the iterative process of filling the shadow unit values ​​has reached sufficient accuracy. It is a "stopping condition" that tells the algorithm when it can end the iteration and use the currently obtained shadow unit values ​​instead of calculating indefinitely.

[0060] Extracting the linear mapping matrix : When the iteration converges, the final shadow unit value It can be represented as an internal unit value Linear functions: ; matrix This is the desired linear mapping matrix. In practice, this matrix can be efficiently generated using the automatic differentiation framework (PyTorch): by automatically differentiating the above iteratively converged mapping relationship, its Jacobian matrix is ​​directly obtained, i.e. .

[0061] In another exemplary embodiment of this application, step S4 above, which generates an initial field based on atmospheric conditions and a cubic sphere grid, mainly includes the following steps.

[0062] An initial field matching a cubic sphere grid is generated from the atmospheric state. In HOPE, atmospheric state point values ​​at the centers of all cells are first generated based on the scene, and the atmospheric state field is then converted into model forecasts. ,in: ; In the formula, Indicates the potential height; u s This refers to the zonal wind speed; v s The meridional wind speed; via u s and v s u and v are calculated.

[0063] The next step is a crucial part of numerical weather prediction—reconstruction. Mathematically, reconstruction involves reconstructing a continuous variable field from discretely distributed cell average variables. Once the continuous variable field is obtained, point values ​​at any location can be acquired. Here, only the values ​​of EP (integration point at the cell boundary) and GQP (Gaussian integration point within the cell) are needed to complete subsequent calculations. At the cell boundary, only a single sampling point is required for integration. Therefore, the reconstruction process does not require explicit polynomial expressions; instead, it only needs to obtain the mapping matrix from the discrete mean field to the function points, thus significantly reducing computational costs.

[0064] Tensor Product Polynomial (TPP) Reconstruction: The average atmospheric state in each cell is obtained through polynomial integration and used as the initial state of the model.

[0065] In this embodiment, the discrete atmospheric variable field is reconstructed using tensor product polynomials. Compared with directly using a two-dimensional polynomial basis, the terms in the tensor product polynomial can be naturally matched with the square template, ensuring the well-posedness of the reconstructed linear system solution.

[0066] Constructing a polynomial expression for space reconstruction: Reconstruction is based on the unit average (i.e., the forecast). As known quantities, construct a polynomial function, and then use the polynomial function to obtain the edge points and Gaussian integration points. The value is obtained by averaging the values ​​of various functional points around the unit.

[0067] In the above text, a cubic sphere mesh has been generated. The area around the target cell is then selected on the cubic sphere mesh. A rectangular area composed of neighboring units serves as a template (e.g., TPP3 corresponds to this). Template, TPP5 equivalent (Template), the position of each cell in the template is calculated using local coordinates. OK. The template is used for solving the coefficients of the tensor product polynomial (TPP) reconstruction.

[0068] Horizontal reconstruction using TPP: The nth-order TPP reconstruction scheme TPPn is represented by a linear combination of two-dimensional polynomials, as follows: ; Where i and j are two-dimensional row and column indices in the template grid, i represents the index value in the x-direction within the template, ranging from 1 to n; j represents the index value in the y-direction within the template, ranging from 1 to n; and n represents the template width. For polynomial coefficients, This is the local coordinate system for the template. (The remaining text appears to be incomplete and possibly contains errors. A more accurate translation would require the full context.) A single neighboring unit can uniquely determine this. The value of , .

[0069] Establish a system of linear equations and solve for the polynomial coefficients: Based on the center point values ​​of all elements in the template ,in The total number of units participating in the reconstruction within the template. Construct a system of linear equations, that is ; Solving this system of linear equations will yield the reconstructed polynomial. coefficient This uniquely determines the polynomial.

[0070] The steps to solve for the coefficients are as follows: ①Establish a linear system: Based on the volume fractional average (VIA) of each cell. Construct a system of linear equations: ; in: It is a basis function integral matrix with a size of ,in This is the total number of cells in the template. It is a VIA vector with elements of , ..., ; ; It is the vector of polynomial coefficients to be found, i.e. .

[0071] ② Matrix inversion and coefficient solving: because Since it is a square matrix and not singular, we can directly find its inverse: ; This step directly yields the reconstructed polynomial coefficient vector. That is, all polynomial coefficients .

[0072] Due to the matrix It relies solely on the mesh geometry, is calculated and stored once during mode initialization, and subsequent calculations only require matrix-vector multiplication.

[0073] ③ Reconstruct point values ​​using coefficients: ; The obtained coefficients Substituting into the polynomial expression, it is possible to calculate any position. The function value.

[0074] Calculate the average volume fraction: Volume fraction average (VIA) The initial atmospheric conditions are obtained by performing a polynomial integration on the cubic sphere grid cells. During model initialization, the initial atmospheric conditions (potential height, wind speed field) are interpolated onto the cubic sphere grid, and the VIA of each cell is calculated as the starting point for the forecast. The VIA of each cell is obtained by integrating the initial conditions during model initialization.

[0075] The polynomial in each unit Integrating within the inner zone yields the average field quantity of the cell: ; in, Represents a cell within the computational space of The side length, in HOPE, is locally reconstructed using computational space, therefore let It can greatly simplify the calculation without affecting the result.

[0076] Specific ways to combine data processing with grids: ① Gridding: The global atmospheric field is mapped onto six panels of a cubic sphere, each panel being a rectangular computational grid.

[0077] ②Integral calculation: Within each cell, the volume average (VIA) of the field quantity is calculated using the polynomial integration method.

[0078] ③ Tensor storage format: To achieve GPU-friendly and efficient computation, this embodiment organizes the forecast field data into a five-dimensional tensor format to be compatible with PyTorch's convolution operations: VIA data is organized into tensor form. This facilitates convolutional reconstruction and integral calculation. =3 represents the number of variables (including mass, x-direction force, and y-direction momentum). =6 represents the number of panels. The number of grids for each panel.

[0079] ④ Time integration: During the time integration process, VIA, as the core forecast variable, is advanced to the next time step through the process of reconstruction → boundary mapping → Riemann solution → flux integration → time update.

[0080] This step yields the initial states of each unit in the model, providing the data foundation for subsequent forecasts, such as... Figure 6 As shown. Figure 6 The red dots represent edge points (EP), and the green dots represent cell center points (CP).

[0081] In another exemplary embodiment of this application, step S5, which involves spatial reconstruction based on the initial field to obtain cell boundary values, mainly includes the following steps.

[0082] The above text completes the basic preparations: (1) generating the initial field of the matching cubic sphere mesh (the volume integral average field value of each element, such as the element center point value of physical quantities such as potential height and wind speed), which is the core input physical quantity for reconstruction; (2) determining The selection range of neighboring units, that is, the set of units participating in WENO2D reconstruction; (3) Calculate the side length of the spatial unit. , ,( (4) Solve for the polynomial coefficients of the TPP reconstruction. As the basis for the accuracy of low-order sub-templates in WENO2D reconstruction, it is used to calculate the reconstruction value of each sub-template; (5) the relative coordinates of the elements. It is used for the construction of WENO2D basis functions and weight calculation.

[0083] Based on the above preparations, this step constructs two-dimensional reconstruction weights using tensor product polynomials to obtain cell boundary values. , Used for subsequent Riemannian solutions and time integration, such as Figure 7As shown, (a) is the local reconstruction coordinate system, with the red dots representing the cell centers; (b) is the quadratic polynomial template; (c) is the TPP3 template; and (d) is the TPP5 template.

[0084] Step S51, Spatial Reconstruction: Spatial reconstruction adopts the Tensor Product Polynomial (TPP) reconstruction scheme.

[0085] Step S51.1, Tensor Product Polynomial (TPP) Reconstruction: The preceding steps have prepared the data and templates. In the current step S51.1, the steps from the previous step are invoked. The neighbor cell template is used to determine the range of neighbor cells corresponding to the current target cell; from the initial field in step two, the average field value of all neighbor cells within this template is extracted. .

[0086] Use the same polynomial expression as in the step of "constructing a polynomial expression for space reconstruction". The polynomial is in each unit Integrating within the inner boundary, we obtain the form of the unit mean field: ; in, Represents a cell within the computational space of The side length is as described above. c is the position coefficient vector. , of which Each element is expressed as , , This is the floor function. Polynomial coefficient vector. 。

[0087] Step S51.2: Construct a linear system: Based on the unit mean values ​​obtained in the "calculate the volume integral average value" step, a system of linear equations is formed by combining them with the expression in step S51.1: ; in Let A be a matrix formed by the integrals of the basis functions (each row of matrix A is the transpose of psi, and ps is a column vector), where the elements are... , Let be the polynomial coefficient vector. This is the unit average field vector.

[0088] in This is the total number of cells in the template. Here... This represents the average field vector of the Nth unit. Let N represent the basis function of the polynomial.

[0089] Step S51.3: Solve for the polynomial coefficients: Solving the system of linear equations using matrix inversion or the least squares method yields the polynomial coefficients. : .

[0090] Step S51.4, Point Value Reconstruction: Using the obtained coefficients, calculate the field values ​​at specific points (including edge point EP and element center point CP): ; in superscript This represents the transpose of the matrix. To reconstruct the matrix, it should be noted that A matrix depends only on the geometric properties of the computational space, which does not change over time. Therefore, in practical programs, only one matrix needs to be computed during initialization. The reconstruction matrix is ​​saved to memory (video memory). When spatial reconstruction is needed, only the reconstruction matrix needs to be saved to memory (video memory). The matrix is ​​then assembled into the convolution operator `torch.nn.functional.conv2d`. `M` represents the total number of target sampling points participating in this point value reconstruction, which determines the row dimension of vector `p` and matrix `C`. `M` represents the variable index. This represents the local coordinates of the Mth reconstructed point.

[0091] Step S52: Obtain cell boundary values , : Based on the above TPP and WENO 2D reconstruction, the cell boundary values ​​are obtained through high-precision polynomial approximation and adaptive weight adjustment. and : in For potential height, The phase velocity of the gravitational wave is indicated by the subscript. These represent the left and right sides of the boundary, respectively. Normal wind speed at the unit boundary: .

[0092] This formula converts the inverter wind u into a physical velocity perpendicular to the cell boundary. , where variable i is used to extract from the contravariant tensor Select the corresponding component from the list.

[0093] In another exemplary embodiment of this application, before step S6, a step of performing an equivalent mapping of the cubic spherical boundary based on the spatial reconstruction result is further included: The essence of this step is to address the discontinuity of cubic sphere mesh coordinates across panels. It does not introduce new physical quantities; instead, it performs a coordinate system-wide physical quantity transformation on the panel boundary data reconstructed above. This allows the subsequent Riemann solution to calculate the flux across the panel boundaries while maintaining physical consistency. This process ensures the algorithm's conservation and high accuracy.

[0094] The core of this process is to utilize the coordinate transformation relationships on different panels to remap the physical quantities of the mass and momentum fields at the panel boundaries, ensuring physical consistency of flux calculations at these boundaries. Through mass and momentum field transformations, combined with Gaussian integration point processing at the cubic sphere boundary, vector uniformity at the cubic sphere mesh boundary is ensured, providing accurate input for subsequent Riemannian solutions and time integration.

[0095] The data required for this step includes: (1) Gaussian integral point (EP) values ​​on the cell boundary: derived from the spatial reconstruction results. On the panel boundary, these EP points are shared by two adjacent panels (e.g., Figure 7 (a) Red dots on the interface between panel 1 and panel 4. (2) Coordinate transformation matrix of adjacent panels. Used to transform physical quantities from local coordinates of one panel to another.

[0096] The key to this step lies in the "conversion of physical quantities between panels," which is specifically divided into two sub-steps: (1) Mass field transformation: For each Gaussian integration point on the panel boundary, the mass field value on a panel (such as panel 1) needs to be calculated. Transform to the coordinate system of an adjacent panel (such as panel 4). The transformation formula is as follows: ; In the formula, the subscript represents the target panel, and the superscript represents the source panel. This is a Jacobian determinant that measures the area of ​​the mass field, varying with the panel and its location. It describes the rate of change of the area mapped from computational space (rectangular grids on a cubic sphere panel) to physical space (the sphere). This scaling ensures consistency in the measurement of the mass field across the panel boundaries. In the formula... This represents the mass field value after transforming from the coordinate system of panel 1 to panel 4; Represents the original mass field value in the coordinate system of panel 1; transformation factor Used to correct differences between two panels caused by different metrics.

[0097] in ; ; It is a covariant metric tensor, a 2×2 tensor used to describe the local geometry in a curvilinear coordinate system, where i,j represent coordinates; For covariant metric tensors The determinant value; X and Y both represent the tangential components of the local isoangular coordinates, where These are the local isoangular coordinates of the panel. , is the distance factor. r represents the Earth's radius, r .

[0098] (2) Momentum field transformation: Momentum is a vector quantity, and its physical direction must remain unchanged during the conversion process. This process includes two steps: ① Transformation matrix through the source panel Inverse variable component , Projected onto the global spherical coordinate system, it is converted into momentum in the spherical coordinate system. and The corresponding conversion formula is: ; ② The obtained spherical momentum components are then converted to the inverse representation corresponding to panel 4 to ensure compatibility with the local coordinate system of the target panel. The corresponding conversion formula is: .

[0099] In another exemplary embodiment of this application, in step S6, the boundary flux is first solved using Riemannian methods based on the obtained atmospheric state at the cell boundary points: After preliminary spatial reconstruction and equivalent mapping operations, the atmospheric states at the edge points and center points of all cells are obtained. At this point, noting the temporal tendency of solving for the cell mean, it is necessary to obtain the flux values ​​at the edge points. However, on a single edge, different reconstruction results from cells on both sides appear simultaneously. This embodiment obtains the unique flux value at the boundary point by approximating the Riemann problem. A Low Mach number Approximate Riemann Solver (LMARS) is used to solve for the edge flux.

[0100] The covariate tensor was mentioned earlier. Its inverse matrix They are called inverse variable tensors, and the two are inverse matrices of each other. .

[0101] by Taking direction as an example, the flux at the element boundary is: ; ; This represents the projection coefficient of the pressure gradient force in the x-direction in the momentum equation in the x-direction; This represents the projection coefficient of the pressure gradient force in the y direction in the momentum equation in the x direction; in For pressure-driven flux, other relevant parameters are calculated as follows: ; This represents the diagonal component of the contravariant metric tensor. When i=1, it represents the x-direction; when i=2, it represents the y-direction. This represents the wind speed in the computational space. It is the normal physical wind speed obtained through the Riemann solution. Indicates wind speed. Indicates the wind speed to the left of the unit boundary. This indicates the wind speed to the right of the unit boundary.

[0102] ; ; ; ; ; Indicates the total potential height. Indicates the height of potential. This represents the potential height obtained by reconstructing the left side of the cell boundary. This represents the potential height obtained by reconstructing the right side of the cell boundary; It represents the elevation of the terrain.

[0103] The direction method is similar. Here, the y-direction is indicated when the relevant variable i=2 for G.

[0104] Then, the average flux at the cell boundary is obtained by integrating the results of the edge flux calculation: The purpose of this step is to numerically integrate the flux values ​​obtained by the Riemann solution at each discrete Gaussian integration point (EP) on the cell boundary along the boundary to calculate the average flux on each boundary.

[0105] The boundary flux integral is expanded using a one-dimensional polynomial function: ; Here, i represents the index value of the summation. express The (i-1)th power; k represents the summation index; This is the total number of cells in the template; These are the coefficients of a one-dimensional polynomial; It is a known polynomial basis. .

[0106] Using edge-proximity A linear relationship is established using the flux points on each edge as known conditions. The method for selecting these points is to take symmetrical points on both sides of each edge point. Points. For example, when using TPP5. If unit is required If you accumulate points on the right-hand side, then take the sequence number. Using the five element edge points as known conditions, a system of linear equations can be established, and solving it yields the results. coefficient This uniquely determines the polynomial expression. This process is similar to the integration process in step S51.2, but the dimension is reduced from two dimensions to one dimension.

[0107] Calculation results: The flux information of discrete points is aggregated into an average flux defined on the edges of the cell. For each cell in the grid, an average flux value is calculated for its four edges. This prepares the groundwork for subsequent calculations of the cell's total cost and expenditure.

[0108] In another exemplary embodiment of this application, time tendency calculation is performed in step S7 based on the edge flux results.

[0109] Based on the average source term within the cell and the net flux at the cell boundary, the temporal tendency of the average forecast value within each cell is calculated using the finite volume method formula.

[0110] The process can be summarized as: average source term + average flux of the cell → finite volume discretization → time tendency of each cell →(Passed to step eight: time integration).

[0111] This step requires two key data inputs: the average source term value within the cell and the average flux on each edge of the cell.

[0112] Average of source items within a cell: for each cell The source terms at all Gaussian integration points (CQPs) within its internal cells The element-averaged source term obtained after Gaussian integration .

[0113] Average flux per edge of a cell: for each cell The average flux along its four sides (left, right, bottom, top) .

[0114] This step executes the core discretization process of the finite volume method, namely, calculating the time tendency of the average forecast value within each cell. The relevant formula is as follows: Based on the finite volume method, the time tendency is calculated as follows: ; ; ; ; This indicates the grid spacing in the x-direction; This indicates the grid spacing in the y-direction.

[0115] in Representation unit The area covered Representation unit Left, right, bottom, and top. Integral of source terms within the cell. The method and steps are the same as in step S51.2.

[0116] In another exemplary embodiment of this application, in step S8, the forecast field is obtained by time integration based on the time tendency calculation result.

[0117] The forecast field at time n+1 is obtained using an explicit Runge-Kutta (RK) time integration scheme. ; ; ; This represents the atmospheric state at time n; through intermediate stages ( and The atmospheric state at time n+1 is obtained through calculation. .

[0118] The output of this step indicates the time elapsed after one time step. The latest predictions of atmospheric states (mass distribution and momentum distribution) after calculation are used as the basis for the next round (from arrive 2) Calculate the initial field, and perform iterative calculations.

[0119] In another exemplary embodiment of this application, step S9 performs shadow unit interpolation filling based on the forecast field.

[0120] This step performs a linear mapping operation, the purpose of which is to efficiently and accurately fill in the values ​​of all shadow cells. The output of this step is a complete set of values ​​containing all internal and shadow cell values. Atmospheric conditions at any given time. The next step can be considered a "bridge" step connecting the two time integration steps. It provides complete global data for high-precision calculations in the next time step through an efficient linear transformation.

[0121] This step transforms the shadow cell filling calculation, which is traditionally performed at each time step and may involve iterations, into a one-time pre-calculation of a sparse linear mapping matrix during the initialization phase. Furthermore, it performs efficient sparse matrix-vector multiplication only once per time step during runtime. This significantly improves computational efficiency and is one of the key technologies for achieving high performance and high scalability in this embodiment.

[0122] To achieve the above objectives, this embodiment proposes a method based on pre-computed linear mapping matrices, replacing the traditional inefficient iterative filling algorithm. The complete implementation of this method consists of two stages: a pattern initialization stage (offline computation) and a pattern runtime stage (online computation).

[0123] The front-end input data for this step includes two parts: (1) Static pre-computation data: linear mapping matrix .

[0124] This matrix was completed during the pattern initialization phase mentioned above; the linear mapping matrix was pre-calculated and stored.

[0125] (2) Dynamic input data: derived from the current time in step S8 ( Internal unit forecast field This data comes from the output of step S8 and contains atmospheric state data for all internal units. For each internal computing unit... Step S8 has already calculated its The average predicted values ​​at each time point (including mass and momentum) will be used in this step, based on the average predicted values ​​of these internal units, combined with a pre-calculated linear mapping matrix. By using sparse matrix-vector multiplication, the approximate state values ​​of all shadow units at the current time can be efficiently calculated, thus providing complete global data for the next step of boundary processing or flux calculation.

[0126] The final step is to output numerical weather prediction results based on the solutions of meteorological elements: This step is the final output of the entire forecasting method. Based on the latest forecast field with complete shadow unit filling output in step S9, the numerical solution is converted into meteorological elements required for weather forecasting and the results are output.

[0127] To ensure the stable and reliable operation of the numerical weather prediction method described in this embodiment and to obtain effective forecast results, the entire calculation process must adhere to a fundamental physical-mathematical constraint called the CFL condition (Courant-Friedrichs-Lewy condition). The CFL condition is a numerical stability criterion that specifies a quantitative relationship between the time step and the spatial grid scale when discretely solving partial differential equations (such as the shallow water equations in this application). Its core principle can be stated as follows: within a computational time step, the distance that physical information travels in the real world cannot exceed the size of the discrete grid used to capture that information in the computational grid.

[0128] Based on a fixed time step that satisfies the CFL condition Repeat steps S1-S9 until the desired time is reached to obtain the solution of meteorological elements of the shallow water numerical weather prediction model, namely the results of wind field, height field, mass, energy, angular momentum, potential vorticity, potential vorticity pseudoenergy, vorticity, divergence, etc. at future time. The relevant variables are all conclusions of numerical weather prediction.

[0129] To facilitate understanding, the meteorological element solutions obtained in this step are compared with the weather elements familiar to the public, as shown in Table 1: Table 1. Relationship between Weather Elements

[0130] Therefore, the core physical quantities output in this embodiment, such as "wind field, height field, vorticity, divergence, and potential vorticity," are the cornerstone and direct product of modern numerical weather prediction. These quantities can be directly or indirectly derived from mature meteorological diagnostic principles: Determine the location and intensity of weather systems (such as high pressure, low pressure, and fronts) to judge the weather situation over a wide area (sunny, rainy, windy).

[0131] Diagnosing atmospheric dynamic processes (such as rising / sinking motion) is a key physical basis for determining whether weather phenomena such as precipitation will occur.

[0132] It provides the essential input driving field for subsequent specialized physical parameterization modules (such as cloud precipitation physics and radiation transfer) to ultimately calculate specific public forecast elements such as temperature, humidity, and precipitation.

[0133] In another exemplary embodiment of this application, a high-precision numerical weather prediction system with automatic differentiation is proposed: HOPE (High Order Prediction Environment) is an automatically differentiable, high-precision, non-oscillating, finite-volume shallow-water dynamic core. It uses the finite-volume method to discretize and solve the equations of motion of the spherical barotropic atmosphere, and features high precision, non-oscillation, flux conservation, and high scalability.

[0134] At the algorithm level, the model is numerically discretized based on a cubic sphere mesh, and spatial reconstruction can be performed using high-precision algorithms such as tensor product polynomials or two-dimensional WENO (Weighted Essential Non-Oscillation). At the same time, an intensive shadow unit reconstruction technology has been developed, which can complete the filling of all shadow units with only a single matrix-vector multiplication. In addition, a physically meaningful equivalent mapping scheme is adopted at the intersection of each face of the cubic sphere to achieve strict flux conservation.

[0135] In terms of engineering implementation, it is innovatively developed based on the artificial intelligence software platform PyTorch. Based on the platform's built-in high-performance operators such as convolution and Einstein summation, it achieves cross-platform heterogeneous acceleration capabilities. At the same time, all algorithms can be automatically differentiated, enabling seamless integration with artificial intelligence models.

[0136] Process Overview: The equations of motion for a spherical barotropic atmosphere are multivariable coupled nonlinear systems in infinite-dimensional space, and cannot be solved analytically. They must be solved numerically through discretization. The core task of the HOPE shallow-water dynamics framework is to discretize and solve the equations of motion for a barotropic atmosphere. The specific steps of the model are as follows: Please see Figure 9 Step 1: Obtain initial field data and initialize the system.

[0137] The shallow water equations based on atmospheric motion (including mass conservation and momentum conservation equations) are initialized once before the model begins integration to drive the numerical weather prediction model and provide a complete computational framework and static data for subsequent numerical solutions.

[0138] Acquire initial field data and parameters: Obtain initial atmospheric state data (global geopotential height field and wind field), topographic data (surface elevation data), and set physical parameters such as gravitational acceleration and Earth's rotation angular velocity for solving the atmospheric motion equations.

[0139] Pre-compute and store the linear mapping matrix: To efficiently and stably handle cross-panel boundary calculations, an automatic differentiation framework is used to pre-compute the linear mapping matrix from interior cell forecasts to shadow cell filling values ​​and store it in sparse matrix format. This matrix is ​​reused in all subsequent time steps, transforming boundary filling into efficient sparse matrix-vector multiplication.

[0140] Generating a spherical cubic mesh: A spherical computational domain is established, specifically by generating an isoangular cubic sphere mesh. The sphere is decomposed into six faces, each mapped to a rectangular computational space for subsequent discretization of the equations. This mesh structure provides a discrete framework for solving systems of partial differential equations, ensuring mesh uniformity and computational efficiency.

[0141] Step 2: Generate the initial field at time n based on the cubic sphere mesh from Step 1.

[0142] Select a scenario corresponding to the shallow water equations, such as geostrophic flow, hill-crossing airflow, Rossby-Haurwitz waves, or actual atmospheric conditions at 500 hPa. Generate atmospheric conditions using a matching cubic sphere grid based on the initial value function or the original data. This step provides the model with initial atmospheric conditions, offering initial conditions for solving the partial differential equations and serving as the starting point for subsequent calculations.

[0143] Step 3: Based on the initial field in Step 2, perform spatial reconstruction to obtain the cell boundary values.

[0144] Spatial reconstruction can be achieved using either the Tensor Product Polynomial (TPP) scheme or the Two-Dimensional Weighted Fundamentally Oscillated (WENO 2D) scheme, thereby obtaining the internal integration point values ​​on the cell edges for subsequent Riemann solutions and source term calculations. The WENO 2D scheme achieves oscillation-free properties, but the computational cost increases due to the presence of nonlinear weights, making it suitable for flow fields with poor continuity. The TPP scheme, while simple and computationally efficient, cannot adaptively suppress non-physical oscillations and is suitable for smooth flow fields. This step ensures computational accuracy and robustness through high-precision reconstruction and an adaptive oscillation-free algorithm.

[0145] Step 4: Perform equivalent mapping of the cubic sphere boundary based on the spatial reconstruction results in Step 3.

[0146] At the integration points on the boundaries of each face of the cubic sphere, the mass and momentum fields are remapped, mapping the information of neighboring faces to the coordinate system of this face, thus unifying the atmospheric state at the boundaries of the cubic sphere mesh. This step ensures the physical equivalence of the atmospheric state at the boundaries, ensuring flux conservation when the flow field crosses the face boundaries, and guaranteeing strict conservation of mass and momentum during subsequent time integration. Combining this step with the element boundary values ​​obtained in step three, the atmospheric state at all element boundary points of the cubic sphere mesh has been obtained.

[0147] Step 5: Solve the Riemann problem to obtain the edge flux based on the atmospheric state of the unit boundary points obtained in Steps 3 and 4.

[0148] The low Mach number approximate Riemann solver (LMARS) is used to calculate the element boundary flux. This step obtains the flux at the element interface through the approximate Riemann solution, providing necessary data for subsequent spatial dip calculations.

[0149] Step 6: Calculate the time tendency based on the edge flux results in Step 5.

[0150] The net flux of the cell is calculated based on the flux values ​​at each interface, and the source terms within the cell are also calculated. The tendency of the cell mean to change over time can then be obtained using the finite volume method. This time tendency is the sum of the negative of the net flux at the cell boundary and the source terms within the cell. This step, by calculating the net flux at the cell interface and the source terms within the cell, provides a basis for time integration to obtain the tendency of the cell mean to change over time.

[0151] Step 7: Based on the time tendency calculation results in Step 6, perform time integration to obtain the forecast field.

[0152] The forecast field at time n+1 is obtained using a third-order explicit Runge-Kutta time integration scheme. This step, through time progression, yields the atmospheric state at the forecast time.

[0153] Step 8: Perform shadow unit interpolation filling based on the forecast field in Step 7.

[0154] A central reconstruction scheme is adopted, using an iterative algorithm to fill in the information of unknown cells. Then, based on the linear characteristics of the mapping algorithm, the interpolation mapping matrix is ​​obtained using the model's automatic differentiability property, and this matrix is ​​saved in CSR sparse matrix format. This process only needs to be performed once during model initialization. During model integration, the entire filling process is directly expressed as the multiplication of the mapping matrix and the atmospheric state vector, eliminating the need for iterative calculations, which greatly improves computational efficiency and reduces the number of communications between distributed computing nodes. This step improves computational efficiency by optimizing the handling of shadow cells while ensuring the high accuracy of the algorithm.

[0155] Step 9: Output weather forecast results based on the solution of meteorological elements Through steps one through eight, the solutions for meteorological elements in the shallow water numerical weather prediction model are obtained, and thus the weather forecast results are obtained.

[0156] HOPE is a high-precision, oscillatory-free, finite-volume shallow-water dynamic core that is automatically differentiable. This application utilizes tensor product polynomials to construct a class of convolution-like reconstruction algorithms that efficiently leverage the heterogeneous computing capabilities of GPUs. It also provides convenient compatibility with nonlinear, oscillatory algorithms such as the two-dimensional WENO algorithm. This efficiently achieves high-precision, oscillatory characteristics for the dynamic framework of numerical models. It pioneers the development of high-precision numerical weather prediction models based on artificial intelligence (AI) software platforms. The automatic differentiation function of the PyTorch platform enables seamless integration of HOPE with AI models. Furthermore, fully utilizing PyTorch's high-performance built-in functions fully leverages the performance advantages of heterogeneous computing platforms and facilitates cross-platform development. Based on the results of this research, it can be easily extended to baroclinic atmospheric models and combined with AI physical process parameterized models to form a fusion-type numerical weather prediction model. This model maintains the high-precision characteristics of traditional numerical methods while possessing the automatic learning capabilities of AI physical processes.

[0157] Problems with existing technology: Numerical weather forecasting models that possess all five characteristics—high accuracy, no oscillation, scalability, mass conservation, and automatic differentiation—are currently rare. Existing numerical models typically only possess a few of these characteristics and cannot possess all of them. Existing numerical models typically adopt an optimization strategy of first developing them on the CPU x86 platform and then porting them to heterogeneous platforms, without taking into account the characteristics of heterogeneous computing in the top-level design. This limits the acceleration performance that the models can achieve on heterogeneous platforms. Most existing numerical weather prediction models are developed based on Fortran, lack automatic differentiability, and are inconvenient to integrate with AI models.

[0158] The HOPE proposed in this application aims to address the following key issues: Improve accuracy and stability: Existing problem: Traditional high-order numerical methods are prone to non-physical oscillations in non-smooth flow fields (such as abrupt terrain changes or regions with strong gradients), leading to simulation distortion.

[0159] Solution: By employing arbitrarily high-order precision algorithms, the trade-off between accuracy and stability inherent in traditional methods is resolved. While ensuring mass conservation, simulation errors are effectively reduced, significantly improving accuracy, especially when dealing with smooth flow fields. Simultaneously, by utilizing schemes such as WENO, the model remains stable when handling non-smooth regions (such as shock waves and discontinuous fields), avoiding the generation of non-physical oscillations.

[0160] Enhance computational efficiency and scalability: Existing problems: Higher-order methods suffer from low parallel efficiency and difficulty in adapting to GPUs on cubic sphere meshes.

[0161] Solution: By optimizing numerical algorithms and data structures, computational efficiency is significantly improved while maintaining accuracy. The algorithm design fully considers the characteristics of GPU architecture, transforming the main operations into convolution and matrix multiplication, facilitating accelerated computation using GPUs and greatly improving computational efficiency. Furthermore, the model architecture possesses good scalability, capable of adapting to increasing computational demands and supporting future higher-resolution numerical simulations.

[0162] Achieving automatic differentiation and flexible coupling: Existing problem: Traditional numerical models lack automatic differentiation capabilities, making it difficult to combine them with neural networks to build hybrid models.

[0163] Solution: HOPE was developed based on PyTorch, giving it inherent automatic differentiation capabilities. This greatly simplifies the development of mode adjoint operators in four-dimensional variational assimilation, facilitating the development of mode assimilation systems. Full-mode automatic differentiation enables end-to-end training, allowing direct optimization of physically parameterized models using multivariate data, eliminating reliance on reanalysis data.

[0164] It is worth noting that the core advantage of this application lies in its significantly simplified integration architecture, reducing nonlinear computation and memory consumption. This application proposes an optimized scheme for high-precision finite-volume single-point integration, which deeply simplifies the integration point architecture: Compared with the original HOPE algorithm, this application utilizes a polynomial basis to construct a numerical integration scheme, overcoming the problem of excessive sampling points in Gaussian integration. Only a single sampling point is needed at both the cell boundary and the center to achieve high-precision integration. This optimized scheme significantly reduces the number of reconstruction points, integration points, and storage overhead while maintaining arbitrary high-order precision, oscillation-free operation, strict mass conservation, and native automatic differentiability. This further improves GPU parallel computing efficiency and is more suitable for ultra-high-resolution numerical weather forecasting scenarios.

[0165] While maintaining the same accuracy, the nonlinear computation requirements, reconstruction computation costs, and storage requirements are all significantly reduced. A 5th-order accuracy scheme (TPP5) was used to simulate the steady-state flow field of geostrophic flow. The errors and convergence rates of the original multi-point algorithm and the single-point algorithm of this application were compared, and the results are shown in Table 2.

[0166] Table 2. Comparison of geostrophic flow accuracy tests between the original multi-point scheme and the single-point scheme of this application at different resolutions. In Table 2, C30, C45, and C90 represent three resolutions of the cubic sphere grid, corresponding to each panel being divided into 30×30, 45×45, and 90×90 grid units, respectively, with the grid ranging from coarse to fine, and C90 being the finest.

[0167] L1 norm error ( L2 norm error L-infinity norm error ( These three errors are standard indicators for evaluating the accuracy of numerical weather prediction models. They measure the deviation between simulation results and the true solution (analytical solution) from different perspectives, each with its own emphasis. The smaller all three errors are, the better. Their calculation formulas are as follows:

[0168] in, Indicated in the panel Above, grid cells The predicted values ​​obtained from numerical simulation; This represents the corresponding analytical solution (true solution); It is a weighted summation operator; The Jacobian determinant, representing the metric of a cubic sphere mesh, is used to transform the computational space into the physical space.

[0169] The numerator of norm error This means that for all grid cells, the difference between the predicted value and the actual solution is multiplied by the corresponding cell area weight. Then sum them; denominator This represents the summation of the true solution by area weight.

[0170] The numerator of norm error , representing the square root of the square of the deviation on each grid cell; denominator This represents the square root of the weighted sum of the squares of the true solution.

[0171] The numerator of norm error , representing the maximum absolute deviation on the grid cell; denominator This represents the maximum absolute value of the true solution.

[0172] Convergence rate is a core metric for measuring the accuracy of numerical algorithms, reflecting the rate at which simulation errors decrease when the computational grid is refined. For the 5th-order accuracy scheme (TPP5) tested in this application, the theoretical convergence rate should be 5.0.

[0173] As shown in Table 2, during the mesh refinement process (C30→C45→C90), the convergence rates of the single-point scheme and the original multi-point scheme of this application are highly consistent: (1) C30 to C45 refinement: The convergence rates of L1, L2 and L infinity of the single-point scheme of this application are 5.0038, 4.9953 and 4.9663, respectively, while those of the original multi-point scheme are 5.0039, 4.9954 and 4.9652, respectively; (2) C45 to C90 refinement: The convergence rates of L1, L2 and L infinity of the single-point scheme of this application are 5.0013, 4.9989 and 4.9922, respectively, while those of the original scheme are 5.0014, 4.9992 and 4.9923, respectively.

[0174] Both convergence rates are close to the theoretical value of 5.0, indicating that the single-point scheme in this application strictly maintains the 5th-order convergence characteristic. (Note: The slight deviation between the convergence rate and the theoretical value of 5.0 is normal and is mainly caused by mesh non-uniformity and asymptotic convergence behavior under finite resolution. It does not affect the qualitative conclusion on the accuracy of the algorithm.) The reduction in storage and computational load of the new HOPE algorithm is shown in Table 3.

[0175] Table 3. Reduction in storage and computational load of the new HOPE algorithm

[0176] It possesses five characteristics: high precision, no oscillation, scalability, mass conservation, and support for automatic differentiation.

[0177] Achieving high precision of arbitrary order: HOPE proposes a high-order interpolation scheme for cubic sphere meshes, which achieves arbitrary odd-order precision through tensor product polynomial (TPP) reconstruction in step 3.

[0178] Oscillation-free: This scheme overcomes the shortcomings of traditional spectral methods (such as NeuralGCM) in terms of mass conservation and numerical oscillation. It introduces two-dimensional WENO (WENO 2D) adaptive reconstruction, which automatically reduces the order in shock wave or steep gradient regions, suppresses non-physical oscillations, achieves essentially oscillation-free operation, and has numerical stability.

[0179] Scalability: This algorithm has high scalability on GPUs, has no global dependencies, requires only local communication, and supports large-scale parallel computing (such as PyTorch version adapted to GPU clusters), breaking through the scalability limitations of traditional spectral methods.

[0180] Strict mass conservation: Based on the finite volume method (FVM), a conservation algorithm is designed to naturally guarantee global mass conservation, thus solving the problem of mass non-conservation caused by global spectral expansion in spectral models.

[0181] Supports automatic differentiation: Implemented based on the PyTorch framework, it natively supports automatic differentiation. In scenarios that require residual backpropagation calculations, such as joint training with machine learning models and four-dimensional variational assimilation, there is no need to manually develop an adjoint model.

[0182] Algorithm design adapted to GPU architecture to achieve performance innovation GPU-friendly design: Spatial reconstruction is equivalent to convolution operation (torch.nn.Functional.conv2d), flux integration is implemented through matrix multiplication, which is fully vectorized, avoids loops and conditional branches, and significantly improves computational efficiency compared to Fortran code.

[0183] Memory and computation optimization: The shadow unit interpolation in step eight reduces the amount of real-time computation by pre-computing the linear mapping matrix, thus avoiding memory pressure during the iteration process.

[0184] An automated differentiable framework for AI convergence has the advantage of ecosystem compatibility. This embodiment is implemented based on the PyTorch framework, supports automatic differentiation, and facilitates the coupling of arbitrary neural network (NN) methods. Using PyTorch as a bridge, it connects numerical models with AI toolchains, laying the foundation for building a fusion "numerical + AI" forecasting system.

[0185] The PyTorch platform offers the following advantages for engineering development: First, it enables cross-platform high-performance computing using PyTorch's built-in high-performance functions; second, its full-mode automatic differentiation makes end-to-end training possible, allowing direct optimization of physically parameterized models using multivariate data, eliminating reliance on reanalysis data; third, its automatic differentiation function greatly simplifies the development of mode adjoint operators in four-dimensional variational assimilation, facilitating the development of mode assimilation systems; and fourth, it allows for seamless integration with AI models without requiring the development of additional cross-language interface modules.

[0186] In one exemplary embodiment, a computer device is provided, which may be a server or a terminal, and its internal structure diagram may be as follows. Figure 10As shown, this computer device includes a processor, memory, input / output (I / O) interfaces, and a communication interface. The processor, memory, and I / O interfaces are connected via a system bus, and the communication interface is also connected to the system bus via the I / O interfaces. The processor provides computational and control capabilities. The memory includes non-volatile storage media and internal memory. The non-volatile storage media stores the operating system, computer programs, and databases. The internal memory provides the environment for the operation of the operating system and computer programs stored in the non-volatile storage media. The I / O interfaces are used for exchanging information between the processor and external devices. The communication interface is used for communication with external terminals via a network connection. When the computer program is executed by the processor, it implements an optimization method for a high-precision numerical weather prediction method with automatic differentiation.

[0187] Those skilled in the art will understand that Figure 10 The structures shown are merely block diagrams of some structures related to the present application and do not constitute a limitation on the computer device to which the present application is applied. Specific computer devices may include more or fewer components than shown in the figures, or combine certain components, or have different component arrangements. In an exemplary embodiment, a computer device is provided, including a memory and a processor. The memory stores a computer program, and the processor executes the computer program to implement the steps in the above-described method embodiments.

[0188] In one exemplary embodiment, a computer-readable storage medium is provided storing a computer program that, when executed by a processor, implements the steps in the above-described method embodiments.

[0189] In one exemplary embodiment, a computer program product is provided, including a computer program that, when executed by a processor, implements the steps in the above-described method embodiments.

[0190] It should be noted that the user information (including but not limited to user device information, user personal information, etc.) and data (including but not limited to data used for analysis, data stored, data displayed, etc.) involved in this application are all information and data authorized by the user or fully authorized by all parties. Moreover, the collection, use and processing of the relevant data are carried out in compliance with the relevant data protection laws and policies of the country where the location is located, and with the authorization granted by the owner of the corresponding device.

[0191] Those skilled in the art will understand that all or part of the processes in the above embodiments can be implemented by a computer program instructing related hardware. The computer program can be stored in a non-volatile computer-readable storage medium. When executed, the computer program can include the processes of the embodiments described above. Any references to memory, databases, or other media used in the embodiments provided in this application can include at least one of non-volatile and volatile memory. Non-volatile memory can include read-only memory (ROM), magnetic tape, floppy disk, flash memory, optical memory, high-density embedded non-volatile memory, resistive random access memory (ReRAM), magnetic random access memory (MRAM), ferroelectric random access memory (FRAM), phase change memory (PCM), graphene memory, etc. Volatile memory can include random access memory (RAM) or external cache memory, etc. By way of illustration and not limitation, RAM can take many forms, such as Static Random Access Memory (SRAM) or Dynamic Random Access Memory (DRAM).

[0192] The databases involved in the embodiments provided in this application may include at least one type of relational database and non-relational database. Non-relational databases may include, but are not limited to, blockchain-based distributed databases. The processors involved in the embodiments provided in this application may be general-purpose processors, central processing units, graphics processing units, digital signal processors, programmable logic devices, quantum computing-based data processing logic devices, etc., and are not limited to these.

[0193] The technical features of the above embodiments can be combined in any way. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0194] This document uses specific examples to illustrate the principles and implementation methods of this application. The descriptions of the above embodiments are only for the purpose of helping to understand the methods and core ideas of this application. Furthermore, those skilled in the art will recognize that, based on the ideas of this application, there will be changes in the specific implementation methods and application scope. Therefore, the content of this specification should not be construed as a limitation of this application.

Claims

1. An optimization method for a high-precision numerical weather prediction method with automatic differentiation, characterized in that, include: Acquire initial atmospheric state data and topographic data, and configure physical parameters based on geophysical constants; Based on the physical parameters, a cubic sphere mesh is generated and the sphere is decomposed into six panels; on each panel, internal cells for solving the forecast quantity and shadow cells located outside the panel boundary for handling boundary conditions are defined. Based on the geometry of the cubic sphere mesh, a linear mapping matrix from the internal cells to the shadow cells is pre-calculated; The linear mapping matrix is ​​stored in a sparse matrix format; The initial atmospheric state data is mapped onto the cubic sphere grid to generate the initial field at time n, and the source term average value of each of the internal cells is calculated based on the terrain data and the physical parameters; the initial field at time n includes the internal cell mean value stored on each of the internal cells; Based on the mean of the internal units, spatial reconstruction is performed on each internal unit. The distribution of physical quantities inside the unit is constructed using tensor product polynomials. The integration point value on the boundary of the internal unit is obtained by using the midpoint rule as a single sampling point. Using the integration point values ​​on the boundary of the inner element as the left and right states of the Riemann problem, the Riemann problem is solved using a low Mach number approximate Riemann solver to obtain the flux value at each integration point on the boundary of the inner element. Numerical integration is then performed along each boundary of the inner element to obtain the average flux on each boundary of the inner element. Based on the average flux on each boundary of the internal unit and the average source term of the internal unit, the time tendency of the mean of the internal unit with time is calculated according to the finite volume method; the average source term is obtained by integration with the center point of the unit as a single sampling point through the midpoint rule. Based on the aforementioned time tendency, an explicit Runge-Kutta time integration scheme is used to perform time integration, and the internal cell prediction field at time n+1 is obtained. Using the linear mapping matrix, for the The internal cell prediction field at time 1 is subjected to sparse matrix-vector multiplication to calculate the filling values ​​of all shadow cells, thus obtaining a prediction field containing complete shadow cells. Time forecast field.

2. The optimization method for the high-precision numerical weather prediction method with automatic differentiation according to claim 1, characterized in that, Based on the mean of the internal units, spatial reconstruction is performed on each internal unit. The distribution of physical quantities within the unit is constructed using a tensor product polynomial. Then, using the midpoint rule, integration is performed with the midpoint of the unit boundary as a single sampling point to obtain the integral point value on the boundary of the internal unit. Specifically, this includes: For any given internal cell, take that internal cell as the target internal cell, and in the computational space of the cubic sphere mesh, select... The neighboring units form a square reconstruction template; where n is determined according to the required reconstruction order, and n is an odd number. Extract the average volume fraction of all cells within the reconstructed template; Based on the volume integral average, a system of linear equations is constructed with tensor product polynomials as basis functions; the tensor product polynomials consist of the tensor products of two one-dimensional polynomials in the x and y directions; the form of the tensor product polynomials is: ;in, and It is a two-dimensional row and column index in the template grid. This represents the index value in the x-direction within the template, ranging from 1 to n. This represents the index value in the y-direction within the template, ranging from 1 to n, where n represents the template width; and Represents a polynomial with the corresponding power; For polynomial coefficients, is the local coordinate system of the template; k is the calculation parameter. ; These are basis functions in tensor product polynomial reconstruction. ; and For calculation parameters; Solving the linear equations yields all the coefficients of the tensor product polynomial, thereby uniquely determining the physical quantity distribution function within the target's internal unit; the physical quantity distribution function is a continuously differentiable function used to describe the continuous changes in potential height and / or momentum within the unit; On each boundary of the target internal unit, the function value at the midpoint of the boundary is calculated using the physical quantity distribution function, and the function value at the midpoint of the boundary is used as the representative value of the physical quantity on the boundary; the representative value is multiplied by the side length of the boundary to obtain the integral value on the boundary; all boundaries of the target internal unit are traversed to obtain the integral point value on each boundary of the target internal unit. Traverse all internal cells to obtain the integral point values ​​on the boundaries of all internal cells.

3. The optimization method for the high-precision numerical weather prediction method with automatic differentiation according to claim 1, characterized in that, Based on the geometry of the cubic sphere mesh, a linear mapping matrix from the internal cells to the shadow cells is pre-calculated, specifically including: Initialize all shadow units to zero; In each iteration, based on the currently known mean of the internal cells and the value of the shadow cells in the current iteration step, the tensor product polynomial interpolation method is used to calculate the value of all Gaussian integral points in each shadow cell. Perform Gaussian integration on the values ​​of all Gaussian integration points within each shadow unit to obtain the new average value of the shadow unit, thereby obtaining the updated values ​​of all shadow units; The convergence state is reached when the norm of the difference between the values ​​of the shadow units obtained from two consecutive iterations is less than a preset threshold. The values ​​of the shadow units in the convergent state are regarded as linear functions of the mean of the inner units. The derivative of the linear function is obtained by using the automatic differentiation framework to obtain the Jacobian matrix from the inner units to all shadow units, and the Jacobian matrix is ​​used as the linear mapping matrix.

4. The optimization method for the high-precision numerical weather prediction method with automatic differentiation according to claim 1, characterized in that, Using the integration point values ​​on the boundary of the inner element as the left and right states of the Riemann problem, a low Mach number approximate Riemann solver is used to solve the Riemann problem, obtaining the flux value at each integration point on the boundary of the inner element, specifically including: From the integral point values ​​on the boundary of the internal unit, extract the mass field value and momentum field value at the integral point located at the boundary of each panel of the cube sphere, and perform coordinate system physical quantity remapping on the mass field value and momentum field value at the boundary to obtain the remapping result; Using the remapping result as the left and right states of the Riemann problem, the Riemann problem is solved using a low Mach number approximate Riemann solver to obtain the flux value at each integration point on the boundary of the inner cell.

5. The optimization method for the high-precision numerical weather prediction method with automatic differentiation according to claim 4, characterized in that, At the boundary, the physical quantities of the mass field value and the momentum field value are remapped in a coordinate system to obtain the remapping result, which specifically includes: For each Gaussian integration point located on the panel boundary, obtain its mass field value and inverse momentum component obtained by spatial reconstruction on the source panel; the inverse momentum component is the momentum component expressed in the local computational coordinate system of the source panel. Based on the source panel coordinate transformation matrix and target panel coordinate transformation matrix pre-calculated during the cubic sphere mesh generation process, and the metric Jacobian determinant of the source panel and target panel at the Gaussian integration point; The mass field value in the target panel coordinate system is obtained by scaling and transforming the mass field value using the formula: Target panel mass field value = (Target panel metric Jacobian ÷ Source panel metric Jacobian) × Source panel mass field value. The inverse momentum components are converted into momentum components in spherical coordinates by matrix multiplication of the source panel coordinate transformation matrix and the inverse momentum components. By multiplying the inverse matrix of the target panel coordinate transformation matrix with the matrix of the spherical momentum components, and then multiplying by the ratio of the metric Jacobian, the momentum components in the spherical coordinate system are transformed to the contravariant representation of the target panel, thus obtaining the remapping result.

6. The optimization method for the high-precision numerical weather prediction method with automatic differentiation according to claim 1, characterized in that, After using the linear mapping matrix to perform sparse matrix-vector multiplication on the internal cell prediction field at time n+1, calculating the filling values ​​of all shadow cells, and obtaining the prediction field at time n+1 containing complete shadow cells, the process further includes: Based on the forecast field at time n+1 containing complete shadow units, the solution of meteorological elements is output as the numerical weather prediction result; the solution of meteorological elements includes one or more of wind field, height field, vorticity, divergence, potential vorticity, mass conservation quantity, and energy conservation quantity.

7. The optimization method for the high-precision numerical weather prediction method with automatic differentiation according to claim 1, characterized in that, When solving this Riemann problem using a low Mach number approximation Riemann solver, the solution... The formula for the flux at the element boundary in the direction is: ； Where P is the pressure-driven flux, m is the calculated parameter, and q L and q R These represent the boundary values ​​on the left and right sides of the cell, respectively.

8. A computer device, comprising: A memory, a processor, and a computer program stored in the memory and capable of running on the processor, characterized in that the processor executes the computer program to implement the optimization method for the high-precision numerical weather forecasting method with automatic differentiation as described in any one of claims 1-7.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that, When executed by a processor, the computer program implements the optimization method for the high-precision numerical weather forecasting method with automatic differentiation as described in any one of claims 1-7.

10. A computer program product, comprising a computer program, characterized in that, When executed by a processor, the computer program implements the optimization method for the high-precision numerical weather forecasting method with automatic differentiation as described in any one of claims 1-7.

Images